Asymptotic, second-order homogenization of linear elastic beam networks

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 $\boldsymbol{\varepsilon}$ and its gradient $\nabla \boldsymbol{\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.


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].
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 selfcontained simulations so far, as discussed in further detail in the conclusion.

GENERATING THE ELASTIC LATTICE
In this section, we describe the process of generating the elastic lattices.Given a small parameter  ≪ 1, we produce a lattice having cells of size () whose elastic and geometric properties vary over the characteristic length L = (1).This sets the stage for homogenization which captures the effective properties of the lattice in the limit  → 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.

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 .Their position is denoted as  ˜ ∈ ℝ 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 ⩽ b ⩽ 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  belongs is denoted as b().

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

• element families
Elements are grouped into n  families indexed by , having similar elastic and geometric properties.The honeycomb lattice in Figure 2.1a, for instance, has n  = 3 element families,  ∈ {I, II, III}, corresponding to the three possible orientations of the beams.The family to which a particular element  belongs is denoted as   .The density of elements of a particular family  per unit surface (if the space dimension is n = 2) or volume (if n = 3) is denoted as  ˜.
The topological lattice being periodic, the vector  ˜ + −  ˜ − joining the endpoints of a particular element  is the same for all elements  in a particular family  and can be written as where  ˜ maps a family index  to a vector in ℝ d .
By a similar argument, the Bravais sub-lattice b(  ± ) of the endpoints ± of an element  depends on the element family   and not on the element  itself: there exists a function b  ± mapping a family  and an endpoint index ± (head or tail) to a Bravais index b, such that b(  ± ) = b  ± . (2.2) The maps  ˜ and b  ± for the topological honeycomb lattice are given in Table 2.1.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.

Curved reference configuration
The reference configuration of the elastic lattice, shown in Figure 2.1b, is produced by applying a diffeomorphism  : ℝ d → ℝ d combined with a contraction with ratio , to the topological lattice: the position in the reference configuration of the node  is   =  (  ˜).
(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 'complexexponential' diffeomorphism where  r  = cos   1 + sin   2 denotes the unit radial vector and R is the curvature radius of the arc of circle that is the image of the curve X ˜2 = 0.The case of a periodic, non-curved honeycomb lattice corresponds to R → ∞, which yields  = id 2 (identity in ℝ 2 ).Each element  is assigned a conventional center   c in reference configuration.Various definitions of the center with (2.6), we find (2.7) The right-hand side vanishes when  → 0 and the smoothness of  proves the estimate   ± (, ) = () stated in (2.6).By way of illustration, consider elements of type  = I in a curved honeycomb, see Figure 2.1(b).Combining ) and  ˜I = a  1 , we have from (2.7) (2.8) The two other functions  II ± (, ) and  III ± (, ) are obtained in a similar way.

Assumed power-series dependence on η
Expanding Equation (2.7) in powers of , we get where the cancellation at leading order leads to the estimate (2.6).The diffeomorphism  being infinitely smooth by assumption, this expansion could be pushed to any other.More generally, we will assume that any function g(, ) depending on the expansion parameter  can be expanded as a power series in : when we write g(, ) with g = ( p ), we imply that g is of the form g(, ) =  p (g ( p) () +  g ( p+1) () + ⋅ ⋅ ⋅), (2.10) where the g (k) 's are smooth functions of  that do not depend on .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 =   ± 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.

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   are denoted as   ∈ ℝ d , and   ∈ ℝ nr , respectively (in 2D, the infinitesimal rotation   =   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   is common to the adjacent elements at any particular node   .The nodal degrees of freedom are therefore where (2.11) For the 2D elements discussed here, d = 2, n r = 1 and n n = 3.Each element family  is characterized by a set of n E element strain measures, which are collected into an element strain vector   ∈ ℝ nE  .In 2D, beam elements uses n E = 3 strain measures, namely the extensional strain   , the bending strain   and the shear strain   , which are defined by z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z (2.12) where ℓ  =    + −   − is the undeformed length,   = (   + −   −) /ℓ  is the unit tangent,   =  3 ×   is the undeformed unit normal and  3 =  1 ×  2 is the unit normal to the Euclidean plane, see Figure 2.2.The quantities (  ,   ,   ) 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   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.We observe that the undeformed element's length ℓ  , unit tangent   and unit normal   can each be expressed as functions of the element family   , of the scale separation parameter  and of the element center   c .Indeed, with help of (2.6), we have (2.13) The dependence of these function on  applies to lattices having graded geometrical properties (see §7.3).The functions ℓ ˘ ,  ˘ and  ˘ have regular expansions in  of the form (2.10) as can be checked using (2.9).In view of this, the element strain   in (2.12) can be rewritten in block-matrix notation as  ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( (  ( ( + ˘(,  )/ℓ ˘(,  ) 1/2 (2.15) The function   (, ) can be expressed as a regular expansion in  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   applicable to 2D beams, and use the generic form (2.14) of the displacement-strain relation instead.

Strain energy of the lattice
A beam  made up of a Hookean material and having a circular cross-section is characterized by a stretching modulus (E A)  and a bending modulus (E I)  .The linear beam theory yields the strain energy w  of the beam as , which can be rewritten as (no implicit sum over ) in terms of the beam's aspect-ratio parameter (2.17) For the slender beam theory to be applicable, the beam's radius r  must be much smaller than its length, r  ≪ ℓ  , implying that the aspect-ratio parameter   ∼ (r  /ℓ  ) 2 ≪ 1 is small.To handle the case of lattices whose elastic properties vary slowly in space, we assume that the beam's elastic constants (E A)  and   are prescribed as smooth functions (E A )  (, ) and  ˜(, ) of the scale separation parameter  and of the element center   c , with an additional dependence on the discrete family index , i.e., In view of (2.16), we define the element stiffness matrix   ∈  The element stiffness is assumed to obey the scaling assumption 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 Φ = ( 0 ), which is convenient, see Equation (3.42) below.Since ℓ ˘ = () by (2.13), the stretching modulus E A  in (2.19) will be viewed as a quantity of order in stretching-dominated lattices.
Remark 2.2.As noted below Equation (2.17), the function  ˘(, ) is small,  ˘(, ) ≪ 1. Mathematically,  ˘ will be defined as a power of the aspect-ratio .Different power laws are considered in the different examples, such as  ˘ = () for the honeycomb, see (5.7), or  ˘ = ( 2 ) for the pantograph, see (7.28).The only assumption we make at this stage is that  ˘ is a function of  with  ˘ = ( 0 )-this includes  ˘ = ( 1 ) and  ˘ = ( 2 ) as particular cases.In any case, the matrix   (, ) appearing in (2.19) may contain terms of different orders in  and may therefore depend on  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 Ω ⊂ ℝ d of the physical lattice that exclude the edges, see Figure 2.3: the strain energy of the lattice Φ d in the region Ω is defined as the sum of the elastic energy w  over all beams  whose center   c (shown as a red dot in Figure 2.3) lies inside the domain Ω, Points where the loading is not smooth (such as points of application of point-like forces) are also excluded from the domain Ω, along with a finite neighborhood around such points.The subscript 'd' in Φ d stands for 'discrete': a continuous variant Φ of the lattice strain energy, bearing no subscript, will be introduced later.Remark 2.3.In Sections 7.1, we treat the case of an elastic truss whose bending modulus E I() 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   depends explicitly on , see Equation (7.22) for the arch.

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 () collecting the degrees of freedom which are fixed during homogenization ( §3.3) and a vector () collecting those that are relaxed during homogenization ( §3.4), • introduce an assembled strain vector () in terms of (), () and their gradients ( §3.5), • express the strain energy of the lattice in terms of () ( §3.6-3.7).

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 : • a macroscopic displacement () ∈ ℝ d representing the local average of the nodal displacement over the different Bravais sub-lattices, • for each Bravais sub-lattice 1 ⩽ b ⩽ n b , a microscopic displacement  b () ∈ ℝ d and a microscopic rotation  b () ∈ ℝ nr .
Their orders of magnitude are taken as (see Remark 3.1 for a justification), All these functions may depend on the expansion parameter  ≪ 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., where all the 'terms' in the series  (k) (),  b (k) () and ( b ) (k) () are smooth functions of .Throughout this paper and in the interest of legibility, we will use the quantities (),  b () and  b () 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., ∇ k  = ( 0 ), ∇ k  b = ( 1 ) and ∇ k  b = ( 0 ), which we summarize by the formal rule ∇ = ( 0 ). (3.3) As usual in continuum mechanics, the macroscopic rotation () ∈ ℝ nr is a vector extracted from the antisymmetric part of ∇(): in our notation, this extraction is denoted as () = − 1 2 ∇() :  = ( 0 ) where  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 (3.5) By averaging both sides of (3.4) 1 with respect to the Bravais sub-lattices b, one can then confirm that () is the average nodal displacement as announced earlier.
Remark 3.1.The scaling assumptions in (3.1) imply that the element strain   in (2.12) has magnitude Considering the stretching strain, for instance, ) when combined with (3.1) 1,2 and (2.6).A similar reasoning holds for   = ( 0 ) and   = ( 0 ).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.Remark 3.3.Consider the particular case of a rigid-body motion of the lattice made up of a translation  0 and rotation  0 , i.e.,   ,   = ( 0 +  0 ×   ,  0 ) at every node .The macroscopic fields are extracted as outlined in Remark 3.2, which yields the macroscopic displacement as () =  0 +  0 × , a constant macroscopic rotation () =  0 given by (3.8) 2 and zero microscopic displacement  b () =  d and rotation  b () =  nr .By including the term (  ) in the right-hand side of (3.4) 2 , we have made the microscopic displacement insensitive to rigid-body motions of the lattice.

Microscopic degrees of freedom y
Next, we define the microscopic degrees of freedom () as the concatenation into a single vector of those relevant to the different Bravais lattice, where n y is the dimension of , By contrast with , the vector  lists all the degrees of freedom that will be relaxed during homogenization: the discrete energy of the lattice depends on  but the homogenized energy does not.The coefficients 1/  in (3.15) warrant The microscopic displacement and rotation at a particular node can be extracted from  using ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( where  b is defined as (nn,ny) .
We have derived in (3.13) the Taylor expansion of the macroscopic contribution to the nodal degrees of freedom about a generic point  c .The counterpart for the microscopic contributions is obtained by expanding (3.18) as ( ( ( ( ( ( ( ( ( ( ( ( ( ( 𝝃 b where (nn,ny,d,d) , . . . (3.20)
The only point of the macroscopic degrees of freedom  being to compute the element strain by Equation (3.22), it is therefore possible to discard entirely the  and  components in , i.e., to amend (3.10) to (alternate, reduced form of ). (3.24) In this case the length of n l in (3.11) must obviously be changed to n l = n  .The library is flexible: the definition of  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  coincides with the macroscopic strain .We prefer, however, to treat these notions as distinct, as we will occasionally include additional degrees of freedom into , such as the pre-strain magnitude in Equation (7.9).In any case,  collects all the variables on which the homogenized energy depends, so that the energy is always of the form (4.5).

Continualized strain energy
The energy Φ 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  by a 'genuine' integral plus boundary terms-we work with an error of order  3 in the present work, where Φ bt = ∮ ∂Ω . . .da is the boundary term, and Φ is the continualized energy, • the integral term thus obtained corresponds to taking the formal limit  → 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 1 2 (, ) ⋅ (, ) ⋅ (, ) announced in (3.38).We do not need the exact expression of the boundary term Φ bt in the present paper, and therefore leave to future work both the details of the proof and the expression of Φ bt .In this future work, we will combine the boundary term Φ bt with an analysis of the boundary layer to account for the boundary in an effective way.
Combining The above scalings for  and Φ are convenient.They motivate the scaling assumption (2.21).

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  in terms of ,  and their gradients, which we can rewrite as  = (, ; (), ∇(), . . .; (), ∇(), . . .), • the expression (3.32) of the kinematical constraint, • the expression (3.38) of the energy.
We will homogenize the continualized energy Φ, ignoring the boundary terms Φ 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.

SECOND-ORDER HOMOGENIZATION
In this section, we summarize the homogenization procedure described in [AL23] which applies readily to the canonical form of lattice energy Φ = Φ[, , ] derived in Section 3. Here, the square brackets refer to a functional dependence on the arguments  and .

Principle of the homogenization method
The macroscopic state vector  is equal to the macroscopic strain  ˇ, see (3.24) and (3.7).However, to make future extensions easier, we treat  and  ˇas distinct entities, viewing  as the list of degrees of freedom that are fixed during homogenization (which may be more than just the macroscopic strain  ˇ).
The homogenization procedure described in our previous work [AL23] can be summarized as follows.With  a prescribed function of , we revisited homogenization by working at the energy level and showed that the microscopic degrees of freedom () are the solution of a stationary-point problem, This stationary-point problem concerns interior points of the domain Ω only, and proceeds order by order in the expansion parameter , the gradient effect being captured as a higher-order correction in .The canonical form derived in Section 3 provided us with the energy functional Φ[, , ] which we can now plug into the homogenization engine developed in our previous work [Aud23].
The stationary point  ⋆ obtained in (4.1) is a functional of , as indicated by the square brackets.It is also a function of  whose values are denoted as  ⋆ [, ]().Equation (4.1) effectively slaves the microscopic degrees of freedom  to the macroscopic state .
A homogenized energy functional Φ ⋆ is obtained by inserting the stationary point  ⋆ into the original energy, The stationary point problem (4.1) is similar to that governing the equilibrium of the lattice, except that: • only  is allowed to vary,  being kept fixed; • the potential of the external loading Ψ is ignored, only the strain energy Φ being considered; • the energy is written in the domain Ω strictly contained in the physical domain of the lattice, see Figure 2.3, and boundary terms Φ bc are ignored.
Once homogenization is complete, the following additional steps should be addressed: the boundary layers taking place outside the domain Ω must be analyzed, the boundary terms Φ bc must obtained, and the total potential energy-including the potential of the external loading Ψ-must be optimized with respect to .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 Φ ⋆ : 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 Ψ during homogenization is as follows.It follows from standard assumptions on the magnitude of the load that the potential Ψ[, ] is a function of  but not of : since  is fixed during homogenization, it is then justified to ignore the external potential Ψ[, ] during homogenization, and to restore it in the final step, when the total energy Φ ⋆ [, ] + Ψ[, ] is relaxed with respect to -see also the discussion in Section 3.1 in our previous work [AL23].
The homogenization method treats the scale separation parameter  ≪ 1 as an expansion parameter and solves the stationary point problem (4.1) order by order in : the optimal microscopic displacement  ⋆ is obtained in the form of an expansion, where the localization tensors  and ′ are an output in symbolic form by the homogenization procedure.They scale as 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 Φ ⋆ in (4.2) is obtained in the form of an expansion as well, which can be written in compact form as

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  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 × 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  l (, ),  l ′(, ),  l ′′(, ),  y (, ),  y ′(, ) and  y ′′(, ) of the strain expansion (3.34); • the constraint extraction matrix  in (3.33) 1 ; • the global stiffness matrix (, ) in (3.39), in symbolic form; • the list  of variable material parameters declared by the user is forwarded to the homogenization engine.This list  declares all the symbols entering in the expressions of ,  l ,  l ′,  l ′′,  y ,  y ′ or  y ′′, that vary in space.For lattices having homogeneous properties,  = {} is set to be an empty list.For lattices having graded properties, one can either work out the explicit dependence of ,  l , . . . on  as in (3.34) and set  = {X 1 , . . ., 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  was implicit in our previous work [AL23] and is now denoted explicitly by the argument .As long as the orders of magnitude are preserved, nothing changes when the input to the homogenization procedure depends on .
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 ( y (, ) ⋅ ) ⋅ (, ) ⋅ ( y (, ) ⋅ ) > 0 must be strictly positive for any  such that  ⋅  y (, ) ⋅  =  (see equation (11) in our previous work [AL23]).We rely on the homogenization engine to check this condition.Whenever present, offending microscopic vectors  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 .

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.

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  = 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.Since a contraction factor  is applied to the underlying topological lattice (whose beam length is a = (1)), the generated lattice has beam length ℓ =  a ≪ 1, see also Figure 2.1b.In view of the value of  ˜ = 2 3 3 a 2 in Table 2.1 and of Equation (3.40), the element density   is given, for each one of the three element families  ∈ {I, II, III}, by (5.1) With  = id 2 , the element centers   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 homogenizeextensible-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  = 3 and the array b  ± characterizing the underlying topological lattice, see Table 2.1, • the reference configuration is characterized by   from (5.1) and by the relative position   ± = ±   ˜ 2 of the endpoints with respect to the element center (midpoints), where  ˜ is available from Table 2.1, • each one of the n  = 3 element type is assigned the stiffness matrix   in (2.19) applicable to extensible beams, with constant material parameters; indeed, for this particular lattice, the stretching modulus E A, the beam length ℓ =  a and the aspect ratio parameter  are all (i) invariant in space and (ii) identical for all three element families  ∈ {I, II, III}, • the variable material parameters  = {} is set up as an empty list since the lattice is homogeneous.
The parameters E A, ℓ and  are treated as symbolic parameters by the library, and the homogenized energy is delivered in the form of a function of E A, ℓ and .
Remark 5.1.We noted earlier that the aspect-ratio parameter  ≪ 1 needs to be small for the beam theory to be applicable, see the discussion below Equation (2.17).Mathematically,  has to be linked to the small parameter , see Remark 2.2.Remarkably, we can carry out homogenization without having to specify how  scales with : the only assumption we made is  = ( 0 ) and this does not exclude more specific assumptions such as  = ( 1 ) or  = ( 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  → 0 in the homogenized energy ( §5.3).

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 Φ ⋆ [] in the form (4.5) that makes use of the strain vector  =  ˇin Mandel representation, see (3.7) and (3.24).We systematically rewrite it in terms of the standard 2D strain tensor  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, (5.2) and its Lamé parameters are given by This matches the previous results of [VDP14; GLTK19].In the slender limit ( → 0), we have  =  E A / ℓ = (E I /ℓ 3 ) ≪  = (E A/ℓ ), from which we conclude that (i) the lattice is bending-dominated and not stretching-dominated (due to  ∝ E I) and (ii) is incompressible in this limit (due to  ∝ E A, see (5.8) below).
The microscopic displacement underlying this leading-order energy prediction (5.2-5.3) is obtained from the localization tensor  appearing in (4.3) as • First order.When pushed to first order, the homogenization procedure delivers the homogenized tensors  =  and  = .In view of (4.5), and as usual with centrosymmetric lattices [ABB10], this implies that the energy has no contribution of order .As a result, the leading-order approximation (5.2) is exact up to ( 2 ) and not just ().
• Second order.In the Mathematica notebook, both the bulk term t:: and the boundary term t: : ( ⊗ ∇ ⊗ ) entering at order  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  → 0 are nevertheless given in Section 5.3.
The first-order correction to the microscopic displacement  b [1] ⋆ and  b [1] ⋆ =  that underlies this second-order energy correction is encoded in the localization tensor ′ appearing in (4.3): interpreting the result, we get where curl  = ∂a 2 /∂X 1 − ∂a 1 /∂X 2 denotes the (scalar) curl of a vector field () in 2D.

Inextensible limit
We now take the slender limit  → 0 in the homogenization results from Section 5.2.Specifically, we address the distinguished limit in which the slenderness parameter  and the scale separation parameter  both go to zero with (E A)  = (). (5.7) The equivalent compression modulus  ∼ 3) then becomes infinite while the shear modulus  ∼  E A ℓ ∼ (1) remains finite.For better legibility of the results, we eliminate E A in favor of E I everywhere using E A = E I /  ℓ 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  → ∞ in (5.3) 1 , implying that the equivalent Cauchy medium is incompressible, tr () = 0, ∀. (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  → 0. The result is obtained in the companion Mathematica notebook in the form where div  = ∂a 1 /∂X 1 + ∂a 2 /∂X 2 denotes the divergence of a 2D vector field ().The successive lines in the righthand side represent the leading-order contribution-of order E I/ℓ 3 ∼ E A /ℓ = (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),  2 denotes the identity matrix in dimension 2 and the constant tensor  has been defined in (5.5).Due to this  term, the equivalent strain-gradient continuum defined by Φ ⋆ 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 ( ( ( ( ( ( ( ( ( ( (−1) b+1 4  :  ℓ 0 + 0 (5.10) Remark 5.2.As mentioned earlier, the homogenized energy Φ ⋆ approximates only the integral term Φ that is produced by the continualization of the discrete energy Φ d but not its companion boundary term Φ bt , see (3.37).The boundary terms in (5.9) must therefore be added up to Φ 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.

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  → 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 E I and ℓ only (no E A or  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  representing the macroscopic degrees of freedom with a fourth slot p = l 4 : using block-matrix notation, we have now where  ˇare the usual, 3 components of the 2D macroscopic strain tensor  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 .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  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 , ,  and  and localization tensors , ′ 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  ∈ {I, II, III}, the inextensibility constraint   = 0 is handled by adding one row to the constraint matrix  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   in , warranting that the generic form of the kinematical constraint  ⋅  = 0 used throughout now includes the inextensibility constraints   = 0.
• The element stiffness in (2.19) is changed to ℓ ˘(, ) diag(0, 1, 12) ∈  (3,3) (2D inextensible beam), (5.12) where E I  (, ) is the prescribed map of elastic (bending) modulus-we set E I  (, ) = (E I) in the present case of a homogeneous lattice.The value (0) appearing in first position of the diagonal matrix in the righthand 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].

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  ¯,  ¯ ∈ ℝ nn 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  ¯() =  ¯1() / ,  ¯1(), . . . in the form of continuous functions of .Here, the overbar marks quantities extracted from the numerical simulation.By differentiating  ¯() symbolically, one can obtain the macroscopic strain  ¯() and its gradient ∇ ¯().The accuracy of the homogenization procedure is assessed by comparing locally the microscopic degrees of freedom  ¯() that are directly extracted from the simulations, to the prediction  =  ⋅  ¯+ ′ ⋅ ∇ ¯in (4.3) based on the localization tensors  and ′ obtained via homogenization, combined with the interpolated macroscopic strain  ¯() and its gradient ∇ ¯().
The verification is therefore focussed on the prediction (4.3) of the microscopic displacement .This is indeed the key prediction of the homogenization method as it slaves the microscopic displacement  to the macroscopic degrees of freedom.The expression of the homogenized energy (4.5) follows straightforwardly by inserting this expression of  into the original energy.

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  = 0.36.The displacements of the nodes at the left vertical boundary are blocked,  ∈lb = .Those at the right boundary are displaced by an amount Δ,  ∈rb = −Δ  1 , where  1 is a unit vector along the long axis of the strip.The rotations are blocked along both short edges,  ∈lb =  ∈rb = 0. We check convergence for finer and finer lattices, when the beam length ℓ = () goes to zero (ℓ → 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 E I = 1 and a stretching modulus E A = 10 12 : this warrants that the aspect-ratio parameter  = E I E A ℓ 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).build the second-order interpolations in dimension d).These interpolations applied to the sub-lattice b are valid in the neighborhood of the element .In Figure 6.2b, we plot the interpolation of the longitudinal displacement u ¯b  (X 1 ) =  ¯b  (X 1 ) ⋅  1 for the particular case of simple compression (uniaxial stress): the interpolations u ¯1  (X 1 ) and u ¯2  (X 1 ) on the two Bravais sub-lattices are different (dashed lines), which points to the existence of a microscopic displacement.
u ¯(X 1 ) u ¯B′  ( ( ( ( ( ( ( ( ( ( ( ( ( ( ) . (6.3) We compare this to the leading order ( ⋆ ) [ where a = ( 0 ) is the conventional length of the beams in the topological lattice and Δ is the conventional magni-tude of the applied load in the numerical simulations.The denominators in right-hand sides in (6.5) make the errors e [⋅]  's independent of both these conventions.An interpolation present in the right-hand sides is evaluated at the center   c of the element .
In Figure 6.3, these errors are painted on the deformed lattice, for two different beam lengths ℓ , using a logarithmic colormap.We observe a significant and consistent reduction on the magnitude of the error when we go from no homogenization (e [×]  ) to leading-order homogenization (e [0]  ), and then to second-order homogenization (e [1]  ).As anticipated earlier, this improvement takes place everywhere in the lattice, except in a layer that is a few cells thick at the boundaries, where it does not decrease noticeably. on the microscopic displacement predicted by the leading-order (i = 0) or second-order (i = 1) homogenization, compared with its magnitude (i = ×) in the raw numerical solution, for two different lattices (coarser one in the left column, finer one in the right column).The error is painted over the deformed lattice as a base-10-logarithmic colormap.An arbitrary magnification factor is applied on the displacement for plotting the deformed lattice.The numerical simulations were produced with an imposed displacement Δ = 0.1 but the results are rescaled in such a way that this value of Δ doest not affect the colors.The consistent decrease of the error in logarithmic scale with the order of homogenization is a sign that the homogenization is correct.
This first visual indication that the homogenization is correct can be confirmed by a convergence test with respect to the beam length (ℓ → 0).We select a point in the interior of the rectangular domain, shown by the blue cross in the inset in Figure 6.4, and we compute the local errors e for ℓ → 0. The fact that the order of the error increases from 1 to 2 when we switch from leading-order to second-order homogenization is a numerical evidence that the latter is asymptotically correct.A similar behavior as in Figure 6.4 is obtained when a different target point is chosen (data not shown).For points close to the boundary, the scaling regime seen in Figure 6.4 does not start until ℓ becomes significantly smaller than the distance to boundary (data not shown). and e [1]  computed at the nearest element  to a fixed target (blue cross in the inset) as a function of the 'microscopic' element length ℓ .The blue dashed lines are linear fits with slope 1 and 2, respectively.

Sheared strip, with or without a hole
Next, we test the rectangular strip when subjected to a shearing displacement Δ, as shown in Figure 6.5.The displacement is still blocked along the short boundary on the left-hand side ( ∈lb = ) and the nodes along the short boundary on the right-hand side are moved vertically ( ∈rb = +Δ  2 ).In addition, the nodes along both short boundaries can now rotate freely.The same testing procedure is applied.The colormaps of the homogenization errors are shown in logarithmic scale in Figure 6.6 for two different beam lengths (a base-10 logarithm is used throughout in the paper).Here again, the higher-order correction significantly increases the accuracy of the prediction, except in a small layer along the boundaries.In the next test, we punch the lattice by removing a disk-like region with radius r = 0.04 in its center.The boundary conditions are identical to those in Figure 6.5.This relatively small hole causes stress concentration, visible on the raw output of the numerical simulations in the first row in Figure 6.7.The typical flower-like stress concentration pattern survives to leading-order homogenization, see the second row.When higher-order homogenization is used (third row), however, it almost disappears from the error map, indicating that second-order homogenization accurately captures the fast variations of the stress in the vicinity of the hole: the error goes to zero quickly away from the hole.The existence of a 'safe zone' surrounding the hole (extended blue region in this third row, where the relative error is as small as ∼10 −3 ) suggests that the punched strip could be very accurately represented by combining effective boundary conditions along the edge of the hole with a higher-order continuum model for the strip.This is probably not possible with leading-order homogenization (second row) as the region of influence of the hole stretches all the way to the lateral boundaries.Honeycomb strip with a hole, subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure 6.3 and 6.5.Same boundary conditions as in Figure 6.5.

Cracked strip in tension
Before moving on to more complex examples, we consider the cracked version of the rectangular honeycomb strip shown in Figure 6.8, having an aspect-ratio / L = 0.78.The crack is obtained by removing the beams that intersect the segment with length a, making an angle π/6 with the direction  1 , and whose midpoint is at the center of the lattice.Crack propagation is not considered.The same boundary conditions are used as for the shearing test, but a longitudinal displacement  ∈rb = Δ  1 is imposed on the right-hand-side boundary.Colormaps of the homogenization error are shown in Figure 6.9 for ℓ = 0.0125.Similar to what has been obtained for a punched lattice, some stress concentration is observed at both crack tips.When second-order homogenization is used, not only does the homogenization error drop significantly, it is also very small except in the intermediate neighborhood of the crack tips and lips and of the edges of the rectangular region.The fact that the perturbation at the crack tips are disconnected from each other suggests that this geometry can be accurately represented by a higher-order continuum model endowed with effective boundary conditions for the crack tips, crack lips and free boundaries.

Inhomogeneous elastic properties
We consider an extension of the honeycomb compression test done in Section 6.1, now using a lattice whose bending rigidity is graded in space: a beam  with midpoint position   c is assigned a stiffness E I(  c ), where ( ( ( ( ( ( ( ( ( ( 1 + 5 tanh 2 ( ( ( ( ( ( ( ( ( ( ( ( 6 2  ( ( ( ( ( ( ( ( ( ( (− ( ( ( ( ( ( ( ( (−curl(:)+ ( : ) ∧ ∇(E I) where the constant tensor  has been defined in (5.5), the curl operator has been introduced in (5.6), and ∧ denotes the 2D wedge product,  ∧  = a 1 b 2 − a 2 b 2 .The only difference with the homogeneous case (5.10) is the presence of the gradient term ∇(E I) in the microscopic rotation  b .The gradient effect has now two contributions, one associated with the strain gradient ∇ and the other one with the gradient of elastic properties ∇(E I).
Remark 7.1.The statement (7.2) informs the library that the parameter E I entering in the stiffness  varies in space, implying that the gradient must be calculated by the chain rule as ∇ = d d(E I) ⋅ ∇E I.The library does not know about the specific profile E I() in (7.1) and treats the gradient ∇(E I) symbolically: the explicit expression of ∇(E I) must be inserted by the user into the output (7.3) of the library.
The validity of the predictions (7.3) are tested using the same procedure as earlier.The results are shown in Figure 7.2.In the discrete simulations, the deformation is concentrated in the soft part of the lattice; by a Poisson's effect, the lattice forms a neck in the center.The rapid variations in bending stiffness creates strong gradients that make the predictions of the leading-order homogenization relatively inaccurate, see Figure 7.2(b).Taking this gradient effect into account improves the accuracy by approximately one order of magnitude, see Figure 7.2(c).Both terms in (7.3) 2 are equally important: dropping either one deteriorates the accuracy significantly (plots not shown).

Kagome lattice with pre-strain
In this section, we change the lattice topology to Kagome and analyze the effect of a non-uniform pre-strain, as produced typically by thermal expansion in the presence of a localized heat source.The rectangular lattice is shown in Figure 7.3a: it is pinned at two adjacent corners and there are no other loading than the pre-strain.It is made up of beams with length ℓ = 0.006579, aspect-ratio parameter  = 0.002, stretching modulus E A = 1 and bending modulus EI = EA ℓ 2 .
In the discrete simulations, the extensional pre-strain p() is accounted for by changing the strain energy (2.20) of the beams to where  0 (p) = (p, 0, 0).(7.4) The pre-strain vector  0 reflects the ordering convention in (2.12): the pre-strain p being extensional, it offsets the first slot of   representing the stretching strain   .In (7.4), the value of the pre-strain p in any particular beam  is picked based on a pre-strain map p() evaluated at the midpoint   c : in the simulations, we use the map plotted in Figure 7.3a, which has a maximum at the center  =  of the upper boundary, and decreases exponentially away from this point.The pre-strain magnitude is set to Δ = 0.3 in the simulations but this value is a matter of convention as all our results are rescaled by Δ.The results of the discrete simulations are shown in Figure 2.20 (deformed lattice) along with the homogenization results (colors).We proceed to homogenize the lattice.The underlying topological Kagome lattice makes use of the conventions shown in Figure 7.3b.It is specified in the homogenization code based on the properties listed in Table 7.1.
Properties of the topological Kagome lattice shown in Figure 7.3b, used as an input to the homogenization code: space dimension d, number of Bravais sub-lattices n b , number of element families n  , density  ˜ per unit area of element belonging to any particular family , connectivity matrix b  ± and end-to-end vector   =   + −   − for a beam of type  in reference configuration.
In the homogenization library, the pre-strain is taken into account by expanding (7.4) as Contrary to what was assumed earlier in (3.29), the element energy (7.6) is no longer homogeneous of degree 2 with respect to the unknown strain   .To work around this, we incorporate the 1's appearing in (7.6) into the macroscopic strain , amending its definition (3.25) as-mind the trailing 1,  = ( ) ) ) ) ) ) ) ) ) ) ) ) ) ).
(7.7) Equation (7.6) can then be implemented by incorporating the square brackets appearing in (7.6) into the element contribution   g (3.31), which we now define as ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( Now, since  is required to depend linearly on ,  and their gradients by (3.34), the updated definition of  in (7.7) requires that we append a trailing entry 1 in the macroscopic strain  as well, modifying (3.24) into and that we propagate it to : to do so, we extend the dimension of the tensors  l , . . .,  y ′′ appropriately and fill the new entries with 0's, except that in the lower-right corner of  l which is set to 1.In the library, the content of the vectors , , , etc. is not fixed once for all but can instead be chosen in a flexible way, as documented in the Implementation notes that are distributed along with the library: this makes the implementation of the above changes straightforward.
From the user's perspective, the pre-strain is dealt with by including an option includeUnitStrain=True that sets up the trailing 1's in both  and , and by providing an optional argument prestrain={"ϵ":p} to the constructor for beam elements.The pre-strain is represented by a pure symbol p that is treated in a similar way as the variable stiffness in section 7.1: it is declared as a variable material parameter by the assignment m={p}, and the homogenization results are given symbolically in terms of  = (p) and its gradient ∇ = (∇p), as well as on  and its gradients ∇, etc.To interpret these results, the user must take into account the fact that (∇l) 4i = ∂(l 4 )/∂X i = ∂1/∂X i = 0 in view of (7.9).
In the limit E I →0 of perfectly flexible beams, the Kagome lattice possesses a number of zero-energy modes [HF06].One of these zero modes has infinite wavelength and therefore survives homogenization: it involves a rigid rotation of adjacent triangles in opposite directions and translates the Bravais sub-lattice having index b by a vector  b , Indeed, the corresponding microscopic displacement  = ( 1 0 2 0 3 0 ) makes all the extensional degrees of freedom in  l ⋅  vanish.The Kagome lattice is homogenized in the Mathematica notebook homogenize-kagome-prestress.nbdistributed along with the library.The lattice made up of extensible beams, corresponding to the element stiffness matrix in (2.19).The microscopic displacement computed by the library takes the form ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( (  ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( (  ( ( ) where  is the Levi-Civita (purely anti-symmetric) symbol in dimension 2,  =  1 ⊗  2 −  2 ⊗  1 and the coefficients c i are rational functions of  that are bounded for  → 0, (7.12) The compact tensorial expressions in (7.11) has been obtained by applying to the raw output of the homogenization code a procedure similar to that yielding the irreducible form of elasticity tensors in the presence of symmetries [AKO17].In Figure 7.4, the expression (7.11) of the microscopic displacement is verified using the usual procedure.The leading-order homogenization fails to give accurate predictions in regions where the gradient of pre-strain is important and the improvement brought about by second-order homogenization is very significant there.If the correction proportional to ∇p is dropped in the second-order homogenization results (7.11) 2 , the agreement with the numerical solution degrades noticeably (data not shown).Remark 7.2.We have considered a Kagome lattice made up of extensible beams.Indeed, it is neither possible to use springs E I = 0 as the lattice then possess a zero-energy mode-the library aborts, printing out the zero mode (7.10)-, nor inextensible beams ( = 0) as this makes the lattice fully rigid-the library aborts, reporting a compatibility condition in the form () =  for all .
Remark 7.3.The  −1 coefficient in factor of c 1 in the right-hand side of (7.11) 1 implies that the amplitude of the soft mode  b becomes larger and larger in the limit of inextensible beams,  → 0. This is not surprising as, by definition, this mode has a vanishing elastic energy in the limit.A detailed analysis of the inextensible limit of the Kagome lattice is beyond the scope of this paper-in principle, a limit energy can be obtained by analyzing the energy of the discrete lattice computed by the library in the limit  → 0, similar to what we did earlier in Section 5.3.

Circular arch
Our next example is an half-annular lattice made up of elastic springs (E I = 0), sketched in Figure 7.5.This lattice features spatial inhomogeneity: both the orientation and the natural length of the springs vary across the annular domain.
The discrete arch is generated as follows.We pick two integers n 1 and n 2 and define In the simulations shown in Figure 7.7, we used (n 1 , n 2 ) = (13, 90).The nodes are indexed by a pair of integers (i, j) with −n 1 ⩽ i ⩽ +n 1 and 0 ⩽ j ⩽ n 2 and their position  (i, j) in reference configuration is set to where  r  = cos   1 + sin   2 . (7.14) The inner and outer radii of the arch are therefore given as exp(± n 1 ) = exp ± .Each node is assigned a Bravais index b with b = 1 if i + j is even (black nodes in the figure) and b = 2 if it is odd (white nodes).As shown in the figure, we set up elastic springs radially between all pairs of adjacent nodes, azimuthally between all pairs of adjacent nodes, and diagonally between pairs of adjacent black nodes (b = 1).
) and (i  + , j  + ), we define the undeformed end-to-end vector as   =  (i  + , j  + ) −  (i − , j − ) , the natural spring length as ℓ  =   , the undeformed unit tangent   =   /ℓ  and the linearized spring elongation z  as z  =   ⋅ ( (i  + , j  + ) −  (i − , j − ) ). (7.15) The spring is assigned an elastic energy where k  is the spring constant which we set up as to represent elastic bars having all the same traction modulus E A. The short edge j = n 2 is blocked and a force is applied on the node (i, j), representing a weight-like force having a constant density (−Δ  2 ) per unit area, where r i = exp ( i) is the radial coordinate of the node.The nodal displacements   at equilibrium are proportional to Δ and independent of the global elasticity constant E A. The invariable forces  (i, j) are taken into account through their potential energy −∑ i=−n1 +n1 ∑ j=0 n2  (i, j) ⋅  (i, j) .
Typical results of the discrete simulations are shown in Figure 7.7 (deformed lattice), along with a comparison to the predictions of homogenization (colors) which we proceed to derive now.
We use homogenization to characterize the limit where both discretization parameters n 1 and n 2 are large: the scale separation parameter  = (1/n 2 ) is then small,  ≪ 1.We limit attention to the case of a stubby arch, and not a slender one, by assuming that the aspect-ratio parameter  = (n 1 /n 2 ) remains finite.
The periodic topological lattice underlying the arch is shown in Figure 7.6a.Its nodes are positioned on a square lattice spanning a rectangular domain,  ˜(i, j) = i  1 + j  2 , with −n 1 ⩽ i ⩽ +n 1 and 0 ⩽ j ⩽ n 2 .There are n  = 6 element families (colors in Figure 7.6a) whose properties are listed in Table 7.2.In this imaginary configuration, the area density of the elements belonging to any particular family  is  ˜ = 1/2.The node positions  (i, j) in (7.14) are produced by applying the map  (i, j) =  (  ˜(i, j) ) in (2.3) to the topological lattice  ˜(i, j) using an exponential-map diffeomorphism  similar to (2.4),  ( ˜) = exp (x ˜1)  r (x ˜2). (7.18) The arch therefore complies with the lattice generation framework presented in Section 2. The Jacobian can be calculated as J  ( ˜) = exp (2 x ˜1), and Equation (3.40) yields the density of springs in a given family  as where the edge vectors  ˜ of the topological lattice are listed in Table 7.2.In the homogenization code, the lattice is initialized with the information listed in Table 7.2, along with the expressions of   in (7.19) and   ± in (7.20).
The following changes are made to the code to handle spring elements, as opposed to the beam elements that have been used so far.The rotational degrees of freedom  b are dropped.There is a single strain (n E = 1) attached to a spring , namely the elongation z  : the element strain in (2.12) becomes   = ( z ).In view of the expression (7.15) of the elongation, we change the displacement-strain matrix in (2.14-2.15) to   , r,  = ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( (   (  ( − where  ˘ is still defined by (2.13).The choice of the spring constants in (7.16) is captured by setting up the element stiffness matrix as   , r,  = ( ( ( ( ( ( ( ( ( E A ℓ ˘ , r,  ) ) ) ) ) ) ) ) )∈ (1,1) (springs), (7.22) so that (2.20) yields the elastic energy of the springs.As anticipated in Section 2, the properties (such as   ,   , etc.) of this curved lattice depend on position  through the polar coordinates r,  .To inform the library that the symbols r and  vary in space, we set up the list of variable material parameters as  = r,  . (7.23) The library furnishes the homogenized tensors in terms  and ∇, treated as symbolic quantities: in the companion Mathematica notebook homogenize-circulararch.nb, we insert the explicit expression of the gradient ∇ =  r  ,    /r in this symbolic output, where    = −sin   1 + cos   2 is the azimuthal unit vector.
The homogenization proceeds automatically as usual, except that the symbolic expressions of  l ; r,  , . . .,  y ′′ ; r,  obtained by (3.23) are so long that they make the procedure hang: to work around this, we replace  l ; r,  , . . .,  y ′′ ; r,  by their Taylor expansions with respect to  to order  3 , which are considerably simpler.These quantities being all of order , expanding to order  3 warrants a relative accuracy of order  2 is preserved, which is consistent with our second-order homogenization scheme.
The second-order homogenization procedure delivers the microscopic displacement in the form  b ⋆ = (−1) b+1 r  ( ( ( ( ( ( ( ( ( ( ( ( ( ( 0  0 +  1 8 2 [ [ [ [ [ [ [ [ [ ( ( ( ( ( ( ( (7.24) where tr  =  rr +   , and  rr ,   and  r are the components of the macroscopic strain in the polar basis ( r ,   ).This prediction is verified by comparing to discrete simulations following the usual procedure in Figure 7.7.The vanishing contribution 0  0 in the parenthesis in (7.24) points to the fact that the leading-order microscopic displacement is zero for this particular truss: it does not help improving the error, e  p +  r ⊗ ( r −  2 ) :  + r ∇tr  +   ⋅ ∇(  : where  p is the polar version of the tensor  introduced in (5.5),  p =  r ⊗  r +   ⊗   with  r = − r ⊗  r +   ⊗   and   =  r ⊗   +   ⊗  r .In the presence of variable properties, the corrective displacement not only includes a term proportional to ∇ but also one proportional to  (see the first term in the square bracket above).The latter is produced by inserting ∇ =  r  ,    /r into the corrective ∇ ⊗  term returned by the library in symbolic form.In Figure 7.9, a shear test is conducted on a block comprising 20 × 10 × 10 unit cells: the displacement of the nodes belonging to the face shown in the left-side is blocked, whereas that of those belonging to the opposite face is set to Δ  3 , with Δ = 0.1 ℓ .In Figure 7.10, a twisting test is conducted on a block comprising 30 × 15 × 11 unit cells: a rotation by an angle Δ is imposed on the second face, with Δ = 0.1 (only the components of the displacement in the perpendicular plane ( 2 ,  3 ) are prescribed, the longitudinal component along  1 remaining unconstrained).lattice is set up using the parameters listed in Table 7.3 (the connectivity matrix b  ± and the edge vectors   are provided to the library but are not included in the table).To represent the springs, the 3D version of the displacementto-strain matrix   in (7.21) is used, along with the element stiffness   in (7.22).For this particular truss, the leading-order microscopic displacement vanishes,  0 = , and the first-order one is returned by the library in the form

Shearing and twisting of a 3D elastic truss
This prediction is verified against numerical simulations in Figures 7.9

A pantograph
Our final example is the 1D elastic truss shown in In all the previous illustrations given in this paper, the gradient effect was a correction to the leading order effect-our homogenization scheme has indeed been derived with this situation in mind.The pantograph is different as the strain-gradient term d 11 /dX 1 ends up contributing to its energy at the same order as the strain  11 , see Equation (7.29) below.Our homogenization scheme can be revisited to deal with this situation, as we show now.A pantograph is a 1D truss made up of springs.All springs have the same constant k, except for the vertical dashed springs whose spring constant is  k.We consider the limit of large elastic contrast ≪ 1.There are n b = 6 Bravais sublattices (symbols and integers) and n  = 13 different types of springs: a unit cell is shown using colored segments.In the limit  = 0 where the dashed vertical springs are absent, the truss possesses a zero-energy mechanism.This is the same lattice as studied in [AS18b], except that beams have been replaced by springs and the weak (dashed vertical) springs have been added.
The pantograph is homogenized in the companion Mathematica notebook homogenize-pantograph-springshighcontrast.nb.In is drawn in the Euclidean plane (d = 2), see Figure 7.11, but the limit model is a one-dimensional bar.We could write specialized code handling 2D lattices that are 1D in the limit but a straightforward work-around is to work in the plane and imagine that we are homogenizing an infinite number of replicas of the pantograph obtained by copying the original one, translating it along  2 and pasting it, repeatedly.Concretely, we set up the dimension as d = 2 in the library and identify the density per area   of each family of elements with its density per unit length   = 1/ℓ , where ℓ = O() is the size of the unit cell.We also check that the homogenized energy computed by the library does not depend on  12 and  22 since the replicas are disconnected.
The homogenized energy returned by the library takes the form where c 5  = 1 + ⋆ is available from the companion notebook.We have not yet specified how  scales with the cell size , and the above result is correct as long as  remains bounded,  = ( 0 ), see Remark 2.2.The limit of large elastic contrast is now addressed by taking the formal limit  → 0 in (7.27).The distinguished limit where both  → 0 and  → 0 with  = (ℓ 2 ) = ( 2 ) (7.28) is particularly interesting, as it makes the two terms balanced in the integral: at leading order, we then have where ℓ / 1/2 = (1) is the macroscopic persistence length of the gradient effect.This energy can be interpreted based on the existence of a soft mode having zero energy in the limit  → 0 (vertical dashed springs are removed).This soft mode involves combined rigid rotations of the triangles, as shown by the light gray circular arrows in Figure 7.11, with an amplitude proportional to the macroscopic strain  11 .The first term in the integral in (7.29) can be rewritten as 1 2 ∫ k  (ℓ  11 ) 2 dX 1 /ℓ , which describes the small energy penalty brought about by the weak springs when the soft mode is activated.Classical homogenization captures only this term, yielding a non-definite energy functional 1 2 ∫ 0 × (ℓ  11 ) 2 dX 1 /ℓ in the limit of large elastic contrast  = 0.The gradient effect solves this issue.The pantograph example shows that the special behavior of high-contrast microstructures, which has been documented in the literature, can be readily recovered by taking the formal limit  → 0 in the output of the homogenization library.The symbolic and second-order capabilities of the library are instrumental here.
Remark 7.5.For consistency with the scaling assumption that the energy is Φ ⋆ = ( 0 ), one should view the spring constant k as quantity of order 1/ℓ 3 = ( −3 ).However, the magnitude of this overall multiplicative constant in the energy is a matter of convention; nothing changes if different scaling assumptions are used.

DISCUSSION AND CONCLUSION
We have proposed a symbolic, second-order homogenization method for elastic lattices.It is distributed as a library, named shoal, implemented in the symbolic calculation language Wolfram Mathematica: given a description of the lattice, the library runs automatically and returns both the localization tensors in (4.3) and the homogenized energy (4.5) in symbolic form.
In the design of the library, emphasis has been placed on versatility: it works in arbitrary dimension d (1 ⩽ d ⩽ 3), can handle lattices featuring graded geometric and/or elastic properties, pre-stress or a large elastic contrast.Elastic lattices made up of springs or beams have been demonstrated and new types of elements can be created easily.We are not aware of previous work addressing both graded properties and pre-stress in the context of higher-order, discrete homogenization.
For a variety of lattices, we have derived in closed analytical form both the nodal displacements and rotations-see Equations (5.10), (7.3), (7.11), (7.25), (7.26)-and the homogenized energy-see Equations (5.9), (7.29)-up to second order.Expressions of this kind are novel to the best of our knowledge.Being symbolic, they fully convey the influence of the various lattice parameters.We hope that this kind of results will be used in the future to gain analytical insights into the effective properties of various lattices.
We have proposed a procedure for verifying the correctness (in an asymptotic sense) of the microscopic displacement predicted by the library, based on comparisons with discrete simulations, see Figure 6.4.We have presented a number of validation examples for which we have plotted physical maps of the homogenization error.The improvement in accuracy brought about by second-order homogenization is significant everywhere except in layers forming at free boundaries, where the assumption of scale separation breaks down.
Until these boundary layers have been analyzed and accounted for in the homogenized energy Φ ⋆ , it is not possible to solve global equilibrium problems for the lattice by making the functional Φ ⋆ stationary.This is one main limitation of the present work, which is just a step towards this goal-a rather ambitious one in the context of higher-order homogenization.The derivation of the other boundary terms in (4.5)-produced by homogenization and not by boundary layers-is original to the best of our knowledge and will need to be included in the homogenized energy functional predicting equilibrium.
There is a second reason that prevents the homogenized energy Φ ⋆ from being used in self-contained simulations, at least in its current form.As noted in Remark 5.3 and in our previous work [AL23], it contains negative straingradient moduli, such as the coefficient −7 appearing in the second line of (5.9), and is therefore non-positive.This is a common finding in higher-order homogenization [LM18;DLSS22].This strong limitation explains the limited popularity of higher-order homogenization so far.In future work, we hope to introduce a regularization strategy that preserves the accuracy of the solution while restoring the desired positivity.
The library is symbolic and does second-order homogenization.When combined together, these two original features are quite powerful.Consider for instance the challenging class of lattice homogenization problems involving a small dimensionless parameter  ≪ 1 (such as the slenderness ratio E I /[E A ℓ 2 ] in beam lattices or the elastic contrast, see §7.5) in addition to the scale separation parameter  ≪ 1. Non-conventional limit behaviors can be obtained in such lattices, which have been approached using specialized methods in the literature.These limit behaviors are often associated with the presence of a large-amplitude, slowly decaying soft mode (Remark 7.3 in §7.2, §7.5).Our approach provides both a unifying perspective and a straightforward approach to these problems: by homogenizing ( → 0) with fixed  first, and taking a distinguished limit in  and  in the homogenized energy next, one can identify the limit energy based on simple dimensional analysis ( §7.5).This approach involves exchanging the limits on  and , which is not mathematically rigorous but yields correct results in all the cases we have tried.
Along the same lines, albeit at a more basic (and less original) level, the homogenization procedure handles bending-dominated and stretching-dominated lattices in a unified way: the class of a particular lattice can be deducted from the symbolic expression of the leading-order homogenized energy, see the discussion immediately after Equation (5.3).
The proposed method is subject to the following practical limitations.(i) The library solves linear algebra problems in symbolic form, which can require a large amount of computing time especially if there are many microscopic degrees of freedom.When this happens, it is possible to revert to numeric homogenization; this is done easily, by assigning numerical values to the lattice parameters.The symbolic computations required to generate the results included in this paper all took from a few seconds to less than a minute on a personal computer.(ii) The homogenized energy and microscopic displacement computed by the library have been cast by hand into compact tensor form, shown for instance in Equation (5.9-5.10): this process could in principle be automated by exploiting the symmetries of the lattice.(iii) We have homogenized discrete lattices made up of 1D beams of springs, without questioning the validity of these 1D models.To address this question, one needs to revert to the equations of bulk elasticity in a periodic domain comprising slender junctions [AB21;DLSS22].
In this work, we have limited attention to linear elastic lattices and thus have not considered material or geometrical nonlinearity.When treated incrementally, homogenization problems for lattices having a more general constitutive behavior (such as non-linearly elastic, visco-elastic or elasto-plastic) yield linearized problems featuring both pre-stress and spatially graded properties (such as the tangent elastic stiffness).This is precisely the kind of problems that our approach can solve.In future work, we plan to apply our homogenization scheme incrementally.By doing so, we hope to gain insights into the homogenizing of lattices featuring a wide range of constitutive behaviors.
where  ⊗2 =  ⊗ ,  ⊗3 =  ⊗  ⊗ .In Section 3, Equation (A.8) is used to write the strain energy of a particular element, which is expressed in terms of the degrees of freedom at the adjacent nodes, into a Taylor expansion at the center of the element.

B.1. Mandel representation ε ˇof the macroscopic strain
We introduce the following constant tensor  ∈  (d,d,n)

B.2. Properties of the tensor 𝓝
The tensor  introduced in (A.6) is similar to  as it delivers an orthogonal (but not orthonormal) basis of antisymmetric tensors in  (d,d) indexed by k in the form ( ⋅  k nr ) where • in dimension d = 2, the only possible rotation index k is limited to k ∈ {3} and we set  3 1 = (1), • in dimension d = 3, 1 ⩽ k ⩽ 3 and the usual definition of  k nr in (B.2) applies.
In view of this and of Equation (A.5) 1 ,  is nothing but the components of  ˆin this basis, up to a sign.The reciprocal formula (A.5) 2 follows from the identity 1 2  T231 :  =  nr , (B.6) which expresses itself the orthogonal character of the underlying basis.

B.3. Elimination of the rotation gradient
In this section, we provide a detailed justification of the identities (A.7), as well as expressions of the geometric tensors ′, ′′, ′, ′′, ′′′ appearing in the right-hand side.These tensors depend on the dimension d of the Euclidean space only.
The rotation in (3.8) 2 can be written in components as  a = −

Figure 2 . 1 .
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  and a diffeomorphism  , see (2.3).(b') Close-up view of a particular element (beam) , belonging to family   = I; note that the conventional element center   c is distinct from its midpoint.The figure was generated for a honeycomb lattice with the diffeomorphism in (2.4) and with parameters a = 1,  = 0.4 and R = 5.

Figure 2 .
Figure 2.2.A 2D beam element, having index  and conventional center   c : (a) geometry in reference configuration and (b) the n E = 3 deformation modes : extension, bending and shearing.

Figure 2 . 3 .
Figure 2.3.To avoid boundary layers, the homogenization is carried out in a sub-region Ω of the lattice (shaded blue region) that excludes the physical boundaries of the lattice.
= (  ) +  b() (  ),   = (  ) +  b() (  ).(3.4) Here, both the nodal displacement   and rotation   are expressed as the sum of a macroscopic quantity (displacement () and rotation ()) and a microscopic one (displacement  b () and rotation  b ()).The macroscopic quantities () and () are functions on  only, and therefore vary slowly, whereas the additional dependence of  b and  b on the Bravais index b = b() 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   in (3.4) 1 is not unique: shifting () by  ˜(), and all the  b ()'s by − ˜() leaves   unchanged.To ensure unicity, we impose the constraint b=1 nb  b () = , ∀.
First, the nodal displacements   in each particular Bravais sub-lattice b ∋  are interpolated by smooth functions  b s (  ) of their undeformed position   : a set of n b such interpolating functions  b s with b = 1, . . ., n b are obtained.Next, () is defined as the average () = 1 nb ∑ b=1 nb  b s ().Finally, the 'continualized' microscopic fields  b () and  b () are obtained by interpolating   − (  ) and   − (  ), respectively, over each Bravais sublattice separately.The constraint (3.5) is then automatically satisfied.

Figure 5 . 1 .
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.

Figure 6 . 1 .
Figure 6.1.Honeycomb strip in compression: the lattice is loaded by bringing the short edges closer to one another by a distance Δ. Sliding and rotation is blocked on both short edges.We check convergence of the discrete numerical solution when the beam length goes to zero, ℓ → 0.

For a given beam
length ℓ , the discrete simulation yields a set of nodal displacements  ¯ and rotations  ¯.For each beam  in the lattice and for every Bravais sub-lattice b ∈ {1, 2}, we pick  b  = 6 nodes from the sub-lattice b in the neighborhood of the beam , and use the corresponding nodal values   and   to build second-order interpolations  ¯b  () and  ¯b () (in general, we need  b  = 1 + d + d (d + 1)/2 = (d + 1) (d + 2) /2 neighbors to

Figure 6 . 2 .
Figure 6.2.Local interpolations u ¯b  () =  ¯b  () ⋅  1 of the longitudinal displacement, illustrated in the particular case of uniaxial compression of a honeycomb lattice.The state of uniaxial compression is obtained by modifying the boundary conditions on the short edges of the lattice shown in Figure 6.1 to let the nodes to slide freely in the transverse direction.(a) Schematic view of the deformed lattice.The link  about which the interpolations are constructed is highlighted in blue.Only the  b  = 6 neighboring nodes from each Bravais sub-lattice b∈ {1,2} that enter into the interpolation are shown.(b) Longitudinal displacement from discrete simulation (u ¯, open and filled disks), interpolated in each Bravais sub-lattice (u ¯b  (X 1 ), dashed lines) and macroscopic displacement (u ¯(X 1 ), solid line).In this simple illustration, the X 2 direction is an invariance direction and all displacements are functions of X 1 only.Following Remark 3.2, we obtain a local interpolation  ¯() of the macroscopic displacement by averaging over the Bravais sub-lattices,  ¯() = 1 n b b=1 nb 0]  and first correction ( ⋆ )[1]    predicted by the homogenized results (5.10) for an inextensible lattice, ( :  ¯)  2 −curl ( :  ¯) .(6.4)In the right-hand sides, we use the estimate  ¯ of the strain, and ∇ ¯ of the strain gradient at the element center  c ,as obtained from the discrete numerical simulation: these values are extracted from the interpolations  ¯(  c ) and ∇ ¯(  c ), based on the definition (3.10) or (3.24) of  in terms of the strain  ˇand the Mandel convention (3.7) for representing the strain  as a flat vector  ˇ.Specifically, we compute at every element  the error e [0]  on the microscopic displacement based on leading-order homogenization, the error e [1]  based on second-order homogenization, as well as its magnitude e [×]  in the raw output of numerical simulations as e − (( ⋆ ) [0]  + ( ⋆ ) [1]  )

Figure 6 . 3 .
Figure 6.3.Honeycomb strip in compression: error e [i] on the microscopic displacement predicted by the leading-order (i = 0) or second-order (i = 1) homogenization, compared with its magnitude (i = ×) in the raw numerical solution, for two different lattices (coarser one in the left column, finer one in the right column).The error is painted over the deformed lattice as a base-10-logarithmic colormap.An arbitrary magnification factor is applied on the displacement for plotting the deformed lattice.The numerical simulations were produced with an imposed displacement Δ = 0.1 but the results are rescaled in such a way that this value of Δ doest not affect the colors.The consistent decrease of the error in logarithmic scale with the order of homogenization is a sign that the homogenization is correct.

Figure 6 . 4 .
Figure 6.4.Honeycomb strip in compression: log-log plot of the homogenization errors e [0] and e [1]  computed at the nearest element  to a fixed target (blue cross in the inset) as a function of the 'microscopic' element length ℓ .The blue dashed lines are linear fits with slope 1 and 2, respectively.

Figure 6 . 5 .
Figure 6.5.Honeycomb strip subjected shear: loading geometry.The nodes on both short edges can now rotate freely.

Figure 6 . 6 .
Figure 6.6.Honeycomb strip subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure6.3.The boundary conditions are represented in a simplified way, the reader is referred to Figure6.5 for details.
Figure6.7.Honeycomb strip with a hole, subjected to shear: homogenization error for two different lattices (columns), using the same conventions as in Figure6.3 and 6.5.Same boundary conditions as in Figure6.5.

Figure 6 . 9 .
Figure6.9.Cracked honeycomb strip in tension: homogenization errors for ℓ = 0.0125, using the same conventions as in Figure6.3.The representation of the boundary conditions is simplified, the reader is referred to Figure6.8 for details.

Figure 7 . 1 .
Figure 7.1.Honeycomb lattice with a graded bending modulus E I() given by (7.1).The bending modulus E I is constant in every beam and is found by evaluating the function E I() at the beam midpoint   c .The boundary conditions are depicted using the same conventions as earlier.The homogenization is based on the inextensible beam model, as earlier.Compared to the homogeneous case ( §5.4), the only change is that the symbolic parameter E I must now be included in the list of variable material parameters,  = {E I}. (7.2)

Figure 7 . 3 .
Figure 7.3.Kagome lattice with pre-strain.(a) Imposed pre-strain map p() and boundary conditions.(b) Underlying topological lattice, made up of n  = 6 beam families and n b = 3 Bravais sub-lattices.In each family, a particular beam is shown using a thick line, with the arrow denoting the conventional orientation.All beams have equal length ℓ = 0.006579, stretching modulus E A = 1 and aspect-ratio parameter  = 0.002.

Figure 7 . 4 .
Figure 7.4.Kagome lattice with spatially graded pre-strain: verification of the microscopic displacement (7.11) based on the discrete numerical solution.Same lattice as in Figure 7.3.

Figure 7 . 5 .
Figure 7.5.Circular arch.This annular lattice is made up of spring elements connecting nodes  (i, j) indexed by two integers, a radial coordinate −n 1 ⩽ i ⩽ +n 1 and an azimuthal coordinate 0 ⩽ j ⩽ n 2 .The short edge j = n 2 is blocked and the other edges are free.A weight-like vertical force  (i, j) is applied on the nodes, with uniform density Δ per unit area.The elasticity of the spring (elements) is set up as follows.Considering a spring  having adjacent nodes (i  − , j  −

Figure 7 . 6 .
Figure 7.6.Circular arch: homogenization conventions.(a) Underlying periodic topological lattice, with the different element families represented by colors.(b) Close-up view of particular element of type  = I in reference configuration, drawn with a large value of the cell size ∼  for the sake of legibility.Note that the element center  I c (target-like symbol, defined as the image of the element midpoint  ˜I c by the diffeomorphism  ) is distinct from the spring midpoint.

Figure 7 . 8 .
Figure 7.8.Unit cell of a 3D lattice made up of springs, with box dimensions ℓ × ℓ × 3 ℓ .The lattice has n b = 2 Bravais sublattices (solid vs. open disks) and n  = 11 element families (colors).The springs are treated as elastic bars having identical stretching modulus E A.

Figure 7 . 9 .
Figure 7.9.3D truss, shearing test.The springs are colored according to the magnitude of the microscopic displacement (no homogenization) or of the error the displacement predicted by first-order or second-order homogenization.To aid visualization, a single layer of springs is shown, the full lattice domain being shown by the black box.

Figure 7 . 10 .
Figure 7.10.3D truss, twisting test.Same conventions as in Figure 7.9.On the face shown in foreground, the rotation with angle Δ is imposed on the  2 and  3 (normal) components of the displacement, the  1 (longitudinal) component remains unconstrained.
Figure 7.11, known as a pantograph.It is made up of springs having identical elastic constants k, except for the vertical springs whose elastic constant is  k ≪ k; we address the case of a small elastic contrast,  ≪ 1.The homogenization of pantograph-like lattices has been studied in a number of work, see for instance [SAd11; AS18b].More broadly, microstructures featuring large elastic contrast have been extensively studied, see for instance the work of [PS97; Bel17; JS20].By recovering known results about the pantograph, we demonstrate how our method can handle the special case of a large elastic contrast.
Figure7.11.A pantograph is a 1D truss made up of springs.All springs have the same constant k, except for the vertical dashed springs whose spring constant is  k.We consider the limit of large elastic contrast ≪ 1.There are n b = 6 Bravais sublattices (symbols and integers) and n  = 13 different types of springs: a unit cell is shown using colored segments.In the limit  = 0 where the dashed vertical springs are absent, the truss possesses a zero-energy mechanism.This is the same lattice as studied in[AS18b], except that beams have been replaced by springs and the weak (dashed vertical) springs have been added.
∞ -smooth functions of  at  = 0, with c 5 (0) = c 6 (0) = 1.The underlying microscopic displacement  b depending only on the dimension d of the Euclidean space, b ∈ ℝ b is the discrete Dirac vector in dimension b, whose components are all zero except for the a-th one which is equal to 1,( a b ) j = { { { { { { { { { { { { { { { { { { 0 if j ≠ a 1 if j = a , with 1 ⩽ j ⩽ b. (B.2)The rationale behind the definition of  is that the quantities in square brackets form an orthonormal basis of the space of symmetric d × d matrices, and that the indices i (first term) and i+ ( j − i) d − j − i − 1 2(second term) provide a sequential numbering of the tensors forming this orthonormal basis, from 1 to n  .The orthonormal basis of symmetric d × d matrices can be therefore be extracted as ( ⋅  k n ) 1⩽k⩽n .The following two identities are a consequence of the orthonormal character of the basis appearing in square brackets in (B.1), T231 :  =  n  ⋅  T231 = (( d ⊗  d ) T1324 ) S12∘S34 , (B.3)where we recall that ∘ stands for the composition of symmetrization operations, see (5.1-2.3).The Mandel representation  ˇ∈ ℝ n of the strain tensor  introduced in (3.7) is nothing but the decomposition of  in this orthonormal basis,  ˇ=  : .(B.4)Conversely, the symmetric strain tensor  can be reconstructed from the vector  ˇby  =  ⋅  ˇ, (B.5) as can be shown using (B.3).

Table 2 .
1. Properties of the topological honeycomb lattice in Figure 2.1a: space dimension d, number n b of Bravais sublattices, number n  of element families, density  ˜ of elements in a given family per unit area, connectivity b  ± , edge vectors  ˜.
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   ∈  (nE  ,nE  ) is always symmetric but not necessarily diagonal.

Table 4 .
1. Lattice specification used as an input to the homogenization procedure.

Table 7 .
2. Properties of the topological square lattice underlying the arch shown Figure 7.6a: space dimension d, number n b of Bravais sub-lattices, number n  of element (spring) families, element density  ˜, connectivity b  ± , and end-to-end vectors  ˜.
are the polar coordinates in reference configuration,  = r  r  .As shown by the target-like symbols in Figures 2.1b' and 7.6b, the link centers   c , conventionally defined as the image by the diffeomorphism of the edge midpoints, see (2.5), do not coincide with the midpoints.The vectors   ± yielding the endpoints ± relative to the center   c (thin black arrows in Figure 7.6b) are given by Equation (2.7) in terms of the family index , the scale separation parameter  and of the polar coordinates r,  of the center   c as

Table 7 .
3. Parameters used to set up the 3D elastic truss in the homogenization code, see Figures 7.8-7.10.The connectivity matrix b  ± and the edge vectors   are omitted as they can be read off easily from Figure 7.8.The homogenization of the 3D truss is carried out in the companion notebook homogenize-3Dtruss.nb.The YANG YE, B. AUDOLY, C. LESTRINGANT