# Asymptotic, second-order homogenization of linear elastic beam networks

Yang Ye, Basile Audoly, Claire Lestringant

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yang Ye, Basile Audoly, Claire Lestringant. Asymptotic, second-order homogenization of linear elastic beam networks. Journal of the Mechanics and Physics of Solids, inPress, 10.1016/j.jmps.2024.105637. hal-04549113v2

HAL Id: hal-04549113

https://hal.science/hal-04549113v2

Submitted on 17 Apr 2024

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.

<figure><img src="image_2.png" /></figure>
# Asymptotic, second-order homogenization of linear elastic beam networks

YANG YE

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

B. AUDOLY

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

C. LESTRINGANT

Institut Jean Le Rond d'Alembert,
Sorbonne Université, CNRS,
75005 Paris, France

April 17, 2024

We propose a general approach to the higher-order homogenization of discrete elastic networks made up of linear elastic beams or springs in dimension 2 or 3. The network may be *nearly* (rather than *exactly*) periodic: its elastic and geometric properties are allowed to vary slowly in space, in addition to being periodic at the scale of the unit cell. The reference configuration may be prestressed. A homogenized strain energy depending on both the macroscopic strain $\varepsilon$ and its gradient $\nabla\varepsilon$ is obtained by means of a two-scale expansion. The homogenized energy is asymptotically exact two orders beyond that obtained by classical homogenization. The homogenization method is implemented in a symbolic calculation language and applied to various types of networks, such as a 2D honeycomb, a 2D Kagome lattice, a 3D truss and a 1D pantograph. It is validated by comparing the predictions of the microscopic displacement to that obtained by full, discrete simulations. This second-order method remains highly accurate even when the strain gradient effects are significant, such as near the lips of a crack tip or in regions where a gradient of pre-strain is imposed.

**Keywords** Elastic lattices, Asymptotic homogenization, Second-order homogenization, Energy methods.

## 1. INTRODUCTION

Man-made materials have evolved over time, achieving a constant improvement in properties such as strength, stiffness, lightweight or wave transmission. Some of the more recent breakthroughs were made possible by designing the architecture at small scale. Progress in additive manufacturing enables the fabrication of cellular materials whose microstructure can be precisely controlled and tailored: a wide range of new architectured materials can be fabricated, using periodic arrangements of thin beams as fundamental structural elements and ceramics, metals, polymers or composites as base materials [FDA10]. The ongoing revolution consists in integrating these materials into structures having optimized properties.

The design of such systems calls for effective *homogenized* models capable of precisely capturing the mechanical behavior of these complex periodic or quasi-periodic microstructures [AG97; CP12; VP12]. Recent work shows that the accuracy of Cauchy-type continuum obtained by classical homogenization techniques is limited. In various cases, the effective models need to include higher-order terms in order to accurately capture physical phenomena such as localization in shear-bands or domain walls [GLTK19; DYF+20], mechanisms generating long-range modes of deformation [NCH20; DLSS22; ] or symmetries that are not present in the low-order models [RA16; RA19]. Besides, higher-order terms improve the accuracy of the effective models, especially when scale separation is poor, for example when modeling a crack in a relatively coarse microstructure [RKD+15; RDK17; SCO+22].

Various strategies have been proposed to derive higher-order effective models: numerical identification procedures [VDP14; GLTK19; RDVB19a; RDVB19b; RAPG19] or analytical approaches based on formal derivations that slave local displacements to macroscopic fields through a Taylor expansion [TB93; BT94; KM04]. The two-scale expansion method offers a rigorous methodology to derive effective models without introducing such ad-hoc kinematic assumptions. It has been applied to continuous periodic composites [Bou96; SC00; Bou19] and discrete microstructures [CMR06; DO11; DG12; Dav13; LR13], including for deriving higher-order contributions [AS18b; AS18a; ASB19; AB21; DLSS22].

In a recent paper [AL23], we revisited higher-order asymptotic homogenization of linear elastic, discrete microstructures. Building upon our previous work on asymptotic dimension reduction [LA20] and following [LM18], we proposed to tackle asymptotic homogenization at the energy level. Starting from an abstract, generic energy formulation of the discrete microstructure, we introduced a series of steps that can be applied in an automated manner to derive a homogenized energy. All the steps are implemented in a symbolic calculation language and distributed as the open-source library shoal (for Second-order HOMogenization Automated in a Library) [Aud23]. This approach can handle microstructures featuring both pre-stress and spatially varying (graded) properties, two aspects that are rarely addressed in the context of high-order, asymptotic homogenization and are very useful for designing and optimizing microstructures.
Asymptotic, second-order homogenization
of linear elastic beam networks

YANG YE

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

B. AUDOLY

Laboratoire de Mécanique des Solides
CNRS, Institut Polytechnique de Paris
91120 Palaiseau, France

C. LESTRINGANT

Institut Jean Le Rond d'Alembert,
Sorbonne Université, CNRS,
75005 Paris, France

April 17, 2024

We propose a general approach to the higher-order homogenization of discrete elastic networks made up of linear elastic beams or springs in dimension 2 or 3. The network may be *nearly* (rather than *exactly*) periodic: its elastic and geometric properties are allowed to vary slowly in space, in addition to being periodic at the scale of the unit cell. The reference configuration may be prestressed. A homogenized strain energy depending on both the macroscopic strain $\varepsilon$ and its gradient $\nabla\varepsilon$ is obtained by means of a two-scale expansion. The homogenized energy is asymptotically exact two orders beyond that obtained by classical homogenization. The homogenization method is implemented in a symbolic calculation language and applied to various types of networks, such as a 2D honeycomb, a 2D Kagome lattice, a 3D truss and a 1D pantograph. It is validated by comparing the predictions of the microscopic displacement to that obtained by full, discrete simulations. This second-order method remains highly accurate even when the strain gradient effects are significant, such as near the lips of a crack tip or in regions where a gradient of pre-strain is imposed.

**Keywords** Elastic lattices, Asymptotic homogenization, Second-order homogenization, Energy methods.

# 1. INTRODUCTION

Man-made materials have evolved over time, achieving a constant improvement in properties such as strength, stiffness, lightweight or wave transmission. Some of the more recent breakthroughs were made possible by designing the architecture at small scale. Progress in additive manufacturing enables the fabrication of cellular materials whose microstructure can be precisely controlled and tailored: a wide range of new architECTured materials can be fabricated, using periodic arrangements of thin beams as fundamental structural elements and ceramics, metals, polymers or composites as base materials [FDA10]. The ongoing revolution consists in integrating these materials into structures having optimized properties.

The design of such systems calls for effective *homogenized* models capable of precisely capturing the mechanical behavior of these complex periodic or quasi-periodic microstructures [AG97; CP12; VP12]. Recent work shows that the accuracy of Cauchy-type continuum obtained by classical homogenization techniques is limited. In various cases, the effective models need to include higher-order terms in order to accurately capture physical phenomena such as localization in shear-bands or domain walls [GLTK19; DYF+20], mechanisms generating long-range modes of deformation [NCH20; DLSS22; ] or symmetries that are not present in the low-order models [RA16; RA19]. Besides, higher-order terms improve the accuracy of the effective models, especially when scale separation is poor, for example when modeling a crack in a relatively coarse microstructure [RKD+15; RDK17; SCO+22].

Various strategies have been proposed to derive higher-order effective models: numerical identification procedures [VDP14; GLTK19; RDVB19a; RDVB19b; RAPG19] or analytical approaches based on formal derivations that slave local displacements to macroscopic fields through a Taylor expansion [TB93; BT94; KM04]. The two-scale expansion method offers a rigorous methodology to derive effective models without introducing such ad-hoc kinematic assumptions. It has been applied to continuous periodic composites [Bou96; SC00; Bou19] and discrete microstructures [CMR06; DO11; DG12; Dav13; LR13], including for deriving higher-order contributions [AS18b; AS18a; ASB19; AB21; DLSS22].

In a recent paper [AL23], we revisited higher-order asymptotic homogenization of linear elastic, discrete microstructures. Building upon our previous work on asymptotic dimension reduction [LA20] and following [LM18], we proposed to tackle asymptotic homogenization at the energy level. Starting from an abstract, generic energy formulation of the discrete microstructure, we introduced a series of steps that can be applied in an automated manner to derive a homogenized energy. All the steps are implemented in a symbolic calculation language and distributed as the open-source library shoal (for Second-order HOMogenization Automated in a Library) [Aud23]. This approach can handle microstructures featuring both pre-stress and spatially varying (graded) properties, two aspects that are rarely addressed in the context of high-order, asymptotic homogenization and are very useful for designing and optimizing microstructures.
The present paper addresses higher-order homogenization of periodic networks made up of slender 1D elements (typically beams or springs), further referred to as *lattice material* or simply *lattice*. The homogenization is done symbolically and in an automated way, and builds up on the general-purpose homogenization engine from our previous work [AL23]. It delivers closed-form analytical expressions of the effective elastic properties without any adjustable internal kinematic variable or parameter. By contrast with earlier work focussing on particular lattice architectures [DO11; Dav13; LR13] or 1D/2D networks of springs [MA02; MA06], our approach is designed to be general and versatile. It works in 1D, 2D or 3D and can handle arbitrary lattice topologies, the presence of prestress, slowly varying geometrical or elastic properties, a variety of member types (such as springs, beams or curved arches), as well as a large elastic contrast.

We illustrate our homogenization method on a comprehensive set of examples: a uniform 2D honeycomb, spatially varying 2D Kagome and honeycomb architectures, a 3D quasi-octet lattice and 1D pantograph. With a view of assessing the accuracy of the resulting models, we propose a numerical verification procedure that compares the microscopic fluctuations extracted from full numerical simulations to those predicted by the homogenization method. We apply this procedure to a range of structural problems including macroscopic samples with cracks, holes, as well as spatially varying pre-stress, elastic and geometric properties. Our numerical results show that the higher-order effective model very accurately captures the fine features of the local displacements and rotations away from the boundaries, where the presence of boundary layers is known to make homogenization break down.

Higher-order homogenization comes with several difficulties, and some remain beyond the scope of the present work. Most notably, we do not address the presence of boundary layers, and accept that the gradient stiffness is often negative (see Remark 5.3 below). These limitations prevent from using the homogenized energy to set up self-contained simulations so far, as discussed in further detail in the conclusion.

## 2. GENERATING THE ELASTIC LATTICE

In this section, we describe the process of generating the elastic lattices. Given a small parameter $\eta \ll 1$, we produce a lattice having cells of size $\mathcal{O}(\eta)$ whose elastic and geometric properties vary over the characteristic length $L = \mathcal{O}(1)$. This sets the stage for homogenization which captures the effective properties of the lattice in the limit $\eta \to 0$ where the scales are well separated.

We introduce a general description of lattices based on the particular example of a *curved honeycomb lattice*. Other lattice geometries are illustrated in the numerical examples in Section 7. Specifically, we consider *networks of elastic beams connected by rigid hinges*. We refer to rigid hinges as *nodes*, and to the elastic beams as *elements*: the degrees of freedom are attached to nodes and the strain energy is stored in elements, as in the finite-element method.

### 2.1. Underlying topological lattice

The lattice topology is defined by an *underlying topological lattice* which is periodic by assumption, see Figure 2.1a. It is a mathematical abstraction: the real elastic lattice may be curved, hence non-periodic, as shown in Figure 2.1b. The following notions are associated with the topological lattice:

* **nodes** The nodes are labelled by an index $\beta$. Their position is denoted as $\tilde{X}_{\beta} \in \mathbb{R}^d$, where $d$ is the dimension of the space in which the lattice is embedded. All quantities relating to the topological lattice are denoted with a tilde.

* **Bravais sub-lattices** The nodes can be grouped into Bravais sub-lattices indexed by an integer $b$ with $1 \le b \le n_b$, where $n_b$ is the number of sub-lattices. The sub-lattices are defined as the subsets of nodes that are mapped to one another through the fundamental translations of the topological lattice. For instance, the honeycomb lattice in Figure 2.1a has $n_b = 2$ Bravais sub-lattices, shown using the open vs. solid disks. The index of the Bravais sub-lattice to which a particular node $\beta$ belongs is denoted as $b(\beta)$.

* **elements** Elements, which are elastic beams in the present illustration, are labelled using an index $a$. An orientation is assigned to every element $a$, and we denote by $(\beta_a^-, \beta_a^+)$ its endpoints, with $\beta_a^-$ as the tail and $\beta_a^+$ as the head, see Figure 2.1b'.

* **element families** Elements are grouped into $n_\varphi$ families indexed by $\varphi$, having similar elastic and geometric properties. The honeycomb lattice in Figure 2.1a, for instance, has $n_\varphi = 3$ element families, $\varphi \in \{I, II, III\}$, corresponding to the three possible orientations of the beams. The family to which a particular element $a$ belongs is denoted as $\varphi_a$. The density of elements of a particular family $\varphi$ per unit surface (if the space dimension is $n=2$) or volume (if $n=3$) is denoted as $\tilde{\rho}_\varphi$.

The topological lattice being periodic, the vector $\tilde{X}_{\beta_a^+} - \tilde{X}_{\beta_a^-}$ joining the endpoints of a particular element $a$ is the same for all elements $a$ in a particular family $\varphi$ and can be written as

$$ \tilde{X}_{\beta_a^+} - \tilde{X}_{\beta_a^-} = \delta_{\varphi_a} \qquad (2.1) $$

where $\delta_\varphi$ maps a family index $\varphi$ to a vector in $\mathbb{R}^d$.
The present paper addresses higher-order homogenization of periodic networks made up of slender 1D elements (typically beams or springs), further referred to as *lattice material* or simply *lattice*. The homogenization is done symbolically and in an automated way, and builds up on the general-purpose homogenization engine from our previous work [AL23]. It delivers closed-form analytical expressions of the effective elastic properties without any adjustable internal kinematic variable or parameter. By contrast with earlier work focussing on particular lattice architectures [DO11; Dav13; LR13] or 1D/2D networks of springs [MA02; MA06], our approach is designed to be general and versatile. It works in 1D, 2D or 3D and can handle arbitrary lattice topologies, the presence of prestress, slowly varying geometrical or elastic properties, a variety of member types (such as springs, beams or curved arches), as well as a large elastic contrast.

We illustrate our homogenization method on a comprehensive set of examples: a uniform 2D honeycomb, spatially varying 2D Kagome and honeycomb architectures, a 3D quasi-octet lattice and 1D pantograph. With a view of assessing the accuracy of the resulting models, we propose a numerical verification procedure that compares the microscopic fluctuations extracted from full numerical simulations to those predicted by the homogenization method. We apply this procedure to a range of structural problems including macroscopic samples with cracks, holes, as well as spatially varying pre-stress, elastic and geometric properties. Our numerical results show that the higher-order effective model very accurately captures the fine features of the local displacements and rotations away from the boundaries, where the presence of boundary layers is known to make homogenization break down.

Higher-order homogenization comes with several difficulties, and some remain beyond the scope of the present work. Most notably, we do not address the presence of boundary layers, and accept that the gradient stiffness is often negative (see Remark 5.3 below). These limitations prevent from using the homogenized energy to set up self-contained simulations so far, as discussed in further detail in the conclusion.

## 2. GENERATING THE ELASTIC LATTICE

In this section, we describe the process of generating the elastic lattices. Given a small parameter $\eta \ll 1$, we produce a lattice having cells of size $\mathcal{O}(\eta)$ whose elastic and geometric properties vary over the characteristic length $L = \mathcal{O}(1)$. This sets the stage for homogenization which captures the effective properties of the lattice in the limit $\eta \to 0$ where the scales are well separated.

We introduce a general description of lattices based on the particular example of a *curved honeycomb lattice*. Other lattice geometries are illustrated in the numerical examples in Section 7. Specifically, we consider *networks of elastic beams connected by rigid hinges*. We refer to rigid hinges as *nodes*, and to the elastic beams as *elements*: the degrees of freedom are attached to nodes and the strain energy is stored in elements, as in the finite-element method.

### 2.1. Underlying topological lattice

The lattice topology is defined by an *underlying topological lattice* which is periodic by assumption, see Figure 2.1a. It is a mathematical abstraction: the real elastic lattice may be curved, hence non-periodic, as shown in Figure 2.1b. The following notions are associated with the topological lattice:

* **nodes** The nodes are labelled by an index $\beta$. Their position is denoted as $\tilde{X}_{\beta} \in \mathbb{R}^d$, where $d$ is the dimension of the space in which the lattice is embedded. All quantities relating to the topological lattice are denoted with a tilde.

* **Bravais sub-lattices** The nodes can be grouped into Bravais sub-lattices indexed by an integer $b$ with $1 \le b \le n_b$, where $n_b$ is the number of sub-lattices. The sub-lattices are defined as the subsets of nodes that are mapped to one another through the fundamental translations of the topological lattice. For instance, the honeycomb lattice in Figure 2.1a has $n_b = 2$ Bravais sub-lattices, shown using the open vs. solid disks. The index of the Bravais sub-lattice to which a particular node $\beta$ belongs is denoted as $b(\beta)$.

* **elements** Elements, which are elastic beams in the present illustration, are labelled using an index $a$. An orientation is assigned to every element $a$, and we denote by $(\beta_a^-, \beta_a^+)$ its endpoints, with $\beta_a^-$ as the tail and $\beta_a^+$ as the head, see Figure 2.1b'.

* **element families** Elements are grouped into $n_\varphi$ families indexed by $\varphi$, having similar elastic and geometric properties. The honeycomb lattice in Figure 2.1a, for instance, has $n_\varphi = 3$ element families, $\varphi \in \{I, II, III\}$, corresponding to the three possible orientations of the beams. The family to which a particular element $a$ belongs is denoted as $\varphi_a$. The density of elements of a particular family $\varphi$ per unit surface (if the space dimension is $n=2$) or volume (if $n=3$) is denoted as $\tilde{\rho}_\varphi$.

The topological lattice being periodic, the vector $\tilde{X}_{\beta_a^+} - \tilde{X}_{\beta_a^-}$ joining the endpoints of a particular element $a$ is the same for all elements $a$ in a particular family $\varphi$ and can be written as

$$ \tilde{X}_{\beta_a^+} - \tilde{X}_{\beta_a^-} = \delta_{\varphi_a} \qquad (2.1) $$

where $\delta_\varphi$ maps a family index $\varphi$ to a vector in $\mathbb{R}^d$.
By a similar argument, the Bravais sub-lattice $b(\beta_a^\pm)$ of the endpoints $\pm$ of an element $\alpha$ depends on the element family $\varphi_\alpha$ and not on the element $\alpha$ itself: there exists a function $b_\varphi^\pm$ mapping a family $\varphi$ and an endpoint index $\pm$ (head or tail) to a Bravais index $b$, such that

$$b(\beta_a^\pm) = b_{\varphi_a^\pm}. \qquad (2.2)$$

The maps $\tilde{\delta}_\varphi$ and $b_\varphi^\pm$ for the topological honeycomb lattice are given in Table 2.1.

<figure><img src="image_1.png" /><figcaption>**Figure 2.1.** Generation of a lattice possessing slowly varying properties. (a) Periodic, underlying topological lattice; the conventional element orientation is indicated by the arrows. (b) The reference configuration of the lattice is produced by applying a contraction factor $\eta$ and a diffeomorphism $f$, see (2.3). (b') Close-up view of a particular element (beam) $\alpha$, belonging to family $\varphi_\alpha = I$; note that the conventional element center $\tilde{X}_\alpha$ is distinct from its midpoint. The figure was generated for a honeycomb lattice with the diffeomorphism in (2.4) and with parameters $\alpha = 1$, $\eta = 0.4$ and $R = 5$.</figcaption></figure>

$$\begin{array}{ccc} d & n_b & n_\varphi & \tilde{\rho}_\varphi \\ \hline 2 & 2 & 3 & \frac{2}{3\sqrt{3}a^2} \end{array}$$

$$\begin{array}{c|ccc} \varphi & I & II & III \\ \hline b_{\tilde{\varphi}} & 2 & 1 & 1 \\ b_{\varphi}^{+} & 1 & 2 & 2 \\ \tilde{\delta}_{\varphi} & a e_1 & a \left( \frac{e_1}{2} + \frac{\sqrt{3}}{2} e_2 \right) & a \left( \frac{e_1}{2} - \frac{\sqrt{3}}{2} e_2 \right) \end{array}$$

<caption>**Table 2.1.** Properties of the topological honeycomb lattice in Figure 2.1a: space dimension $d$, number $n_b$ of Bravais sub-lattices, number $n_\varphi$ of element families, density $\tilde{\rho}_\varphi$ of elements in a given family per unit area, connectivity $b_\varphi^\pm$, edge vectors $\tilde{\delta}_\varphi$.</caption>

**Remark 2.1.** By contrast with most of the existing work on periodic homogenization, our approach does not require a unit cell to be defined: the information contained in a unit cell is entirely captured by the quantities shown in Table 2.1.

## 2.2. Curved reference configuration

The reference configuration of the elastic lattice, shown in Figure 2.1b, is produced by applying a diffeomorphism $f: \mathbb{R}^d \to \mathbb{R}^d$ combined with a contraction with ratio $\eta$, to the topological lattice: the position in the reference configuration of the node $\beta$ is

$$\mathbf{X}_{\beta} = \mathbf{f}(\eta \tilde{\mathbf{X}}_{\beta}). \qquad (2.3)$$

This generation procedure allows for a variety of lattice designs for a given underlying topological lattice.

By way of illustration, the curved honeycomb lattice shown in Figure 2.1b has been produced using the 'complex-exponential' diffeomorphism

$$f(\tilde{\mathbf{x}}) = R \exp\left(\frac{\tilde{\mathbf{x}}_2}{R}\right) e_r\left(\frac{\pi}{2} - \frac{\tilde{\mathbf{x}}_1}{R}\right) - R e_2, \quad \text{(curved honeycomb lattice)} \qquad (2.4)$$

where $e_r(\theta) = \cos\theta e_1 + \sin\theta e_2$ denotes the unit radial vector and $R$ is the curvature radius of the arc of circle that is the image of the curve $\tilde{X}_2 = 0$. The case of a periodic, non-curved honeycomb lattice corresponds to $R \to \infty$, which yields $f = \text{id}_2$ (identity in $\mathbb{R}^2$).

Each element $\alpha$ is assigned a conventional center $\mathbf{X}_\alpha^c$ in reference configuration. Various definitions of the center are possible. We choose to define it as the image by the diffeomorphism of the edge's midpoint of the underlying topological lattice,

$$\mathbf{X}_\alpha^c = f\left(\frac{\eta}{2}\left(\tilde{\mathbf{X}}_{\beta_\alpha} + \tilde{\mathbf{X}}_{\beta_\alpha^+}\right)\right). \qquad (2.5)$$

A particular center $\mathbf{X}_\alpha^c$ is shown in Figure 2.1b'.
By a similar argument, the Bravais sub-lattice $b(\beta_a^\pm)$ of the endpoints $\pm$ of an element $\alpha$ depends on the element family $\varphi_\alpha$ and not on the element $\alpha$ itself: there exists a function $b_\varphi^\pm$ mapping a family $\varphi$ and an endpoint index $\pm$ (head or tail) to a Bravais index $b$, such that

$$b(\beta_a^\pm) = b_{\varphi_a^\pm}. \qquad (2.2)$$

The maps $\tilde{\delta}_\varphi$ and $b_\varphi^\pm$ for the topological honeycomb lattice are given in Table 2.1.

<figure><img src="image_3.png" /><figcaption>**Figure 2.1.** Generation of a lattice possessing slowly varying properties. (a) Periodic, underlying topological lattice; the conventional element orientation is indicated by the arrows. (b) The reference configuration of the lattice is produced by applying a contraction factor $\eta$ and a diffeomorphism $f$, see (2.3). (b') Close-up view of a particular element (beam) $\alpha$, belonging to family $\varphi_\alpha = I$; note that the conventional element center $\tilde{X}_\alpha$ is distinct from its midpoint. The figure was generated for a honeycomb lattice with the diffeomorphism in (2.4) and with parameters $\alpha = 1$, $\eta = 0.4$ and $R = 5$.</figcaption></figure>

$$\begin{array}{cccc}
d & n_b & n_\varphi & \tilde{\rho}_\varphi \\
\hline
2 & 2 & 3 & \frac{2}{3\sqrt{3}a^2}
\end{array}$$

$$\begin{array}{c|ccc}
\varphi & I & II & III \\
\hline
b_{\tilde{\varphi}} & 2 & 1 & 1 \\
b_{\varphi}^{+} & 1 & 2 & 2 \\
\tilde{\delta}_{\varphi} & a e_1 & a \left(\frac{e_1}{2} + \frac{\sqrt{3}}{2} e_2\right) & a \left(\frac{e_1}{2} - \frac{\sqrt{3}}{2} e_2\right)
\end{array}$$

<caption>**Table 2.1.** Properties of the topological honeycomb lattice in Figure 2.1a: space dimension $d$, number $n_b$ of Bravais sub-lattices, number $n_\varphi$ of element families, density $\tilde{\rho}_\varphi$ of elements in a given family per unit area, connectivity $b_\varphi^\pm$, edge vectors $\tilde{\delta}_\varphi$.</caption>

**Remark 2.1.** By contrast with most of the existing work on periodic homogenization, our approach does not require a unit cell to be defined: the information contained in a unit cell is entirely captured by the quantities shown in Table 2.1.

## 2.2. Curved reference configuration

The reference configuration of the elastic lattice, shown in Figure 2.1b, is produced by applying a diffeomorphism $f: \mathbb{R}^d \to \mathbb{R}^d$ combined with a contraction with ratio $\eta$, to the topological lattice: the position in the reference configuration of the node $\beta$ is

$$\mathbf{X}_{\beta} = f(\eta \tilde{\mathbf{X}}_{\beta}). \qquad (2.3)$$

This generation procedure allows for a variety of lattice designs for a given underlying topological lattice.

By way of illustration, the curved honeycomb lattice shown in Figure 2.1b has been produced using the 'complex-exponential' diffeomorphism

$$f(\tilde{\mathbf{x}}) = R \exp\left(\frac{\tilde{\mathbf{x}}_2}{R}\right) e_r\left(\frac{\pi}{2} - \frac{\tilde{\mathbf{x}}_1}{R}\right) - R e_2, \quad \text{(curved honeycomb lattice)} \qquad (2.4)$$

where $e_r(\theta) = \cos\theta e_1 + \sin\theta e_2$ denotes the unit radial vector and $R$ is the curvature radius of the arc of circle that is the image of the curve $\tilde{X}_2 = 0$. The case of a periodic, non-curved honeycomb lattice corresponds to $R \to \infty$, which yields $f = \text{id}_2$ (identity in $\mathbb{R}^2$).

Each element $\alpha$ is assigned a conventional center $\mathbf{X}_\alpha^c$ in reference configuration. Various definitions of the center are possible. We choose to define it as the image by the diffeomorphism of the edge's midpoint of the underlying topological lattice,

$$\mathbf{X}_\alpha^c = f\left(\frac{\eta}{2}\left(\tilde{\mathbf{X}}_{\beta_\alpha} + \tilde{\mathbf{X}}_{\beta_\alpha^+}\right)\right). \qquad (2.5)$$

A particular center $\mathbf{X}_\alpha^c$ is shown in Figure 2.1b'.
For any particular element, we now consider the vector $\mathbf{X}_{\beta_{\alpha}^{\pm}} - \mathbf{X}_{\alpha}^{c}$ joining (in reference configuration) the conventional center $\mathbf{X}_{\alpha}^{c}$ to either one of the endpoints $\mathbf{X}_{\beta_{\alpha}^{\pm}}$. We claim that this vector can be written as

$$ \mathbf{X}_{\beta_{\alpha}^{\pm}} - \mathbf{X}_{\alpha}^{c} = \Delta_{\phi}^{\pm}(\eta, \mathbf{X}_{\alpha}^{c}), \quad \text{with} \quad \Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = \mathcal{O}(\eta). \qquad (2.6) $$

Here, the $\Delta_{\phi}^{\pm}(\eta, \mathbf{X})$'s are smooth functions of their continuous arguments $(\eta, \mathbf{X})$ that further depend on the discrete variables in subscript (element family $\phi$) and superscript ($\pm$ labelling endpoints). The dependence on the element $\alpha$ in the left-hand side of (2.6) has been turned in the right-hand side into a dependence on the element family $\phi_{\alpha}$ and on the center position $\mathbf{X}_{\alpha}^{c}$. Equation (2.6) captures the fact that $\mathbf{X}_{\beta_{\alpha}^{\pm}} - \mathbf{X}_{\alpha}^{c}$ varies slowly across neighboring elements belonging to the same family.

A constructive proof of (2.6) is as follows. Observe first that (2.1) and (2.5) can be combined into $\tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}} = \frac{1}{2}(\tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}} + \tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}}^{\prime}) \pm \frac{1}{2}(\tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}} - \tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}}^{\prime}) = \frac{1}{\eta}f^{-1}(\mathbf{X}_{\alpha}^{c}) \pm \frac{1}{2}\tilde{\delta}_{\phi_{\alpha}}$, so that $\mathbf{X}_{\beta_{\alpha}^{\pm}} - \mathbf{X}_{\alpha}^{c} = f(\eta(\tilde{\mathbf{X}}_{\beta_{\alpha}^{\pm}}^{\prime}) - \mathbf{X}_{\alpha}^{c}) = f(f^{-1}(\mathbf{X}_{\alpha}^{c}) \pm \frac{\eta}{2}\tilde{\delta}_{\phi_{\alpha}}) - \mathbf{X}_{\alpha}^{c}$. Identifying with (2.6), we find

$$ \Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = f(f^{-1}(\mathbf{X}) \pm \frac{\eta}{2}\tilde{\delta}_{\phi}) - \mathbf{X}. \qquad (2.7) $$

The right-hand side vanishes when $\eta \to 0$ and the smoothness of $f$ proves the estimate $\Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = \mathcal{O}(\eta)$ stated in (2.6).

By way of illustration, consider elements of type $\phi = I$ in a curved honeycomb, see Figure 2.1(b). Combining $f^{-1}(\mathbf{X}) = R(\tan^{-1}(\frac{X_1}{R+X_2})e_1 + \frac{1}{R}\ln\sqrt{X_1^2+(R+X_2)^2}e_2)$ from (2.4) and $\tilde{\delta}_I = ae_1$, we have from (2.7)

$$ \Delta_I^{\pm}(\eta, \mathbf{X}) = \sqrt{X_1^2 + (R+X_2)^2} e_r \left(\frac{\pi}{2} - \tan^{-1}\left(\frac{X_1}{R+X_2}\right) \mp \eta \frac{a}{2R}\right) - (X_1 e_1 + (R+X_2) e_2) \quad \text{(curved honeycomb lattice).} \qquad (2.8) $$

The two other functions $\Delta_{II}^{\pm}(\eta, \mathbf{X})$ and $\Delta_{III}^{\pm}(\eta, \mathbf{X})$ are obtained in a similar way.

## 2.3. Assumed power-series dependence on $\eta$

Expanding Equation (2.7) in powers of $\eta$, we get

$$ \Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = \overbrace{f(f^{-1}(\mathbf{X})) - \mathbf{X}}^{0} \pm \frac{\eta}{2} \nabla f(f^{-1}(\mathbf{X})) \cdot \tilde{\delta}_{\phi} + \frac{\eta^2}{4} \nabla^2 f(f^{-1}(\mathbf{X})) : (\tilde{\delta}_{\phi} \otimes \tilde{\delta}_{\phi}) + \dots \qquad (2.9) $$

where the cancellation at leading order leads to the estimate (2.6). The diffeomorphism $f$ being infinitely smooth by assumption, this expansion could be pushed to any other.

More generally, we will assume that any function $g(\eta, \mathbf{X})$ depending on the expansion parameter $\eta$ can be expanded as a power series in $\eta$: when we write $g(\eta, \mathbf{X})$ with $g = \mathcal{O}(\eta^p)$, we imply that $g$ is of the form

$$ g(\eta, \mathbf{X}) = \eta^p (g_{(p)}(\mathbf{X}) + \eta g_{(p+1)}(\mathbf{X}) + \dots), \qquad (2.10) $$

where the $g_{(k)}$'s are smooth functions of $\mathbf{X}$ that do not depend on $\eta$. Expansions of the form (2.10) are known as regular. For any lattice property, regular expansions can be obtained in explicit form as we just did for $g = \Delta_{\phi}^{\pm}$ in (2.9). For all the elastic fields, the regular expansions will be assumed, see Equation (3.2) below, as we approach homogenization based on formal expansions. At boundary layers, the regular expansions break down and the elastic fields are given by singular expansions instead: their analysis is beyond the scope of this paper.

## 2.4. Kinematic analysis of the beam-elements

For the moment, we limit attention to naturally straight, 2D beam elements: this applies to the lattice shown in Figure 2.1b. The extension to elastic springs is discussed later in Section 7. The extension to 3D beams or arches in both 2D and 3D (i.e., of beams possessing natural curvature) is straightforward and is left for future work.

The lattice is deformed by applying forces. The displacement and rotation of the node $\mathbf{X}_{\beta}$ are denoted as $v_{\beta} \in \mathbb{R}^d$, and $\theta_{\beta} \in \mathbb{R}^{n_r}$, respectively (in 2D, the infinitesimal rotation $\theta_{\beta} = (\theta_{\beta})$ is a scalar, $n_r = 1$; see Section A.2 in the Appendix for the extension to an arbitrary dimension $d$). The assumption of rigid nodes implies that the infinitesimal rotation $\theta_{\beta}$ is common to the adjacent elements at any particular node $\mathbf{X}_{\beta}$. The nodal degrees of freedom are therefore

$$ (v_{\beta}, \theta_{\beta}) \in \mathbb{R}^{n_n}, \quad \text{where } n_n = d + n_r = \begin{cases} 1 & \text{if } d=1 \\ 3 & \text{if } d=2. \\ 6 & \text{if } d=3 \end{cases}. \qquad (2.11) $$

For the 2D elements discussed here, $d=2$, $n_r = 1$ and $n_n = 3$.

Each element family $\phi$ is characterized by a set of $n_{E_{\phi}}$ element strain measures, which are collected into an element strain vector $\mathbf{E}_{\alpha} \in \mathbb{R}^{n_{E_{\phi}}}$ in 2D, beam elements uses $n_{E_{\phi}} = 3$ strain measures, namely the extensional strain $\epsilon_{\alpha}$, the bending strain $\kappa_{\alpha}$ and the shear strain $\tau_{\alpha}$, which are defined by

$$ \mathbf{E}_{\alpha} = \left( \overbrace{(\mathbf{v}_{\beta_{\alpha}^{\pm} - \mathbf{v}_{\beta_{\alpha}^{\pm}}) \cdot \mathbf{t}_{\alpha}}^{\kappa_{\alpha}}, \quad \overbrace{\theta_{\beta_{\alpha}^{\pm} - \theta_{\beta_{\alpha}^{\pm}}}^{\tau_{\alpha}}, \quad \overbrace{\frac{\theta_{\beta_{\alpha}^{\pm} + \theta_{\beta_{\alpha}^{\pm}}}{2} - \frac{(\mathbf{v}_{\beta_{\alpha}^{\pm} - \mathbf{v}_{\beta_{\alpha}^{\pm}}) \cdot \mathbf{n}_{\alpha}}{\ell_{\alpha}}}^{\tau_{\alpha}}} \right) \in \mathbb{R}^3, \quad n_{E_{\phi}} = 3, \quad (2\text{D beam}) \qquad (2.12) $$
For any particular element, we now consider the vector $X_{\beta_a}^c - X_{\alpha_a}^c$ joining (in reference configuration) the conventional center $X_{\alpha_a}^c$ to either one of the endpoints $X_{\beta_a}^c$. We claim that this vector can be written as

$$X_{\beta_a}^c - X_{\alpha_a}^c = \Delta_{\phi_a}^{\pm}(\eta, X_{\alpha_a}^c), \quad \text{with } \Delta_{\phi_a}^{\pm}(\eta, X) = \mathcal{O}(\eta). \qquad (2.6)$$

Here, the $\Delta_{\phi_a}^{\pm}(\eta, X)$'s are smooth functions of their continuous arguments $(\eta, X)$ that further depend on the discrete variables in subscript (element family $\phi$) and superscript ($\pm$ labelling endpoints). The dependence on the element $\alpha$ in the left-hand side of (2.6) has been turned in the right-hand side into a dependence on the element family $\phi_a$ and on the center position $X_{\alpha_a}^c$: Equation (2.6) captures the fact that $X_{\beta_a}^c - X_{\alpha_a}^c$ varies slowly across neighboring elements belonging to the same family.

A constructive proof of (2.6) is as follows. Observe first that (2.1) and (2.5) can be combined into $\tilde{X}_{\beta_a}^c = \frac{1}{2}(\tilde{X}_{\beta_a}^c + \tilde{X}_{\beta_a^+}^c) \pm \frac{1}{2}(\tilde{X}_{\beta_a^+} - \tilde{X}_{\beta_a}) = \frac{1}{\eta}f^{-1}(X_{\alpha_a}^c) \pm \frac{1}{2}\tilde{\delta}_{\phi_a}$, so that $X_{\beta_a^+} - X_{\alpha_a}^c = f(\eta\tilde{X}_{\beta_a^+}^c) - X_{\alpha_a}^c = f(f^{-1}(X_{\alpha_a}^c) \pm \frac{\eta}{2}\tilde{\delta}_{\phi_a}) - X_{\alpha_a}^c$. Identifying with (2.6), we find

$$\Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = f(f^{-1}(\mathbf{X}) \pm \frac{\eta}{2}\tilde{\delta}_{\phi}) - \mathbf{X}. \qquad (2.7)$$

The right-hand side vanishes when $\eta \to 0$ and the smoothness of $f$ proves the estimate $\Delta_{\phi}^{\pm}(\eta, X) = \mathcal{O}(\eta)$ stated in (2.6).

By way of illustration, consider elements of type $\phi = I$ in a curved honeycomb, see Figure 2.1(b). Combining $f^{-1}(X) = R(\tan^{-1}(\frac{X_1}{R+X_2})e_1 + \frac{1}{R}\ln\sqrt{X_1^2+(R+X_2)^2}e_2)$ from (2.4) and $\tilde{\delta}_I = ae_1$, we have from (2.7)

$$\Delta_I^{\pm}(\eta, \mathbf{X}) = \sqrt{X_1^2 + (R+X_2)^2} e_r \left(\frac{\pi}{2} - \tan^{-1}\left(\frac{X_1}{R+X_2}\right) \mp \eta \frac{a}{2R}\right) - (X_1 e_1 + (R+X_2) e_2) \quad \text{(curved honeycomb lattice).} \qquad (2.8)$$

The two other functions $\Delta_{II}^{\pm}(\eta, X)$ and $\Delta_{III}^{\pm}(\eta, X)$ are obtained in a similar way.

## 2.3. Assumed power-series dependence on $\eta$

Expanding Equation (2.7) in powers of $\eta$, we get

$$\Delta_{\phi}^{\pm}(\eta, \mathbf{X}) = \overbrace{f(f^{-1}(\mathbf{X})) - \mathbf{X}}^{0} \pm \frac{\eta}{2} \nabla f(f^{-1}(\mathbf{X})) \cdot \tilde{\delta}_{\phi} + \frac{\eta^2}{4} \nabla^2 f(f^{-1}(\mathbf{X})) : (\tilde{\delta}_{\phi} \otimes \tilde{\delta}_{\phi}) + \dots \qquad (2.9)$$

where the cancellation at leading order leads to the estimate (2.6). The diffeomorphism $f$ being infinitely smooth by assumption, this expansion could be pushed to any other.

More generally, we will assume that any function $g(\eta, X)$ depending on the expansion parameter $\eta$ can be expanded as a power series in $\eta$: when we write $g(\eta, X)$ with $g = \mathcal{O}(\eta^p)$, we imply that $g$ is of the form

$$g(\eta, \mathbf{X}) = \eta^p (g_{(p)}(\mathbf{X}) + \eta g_{(p+1)}(\mathbf{X}) + \dots), \qquad (2.10)$$

where the $g_{(k)}$'s are smooth functions of $X$ that do not depend on $\eta$. Expansions of the form (2.10) are known as regular. For any lattice property, regular expansions can be obtained in explicit form as we just did for $g = \Delta_{\phi}^{\pm}$ in (2.9). For all the elastic fields, the regular expansions will be assumed, see Equation (3.2) below, as we approach homogenization based on formal expansions. At boundary layers, the regular expansions break down and the elastic fields are given by singular expansions instead: their analysis is beyond the scope of this paper.

## 2.4. Kinematic analysis of the beam-elements

For the moment, we limit attention to naturally straight, 2D beam elements: this applies to the lattice shown in Figure 2.1b. The extension to elastic springs is discussed later in Section 7. The extension to 3D beams or arches in both 2D and 3D (i.e., of beams possessing natural curvature) is straightforward and is left for future work.

The lattice is deformed by applying forces. The displacement and rotation of the node $X_{\beta}$ are denoted as $v_{\beta} \in \mathbb{R}^d$, and $\theta_{\beta} \in \mathbb{R}^{n_r}$, respectively (in 2D, the infinitesimal rotation $\theta_{\beta} = (\theta_{\beta})$ is a scalar, $n_r = 1$; see Section A.2 in the Appendix for the extension to an arbitrary dimension $d$). The assumption of rigid nodes implies that the infinitesimal rotation $\theta_{\beta}$ is common to the adjacent elements at any particular node $X_{\beta}$. The nodal degrees of freedom are therefore

$$(v_{\beta}, \theta_{\beta}) \in \mathbb{R}^{n_n}, \quad \text{where } n_n = d + n_r = \begin{cases} 1 & \text{if } d=1 \\ 3 & \text{if } d=2. \\ 6 & \text{if } d=3 \end{cases}. \qquad (2.11)$$

For the 2D elements discussed here, $d=2$, $n_r = 1$ and $n_n = 3$.

Each element family $\phi$ is characterized by a set of $n_{E_{\phi}}$ element strain measures, which are collected into an element strain vector $E_{\alpha} \in \mathbb{R}^{n_{E_{\phi}}}$, in 2D, beam elements uses $n_{E_{\phi}} = 3$ strain measures, namely the extensional strain $\epsilon_{\alpha}$, the bending strain $\kappa_{\alpha}$ and the shear strain $\tau_{\alpha}$, which are defined by

$$E_{\alpha} = \left( \overbrace{\left( \frac{\mathbf{v}_{\beta_a} + \mathbf{v}_{\beta_a^+} \cdot \mathbf{t}_a}{\ell_a} \right)}^{\kappa_a}, \quad \overbrace{\theta_{\beta_a^+} - \theta_{\beta_a}}^{\kappa_a}, \quad \overbrace{\frac{\theta_{\beta_a}^+ + \theta_{\beta_a^+} - \frac{(\mathbf{v}_{\beta_a}^+ - \mathbf{v}_{\beta_a^+}) \cdot \mathbf{n}_a}{\ell_a}}^{\tau_a}}^{\kappa_a} \right) \in \mathbb{R}^3, \quad n_{E_{\phi}} = 3, \quad (2\text{D beam}) \qquad (2.12)$$
where $\ell_a = |\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-}|$ is the undeformed length, $\mathbf{t}_a = (\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-}) / \ell_a$ is the unit tangent, $\mathbf{n}_a = \mathbf{e}_3 \times \mathbf{t}_a$ is the undeformed unit normal and $\mathbf{e}_3 = \mathbf{e}_1 \times \mathbf{e}_2$ is the unit normal to the Euclidean plane, see Figure 2.2. The quantities $(\epsilon_a, \kappa_a, \tau_a)$ are effective measures characterizing the strain in a beam-element globally, which do not depend on the arc-length parameter. As shown in Figure 2.2(b), for instance, the shearing mode $\tau_a$ defined at the element scale yields ‘locally’ non-uniform curvature as a function of arc-length—this local response will captured by the element stiffness derived in the forthcoming Section 2.5.

<figure><img src="image_2.png" /><figcaption>**Figure 2.2.** A 2D beam element, having index $\alpha$ and conventional center $\mathbf{X}_\alpha^c$: (a) geometry in reference configuration and (b) the $n_{E_\phi} = 3$ deformation modes: extension, bending and shearing.</figcaption></figure>

We observe that the undeformed element's length $\ell_a$, unit tangent $\mathbf{t}_a$ and unit normal $\mathbf{n}_a$ can each be expressed as functions of the element family $\varphi_\alpha$, of the scale separation parameter $\eta$ and of the element center $\mathbf{X}_\alpha^c$. Indeed, with help of (2.6), we have $\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-} = \Delta_{\varphi_\alpha}^+(\eta, \mathbf{X}_\alpha^c) - \Delta_{\varphi_\alpha}^-(\eta, \mathbf{X}_\alpha^c)$, which yields

$$ 
\begin{align} 
\ell_a &= \hat{\ell}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\ell}_{\varphi}(\eta, \mathbf{X}) = |\Delta_{\varphi}^+(\eta, \mathbf{X}) - \Delta_{\varphi}^-(\eta, \mathbf{X})| = \mathcal{O}(\eta^1), \nonumber \\ 
\mathbf{t}_a &= \hat{\mathbf{t}}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\mathbf{t}}_{\varphi}(\eta, \mathbf{X}) = \frac{\Delta_{\varphi}^+(\eta, \mathbf{X}) - \Delta_{\varphi}^-(\eta, \mathbf{X})}{\hat{\ell}_{\varphi}(\eta, \mathbf{X})} = \mathcal{O}(\eta^0) \tag{2.13} \\ 
\mathbf{n}_a &= \hat{\mathbf{n}}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\mathbf{n}}_{\varphi}(\eta, \mathbf{X}) = \mathbf{e}_3 \times \hat{\mathbf{t}}_{\varphi}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0). \nonumber 
\end{align} 
$$

The dependence of these function on $\mathbf{X}$ applies to lattices having graded geometrical properties (see §7.3). The functions $\hat{\ell}_{\varphi_\alpha}$, $\hat{\mathbf{t}}_{\varphi_\alpha}$ and $\hat{\mathbf{n}}_{\varphi_\alpha}$ have regular expansions in $\eta$ of the form (2.10) as can be checked using (2.9).

In view of this, the element strain $\mathbf{E}_\alpha$ in (2.12) can be rewritten in block-matrix notation as

$$ \mathbf{E}_\alpha = \mathcal{D}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) : \begin{pmatrix} \mathbf{v}_{\beta_\alpha^-} \\ \theta_{\beta_\alpha^-} \end{pmatrix} \begin{pmatrix} \mathbf{v}_{\beta_\alpha^+} \\ \theta_{\beta_\alpha^+} \end{pmatrix}, \quad (2.14) $$

where: denotes the double contraction, see Appendix A, $(\mathbf{v}_\beta, \theta_\beta) \in \mathbb{R}^{n_\alpha}$ are the nodal degrees of freedom, see (2.11), and $\mathcal{D}_\varphi \in \mathbb{T}^{(n_{E_\varphi}, 2, n_\alpha)}$ is the displacement-to-strain tensor,

$$ \mathcal{D}_\varphi(\eta, \mathbf{X}) = \begin{pmatrix} \left( \frac{-\hat{\mathbf{t}}_\varphi(\eta, \mathbf{X})}{\hat{\ell}_\varphi(\eta, \mathbf{X})} \begin{pmatrix} 0 \\ \mathbf{v}_2 \end{pmatrix} \right) & \left( \frac{+\hat{\mathbf{n}}_\varphi(\eta, \mathbf{X})}{\hat{\ell}_\varphi(\eta, \mathbf{X})} \begin{pmatrix} 1/2 \\ \mathbf{v}_2 \end{pmatrix} \right) \\ \left( \frac{+\hat{\mathbf{t}}_\varphi(\eta, \mathbf{X})}{\hat{\ell}_\varphi(\eta, \mathbf{X})} \begin{pmatrix} 0 \\ \mathbf{v}_2 \end{pmatrix} \right) & \left( \frac{-\hat{\mathbf{n}}_\varphi(\eta, \mathbf{X})}{\hat{\ell}_\varphi(\eta, \mathbf{X})} \begin{pmatrix} 1/2 \\ \mathbf{v}_2 \end{pmatrix} \right) \end{pmatrix} \in \mathbb{T}^{(3,2,3)} \quad (2D \text{ beam}). \quad (2.15) $$

The function $\mathcal{D}_\varphi(\eta, \mathbf{X})$ can be expressed as a regular expansion in $\eta$ of the form (2.10), as can be checked using (2.9) and (2.13).

To keep the presentation general, we will refrain from using the specific expression (2.15) of $\mathcal{D}_\varphi$ applicable to 2D beams, and use the generic form (2.14) of the displacement-strain relation instead.

## 2.5. Strain energy of the lattice

A beam $\alpha$ made up of a Hookean material and having a circular cross-section is characterized by a stretching modulus $(EA)_\alpha$ and a bending modulus $(EI)_\alpha$. The linear beam theory yields the strain energy $w_\alpha$ of the beam as $w_\alpha = \frac{1}{2}((EA)_\alpha \ell_\alpha \epsilon_\alpha^2 + \frac{(EI)_\alpha}{\ell_\alpha}(\kappa_\alpha^2 + 12\tau_\alpha^2))$, which can be rewritten as

$$ w_\alpha = \frac{(EA)_\alpha}{2} \left( \epsilon_\alpha^2 + \chi_\alpha (\kappa_\alpha^2 + 12 \tau_\alpha^2) \right) \quad (2.16) $$

(no implicit sum over $\alpha$) in terms of the beam's aspect-ratio parameter

$$ \chi_\alpha = \frac{(EI)_\alpha}{\ell_\alpha^2 (EA)_\alpha}. \quad (2.17) $$

For the slender beam theory to be applicable, the beam's radius $r_\alpha$ must be much smaller than its length, $r_\alpha \ll \ell_\alpha$, implying that the aspect-ratio parameter $\chi_\alpha \sim (r_\alpha / \ell_\alpha)^2 \ll 1$ is small.

To handle the case of lattices whose elastic properties vary slowly in space, we assume that the beam's elastic constants $(EA)_\alpha$ and $\chi_\alpha$ are prescribed as smooth functions $(\tilde{E}A)_\varphi(\eta, \mathbf{X})$ and $\tilde{\chi}_\varphi(\eta, \mathbf{X})$ of the scale separation parameter $\eta$ and of the element center $\mathbf{X}_\alpha^c$, with an additional dependence on the discrete family index $\varphi$, i.e.,

$$ (EA)_\alpha = \tilde{E}A_{(\varphi_\alpha)}(\eta, \mathbf{X}_\alpha^c), \quad \chi_\alpha = \tilde{\chi}_{(\varphi_\alpha)}(\eta, \mathbf{X}_\alpha^c). \quad (2.18) $$

In view of (2.16), we define the element stiffness matrix $\mathcal{H}_\varphi \in \mathbb{T}^{(n_{E_\varphi}, n_{E_\varphi})}$ as

$$ \mathcal{H}_\varphi(\eta, \mathbf{X}) = \tilde{E}A_\varphi(\eta, \mathbf{X}) \hat{\ell}_\varphi(\eta, \mathbf{X}) \operatorname{diag}(1, \tilde{\chi}_\varphi(\eta, \mathbf{X}), 12 \tilde{\chi}_\varphi(\eta, \mathbf{X})) \in \mathbb{T}^{(3,3)} \quad (2D \text{ beam}), \quad (2.19) $$
where $\ell_a = |\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-}|$ is the undeformed length, $\mathbf{t}_a = (\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-}) / \ell_a$ is the unit tangent, $\mathbf{n}_a = \mathbf{e}_3 \times \mathbf{t}_a$ is the undeformed unit normal and $\mathbf{e}_3 = \mathbf{e}_1 \times \mathbf{e}_2$ is the unit normal to the Euclidean plane, see Figure 2.2. The quantities $(\varepsilon_a, \kappa_a, \tau_a)$ are effective measures characterizing the strain in a beam-element globally, which do not depend on the arc-length parameter. As shown in Figure 2.2(b), for instance, the shearing mode $\tau_a$ defined at the element scale yields ‘locally’ non-uniform curvature as a function of arc-length—this local response will captured by the element stiffness derived in the forthcoming Section 2.5.

<figure><img src="image_4.png" /><figcaption>**Figure 2.2.** A 2D beam element, having index $\alpha$ and conventional center $\mathbf{X}_\alpha^c$: (a) geometry in reference configuration and (b) the $n_{E_\varphi} = 3$ deformation modes: extension, bending and shearing.</figcaption></figure>

We observe that the undeformed element's length $\ell_a$, unit tangent $\mathbf{t}_a$ and unit normal $\mathbf{n}_a$ can each be expressed as functions of the element family $\varphi_\alpha$, of the scale separation parameter $\eta$ and of the element center $\mathbf{X}_\alpha^c$. Indeed, with help of (2.6), we have $\mathbf{X}_{\beta_a^+} - \mathbf{X}_{\beta_a^-} = \Delta_{\varphi_\alpha}^+(\eta, \mathbf{X}_\alpha^c) - \Delta_{\varphi_\alpha}^-(\eta, \mathbf{X}_\alpha^c)$, which yields

$$ 
\begin{align} 
\ell_a &= \hat{\ell}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\ell}_{\varphi}(\eta, \mathbf{X}) = |\Delta_{\varphi}^+(\eta, \mathbf{X}) - \Delta_{\varphi}^-(\eta, \mathbf{X})| = \mathcal{O}(\eta^1), \nonumber \\ 
\mathbf{t}_a &= \hat{\mathbf{t}}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\mathbf{t}}_{\varphi}(\eta, \mathbf{X}) = \frac{\Delta_{\varphi}^+(\eta, \mathbf{X}) - \Delta_{\varphi}^-(\eta, \mathbf{X})}{\hat{\ell}_{\varphi}(\eta, \mathbf{X})} = \mathcal{O}(\eta^0) \tag{2.13} \\ 
\mathbf{n}_a &= \hat{\mathbf{n}}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) && \text{where } \hat{\mathbf{n}}_{\varphi}(\eta, \mathbf{X}) = \mathbf{e}_3 \times \hat{\mathbf{t}}_{\varphi}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0). \nonumber 
\end{align} 
$$

The dependence of these function on $\mathbf{X}$ applies to lattices having graded geometrical properties (see §7.3). The functions $\hat{\ell}_{\varphi_\alpha}$, $\hat{\mathbf{t}}_{\varphi_\alpha}$ and $\hat{\mathbf{n}}_{\varphi_\alpha}$ have regular expansions in $\eta$ of the form (2.10) as can be checked using (2.9).

In view of this, the element strain $\mathbf{E}_\alpha$ in (2.12) can be rewritten in block-matrix notation as

$$ \mathbf{E}_\alpha = \mathcal{D}_{\varphi_\alpha}(\eta, \mathbf{X}_\alpha^c) : \begin{pmatrix} \mathbf{v}_{\beta_\alpha^-} \\ \theta_{\beta_\alpha^-} \end{pmatrix} \begin{pmatrix} \mathbf{v}_{\beta_\alpha^+} \\ \theta_{\beta_\alpha^+} \end{pmatrix}, \quad (2.14) $$

where: denotes the double contraction, see Appendix A, $(\mathbf{v}_\beta, \theta_\beta) \in \mathbb{R}^{n_\alpha}$ are the nodal degrees of freedom, see (2.11), and $\mathcal{D}_\varphi \in \mathbb{T}^{(n_{E_\varphi}, 2, n_\alpha)}$ is the displacement-to-strain tensor,

$$ \mathcal{D}_{\varphi}(\eta, \mathbf{X}) = \begin{pmatrix} \left( \frac{-\hat{\mathbf{t}}_{\varphi}(\eta, \mathbf{X})}{\hat{\ell}_{\varphi}(\eta, \mathbf{X})} \quad 0 \right) & \left( \frac{\mathbf{v}_{\beta_\alpha^-}}{\theta_{\beta_\alpha^-}} \quad \frac{\mathbf{v}_{\beta_\alpha^+}}{\theta_{\beta_\alpha^+}} \right) \\ \left( \frac{\mathbf{v}_{\beta_\alpha^-}}{\theta_{\beta_\alpha^-}} \quad -1 \right) & \left( \frac{-\hat{\mathbf{n}}_{\varphi}(\eta, \mathbf{X})}{\hat{\ell}_{\varphi}(\eta, \mathbf{X})} \quad \frac{\mathbf{v}_{\beta_\alpha^-}}{\theta_{\beta_\alpha^-}} \right) & \frac{\mathbf{v}_{\beta_\alpha^+}}{\theta_{\beta_\alpha^+}} \quad \frac{1}{2} \end{pmatrix} \in \mathbb{T}^{(3,2,3)} \quad (2\text{D beam}). \quad (2.15) $$

The function $\mathcal{D}_\varphi(\eta, \mathbf{X})$ can be expressed as a regular expansion in $\eta$ of the form (2.10), as can be checked using (2.9) and (2.13).

To keep the presentation general, we will refrain from using the specific expression (2.15) of $\mathcal{D}_\varphi$ applicable to 2D beams, and use the generic form (2.14) of the displacement-strain relation instead.

## 2.5. Strain energy of the lattice

A beam $\alpha$ made up of a Hookean material and having a circular cross-section is characterized by a stretching modulus $(EA)_\alpha$ and a bending modulus $(EI)_\alpha$. The linear beam theory yields the strain energy $w_\alpha$ of the beam as $w_\alpha = \frac{1}{2}((EA)_\alpha \ell_\alpha \varepsilon_\alpha^2 + \frac{(EI)_\alpha}{\ell_\alpha}(\kappa_\alpha^2 + 12\tau_\alpha^2))$, which can be rewritten as

$$ w_\alpha = \frac{(EA)_\alpha}{2} \left( \varepsilon_\alpha^2 + \chi_\alpha (\kappa_\alpha^2 + 12\tau_\alpha^2) \right) \quad (2.16) $$

(no implicit sum over $\alpha$) in terms of the beam's aspect-ratio parameter

$$ \chi_\alpha = \frac{(EI)_\alpha}{\ell_\alpha^2 (EA)_\alpha}. \quad (2.17) $$

For the slender beam theory to be applicable, the beam's radius $r_\alpha$ must be much smaller than its length, $r_\alpha \ll \ell_\alpha$, implying that the aspect-ratio parameter $\chi_\alpha \sim (r_\alpha / \ell_\alpha)^2 \ll 1$ is small.

To handle the case of lattices whose elastic properties vary slowly in space, we assume that the beam's elastic constants $(EA)_\alpha$ and $\chi_\alpha$ are prescribed as smooth functions $(\tilde{EA})_\varphi(\eta, \mathbf{X})$ and $\tilde{\chi}_\varphi(\eta, \mathbf{X})$ of the scale separation parameter $\eta$ and of the element center $\mathbf{X}_\alpha^c$, with an additional dependence on the discrete family index $\varphi$, i.e.,

$$ (EA)_\alpha = \tilde{EA}_{(\varphi_\alpha)}(\eta, \mathbf{X}_\alpha^c), \quad \chi_\alpha = \tilde{\chi}_{(\varphi_\alpha)}(\eta, \mathbf{X}_\alpha^c). \quad (2.18) $$

In view of (2.16), we define the element stiffness matrix $\mathcal{H}_\varphi \in \mathbb{T}^{(n_{E_\varphi}, n_{E_\varphi})}$ as

$$ \mathcal{H}_\varphi(\eta, \mathbf{X}) = \tilde{EA}_\varphi(\eta, \mathbf{X}) \hat{\ell}_\varphi(\eta, \mathbf{X}) \operatorname{diag}(1, \tilde{\chi}_\varphi(\eta, \mathbf{X}), 12 \tilde{\chi}_\varphi(\eta, \mathbf{X})) \in \mathbb{T}^{(3,3)} \quad (2\text{D beam}), \quad (2.19) $$
so that the strain energy of an element takes the form

$$w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2} \mathbf{E}_{\alpha} \cdot \mathcal{H}_{\varphi_{\alpha}}(\eta, \mathbf{X}_{\alpha}^{c}) \cdot \mathbf{E}_{\alpha} \quad (2.20)$$

(no implicit sum over $\alpha$). To keep the presentation general, we will use this expression (2.20) of the strain energy in the following, as it is applies to linearly elastic element of various kinds such as springs, beams or arches in both 2D and 3D. The element stiffness matrix $\mathcal{H}_{\varphi} \in T^{(n_{E\varphi}, n_{E\bar{\varphi}})}$ is always symmetric but not necessarily diagonal.

The element stiffness is assumed to obey the scaling assumption

$$\mathcal{H}_{\varphi} = \mathcal{O}(\eta^d), \quad (2.21)$$

where $d$ is the space dimension. This assumption is a matter of convention: applying a global factor to all elements' stiffness scales the homogenized energy by the same factor. The convention (2.21) warrants that the energy is $\Phi = \mathcal{O}(\eta^0)$, which is convenient, see Equation (3.42) below. Since $\hat{\ell}_{\varphi} = \mathcal{O}(\eta)$ by (2.13), the stretching modulus $\tilde{E}\tilde{A}_{\varphi}$ in (2.19) will be viewed as a quantity of order

$$EA = \mathcal{O}(\eta^{d-1}) \quad (2.22)$$

in stretching-dominated lattices.

**Remark 2.2.** As noted below Equation (2.17), the function $\hat{\chi}_{\varphi}(\eta, \mathbf{X})$ is small, $\hat{\chi}_{\varphi}(\eta, \mathbf{X}) \ll 1$. Mathematically, $\hat{\chi}_{\varphi}$ will be defined as a power of the aspect-ratio $\eta$. Different power laws are considered in the different examples, such as $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^2)$ for the honeycomb, see (5.7), or $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^0)$ for the pantograph, see (7.28). The only assumption we make at this stage is that $\hat{\chi}_{\varphi}$ is a function of $\eta$ with $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^0)$—this includes $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^1)$ and $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^2)$ as particular cases. In any case, the matrix $\mathcal{H}_{\varphi}(\eta, \mathbf{X})$ appearing in (2.19) may contain terms of different orders in $\eta$ and may therefore depend on $\eta$ through more than just an overall scaling factor.

Anticipating that the homogenization will fail near the edges of the lattice due to the presence of boundary layers, we limit attention to a sub-region $\Omega \subset \mathbb{R}^d$ of the physical lattice that exclude the edges, see Figure 2.3: the strain energy of the lattice $\Phi_d$ in the region $\Omega$ is defined as the sum of the elastic energy $w_\alpha$ over all beams $\alpha$ whose center $\mathbf{X}_\alpha^c$ (shown as a red dot in Figure 2.3) lies inside the domain $\Omega$,

$$\Phi_d = \sum_{\alpha \text{ such that } \mathbf{X}_\alpha^c \in \Omega} w_\alpha(\mathbf{E}_\alpha). \quad (2.23)$$

Points where the loading is not smooth (such as points of application of point-like forces) are also excluded from the domain $\Omega$, along with a finite neighborhood around such points. The subscript ‘d’ in $\Phi_d$ stands for ‘discrete’: a continuous variant $\Phi$ of the lattice strain energy, bearing no subscript, will be introduced later.

<figure><img src="image_3.png" /><figcaption>**Figure 2.3.** To avoid boundary layers, the homogenization is carried out in a sub-region $\Omega$ of the lattice (shaded blue region) that excludes the physical boundaries of the lattice.</figcaption></figure>

**Remark 2.3.** In Sections 7.1, we treat the case of an elastic truss whose bending modulus $EI(\mathbf{X})$ varies slowly in space. In Section 7.3, we treat the case of circular arch, whose geometric properties vary slowly in space. In both case, the element stiffness matrix $\mathcal{H}_{\varphi}$ depends explicitly on $\mathbf{X}$, see Equation (7.22) for the arch.

### 3. HOMOGENIZATION PROBLEM IN CANONICAL FORM

In this Section, we cast the energy formulation of the discrete lattice problem obtained in Section 2 into a standard form, the canonical form, which we will be able to feed into the symbolic homogenization tool [Aud23] available from our previous work [AL23]. The following steps are involved:

* describe the discrete solution in terms of continuous functions (continualization step, §3.1),

* introduce a vector $l(\mathbf{X})$ collecting the degrees of freedom which are fixed during homogenization (§3.3) and a vector $y(\mathbf{X})$ collecting those that are relaxed during homogenization (§3.4),

* introduce an assembled strain vector $\mathbf{E}(\mathbf{X})$ in terms of $l(\mathbf{X})$, $y(\mathbf{X})$ and their gradients (§3.5),

* express the strain energy of the lattice in terms of $\mathbf{E}(\mathbf{X})$ (§3.6–3.7).
so that the strain energy of an element takes the form

$$w_{\alpha}(\mathbf{E}_{\alpha}) = \frac{1}{2} \mathbf{E}_{\alpha} \cdot \mathcal{H}_{\varphi_{\alpha}}(\eta, \mathbf{X}_{\alpha}^{c}) \cdot \mathbf{E}_{\alpha} \quad (2.20)$$

(no implicit sum over $\alpha$). To keep the presentation general, we will use this expression (2.20) of the strain energy in the following, as it is applies to linearly elastic element of various kinds such as springs, beams or arches in both 2D and 3D. The element stiffness matrix $\mathcal{H}_{\varphi} \in T^{(n_{E\varphi}, n_{E\varphi})}$ is always symmetric but not necessarily diagonal.

The element stiffness is assumed to obey the scaling assumption

$$\mathcal{H}_{\varphi} = \mathcal{O}(\eta^d), \quad (2.21)$$

where $d$ is the space dimension. This assumption is a matter of convention: applying a global factor to all elements' stiffness scales the homogenized energy by the same factor. The convention (2.21) warrants that the energy is $\Phi = \mathcal{O}(\eta^0)$, which is convenient, see Equation (3.42) below. Since $\hat{\ell}_{\varphi} = \mathcal{O}(\eta)$ by (2.13), the stretching modulus $\tilde{E}\tilde{A}_{\varphi}$ in (2.19) will be viewed as a quantity of order

$$EA = \mathcal{O}(\eta^{d-1}) \quad (2.22)$$

in stretching-dominated lattices.

**Remark 2.2.** As noted below Equation (2.17), the function $\hat{\chi}_{\varphi}(\eta, \mathbf{X})$ is small, $\hat{\chi}_{\varphi}(\eta, \mathbf{X}) \ll 1$. Mathematically, $\hat{\chi}_{\varphi}$ will be defined as a power of the aspect-ratio $\eta$. Different power laws are considered in the different examples, such as $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^2)$ for the honeycomb, see (5.7), or $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^0)$ for the pantograph, see (7.28). The only assumption we make at this stage is that $\hat{\chi}_{\varphi}$ is a function of $\eta$ with $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^0)$—this includes $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^1)$ and $\hat{\chi}_{\varphi} = \mathcal{O}(\eta^2)$ as particular cases. In any case, the matrix $\mathcal{H}_{\varphi}(\eta, \mathbf{X})$ appearing in (2.19) may contain terms of different orders in $\eta$ and may therefore depend on $\eta$ through more than just an overall scaling factor.

Anticipating that the homogenization will fail near the edges of the lattice due to the presence of boundary layers, we limit attention to a sub-region $\Omega \subset \mathbb{R}^d$ of the physical lattice that exclude the edges, see Figure 2.3: the strain energy of the lattice $\Phi_d$ in the region $\Omega$ is defined as the sum of the elastic energy $w_\alpha$ over all beams $\alpha$ whose center $\mathbf{X}_\alpha^c$ (shown as a red dot in Figure 2.3) lies inside the domain $\Omega$,

$$\Phi_d = \sum_{\alpha \text{ such that } \mathbf{X}_\alpha^c \in \Omega} w_\alpha(\mathbf{E}_\alpha). \quad (2.23)$$

Points where the loading is not smooth (such as points of application of point-like forces) are also excluded from the domain $\Omega$, along with a finite neighborhood around such points. The subscript ‘d’ in $\Phi_d$ stands for ‘discrete’: a continuous variant $\Phi$ of the lattice strain energy, bearing no subscript, will be introduced later.

<figure><img src="image_5.png" /><figcaption>**Figure 2.3.** To avoid boundary layers, the homogenization is carried out in a sub-region $\Omega$ of the lattice (shaded blue region) that excludes the physical boundaries of the lattice.</figcaption></figure>

**Remark 2.3.** In Sections 7.1, we treat the case of an elastic truss whose bending modulus $EI(\mathbf{X})$ varies slowly in space. In Section 7.3, we treat the case of circular arch, whose geometric properties vary slowly in space. In both case, the element stiffness matrix $\mathcal{H}_{\varphi}$ depends explicitly on $\mathbf{X}$, see Equation (7.22) for the arch.

### 3. HOMOGENIZATION PROBLEM IN CANONICAL FORM

In this Section, we cast the energy formulation of the discrete lattice problem obtained in Section 2 into a standard form, the canonical form, which we will be able to feed into the symbolic homogenization tool [Aud23] available from our previous work [AL23]. The following steps are involved:

* describe the discrete solution in terms of continuous functions (continualization step, §3.1),

* introduce a vector $l(\mathbf{X})$ collecting the degrees of freedom which are fixed during homogenization (§3.3) and a vector $y(\mathbf{X})$ collecting those that are relaxed during homogenization (§3.4),

* introduce an assembled strain vector $\mathbf{E}(\mathbf{X})$ in terms of $l(\mathbf{X})$, $y(\mathbf{X})$ and their gradients (§3.5),

* express the strain energy of the lattice in terms of $\mathbf{E}(\mathbf{X})$ (§3.6–3.7).
## 3.1. Parameterization of discrete solution using continuous functions

A key assumption in homogenization is that the equilibrium solution of the discrete lattice can be captured by the following functions of the continuous variable $\mathbf{X}$:

* a macroscopic displacement $\mathbf{u}(\mathbf{X}) \in \mathbb{R}^d$ representing the local average of the nodal displacement over the different Bravais sub-lattices,

* for each Bravais sub-lattice $1 \le b \le n_b$, a microscopic displacement $\xi_b(\mathbf{X}) \in \mathbb{R}^d$ and a microscopic rotation $\psi_b(\mathbf{X}) \in \mathbb{R}^{n_b}$.

Their orders of magnitude are taken as (see Remark 3.1 for a justification),

$$ \mathbf{u}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad \xi_b(\mathbf{X}) = \mathcal{O}(\eta^1), \quad \psi_b(\mathbf{X}) = \mathcal{O}(\eta^0). \qquad (3.1) $$

All these functions may depend on the expansion parameter $\eta \ll 1$, even though this dependence is implicit in our notation. As earlier in (2.10), we assume that this dependence takes place through regular expansions, with a leading term set by (3.1), i.e.,

$$ \begin{aligned} \mathbf{u}(\mathbf{X}) &= \mathbf{u}_{(0)}(\mathbf{X}) + \eta \mathbf{u}_{(1)}(\mathbf{X}) + \dots \\ \xi_b(\mathbf{X}) &= \eta(\xi_b)_{(1)}(\mathbf{X}) + \dots \\ \psi_b(\mathbf{X}) &= (\psi_b)_{(0)}(\mathbf{X}) + \eta(\psi_b)_{(1)}(\mathbf{X}) + \dots \end{aligned} \qquad (3.2) $$

where all the 'terms' in the series $\mathbf{u}_{(k)}(\mathbf{X})$, $(\xi_b)_{(k)}(\mathbf{X})$ and $(\psi_b)_{(k)}(\mathbf{X})$ are smooth functions of $\mathbf{X}$. Throughout this paper and in the interest of legibility, we will use the quantities $\mathbf{u}(\mathbf{X})$, $\xi_b(\mathbf{X})$ and $\psi_b(\mathbf{X})$ as blocks, refraining as much as possible from exposing the underlying expansions. The only concrete consequences of (3.2) are that (i) we determine these fields perturbatively, see (4.3–4.4), at the successive orders of the homogenization procedure and (ii) the gradients scale the same way as the undifferentiated quantities, i.e., $\nabla^k \mathbf{u} = \mathcal{O}(\eta^0)$, $\nabla^k \xi_b = \mathcal{O}(\eta^1)$ and $\nabla^k \psi_b = \mathcal{O}(\eta^0)$, which we summarize by the formal rule

$$ \nabla = \mathcal{O}(\eta^0). \qquad (3.3) $$

As usual in continuum mechanics, the macroscopic rotation $\gamma(\mathbf{X}) \in \mathbb{R}^{n_r}$ is a vector extracted from the antisymmetric part of $\nabla \mathbf{u}(\mathbf{X})$: in our notation, this extraction is denoted as $\gamma(\mathbf{X}) = -\frac{1}{2} \nabla \mathbf{u}(\mathbf{X}) : \mathcal{N}' = \mathcal{O}(\eta^0)$ where $\mathcal{N}'$ is a constant tensor that is antisymmetric with respect to its first pair of indices, see Equation (3.8)$_2$ below.

Following the Cauchy-Born 'hypothesis', the nodal degrees of freedom are postulated in terms of the continuous fields as

$$ \mathbf{v}_{\beta} = \mathbf{u}(\mathbf{X}_{\beta}) + \xi_{b(\beta)}(\mathbf{X}_{\beta}), \quad \theta_{\beta} = \gamma(\mathbf{X}_{\beta}) + \psi_{b(\beta)}(\mathbf{X}_{\beta}). \qquad (3.4) $$

Here, both the nodal displacement $\mathbf{v}_{\beta}$ and rotation $\theta_{\beta}$ are expressed as the sum of a macroscopic quantity (displacement $\mathbf{u}(\mathbf{X})$ and rotation $\gamma(\mathbf{X})$) and a microscopic one (displacement $\xi_b(\mathbf{X})$ and rotation $\psi_b(\mathbf{X})$). The macroscopic quantities $\mathbf{u}(\mathbf{X})$ and $\gamma(\mathbf{X})$ are functions on $\mathbf{X}$ only, and therefore vary slowly, whereas the additional dependence of $\xi_b$ and $\psi_b$ on the Bravais index $b = b(\beta)$ of the node causes 'fast' variations at the scale of the cell. Equation (3.4) is a variant of the double scale expansion used as a starting point for the homogenization of periodic continua, see for instance [San80; Bou96], which has been modified to account for the fact that $b$ is a discrete variable here.

The decomposition of the nodal displacement $\mathbf{v}_{\beta}$ in (3.4)$_1$ is not unique: shifting $\mathbf{u}(\mathbf{X})$ by $\tilde{\mathbf{u}}(\mathbf{X})$, and all the $\xi_b(\mathbf{X})$'s by $-\tilde{\mathbf{u}}(\mathbf{X})$ leaves $\mathbf{v}_{\beta}$ unchanged. To ensure unicity, we impose the constraint

$$ \sum_{b=1}^{n_b} \xi_b(\mathbf{X}) = 0, \quad \forall \mathbf{X}. \qquad (3.5) $$

By averaging both sides of (3.4)$_1$ with respect to the Bravais sub-lattices $b$, one can then confirm that $\mathbf{u}(\mathbf{X})$ is the average nodal displacement as announced earlier.

**Remark 3.1.** The scaling assumptions in (3.1) imply that the element strain $E_{\alpha}$ in (2.12) has magnitude

$$ E_{\alpha} = \mathcal{O}(\eta^0). \qquad (3.6) $$

Considering the stretching strain, for instance, $|\epsilon_{\alpha}| = \left| \frac{(\mathbf{v}_{\beta\alpha}^{+} - \mathbf{v}_{\beta\alpha}) \cdot \mathbf{t}_{\alpha}}{\ell_{\alpha}} \right| = \left| \frac{(u(\mathbf{X}_{\beta\alpha}^{+}) + \xi_{b(\beta\alpha)}^{+})(\mathbf{X}_{\beta\alpha}^{+}) - u(\mathbf{X}_{\beta\alpha}) - \xi_{b(\beta\alpha)}^{+})(\mathbf{X}_{\beta\alpha})}{\ell_{\alpha}} \right| \le \frac{|\nabla \mathbf{u} \cdot (\Delta_{\psi_{\alpha}} - \Delta_{\psi_{\alpha}})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} + \frac{|\xi_{b(\beta\alpha)}(\mathbf{X}_{\beta\alpha})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} + \frac{|\xi_{b(\beta\alpha)}(\mathbf{X}_{\beta\alpha})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} \quad \text{which yields } |\epsilon_{\alpha}| = \mathcal{O}(\eta^0) \text{ when combined with (3.1)$_{1,2}$ and (2.6). A similar reasoning holds for } \kappa_{\alpha} = \mathcal{O}(\eta^0) \text{ and } \tau_{\alpha} = \mathcal{O}(\eta^0)\text{. In situations where the strain scales differently than in (3.6), one can simply shift the orders in the scaling assumptions (3.1).}$$

**Remark 3.2.** Given an equilibrium configuration of the lattice, the continuous fields can be explicitly constructed as follows. First, the nodal displacements $\mathbf{v}_{\beta}$ in each particular Bravais sub-lattice $b \ni \beta$ are interpolated by smooth functions $\mathbf{v}_{\beta}^s(\mathbf{X}_{\beta})$ of their undeformed position $\mathbf{X}_{\beta}$: a set of $n_b$ such interpolating functions $\mathbf{v}_{\beta}^s$ with $b = 1, \dots, n_b$ are obtained. Next, $\mathbf{u}(\mathbf{X})$ is defined as the average $\mathbf{u}(\mathbf{X}) = \frac{1}{n_b} \sum_{b=1}^{n_b} \mathbf{v}_{\beta}^s(\mathbf{X})$. Finally, the 'continualized' microscopic fields $\xi_b(\mathbf{X})$ and $\psi_b(\mathbf{X})$ are obtained by interpolating $\mathbf{v}_{\beta} - \mathbf{u}(\mathbf{X}_{\beta})$ and $\theta_{\beta} - \gamma(\mathbf{X}_{\beta})$, respectively, over each Bravais sub-lattice separately. The constraint (3.5) is then automatically satisfied.
## 3.1. Parameterization of discrete solution using continuous functions

A key assumption in homogenization is that the equilibrium solution of the discrete lattice can be captured by the following functions of the continuous variable $\mathbf{X}$:

* a macroscopic displacement $\mathbf{u}(\mathbf{X}) \in \mathbb{R}^d$ representing the local average of the nodal displacement over the different Bravais sub-lattices,

* for each Bravais sub-lattice $1 \le b \le n_b$, a microscopic displacement $\xi_b(\mathbf{X}) \in \mathbb{R}^d$ and a microscopic rotation $\psi_b(\mathbf{X}) \in \mathbb{R}^{n_b}$.

Their orders of magnitude are taken as (see Remark 3.1 for a justification),

$$ \mathbf{u}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad \xi_b(\mathbf{X}) = \mathcal{O}(\eta^1), \quad \psi_b(\mathbf{X}) = \mathcal{O}(\eta^0). \qquad (3.1) $$

All these functions may depend on the expansion parameter $\eta \ll 1$, even though this dependence is implicit in our notation. As earlier in (2.10), we assume that this dependence takes place through regular expansions, with a leading term set by (3.1), i.e.,

$$ \begin{aligned} \mathbf{u}(\mathbf{X}) &= \mathbf{u}_{(0)}(\mathbf{X}) + \eta \mathbf{u}_{(1)}(\mathbf{X}) + \dots \\ \xi_b(\mathbf{X}) &= \eta(\xi_b)_{(1)}(\mathbf{X}) + \dots \\ \psi_b(\mathbf{X}) &= (\psi_b)_{(0)}(\mathbf{X}) + \eta(\psi_b)_{(1)}(\mathbf{X}) + \dots \end{aligned} \qquad (3.2) $$

where all the 'terms' in the series $\mathbf{u}_{(k)}(\mathbf{X})$, $(\xi_b)_{(k)}(\mathbf{X})$ and $(\psi_b)_{(k)}(\mathbf{X})$ are smooth functions of $\mathbf{X}$. Throughout this paper and in the interest of legibility, we will use the quantities $\mathbf{u}(\mathbf{X})$, $\xi_b(\mathbf{X})$ and $\psi_b(\mathbf{X})$ as blocks, refraining as much as possible from exposing the underlying expansions. The only concrete consequences of (3.2) are that (i) we determine these fields perturbatively, see (4.3–4.4), at the successive orders of the homogenization procedure and (ii) the gradients scale the same way as the undifferentiated quantities, i.e., $\nabla^k \mathbf{u} = \mathcal{O}(\eta^0)$, $\nabla^k \xi_b = \mathcal{O}(\eta^1)$ and $\nabla^k \psi_b = \mathcal{O}(\eta^0)$, which we summarize by the formal rule

$$ \nabla = \mathcal{O}(\eta^0). \qquad (3.3) $$

As usual in continuum mechanics, the macroscopic rotation $\gamma(\mathbf{X}) \in \mathbb{R}^{n_r}$ is a vector extracted from the antisymmetric part of $\nabla \mathbf{u}(\mathbf{X})$: in our notation, this extraction is denoted as $\gamma(\mathbf{X}) = -\frac{1}{2} \nabla \mathbf{u}(\mathbf{X}) : \mathcal{N}' = \mathcal{O}(\eta^0)$ where $\mathcal{N}'$ is a constant tensor that is antisymmetric with respect to its first pair of indices, see Equation (3.8)$_2$ below.

Following the Cauchy-Born 'hypothesis', the nodal degrees of freedom are postulated in terms of the continuous fields as

$$ \mathbf{v}_{\beta} = \mathbf{u}(\mathbf{X}_{\beta}) + \xi_{b(\beta)}(\mathbf{X}_{\beta}), \quad \theta_{\beta} = \gamma(\mathbf{X}_{\beta}) + \psi_{b(\beta)}(\mathbf{X}_{\beta}). \qquad (3.4) $$

Here, both the nodal displacement $\mathbf{v}_{\beta}$ and rotation $\theta_{\beta}$ are expressed as the sum of a macroscopic quantity (displacement $\mathbf{u}(\mathbf{X})$ and rotation $\gamma(\mathbf{X})$) and a microscopic one (displacement $\xi_b(\mathbf{X})$ and rotation $\psi_b(\mathbf{X})$). The macroscopic quantities $\mathbf{u}(\mathbf{X})$ and $\gamma(\mathbf{X})$ are functions on $\mathbf{X}$ only, and therefore vary slowly, whereas the additional dependence of $\xi_b$ and $\psi_b$ on the Bravais index $b = b(\beta)$ of the node causes 'fast' variations at the scale of the cell. Equation (3.4) is a variant of the double scale expansion used as a starting point for the homogenization of periodic continua, see for instance [San80; Bou96], which has been modified to account for the fact that $b$ is a discrete variable here.

The decomposition of the nodal displacement $\mathbf{v}_{\beta}$ in (3.4)$_1$ is not unique: shifting $\mathbf{u}(\mathbf{X})$ by $\tilde{\mathbf{u}}(\mathbf{X})$, and all the $\xi_b(\mathbf{X})$'s by $-\tilde{\mathbf{u}}(\mathbf{X})$ leaves $\mathbf{v}_{\beta}$ unchanged. To ensure unicity, we impose the constraint

$$ \sum_{b=1}^{n_b} \xi_b(\mathbf{X}) = 0, \quad \forall \mathbf{X}. \qquad (3.5) $$

By averaging both sides of (3.4)$_1$ with respect to the Bravais sub-lattices $b$, one can then confirm that $\mathbf{u}(\mathbf{X})$ is the average nodal displacement as announced earlier.

**Remark 3.1.** The scaling assumptions in (3.1) imply that the element strain $E_{\alpha}$ in (2.12) has magnitude

$$ E_{\alpha} = \mathcal{O}(\eta^0). \qquad (3.6) $$

Considering the stretching strain, for instance, $|\varepsilon_{\alpha}| = \left| \frac{(\mathbf{v}_{\beta\alpha}^{+} - \mathbf{v}_{\beta\alpha}) \cdot \mathbf{t}_{\alpha}}{\ell_{\alpha}} \right| = \left| \frac{(u(\mathbf{X}_{\beta\alpha}^{+}) + \xi_{b(\beta\alpha)}^{+})(\mathbf{X}_{\beta\alpha}^{+}) - u(\mathbf{X}_{\beta\alpha}) - \xi_{b(\beta\alpha)}^{+})(\mathbf{X}_{\beta\alpha})}{\ell_{\alpha}} \right| \le \frac{|\nabla \mathbf{u} \cdot (\Delta_{\psi_{\alpha}} - \Delta_{\psi_{\alpha}})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} + \frac{|\xi_{b(\beta\alpha)}(\mathbf{X}_{\beta\alpha})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} + \frac{|\xi_{b(\beta\alpha)}(\mathbf{X}_{\beta\alpha})|}{|\Delta_{\psi_{\alpha}}^{+} - \Delta_{\psi_{\alpha}}|} \quad \text{which yields } |\varepsilon_{\alpha}| = \mathcal{O}(\eta^0) \text{ when combined with (3.1)$_{1,2}$ and (2.6). A similar reasoning holds for } \kappa_{\alpha} = \mathcal{O}(\eta^0) \text{ and } \tau_{\alpha} = \mathcal{O}(\eta^0)\text{. In situations where the strain scales differently than in (3.6), one can simply shift the orders in the scaling assumptions (3.1).}$$

**Remark 3.2.** Given an equilibrium configuration of the lattice, the continuous fields can be explicitly constructed as follows. First, the nodal displacements $\mathbf{v}_{\beta}$ in each particular Bravais sub-lattice $b \ni \beta$ are interpolated by smooth functions $\mathbf{v}_{\beta}^s(\mathbf{X}_{\beta})$ of their undeformed position $\mathbf{X}_{\beta}$: a set of $n_b$ such interpolating functions $\mathbf{v}_{\beta}^s$ with $b = 1, \dots, n_b$ are obtained. Next, $\mathbf{u}(\mathbf{X})$ is defined as the average $\mathbf{u}(\mathbf{X}) = \frac{1}{n_b} \sum_{b=1}^{n_b} \mathbf{v}_{\beta}^s(\mathbf{X})$. Finally, the 'continualized' microscopic fields $\xi_b(\mathbf{X})$ and $\psi_b(\mathbf{X})$ are obtained by interpolating $\mathbf{v}_{\beta} - \mathbf{u}(\mathbf{X}_{\beta})$ and $\theta_{\beta} - \gamma(\mathbf{X}_{\beta})$, respectively, over each Bravais sub-lattice separately. The constraint (3.5) is then automatically satisfied.
**Remark 3.3.** Consider the particular case of a rigid-body motion of the lattice made up of a translation $\mathbf{u}_0$ and rotation $\gamma_0$, i.e., $(\mathbf{v}_\beta, \theta_\beta) = (\mathbf{u}_0 + \gamma_0 \times \mathbf{X}_\beta, \gamma_0)$ at every node $\beta$. The macroscopic fields are extracted as outlined in Remark 3.2, which yields the macroscopic displacement as $\mathbf{u}(\mathbf{X}) = \mathbf{u}_0 + \gamma_0 \times \mathbf{X}$, a constant macroscopic rotation $\gamma(\mathbf{X}) = \gamma_0$ given by (3.8)$_2$ and zero microscopic displacement $\xi_b(\mathbf{X}) = \mathbf{0}_d$ and rotation $\psi_b(\mathbf{X}) = \mathbf{0}_{n_r}$. By including the term $\gamma(\mathbf{X}_\beta)$ in the right-hand side of (3.4)$_2$, we have made the microscopic displacement insensitive to rigid-body motions of the lattice.

## 3.2. Macroscopic strain $\tilde{\mathbf{\epsilon}}$ and rotation $\gamma$

In linear elasticity, the displacement gradient $\nabla \mathbf{u}(\mathbf{X})$ is classically decomposed into:

* a symmetric strain tensor $\varepsilon(\mathbf{X}) = \frac{1}{2}(\nabla \mathbf{u}(\mathbf{X}) + \nabla \mathbf{u}^T(\mathbf{X})) \in \mathbb{T}^{(d,d)}$, whose components we arrange into a vector $\tilde{\mathbf{\epsilon}}(\mathbf{X}) \in \mathbb{R}^{n_r}$, following Mandel's notation,

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \begin{cases} (\varepsilon_{11}(\mathbf{X}), \varepsilon_{22}(\mathbf{X}), \sqrt{2}\varepsilon_{12}(\mathbf{X})) & (d=2) \\ (\varepsilon_{11}(\mathbf{X}), \varepsilon_{22}(\mathbf{X}), \varepsilon_{33}(\mathbf{X}), \sqrt{2}\varepsilon_{12}(\mathbf{X}), \sqrt{2}\varepsilon_{23}(\mathbf{X}), \sqrt{2}\varepsilon_{13}(\mathbf{X})) & (d=3), \end{cases} \quad (3.7) $$

where $n_\varepsilon = \frac{d(d+1)}{2}$ is the number of independent strain components,

* and an antisymmetric part $\hat{\gamma}(\mathbf{X}) = \frac{1}{2}(\nabla \mathbf{u}(\mathbf{X}) - \nabla \mathbf{u}^T(\mathbf{X}))$ representing an infinitesimal rotation.

With the help of the geometric tensors $\mathcal{M} \in \mathbb{T}^{(d,d,n_r)}$ and $\mathcal{N}^* \in \mathbb{T}^{(d,d,n_r)}$, which depend on the space dimension $d$ only and are defined in Equation (B.1) in the Appendix B and Equation (A.6) in the Appendix A.2, respectively, the macroscopic strain $\tilde{\mathbf{\epsilon}}$ and the pseudo-vector $\gamma$ associated with the infinitesimal rotation $\hat{\gamma}$ are obtained from $\nabla \mathbf{u}$ as

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \nabla \mathbf{u}(\mathbf{X}) : \mathcal{M}, \quad \gamma(\mathbf{X}) = -\frac{1}{2} \nabla \mathbf{u}(\mathbf{X}) : \mathcal{N}. \quad (3.8) $$

The first equality follows from Equation (B.4) in Appendix B, and from the fact that $\mathcal{M}$ is symmetric with respect to its first pair of indices. The second equality follows by setting $\gamma = \theta$ in (A.5)$_2$ and from the fact that $\mathcal{N}^*$ is antisymmetric with respect to its first pair of indices.

By the scaling assumptions (3.1), $\nabla \mathbf{u} = \mathcal{O}(\eta^0)$, and we have the order-of-magnitude estimates

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad \gamma(\mathbf{X}) = \mathcal{O}(\eta^0). \quad (3.9) $$

## 3.3. Macroscopic state vector $l$

Next, we define the macroscopic degrees of freedom vector $l(\mathbf{X}) \in \mathbb{R}^{n_l}$ by concatenating the macroscopic displacement $\mathbf{u}(\mathbf{X})$, rotation $\gamma(\mathbf{X})$ and strain $\tilde{\mathbf{\epsilon}}(\mathbf{X})$: in block-vector notation,

$$ l(\mathbf{X}) = (\mathbf{u}(\mathbf{X}) \quad \gamma(\mathbf{X}) \quad \tilde{\mathbf{\epsilon}}(\mathbf{X})) \in \mathbb{R}^{n_l}, \quad (3.10) $$

where

$$ n_l = d + n_r + n_\varepsilon. \quad (3.11) $$

By design, the macroscopic state vector is a quantity of order 1, see (3.1) and (3.8),

$$ l(\mathbf{X}) = \mathcal{O}(\eta^0). \quad (3.12) $$

The macroscopic state vector $l(\mathbf{X})$ collects the quantities that are fixed during homogenization.

As shown in Equation (A.8) in Appendix A.3, the Taylor expansion about a point $\mathbf{X}_c$ of the macroscopic contribution $(\mathbf{u}, \gamma)$ to the nodal degrees of freedom, see (3.4), can be written in terms of the macroscopic state vector $l(\mathbf{X})$ as

$$ \left( \begin{array}{c} \mathbf{u}(\mathbf{X}) \\ \gamma(\mathbf{X}) \end{array} \right) = \mathcal{T}_l(\mathbf{X}-\mathbf{X}_c) \cdot l(\mathbf{X}_c) + \mathcal{T}_{l'}(\mathbf{X}-\mathbf{X}_c) : \nabla l(\mathbf{X}_c) + \mathcal{T}_{l''}(\mathbf{X}-\mathbf{X}_c) : \nabla^2 l(\mathbf{X}_c) + \dots, \quad (3.13) $$

where the tensors $\mathcal{T}_{l(\varepsilon)}$ are expressed in block-matrix notation as

$$ \begin{align} \mathcal{T}_l(\Delta) &= \begin{pmatrix} I_d & -\mathcal{N}^{T_{132}} \cdot \Delta & \mathcal{U}'^{T_{132}} \cdot \Delta \\ \mathbf{0}_{n_r \times d} & I_{n_r} & \mathbf{0}_{n_r \times n_r} \end{pmatrix} \in \mathbb{T}^{(n_n, n_l)}, \notag \\ \mathcal{T}_{l'}(\Delta) &= \begin{pmatrix} \mathbf{0}_{d \times d \times d} & \mathbf{0}_{d \times n_r \times d} & \mathcal{U}''^{T_{14523}} : \frac{\Delta^{\otimes 2}}{2} \\ \mathbf{0}_{n_r \times d \times d} & \mathbf{0}_{n_r \times n_r \times d} & \mathcal{G}''^{T_{1423}} \cdot \Delta \end{pmatrix} \in \mathbb{T}^{(n_n, n_l, d)}, \tag{3.14} \\ \mathcal{T}_{l''}(\Delta) &= \begin{pmatrix} \mathbf{0}_{d \times d \times d \times d} & \mathbf{0}_{d \times n_r \times d \times d} & \mathcal{U}'''^{T_{1567234}} : \frac{\Delta^{\otimes 3}}{6} \\ \mathbf{0}_{n_r \times d \times d \times d} & \mathbf{0}_{n_r \times n_r \times d \times d} & \mathcal{G}'''^{T_{156234}} : \frac{\Delta^{\otimes 2}}{2} \end{pmatrix} \in \mathbb{T}^{(n_n, n_l, d, d)}, \notag \end{align} $$

Equations (3.13–3.14) follow from the geometric analysis in Sections A.3 and B.3 in the Appendix, and make use of the geometric tensors $\mathcal{N}$, $\mathcal{U}'$, $\mathcal{U}''$, $\mathcal{U}'''$, $\mathcal{G}'$ and $\mathcal{G}''$ defined in Equations (B.7–B.15): they are a compact rewriting of the Taylor expansion (A.8) derived in the Appendix. In the right-hand sides of (3.14), the vertical blocks correspond to the nodal displacement $(\mathbf{u})$ versus rotation $(\gamma)$, and the horizontal blocks to the $\nabla^i \mathbf{u}$, $\nabla^i \gamma$ and $\nabla^i \tilde{\mathbf{\epsilon}}$ parts of $\nabla^i l$, see (3.10). The generalized-transpose notation is defined in Appendix A.
**Remark 3.3.** Consider the particular case of a rigid-body motion of the lattice made up of a translation $\mathbf{u}_0$ and rotation $\gamma_0$, i.e., $(\mathbf{v}_\beta, \theta_\beta) = (\mathbf{u}_0 + \gamma_0 \times \mathbf{X}_\beta, \gamma_0)$ at every node $\beta$. The macroscopic fields are extracted as outlined in Remark 3.2, which yields the macroscopic displacement as $\mathbf{u}(\mathbf{X}) = \mathbf{u}_0 + \gamma_0 \times \mathbf{X}$, a constant macroscopic rotation $\gamma(\mathbf{X}) = \gamma_0$ given by (3.8)$_2$ and zero microscopic displacement $\xi_b(\mathbf{X}) = \mathbf{0}_d$ and rotation $\psi_b(\mathbf{X}) = \mathbf{0}_{n_r}$. By including the term $\gamma(\mathbf{X}_\beta)$ in the right-hand side of (3.4)$_2$, we have made the microscopic displacement insensitive to rigid-body motions of the lattice.

## 3.2. Macroscopic strain $\tilde{\mathbf{\epsilon}}$ and rotation $\gamma$

In linear elasticity, the displacement gradient $\nabla \mathbf{u}(\mathbf{X})$ is classically decomposed into:

* a symmetric strain tensor $\varepsilon(\mathbf{X}) = \frac{1}{2}(\nabla \mathbf{u}(\mathbf{X}) + \nabla \mathbf{u}^T(\mathbf{X})) \in \mathbb{T}^{(d,d)}$, whose components we arrange into a vector $\tilde{\mathbf{\epsilon}}(\mathbf{X}) \in \mathbb{R}^{n_r}$, following Mandel's notation,

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \begin{cases} (\varepsilon_{11}(\mathbf{X}), \varepsilon_{22}(\mathbf{X}), \sqrt{2}\varepsilon_{12}(\mathbf{X})) & (d=2) \\ (\varepsilon_{11}(\mathbf{X}), \varepsilon_{22}(\mathbf{X}), \varepsilon_{33}(\mathbf{X}), \sqrt{2}\varepsilon_{12}(\mathbf{X}), \sqrt{2}\varepsilon_{23}(\mathbf{X}), \sqrt{2}\varepsilon_{13}(\mathbf{X})) & (d=3), \end{cases} \quad (3.7) $$

where $n_\varepsilon = \frac{d(d+1)}{2}$ is the number of independent strain components,

* and an antisymmetric part $\hat{\gamma}(\mathbf{X}) = \frac{1}{2}(\nabla \mathbf{u}(\mathbf{X}) - \nabla \mathbf{u}^T(\mathbf{X}))$ representing an infinitesimal rotation.

With the help of the geometric tensors $\mathcal{M} \in \mathbb{T}^{(d,d,n_r)}$ and $\mathcal{N}^* \in \mathbb{T}^{(d,d,n_r)}$, which depend on the space dimension $d$ only and are defined in Equation (B.1) in the Appendix B and Equation (A.6) in the Appendix A.2, respectively, the macroscopic strain $\tilde{\mathbf{\epsilon}}$ and the pseudo-vector $\gamma$ associated with the infinitesimal rotation $\hat{\gamma}$ are obtained from $\nabla \mathbf{u}$ as

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \nabla \mathbf{u}(\mathbf{X}) : \mathcal{M}, \quad \gamma(\mathbf{X}) = -\frac{1}{2} \nabla \mathbf{u}(\mathbf{X}) : \mathcal{N}. \quad (3.8) $$

The first equality follows from Equation (B.4) in Appendix B, and from the fact that $\mathcal{M}$ is symmetric with respect to its first pair of indices. The second equality follows by setting $\gamma = \theta$ in (A.5)$_2$ and from the fact that $\mathcal{N}^*$ is antisymmetric with respect to its first pair of indices.

By the scaling assumptions (3.1), $\nabla \mathbf{u} = \mathcal{O}(\eta^0)$, and we have the order-of-magnitude estimates

$$ \tilde{\mathbf{\epsilon}}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad \gamma(\mathbf{X}) = \mathcal{O}(\eta^0). \quad (3.9) $$

## 3.3. Macroscopic state vector $\mathbf{l}$

Next, we define the macroscopic degrees of freedom vector $\mathbf{l}(\mathbf{X}) \in \mathbb{R}^{n_l}$ by concatenating the macroscopic displacement $\mathbf{u}(\mathbf{X})$, rotation $\gamma(\mathbf{X})$ and strain $\tilde{\mathbf{\epsilon}}(\mathbf{X})$: in block-vector notation,

$$ \mathbf{l}(\mathbf{X}) = (\mathbf{u}(\mathbf{X}) \quad \gamma(\mathbf{X}) \quad \tilde{\mathbf{\epsilon}}(\mathbf{X})) \in \mathbb{R}^{n_l}, \quad (3.10) $$

where

$$ n_l = d + n_r + n_\varepsilon. \quad (3.11) $$

By design, the macroscopic state vector is a quantity of order 1, see (3.1) and (3.8),

$$ \mathbf{l}(\mathbf{X}) = \mathcal{O}(\eta^0). \quad (3.12) $$

The macroscopic state vector $\mathbf{l}(\mathbf{X})$ collects the quantities that are fixed during homogenization.

As shown in Equation (A.8) in Appendix A.3, the Taylor expansion about a point $\mathbf{X}_c$ of the macroscopic contribution $(\mathbf{u}, \gamma)$ to the nodal degrees of freedom, see (3.4), can be written in terms of the macroscopic state vector $\mathbf{l}(\mathbf{X})$ as

$$ \left( \begin{array}{c} \mathbf{u}(\mathbf{X}) \\ \gamma(\mathbf{X}) \end{array} \right) = \mathcal{T}_l(\mathbf{X}-\mathbf{X}_c) \cdot \mathbf{l}(\mathbf{X}_c) + \mathcal{T}_{l'}(\mathbf{X}-\mathbf{X}_c) : \nabla \mathbf{l}(\mathbf{X}_c) + \mathcal{T}_{l''}(\mathbf{X}-\mathbf{X}_c) : \nabla^2 \mathbf{l}(\mathbf{X}_c) + \dots, \quad (3.13) $$

where the tensors $\mathcal{T}_{l(i)}$ are expressed in block-matrix notation as

$$ \begin{align} \mathcal{T}_l(\Delta) &= \begin{pmatrix} \mathbf{I}_d & -\mathcal{N}^{T_{132}} \cdot \Delta & \mathcal{U}'^{T_{132}} \cdot \Delta \\ \mathbf{0}_{n_r \times d} & \mathbf{I}_{n_r} & \mathbf{0}_{n_r \times n_r} \end{pmatrix} \in \mathbb{T}^{(n_n, n_l)}, \notag \\ \mathcal{T}_{l'}(\Delta) &= \begin{pmatrix} \mathbf{0}_{d \times d \times d} & \mathbf{0}_{d \times n_r \times d} & \mathcal{U}''^{T_{14523}} : \frac{\Delta^{\otimes 2}}{2} \\ \mathbf{0}_{n_r \times d \times d} & \mathbf{0}_{n_r \times n_r \times d} & \mathcal{G}'^{T_{1423}} \cdot \Delta \end{pmatrix} \in \mathbb{T}^{(n_n, n_l, d)}, \tag{3.14} \\ \mathcal{T}_{l''}(\Delta) &= \begin{pmatrix} \mathbf{0}_{d \times d \times d \times d} & \mathbf{0}_{d \times n_r \times d \times d} & \mathcal{U}'''^{T_{1567234}} : \frac{\Delta^{\otimes 3}}{6} \\ \mathbf{0}_{n_r \times d \times d \times d} & \mathbf{0}_{n_r \times n_r \times d \times d} & \mathcal{G}''^{T_{156234}} : \frac{\Delta^{\otimes 2}}{2} \end{pmatrix} \in \mathbb{T}^{(n_n, n_l, d, d)}, \notag \\ \dots & \end{align} $$

Equations (3.13–3.14) follow from the geometric analysis in Sections A.3 and B.3 in the Appendix, and make use of the geometric tensors $\mathcal{N}$, $\mathcal{U}'$, $\mathcal{U}''$, $\mathcal{U}''$, $\mathcal{G}'$ and $\mathcal{G}''$ defined in Equations (B.7–B.15): they are a compact rewriting of the Taylor expansion (A.8) derived in the Appendix. In the right-hand sides of (3.14), the vertical blocks correspond to the nodal displacement $(\mathbf{u})$ versus rotation $(\gamma)$, and the horizontal blocks to the $\nabla^i \mathbf{u}$, $\nabla^i \gamma$ and $\nabla^i \tilde{\mathbf{\epsilon}}$ parts of $\nabla^i \mathbf{l}$, see (3.10). The generalized-transpose notation is defined in Appendix A.
### 3.4. Microscopic degrees of freedom y

Next, we define the microscopic degrees of freedom $y(\mathbf{X})$ as the concatenation into a single vector of those relevant to the different Bravais lattice,

$$ \mathbf{y}(\mathbf{X}) = \left( \frac{\xi_1(\mathbf{X})}{\eta} \ \psi_1(\mathbf{X}) \ \dots \ \frac{\xi_b(\mathbf{X})}{\eta} \ \psi_b(\mathbf{X}) \ \dots \ \frac{\xi_{n_b}(\mathbf{X})}{\eta} \ \psi_{n_b}(\mathbf{X}) \right) \in \mathbb{R}^{n_y}, \quad (3.15) $$

where $n_y$ is the dimension of $y$,

$$ n_y = n_b \times n_n. \quad (3.16) $$

By contrast with $l$, the vector $y$ lists all the degrees of freedom that will be relaxed during homogenization: the discrete energy of the lattice depends on $y$ but the homogenized energy does not. The coefficients $1/\eta$ in (3.15) warrant

$$ \mathbf{y}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad (3.17) $$

see (3.1).

The microscopic displacement and rotation at a particular node can be extracted from $y$ using

$$ \begin{pmatrix} \xi_b(\beta)(\mathbf{X}_\beta) \\ \psi_b(\beta)(\mathbf{X}_\beta) \end{pmatrix} = \mathcal{B}_{b(\beta)} \cdot \mathbf{y}(\mathbf{X}_\beta), \quad (3.18) $$

where $\mathcal{B}_b$ is defined as

$$ \mathcal{B}_b = \left( \mathbf{0}_{n_n \times (b-1)n_n} \begin{pmatrix} \eta I_d & \mathbf{0} \\ \mathbf{0} & I_{n_n} \end{pmatrix} \mathbf{0}_{n_n \times (n_b-b)n_n} \right) \in \mathbb{T}^{(n_n, n_y)}. $$

We have derived in (3.13) the Taylor expansion of the macroscopic contribution to the nodal degrees of freedom about a generic point $\mathbf{X}_c$. The counterpart for the microscopic contributions is obtained by expanding (3.18) as

$$ \begin{pmatrix} \xi_b(\mathbf{X}) \\ \psi_b(\mathbf{X}) \end{pmatrix} = \mathcal{T}_y^b(\mathbf{X}-\mathbf{X}_c) \cdot \mathbf{y}(\mathbf{X}_c) + \mathcal{T}_{y'}^b(\mathbf{X}-\mathbf{X}_c) : \nabla \mathbf{y}(\mathbf{X}_c) + \mathcal{T}_{y''}^b(\mathbf{X}-\mathbf{X}_c) : \nabla^2 \mathbf{y}(\mathbf{X}_c) + \dots \quad (3.19) $$

where

$$ \begin{align} \mathcal{T}_{y'}^b(\Delta) &= \mathcal{B}_b \quad \in \mathbb{T}^{(n_n, n_y)}, \nonumber \\ \mathcal{T}_{y'}^b(\Delta) &= \mathcal{B}_b \otimes \Delta \quad \in \mathbb{T}^{(n_n, n_y, d)}, \tag{3.20} \\ \mathcal{T}_{y''}^b(\Delta) &= \mathcal{B}_b \otimes \frac{\Delta^{\otimes 2}}{2} \quad \in \mathbb{T}^{(n_n, n_y, d, d)}, \nonumber \\ \dots \nonumber \end{align} $$

### 3.5. Element strain expansion

Combining (2.14), (3.4), (3.13), (3.19), (2.6) and (2.2), one can rewrite the strain $\mathbf{E}_\alpha$ in an element $\alpha$ as

$$ \mathbf{E}_\alpha = \mathbf{E}_{\phi_\alpha}(\eta, \mathbf{X}_\alpha^c) \quad (3.21) $$

where $\mathbf{E}_\phi(\eta, \mathbf{X})$ is a smooth function of $\mathbf{X} = \mathbf{X}_\alpha^c$ depending on the element family $\phi = \phi_\alpha$,

$$ \mathbf{E}_\phi(\eta, \mathbf{X}) = \mathcal{D}_\phi(\eta, \mathbf{X}) : \left( \left[ \sum_{i=0}^{\infty} \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^-(\eta, \mathbf{X})) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) \cdots \right] \left[ \sum_{i=0}^{\infty} \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) \cdots \right] \right), $$

and $\mathbf{I}^{(i)}$ is a shorthand for the contraction of order $i$ defined in Section A.1 of the Appendix, that is $\mathbf{I}^{(1)}=\mathbf{I}, \mathbf{I}^{(2)}=\mathbf{I}, \mathbf{I}^{(3)}=\mathbf{I}, \dots$.

This can be rewritten in compact form as

$$ \mathbf{E}_{\phi}(\eta, \mathbf{X}) = \sum_{i=0}^{\infty} \mathbf{E}_{\phi}^{l^{(i)}}(\eta, \mathbf{X}) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) + \sum_{i=0}^{\infty} \mathbf{E}_{\phi}^{y^{(i)}}(\eta, \mathbf{X}) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{y}(\mathbf{X}), \quad (3.22) $$

where

$$ \begin{align} \mathbf{E}_{\phi}^{l^{(i)}}(\eta, \mathbf{X}) &= \mathcal{D}_{\phi}(\eta, \mathbf{X}) : \left( \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^-(\eta, \mathbf{X})) \cdot \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \right) \nonumber \\ \mathbf{E}_{\phi}^{y^{(i)}}(\eta, \mathbf{X}) &= \mathcal{D}_{\phi}(\eta, \mathbf{X}) : \left( \mathcal{T}_{y^{(i)}}^{b\bar{\phi}}(\Delta_{\bar{\phi}}(\eta, \mathbf{X})) \cdot \mathcal{T}_{y^{(i)}}^{b\phi}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \right). \tag{3.23} \end{align} $$

In (3.21–3.23), the element strain $\mathbf{E}_\alpha$ has been obtained as a series expansion about the element center $\mathbf{X} = \mathbf{X}_\alpha^c$ in terms of the macroscopic state vector $\mathbf{l}$, of the microscopic degrees of freedom $y$ and of their successive derivatives. Using explicit formulas for the $\mathcal{T}$'s in (3.14) and (3.20), we can use (3.23) to calculate symbolically and in closed form the coefficients $\mathbf{E}_\phi^l(\eta, \mathbf{X}), \mathbf{E}_\phi^l'(\eta, \mathbf{X}), \mathbf{E}_\phi^y(\eta, \mathbf{X}), \mathbf{E}_\phi^y'(\eta, \mathbf{X}), \dots$ entering in the series expansion (3.22) of the element strain. These coefficients depend on both $\eta$ and $\mathbf{X}$ when the lattice has non-uniform geometrical properties in space, as happens for instance with the curved lattice shown in Figure 2.1b.
### 3.4. Microscopic degrees of freedom y

Next, we define the microscopic degrees of freedom $y(\mathbf{X})$ as the concatenation into a single vector of those relevant to the different Bravais lattice,

$$ \mathbf{y}(\mathbf{X}) = \left( \frac{\xi_1(\mathbf{X})}{\eta} \ \psi_1(\mathbf{X}) \ \dots \ \frac{\xi_b(\mathbf{X})}{\eta} \ \psi_b(\mathbf{X}) \ \dots \ \frac{\xi_{n_b}(\mathbf{X})}{\eta} \ \psi_{n_b}(\mathbf{X}) \right) \in \mathbb{R}^{n_y}, \quad (3.15) $$

where $n_y$ is the dimension of $y$,

$$ n_y = n_b \times n_n. \quad (3.16) $$

By contrast with $l$, the vector $y$ lists all the degrees of freedom that will be relaxed during homogenization: the discrete energy of the lattice depends on $y$ but the homogenized energy does not. The coefficients $1/\eta$ in (3.15) warrant

$$ \mathbf{y}(\mathbf{X}) = \mathcal{O}(\eta^0), \quad (3.17) $$

see (3.1).

The microscopic displacement and rotation at a particular node can be extracted from $y$ using

$$ \begin{pmatrix} \xi_b(\beta)(\mathbf{X}_\beta) \\ \psi_b(\beta)(\mathbf{X}_\beta) \end{pmatrix} = \mathcal{B}_{b(\beta)} \cdot \mathbf{y}(\mathbf{X}_\beta), \quad (3.18) $$

where $\mathcal{B}_b$ is defined as

$$ \mathcal{B}_b = \left( \mathbf{0}_{n_n \times (b-1)n_n} \begin{pmatrix} \eta \mathbf{I}_d & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{n_n} \end{pmatrix} \mathbf{0}_{n_n \times (n_b-b)n_n} \right) \in \mathbb{T}^{(n_n, n_y)}. $$

We have derived in (3.13) the Taylor expansion of the macroscopic contribution to the nodal degrees of freedom about a generic point $\mathbf{X}_c$. The counterpart for the microscopic contributions is obtained by expanding (3.18) as

$$ \begin{pmatrix} \xi_b(\mathbf{X}) \\ \psi_b(\mathbf{X}) \end{pmatrix} = \mathcal{T}_y^b(\mathbf{X}-\mathbf{X}_c) \cdot \mathbf{y}(\mathbf{X}_c) + \mathcal{T}_{y'}^b(\mathbf{X}-\mathbf{X}_c) : \nabla \mathbf{y}(\mathbf{X}_c) + \mathcal{T}_{y''}^b(\mathbf{X}-\mathbf{X}_c) : \nabla^2 \mathbf{y}(\mathbf{X}_c) + \dots \quad (3.19) $$

where

$$ \begin{align} \mathcal{T}_{y'}^b(\Delta) &= \mathcal{B}_b \quad \in \mathbb{T}^{(n_n, n_y)}, \notag \\ \mathcal{T}_{y'}^b(\Delta) &= \mathcal{B}_b \otimes \Delta \quad \in \mathbb{T}^{(n_n, n_y, d)}, \tag{3.20} \\ \mathcal{T}_{y''}^b(\Delta) &= \mathcal{B}_b \otimes \frac{\Delta^{\otimes 2}}{2} \quad \in \mathbb{T}^{(n_n, n_y, d, d)}, \notag \\ \dots \notag \end{align} $$

### 3.5. Element strain expansion

Combining (2.14), (3.4), (3.13), (3.19), (2.6) and (2.2), one can rewrite the strain $\mathbf{E}_\alpha$ in an element $\alpha$ as

$$ \mathbf{E}_\alpha = \mathbf{E}_{\phi_\alpha}(\eta, \mathbf{X}_\alpha^c) \quad (3.21) $$

where $\mathbf{E}_\phi(\eta, \mathbf{X})$ is a smooth function of $\mathbf{X} = \mathbf{X}_\alpha^c$ depending on the element family $\phi = \phi_\alpha$,

$$ \mathbf{E}_\phi(\eta, \mathbf{X}) = \mathcal{D}_\phi(\eta, \mathbf{X}) : \left( \left[ \sum_{i=0}^{\infty} \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^-(\eta, \mathbf{X})) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) \cdots \right] \left[ \sum_{i=0}^{\infty} \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) \cdots \right] \right), $$

and $\mathbf{I}^{(i)}$ is a shorthand for the contraction of order $i$ defined in Section A.1 of the Appendix, that is $\mathbf{I}^{(1)} = \mathbf{I}, \mathbf{I}^{(2)} = \mathbf{I}, \dots, \mathbf{I}^{(3)} = \mathbf{I}, \dots$.

This can be rewritten in compact form as

$$ \mathbf{E}_{\phi}(\eta, \mathbf{X}) = \sum_{i=0}^{\infty} \mathbf{E}_{\phi}^{l^{(i)}}(\eta, \mathbf{X}) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{l}(\mathbf{X}) + \sum_{i=0}^{\infty} \mathbf{E}_{\phi}^{y^{(i)}}(\eta, \mathbf{X}) \cdot \mathbf{I}^{(i+1)} \nabla^i \mathbf{y}(\mathbf{X}), \quad (3.22) $$

where

$$ \begin{align} \mathbf{E}_{\phi}^{l^{(i)}}(\eta, \mathbf{X}) &= \mathcal{D}_{\phi}(\eta, \mathbf{X}) : \left( \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^-(\eta, \mathbf{X})) \cdot \mathcal{T}_{l^{(i)}}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \right) \notag \\ \mathbf{E}_{\phi}^{y^{(i)}}(\eta, \mathbf{X}) &= \mathcal{D}_{\phi}(\eta, \mathbf{X}) : \left( \mathcal{T}_{y^{(i)}}^{b\bar{\phi}}(\Delta_{\bar{\phi}}^-(\eta, \mathbf{X})) \cdot \mathcal{T}_{y^{(i)}}^{b\phi}(\Delta_{\bar{\phi}}^+(\eta, \mathbf{X})) \right). \tag{3.23} \end{align} $$

In (3.21–3.23), the element strain $\mathbf{E}_\alpha$ has been obtained as a series expansion about the element center $\mathbf{X} = \mathbf{X}_\alpha^c$ in terms of the macroscopic state vector $\mathbf{l}$, of the microscopic degrees of freedom $y$ and of their successive derivatives. Using explicit formulas for the $\mathcal{T}$'s in (3.14) and (3.20), we can use (3.23) to calculate symbolically and in closed form the coefficients $\mathbf{E}_\phi^l(\eta, \mathbf{X}), \mathbf{E}_\phi^l'(\eta, \mathbf{X}), \mathbf{E}_\phi^y(\eta, \mathbf{X}), \mathbf{E}_\phi^y'(\eta, \mathbf{X}), \dots$ entering in the series expansion (3.22) of the element strain. These coefficients depend on both $\eta$ and $\mathbf{X}$ when the lattice has non-uniform geometrical properties in space, as happens for instance with the curved lattice shown in Figure 2.1b.
**Remark 3.4.** By (3.10), the first slots in $l$ correspond to a macroscopic rigid-body displacement (translation $u$ and rotation $\gamma$). The corresponding subblocks in $E_{\varphi}^l(\eta, \mathbf{X})$, $E_{\varphi}^{l'}(\eta, \mathbf{X})$ and $E_{\varphi}^{l''}(\eta, \mathbf{X})$ are actually zero, reflecting the fact that the lattice elements (such as beams in the present case) are unaffected by rigid-body displacement:

* because of the zero blocks in (3.14)$_{2,3}$, the subblocks in $E_{\varphi}^{l'}(\eta, \mathbf{X})$ and $E_{\varphi}^{l''}(\eta, \mathbf{X})$ corresponding to $u$ and $\gamma$ are zero, see (3.23);

* in view of (3.14) and (3.23), the $u$ subblock in $E_{\varphi}^l(\eta, \mathbf{X})$ is $\mathcal{D}_{\varphi}(\eta, \mathbf{X}) : (\mathbf{I}_d \quad \mathbf{I}_d)$, which is zero when the expression of $\mathcal{D}_{\varphi}$ in (2.15) is used (invariance of the element strain $E_{\alpha}$ by a rigid translation);

* the $\gamma$ subblock in $E_{\varphi}^l(\eta, \mathbf{X})$ is $\mathcal{D}_{\varphi}(\eta, \mathbf{X}) : (-\mathcal{N}^{T132} \cdot \Delta_{\varphi}^{-} - \mathcal{N}^{T132} \cdot \Delta_{\varphi}^{+})$, which can be shown to be zero by a similar argument (invariance of the element strain $E_{\alpha}$ by a rigid rotation).

The only point of the macroscopic degrees of freedom $l$ being to compute the element strain by Equation (3.22), it is therefore possible to discard entirely the $u$ and $\gamma$ components in $l$, i.e., to amend (3.10) to

$$ l(\mathbf{X}) = (\tilde{\mathcal{E}}(\mathbf{X})) \quad \text{(alternate, reduced form of } l). \qquad (3.24) $$

In this case the length of $n_l$ in (3.11) must obviously be changed to $n_l=n_{\epsilon}$. The library is flexible: the definition of $l$ can be selected using an option and we check that the expected subblocks are zero in any case.

**Remark 3.5.** In Equation (3.24), the macroscopic degrees of freedom $l$ coincides with the macroscopic strain $\epsilon$. We prefer, however, to treat these notions as distinct, as we will occasionally include additional degrees of freedom into $l$, such as the pre-strain magnitude in Equation (7.9). In any case, $l$ collects all the variables on which the homogenized energy depends, so that the energy is always of the form (4.5).

## 3.6. Assembly

Next, we concatenate the elements' strain measures $E_{\varphi}$ corresponding to the different element families $1 \le \varphi \le n_{\varphi}$, along the left-hand side of the constraint (3.5) into a global strain vector $E(\mathbf{X})$,

$$ E(\eta, \mathbf{X}) = (E_I(\eta, \mathbf{X}) \quad \dots \quad E_{\varphi}(\eta, \mathbf{X}) \quad \dots \quad E_{n_{\varphi}}(\eta, \mathbf{X}) \quad \frac{1}{\eta} \sum_{b=1}^{n_b} \xi_b(\mathbf{X})) \in \mathbb{R}^{n_{\mathcal{E}}}, \qquad (3.25) $$

whose dimension is

$$ n_E = \left( \sum_{\varphi=I}^{n_{\varphi}} n_{E_{\varphi}} \right) + d. \qquad (3.26) $$

By (3.1) and (3.6), it scales as

$$ E(\eta, \mathbf{X}) = O(\eta^0). \qquad (3.27) $$

The honeycomb lattice shown in Figure 2.1, for instance, has $E = (\epsilon_I \quad \kappa_I \quad \epsilon_{II} \quad \kappa_{II} \quad \epsilon_{III} \quad \kappa_{III} \quad \tau_{II} \quad \tau_{III} \quad \frac{\xi_1+\xi_2}{\eta} \cdot \mathbf{e}_1 \quad \frac{\xi_1+\xi_2}{\eta} \cdot \mathbf{e}_2)$ and $n_E = 11$.

The global strain vector $E(\mathbf{X})$ plays a central role in the formulation of the homogenized problem:

* the strain $E_{\alpha}$ in a particular element can be found by evaluating $E(\eta, \mathbf{X})$ at the element center $\mathbf{X}=\mathbf{X}_{\alpha}^c$, and extracting the range of indices corresponding to element family $\varphi_{\alpha}$, see (3.25) and (3.21),

$$ E_{\alpha} = (\mathbf{0}_{n_{E}(\varphi_{\alpha})} \times p_{\varphi} \quad \mathbf{I}_{n_{E}(\varphi_{\alpha})} \quad \mathbf{0}_{n_{E}(\varphi_{\alpha})} \times p'_{\varphi}) \cdot E(\eta, \mathbf{X}_{\alpha}^c) \qquad (3.28) $$

where $p_{\varphi} = \sum_{\varphi'=I}^{\varphi-1} n_{E_{\varphi'}}$ is the number of slots before the $E_{\varphi}$ block in (3.25) and $p'_{\varphi} = \sum_{\varphi'=\varphi+1}^{n_{\varphi}} n_{E_{\varphi'}} + d$ is the number of slots after the $E_{\varphi}$ block;

* in view of (3.28), the elastic energy in (2.20–2.23) can be rewritten as

$$ \Phi_d = \int_{\Omega} \frac{1}{2} \mathcal{E}(\eta, \mathbf{X}) \cdot \mathcal{H}_d(\eta, \mathbf{X}) \cdot \mathcal{E}(\eta, \mathbf{X}) \, d\mathbf{X} \qquad (3.29) $$

where $\mathcal{H}_d(\eta, \mathbf{X})$ is the global stiffness matrix,

$$ \mathcal{H}_d(\eta, \mathbf{X}) = \sum_{\varphi=I}^{n_{\varphi}} \left( \sum_{\substack{\alpha \in \varphi \\ \mathbf{X}_{\alpha}^c \in \Omega}} \delta_D(\mathbf{X} - \mathbf{X}_{\alpha}^c) \right) \mathcal{H}_{\varphi}^g(\eta, \mathbf{X}) \in \mathbb{T}^{(n_E, n_E)}, \qquad (3.30) $$

which is symmetric, $\delta_D$ denotes the Dirac distribution, and $\mathcal{H}_{\varphi}^g(\eta, \mathbf{X})$ represents the contribution from the element family $\varphi$ to the global stiffness matrix,

$$ \mathcal{H}_{\varphi}^{g}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{0}_{p_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p_{\varphi} \times p'_{\varphi}} \\ \mathbf{0}_{n_{E_{\varphi}} \times p_{\varphi}} & \mathcal{H}_{\varphi}(\eta, \mathbf{X}) & \mathbf{0}_{n_{E_{\varphi}} \times p'_{\varphi}} \\ \mathbf{0}_{p'_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p'_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p'_{\varphi} \times p'_{\varphi}} \end{pmatrix} \in \mathbb{T}^{(n_E, n_E)}; \qquad (3.31) $$

* the $n_c=d=2$ kinematic constraints in (3.5) can be rewritten as

$$ \mathcal{Q} \cdot \mathcal{E}(\eta, \mathbf{X}) = \mathbf{0} \quad \forall \mathbf{X} \in \Omega \qquad (3.32) $$
**Remark 3.4.** By (3.10), the first slots in $l$ correspond to a macroscopic rigid-body displacement (translation $u$ and rotation $\gamma$). The corresponding subblocks in $E_{\varphi}^l(\eta, \mathbf{X})$, $E_{\varphi}^{l'}(\eta, \mathbf{X})$ and $E_{\varphi}^{l''}(\eta, \mathbf{X})$ are actually zero, reflecting the fact that the lattice elements (such as beams in the present case) are unaffected by rigid-body displacement:

* because of the zero blocks in (3.14)$_{2,3}$, the subblocks in $E_{\varphi}^{l'}(\eta, \mathbf{X})$ and $E_{\varphi}^{l''}(\eta, \mathbf{X})$ corresponding to $u$ and $\gamma$ are zero, see (3.23);

* in view of (3.14) and (3.23), the $u$ subblock in $E_{\varphi}^l(\eta, \mathbf{X})$ is $\mathcal{D}_{\varphi}(\eta, \mathbf{X}) : (\mathbf{I}_d \quad \mathbf{I}_d)$, which is zero when the expression of $\mathcal{D}_{\varphi}$ in (2.15) is used (invariance of the element strain $E_{\alpha}$ by a rigid translation);

* the $\gamma$ subblock in $E_{\varphi}^l(\eta, \mathbf{X})$ is $\mathcal{D}_{\varphi}(\eta, \mathbf{X}) : (-\mathcal{N}^{T132} \cdot \Delta_{\varphi}^{-} - \mathcal{N}^{T132} \cdot \Delta_{\varphi}^{+})$, which can be shown to be zero by a similar argument (invariance of the element strain $E_{\alpha}$ by a rigid rotation).

The only point of the macroscopic degrees of freedom $l$ being to compute the element strain by Equation (3.22), it is therefore possible to discard entirely the $u$ and $\gamma$ components in $l$, i.e., to amend (3.10) to

$$ l(\mathbf{X}) = (\tilde{\mathcal{E}}(\mathbf{X})) \quad \text{(alternate, reduced form of } l). \qquad (3.24) $$

In this case the length of $n_l$ in (3.11) must obviously be changed to $n_l=n_{\epsilon}$. The library is flexible: the definition of $l$ can be selected using an option and we check that the expected subblocks are zero in any case.

**Remark 3.5.** In Equation (3.24), the macroscopic degrees of freedom $l$ coincides with the macroscopic strain $\epsilon$. We prefer, however, to treat these notions as distinct, as we will occasionally include additional degrees of freedom into $l$, such as the pre-strain magnitude in Equation (7.9). In any case, $l$ collects all the variables on which the homogenized energy depends, so that the energy is always of the form (4.5).

## 3.6. Assembly

Next, we concatenate the elements' strain measures $E_{\varphi}$ corresponding to the different element families $1 \le \varphi \le n_{\varphi}$, along the left-hand side of the constraint (3.5) into a global strain vector $E(\mathbf{X})$,

$$ E(\eta, \mathbf{X}) = (E_I(\eta, \mathbf{X}) \quad \dots \quad E_{\varphi}(\eta, \mathbf{X}) \quad \dots \quad E_{n_{\varphi}}(\eta, \mathbf{X}) \quad \frac{1}{\eta} \sum_{b=1}^{n_b} \xi_b(\mathbf{X})) \in \mathbb{R}^{n_{\mathcal{E}}}, \qquad (3.25) $$

whose dimension is

$$ n_E = \left( \sum_{\varphi=I}^{n_{\varphi}} n_{E_{\varphi}} \right) + d. \qquad (3.26) $$

By (3.1) and (3.6), it scales as

$$ E(\eta, \mathbf{X}) = O(\eta^0). \qquad (3.27) $$

The honeycomb lattice shown in Figure 2.1, for instance, has $E = (\epsilon_I \quad \kappa_I \quad \epsilon_{II} \quad \kappa_{II} \quad \epsilon_{III} \quad \kappa_{III} \quad \tau_{II} \quad \tau_{III} \quad \frac{\xi_1+\xi_2}{\eta} \cdot \mathbf{e}_1 \quad \frac{\xi_1+\xi_2}{\eta} \cdot \mathbf{e}_2)$ and $n_E = 11$.

The global strain vector $E(\mathbf{X})$ plays a central role in the formulation of the homogenized problem:

* the strain $E_{\alpha}$ in a particular element can be found by evaluating $E(\eta, \mathbf{X})$ at the element center $\mathbf{X}=\mathbf{X}_{\alpha}^c$, and extracting the range of indices corresponding to element family $\varphi_{\alpha}$, see (3.25) and (3.21),

$$ E_{\alpha} = (\mathbf{0}_{n_{E}(\varphi_{\alpha})} \times p_{\varphi} \quad \mathbf{I}_{n_{E}(\varphi_{\alpha})} \quad \mathbf{0}_{n_{E}(\varphi_{\alpha})} \times p'_{\varphi}) \cdot E(\eta, \mathbf{X}_{\alpha}^c) \qquad (3.28) $$

where $p_{\varphi} = \sum_{\varphi'=I}^{\varphi-1} n_{E_{\varphi'}}$ is the number of slots before the $E_{\varphi}$ block in (3.25) and $p'_{\varphi} = \sum_{\varphi'=\varphi+1}^{n_{\varphi}} n_{E_{\varphi'}} + d$ is the number of slots after the $E_{\varphi}$ block;

* in view of (3.28), the elastic energy in (2.20–2.23) can be rewritten as

$$ \Phi_d = \int_{\Omega} \frac{1}{2} \mathcal{E}(\eta, \mathbf{X}) \cdot \mathcal{H}_d(\eta, \mathbf{X}) \cdot \mathcal{E}(\eta, \mathbf{X}) \, d\mathbf{X} \qquad (3.29) $$

where $\mathcal{H}_d(\eta, \mathbf{X})$ is the global stiffness matrix,

$$ \mathcal{H}_d(\eta, \mathbf{X}) = \sum_{\varphi=I}^{n_{\varphi}} \left( \sum_{\substack{\alpha \in \varphi \\ \mathbf{X}_{\alpha}^c \in \Omega}} \delta_D(\mathbf{X} - \mathbf{X}_{\alpha}^c) \right) \mathcal{H}_{\varphi}^g(\eta, \mathbf{X}) \in \mathbb{T}^{(n_E, n_E)}, \qquad (3.30) $$

which is symmetric, $\delta_D$ denotes the Dirac distribution, and $\mathcal{H}_{\varphi}^g(\eta, \mathbf{X})$ represents the contribution from the element family $\varphi$ to the global stiffness matrix,

$$ \mathcal{H}_{\varphi}^{g}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{0}_{p_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p_{\varphi} \times p'_{\varphi}} \\ \mathbf{0}_{n_{E_{\varphi}} \times p_{\varphi}} & \mathcal{H}_{\varphi}(\eta, \mathbf{X}) & \mathbf{0}_{n_{E_{\varphi}} \times p'_{\varphi}} \\ \mathbf{0}_{p'_{\varphi} \times p_{\varphi}} & \mathbf{0}_{p'_{\varphi} \times n_{E_{\varphi}}} & \mathbf{0}_{p'_{\varphi} \times p'_{\varphi}} \end{pmatrix} \in \mathbb{T}^{(n_E, n_E)}; \qquad (3.31) $$

* the $n_c=d=2$ kinematic constraints in (3.5) can be rewritten as

$$ \mathcal{Q} \cdot \mathcal{E}(\eta, \mathbf{X}) = \mathbf{0} \quad \forall \mathbf{X} \in \Omega \qquad (3.32) $$
where the constraint extraction matrix $\mathcal{Q}$ extracts the trailing $n_c$ in the expression (3.25) of $\mathbf{E}$,

$$ \mathcal{Q} = ( \mathbf{0}_{(n_E-n_c) \times n_c} \mathbf{I}_{n_c} ) \in \mathbb{T}^{(n_c, n_E)}; \qquad (3.33) $$

• by (3.22), (3.25) the assembled strain vector can be expressed in terms of $\mathbf{l}, \mathbf{y}$ and their gradients as

$$ \begin{aligned} \mathbf{E}(\eta, \mathbf{X}) ={}& \mathbf{E}_l(\eta, \mathbf{X}) \cdot \mathbf{l}(\mathbf{X}) + \mathbf{E}'_l(\eta, \mathbf{X}) : \nabla \mathbf{l}(\mathbf{X}) + \mathbf{E}''_l(\eta, \mathbf{X}) : \nabla^2 \mathbf{l}(\mathbf{X}) + \dots \\ & + \mathbf{E}_y(\eta, \mathbf{X}) \cdot \mathbf{y}(\mathbf{X}) + \mathbf{E}'_y(\eta, \mathbf{X}) : \nabla \mathbf{y}(\mathbf{X}) + \mathbf{E}''_y(\eta, \mathbf{X}) : \nabla^2 \mathbf{y}(\mathbf{X}) + \dots \end{aligned} \qquad (3.34) $$

where the coefficients $\mathbf{E}_l, \mathbf{E}'_l, \mathbf{E}''_l, \mathbf{E}_y, \mathbf{E}'_y$ and $\mathbf{E}''_y$ are obtained by a simple assembly process, see (3.25) and (3.15), as

$$ \mathbf{E}_l^{(i)}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{E}_l^{(i)}(\eta, \mathbf{X}) \\ \vdots \\ \mathbf{E}_{n_\psi}^{(i)}(\eta, \mathbf{X}) \\ \mathbf{0}_{n_c \times n_l \times d^i} \end{pmatrix}, \quad \mathbf{E}_y^{(i)}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{E}_l^{y(i)}(\eta, \mathbf{X}) \\ \vdots \\ \mathbf{E}_{n_\psi}^{y(i)}(\eta, \mathbf{X}) \\ \begin{cases} I_d \mathbf{0}_{d \times n_r} \dots I_d \mathbf{0}_{d \times n_r} & \text{if } i=0 \\ \mathbf{0}_{n_c \times n_l \times d^i} & \text{if } i>0 \end{cases} \end{pmatrix}. \qquad (3.35) $$

**Remark 3.6.** Both $\mathbf{l}(\mathbf{X})$ and $\mathbf{y}(\mathbf{X})$ are of order $\eta^0$ by (3.12) and (3.17), and so is $\mathbf{E}(\eta, \mathbf{X})$ as well by (3.27). As a result, the leading-order coefficients in (3.34) are of order $\mathbf{E}_l = \mathcal{O}(\eta^0)$ and $\mathbf{E}_y = \mathcal{O}(\eta^0)$. A more detailed analysis shows that

$$ \begin{aligned} \mathbf{E}_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^0), & \mathbf{E}'_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^1), & \mathbf{E}''_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^2), \\ \mathbf{E}_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^0), & \mathbf{E}'_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^1), & \mathbf{E}''_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^2). \end{aligned} \qquad (3.36) $$

The increasing powers of $\eta$ from left to right in this table can be traced back to the presence of increasing powers of $\Delta = \mathcal{O}(\eta)$ in the successive $\mathcal{T}$ operators appearing in (3.14) and (3.20), see (2.6).

## 3.7. Continualized strain energy

The energy $\Phi_d$ has been artificially rewritten as an integral in (3.29) but, because of the Dirac weights in (3.30), it is in fact a discrete sum over the element centers. By the Euler–MacLaurin formula [LP16], this discrete sum can be approximated to any desired order in $\eta$ by a ‘genuine’ integral plus boundary terms—we work with an error of order $\eta^3$ in the present work,

$$ \Phi_d = \Phi + \Phi_{bt} + \mathcal{O}(\eta^3) \qquad (3.37) $$

where $\Phi_{bt} = \int_{\partial\Omega} \dots da$ is the boundary term, and $\Phi$ is the continualized energy,

$$ \Phi = \int_{\Omega} \frac{1}{2} \mathbf{E}(\eta, \mathbf{X}) \cdot \mathcal{H}(\eta, \mathbf{X}) \cdot \mathbf{E}(\eta, \mathbf{X}) \, d\mathbf{X}, \qquad (3.38) $$

obtained by taking formally the limit $\eta \to 0$ in (3.29–3.30),

$$ \mathcal{H}(\eta, \mathbf{X}) = \sum_{\varphi=I}^{n_\varphi} \rho_{\varphi}(\eta, \mathbf{X}) \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, \mathbf{X}). \qquad (3.39) $$

Here, $\rho_{\varphi}(\eta, \mathbf{X})$ is the density of element centers for family $\varphi$ in the reference configuration space $\mathbf{X}$, see Figure 2.1b. It is given in terms of the known density $\tilde{\rho}_{\varphi}$ of element centers in the topological lattice space $\tilde{\mathbf{X}}$ (see Table 2.1) as

$$ \rho_{\varphi}(\eta, \mathbf{X}) = \frac{\tilde{\rho}_{\varphi}}{\eta^d |J_f(f^{-1}(\mathbf{X}))|}, \qquad (3.40) $$

where $J_f = \det \nabla f$ is the Jacobian of the diffeomorphism (2.3). Like $\tilde{\rho}_{\varphi}$, the density $\rho_{\varphi}$ is expressed per unit area if $d=2$, and per unit volume if $d=3$.

The proof of Equations (3.37–3.40) can be outlined as follows:

1. the energy in (3.29) is a sum by (3.30) of terms concerning a particular family $\varphi$;

2. each of the terms writes as $\sum_{\alpha \in \varphi} \mathbf{X}_\alpha^c \in \Omega \frac{1}{2} \mathbf{E}(\eta, \mathbf{X}_\alpha^c) \cdot \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, \mathbf{X}_\alpha^c) \cdot \mathbf{E}(\eta, \mathbf{X}_\alpha^c)$ and is continualized separately;

3. with the help of the diffeomorphism, each term is rewritten as a sum over the (periodic) underlying topological lattice, namely $\sum_{\alpha \in \varphi} \tilde{\mathbf{X}}_\alpha^c \in \eta^{-1}f^{-1}(\Omega) \frac{1}{2} \mathbf{E}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c)) \cdot \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c)) \cdot \mathbf{E}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c))$;

4. this is now a sum over a large domain in a periodic lattice, whose term is given by a slowly-variable function: by the Euler–MacLaurin formula, see Theorem 2.1(ii) by Le Floch and Pelayo [LP16], it can be rewritten as an integral plus boundary terms—this proof is constructive and explicit expressions are available:

• the integral term thus obtained corresponds to taking the formal limit $\eta \to 0$ in (3.29–3.30) where the Dirac weights are replaced by their density;

• the boundary terms are available as well but are not needed in the present paper;
where the constraint extraction matrix $\mathcal{Q}$ extracts the trailing $n_c$ in the expression (3.25) of $\mathbf{E}$,

$$ \mathcal{Q} = ( \mathbf{0}_{(n_E-n_c) \times n_c} \mathbf{I}_{n_c} ) \in \mathbb{T}^{(n_c, n_E)}; \qquad (3.33) $$

• by (3.22), (3.25) the assembled strain vector can be expressed in terms of $\mathbf{l}, \mathbf{y}$ and their gradients as

$$ \begin{aligned} \mathbf{E}(\eta, \mathbf{X}) ={}& \mathbf{E}_l(\eta, \mathbf{X}) \cdot \mathbf{l}(\mathbf{X}) + \mathbf{E}'_l(\eta, \mathbf{X}) : \nabla \mathbf{l}(\mathbf{X}) + \mathbf{E}''_l(\eta, \mathbf{X}) : \nabla^2 \mathbf{l}(\mathbf{X}) + \dots \\ & + \mathbf{E}_y(\eta, \mathbf{X}) \cdot \mathbf{y}(\mathbf{X}) + \mathbf{E}'_y(\eta, \mathbf{X}) : \nabla \mathbf{y}(\mathbf{X}) + \mathbf{E}''_y(\eta, \mathbf{X}) : \nabla^2 \mathbf{y}(\mathbf{X}) + \dots \end{aligned} \qquad (3.34) $$

where the coefficients $\mathbf{E}_l, \mathbf{E}'_l, \mathbf{E}''_l, \mathbf{E}_y, \mathbf{E}'_y$ and $\mathbf{E}''_y$ are obtained by a simple assembly process, see (3.25) and (3.15), as

$$ \mathbf{E}_l^{(i)}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{E}_l^{(i)}(\eta, \mathbf{X}) \\ \vdots \\ \mathbf{E}_{n_{\psi}}^{(i)}(\eta, \mathbf{X}) \\ \mathbf{0}_{n_c \times n_l \times d^i} \end{pmatrix}, \quad \mathbf{E}_y^{(i)}(\eta, \mathbf{X}) = \begin{pmatrix} \mathbf{E}_l^{y(i)}(\eta, \mathbf{X}) \\ \vdots \\ \mathbf{E}_{n_{\psi}}^{y(i)}(\eta, \mathbf{X}) \\ \begin{cases} I_d & \mathbf{0}_{d \times n_r} \dots I_d \mathbf{0}_{d \times n_r} & \text{if } i=0 \\ \mathbf{0}_{n_c \times n_l \times d^i} & \mathbf{0}_{d \times n_r} & \text{if } i>0 \end{cases} \end{pmatrix}. \qquad (3.35) $$

**Remark 3.6.** Both $\mathbf{l}(\mathbf{X})$ and $\mathbf{y}(\mathbf{X})$ are of order $\eta^0$ by (3.12) and (3.17), and so is $\mathbf{E}(\eta, \mathbf{X})$ as well by (3.27). As a result, the leading-order coefficients in (3.34) are of order $\mathbf{E}_l = \mathcal{O}(\eta^0)$ and $\mathbf{E}_y = \mathcal{O}(\eta^0)$. A more detailed analysis shows that

$$ \begin{aligned} \mathbf{E}_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^0), & \mathbf{E}'_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^1), & \mathbf{E}''_l(\eta, \mathbf{X}) &= \mathcal{O}(\eta^2), \\ \mathbf{E}_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^0), & \mathbf{E}'_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^1), & \mathbf{E}''_y(\eta, \mathbf{X}) &= \mathcal{O}(\eta^2). \end{aligned} \qquad (3.36) $$

The increasing powers of $\eta$ from left to right in this table can be traced back to the presence of increasing powers of $\Delta = \mathcal{O}(\eta)$ in the successive $\mathcal{T}$ operators appearing in (3.14) and (3.20), see (2.6).

## 3.7. Continualized strain energy

The energy $\Phi_d$ has been artificially rewritten as an integral in (3.29) but, because of the Dirac weights in (3.30), it is in fact a discrete sum over the element centers. By the Euler–MacLaurin formula [LP16], this discrete sum can be approximated to any desired order in $\eta$ by a ‘genuine’ integral plus boundary terms—we work with an error of order $\eta^3$ in the present work,

$$ \Phi_d = \Phi + \Phi_{bt} + \mathcal{O}(\eta^3) \qquad (3.37) $$

where $\Phi_{bt} = \int_{\partial\Omega} \dots da$ is the boundary term, and $\Phi$ is the continualized energy,

$$ \Phi = \int_{\Omega} \frac{1}{2} \mathbf{E}(\eta, \mathbf{X}) \cdot \mathcal{H}(\eta, \mathbf{X}) \cdot \mathbf{E}(\eta, \mathbf{X}) \, d\mathbf{X}, \qquad (3.38) $$

obtained by taking formally the limit $\eta \to 0$ in (3.29–3.30),

$$ \mathcal{H}(\eta, \mathbf{X}) = \sum_{\varphi=I}^{n_{\varphi}} \rho_{\varphi}(\eta, \mathbf{X}) \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, \mathbf{X}). \qquad (3.39) $$

Here, $\rho_{\varphi}(\eta, \mathbf{X})$ is the density of element centers for family $\varphi$ in the reference configuration space $\mathbf{X}$, see Figure 2.1b. It is given in terms of the known density $\tilde{\rho}_{\varphi}$ of element centers in the topological lattice space $\tilde{\mathbf{X}}$ (see Table 2.1) as

$$ \rho_{\varphi}(\eta, \mathbf{X}) = \frac{\tilde{\rho}_{\varphi}}{\eta^d |J_f(f^{-1}(\mathbf{X}))|}, \qquad (3.40) $$

where $J_f = \det \nabla f$ is the Jacobian of the diffeomorphism (2.3). Like $\tilde{\rho}_{\varphi}$, the density $\rho_{\varphi}$ is expressed per unit area if $d=2$, and per unit volume if $d=3$.

The proof of Equations (3.37–3.40) can be outlined as follows:

1. the energy in (3.29) is a sum by (3.30) of terms concerning a particular family $\varphi$;

2. each of the terms writes as $\sum_{\alpha \in \varphi} \mathbf{X}_\alpha^c \in \Omega \frac{1}{2} \mathbf{E}(\eta, \mathbf{X}_\alpha^c) \cdot \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, \mathbf{X}_\alpha^c) \cdot \mathbf{E}(\eta, \mathbf{X}_\alpha^c)$ and is continualized separately;

3. with the help of the diffeomorphism, each term is rewritten as a sum over the (periodic) underlying topological lattice, namely $\sum_{\alpha \in \varphi} \tilde{\mathbf{X}}_\alpha^c \in \eta^{-1}f^{-1}(\Omega) \frac{1}{2} \mathbf{E}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c)) \cdot \mathcal{H}_{\varphi}^{\mathbf{g}}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c)) \cdot \mathbf{E}(\eta, f(\eta \tilde{\mathbf{X}}_\alpha^c))$;

4. this is now a sum over a large domain in a periodic lattice, whose term is given by a slowly-variable function: by the Euler–MacLaurin formula, see Theorem 2.1(ii) by Le Floch and Pelayo [LP16], it can be rewritten as an integral plus boundary terms—this proof is constructive and explicit expressions are available:

• the integral term thus obtained corresponds to taking the formal limit $\eta \to 0$ in (3.29–3.30) where the Dirac weights are replaced by their density;

• the boundary terms are available as well but are not needed in the present paper;
5. the proof of Equations (3.37–3.40) follows by summing up the integrals and boundary terms coming from the different families.

This yields the integrand $\frac{1}{2}\mathbf{E}(\eta, \mathbf{X}) \cdot \mathcal{H}(\eta, \mathbf{X}) \cdot \mathbf{E}(\eta, \mathbf{X})$ announced in (3.38). We do not need the exact expression of the boundary term $\Phi_{bt}$ in the present paper, and therefore leave to future work both the details of the proof and the expression of $\Phi_{bt}$. In this future work, we will combine the boundary term $\Phi_{bt}$ with an analysis of the boundary layer to account for the boundary in an effective way.

Combining (2.21), (3.31), (3.39) and (3.40), we obtain an order-of-magnitude estimate for the continuous stiffness,

$$ \mathcal{H}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0). \qquad (3.41) $$

Since $\mathbf{E}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0)$ by (3.27), the continualized energy in (3.38) scales as

$$ \Phi = \mathcal{O}(\eta^0). \qquad (3.42) $$

The above scalings for $\mathcal{H}$ and $\Phi$ are convenient. They motivate the scaling assumption (2.21).

## 3.8. Summary: canonical form

In our previous work, we carried out homogenization of a generic elasticity problem, formulated in a special form called the canonical form, see Section 2.1 in [AL23]. In Section 3 above, we have rewritten the equilibrium of the lattice in canonical form, which will allow us to apply the homogenization procedure readily. The canonical form is nothing but

* the expansion (3.34) of the strain $\mathbf{E}$ in terms of $\mathbf{l}$, $\mathbf{y}$ and their gradients, which we can rewrite as $\mathbf{E} = \mathbf{E}(\eta, \mathbf{X}; \mathbf{l}(\mathbf{X}), \nabla\mathbf{l}(\mathbf{X}), \dots; \mathbf{y}(\mathbf{X}), \nabla\mathbf{y}(\mathbf{X}), \dots)$,

* the expression (3.32) of the kinematical constraint,

* the expression (3.38) of the energy.

We will homogenize the continualized energy $\Phi$, ignoring the boundary terms $\Phi_{bt}$ appearing in (3.37) during homogenization: the latter can only be interpreted in combination with an analysis of the boundary layers that form at edges of the lattice, which we leave for future work.

# 4. SECOND-ORDER HOMOGENIZATION

In this section, we summarize the homogenization procedure described in [AL23] which applies readily to the canonical form of lattice energy $\Phi = \Phi[\eta, \mathbf{l}, \mathbf{y}]$ derived in Section 3. Here, the square brackets refer to a functional dependence on the arguments $\mathbf{l}$ and $\mathbf{y}$.

## 4.1. Principle of the homogenization method

The macroscopic state vector $\mathbf{l}$ is equal to the macroscopic strain $\tilde{\mathbf{\epsilon}}$, see (3.24) and (3.7). However, to make future extensions easier, we treat $\mathbf{l}$ and $\tilde{\mathbf{\epsilon}}$ as distinct entities, viewing $\mathbf{l}$ as the list of degrees of freedom that are fixed during homogenization (which may be more than just the macroscopic strain $\tilde{\mathbf{\epsilon}}$).

The homogenization procedure described in our previous work [AL23] can be summarized as follows. With $\mathbf{l}$ a prescribed function of $\mathbf{X}$, we revisited homogenization by working at the energy level and showed that the microscopic degrees of freedom $\mathbf{y}(\mathbf{X})$ are the solution of a stationary-point problem,

$$ \mathbf{y}^*[\eta, \mathbf{l}] = \underset{\mathbf{y} \text{ such that } \mathcal{Q}\cdot\mathbf{E}=0 \forall \mathbf{X}}{\text{stpt}} \Phi[\eta, \mathbf{l}, \mathbf{y}]. \qquad (4.1) $$

This stationary-point problem concerns interior points of the domain $\Omega$ only, and proceeds order by order in the expansion parameter $\eta$, the gradient effect being captured as a higher-order correction in $\eta$. The canonical form derived in Section 3 provided us with the energy functional $\Phi[\eta, \mathbf{l}, \mathbf{y}]$ which we can now plug into the homogenization engine developed in our previous work [Aud23].

The stationary point $\mathbf{y}^*$ obtained in (4.1) is a functional of $\mathbf{l}$, as indicated by the square brackets. It is also a function of $\mathbf{X}$ whose values are denoted as $\mathbf{y}^*[\eta, \mathbf{l}](\mathbf{X})$. Equation (4.1) effectively slaves the microscopic degrees of freedom $\mathbf{y}$ to the macroscopic state $\mathbf{l}$.

A homogenized energy functional $\Phi^*$ is obtained by inserting the stationary point $\mathbf{y}^*$ into the original energy,

$$ \Phi^*[\eta, \mathbf{l}] = \Phi[\eta, \mathbf{l}, \mathbf{y}^*[\eta, \mathbf{l}]]. \qquad (4.2) $$

The stationary point problem (4.1) is similar to that governing the equilibrium of the lattice, except that:

* only $\mathbf{y}$ is allowed to vary, $\mathbf{l}$ being kept fixed;

* the potential of the external loading $\Psi$ is ignored, only the strain energy $\Phi$ being considered;

* the energy is written in the domain $\Omega$ strictly contained in the physical domain of the lattice, see Figure 2.3, and boundary terms $\Phi_{bc}$ are ignored.
5. the proof of Equations (3.37–3.40) follows by summing up the integrals and boundary terms coming from the different families.

This yields the integrand $\frac{1}{2}\mathbf{E}(\eta, \mathbf{X}) \cdot \mathcal{H}(\eta, \mathbf{X}) \cdot \mathbf{E}(\eta, \mathbf{X})$ announced in (3.38). We do not need the exact expression of the boundary term $\Phi_{bt}$ in the present paper, and therefore leave to future work both the details of the proof and the expression of $\Phi_{bt}$. In this future work, we will combine the boundary term $\Phi_{bt}$ with an analysis of the boundary layer to account for the boundary in an effective way.

Combining (2.21), (3.31), (3.39) and (3.40), we obtain an order-of-magnitude estimate for the continuous stiffness,

$$ \mathcal{H}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0). \qquad (3.41) $$

Since $\mathbf{E}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0)$ by (3.27), the continualized energy in (3.38) scales as

$$ \Phi = \mathcal{O}(\eta^0). \qquad (3.42) $$

The above scalings for $\mathcal{H}$ and $\Phi$ are convenient. They motivate the scaling assumption (2.21).

## 3.8. Summary: canonical form

In our previous work, we carried out homogenization of a generic elasticity problem, formulated in a special form called the canonical form, see Section 2.1 in [AL23]. In Section 3 above, we have rewritten the equilibrium of the lattice in canonical form, which will allow us to apply the homogenization procedure readily. The canonical form is nothing but

* the expansion (3.34) of the strain $\mathbf{E}$ in terms of $\mathbf{l}$, $\mathbf{y}$ and their gradients, which we can rewrite as $\mathbf{E} = \mathbf{E}(\eta, \mathbf{X}; \mathbf{l}(\mathbf{X}), \nabla\mathbf{l}(\mathbf{X}), \dots; \mathbf{y}(\mathbf{X}), \nabla\mathbf{y}(\mathbf{X}), \dots)$,

* the expression (3.32) of the kinematical constraint,

* the expression (3.38) of the energy.

We will homogenize the continualized energy $\Phi$, ignoring the boundary terms $\Phi_{bt}$ appearing in (3.37) during homogenization: the latter can only be interpreted in combination with an analysis of the boundary layers that form at edges of the lattice, which we leave for future work.

# 4. SECOND-ORDER HOMOGENIZATION

In this section, we summarize the homogenization procedure described in [AL23] which applies readily to the canonical form of lattice energy $\Phi = \Phi[\eta, \mathbf{l}, \mathbf{y}]$ derived in Section 3. Here, the square brackets refer to a functional dependence on the arguments $\mathbf{l}$ and $\mathbf{y}$.

## 4.1. Principle of the homogenization method

The macroscopic state vector $\mathbf{l}$ is equal to the macroscopic strain $\tilde{\mathbf{\epsilon}}$, see (3.24) and (3.7). However, to make future extensions easier, we treat $\mathbf{l}$ and $\tilde{\mathbf{\epsilon}}$ as distinct entities, viewing $\mathbf{l}$ as the list of degrees of freedom that are fixed during homogenization (which may be more than just the macroscopic strain $\tilde{\mathbf{\epsilon}}$).

The homogenization procedure described in our previous work [AL23] can be summarized as follows. With $\mathbf{l}$ a prescribed function of $\mathbf{X}$, we revisited homogenization by working at the energy level and showed that the microscopic degrees of freedom $\mathbf{y}(\mathbf{X})$ are the solution of a stationary-point problem,

$$ \mathbf{y}^*[\eta, \mathbf{l}] = \underset{\mathbf{y} \text{ such that } \mathcal{Q}\cdot\mathbf{E}=0 \forall \mathbf{X}}{\text{stpt}} \Phi[\eta, \mathbf{l}, \mathbf{y}]. \qquad (4.1) $$

This stationary-point problem concerns interior points of the domain $\Omega$ only, and proceeds order by order in the expansion parameter $\eta$, the gradient effect being captured as a higher-order correction in $\eta$. The canonical form derived in Section 3 provided us with the energy functional $\Phi[\eta, \mathbf{l}, \mathbf{y}]$ which we can now plug into the homogenization engine developed in our previous work [Aud23].

The stationary point $\mathbf{y}^*$ obtained in (4.1) is a functional of $\mathbf{l}$, as indicated by the square brackets. It is also a function of $\mathbf{X}$ whose values are denoted as $\mathbf{y}^*[\eta, \mathbf{l}](\mathbf{X})$. Equation (4.1) effectively slaves the microscopic degrees of freedom $\mathbf{y}$ to the macroscopic state $\mathbf{l}$.

A homogenized energy functional $\Phi^*$ is obtained by inserting the stationary point $\mathbf{y}^*$ into the original energy,

$$ \Phi^*[\eta, \mathbf{l}] = \Phi[\eta, \mathbf{l}, \mathbf{y}^*[\eta, \mathbf{l}]]. \qquad (4.2) $$

The stationary point problem (4.1) is similar to that governing the equilibrium of the lattice, except that:

* only $\mathbf{y}$ is allowed to vary, $\mathbf{l}$ being kept fixed;

* the potential of the external loading $\Psi$ is ignored, only the strain energy $\Phi$ being considered;

* the energy is written in the domain $\Omega$ strictly contained in the physical domain of the lattice, see Figure 2.3, and boundary terms $\Phi_{bc}$ are ignored.
Once homogenization is complete, the following additional steps should be addressed: the boundary layers taking place outside the domain $\Omega$ must be analyzed, the boundary terms $\Phi_{bc}$ must be obtained, and the total potential energy—including the potential of the external loading $\Psi$—must be optimized with respect to $l$. When this is done, the proper equilibrium of the lattice is recovered, warranting that the homogenization procedure summarized in (4.1–4.2) is asymptotically consistent. This paper focusses on the derivation of the homogenized model $\Phi^*$: its application to the solution of complete structural problems will be treated in future work.

**Remark 4.1.** The motivation for ignoring the potential of the external loading $\Psi$ during homogenization is as follows. It follows from standard assumptions on the magnitude of the load that the potential $\Psi[\eta, l]$ is a function of $l$ but not of $y$: since $l$ is fixed during homogenization, it is then justified to ignore the external potential $\Psi[\eta, l]$ during homogenization, and to restore it in the final step, when the total energy $\Phi^*[\eta, l] + \Psi[\eta, l]$ is relaxed with respect to $l$—see also the discussion in Section 3.1 in our previous work [AL23].

The homogenization method treats the scale separation parameter $\eta \ll 1$ as an expansion parameter and solves the stationary point problem (4.1) order by order in $\eta$: the optimal microscopic displacement $y^*$ is obtained in the form of an expansion,

$$y^*[\eta, l](\mathbf{X}) = \mathbf{Y}(\eta, \mathbf{X}) \cdot l(\mathbf{X}) + \mathbf{Y}'(\eta, \mathbf{X}) : \nabla l(\mathbf{X}) + \mathcal{O}(\eta^2), \quad (4.3)$$

where the localization tensors $\mathbf{Y}$ and $\mathbf{Y}'$ are an output in symbolic form by the homogenization procedure. They scale as

$$\mathbf{Y}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0), \quad \mathbf{Y}'(\eta, \mathbf{X}) = \mathcal{O}(\eta^1). \quad (4.4)$$

implying that the first term in the right-hand side of (4.3) is the leading-order term, whose order of magnitude agrees with (3.17), while the second term is a corrector capturing the gradient effect.

The homogenized energy $\Phi^*$ in (4.2) is obtained in the form of an expansion as well, which can be written in compact form as

$$\Phi^*[l] = \int_{\Omega} \left( \mathbf{K}(\eta, \mathbf{X}) : \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} + \mathbf{A}(\eta, \mathbf{X}) : (l(\mathbf{X}) \otimes \nabla l(\mathbf{X})) + \mathbf{B}(\eta, \mathbf{X}) : \frac{\nabla l(\mathbf{X}) \otimes \nabla l(\mathbf{X})}{2} \right) d\mathbf{X} + \\ \oint_{\partial\Omega} \left[ \mathbf{k}(\eta, \mathbf{X}) : \left( \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} \otimes \mathbf{n}(\mathbf{X}) \right) + \mathbf{a}(\eta, \mathbf{X}) : (l(\mathbf{X}) \otimes \nabla l(\mathbf{X}) \otimes \mathbf{n}(\mathbf{X})) \right] d\mathbf{X} + \mathcal{O}(\eta^3). \quad (4.5)$$

The boundary integral $\oint_{\partial\Omega} \dots d\mathbf{x}$ in the second line takes place over the boundary $\partial\Omega$ of the domain, with $\mathbf{n}$ denoting the unit outward normal and $d\mathbf{x}$ the area (if $d=3$) or the length (if $d=2$) of a boundary element, see Figure 2.3. The special integral sign $\oint$ will be used throughout for boundary integrals. Note that the boundary terms appearing in (4.5) as a result of homogenization are different from those that appeared earlier in (3.37) as a result of continuing the energy: they must be summed up.

The homogenization procedure delivers the tensors characterizing the homogenized energy in symbolic form. Their orders of magnitude are given by

$$\mathbf{K}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0), \quad \mathbf{A}(\eta, \mathbf{X}) = \mathcal{O}(\eta^1), \quad \mathbf{B}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2), \quad \mathbf{k}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2), \quad \mathbf{a}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2). \quad (4.6)$$

**Remark 4.2.** When the microscopic degrees of freedom and energy $\Phi^*$ are truncated to leading-order $\mathcal{O}(\eta^0)$, we obtain the approximations

$$y^*[\eta, l](\mathbf{X}) = \mathbf{Y}(\eta, \mathbf{X}) \cdot l(\mathbf{X}) + \mathcal{O}(\eta^1), \quad \Phi^*[\eta, l] = \int_{\Omega} \mathbf{K}(\eta, \mathbf{X}) : \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} d\mathbf{X} + \mathcal{O}(\eta^1)$$

that correspond to classical homogenization, in which the effect of the gradient $\nabla l$ is ignored [CMR06; DO11; CP12; Dav13; LR13].

**Remark 4.3.** The homogenization method producing Equations (4.3–4.5) has been presented in our previous work [AL23] using different scaling conventions: the microscopic cell size, now $\mathcal{O}(\eta^1)$, was then $\mathcal{O}(\eta^0)$ and the lattice size, now $\mathcal{O}(\eta^0)$, was then $\mathcal{O}(\eta^{-1})$. The two approaches are ultimately equivalent. The scaling relations announced in (4.4) and (4.6) are obtained by inserting the estimates (3.36) for $\mathbf{E}_l, \mathbf{E}'_l$, etc., into the homogenization formula from Appendix C in [AL23]: in Equation [C.9] in this Appendix, for instance, $\mathcal{F}^1 = \mathcal{O}(\mathbf{E}'_l, \mathbf{E}'_l) = \mathcal{O}(\eta^1)$ while $\mathcal{F}^{11} = \mathcal{O}(\mathbf{E}'_l)$ and similarly $\mathcal{F}^2 = \mathcal{O}(\mathbf{E}'_l, \mathbf{E}'_l) = \mathcal{O}(\eta^2)$. The tensor $\mathcal{A}$ in [C.19] capturing the first-order gradient effect is then $\mathcal{O}(\eta)$, like the quantity $\mathbf{A}_0 \sim \mathbf{A}$ that derives from it, see [C.21]. A similar analysis shows that the quantities derived at second order, such as $\mathbf{Y}', \mathbf{B}$, etc., are all $\mathcal{O}(\eta^2)$.

## 4.2. Connecting with the discrete homogenization library

The homogenization of elastic lattices is implemented and distributed as an extension of the general-purpose homogenization library shoal [Aud23] described in Section 4.1. Both components of the library are written in the symbolic calculation language Wolfram Mathematica. The extension is implemented as follows.

In terms of the parameters describing the lattice listed in Table 4.1, we calculate the following quantities in symbolic form using Wolfram Mathematica [21]:

* the tensors $\mathcal{M}', \mathcal{N}', \mathcal{U}', \mathcal{U}'', \mathcal{G}'$ and $\mathcal{G}''$ depending on space dimension $d$, from Appendix B;
Once homogenization is complete, the following additional steps should be addressed: the boundary layers taking place outside the domain $\Omega$ must be analyzed, the boundary terms $\Phi_{bc}$ must be obtained, and the total potential energy—including the potential of the external loading $\Psi$—must be optimized with respect to $l$. When this is done, the proper equilibrium of the lattice is recovered, warranting that the homogenization procedure summarized in (4.1–4.2) is asymptotically consistent. This paper focusses on the derivation of the homogenized model $\Phi^*$: its application to the solution of complete structural problems will be treated in future work.

**Remark 4.1.** The motivation for ignoring the potential of the external loading $\Psi$ during homogenization is as follows. It follows from standard assumptions on the magnitude of the load that the potential $\Psi[\eta, l]$ is a function of $l$ but not of $y$: since $l$ is fixed during homogenization, it is then justified to ignore the external potential $\Psi[\eta, l]$ during homogenization, and to restore it in the final step, when the total energy $\Phi^*[\eta, l] + \Psi[\eta, l]$ is relaxed with respect to $l$—see also the discussion in Section 3.1 in our previous work [AL23].

The homogenization method treats the scale separation parameter $\eta \ll 1$ as an expansion parameter and solves the stationary point problem (4.1) order by order in $\eta$: the optimal microscopic displacement $y^*$ is obtained in the form of an expansion,

$$y^*[\eta, l](\mathbf{X}) = \mathbf{Y}(\eta, \mathbf{X}) \cdot l(\mathbf{X}) + \mathbf{Y}'(\eta, \mathbf{X}) : \nabla l(\mathbf{X}) + \mathcal{O}(\eta^2), \quad (4.3)$$

where the localization tensors $\mathbf{Y}$ and $\mathbf{Y}'$ are an output in symbolic form by the homogenization procedure. They scale as

$$\mathbf{Y}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0), \quad \mathbf{Y}'(\eta, \mathbf{X}) = \mathcal{O}(\eta^1). \quad (4.4)$$

implying that the first term in the right-hand side of (4.3) is the leading-order term, whose order of magnitude agrees with (3.17), while the second term is a corrector capturing the gradient effect.

The homogenized energy $\Phi^*$ in (4.2) is obtained in the form of an expansion as well, which can be written in compact form as

$$\Phi^*[l] = \int_{\Omega} \left( \mathbf{K}(\eta, \mathbf{X}) : \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} + \mathbf{A}(\eta, \mathbf{X}) : (l(\mathbf{X}) \otimes \nabla l(\mathbf{X})) + \mathbf{B}(\eta, \mathbf{X}) : \frac{\nabla l(\mathbf{X}) \otimes \nabla l(\mathbf{X})}{2} \right) d\mathbf{X} + \int_{\partial\Omega} \left[ \mathbf{k}(\eta, \mathbf{X}) : \left( \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} \otimes \mathbf{n}(\mathbf{X}) \right) + \mathbf{a}(\eta, \mathbf{X}) : (l(\mathbf{X}) \otimes \nabla l(\mathbf{X}) \otimes \mathbf{n}(\mathbf{X})) \right] d\mathbf{X} + \mathcal{O}(\eta^3). \quad (4.5)$$

The boundary integral $\oint_{\partial\Omega} \dots da$ in the second line takes place over the boundary $\partial\Omega$ of the domain, with $\mathbf{n}$ denoting the unit outward normal and $da$ the area (if $d=3$) or the length (if $d=2$) of a boundary element, see Figure 2.3. The special integral sign $\oint$ will be used throughout for boundary integrals. Note that the boundary terms appearing in (4.5) as a result of homogenization are different from those that appeared earlier in (3.37) as a result of continuing the energy: they must be summed up.

The homogenization procedure delivers the tensors characterizing the homogenized energy in symbolic form. Their orders of magnitude are given by

$$\mathbf{K}(\eta, \mathbf{X}) = \mathcal{O}(\eta^0), \quad \mathbf{A}(\eta, \mathbf{X}) = \mathcal{O}(\eta^1), \quad \mathbf{B}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2), \quad \mathbf{k}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2), \quad \mathbf{a}(\eta, \mathbf{X}) = \mathcal{O}(\eta^2). \quad (4.6)$$

**Remark 4.2.** When the microscopic degrees of freedom and energy $\Phi^*$ are truncated to leading-order $\mathcal{O}(\eta^0)$, we obtain the approximations

$$y^*[\eta, l](\mathbf{X}) = \mathbf{Y}(\eta, \mathbf{X}) \cdot l(\mathbf{X}) + \mathcal{O}(\eta^1), \quad \Phi^*[\eta, l] = \int_{\Omega} \mathbf{K}(\eta, \mathbf{X}) : \frac{l(\mathbf{X}) \otimes l(\mathbf{X})}{2} d\mathbf{X} + \mathcal{O}(\eta^1)$$

that correspond to classical homogenization, in which the effect of the gradient $\nabla l$ is ignored [CMR06; DO11; CP12; Dav13; LR13].

**Remark 4.3.** The homogenization method producing Equations (4.3–4.5) has been presented in our previous work [AL23] using different scaling conventions: the microscopic cell size, now $\mathcal{O}(\eta^1)$, was then $\mathcal{O}(\eta^0)$ and the lattice size, now $\mathcal{O}(\eta^0)$, was then $\mathcal{O}(\eta^{-1})$. The two approaches are ultimately equivalent. The scaling relations announced in (4.4) and (4.6) are obtained by inserting the estimates (3.36) for $\mathbf{E}_l, \mathbf{E}'_l$, etc., into the homogenization formula from Appendix C in [AL23]: in Equation [C.9] in this Appendix, for instance, $\mathcal{F}^1 = \mathcal{O}(\mathbf{E}'_l, \mathbf{E}'_l) = \mathcal{O}(\eta^1)$ while $\mathcal{F}^{11} = \mathcal{O}(\mathbf{E}'_l) = \mathcal{O}(\eta^2)$ and similarly $\mathcal{F}^2 = \mathcal{O}(\mathbf{E}'_l, \mathbf{E}'_l) = \mathcal{O}(\eta^2)$. The tensor $\mathcal{A}$ in [C.19] capturing the first-order gradient effect is then $\mathcal{O}(\eta)$, like the quantity $\mathbf{A}_0 \sim \mathbf{A}$ that derives from it, see [C.21]. A similar analysis shows that the quantities derived at second order, such as $\mathbf{Y}', \mathbf{B}$, etc., are all $\mathcal{O}(\eta^2)$.

## 4.2. Connecting with the discrete homogenization library

The homogenization of elastic lattices is implemented and distributed as an extension of the general-purpose homogenization library shoal [Aud23] described in Section 4.1. Both components of the library are written in the symbolic calculation language Wolfram Mathematica. The extension is implemented as follows.

In terms of the parameters describing the lattice listed in Table 4.1, we calculate the following quantities in symbolic form using Wolfram Mathematica [21]:

* the tensors $\mathcal{M}', \mathcal{N}', \mathcal{U}', \mathcal{U}'', \mathcal{G}'$ and $\mathcal{G}''$ depending on space dimension $d$, from Appendix B;
* the tensors $\mathcal{T}_{l^{(i)}}$, $\mathcal{T}_{y^{(i)}}^b$ and $\mathcal{D}_{\varphi}$ using (3.14), (3.20) and (2.15), respectively;

* the tensors $\mathbf{E}_{\varphi}^{l^{(i)}}$ and $\mathbf{E}_{\varphi}^{y^{(i)}}$ using (3.23) and then the tensors $\mathbf{E}_{l}^{(i)}$ and $\mathbf{E}_{y}^{(i)}$ using (3.35);

* the element stiffness matrix $\mathcal{H}_{\varphi}$, its contribution $\mathcal{H}_{\varphi}^g$ to the global stiffness and the global stiffness $\mathcal{H}$ itself using (2.19), (3.31) and (3.39), respectively.

<table><caption>**Table 4.1.** Lattice specification used as an input to the homogenization procedure.</caption><tbody><tr><td rowspan="4">underlying topological lattice<br>(see Table 2.1)</td><td>d</td><td>space dimension</td></tr><tr><td>n<sub>b</sub></td><td>number of Bravais sub-lattices</td></tr><tr><td>n<sub>ϕ</sub></td><td>number of element families</td></tr><tr><td>b<sub>ϕ</sub><sup>±</sup></td><td>elements connectivity</td></tr><tr><td rowspan="2">reference configuration</td><td>ρ<sub>ϕ</sub>(η, X)</td><td>density of element family</td></tr><tr><td>Δ<sub>ϕ</sub><sup>±</sup>(η, X)</td><td>endpoint position with respect to element center</td></tr><tr><td rowspan="2">elastic properties</td><td>ẼA<sub>ϕ</sub>(η, X)</td><td>traction modulus</td></tr><tr><td>χ̃<sub>ϕ</sub>(η, X)</td><td>beam aspect-ratio parameter</td></tr><tr><td>variable material parameters</td><td>m</td><td>(see Remark 7.1)</td></tr></tbody></table>

Next, we feed these quantities to the homogenization engine proposed in our previous work. Specifically, we process and wrap up the data provided in Table 4.1 to produce the following list of arguments that are passed to the homogenization engine—we are still considering the generic beam lattice from Figure 2.1 for the purpose of illustration:

* the space dimension $d=2$, the dimension $n_l=3$ of the macroscopic state vector, the number $n_y=2 \times 3=6$ of microscopic degrees of freedom, the dimension $n_E=11$ of the global strain vector and the numbers $n_c=2$ of kinematical constraint, see (3.11), (3.16), (3.26) and (3.33)$_2$;

* the tensorial functions $\mathbf{E}_l(\eta, X)$, $\mathbf{E}'_l(\eta, X)$, $\mathbf{E}''_l(\eta, X)$, $\mathbf{E}_y(\eta, X)$, $\mathbf{E}'_y(\eta, X)$ and $\mathbf{E}''_y(\eta, X)$ of the strain expansion (3.34);

* the constraint extraction matrix $Q$ in (3.33)$_1$;

* the global stiffness matrix $\mathcal{H}(\eta, X)$ in (3.39), in symbolic form;

* the list $m$ of variable material parameters declared by the user is forwarded to the homogenization engine.

This list $m$ declares all the symbols entering in the expressions of $\mathcal{H}$, $\mathbf{E}_l$, $\mathbf{E}'_l$, $\mathbf{E}''_l$, $\mathbf{E}_y$, $\mathbf{E}'_y$ or $\mathbf{E}''_y$, that vary in space. For lattices having homogeneous properties, $m=\{\}$ is set to be an empty list. For lattices having graded properties, one can either work out the explicit dependence of $\mathcal{H}$, $\mathbf{E}_l$, ... on $X$ as in (3.34) and set $m=\{X_1, \dots, X_d\}$, or use the simpler approach discussed in §7.1–7.3.

**Remark 4.4.** The explicit dependence of the various quantities on the scale-separation parameter $\eta$ was implicit in our previous work [AL23] and is now denoted explicitly by the argument $\eta$. As long as the orders of magnitude are preserved, nothing changes when the input to the homogenization procedure depends on $\eta$.

**Remark 4.5.** The homogenization procedure requires that the energy is positive-definite over the subspace of admissible microscopic degrees of freedom, in the sense that $(\mathbf{E}_y(\eta, X) \cdot \mathbf{y}) \cdot \mathcal{H}(\eta, X) \cdot (\mathbf{E}_y(\eta, X) \cdot \mathbf{y}) > 0$ must be strictly positive for any $\mathbf{y}$ such that $\mathbf{Q} \cdot \mathbf{E}_y(\eta, X) \cdot \mathbf{y} = 0$ (see equation (11) in our previous work [AL23]). We rely on the homogenization engine to check this condition. Whenever present, offending microscopic vectors $\mathbf{y}$ are reported as ‘soft modes’ and the engine aborts with an error message. In future work we will cover how soft modes can be dealt with: this requires treating them as macroscopic degrees of freedom, i.e., incorporating them into the vector $l$.

## 5. SYMBOLIC HOMOGENIZATION OF A PERIODIC HONEYCOMB LATTICE

In this section, a first concrete application of the homogenization scheme is presented. Leaving more advanced examples for the following sections, we start with a simple lattice geometry: we address the *periodic* (non-curved) variant of the honeycomb lattice shown in Figure 2.1, limiting attention to the case of *uniform* elastic properties in space. We derive a simple closed-form expression for its second-order homogenized energy, see Equation (5.9) below. Its leading-order is consistent with existing results from the literature but we also include a higher-order contribution which is novel to the best of our knowledge.

### 5.1. Setting up the lattice in the extensible case

The *periodic* honeycomb lattice is generated using the method described in Section 2, with the diffeomorphism set to $f = \mathrm{id}_2$ (identity in the plane). The properties of topological honeycomb lattice listed in Table 2.1 are still used. The resulting lattice is shown in Figure 5.1.
* the tensors $\mathcal{T}_l^{(i)}$, $\mathcal{T}_y^{(i)}$ and $\mathcal{D}_\varphi$ using (3.14), (3.20) and (2.15), respectively;

* the tensors $\mathcal{E}_\varphi^{(i)}$ and $\mathcal{E}_\varphi^{(i)}$ using (3.23) and then the tensors $\mathcal{E}_l^{(i)}$ and $\mathcal{E}_y^{(i)}$ using (3.35);

* the element stiffness matrix $\mathcal{H}_\varphi$, its contribution $\mathcal{H}_\varphi^g$ to the global stiffness and the global stiffness $\mathcal{H}$ itself using (2.19), (3.31) and (3.39), respectively.

<table><caption>**Table 4.1.** Lattice specification used as an input to the homogenization procedure.</caption><tbody><tr><td rowspan="4">underlying topological lattice<br>(see Table 2.1)</td><td>d</td><td>space dimension</td></tr><tr><td>n<sub>b</sub></td><td>number of Bravais sub-lattices</td></tr><tr><td>n<sub>ϕ</sub></td><td>number of element families</td></tr><tr><td>b<sub>ϕ</sub><sup>±</sup></td><td>elements connectivity</td></tr><tr><td rowspan="2">reference configuration</td><td>ρ<sub>ϕ</sub>(η, X)</td><td>density of element family</td></tr><tr><td>Δ<sub>ϕ</sub><sup>±</sup>(η, X)</td><td>endpoint position with respect to element center</td></tr><tr><td rowspan="2">elastic properties</td><td>ẼA<sub>ϕ</sub>(η, X)</td><td>traction modulus</td></tr><tr><td>χ̃<sub>ϕ</sub>(η, X)</td><td>beam aspect-ratio parameter</td></tr><tr><td>variable material parameters</td><td>m</td><td>(see Remark 7.1)</td></tr></tbody></table>

Next, we feed these quantities to the homogenization engine proposed in our previous work. Specifically, we process and wrap up the data provided in Table 4.1 to produce the following list of arguments that are passed to the homogenization engine—we are still considering the generic beam lattice from Figure 2.1 for the purpose of illustration:

* the space dimension $d=2$, the dimension $n_l=3$ of the macroscopic state vector, the number $n_y=2 \times 3 = 6$ of microscopic degrees of freedom, the dimension $n_E=11$ of the global strain vector and the numbers $n_c=2$ of kinematical constraint, see (3.11), (3.16), (3.26) and (3.33)$_2$;

* the tensorial functions $\mathcal{E}_l(\eta, X)$, $\mathcal{E}'_l(\eta, X)$, $\mathcal{E}''_l(\eta, X)$, $\mathcal{E}_y(\eta, X)$, $\mathcal{E}'_y(\eta, X)$ and $\mathcal{E}''_y(\eta, X)$ of the strain expansion (3.34);

* the constraint extraction matrix $Q$ in (3.33)$_1$;

* the global stiffness matrix $\mathcal{H}(\eta, X)$ in (3.39), in symbolic form;

* the list $m$ of variable material parameters declared by the user is forwarded to the homogenization engine.

This list $m$ declares all the symbols entering in the expressions of $\mathcal{H}$, $\mathcal{E}_l$, $\mathcal{E}'_l$, $\mathcal{E}''_l$, $\mathcal{E}_y$, $\mathcal{E}'_y$ or $\mathcal{E}''_y$, that vary in space. For lattices having homogeneous properties, $m=\{\}$ is set to be an empty list. For lattices having graded properties, one can either work out the explicit dependence of $\mathcal{H}$, $\mathcal{E}_l$, ... on $X$ as in (3.34) and set $m=\{X_1, \dots, X_d\}$, or use the simpler approach discussed in §7.1–7.3.

**Remark 4.4.** The explicit dependence of the various quantities on the scale-separation parameter $\eta$ was implicit in our previous work [AL23] and is now denoted explicitly by the argument $\eta$. As long as the orders of magnitude are preserved, nothing changes when the input to the homogenization procedure depends on $\eta$.

**Remark 4.5.** The homogenization procedure requires that the energy is positive-definite over the subspace of admissible microscopic degrees of freedom, in the sense that $(\mathcal{E}_y(\eta, X) \cdot y) \cdot \mathcal{H}(\eta, X) \cdot (\mathcal{E}_y(\eta, X) \cdot y) > 0$ must be strictly positive for any $y$ such that $Q \cdot \mathcal{E}_y(\eta, X) \cdot y = 0$ (see equation (11) in our previous work [AL23]). We rely on the homogenization engine to check this condition. Whenever present, offending microscopic vectors $y$ are reported as ‘soft modes’ and the engine aborts with an error message. In future work we will cover how soft modes can be dealt with: this requires treating them as macroscopic degrees of freedom, i.e., incorporating them into the vector $l$.

# 5. SYMBOLIC HOMOGENIZATION OF A PERIODIC HONEYCOMB LATTICE

In this section, a first concrete application of the homogenization scheme is presented. Leaving more advanced examples for the following sections, we start with a simple lattice geometry: we address the *periodic* (non-curved) variant of the honeycomb lattice shown in Figure 2.1, limiting attention to the case of *uniform* elastic properties in space. We derive a simple closed-form expression for its second-order homogenized energy, see Equation (5.9) below. Its leading-order is consistent with existing results from the literature but we also include a higher-order contribution which is novel to the best of our knowledge.

## 5.1. Setting up the lattice in the extensible case

The *periodic* honeycomb lattice is generated using the method described in Section 2, with the diffeomorphism set to $f = \mathrm{id}_2$ (identity in the plane). The properties of topological honeycomb lattice listed in Table 2.1 are still used. The resulting lattice is shown in Figure 5.1.
<figure><img src="image_4.png" /><figcaption>**Figure 5.1.** Undeformed configuration of the periodic honeycomb lattice having uniform properties, used as an input to the homogenization procedure. For this un-curved lattice, the elements centers (red dots) coincide with the midpoints.</figcaption></figure>

Since a contraction factor $\eta$ is applied to the underlying topological lattice (whose beam length is $a = \mathcal{O}(1)$), the generated lattice has beam length $\ell = \eta a \ll 1$, see also Figure 2.1b. In view of the value of $\tilde{\rho}_{\varphi} = \frac{2}{3\sqrt{3}a^2}$ in Table 2.1 and of Equation (3.40), the element density $\rho_{\varphi}$ is given, for each one of the three element families $\varphi \in \{I, II, III\}$, by

$$ \rho_{\varphi} = \frac{2}{3\sqrt{3}\ell^2}. \qquad (5.1) $$

With $f = \mathbf{id}_2$, the element centers $X_a^c$ coincide with the midpoints (red dots in the figure), see (2.5).

The homogenization of the honeycomb lattice is carried out in the companion Mathematica notebook homogenize-extensible-honeycomb nb which is distributed along with the library. The quantities appearing in Table 4.1 are passed on input:

* the integers $d=2$, $n_b=2$, $n_{\varphi}=3$ and the array $b_{\varphi}^{\pm}$ characterizing the underlying topological lattice, see Table 2.1,

* the reference configuration is characterized by $\rho_{\varphi}$ from (5.1) and by the relative position $\Delta_{\varphi}^{\pm} = \pm \frac{\eta \delta_{\varphi}}{2}$ of the endpoints with respect to the element center (midpoints), where $\delta_{\varphi}$ is available from Table 2.1,

* each one of the $n_{\varphi}=3$ element type is assigned the stiffness matrix $\mathcal{H}_{\varphi}$ in (2.19) applicable to extensible beams, with constant material parameters; indeed, for this particular lattice, the stretching modulus EA, the beam length $\ell = \eta a$ and the aspect ratio parameter $\chi$ are all (i) invariant in space and (ii) identical for all three element families $\varphi \in \{I, II, III\}$,

* the variable material parameters $\mathbf{m} = \{\}$ is set up as an empty list since the lattice is homogeneous.

The parameters EA, $\ell$ and $\chi$ are treated as symbolic parameters by the library, and the homogenized energy is delivered in the form of a function of EA, $\ell$ and $\chi$.

**Remark 5.1.** We noted earlier that the aspect-ratio parameter $\chi \ll 1$ needs to be small for the beam theory to be applicable, see the discussion below Equation (2.17). Mathematically, $\chi$ has to be linked to the small parameter $\eta$, see Remark 2.2. Remarkably, we can carry out homogenization without having to specify how $\chi$ scales with $\eta$: the only assumption we made is $\chi = \mathcal{O}(\eta^0)$ and this does not exclude more specific assumptions such as $\chi = \mathcal{O}(\eta^1)$ or $\chi = \mathcal{O}(\eta^2)$. Delaying the choice of this scaling law until after homogenization is complete enables us to carry out homogenization once for all, and to identify the interesting scaling regimes based on the homogenized energy, see for instance the pantograph application example in Section 7.5. In bending-dominated lattices such as the honeycomb, the inextensible limit can be obtained simply by taking the limit $\chi \to 0$ in the homogenized energy (§5.3).

## 5.2. Homogenization results

This section is a summary of the symbolic homogenization results obtained in the companion Mathematica notebook homogenize-extensible-honeycomb nb. The library yields the homogenized energy $\Phi^*[I]$ in the form (4.5) that makes use of the strain vector $\mathbf{l} = \tilde{\mathbf{e}}$ in Mandel representation, see (3.7) and (3.24). We systematically rewrite it in terms of the standard 2D strain tensor $\boldsymbol{\epsilon}$ with the help of (B.4).

The homogenization results can be summarized as follows, order by order:

* **Leading order.** The equivalent elastic continuum is isotropic,

$$ \begin{aligned} \Phi^* &= \frac{1}{2} \int_{\Omega} \frac{1}{2} \dot{\tilde{\mathbf{e}}} \cdot \mathbf{K} \cdot \dot{\tilde{\mathbf{e}}} d\mathbf{X} + O(\eta) \\ &= \frac{1}{2} \int_{\Omega} (\lambda \operatorname{tr}^2 \boldsymbol{\epsilon} + 2\mu \boldsymbol{\epsilon} : \boldsymbol{\epsilon}) d\mathbf{X} + O(\eta), \end{aligned} \qquad (5.2) $$

and its Lamé parameters are given by

$$ \lambda = \frac{\sqrt{3} EA}{6 \ell} \frac{1-12 \chi}{1+12 \chi}, \qquad \mu = \frac{4 \sqrt{3} EA}{\ell} \frac{\chi}{1+12 \chi}. \qquad (5.3) $$

This matches the previous results of [VDP14; GLTK19]. In the slender limit ($\chi \to 0$), we have $\mu = \mathcal{O}(EA\chi/\ell) = \mathcal{O}(EI/\ell^3) \ll \lambda = \mathcal{O}(EA/\ell)$, from which we conclude that (i) the lattice is bending-dominated and not stretching-dominated (due to $\mu \propto EI$) and (ii) is incompressible in this limit (due to $\lambda \propto EA$, see (5.8) below).
<figure><img src="image_6.png" /><figcaption>**Figure 5.1.** Undeformed configuration of the periodic honeycomb lattice having uniform properties, used as an input to the homogenization procedure. For this un-curved lattice, the elements centers (red dots) coincide with the midpoints.</figcaption></figure>

Since a contraction factor $\eta$ is applied to the underlying topological lattice (whose beam length is $a = \mathcal{O}(1)$), the generated lattice has beam length $\ell = \eta a \ll 1$, see also Figure 2.1b. In view of the value of $\tilde{\rho}_{\varphi} = \frac{2}{3\sqrt{3}a^2}$ in Table 2.1 and of Equation (3.40), the element density $\rho_{\varphi}$ is given, for each one of the three element families $\varphi \in \{I, II, III\}$, by

$$ \rho_{\varphi} = \frac{2}{3\sqrt{3}\ell^2}. \qquad (5.1) $$

With $f = \mathbf{id}_2$, the element centers $X_a^c$ coincide with the midpoints (red dots in the figure), see (2.5).

The homogenization of the honeycomb lattice is carried out in the companion Mathematica notebook homogenize-extensible-honeycomb nb which is distributed along with the library. The quantities appearing in Table 4.1 are passed on input:

* the integers $d=2$, $n_b=2$, $n_{\varphi}=3$ and the array $b_{\varphi}^{\pm}$ characterizing the underlying topological lattice, see Table 2.1,

* the reference configuration is characterized by $\rho_{\varphi}$ from (5.1) and by the relative position $\Delta_{\varphi}^{\pm} = \pm \frac{\eta \delta_{\varphi}}{2}$ of the endpoints with respect to the element center (midpoints), where $\delta_{\varphi}$ is available from Table 2.1,

* each one of the $n_{\varphi}=3$ element type is assigned the stiffness matrix $\mathcal{H}_{\varphi}$ in (2.19) applicable to extensible beams, with constant material parameters; indeed, for this particular lattice, the stretching modulus EA, the beam length $\ell = \eta a$ and the aspect ratio parameter $\chi$ are all (i) invariant in space and (ii) identical for all three element families $\varphi \in \{I, II, III\}$,

* the variable material parameters $\mathbf{m} = \{\}$ is set up as an empty list since the lattice is homogeneous.

The parameters EA, $\ell$ and $\chi$ are treated as symbolic parameters by the library, and the homogenized energy is delivered in the form of a function of EA, $\ell$ and $\chi$.

**Remark 5.1.** We noted earlier that the aspect-ratio parameter $\chi \ll 1$ needs to be small for the beam theory to be applicable, see the discussion below Equation (2.17). Mathematically, $\chi$ has to be linked to the small parameter $\eta$, see Remark 2.2. Remarkably, we can carry out homogenization without having to specify how $\chi$ scales with $\eta$: the only assumption we made is $\chi = \mathcal{O}(\eta^0)$ and this does not exclude more specific assumptions such as $\chi = \mathcal{O}(\eta^1)$ or $\chi = \mathcal{O}(\eta^2)$. Delaying the choice of this scaling law until after homogenization is complete enables us to carry out homogenization once for all, and to identify the interesting scaling regimes based on the homogenized energy, see for instance the pantograph application example in Section 7.5. In bending-dominated lattices such as the honeycomb, the inextensible limit can be obtained simply by taking the limit $\chi \to 0$ in the homogenized energy (§5.3).

## 5.2. Homogenization results

This section is a summary of the symbolic homogenization results obtained in the companion Mathematica notebook homogenize-extensible-honeycomb nb. The library yields the homogenized energy $\Phi^*[I]$ in the form (4.5) that makes use of the strain vector $\mathbf{l} = \tilde{\mathbf{e}}$ in Mandel representation, see (3.7) and (3.24). We systematically rewrite it in terms of the standard 2D strain tensor $\boldsymbol{\varepsilon}$ with the help of (B.4).

The homogenization results can be summarized as follows, order by order:

* **Leading order.** The equivalent elastic continuum is isotropic,

$$ \begin{aligned} \Phi^* &= \frac{1}{2} \int_{\Omega} \frac{1}{2} \dot{\tilde{\mathbf{e}}} \cdot \mathbf{K} \cdot \dot{\tilde{\mathbf{e}}} d\mathbf{X} + O(\eta) \\ &= \frac{1}{2} \int_{\Omega} (\lambda \operatorname{tr}^2 \boldsymbol{\varepsilon} + 2\mu \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon}) d\mathbf{X} + O(\eta), \end{aligned} \qquad (5.2) $$

and its Lamé parameters are given by

$$ \lambda = \frac{\sqrt{3} EA}{6 \ell} \frac{1-12 \chi}{1+12 \chi}, \qquad \mu = \frac{4 \sqrt{3} EA}{\ell} \frac{\chi}{1+12 \chi}. \qquad (5.3) $$

This matches the previous results of [VDP14; GLTK19]. In the slender limit ($\chi \to 0$), we have $\mu = \mathcal{O}(EA\chi/\ell) = \mathcal{O}(EI/\ell^3) \ll \lambda = \mathcal{O}(EA/\ell)$, from which we conclude that (i) the lattice is bending-dominated and not stretching-dominated (due to $\mu \propto EI$) and (ii) is incompressible in this limit (due to $\lambda \propto EA$, see (5.8) below).
The microscopic displacement underlying this leading-order energy prediction (5.2–5.3) is obtained from the localization tensor $\mathbf{Y}$ appearing in (4.3) as

$$ (\xi_b)^*_{[0]} = (-1)^{b+1} \frac{\ell}{4} \frac{1-12\chi}{1+12\chi} \mathcal{J}:\boldsymbol{\varepsilon}, \quad (\psi_b)^*_{[0]} = 0. \qquad (5.4) $$

where $(\xi_b)^*_{[0]}$ and $(\psi_b)^*_{[0]}$ are the leading-order microscopic displacement and rotation, respectively, of the Bravais sub-lattice $b \in \{1, 2\}$, and $\mathcal{J}$ is a constant tensor associated with the $D_6$ symmetry,

$$ \mathcal{J} = e_1 \otimes (-e_1 \otimes e_1 + e_2 \otimes e_2) + e_2 \otimes (e_1 \otimes e_2 + e_2 \otimes e_1). \qquad (5.5) $$

In (5.4), $(\psi_b)^*_{[0]} = 0$ means that the nodal rotation $\theta_b = \gamma$ matches the macroscopic rotation $\gamma = -\frac{1}{2} \nabla \mathbf{u}: \mathcal{N}$ set by the macroscopic displacement $\mathbf{u}(\mathbf{X})$, and not that the node rotation is zero—see (3.4)$_2$.

* **First order.** When pushed to first order, the homogenization procedure delivers the homogenized tensors $\mathbf{A}=0$ and $\mathbf{k}=0$. In view of (4.5), and as usual with centrosymmetric lattices [ABB10], this implies that the energy has no contribution of order $\eta$. As a result, the leading-order approximation (5.2) is exact up to $\mathcal{O}(\eta^2)$ and not just $\mathcal{O}(\eta)$.

* **Second order.** In the Mathematica notebook, both the bulk term $\mathbf{B}_t := \frac{\nabla l \otimes \nabla l}{2}$ and the boundary term $\mathbf{a}_t := (l \otimes \nabla l \otimes n)$ entering at order $\eta^2$ in the energy (4.5) have been worked out in terms of invariants relevant to the $D_6$ symmetry of the hexagonal network: there are 6 invariants for the surface integral, and 5 for the boundary integral. The full expressions are available in the Mathematica notebook but they are cumbersome and are not included here—their values in the inextensible limit $\chi \to 0$ are nevertheless given in Section 5.3.

The first-order correction to the microscopic displacement $(\xi_b)^*_{[1]}$ and $(\psi_b)^*_{[1]} = 0$ that underlies this second-order energy correction is encoded in the localization tensor $\mathbf{Y}'$ appearing in (4.3): interpreting the result, we get

$$ (\xi_b)^*_{[1]} = 0, \quad (\psi_b)^*_{[1]} = (-1)^{b+1} \frac{3\ell}{4(1+12\chi)} \mathrm{curl}(\mathcal{J}:\boldsymbol{\varepsilon}), \qquad (5.6) $$

where $\mathrm{curl} \mathbf{a} = \partial a_2 / \partial X_1 - \partial a_1 / \partial X_2$ denotes the (scalar) curl of a vector field $\mathbf{a}(\mathbf{X})$ in 2D.

## 5.3. Inextensible limit

We now take the slender limit $\chi \to 0$ in the homogenization results from Section 5.2. Specifically, we address the distinguished limit in which the slenderness parameter $\chi$ and the scale separation parameter $\eta$ both go to zero with

$$ (EA) \chi = \mathcal{O}(\eta). \qquad (5.7) $$

The equivalent compression modulus $\lambda \sim \frac{EA}{\ell} = \mathcal{O}(\chi^{-1}) \to \infty$ in (5.3) then becomes infinite while the shear modulus $\mu \sim \frac{\chi EA}{\ell} \sim \mathcal{O}(1)$ remains finite. For better legibility of the results, we eliminate $EA$ in favor of $EI$ everywhere using $EA = EI / (\chi \ell^2)$, see (2.17), anticipating that in the limit, the beam elements are effectively inextensible and the lattice deforms by pure bending.

In this limit $\lambda \to \infty$ in (5.3)$_1$, implying that the equivalent Cauchy medium is incompressible,

$$ \mathrm{tr} \boldsymbol{\varepsilon}(\mathbf{X}) = 0, \quad \forall \mathbf{X}. \qquad (5.8) $$

The homogenized energy is found by combining the shearing term in leading-order correction (5.2–5.3) with the second-order correction described in words in Section 5.2, taking the kinematic constraint of incompressibility (5.8) into account to remove the apparent divergences that arise when $\chi \to 0$. The result is obtained in the companion Mathematica notebook in the form

$$ \Phi^* = \frac{EI}{\ell^3} \left( \frac{1}{2} \int_{\Omega} 2 \times 4 \sqrt{3} \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon} \, d\mathbf{X} \right. \\ \left. + \frac{\ell^2}{2} \int_{\Omega} \sqrt{3} \left( -7 \nabla \boldsymbol{\varepsilon} : \nabla \boldsymbol{\varepsilon} + \frac{22}{3} \| \nabla \boldsymbol{\varepsilon} : \mathbf{I}_2 \|^2 + 2 (\mathrm{div}(\mathcal{J}:\boldsymbol{\varepsilon}))^2 \right) d\mathbf{X} \right. \\ \left. + \ell^2 \oint_{\partial\Omega} \frac{\sqrt{3}}{2} \left( 13 \nabla \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon} - \frac{7}{2} \mathrm{div}(\mathcal{J}:\boldsymbol{\varepsilon}) \mathcal{J}:\boldsymbol{\varepsilon} \right) \cdot \mathbf{n} \, d\mathbf{S} + \mathcal{O}(\eta^3), \qquad (5.9) $$

where $\mathrm{div} \mathbf{a} = \partial a_1 / \partial X_1 + \partial a_2 / \partial X_2$ denotes the divergence of a 2D vector field $\mathbf{a}(\mathbf{X})$. The successive lines in the right-hand side represent the leading-order contribution—of order $EI / \ell^3 \sim EA \chi / \ell = \mathcal{O}(1)$ by (5.7)—, the surface integral for the second-order correction, and the boundary term for the second-order correction, respectively. In (5.9), $\mathbf{I}_2$ denotes the identity matrix in dimension 2 and the constant tensor $\mathcal{J}$ has been defined in (5.5). Due to this $\mathcal{J}$ term, the equivalent strain-gradient continuum defined by $\Phi^*$ is $D_6$-symmetric, like the original honeycomb itself, and the effective isotropy obtained at the leading order is no longer applicable [RA16; RA19].

The analytical, second-order expansion (5.9) of the honeycomb energy is novel to the best of our knowledge.

In the inextensible limit, the microscopic displacement found in (5.4–5.6) becomes

$$ \xi_b^* = \ell \left( \frac{(-1)^{b+1}}{4} \mathcal{J}:\boldsymbol{\varepsilon} \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2) \right), \qquad (5.10) \\ \psi_b^* = 0 \ell^0 - \frac{3(-1)^{b+1}}{4} \mathrm{curl}(\mathcal{J}:\boldsymbol{\varepsilon}) \ell^1 + \mathcal{O}(\ell^2). $$
The microscopic displacement underlying this leading-order energy prediction (5.2–5.3) is obtained from the localization tensor $\mathbf{Y}$ appearing in (4.3) as

$$ (\xi_b)^*_{[0]} = (-1)^{b+1} \frac{\ell}{4} \frac{1-12\chi}{1+12\chi} \mathcal{J}:\boldsymbol{\varepsilon}, \quad (\psi_b)^*_{[0]} = 0. \qquad (5.4) $$

where $(\xi_b)^*_{[0]}$ and $(\psi_b)^*_{[0]}$ are the leading-order microscopic displacement and rotation, respectively, of the Bravais sub-lattice $b \in \{1, 2\}$, and $\mathcal{J}$ is a constant tensor associated with the $D_6$ symmetry,

$$ \mathcal{J} = e_1 \otimes (-e_1 \otimes e_1 + e_2 \otimes e_2) + e_2 \otimes (e_1 \otimes e_2 + e_2 \otimes e_1). \qquad (5.5) $$

In (5.4), $(\psi_b)^*_{[0]} = 0$ means that the nodal rotation $\theta_b = \gamma$ matches the macroscopic rotation $\gamma = -\frac{1}{2} \nabla \mathbf{u}: \mathcal{N}$ set by the macroscopic displacement $\mathbf{u}(\mathbf{X})$, and not that the node rotation is zero—see (3.4)$_2$.

* **First order.** When pushed to first order, the homogenization procedure delivers the homogenized tensors $\mathbf{A}=\mathbf{0}$ and $\mathbf{k}=\mathbf{0}$. In view of (4.5), and as usual with centrosymmetric lattices [ABB10], this implies that the energy has no contribution of order $\eta$. As a result, the leading-order approximation (5.2) is exact up to $\mathcal{O}(\eta^2)$ and not just $\mathcal{O}(\eta)$.

* **Second order.** In the Mathematica notebook, both the bulk term $\mathbf{B}_t := \frac{\nabla l \otimes \nabla l}{2}$ and the boundary term $\mathbf{a}_t := (l \otimes \nabla l \otimes n)$ entering at order $\eta^2$ in the energy (4.5) have been worked out in terms of invariants relevant to the $D_6$ symmetry of the hexagonal network: there are 6 invariants for the surface integral, and 5 for the boundary integral. The full expressions are available in the Mathematica notebook but they are cumbersome and are not included here—their values in the inextensible limit $\chi \to 0$ are nevertheless given in Section 5.3.

The first-order correction to the microscopic displacement $(\xi_b)^*_{[1]}$ and $(\psi_b)^*_{[1]} = \mathbf{0}$ that underlies this second-order energy correction is encoded in the localization tensor $\mathbf{Y}'$ appearing in (4.3): interpreting the result, we get

$$ (\xi_b)^*_{[1]} = \mathbf{0}, \quad (\psi_b)^*_{[1]} = (-1)^{b+1} \frac{3\ell}{4(1+12\chi)} \mathrm{curl}(\mathcal{J}:\boldsymbol{\varepsilon}), \qquad (5.6) $$

where $\mathrm{curl} \mathbf{a} = \partial a_2 / \partial X_1 - \partial a_1 / \partial X_2$ denotes the (scalar) curl of a vector field $\mathbf{a}(\mathbf{X})$ in 2D.

## 5.3. Inextensible limit

We now take the slender limit $\chi \to 0$ in the homogenization results from Section 5.2. Specifically, we address the distinguished limit in which the slenderness parameter $\chi$ and the scale separation parameter $\eta$ both go to zero with

$$ (EA) \quad \chi = \mathcal{O}(\eta). \qquad (5.7) $$

The equivalent compression modulus $\lambda \sim \frac{EA}{\ell} = \mathcal{O}(\chi^{-1}) \to \infty$ in (5.3) then becomes infinite while the shear modulus $\mu \sim \frac{\chi EA}{\ell} \sim \mathcal{O}(1)$ remains finite. For better legibility of the results, we eliminate $EA$ in favor of $EI$ everywhere using $EA = EI / (\chi \ell^2)$, see (2.17), anticipating that in the limit, the beam elements are effectively inextensible and the lattice deforms by pure bending.

In this limit $\lambda \to \infty$ in (5.3)$_1$, implying that the equivalent Cauchy medium is incompressible,

$$ \mathrm{tr} \, \boldsymbol{\varepsilon}(\mathbf{X}) = 0, \quad \forall \mathbf{X}. \qquad (5.8) $$

The homogenized energy is found by combining the shearing term in leading-order correction (5.2–5.3) with the second-order correction described in words in Section 5.2, taking the kinematic constraint of incompressibility (5.8) into account to remove the apparent divergences that arise when $\chi \to 0$. The result is obtained in the companion Mathematica notebook in the form

$$ \Phi^* = \frac{EI}{\ell^3} \left( \frac{1}{2} \int_{\Omega} 2 \times 4 \sqrt{3} \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon} \, d\mathbf{X} \right. \\ \left. + \frac{\ell^2}{2} \int_{\Omega} \sqrt{3} \left( -7 \nabla \boldsymbol{\varepsilon} : \nabla \boldsymbol{\varepsilon} + \frac{22}{3} \| \nabla \boldsymbol{\varepsilon} : \mathbf{I}_2 \|^2 + 2 (\mathrm{div}(\mathcal{J}:\boldsymbol{\varepsilon}))^2 \right) d\mathbf{X} \right. \\ \left. + \ell^2 \oint_{\partial\Omega} \frac{\sqrt{3}}{2} \left( 13 \nabla \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon} - \frac{7}{2} \mathrm{div}(\mathcal{J}:\boldsymbol{\varepsilon}) \mathcal{J}:\boldsymbol{\varepsilon} \right) \cdot \mathbf{n} \, d\mathbf{S} + \mathcal{O}(\eta^3), \right) \qquad (5.9) $$

where $\mathrm{div} \mathbf{a} = \partial a_1 / \partial X_1 + \partial a_2 / \partial X_2$ denotes the divergence of a 2D vector field $\mathbf{a}(\mathbf{X})$. The successive lines in the right-hand side represent the leading-order contribution—of order $EI / \ell^3 \sim EA \chi / \ell = \mathcal{O}(1)$ by (5.7)—, the surface integral for the second-order correction, and the boundary term for the second-order correction, respectively. In (5.9), $\mathbf{I}_2$ denotes the identity matrix in dimension 2 and the constant tensor $\mathcal{J}$ has been defined in (5.5). Due to this $\mathcal{J}$ term, the equivalent strain-gradient continuum defined by $\Phi^*$ is $D_6$-symmetric, like the original honeycomb itself, and the effective isotropy obtained at the leading order is no longer applicable [RA16; RA19].

The analytical, second-order expansion (5.9) of the honeycomb energy is novel to the best of our knowledge.

In the inextensible limit, the microscopic displacement found in (5.4–5.6) becomes

$$ \xi_b^* = \ell \left( \frac{(-1)^{b+1}}{4} \mathcal{J}:\boldsymbol{\varepsilon} \ell^0 + 0 \ell^1 + \mathcal{O}(\ell^2) \right), \qquad (5.10) \\ \psi_b^* = 0 \ell^0 - \frac{3(-1)^{b+1}}{4} \mathrm{curl}(\mathcal{J}:\boldsymbol{\varepsilon}) \ell^1 + \mathcal{O}(\ell^2). $$
**Remark 5.2.** As mentioned earlier, the homogenized energy $\Phi^*$ approximates only the integral term $\Phi$ that is produced by the continualization of the discrete energy $\Phi_d$ but not its companion boundary term $\Phi_{bt}$, see (3.37). The boundary terms in (5.9) must therefore be added up to $\Phi_{bt}$. This will be covered in future work, along with an analysis of the boundary layers.

**Remark 5.3.** The strain-gradient energy appearing in the second line of (5.9) is not positive, as is often the case in higher-order homogenization. As discussed at the end of Section 4.1 of our previous paper [AL23], this lack of positivity calls for a regularization of the functional, a point which we will address in future work.

## 5.4. Short track: homogenizing the inextensible lattice

In the previous sections (§5.1–5.3), we have homogenized a honeycomb lattice made up of extensible beams first, and taken the inextensible limit $\chi \to 0$ in the homogenized energy next. Here, we show how it is possible to directly homogenize a honeycomb lattice made up of inextensible beams. This yields the same results more elegantly.

The inextensible case is treated in a second Mathematica notebook, also distributed with the library, named homogenize-inextensible-honeycomb nb. The input to the homogenization is identical to that for the extensible lattice (§5.1), except that:

* the beam elements are allocated with the keyword “inextensibleBeam” (as opposed to “beam” in the extensible case) and make use of the properties EI and $\ell$ only (no EA or $\chi$ parameter),

* an option rankDeficiency=1 is passed to the library.

If the rankDeficiency option is overlooked, the solution of the leading-order problem fails as the underlying system of linear equations is deficient (with rank 1) in the inextensible case: the library then stops with an error, identifies the incompressibility condition (5.8) as the root of the deficiency, and suggests that the procedure is run again with rankDeficiency=1.

With the option rankDeficiency=1, the library extends the vector $l$ representing the macroscopic degrees of freedom with a fourth slot $p=l_4$: using block-matrix notation, we have now

$$l = (\tilde{\epsilon} p), \qquad (5.11)$$

where $\tilde{\epsilon}$ are the usual, 3 components of the 2D macroscopic strain tensor $\epsilon$ arranged in a vector, see (3.7). The homogenized energy (4.5) produced by the library must be interpreted in the light of this new definition of $l$. The new variable $p$ is the Lagrange multiplier associated with the incompressibility condition (5.8) and can be interpreted as a pressure.

The output for the inextensible honeycomb lattice, as worked out in the notebook, can be summarized as follows:

* A solvability condition on $l$ is issued by the homogenization engine, which, upon examination, is nothing but the incompressibility condition (5.8). This solvability condition warrants that the right-hand sides of the linear problems involving the rank-deficient matrix are indeed part of its image, so that the inversion of the rank-deficient problem is mathematically possible: see Equation [E.11] in [AL23].

* The elasticity tensors $K$, $A$, $B$ and $k$ and localization tensors $Y$, $Y'$ are calculated. When interpreted with the help of (4.3) and (4.5), they yield the same energy (5.9) and microscopic displacement (5.10) as those from Section 5.3.

In the homogenization library, inextensible lattices are implemented as follows:

* For each element family $\varphi \in \{I, II, III\}$, the inextensibility constraint $\epsilon_\varphi=0$ is handled by adding one row to the constraint matrix $Q$ appearing in (3.32–3.33). Each one of these new rows is filled with 0's, except for a single entry equal to 1 in the column whose index corresponds to the position of the element strain $\epsilon_\varphi$ in $E$, warranting that the generic form of the kinematical constraint $Q \cdot E = 0$ used throughout now includes the inextensibility constraints $\epsilon_\varphi=0$.

* The element stiffness in (2.19) is changed to

$$\mathcal{H}_\varphi(\eta, \mathbf{X}) = \frac{\overline{EI}_\varphi(\eta, \mathbf{X})}{\overline{\ell}_\varphi(\eta, \mathbf{X})} \mathrm{diag}(0, 1, 12) \in \mathbb{T}^{(3,3)} \quad (2D \text{ inextensible beam}), \qquad (5.12)$$

where $\overline{EI}_\varphi(\eta, \mathbf{X})$ is the prescribed map of elastic (bending) modulus—we set $\overline{EI}_\varphi(\eta, \mathbf{X}) = (EI)$ in the present case of a homogeneous lattice. The value (0) appearing in first position of the diagonal matrix in the right-hand side is irrelevant: this is the modulus of the stretching mode which is inhibited by the inextensibility constraint anyway.

The rest of the homogenization procedure is unchanged. Inextensible lattices make use the ability of the homogenization engine to deal with rank-deficient linear problems, as documented in Appendix E of our previous paper [AL23].

# 6. NUMERICAL VERIFICATION

In this section, we set up a numerical procedure to test both the asymptotic validity of our homogenization scheme in the limit of vanishingly small cell size, and its accuracy for small but fixed cell size. We start in this section with honeycomb lattices subjected to various boundary conditions and loading types. More advanced illustrations will be considered in Section 7.
**Remark 5.2.** As mentioned earlier, the homogenized energy $\Phi^*$ approximates only the integral term $\Phi$ that is produced by the continualization of the discrete energy $\Phi_d$ but not its companion boundary term $\Phi_{bt}$, see (3.37). The boundary terms in (5.9) must therefore be added up to $\Phi_{bt}$. This will be covered in future work, along with an analysis of the boundary layers.

**Remark 5.3.** The strain-gradient energy appearing in the second line of (5.9) is not positive, as is often the case in higher-order homogenization. As discussed at the end of Section 4.1 of our previous paper [AL23], this lack of positivity calls for a regularization of the functional, a point which we will address in future work.

## 5.4. Short track: homogenizing the inextensible lattice

In the previous sections (§5.1–5.3), we have homogenized a honeycomb lattice made up of extensible beams first, and taken the inextensible limit $\chi \to 0$ in the homogenized energy next. Here, we show how it is possible to directly homogenize a honeycomb lattice made up of inextensible beams. This yields the same results more elegantly.

The inextensible case is treated in a second Mathematica notebook, also distributed with the library, named homogenize-inextensible-honeycomb nb. The input to the homogenization is identical to that for the extensible lattice (§5.1), except that:

* the beam elements are allocated with the keyword “inextensibleBeam” (as opposed to “beam” in the extensible case) and make use of the properties EI and $\ell$ only (no EA or $\chi$ parameter),

* an option rankDeficiency=1 is passed to the library.

If the rankDeficiency option is overlooked, the solution of the leading-order problem fails as the underlying system of linear equations is deficient (with rank 1) in the inextensible case: the library then stops with an error, identifies the incompressibility condition (5.8) as the root of the deficiency, and suggests that the procedure is run again with rankDeficiency=1.

With the option rankDeficiency=1, the library extends the vector $l$ representing the macroscopic degrees of freedom with a fourth slot $p=l_4$: using block-matrix notation, we have now

$$l = (\tilde{\epsilon} p), \qquad (5.11)$$

where $\tilde{\epsilon}$ are the usual, 3 components of the 2D macroscopic strain tensor $\epsilon$ arranged in a vector, see (3.7). The homogenized energy (4.5) produced by the library must be interpreted in the light of this new definition of $l$. The new variable $p$ is the Lagrange multiplier associated with the incompressibility condition (5.8) and can be interpreted as a pressure.

The output for the inextensible honeycomb lattice, as worked out in the notebook, can be summarized as follows:

* A solvability condition on $l$ is issued by the homogenization engine, which, upon examination, is nothing but the incompressibility condition (5.8). This solvability condition warrants that the right-hand sides of the linear problems involving the rank-deficient matrix are indeed part of its image, so that the inversion of the rank-deficient problem is mathematically possible: see Equation [E.11] in [AL23].

* The elasticity tensors $K$, $A$, $B$ and $k$ and localization tensors $Y$, $Y'$ are calculated. When interpreted with the help of (4.3) and (4.5), they yield the same energy (5.9) and microscopic displacement (5.10) as those from Section 5.3.

In the homogenization library, inextensible lattices are implemented as follows:

* For each element family $\varphi \in \{I, II, III\}$, the inextensibility constraint $\epsilon_\varphi=0$ is handled by adding one row to the constraint matrix $Q$ appearing in (3.32–3.33). Each one of these new rows is filled with 0's, except for a single entry equal to 1 in the column whose index corresponds to the position of the element strain $\epsilon_\varphi$ in $E$, warranting that the generic form of the kinematical constraint $Q \cdot E = 0$ used throughout now includes the inextensibility constraints $\epsilon_\varphi=0$.

* The element stiffness in (2.19) is changed to

$$\mathcal{H}_\varphi(\eta, \mathbf{X}) = \frac{\overline{EI}_\varphi(\eta, \mathbf{X})}{\overline{\ell}_\varphi(\eta, \mathbf{X})} \mathrm{diag}(0, 1, 12) \in \mathbb{T}^{(3,3)} \quad (2D \text{ inextensible beam}), \qquad (5.12)$$

where $\overline{EI}_\varphi(\eta, \mathbf{X})$ is the prescribed map of elastic (bending) modulus—we set $\overline{EI}_\varphi(\eta, \mathbf{X}) = (EI)$ in the present case of a homogeneous lattice. The value (0) appearing in first position of the diagonal matrix in the right-hand side is irrelevant: this is the modulus of the stretching mode which is inhibited by the inextensibility constraint anyway.

The rest of the homogenization procedure is unchanged. Inextensible lattices make use the ability of the homogenization engine to deal with rank-deficient linear problems, as documented in Appendix E of our previous paper [AL23].

# 6. NUMERICAL VERIFICATION

In this section, we set up a numerical procedure to test both the asymptotic validity of our homogenization scheme in the limit of vanishingly small cell size, and its accuracy for small but fixed cell size. We start in this section with honeycomb lattices subjected to various boundary conditions and loading types. More advanced illustrations will be considered in Section 7.
Near the physical edges of the lattice, boundary layers form and the homogenization breaks down [Dum86; LM18]. Until we properly address these layers, they prevent us from making accurate global predictions of the equilibrium solution based on the homogenized theory (4.5)—the analysis of the boundary layers is a significant piece of work and we leave it to a follow-up paper. In the meantime, we give up on making predictions based on the stand-alone homogenized theory (4.5) and keep it tethered to numerical simulations of the full lattice. Our verification method proceeds by extracting the microscopic fluctuations from full numerical simulations of the lattice, and by comparing them to the predictions of the homogenization procedure. This approach is purely local, warranting that the systematic errors in the boundary layers do not interfere with the verification in the bulk. Concretely, our verification method produces color maps of the lattice showing where the homogenization method is accurate and where it is not: it is considered to be valid if the accuracy is good everywhere but in a neighborhood of the boundaries, in the sense that the residual error goes sufficiently rapidly to zero with the cell size (see the convergence plot in Figure 6.4).

The verification procedure builds up on the explicit construction of the continuous fields proposed in Remark 3.2. We start from the nodal displacement and rotation ($\bar{v}_\beta, \theta_\beta) \in \mathbb{R}^{n_a}$ obtained numerically by running discrete lattice simulations based on the energy (2.23). The procedure outlined in Remark 3.2 is implemented through a series of interpolation and averaging steps detailed in Section 6.1 below, delivering interpolations for the macroscopic displacement and the microscopic degrees of freedom $\bar{y}(\mathbf{X}) = (\xi_1(\mathbf{X}) / \eta, \bar{\psi}_1(\mathbf{X}), \dots)$ in the form of continuous functions of $\mathbf{X}$. Here, the overbar marks quantities extracted from the numerical simulation. By differentiating $\bar{\dot{u}}(\mathbf{X})$ symbolically, one can obtain the macroscopic strain $\bar{l}(\mathbf{X})$ and its gradient $\nabla \bar{l}(\mathbf{X})$. The accuracy of the homogenization procedure is assessed by comparing locally the microscopic degrees of freedom $\bar{y}(\mathbf{X})$ that are directly extracted from the simulations, to the prediction $\mathbf{y} = \mathbf{Y} \cdot \bar{l} + \mathbf{Y}' \cdot \nabla \bar{l}$ in (4.3) based on the localization tensors $\mathbf{Y}$ and $\mathbf{Y}'$ obtained via homogenization, combined with the interpolated macroscopic strain $\bar{l}(\mathbf{X})$ and its gradient $\nabla \bar{l}(\mathbf{X})$.

The verification is therefore focussed on the prediction (4.3) of the microscopic displacement $\mathbf{y}$. This is indeed the key prediction of the homogenization method as it slaves the microscopic displacement $\mathbf{y}$ to the macroscopic degrees of freedom. The expression of the homogenized energy (4.5) follows straightforwardly by inserting this expression of $\mathbf{y}$ into the original energy.

## 6.1. Verification procedure

As a first illustration of the verification procedure, we consider a 2D rectangular strip made of a honeycomb lattice subject to compression, as shown in Figure 6.1. The length of the strip is $L=1$ and its width is $\delta = 0.36$. The displacements of the nodes at the left vertical boundary are blocked, $v_{\beta \in lb} = 0$. Those at the right boundary are displaced by an amount $\Delta$, $v_{\beta \in rb} = -\Delta e_1$, where $e_1$ is a unit vector along the long axis of the strip. The rotations are blocked along both short edges, $\theta_{\beta \in lb} = \theta_{\beta \in rb} = 0$. We check convergence for finer and finer lattices, when the beam length $\ell = O(\eta)$ goes to zero ($\ell \to 0$).

The lattice is made up of extensible 2D beams, modeled with the discrete energy (2.19–2.20). All of them are assigned a bending modulus $EI = 1$ and a stretching modulus $EA = 10^{12}$; this warrants that the aspect-ratio parameter $\chi = EI/EA\ell^2$ remains at most of order $10^{-10}$ in all our simulations, so that all beams are effectively in the inextensible regime. The simulations results are compared to the inextensible homogenization theory (§5.3–5.4).

<figure><img src="image_5.png" /><figcaption>**Figure 6.1.** Honeycomb strip in compression: the lattice is loaded by bringing the short edges closer to one another by a distance $\Delta$. Sliding and rotation is blocked on both short edges. We check convergence of the discrete numerical solution when the beam length goes to zero, $\ell \to 0$.</figcaption></figure>

For a given beam length $\ell$, the discrete simulation yields a set of nodal displacements $\bar{v}_\beta$ and rotations $\bar{\theta}_\beta$. For each beam $\alpha$ in the lattice and for every Bravais sub-lattice $b \in \{1, 2\}$, we pick $N_b^\alpha = 6$ nodes from the sub-lattice $b$ in the neighborhood of the beam $\alpha$, and use the corresponding nodal values $\bar{v}_\beta$ and $\bar{\theta}_\beta$ to build second-order interpolations $\bar{v}_b^\alpha(\mathbf{X})$ and $\bar{\theta}_b^\alpha(\mathbf{X})$ (in general, we need $N_b^\alpha = 1 + d + d(d+1)/2 = (d+1)(d+2)/2$ neighbors to build the second-order interpolations in dimension $d$). These interpolations applied to the sub-lattice $b$ are valid in the neighborhood of the element $\alpha$. In Figure 6.2b, we plot the interpolation of the longitudinal displacement $\bar{u}_b^\alpha(X_1) = \bar{v}_b^\alpha(X_1) \cdot e_1$ for the particular case of simple compression (uniaxial stress): the interpolations $\bar{u}_1^\alpha(X_1)$ and