An energy approach to asymptotic, higher-order, linear homogenization

A higher-order homogenization method for linear elastic structures is proposed. While most existing approaches to homogenization start from the equations of equilibrium, this one works at the energy level. We start from an energy functional depending on microscopic degrees of freedom on the one hand and on macroscopic variables on the other hand; the homogenized energy functional is derived by relaxing the microscopic degrees of freedom and applying a formal two-scale expansion. This method delivers the energy functional of the homogenized model directly, including boundary terms that have not been discussed in previous work. Our method is formulated in a generic setting which makes it applicable to a variety of geometries in dimension 1, 2 or 3. An implementation using a symbolic calculation language is proposed and it is distributed as an open-source library. Simple illustrations to elastic trusses having pre-stress or graded elastic properties are presented. The approach is presented in the context of discrete elastic structures and the connection with previous work on the higher-order homogenization of periodic continua is discussed.


Introduction
This work addresses the growing need for effective models capturing accurately the mechanical response of architected materials produced, e.g., by additive manufacturing techniques.In these materials, finite-size effects are not captured by standard (Cauchy) continuum models.Generalized continua theories, including higher-order (gradient) terms or micropolar fields have been introduced to overcome this limitation.Among various strategies proposed to obtain such models, asymptotic analysis offers a rigorous and fully predictive tool to derive exact higher-order or micropolar contributions without resorting to ad-hoc kinematic assumptions (Boutin 2019).
Asymptotic homogenization is a well-established technique for either discrete or continuous periodic microstructures.It aims at deriving an equivalent macroscopic set of equilibrium equations by means of a formal two-scale expansion, see (Sanchez-Palencia 1980;Bakhvalov and Panasenko 1989;Cioranescu and Donato 1999;Cioranescu and Paulin 1999) among others.The results of classical periodic homogenization are recovered at leading order.Later contributions focussed on pushing the asymptotic expansion to higher-orders, with the aim of deriving gradient contributions, either as small correctors to the leading-order prediction (Gambin and Kröner 1989;Boutin 1996;Smyshlyaev and Cherednichenko 2000;Hans and Boutin 2008;Bacigalupo 2014;Le and Marigo 2018;Abali and Barchiesi 2021) or as a non-local leading-order model when the standard equivalent medium is degenerate (Boutin and Soubestre 2011;Abdoul-Anziz and Seppecher 2018a;Abdoul-Anziz and Seppecher 2018b;Abdoul-Anziz et al. 2019;Durand et al. 2022).Asymptotic analysis has also been used to derive micropolar effective theories (Bažant and Christensen 1972;Dos Reis and Ganghoffer 2012;Nassar et al. 2020).
In spite of a large body of existing work, some aspects of higher-order homogenization remain elusive.Besides the question of applicable boundary conditions which is largely overlooked in existing work, most contributions typically start from the strong form of the equilibrium equations thus loosing the variational structure of the initial problem (some try to reconstruct the effective energy a posteriori, but systematically ignore boundary terms).Besides, natural extensions such as pre-stress or slow spatial variations of the elastic or geometric properties of the microstructure are rarely addressed, especially in higher-order contributions.
In this paper, we adapt our earlier work on higher-order, nonlinear asymptotic dimension reduction (Lestringant and Audoly 2020) to linear homogenization of discrete elastic structures, such as elastic trusses or networks of elastic beams.Based on formal arguments, we propose a homogenization method that is directly applicable to any given discrete microstructure (the connection with the homogenization of periodic continua is discussed in Section 7. We carry out homogenization at the energy level, thereby following the pioneering work of (Le and Marigo 2018).As a benefit, important simplifications are made during the homogenization procedure, allowing the higher-order homogenized model to be obtained in compact form via a formal, systematic expansion.In addition, our method is versatile and allows for natural extensions such as the case of pre-stressed structures.We also keep track of boundary terms that have been ignored in previous work.Lastly, our method is implemented in a symbolic calculation language and distributed as an open-source library named shoal, for Second-order HOmogenization Automated in a Library (Audoly 2023).
Our method starts from an energy formulation derived by applying a continuum approximation to a linear discrete microstructure.This preliminary continualization step is not strictly part of this paper; it is briefly illustrated based on an example.We apply to this continualized energy a two-scale expansion, assuming slow variations for macroscopic quantities such as strains and lattice properties, and rapid variations for the degrees of freedom at the microstructural scale.
Our method aims at condensing out these rapidly varying fields, by solving an energy stationarity problem order by order.Constraints are included from the onset in the energy formulation, which helps making the derivation compact.The procedure involves solving a set of elementary linear-algebra problems, and is implemented in a formal calculation language.
2 Input to the homogenization procedure In this section, we present the elastic model serving as the starting point of the homogenization procedure.It is formulated in a continuous domain Ω ⊂ R  of the -dimensional Euclidean space.The elastic model is specified in an abstract and generic form, making the homogenization procedure applicable to a broad range of situations including both periodic continua and discrete structures-such as one-dimensional (beam-like), two-dimensional (plate-like) and three-dimensional (bulk) elastic trusses.It is also versatile enough to handle structures possessing slowly modulated properties or pre-stress, as illustrated in Sections 4.1 and 4.2, respectively.
The task of casting a given problem into the generic form proposed in this section is not particularly difficult but has to be carried out on a case-by-case basis: this aspect is barely touched upon in Section 4 and will be illustrated in follow-up papers.When the structure under study is a discrete elastic truss, this preliminary step involves a so-called continualization assumption that turns the original discrete energy into the continuous energy functional used as a starting point in this paper, see Equation ( 4) and the discussion in Section 4.1.
We limit attention to linear homogenization problems.The extension to nonlinear homogenization can be done by adapting our previous work on nonlinear dimension reduction (Lestringant and Audoly 2020;Audoly and Lestringant 2021).

Energy formulation of the input model
We proceed to specify the continuous elastic model used as an input to the homogenization procedure.The presentation is intentionally abstract: illustrations will be provided in Section 4.
The model is formulated over a continuous domain Ω ⊂ R  , and we denote by  ∈ Ω the space variable, see Figure 1.A deformed configuration of the elastic body is parameterized by three vector fields ( ),  ( ) and ( ) defined over Ω: • microscopic degrees of freedom ( ) ∈ R   which we seek to eliminate using the homogenization procedure, B. Audoly and C. Lestringant An energy approach to higher-order homogenization Typical application of the method: homogenization of a beam lattice in a domain Ω ⊂ R  , here with  = 2.In this stylized representation of the beam lattice, the gradient of material parameters ∇ includes both gradients of elastic properties (varying thickness of beams) and geometry (mesh curvature).The mathematical domain Ω where we carry out homogenization is strictly included in the physical domain of the lattice to avoid boundary layers.The elastic lattice depicted above serves as an illustration: the homogenization method is presented in an abstract setting that does not make any reference to a lattice.
• macroscopic variables which are held fixed during homogenization, namely: the macroscopic strain  ( ) ∈ R   -variable material parameters ( ) ∈ R   .The integers   ,   and   are input parameters of the homogenization procedure.The goal of the procedure is to slave the microscopic degrees of freedom  to the macroscopic variables  and , thereby delivering a homogenized model depending on  and  only.The difference between the macroscopic variables  and  is that  captures the slowly variable properties of the elastic structure which are prescribed once for all, although  is considered fixed during the homogenization procedure but is actually an unknown of the structural problem that the homogenized energy helps solving.
Remark 1 We work in the discrete case, i.e., assume a finite number (  < ∞) of microscopic degrees of freedom.This is typically relevant to elastic truss structures, see Section 4. The extension of the method to infinitely many microscopic degrees of freedom (  = ∞), relevant to periodic continua, is discussed in Section 7.
Remark 2 We will present the homogenization method without providing definitions of the macroscopic strain , variable material parameters  and microscopic degrees of freedom .These quantities can be anything as long as they match the postulated forms of the strain, constraints and energy in Equations ( 2) to ( 5).This abstract presentation makes the method quite versatile, and opens up the way for extensions that we will cover in future work.Specific choices of ,  and  are proposed in the illustration examples (Section 4).
The model used on input may make use of constraints and we denote by  c ⩾ 0 the number of applicable (scalar) kinematic constraints.By convention, the left-hand sides of these kinematic constraints are incorporated into the microscopic strain vector  ∈ R   , and are extracted from  using an appropriate matrix Q ∈ T ( c ,  ) .In view of this, the kinematic constraints are written in the form Q •  (( );  ( ), ∇ ( ), . . .; ( ), ∇( ), . ..) = 0  c ∀ . (3) Section 4 provides a specific illustration on how  and Q can be set up to conform to Equation (3).
The constant tensor Q is provided as an input of the homogenization procedure, see Table 1.
In the input model, the strain energy is assumed to be of the form where the strain energy density in the bulk term is given in the context of linear elasticity by Explicit expressions of the strain , energy Φ and elasticity matrix K appearing in (4-5) will be provided in the illustration Section 4, see for example (33)(34) The elasticity tensor K () ∈ T (  ,  ) is provided as an input to the homogenization method in the form of a tensor-valued, symbolic function of , see Table 1.
The square brackets around the arguments of Φ[, ; ] in the left-hand side of Equation ( 4) denote a functional dependence: the strain energy Φ depends on the functions ,  and  over the entire domain Ω.
The list of parameters passed as an input to the homogenization procedure is recapitulated in stiffness matrix, see Equation ( 5) Table 1 List of parameters passed as an input to the homogenization procedure.The notation used in the last two columns is defined in Appendix A.
Remark 3 This formulation of the input model is designed to be versatile.For instance, the presence of pre-strain can be accommodated by adding a coefficient capturing the pre-strain intensity as an additional entry in the vector  (whose definition is up to the user), and by propagating it to  by an appropriate definition of   (), as illustrated in Section 4.2.Similarly, the presence of pre-stress can be accommodated by including a constant entry with value 1 in , propagating it to , and inserting the pre-stress into the corresponding rows and column of K.

Assumption of slow variations
One of the key assumptions of homogenization is that there is a separation of scales between a microscopic length ℓ (typically the spatial period of the underlying discrete lattice or periodic continuum) and the size  of the domain, ℓ ≪ , as sketched in Figure 1.The goal of homogenization is to deliver an effective model applicable at the macroscopic scale , by hiding the 'details' taking place at the microscopic scale ℓ.
Mathematically, this separation of scale is captured by the small parameter The various fields  ( ), such as  = ,  =  or  = , are assumed to depend on the variable  evolving on the slow scale  = ℓ/, implying that their successive gradients scale as In the following, the gradient ∇ = / will therefore be treated implicitly as a small quantity of order .This implicit notation has the benefit of avoiding a large number of predictable occurrences of the parameter , as discussed in Remark 4 below.The scaling assumption Equation ( 7) is not applicable in the layers that are present near the boundaries or near the point of application of point-like force.The homogenization domain Ω therefore needs to be slightly smaller than the actual physical domain of the elastic body, as shown in Figure 1 (see also Equation [1] in the work of (Abdoul-Anziz and Seppecher 2018b) for an accurate description of how Ω can be shrunk).Alternatively, boundary layers can be solved rigorously and represented in the homogenized model by means of effective boundary terms, see for example (David et al. 2012), but this is beyond the scope of the present work.

Remark 4
The scaling assumption Equation ( 7) can be motivated as follows.By convention, we consider the microscopic length ℓ to be ℓ = O (1) and the macroscopic length to be  = O ( −1 ).In our notation, any macroscopic field  such as ,  or  is implicitly a function of the slow variable  =  , i.e., what we write as  ( ) should be spelled out as  ( ) =  ( ), where  is a dimensionless function, independent of .The gradients can then be obtained as where denotes the gradient with respect to the slow variable: this implies the scaling assumption ∇   ( ) = O (  ) in Equation ( 7).The formal rule (7) dispenses with a notation for the slow variable.
Remark 5 The dependence on ( ) of the strain  in Equations (1) to (2) and of the energy density  in Equations ( 4) to (5) allows one to handle the case of structures whose elastic properties vary over the large scale  = ℓ/, as conveyed by the variations in thickness of the microstructure sketched in Figure 1.The definition of ( ) is entirely up to the user.For a 2D elastic truss possessing rotational symmetry, for instance, one could define ( ) = ( 2 1 +  2 2 ) 1/2 and   = 1 to capture the dependence of the local truss properties on the distance to the center of symmetry.In the case of variable properties without any particular symmetry, one should set ( ) =  and   = .When specifying the input model, one should ensure that any dependence of the properties of the elastic medium on the slow variable  takes place through the quantity ( ), as illustrated in Section 4. For structures having uniform properties over the large scale, one can ignore any dependence on  and set   = 0, see Appendix D.
B. Audoly and C. Lestringant An energy approach to higher-order homogenization 3 Summary of the main results

Homogenization as a partial energy relaxation
In this paper, we identify an equivalent continuum by making stationary the functional Φ in Equation (4) over the microscopic degrees of freedom ( ) for a fixed distribution of the macroscopic variables ( ) and  ( ).The stationary point is denoted as  ★ [, ], and will be assumed to be unique-as indicated by the square brackets, the stationary point  ★ is a functional of  and ; it is also a function of  , whose values are denoted as  ★ [, ] ( ).
In view of (2-3), the kinematic constraint can be written as (( )) : ∇ ( ) + . . .= 0 ∀ .It is treated using Lagrange multipliers ( ) ∈ R  c .The variational problem takes the following form: for given  and , we seek the solution (10) The main goal of this paper is to derive an explicit expression of the homogenized functional Φ ★ .The stationary point problem in (9) can be written formally as Equations ( 10) and ( 11) are at the heart of our variational approach to homogenization.The homogenization works under the assumption that the energy is positive-definite in the subspace of admissible microscopic degrees of freedom, i.e., As we will see, this is a necessary condition for the variational problem (9) to have a unique solution at leading order.It is also a sufficient condition for the homogenization procedure to produce a result up to second-order.

Homogenization results in compact form
The variational problem ( 9) is impossible to solve in closed form in general but thanks to the assumption of scale separation (Section 2.2), we can derive the following approximation of In the boundary terms ∮ Ω . . .d in the second line, Ω denotes the boundary of the domain,  is the unit outward normal, and d the area (if  = 3) or the length (if  = 2) of a boundary element, see Figure 1.The typographical variant of the integral sign ∮ will be used throughout for boundary integrals.
The dimension and symmetries of the homogenized tensors  , , ,  and  are specified in Table 2. Our main result is to derive their expansion in the successive gradients The tensors   (),   (),  0 (),  1 () and  0 () appearing in the right-hand sides of ( 14) are obtained in explicit form in terms of the local material parameters  in Appendix B, see (B.12), and in Appendix C, see Appendices C.5 and C.10.Their properties are listed in Table 3.This completes the specification of the homogenized model (13) up to order  2 included.
Table 2 Dimensions and symmetries of the tensors appearing in the homogenized energy functional in (13).
tensor space symmetry The expansion (13-14) is established in Section 5 by solving the variational problem (9) for  order by order in .The solution is found in the form where  [] and  ′ [] are given as expansions in the successive gradients of , B. Audoly and C. Lestringant An energy approach to higher-order homogenization tensor space symmetry usage Table 3 Tensors delivered by the homogenization procedure, defining the homogenized model in Equations ( 13) and ( 14).
and the localization tensors  0 (),  1 () and  ′ 0 () are derived in explicit form in terms of the variable material parameters  in the Appendix, see (B.8) and (C.48).

Homogenization results in the form of a systematic expansion
The various contributions to Φ ★ in (13) can be grouped order by order as follows, by inserting ( 14) into (13) and using (7): • The leading-order contribution Φ ★ [0] =  (   0 ) is given by and characterizes an equivalent Cauchy medium through a homogenized stiffness tensor  0 () depending only on the local material parameters : this homogenized stiffness  0 () matches that predicted by classical homogenization.• The first correction Φ ★ [1] =  (   1 ) is given by • The second correction The homogenized energy Φ ★ [, ] in (13) is nothing but the sum and it is asymptotically exact up to a higher-order contribution Φ ★ [3] = O (   3 ) which we do not attempt to resolve.

Remark 6
The actual derivation of the homogenized model proceeds in the reverse order than the high-level presentation above: the order-by-order expansion (17-20) is derived first, and the compact form (13-14) is obtained next by rearranging the terms.
The solution for  in (15-16) is derived based on the assumption that the microscopic variables  =  ★ [, ] can be expanded in powers of , where  [ ] ( ) = O (  ) denotes the contribution of order   to .Specifically, the microscopic solution which yields (15-16) by rearranging the terms.
Remark 7 As discussed in Section 2.2, there are implicit scaling factors   in all our formulas.Their consistency can be checked as follows.Take Equation (22b), for instance: the subscript '[1]' in the left-hand side indicates that this is a quantity of order ; this is consistent with the fact that the right-hand side is homogeneous of degree 1 with respect to the symbol ∇ = O ().When checking homogeneity, the boundary terms must be treated with special care: in Equation ( 19), for instance, the integrand of the bulk integral is as quantity of order  2 , in line with the subscript '[2]' appearing in the left-hand.The integrand of the boundary integral is however a quantity of order  1 ; the paradox is resolved by noting that the measure of the domain O (  ) for the bulk integral, but O ( −1 ) = O (  ) for the boundary integral-recall that ℓ = O (1) and  = O ( −1 ).
Ultimately, both integrals are of order O (   2 ).

Illustrations
In this section, we provide simple illustrations of the homogenization method.Equivalent high-order beam models are derived for various truss lattices, in the same line as a number of earlier works on periodic 1D structures, including (Hans and Boutin 2008;Abdoul-Anziz and Seppecher 2018b).This is not a fundamentally new contribution, our main goal being to illustrate how the abstract method can be applied to specific problems.The first two examples demonstrate that the homogenization can naturally handle elastic structures possessing graded properties (Section 4.1) and pre-strain (Section 4.2), two features that are not commonly addressed in the literature.A variant of this truss, this time including rigid bars arranged in a way that the macroscopic strain is constrained, is proposed in the Appendix (Appendix E.9).We also demonstrate how the method can be extended to a frame made up of beams (rather than springs), and show that the homogenization method can be adapted to deliver a Timoshenko beam model (Section 4.3).

A truss lattice having slowly variable elastic properties
We consider a truss lattice comprising elastic bars connected by perfect hinges, arranged in rectangular geometry made up of square cells with side length , see Figure 2. Discrete model.Let ( 1 ,  2 ) be a frame of orthonormal directors with  1 aligned with the longitudinal direction of the truss.The nodes are labelled with indices  ± = (, ±) where  in an integer spanning the longitudinal direction, and the symbol ± is a discrete transverse coordinate labelling the lower versus upper layer of nodes.We denote by   ± = S  1 ±  2  2 with S =  the (undeformed) position of node  ± , and by   ± its infinitesimal displacement.
Given a bar labelled by , we assign to it an orientation and denote by  1  and  2  its ordered end-nodes in reference configuration, by   = ∥  2  −   1  ∥ its undeformed length and by   = (  2  −   1  )/  its undeformed unit tangent.The (dimensionless) axial strain in the bar is given by The discrete elastic energy in the lattice is written in the form The lattice has graded properties as the elastic constant  of the bars depends on the midpoint coordinate S  = 1 2 ( S 1  + S 2  ).In order to keep the homogenization results as simple as possible, we make the simplifying assumption that the different types of bars have identical elastic constants : the more natural assumption that all bars have identical cross-sections (and thus that their elastic constants  is proportional to the length   ) could be addressed by making  a function of not only S  but also of the type of the bar , which does not raise any particular difficulty.
Continualization, scaling assumptions Our continuous model is one-dimensional ( = 1) and involves macroscopic fields that are functions of the longitudinal coordinate S. As part of the continualization step, the nodal displacement   ± is sought in terms of continuous fields in the longitudinal direction as where  ( S) and  ( S) denote the macroscopic longitudinal and transverse displacement of the equivalent rod, respectively, and  ±  ( S) are the components of the microscopic displacement for either row of nodes (±).The special case  ± 1 =  ± 2 = 0 corresponds to the (asymptotically incorrect) assumption of an unshearable model having rigid cross-sections-note that the term ∓ ′ /2 represents the rigid rotation of the cross-section imposed by the centerline.We do not impose  ± 1 and  ± 2 to be zero.We impose the kinematic constraint where ⟨•⟩ denotes the average over the top and bottom rows, i.e., ⟨ ±  ⟩ = 1 2 ( −  +  +  ).This warrants that ( ( S),  ( S)) capture the average nodal displacement at coordinate : indeed, by combining (25-26), we have ⟨  ± ⟩ =  ( S ) 1 +  ( S ) 2 .As a result, the equivalent rod passes through the midpoints of the sections [  −   + ].This choice is somewhat arbitrary: by using different weights in the average (26), we could introduce a lateral offset in the definition of the centerline.
We define the macroscopic length to be /, with  ≪ 1 the scale separation parameter; see Equation ( 6).A standard scaling analysis yields the macroscopic stretching strain as  ∼  ′ ( S), the macroscopic rotation as  ( S) ∼  ′ ( S), the relative rotation between successive transverse links as  ′ ( S), and thus the differential stretching strain of bars located on the inner and outer sides (curvature effect) as  ∼  ′ ( S). Natural scales are found by balancing the two sources of stretching strain , yielding  ∼  ′ ∼  ′ ∼  ′′ , which suggests introducing dimensionless B. Audoly and C. Lestringant An energy approach to higher-order homogenization quantities in the form where ,  and  ±  are dimensionless unknowns and  = S/(/) = S/ is a slow variable, i.e., the arclength scaled by the macroscopic length /.
In addition, we assume that the variations of elastic properties take place on the macroscopic scale, i.e., In what follows, we eliminate the original quantities in favor of the dimensionless ones,  (),  (),  ±  () and  () everywhere.By contrast with the rest of the paper, we keep explicit track of the small coefficient  in the present section.
Setting up the input of the homogenization procedure In the classical theory of linear, planar beams, the two relevant strain measures are the stretching strain  () :=  ′ () and the bending strain  () :=  ′′ ().We thus anticipate that the homogenized energy will depend on the macroscopic strain  (  = 2) which we define as In terms of the unscaled displacement ( ,  ), the dimensionless strain measures are given by  () =  ′ (/) and  () =  ′′ (/).
In addition, we define the vector of microscopic degrees of freedom (  = 4) as These quantities  and  were chosen in such a way that the strain of a bar  can be expressed in terms of ,  and their successive gradients at the midpoint coordinate    , see Equation (32) below.
In our discrete model ( 24), the elastic constants  (   ) depend on the midpoint coordinate    .This makes the energy density  in (34) depend explicitly on .Since  is required by design to depend on  and  only, see (5), we pack up the coordinate  into the list of material parameters (), An alternative (and ultimately equivalent) approach would be to define () as  ().
Next, we define the strain  as the concatenation of (i) the discrete strains   =   given in ( 23), in each of the 5 types  of bars that make up the lattice, see Figure 2, together with (ii) the left-hand sides of the two kinematic constraints appearing in (26), We therefore have   = 7 strain variables.
Having included the left-hand sides of the kinematic constraint (26) at positions 6 and 7 in , we can easily express the  c = 2 constraints in the form Q •  = 0 expected in Equation (3), by defining the constraint extraction matrix as Q = 0 2×5 1 2×2 using block-matrix notation.
Inserting the expression (25) of the nodal displacement in the expression of , using the scaled quantities introduced in ( 27), and performing Taylor expansions about the midpoints of the bars, we get which is of the form  =  (; ,  ′ ,  ′′ , . . .; ,  ′ ,  ′′ , . ..) expected in Equation (2).In ( 33), the argument of the functions ,  ′ ,  ′′ , ,  ′ ,  ′′ is implicitly assumed to be the midpoint    of each bar.The details of the calculations can be found in the companion Mathematica notebook.
The discrete energy ( 24) is finally continualized in the canonical form (4-5) as The notation  ( 1 ) conveys the fact that the argument  1 =  to be passed to the stiffness distribution  () is the first (and only) component  1 of , see (31).In (34), the discrete energy (24) has been continualized by using the formal rule where  is an index running over the 5 different types of links,  ∈  is included in the partial sum of all links  belonging to a particular family , and   denotes the expression of the energy   relevant to the links  belonging to a particular family .The coefficient 1/ appearing in the definition of  in ( 34) is nothing but the lineic density of links of each type per unit dimensionless length .Thanks to the assumed periodicity of the lattice, we have been able to rewrite the discrete sum in ( 24) into an integral in (34).
Remark 8 In the argument sketched above, the continualization is based on the formal approximation rule for a discrete sum in terms of an integral,  ∈   ≈ ∫   d  , which has been used in a number of earlier work.A rigorous justification of this approximation is known as the Euler-MacLaurin formula: it is correct to order  2 in an infinite domain, but needs to be corrected by boundary terms in a finite domain.
At this point, we can identify the quantities that are required on input of the homogenization method, as specified in Table 1: they are listed in Table 4.The integer constants appearing in the left column of Table 4 and the constraint extraction matrix Q have been collected from the above discussion.The tensors   , ...,  ′′  collect the numerical coefficients appearing in (33) and are found by identification with Equation (2).The elastic stiffness tensor K () is found by identifying (34) with (5).Expressions in Table 4 make use of the notation    = (0 1 . . .0 −1 1  0 +1 . . .0  ) ∈ R  for the -th Kronecker vector with length ; for example,  5 2 = (0 1 0 0 0).

B. Audoly and C. Lestringant
An energy approach to higher-order homogenization

Homogenization results
We propose an implementation of the general homogenization method described in this paper (Sections 5 and 6) in the form of an open-source library named shoal (Audoly 2023).The quantities listed in Table 4 are passed to the library in symbolic form, see the input file inhomogeneous-truss.nb1 .The code executes without any further input from the user.It delivers the tensors listed in Table 3, which we now proceed to interpret using (13-16) and (29-30).We note that ∇ =  1 1 ⊗  1 1 in view of the definition (31) of .The homogenization method returns In view of (15-16), the microscopic displacement  =  0 •  +  ′ 0 : ∇ + ( 1 : ∇) •  + • • • is found with the help of (29-30) as The quantities ± ()/6 of order  0 in the right-hand side represent the predictions of classical (leading-order) homogenization.Note that the correction of order  1 includes not only a contribution proportional to the strain gradient  ′ but also one proportional to the gradient  ′ of elastic properties-this corresponds, respectively, to the ∇ and ∇ contributions appearing in the right-hand side of (22b).
The homogenized energy functional is obtained by interpreting the output of the code similarly using Table 3, see also (13)(14), The detailed expressions of the tensors  0 ,  1 , ... underlying the above expression of Φ ★ can be found in the Mathematica notebook inhomogeneous-truss.nb included in the library.The energy Φ ★ in (38) depends on the stretching strain  (), on the dimensionless curvature strain  (), as well as on their gradients  ′ () and  ′ (), and on the gradients of elastic properties  ′ () (the dependence of ,  and  on  implicit in the above integral).This result is valid for any given slowly varying distribution of spring stiffness  () since  () is treated as a symbolic function.As we consider an infinite structure, we have ignored in (38) the boundary terms that are returned by the homogenization library.
In view of the negative coefficients −1/6 and −11/24 appearing in the last term in the integral, the energy Φ ★ can be made arbitrarily large and negative by incorporating oscillations into the unknowns  () and  (), having both small amplitude and small wavelength.This holds even in the simple case where the elastic properties are uniform ( ′ = 0).This points to the lack of lower semi-continuity of the homogenized energy Φ ★ , an undesirable property of higher-order gradient models that has been documented by several authors for various elastic structures, see for instance (Le and Marigo 2018).It calls for a regularization of the functional Φ ★ , a point which we address in future work.
Remark 9 Alternatively, the homogenized functional Φ ★ can be expressed in terms of the original unscaled quantities, namely the elastic constant  ( S), the axial strain  ( S) =  ′ ( S) =  S/ and the curvature  ( S) =  ′′ ( S) = 1   S/ .The result is free of the parameter , 4.2 A truss lattice with pre-strain In order to illustrate the ability of our method to handle pre-strain, we now include an additional tensile pre-strain + ()/2 in the lower layer of bars, and a contractile pre-strain − ()/2 in the upper layer of bars.For the sake of simplicity, we focus on the case of uniform elastic properties across the length, taking  to be independent of .Accordingly, we change the strain definition for the two diagonal bars to  2 =  0 2 −  ()/2 and  4 =  0 4 − (− ()/2), where the quantities  0  refer to the expressions of   in the absence of pre-strain given in the right-hand sides of (33) (Section 4.1).
Since  is required to be a function of  and  and their derivatives, we add a third component  3 =  to the macroscopic strain  and rewrite  2 =  0 2 −  3 ()/2 and  4 =  0 4 +  3 ()/2.Note that the pre-strain is subtracted from the original element strain, so that the energy in the bar labelled 2 in Figure 2, for instance, is given by  2 ( 0 2 −  ()/2) 2 and is minimum when  0 2 =  ()/2.To sum up, we make the following changes in the specification of the input problem: The tensors   ,  ′  and  ′′  are unchanged.The dimension of the tensors   ,  ′  and  ′′  is increased to reflect the new dimension   = 3 of the macroscopic strain vector .This change of dimension takes place simply by adding zero entries to  ′  and  ′′  , and by including a new contribution   = • • • + 1 2 (− 7 2 +  7 4 ) ⊗  3 3 to   , to reflect the new  3 -terms in  2 and  4 .The homogenization code is run again with the modified set of input parameters, see the notebook prestrained-truss.nbincluded in the library.The homogenized energy is obtained as In the leading-order term, the pre-strain sets the natural curvature of the equivalent beam as  0 () =  (), as could be anticipated.Note that the gradient of pre-strain  ′ () contributes to the energy at order  2 .
B. Audoly and C. Lestringant An energy approach to higher-order homogenization Figure 3 An elastic frame made up of identical beams having length , stretching modulus  and bending modulus  .(a) General view, (b) the three modes of deformation of a linear beam, illustrated here for a beam  in the family  = 1: stretching   , bending   and shearing   .

Homogenizing an elastic frame into a Timoshenko model
We now analyze the frame made up of linear elastic beams, shown in Figure 3. Specifically, we show that the homogenization procedure can deliver a Timoshenko beam model, when the rotation of the cross-sections is added to the list of macroscopic parameter , implying that this parameter is kept fixed during the energy relaxation.This example has been treated in (Abdoul-Anziz and Seppecher 2018b) and we use their results to verify our method.
Changes to the homogenization procedure To address the elastic frame, we start over from the model in Section 4.1, with the following changes.• There are 3 families of beams (and not 5) as indicated by the labels  ∈ {1, 2, 3} in Figure 3(a).
• We revert to the case of uniform properties by discarding the parameter .
• We change the scaling assumption on the transverse displacement  in (27b) to in such a way that the rotation  ′ (and not the curvature  ′′ ) is of order  0 : this delivers the Timoshenko model is a more convenient form.• To model the stiff junctions, we introduce the nodal rotation  ± at a node  ± (in addition to the nodal displacement   ± from (25)):  ± is given in terms of two additional unknown continuous functions  −  ( S) and  +  ( S) as Including the centerline rotation  ′ as the first term in the right-hand side makes the microscopic displacements  ±  invariant by rigid-body rotations, so that they can be expressed in terms of the strain as assumed in our procedure.
• For a beam of type , the classical theory of linear (Euler-Bernoulli) beams can be summarized as follows.In addition to the stretching strain   in ( 23), one has to consider a curvature strain   and a shearing strain   , as sketched in Figure 3(b), Here,   is an apparent shearing strain at the scale of unit cells, which actually resolves into bending microscopically, see Figure 3(b): our Euler-Bernoulli 'microscopic' beam model is shearless.The elastic energy of the beam is then of the form • We define the (apparent) rotation of the cross-sections  as the rotation of the line passing through the nodes  − and  + facing each other, • Next, we consider the following scaled macroscopic parameters: centerline extensional strain  () =  ′ () =  ′ (/), centerline curvature  () =  ′′ () =    ′′ (/), rotation of cross-sections  () =  ′ () − ( +  () −  −  ()) =  (/) and shear angle • The vector of microscopic degrees of freedom  is extended to reflect the presence of the two additional fields capturing nodal rotation: Equation ( 30) is changed to • A key modification to the homogenization procedure is that we include the shear angle in the list of macroscopic parameters, by contrast with (29).Equation ( 47) can then be rewritten as This kinematic constraint is taken care of by including the left-hand side of the above equation as a component of strain (specifically,  12 ) and by extending the constraint matrix Q from Table 4 to include a third row filled with zeros, except for an entry equal 1 in column 12, see Equation ( 3).• The strain  ∈ R 12 is now a vector of length 12, made up of the three strain measures (  ,   ,   ) for each of the three types of beams,  ∈ {1, 2, 3}, the two left-hand sides in the zero-averagedisplacement constraints (26), and the left-hand side of (50).• The discrete energy is finally approximated by a continuous integral Φ = ∫ +∞ −∞  d capturing the elastic potentials of the beams ( 45) The coefficient  in the denominators is produced when the discrete sum over beams is approximated by an integral, as earlier in (34).
The scale separation parameter  ≪ 1 is unspecified so far, and we set it to be the small aspect-ratio of the beams, i.e., we set Doing so, we are anticipating that, in the forthcoming Timoshenko model, the shearing mode will relax over a typical macroscopic length / (this macroscopic length / has been introduced earlier in the definition of the scaled arclength S = /).The validity of this assumption will be checked later, by verifying that the the various terms in the Timoshenko model are produced with consistent powers of , see Equation ( 55) below.Equation ( 52) is used to eliminate  in favor of  and  in the elastic potential (51): in the homogenization procedure, the energy Φ provided on input may depend explicitly on .
Results of the homogenization procedure The homogenization procedure is carried out in the Mathematica notebook frame-timoshenko.nbincluded in the library.With the boundary terms omitted (case of an infinitely long truss), the result of the homogenization procedure are interpreted as where  ★ ext (,  ′ ) is a higher-order gradient bar model applicable to the longitudinal displacement , and  ★ Tm (, ,  ′ ) is a Timoshenko beam model applicable to the transverse displacement , The details can be found in the Mathematica notebook.
Remark 10 The square term  2 ( +  ′ ) 2 appearing in  ★ Tm is first obtained in expanded form, with the different terms  2  2 , 2 2  ′ and  2  ′2 appearing respectively in the successive contributions to the energy.The factored form in the equation above emerges when , the stretching energy is dominant: this agrees with the fact that Abdoul-Anziz and Seppecher (Abdoul-Anziz and Seppecher 2018b) reported an inextensible model ( ≡ 0).
Turning therefore attention to the transverse displacement, we eliminate the shear angle  in favor of the apparent rotation  =  ′ +  of cross-sections.Recalling the definition of the scaled curvature  =  ′′ , and rearranging using (52), we obtain the Timoshenko model in standard form, and the elastic moduli are identified as The same elastic frame has been homogenized in (Abdoul-Anziz and Seppecher 2018b) and the authors reported homogenized moduli   = 1/2 and   = 2.These values are a special case of (57), implying that our homogenization results are consistent.Indeed, a careful analysis of the beam model used by these authors reveals that they limited attention to the special case  = 1/ and  = /4 (when setting both their elastic moduli  AS and  AS to 1, see the note frame-timoshenko-AS18-beam-model.pdf in the library).
Remark 11 In terms of the original (un-scaled) deflection  ( S) and cross-section rotation  ( S) =  S/ , the Timoshenko model appearing in (53-56) can be rewritten as (58)

Remark 12
We obtained a Timoshenko model as a consequence of the fact that the shear angle  =  3 is treated as a macroscopic parameter, and not just because the lattice is made up of beams.Indeed, lattices made up of bars ( = 0) can yield Timoshenko models as well when homogenized.The motivation for using beams in this example was to compare with the earlier work of (Abdoul-Anziz and Seppecher 2018b).
5 Derivation of the homogenized model
Inserting ( 64) into (61), we derive in Appendix B the dominant contribution to the energy ] that was announced in (17), namely The expression of the leading-order stiffness tensor  0 () is given in Equation (B.12) in the Appendix.

Analysis of the gradient effect
We now proceed to solve the next orders in the microscopic displacement, see ( 21) and (64), Inserting this into (4), we derive a Taylor expansion of the energy as where Φ [ ] denotes the term of order   , and the star is used to mark energy contributions that do not depend on the yet-unknown corrector  [1] .The dominant term Φ ★ [0] is the functional found earlier in (65), while the next-order terms Φ [1] and Φ [2] are obtained in Equations (C.20) and (C.24) in Appendix C as B. Audoly and C. Lestringant An energy approach to higher-order homogenization 6 Symbolic implementation: The shoal library We have implemented the method in the symbolic calculation language Wolfram Mathematica (Wolfram Research, Inc. 2021), and we distribute it as an open-source library named shoal (for Second-order HOmogenization Automated using a Library) (Audoly 2023).
The implementation makes use of standard tensor algebra operations on symbolic tensors, including general transpositions and multiple contractions, as well as symbolic differentiation with respect to , see (C.28).Note that the vector  never appears explicitly in the implementation.
At leading order and in the non-deficient case, the procedure is implemented by the equations listed in Table 5, corresponding to the first three bullet points above.The homogenization at the two following orders makes use of the subsequent bullet points.The special case of uniform properties, when the parameter  is absent, is worked out in Appendix D, see Table D.2 in particular.
Table 5 Implementation of the leading-order procedure in the non-deficient case ( d = 0), based on the formulas referenced in the first three items in the bullet list from Section 6.
7 Connection with the second-order homogenization of periodic continua In this section, we show that the extension of the proposed homogenization method to the case of periodic elastic continuum gives similar results to the classical approach (Smyshlyaev and Cherednichenko 2000;Boutin 2019;Durand et al. 2022).Similarly to the discrete example discussed earlier, we consider an infinite medium, so that boundary terms do not matter.We consider a periodic microstructure: the extension to slowly varying material properties or geometry could be considered as well, and would yield a derivation similar to the one proposed in (Le and Marigo 2018).
Figure 4 A continuum with a periodic microstructure.

𝒰
It is possible to derive the homogenization for periodic continua from the homogenization results for discrete systems established in the previous sections but this is quite technical.We thus prefer to carry out homogenization from scratch, following the same sequence of steps again.In addition, we do not attempt to link the various quantities relevant to the periodic continua to those relevant to discrete systems in a detailed way: we will content ourselves with formal analogies between them.

Canonical form
We consider a linear elastic continuum of dimension  = 2 or  = 3 with a periodic microstructure of unit cell U, as illustrated in Figure 4 for  = 2.The displacement field in the composite depends on a slow coordinate  and a fast coordinate , with  =  where  is a small parameter, see ( 6), as where (, ) ∈ R  denotes a rapidly fluctuating field, U-periodic with respect to its second argument and subject to a zero-average constraint over one unit-cell, This constraint ensures that  ( ) captures the average displacement, ⟨ (, •)⟩ =  ( ).We define the macroscopic strain  as the symmetrized gradient of  ( ) where the symbol ∇ s denotes the symmetrized gradient with respect to the slow coordinate  .The microscopic strain is the pair  =  (, ), ⟨(, •)⟩ with  ∈ T (,) , where  s denotes the symmetrized gradient with respect to the fast coordinate .This definition is formally similar to (2): we can identify   •  ∼  ( ), 0 ,   •  ∼  s (, ), 0 and   ′ • ∇ ∼ ∇ s (, ), 0 and cast the constraint (75) in a form similar to (3): The energy is postulated in the form and where  (, ) = 1, 0 •  (, ).In (78), C() ∈ T (,,,) is a U-periodic elasticity tensor (slowly varying properties could be considered here by including an additional dependency on  , i.e., by considering C(, ), see (Le and Marigo 2018)).Note that no boundary integral is present in the above expression of Φ. Boundary integrals do emerge, however, when one attempts to justify rigorously the hierarchy in the integrations postulated in (78)-over  first and over  next-, starting from the elastic energy of a periodic continuum in the limit of scale separation  → 0. In addition, boundary terms make it possible to account for the incomplete unit cells located at the boundary of the structure.The analysis of these boundary terms is however beyond the scope of this paper.

Leading order homogenization
At leading order, we neglect slow gradients by setting ∇ s (, ) = 0, the energy writes in the form (61) with For a given symmetric strain tensor  ∈ T (,) , we seek the stationary point () =  [0] () subject to the constraint ⟨⟩ = 0.This writes, using Lagrange multiplier where () is a test function.This problem is similar to (63).
Replacing by a constant  shows that  = 0. Integration by parts provides the strong form of equilibrium.
where we have defined the microscopic stress  (0) (, , ) = C() : ( +  s ()) . (82) Note that  (0) is automatically periodic in , which warrants the equilibrium at the boundary of a cell U.In (81), we use the notation div  for the divergence with respect to the microscopic (fast) coordinate .

B. Audoly and C. Lestringant
An energy approach to higher-order homogenization

Gradient correction
As done earlier in (C.12), we further expand the microscopic strain  order by order: with Following (C.17), (C.18), we eliminate  [1] from the energy at order 1 using ( 84) and obtain In (Durand et al. 2022), odd-order tensors are considered to be zero due to centro-symmetry assumption: here we work in a more general setting.A similar order-1 energy contribution is derived in (Le and Marigo 2018), see their equation ( 36), with an additional contribution coming from the macroscopic gradient of elastic properties.The energy at order 2 further writes, see (C.17), Using relation ( 84), we can eliminate  [2] from the second term Next, we identify the terms (C.26) and (.27) that require an integration by parts After integration by parts, the energy (89) writes where  is the unit normal to the domain boundary.In (92), we use the notation div  for the divergence with respect to the slow coordinate  .The domain being infinite, we can ignore the boundary term and the associated stationarity condition on the boundary.
The corrector  =  [1] (∇, ) is eventually found as a solution to the local variational problem, similar to (72), where we have identified the first order correction to the microscopic stress  (1) (∇, , , ) = C() : ∇ s ( 0 () : ) +  s ().After integrating by parts, we obtain The microscopic stress  (1) is automatically periodic in , which warrants the equilibrium at the boundary of a cell U.This microscopic problem ( 94) is similar to equation ( 16) in (Durand et al. 2022).
The correction  ★ [1] (∇, ) to the microscopic degrees of freedom writes, following (73),  ★ [1] (∇, ) =  ′ 0 () ∴ ∇.We can thus write the first order correction to strain as  [1] =  1 (, ∇, ) where Inserting into (92) allows us to identify the second order contribution to the homogenized energy in the form (19) as with Boundary integrals have been removed in ( 96) due to the fact that we consider an infinite domain.This result is similar to equation ( 35) in (Durand et al. 2022) and equation ( 50) in (Le and Marigo 2018) (this latter work reports additional terms that capture the effect of a gradient of elastic properties).Note that it is possible to identify the tensors  0 ,  0 and  0 appearing in (13-14), by factoring out , ∇ and  in equations ( 85), ( 88), ( 96) and ( 97).However, this introduces cumbersome expressions and we prefer to ignore this step.

Discussion and conclusion
We have proposed an asymptotically exact, second-order homogenization procedure for linear, discrete elastic structures, such as elastic trusses or networks of elastic beams.Our homogenization method works at the energy level, see Equations ( 10) and ( 11), and is similar in this respect to the approach to dimension reduction developed by (Berdichevskii 1981;Hodges 2006;Lestringant and Audoly 2020).It uses an abstract energy formulation as a starting point, see Section 2, and as a result can equally cover two-dimensional or three-dimensional lattices made up of beams or springs.It is designed to be generic, and addresses the case of pre-stress or pre-strain as well as slowly modulated elastic or geometric properties-in their work on the homogenization of periodic continua, (Le and Marigo 2018) address the case of slowly modulated elastic properties but assume that the geometry of the unit cell is invariant.The method can also account for kinematic constraints, and can be readily applied to, e.g., a lattice made up of inextensible beams, see Appendix E. Besides, the connection with existing approaches on the homogenization of periodic continua has been pointed out in Section 7.
The homogenization procedure involves a series of linear algebra calculations that have been implemented once for all in the form of an open-source library, named shoal and based on the symbolic calculation language Wolfram Mathematica.This work will serve as a foundation for a series of applications which we will cover in follow-up papers.In one particular follow-up paper, we will analyze more advanced beam lattices than those treated in the illustration section (Section 4).
The uniqueness of the solutions  ★ [0] and  ★  [1] to the variational problems (B.4) and (C.33) follows from the assumption (12) that W  () is positive definite on the subspace of kinematically admissible microscopic degrees of freedom  defined by Q •   (( )) •  = 0.In these circumstances,  ★ [0] and  ★ [1] do not only make the strain energy stationary, as mentioned in the paper, but also minimum.
Classical work on homogenization postulates the displacement in the form of an expansion  (, ) =      (, ), where (, ) are the slow and fast variables, respectively, and identifies a potential Φ★ [⟨ 0 ⟩, ⟨ 1 ⟩, . ..] depending separately on each one of the macroscopic averages ⟨  ⟩ = ⟨  (, •)⟩, see (Le and Marigo 2018) for instance.As a result, the stationarity condition for Φ★ takes the form of separate problems for ⟨ 0 ⟩, ⟨ 1 ⟩, etc.It has apparently not been appreciated that the homogenized energy has to be a functional warrants that its value is unaffected by a re-parameterization of  leaving the physical displacement  unchanged.Our approach works differently: with the notation of the continuous case from Section 7, we postulate  (, ) =  ( ) +      (, ), subject to the constraint ⟨  (, •)⟩ = 0 (for any  and any ), see (74-75), compute the strain  in terms of the (total) macroscopic displacement  = ⟨⟩, and obtain a homogenized energy Φ ★ [] depending on the total macroscopic strain  (and not on each of the various contributions   to the strain separately).This is simpler than the traditional approach, and ultimately equivalent.
Among the perspectives opened up by the present work, we can mention the extension to non-linear elastic structures (which can be addressed by adapting our previous work on non-linear dimension reduction (Lestringant and Audoly 2020)) and a careful treatment of the boundary layers (which we could incorporate into the homogenized model by means of effective boundary terms, along the lines of what has been done by (David et al. 2012)).
Transposing will allow us to reorder the indices of a tensor in any desired order.Suppose for instance that we wish to rewrite an expression    as the component  ′   of another tensor whose indices must appear in alphabetical order:  ′ is clearly a transpose of , and the permutation is found by noting that the levels (1, 2, 3, 4) in the original tensor , corresponding to the indices (, , , ), become respectively the levels (1, 3, 4, 2) = ( 1 ,  2 ,  3 ,  4 ) in  ′ .This yields (A.3) Index reordering using transposition will be routinely used in combination with contractions to remove indices in tensor algebra, as in     ′   =   1324 ::  ′ .The transpose   of a matrix  ∈ T ( 1 , 2 ) is a particular case of the generalized transpose, The symmetrization of a tensor  with respect to a pair of indices (, ) is denoted as     .For a tensor  of rank  = 4, for instance, the symmetrization with respect to the first and third indices is given by where   3214 is obtained from  by permuting the first and third levels.
A tensor invariant by a permutation of its levels  and  will be said to satisfy the    symmetry; for instance,   13 is  13 symmetric, by construction.More generally, a tensor will be said to satisfy the  {  } { } symmetry if it is symmetric by the combined permutation of indices  ↔  and  ↔ .For instance,  is  {23} {45} symmetric ⇔    =    . (A.5) The symmetrization with respect to pairs of indices works similarly, The composition of symmetrizations is represented by the symbol •.In Equation (C.7), for instance, it stands for Given a tensor  () ∈ T ( 1 , 2 ,...,  ) taking an argument  ∈ R   , we denote as d/d ∈ T ( 1 , 2 ,...,  ,  ) its gradient, By a standard convention, the index  corresponding to differentiation appears last in the gradient.
When the parameter  coincides with the spatial variable  ∈ R  , we use the nabla notation, ∇ = d d . (A.9) The alternate notation ∇ has the advantage of respecting the order of indices but is also less standard.
B Detailed analysis of leading order (classical homogenization) Inserting the microscopic strain  (0) (, , ) given in (59) into the expression (60) of the strain energy density  (0) (, ), we have Using the expression of  in ( 5) and expanding, we rewrite this in block-matrix notation as and the tensors W () ∈ T ( +  ,  +  ) , W  () ∈ T (  ,  ) , W  () ∈ T (, ) and W  () ∈ T (,) are given by Using (B.2), the optimality condition (63) for  ★ [0] and  ★ [0] takes the form where We focus on the case where  () is invertible: the non-invertible (rank-deficient) case is treated in Appendix E. The solution ( ★ [0] ( ),  ★ [0] ( )) is then found by inverting this linear system as In the code, we implemented a more general expression of R that applies to rank-deficient matrices, see Appendix E and Equation (E.19) in particular.Equation (B.6) matches the form of the solution  ★ [0] ( ) =  0 () •  announced in (64) and the localization tensor is identified as The solution for the Lagrange multipliers is B. Audoly and C. Lestringant An energy approach to higher-order homogenization We proceed to introduce important additional quantities that characterize the leading-order solution.
The strain  [0] (, ) =  (0) (, ,  =  0 () • ) is found using (59) as where the strain localization tensor  0 () is given by The leading-order strain energy  ★ [0] (, ) =  (0) (, ,  =  0 () • ) can then be written with the help of (B.2-B.3) and (B.11) in a form that matches that announced in (65), namely  ★ [0] (, ) = 1 2  •  0 () • , where the elasticity tensor  0 () characterizing the equivalent Cauchy-type elastic continuum at order  0 is identified as To complete the analysis of solutions at order  0 , we derive a useful identity that will help simplify the higher orders in the energy expansion.Inserting the solution  ★ [0] ( ) =  0 () •  and  ★ [0] ( ) =  0 () •  into the stationarity condition (B.4) and using (B.3) and identifying  0 from (B.11), we get •  appearing in the first equation can be identified as the leading-order stress, consisting of the elastic stress •  enforcing the constraint.We therefore introduce the stress localization tensor as and rewrite Equation (B.13), after simplification by the arbitrary factor , as ) is the principle of virtual work at leading order: multiplying by a virtual displacement  on the left-hand side and by  on the right-hand side, and rearranging, it takes the usual form ( 0 () • ) • (  () • ) = 0, where the left-hand side is the stress contracted with the virtual increment of strain.
C Detailed analysis of the gradient effect

B. Audoly and C. Lestringant
An energy approach to higher-order homogenization on can rewrite the definition of  in (C.1) as The converse (unpacking) operation is implemented as

C.2 Structure coefficients
The leading-order prediction (22a) for the microscopic degrees of freedom,  ★ [0] =  0 (( )) • ( ) can be expressed in terms of ( ) with the help of (C.4) as The successive gradients of ( ) in (66) can then be calculated as where the symbol O ( 3 ) stands for terms of order  3 and higher, such as By design, the tensor L 1 (), L 11 () and L 2 () capture the successive gradients of  ★ They are identified by differentiating (C.5) with respect to  , which yields Table C.1 lists the properties of all the tensors used in this appendix, starting with the tensors L, L 1 , L 11 and L 2 just defined.The symmetrization operations outside the square brackets in (C.7b-C.7c)are a matter of convention.They reflect the symmetries of the tensors with which the operators L are contracted.
The 'content' column in Table C.1 can be explained as follows.By design,  ★ [0] = L() •  can be 'unpacked' (i.e., expressed in terms of  and ) as  ★ [0] =  0 () • , which is a function of  contracted with : the dependence on  will be treated implicitly, and we express this by writing that the content of  ★ [0] = L() •  is , hence the symbol  appearing in the 'content' column for the row L. Similarly, ∇ ★ [0] = (L 1 () • ) : ∇ can be unpacked as this is the sum of two terms, one being a function of  contracted with  ⊗ ∇, the other one being a function of  contracted with ∇: this is conveyed by the content column in the table, which shows ( ⊗ ∇, ∇) for the row labelled L 1 () used for reconstructing ∇ ★ [0] = (L 1 () • ) : ∇.The point of the  notation is to deal in a simple way with the multiplicity of terms appearing in the last column of Table 1.The remainder of the appendix will make use of this higher-level  notation.On the other hand, we use Table 1 to keep track of the actual content of the various tensors.Now, we proceed to represent the strain in terms of  and its successive gradients.Inserting (C.6) into the strain expression  in (2), identifying  [0] (, ) using (B.10), and rearranging the other terms, we get  =  [0] (, ) + (J 1 (( )) • ( )) : ∇( ) + (J 11 (( )) • ( )) :: (∇( ) ⊗ ∇( )) where the so-called structure coefficients are identified by As indicated in Table 1, J 11 is  {23} {45} -symmetric: this is a consequence of the fact that L 11 is  {45} {67} -symmetric.The  34 -symmetry of J 2 can be justified by a similar argument.

C.3 Strain energy expansion in terms of corrective displacement
Inserting (C.11) into the energy (4), and using the quadratic expression (5) of the energy density, we obtain a Taylor expansion of the energy as • [0] (, ) + K () •  [0] (, ) • ( [1] +  [2] ) Grouping the terms in the integrand order by order, identifying the term of order  0 as 1 2  •  0 •  using (B.12), and the quantity K () • [0] (, ) = K () • [0] () • = ( 0 () − Q  • 0 ()) • by (B.10) and (B.14), we obtain the energy expansion as where Φ [ ] = O (  ) are the successive terms in the expansion, the leading term is the quantity Φ ★ [0] [, ] identified earlier in ( 65) and the higher-order terms are given by (C.16) The Q  terms appearing in both bulk integrals can be removed, as [ ] ) = 0 by (C.13).This yields Inserting the expression of  [1] from (C.12a) in Φ [1] and using the principle of virtual work at dominant order in (B.15a), we have By a similar argument, the expression of  [2] in (C.12b) can be inserted in the bulk integral appearing in Φ [2] , which shows that the term is zero.This yields d .C.5 Extraction of  0 and  1 As announced in Table C.1, the content of A() is  ⊗  ⊗ ∇ and  ⊗ ∇, which means that the right-hand side of Equation (C.20) can be unpacked using (C.3) as where the tensors  0 () and  1 () are extracted from A() as No other term can be present in the right-hand side of (C.22): a term such as  0 () :: ( ⊗  ⊗ ∇), for instance, would be inconsistent with the fact that the energy is homogeneous with degree 2 in  and ∇, see ( 2) and ( 5).
The properties of the tensors  0 () and  1 () are listed in Table 3.
The expression of the first correction Φ ★ [1] to the energy in (C.22) was announced in (18).
C.6 Energy correction at order  2 We now turn attention to the correction Φ [2] in (C.19b).Inserting the expression of  [1] in (C.12a), we get after rearranging the terms where W  () is the tensor introduced in the analysis of the leading order, see (B.3c), and As indicated by the 'delayed' keyword in Table C.1, we do not yet enforce the natural symmetries of " B (0) , which reflect the of the 1 2 (∇ ⊗ ∇) ⊗ ( ⊗ ) with which it gets contracted: they will be enforced later on the children of " B (0) .In the code, we implemented an extension of (C.25) that covers the rank-deficient case as well, see Equation (E.24) in Appendix E.
Note that Φ [2] no longer depends on  [2] thanks to the work done in Appendix C.3.It still depends on  [1] ( ) and its gradient, however.We proceed to remove the dependence on the gradient ∇ [1] ( ) by integrating by parts.

C.7 Integration by parts
The C () terms appearing in (C.24) can be integrated by parts as where Audoly and C. Lestringant An energy approach to higher-order homogenization where With the help of the 'content tracking' done in Table C.1, we obtain the unpacked form of these tensors as where the sub-tensors  1 (),  ′ 0 (),  2 (),  1 (),  0 (),  1 () and  0 () are identified as The expression of

D Special case of homogeneous properties
The special case of homogeneous properties (applicable to a perfectly periodic elastic truss for instance) is considered here.In this special case, the parameter ( ) goes away (  = 0).
The analysis of leading order is unchanged: it delivers  0 and  0 , which are no longer functions of , as illustrated in Table 5 for the non-deficient case.There are many simplifications at the next orders and the corresponding, specialized formulas are provided in Table D.2.Note that the tensor dimensions are changed compared to the general case, see also the 'usage' column in the table.As a consequence, the indices used in the transpose operations are affected.Although it is straightforward in principle, the specialization of the general formulas to this special case is cumbersome-to a point that we found it easier to re-derive them from scratch.Their consistency with the general formulas is checked in a dedicated Mathematica notebook2 .

E Extension to a rank-deficient matrix E.1 Special form of null vectors
Assuming that they exist, let us first characterize the null vectors of the symmetric matrix  () introduced in (B.5), entering in both the leading order problem (B.4) and in the determination of the corrective displacement (C.33).
For any  = (  ,   ) ∈ R   + c such that  () •  = 0, we have Multiplying the first equation by   , and using the second equation, we get   • W  () •   = 0. We observe that the assumption (12) (positive-definiteness of the energy on the subspace of admissible microscopic degrees of freedom) can be rewritten as: Therefore, we have   = 0, which then yields (Q•  ())  •  = 0, i.e., the   block is a null vector of (Q •   ())  .We have just shown The only way that the matrix  () can be singular is because of the Q •   () block.With  d denoting the rank deficiency of the matrix  or (Q •   ())  (both are the same by the argument above), we denote as  () ∈ T ( d , c ) a list of null vectors of (Q •   ())  , arranged in rows.Equation (E.2) then shows that the null vectors of  () are the rows of

E.2 Solutions of the linear equation
We consider the linear equation for a vector  ∈ R   + c ,  usage for an arbitrary choice of the coefficients l ∈ R  d .In the right-hand side, the first term is a particular solution furnished by the pseudo-inverse  † (), and the second term is a linear combination of the column-vectors in    () forming a basis of ker  (), with arbitrary coefficients ( l ) 1⩽ ⩽ d .

E.3 Extended macroscopic strain vector
When the matrix  is rank deficient, we append the  d coefficients l appearing in (E.6) to the macroscopic strain vector , and write where ľ are the usual macroscopic strain vector that defines the microscopic strain , see (2), referred to as  in the main body of the paper, while l are the additional parameters parametrizing the solution  of the rank-deficient linear problem.The dimension of  is now   = ň +  d .We denote the injection matrices Ǐ and Î of ľ and l into , respectively, which enable us to rewrite  as  = Ǐ • ľ + Î • l.Since Ǐ • Ǐ and Î • Î are orthogonal projections from the space R   in which  lives onto the subspaces with equations l = 0 and ľ = 0, respectively, the following identity holds, To capture the fact that the macroscopic strain  in (2) is a function of the original set of macroscopic degrees of freedom ľ, but not of the added l part, we require that any sub-block in   (),  ′  () or  ′′  () corresponding to a range of indices l vanishes, i.e., Indeed, Equation (E.10) warrants that the strain in (2) can be rewritten in terms of ľ and its gradients as where . (E.12) The proof of (E.11) is left to the reader.
In view of Equation (E.10) and Table 1, l does not appear anywhere in the specification of the problem: it is a set of free parameters that are reserved for parameterizing the solution (E.6) of the rank-deficient linear problem.An energy approach to higher-order homogenization

E.4 Changes to leading-order analysis
The solvability condition (E.5) yields •  () = 0. (E.14) The vector in square brackets is an output of the homogenization procedure, representing  d conditions that are linear in the macroscopic strain .
When (E.14) is satisfied, the solution is given by (E.6) as • Î = 0, (E.17) i.e., the operators W  and   do not sense the added degrees of freedom l.Combining with (E.9), this shows that Inserting into the expression of R(), we obtain We use this expression of R() in the code and not that proposed earlier in (B.7).Indeed, the latter can be recovered as a particular case: when the matrix  is invertible,  d = 0, implying that  † () =  −1 () and that   and Î are zero-dimension array, and the last term in (E.19) vanishes.
The definitions (B.8-B.15) of the other leading-order quantities such as  0 ,  0 , etc. are unchanged.
The following identity can be established using (B.8), (E.19), (E.3) and the orthogonality of the The stress  0 • , however, can depend on the l-block of  as well: by adapting the calculation in (E.20), one can show that  0 • Î =   () is non-zero in the rank-deficient case.The components of l can therefore be interpreted as the stress associated with the macroscopic kinematic constraint (E.14); this stress is akin to a Lagrange multiplier, i.e., is not set by any constitutive law.
Combining (E.20) with (B.11) and (B.12), one can show that the strain localization tensor  0 and the equivalent stiffness  0 are also uncoupled to l, implying zero l-sub-blocks, (E.22) The reason we prefer the expression of R in (E.19) to that derived first in (E.16) is that it makes it much more evident  0 ,  0 and  0 are insensitive to the added degrees of freedom l, see the identities (E.20-E.21).
B. Audoly and C. Lestringant An energy approach to higher-order homogenization

E.5 Changes to the energy expansion
With the help of (E.9), (E.22), (C.13), (C.12b) and (B.15a), the second term in the integrand in (C.17b) can be written as where the terms that have been crossed out are zero.In view of this, the definition of the operators " B (0) , " B (1) , C (0) , C (1) can be modified by including the projector Ǐ • Ǐ onto the subspace l = 0, to the right of  0 (), (E.24) The original and amended definitions in (C.25) and (E.24) are equally valid, but the latter has the advantage that it yields final tensors  1 and  0 having zero l-sub-blocks (the sub-blocks obtained with the former set of definitions do evaluate to zero when the solvability constraints are considered but this is much less evident, and potentially confusing).

E.6 Changes to the corrective displacement
The linear problem (C.33) makes use of the same matrix  as the leading-order problem.The solvability condition for  [1] is furnished by (E.5) as The quantity in square brackets is an output of the homogenization procedure that encodes  d conditions depending linearly on  and ∇.
As we did earlier at the leading order, the solution ( ★ having both a leading order contribution ( l) and an order  contribution ( l).A simple way to deal with this complication is (i) to discard the    () • l contribution in (E.26), and (ii) agree that l is a series in .This avoids extending the vector  with  d new entries for every homogenization order.
Concretely, we simply replace the inverse appearing in (C.36) by the Moore-Penrose inverse: our implementation uses . (E.27)

E.7 Solvability condition for 𝒚 [2]
The quantity  [2] entering in  [2] = (J 11 • ) :: (E.28) has been eliminated from (E.23) using the constraint Q •  [2] = 0.For  [2] to exist, one must have . In the code, a basis of vectors perpendicular to Im(Q •    ) is produced using a row-reduction algorithm, and we print out the conditions that each of these vectors is perpendicular to the vector Q • [. ..] appearing in the left-hand side above: this yields conditions depending linearly on ∇ ⊗ ∇ ⊗  and ∇ 2  ⊗ .

E.8 Summary: extension to rank-deficient problems
The following extension of the code enables us to deal with a rank-deficient matrix  (): • provide integers ň and  d and the injection matrix Î as a optional arguments to the homogenization procedure and check the condition (E.10) on the tensors   ,  ′  and  ′′  passed in argument; • compute a set of null vectors of the symmetric matrix  (), check that there are  d such vectors and that they are of the form (E.2), compute the Moore-Penrose inverse  † () if  d > 0; • return the solvability conditions (E.14), (E.25) and (E.30) whenever  d > 0; • replace Equations (B.7), (C.25) and (C.36) yielding R, " B (0) , " B (1) , C (0) , C (1) and R ′ with their extensions (E.19), (E.24) and (E.27) E.9 Illustration: a truss lattice with inextensible beams We consider again the truss lattice shown in Figure 2, but assume this time that the bars on the upper side (+) are inextensible: the corresponding stretching strain is constrained to be zero,  4 = 0.In view of ( 33) and ( 29), the inextensibility condition for the upper beams can be rewritten as () = 0 + O (). (E.31) (9b) where D  Φ[, , ; ] = lim →0 (Φ[,  + ] − Φ c [, ]) / denotes the directional derivative and ( ) is an arbitrary perturbation.Equation (9a) warrants that the stationary point  =  ★ [, ] satisfies the kinematic constraint, while Equation (9b) warrants that it is indeed a stationary point among the 's satisfying the kinematic constraints-the integral term takes care of these constraints by multiplying the Lagrange multipliers by the incremental form of the constraint, as usual in the calculus of variations.Having slaved the microscopic degrees of freedom  =  ★ [, ] to the macroscopic variables, we can define a homogenized energy functional Φ ★ [, ] by inserting the stationary point  ★ [, ] into the original Φ, Φ ★ [, ] = Φ[, ,  ★ [, ]].

Figure 2
Figure 2 Truss lattice example.(a) General view, (b) a specific bar.
) •  =  .(E.4) Multiplying by any null vector  of  and using the symmetry   =  , we obtain  •  = 0. Repeating this argument with all the null vectors that have been arranged into   (), we obtain  d solvability conditions   () •  = 0. (E.5)When (E.5) is satisfied, the solutions  of (E.4) can be expressed with the help of the Moore-Penrose inverse  † () of  () as  =  † () •  +    () • l. (E.6) The leading-order problem (B.4) is of the form (E.4) with = − W  (( )) Q •   (( ))•  ( ).(E.13)B.Audoly and C. Lestringant definition of the extended macroscopic strain vector  in (E.7), this solutions matches the form (B.6) used in the non-deficient case, provided we replace the inverse by the Moore-Penrose inverse and include a new term in the definition of R in (B.7),R() = − † () • W  () Q •   () +    () • Î .(E.16)It is convenient to rewrite this equation in a slightly different form, for a reason that will be discussed later.Using (E.10), one can show that W  ()Q •   ()

Table 1 .
next.Indeed, the stationarity condition of Ψ[, ] + Φ[, , ] with respect to  is nothing but that considered in (11), given that Ψ[, ] does not depend on .Next, it can be checked that the stationarity condition with respect to  of Ψ[, ] + Φ[, , ] is equivalent to the stationarity condition of the modified functional Ψ[, ] + Φ ★ [, ] based on the homogenized strain energy Φ ★ introduced in (10).This argument not only provides a justification to Equations (10) and (11), it also explains why the homogenized energy Φ ★ [, ] can be used as a substitute for the original energy Φ[, , ] in the analysis of the structural problem.

Table 4
(Audoly 2023ters used for the homogenization of the truss lattice.The spring constant  () is set up as a symbolic function of .For the sake of brevity, the terms denoted by ellipses are omitted: full expressions are available in the input file inhomogeneous-truss.nb included in the library(Audoly 2023).

Table D .
2 Higher-order homogenization in the special case of homogeneous properties.A complete implementation of the method in this special case is possible based on Table5(leading order, ignoring any dependence on ) and on the definitions appearing in the first column of the table above.The quantities appearing in the grey rows are the main results (localization tensors for corrective displacement  ′ 0 , Lagrange multipliers  ′ 0 , bulk energy contributions  0 and  0 and boundary contribution  0 ).