Homogenization and continuum limit of mechanical metamaterials

Mechanics of Materials 196 (2024) 105073 Available online 18 June 2024 0167-6636/© 2024 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/bync/4.0/). Contents lists available at ScienceDirect Mechanics of Materials journal homepage: www.elsevier.com/locate/mecmat Research paper Homogenization and continuum limit of mechanical metamaterials M.P. Ariza a,∗, S. Conti b, M. Ortizb,c aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los descubrimientos, s.n., 41092, Sevilla, Spain bInstitut für Angewandte Mathematik and Hausdorff Center for Mathematics, Universität Bonn, Endenicher Allee 60, 53115, Bonn, Germany cDivision of Engineering and Applied Science, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA ARTICLE INFO Dedicated to Alan Needleman, ab imo pectore, on the occasion of his 80th birthday. Keywords: Mechanical metamaterials Structural mechanics Homogenization Discrete-to-continuum Generalized elasticity ABSTRACT When used in bulk applications, mechanical metamaterials set forth a multiscale problem with many orders of magnitude in scale separation between the micro and macro scales. However, mechanical metamaterials fall outside conventional homogenization theory on account of the flexural, or bending, response of their members, including torsion. We show that homogenization theory, based on calculus of variations and notions of Gamma-convergence, can be extended to account for bending. The resulting homogenized metamaterials exhibit intrinsic generalized elasticity in the continuum limit. We illustrate these properties in specific examples including two-dimensional honeycomb and three-dimensional octet-truss metamaterials. 1. Introduction Mechanical metamaterials (Montemayor et al.,2015;Greer and Deshpande,2019;Lu et al.,2022), or architected materials, are characterized by a reticular structure that affords great flexibility for tailoring and shaping their mechanical properties (cf., e.g., Fleck et al. (2010), for a review of stiffness, strength and damage tolerance; Benedetti et al. (2021), for a review of fatigue performance; and Wu et al. (2021), for a review of dynamic properties). Advances in fabrication techniques (Greer and Deshpande,2019;Xia et al.,2022;Jin and Espinosa, 2024;Hamzehei et al.,2024) have brought within the realm of practicality new metamaterial designs that populate previously inaccessible areas of material-property space (Ashby,2006), including lightweight structures with high strength and elastic modulus (Lee et al.,2012; Jang et al.,2013;Bauer et al.,2014;Meza et al.,2014;Rys et al., 2014;Zheng et al.,2014;Meza et al.,2015;Gu and Greer,2015;Bauer et al.,2015;Montemayor et al.,2015;Bauer et al.,2016). Common designs include honeycomb structures (Bauer et al.,2014,2016) and the octet-truss (Deshpande et al.,2001b;Jang et al.,2013;Bauer et al.,2014;Zheng et al.,2014;Meza et al.,2015;Gu and Greer, 2015;Bauer et al.,2015), among others (Lee et al.,2012;Bauer et al., 2014;Rys et al.,2014;Rosário et al.,2015;Bauer et al.,2016;do Rosário et al.,2017). Metamaterials often exhibit scaling properties and size effects (Deshpande et al.,2001a,b) that exploit the ‘‘smaller is stronger’’ principle to beneficial effect, especially as regards fracture resistance (Chen et al.,1998;Shaikeea et al.,2022;Maurizi et al.,2022; Luan et al.,2022) and energy absorption (Zhang et al.,2023). ∗Corresponding author. E-mail addresses: [email protected] (M.P. Ariza), [email protected] (S. Conti), [email protected] (M. Ortiz). Mechanical nanostructured metamaterials, when used in bulk applications, define a quintessential multiscale problem, with many orders of magnitude in scale separation between the micro and macro domains. Under such conditions, models — computational or otherwise — predicated on the explicit tracking of every individual member of the metamaterial are impractical and uncalled for, since locally, away from defects and stress concentrations, large numbers of members deform collectively under applied loading and jointly exhibit well-defined effective behavior that can be characterized by means of a continuum material law. Coarse-graining approaches (Phlipot and Kochmann,2019; Kochmann et al.,2019) based on the quasicontinuum method (Tadmor et al.,1996a,b;Español et al.,2013), data-driven approaches (Korzeniowski and Weinberg,2022;Weinberg,2023;Jiao et al.,2023) and machine-learning approaches (Karapiperis and Kochmann,2023) have been proposed in order to bridge length scales while providing fullresolution of the lattice where necessary (see Fig. 1). In this work, we seek to characterize analytically the effective behavior of metamaterials at the macroscale under conditions of strict separation between the structural and lattice scales. The identification of such effective material behavior is the aim of homogenization theory (Cioranescu and Donato,1999) and discrete-to-continuum techniques (Cicalese et al., 2009;Braides and Gelli,2006). However, mechanical metamaterials fall outside conventional homogenization theory on account of the flexural, or bending, response of their members. We show that homogenization theory can be extended to account for bending and that the resulting homogenized metamaterials exhibit generalized elasticity of the micropolar type (Eringen,1966). https://doi.org/10.1016/j.mechmat.2024.105073 Received 14 March 2024; Received in revised form 12 June 2024; Accepted 14 June 2024 Mechanics of Materials 196 (2024) 105073 2 M.P. Ariza et al. Fig. 1. Octet-truss metamaterial (Deshpande et al.,2001a,b). Such theories are often used to introduce an ad hoc micromechanical length into the material law (cf., e.g., Askes and Aifantis (2011) for a review and discussion). By contrast, in the present context the micropolar terms are physical and are introduced by bending. In addition, the form of the limiting continuum energy is identified uniquely and unequivocally by the homogenization analysis. We additionally show that, for linear-elastic matematerials, the effective moduli of the limiting continuum energy can be computed explicitly using the discrete Fourier transform. We illustrate these properties in specific examples including two-dimensional honeycomb and three-dimensional octet metamaterials. 2. Mechanical metamaterials We consider throughout open-cell mechanical metamaterials comprised of bars, possessing both axial and bending stiffness, connected at rigid joints capable of transmitting forces and moments. 2.1. Notational conventions We denote by Zthe set of integers; by Rthe set of real numbers; and by Cthe set of complex numbers. The characteristic function of a subset 𝐴⊂𝑋is defined as 𝜒𝐴(𝑥) = {1,if 𝑥∈𝐴, 0,if 𝑥∈𝑋∖𝐴. (1) For a Lebesgue-measurable set 𝐴 ⊂ R𝑛, we denote by |𝐴|its Lebesgue measure, i.e., its length if 𝑛= 1; its area if 𝑛= 2; and its volume if 𝑛= 3. We denote by 𝑥⋅𝑦∈Rthe dot product of two vectors 𝑥,𝑦∈R𝑛; by 𝑦 ⊗ 𝑥 ∈R𝑛×𝑚the tensor product of two vectors 𝑥∈R𝑚and 𝑦∈R𝑛, i.e., (𝑦⊗𝑥)𝑗𝑖 =𝑦𝑗𝑥𝑖; and by 𝑦=𝐴𝑥 ∈R𝑛the product between a matrix 𝐴∈R𝑛×𝑚and a vector 𝑥∈R𝑚. We use throughout the abbreviation 𝐴𝑥 ⋅𝑦∶= (𝐴𝑥)⋅𝑦= 𝑛 ∑ 𝑗=1 𝑚 ∑ 𝑖=1 𝐴𝑗𝑖𝑦𝑗𝑥𝑖,(2) corresponding to assigning precedence to matrix multiplication over dot product. For a matrix 𝐴∈R𝑛×𝑛we denote by 𝐴𝑇its transpose and by sym 𝐴=1 2(𝐴+𝐴𝑇),skw 𝐴=1 2(𝐴−𝐴𝑇),(3) its symmetric and skewsymmetric parts. Given two vectors 𝑥∈R𝑚and 𝑦∈R𝑛, we denote by 𝑧= (𝑥;𝑦) ∈ R𝑚×R𝑛the cartesian product of 𝑥and 𝑦, i.e., 𝑧𝑖=𝑥𝑖, for 𝑖= 1,…, 𝑚 and 𝑧𝑗+𝑚=𝑦𝑗, for 𝑗= 1,…, 𝑛. Given two points 𝑥,𝑦∈R𝑛, we denote by [𝑥, 𝑦]the segment of straight line joining the two points in R𝑛. Given a function 𝑓∶R𝑛→R𝑚, we denote by 𝜕𝑖𝑓(𝑥)its partial derivative with respect to 𝑥𝑖. For functions of one variable, we also write 𝑓′(𝑥) = 𝜕𝑓(𝑥). We denote by 𝐷𝑓 ∶R𝑛→R𝑚×𝑛the matrix of partial derivatives of 𝑓, i.e., (𝐷𝑓)𝑖𝑗 =𝜕𝑗𝑓𝑖. Given a Lebesgue-measurable subset 𝐴 ⊂ R𝑛and an integrable function 𝑓∶𝐴→R, we denote by 𝐼=∫𝐴 𝑓(𝑥)𝑑𝑥 (4) the Lebesgue integral of 𝑓(𝑥)over 𝐴. 2.2. Geometry of metamaterials Mechanical metamaterials are locally periodic frame structures consisting of bars connected at a set of joints (cf., e.g., Fleck et al. (2010), Montemayor et al. (2015)). Infinite metamaterials are comprised of 𝑀types, or classes, of oriented bars, all the bars in a class being identical, including orientation, modulo translations. For instance, the bar classes of a honeycomb metamaterial are {⟍,—,⟋}, or 𝑀= 3, Fig. 10. The local environment of a joint is the set of bars incident on the joint. Metamaterials contain 𝑁different types, or classes, of local environments, all local environments in a class being identical modulo translations. Correspondingly, the joints of the metamaterial are classified into 𝑁types, according to their local environment. For instance, the joint classes of a honeycomb metamaterial are { Y , Y }, or 𝑁= 2,Fig. 10. The joints and bars of an infinite metamaterial are embedded periodically in space so that joints and bars of the same class span shifted Bravais lattices. Thus, the positions of the joints in each class span a point set of the form 𝑥(𝑙, 𝛼) = 𝑏𝛼+ 𝑛 ∑ 𝑖=1 𝑙𝑖𝑎𝑖, 𝑙 ∈Z𝑛, 𝛼 = 1,…, 𝑁, (5) where 𝑛is the dimension of space, e.g., 𝑛= 2 for two-dimensional lattices and 𝑛= 3 for three-dimensional lattices; (𝑎𝑖)𝑛 𝑖=1 is a basis of R𝑛;𝑙∶= (𝑙𝑖)𝑛 𝑖=1 is a corresponding array of integer lattice coordinates; and (𝑏𝛼)𝑁 𝛼=1 are translation vectors, or shifts. We additionally denote by 𝑉the volume of the unit cell of the spanning Bravais lattice of the metamaterial. The joints of a metamaterial can therefore be indexed by a pair (𝑙, 𝛼), where 𝛼∈ {1,…, 𝑁}designates the joint class and 𝑙∈Z𝑛 its lattice coordinate array. A full description of an infinite metamaterial additionally requires a scheme for indexing its bars by type and position, and for tabulating their connections to joints. By periodicity and translation invariance, the bars of a given bar class span the same Bravais lattice as the joint classes, modulo translations. Thus, bars can be indexed by pairs (𝑚, 𝛽), where 𝛽∈ {1,…, 𝑀}designates the bar class and 𝑚∈Z𝑛is a lattice coordinate array. In addition, an orientation of the bar classes identifies beginning and end joints in every bar, which we designate by the symbols ±, respectively (cf. Section 4for examples). The adjacency relations between bars and joints, or connectivity, can be defined by designating, for every bar (𝑚, 𝛽), the corresponding beginning and end joints, (𝑚−, 𝛽−)and (𝑚+, 𝛽+), respectively. We note that, by periodicity, ((𝑚+𝑡)±, 𝛽±)=(𝑚±+𝑡, 𝛽±),for all 𝑡∈Z𝑛,(6) i.e., the adjacency relations are translation-invariant. The corresponding joint coordinates follow from (5) as 𝑥±(𝑚, 𝛽) ∶= 𝑥(𝑚±, 𝛽±) = 𝑏𝛽±+ 𝑛 ∑ 𝑖=1 (𝑚±)𝑖𝑎𝑖,(7) Mechanics of Materials 196 (2024) 105073 3 M.P. Ariza et al. Fig. 2. Straight two-dimensional beam of constant cross section. Schematic of reference configuration and local deformation degrees of freedom. and the spanning vectors of the bars follow as 𝑑𝑥𝛽∶= 𝑥+(𝑚, 𝛽) − 𝑥−(𝑚, 𝛽),(8) which are independent of 𝑚by the translation invariance property (6). We also denote by ℎ(𝑚, 𝛽)=[𝑥−(𝑚, 𝛽), 𝑥+(𝑚, 𝛽)] (9) the spanning segments of the bars, regarded as sets. We refer to subsets of the infinite metamaterials just described as metastructures. The extent of a metastructure, i.e., its set of bars, may be characterized by means of an index set 𝐽={(𝑚, 𝛽) ∈ Z𝑛× {1,…, 𝑀} ∶ bar (𝑚, 𝛽)in metastructure}.(10) For an infinite metamaterial, 𝐽=Z𝑛× {1,…, 𝑀}. 2.3. Elastic energy of metastructures Within the framework of linearized kinematics, the joints of a metastructure are endowed with deflection and rotation-angle degrees of freedom, 𝑣(𝑙, 𝛼) ∈ R𝑛and 𝜃(𝑙, 𝛼) ∈ R𝑛(𝑛−1)∕2, respectively. We further denote by 𝑢(𝑙, 𝛼) ∶= (𝑣(𝑙, 𝛼); 𝜃(𝑙, 𝛼)) (11) the set of degrees of freedom of joint (𝑙, 𝛼). At the local level, we denote by 𝑈(𝑚, 𝛽) ∶= (𝑢−(𝑚, 𝛽); 𝑢+(𝑚, 𝛽)) (12) the set of local degrees of freedom of bar (𝑚, 𝛽), where 𝑢±(𝑚, 𝛽)are the degree-of-freedom arrays of the joints of the bar. The total energy of a linear-elastic metastructure of extent 𝐽has the form 𝐸(𝑢) = ∑ (𝑚,𝛽)∈𝐽 1 2S𝛽𝑈(𝑚, 𝛽)⋅𝑈(𝑚, 𝛽),(13) where S𝛽is the elastic stiffness of a bar of class 𝛽. Throughout this work, we consider straight bars of constant cross section. Other cases, such as bars in the form of curved beams, tapered beams, and others, can be treated likewise. 2.4. Straight two-dimensional beam of constant cross section Two-dimensional metamaterials are characterized by planar geometries with joint coordinates 𝑥(𝑙, 𝛼) ∈ R2,𝑛= 2. In addition, we assume that the metamaterial remains planar after deformation, so that the joint degrees of freedom 𝑢(𝑙, 𝛼) ∶= (𝑣(𝑙, 𝛼); 𝜃(𝑙, 𝛼)) of the joints comprise in-plane deflections 𝑣(𝑙, 𝛼) ∈ R2and a single rotation angle 𝜃(𝑙, 𝛼) ∈ R. In order to represent the energy of the metamaterial, we consider a straight beam of length 𝐿and constant cross section deforming in the (𝑥1, 𝑥2)-plane, in a reference configuration spanning the domain [0, 𝐿] on the 𝑥1axis, Fig. 2. We denote by 𝑢−∶= (𝑣−;𝜃−)and 𝑢+∶= (𝑣+;𝜃+) the degrees of freedom of the beginning and end joint of the beam respectively. The corresponding linearized rotation matrices are 𝑤±=(0 −𝜃± 𝜃±0)=∶∗ 𝜃±, 𝜃±=∗ 𝑤±,(14) where ∗denotes the Hodge-* operator. We additionally denote by 𝑈= (𝑢−;𝑢+)the array of local degrees of freedom. In this representation, an application of elementary beam theory (Timoshenko and Young,1965) gives the energy of the beam as 𝐸ref (𝑈) = 𝐸𝐴 2𝐿(𝑣+ 1−𝑣− 1)2+𝐸𝐼 2𝐿(𝜃+−𝜃−)2 +6𝐸𝐼 𝐿(𝑣+ 2−𝑣− 2 𝐿−𝜃−+𝜃+ 2)2, (15) where 𝐸is the Young’s modulus of the beam, 𝐴the cross-sectional area, and 𝐼the moment of inertia of the cross section with principal axis of inertia perpendicular to the plane. The first term in (15) is the energy due to the axial deformation of the beam, whereas the second and third terms account for the energy due to bending. The quadratic energy (15) can also be expressed in matrix form in terms of a 6 ×6 stiffness matrix Sref , whence the stiffness matrix S𝛽of the bars of class 𝛽in (13) follows as S𝛽=𝑇𝑇 𝛽Sref 𝑇𝛽,(16) in terms of some suitable, class-dependent, orthonormal change-of-axes matrix 𝑇𝛽. By direct inspection of (15), we verify that the global energy (13) is material-frame indifferent, i.e., invariant under transformations of the form 𝑣(𝑙, 𝛼)↦𝑣(𝑙, 𝛼) + 𝑤 𝑥(𝑙, 𝛼) + 𝑐, (17a) 𝜃(𝑙, 𝛼)↦𝜃(𝑙, 𝛼)+ ∗ 𝑤, (17b) where 𝑤∈𝑠𝑜(2) represents a rigid linearized rotation and 𝑐∈R2a rigid translation in the plane. From these properties, the global energy (13) follows directly in global coordinates as 𝐸(𝑢) = ∑ (𝑚,𝛽)∈𝐽{𝐸𝐴𝛽 2𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑1,𝛽)2 +𝐸𝐼𝛽 2𝐿𝛽 𝑑𝜃2(𝑚, 𝛽) +6𝐸𝐼𝛽 𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑2,𝛽 𝐿𝛽 − 𝜃(𝑚, 𝛽))2}, (18) where we write 𝑑𝑣(𝑚, 𝛽) = 𝑣(𝑚+, 𝛽+) − 𝑣(𝑚−, 𝛽−),(19a) 𝑑𝜃(𝑚, 𝛽) = 𝜃(𝑚+, 𝛽+) − 𝜃(𝑚−, 𝛽−),(19b)  𝜃(𝑚, 𝛽) = 𝜃(𝑚+, 𝛽+) + 𝜃(𝑚−, 𝛽−) 2,(19c) and 𝑑1,𝛽 ∶= 1 𝐿𝛽 𝑑𝑥𝛽, 𝑑2,𝛽 ∶= 𝑑⟂ 1,𝛽,(20) are orthonormal directors for each bar class, with 𝑑1,𝛽 a unit vector aligned with the oriented axis and 𝑑2,𝛽 a unit vector perpendicular to 𝑑1,𝛽 with orientation chosen such that (𝑑1,𝛽, 𝑑2,𝛽)defines a right-handed orthonormal basis together with the unit vector normal to the plane. The variables 𝜖(𝑚, 𝛽) = 1 𝐿𝛽 𝑑𝑣(𝑚, 𝛽)⋅𝑑1,𝛽,(21a) 𝜒(𝑚, 𝛽) = 1 𝐿𝛽 𝑑𝜃(𝑚, 𝛽),(21b) in (18) may be interpreted as local axial and bending strains, the first and second terms in (18) then supplying the corresponding axial and bending energies, respectively. Mechanics of Materials 196 (2024) 105073 4 M.P. Ariza et al. Fig. 3. Straight three-dimensional beam of constant cross section. (a) Schematic of reference configuration and local deformation degrees of freedom. (b) Cross section and directors. 2.5. Straight three-dimensional beam of constant cross section Three-dimensional metamaterials are characterized by spatial geometries with joint coordinates 𝑥(𝑙, 𝛼) ∈ R3,𝑛= 3, and joint degrees of freedom 𝑢(𝑙, 𝛼) ∶= (𝑣(𝑙, 𝛼); 𝜃(𝑙, 𝛼)) comprising deflections 𝑣(𝑙, 𝛼) ∈ R3and rotation angles 𝜃(𝑙, 𝛼) ∈ R3. With this parametrization, the corresponding linearized rotation matrices are 𝑤±=⎛⎜⎜⎜⎝ 0 −𝜃± 3𝜃± 2 𝜃± 30 −𝜃± 1 −𝜃± 2𝜃± 10 ⎞⎟⎟⎟⎠ ∶=∗ 𝜃±,(22a) 𝜃±=∗ 𝑤±,(22b) where we use the Hodge-* operator. With 𝑈= (𝑢−;𝑢+)as the array of local degrees of freedom, elementary beam theory (Timoshenko and Young,1965) now gives the energy of the beam as 𝐸ref (𝑈) = 𝐸𝐴 2𝐿(𝑣+ 1−𝑣− 1)2+𝐺𝐼1 2𝐿(𝜃+ 1−𝜃− 1)2 +𝐸𝐼2 2𝐿(𝜃+ 2−𝜃− 2)2+𝐸𝐼3 2𝐿(𝜃+ 3−𝜃− 3)2 +6𝐸𝐼2 𝐿(𝑣+ 3−𝑣− 3 𝐿+𝜃− 2+𝜃+ 2 2)2 +6𝐸𝐼3 𝐿(𝑣+ 2−𝑣− 2 𝐿−𝜃− 3+𝜃+ 3 2)2, (23) where 𝐴its cross-sectional area, 𝐼1its torsional moment of inertia accounting for warping of the cross section, 𝐼2and 𝐼3its principal moments of inertia about axes 𝑥2and 𝑥3, respectively, 𝐸is the Young’s modulus and 𝐺is the shear modulus (see Fig. 3). As in the two-dimensional case, the quadratic energy (23) can be expressed in matrix form in terms of a 12 ×12 stiffness matrix Sref , and the stiffness matrices S𝛽of the bars of class 𝛽in cf. (13) then follow as in (16) in terms of orthonormal transformation matrices 𝑇𝛽. The resulting global energy (13) is again material-frame-indifferent, i.e., invariant under transformations of the form (17), with 𝑤∈𝑠𝑜(3) now representing a rigid linearized rotation and 𝑐∈R3a rigid translation in three-dimensional space. Alternatively, from these properties the global energy (13) follows directly in global coordinates as 𝐸(𝑢) = ∑ (𝑚,𝛽)∈𝐽{𝐸𝐴𝛽 2𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑1,𝛽)2 + 𝐺𝐼1,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑1,𝛽)2 + 𝐸𝐼2,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑2,𝛽)2 + 𝐸𝐼3,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑3,𝛽)2 + 6𝐸𝐼2,𝛽 𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑3,𝛽 𝐿𝛽 + 𝜃(𝑚, 𝛽)⋅𝑑2,𝛽 )2 + 6𝐸𝐼3,𝛽 𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑2,𝛽 𝐿𝛽 − 𝜃(𝑚, 𝛽)⋅𝑑3,𝛽 )2} (24) where 𝑑𝑣(𝑚, 𝛽),𝑑𝜃(𝑚, 𝛽)and  𝜃(𝑚, 𝛽)are as in (19) and (𝑑1,𝛽 , 𝑑2,𝛽 , 𝑑3,𝛽 )is an orthonormal director triad in which 𝑑1,𝛽 is the direction of the bar axis, Eq. (20), and 𝑑2,𝛽 and 𝑑3,𝛽 are the principal directions of inertia of the cross section. As in the two-dimensional case, the variables 𝜖(𝑚, 𝛽) = 1 𝐿𝛽 𝑑𝑣(𝑚, 𝛽)⋅𝑑1,𝛽,(25a) 𝜒(𝑚, 𝛽) = 1 𝐿𝛽 𝑑𝜃(𝑚, 𝛽),(25b) may be interpreted as local axial and bending strains, the axial and bending energies in (24) then being quadratic in the corresponding strains. 3. The continuum limit of metamaterials Next, we concern ourselves with the limit in which the metamaterial cell-size is much smaller than the domain size, or continuum limit. To this end, we appeal to tools from calculus of variations, especially methods of 𝛤-convergence (Dal Maso,1993), as they bear on discrete- to-continuum problems (Alicandro and Cicalese,2004;Braides and Gelli,2006). The fundamental property of variational methods is that they ensure convergence of the solutions of the finite cell-size problem to the solutions of the limiting continuum model. However, we note that the local energy of the bars in (18) and (24) includes terms that penalize the difference between the deflection-compatible rotations and the mean rotations of the bars, e.g., 𝑑𝑣(𝑚, 𝛽)⋅𝑑2,𝛽∕𝐿𝛽and  𝜃(𝑚, 𝛽)in (18) and similar terms in (24), respectively. This local energy structure renders the attendant homogenization analysis non-standard and outside the scope of conventional discrete-to-continuum results (cf. Alicandro and Cicalese (2004), Braides and Gelli (2006)). In order to characterize the continuum limit, we consider throughout metastructures whose extent is defined by a domain 𝛺 ⊂ R𝑛. Specifically, the metastructure maximally contained in 𝛺is defined by the index set 𝐽𝛺= {(𝑚, 𝛽) ∈ Z𝑛× {1,…, 𝑀} ∶ ℎ(𝑚, 𝛽)⊂ 𝛺}(26) with 𝑥±(𝑚, 𝛽)as in (7) and ℎ(𝑚, 𝛽)as in (9). Thus, (26) records all the bars that are fully contained in 𝛺. It also supplies a criterion for terminating the metastructure at the boundary. We additionally denote by 𝐼𝛺the index set of the joints in the metastructure, i.e., 𝐼𝛺={(𝑙, 𝛼) ∈ Z𝑛× {1,…, 𝑁} ∶ ∃(𝑚, 𝛽) ∈ 𝐽𝛺 s. t. (𝑙, 𝛼)=(𝑚−, 𝛽−)or (𝑙, 𝛼)=(𝑚+, 𝛽+)}.(27) Henceforth, we denote by 𝐸(𝑢;𝛺)the energy (24) of the metastructure maximally contained within 𝛺. For ease of handling, we decompose the energy as 𝐸(𝑢;𝛺) ∶= 𝐴(𝑢;𝛺) + 𝐵(𝑢;𝛺) + 𝐶(𝑢;𝛺),(28) Mechanics of Materials 196 (2024) 105073 5 M.P. Ariza et al. Fig. 4. Bar-wise Voronoi cells of a honeycomb metamaterial relative to its hexagonal cells. Fig. 5. Octet metamaterial. (a) Space-filling volume cells (Ariza and Ortiz,2005, Fig. 8). (b) Voronoi cell of a bar relative to the volume cells. in terms of an axial energy 𝐴(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺 𝐸𝐴𝛽 2𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑1,𝛽)2,(29) abending energy, 𝐵(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺{𝐺𝐼1,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑1,𝛽)2 +(30) 𝐸𝐼2,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑2,𝛽)2 +𝐸𝐼3,𝛽 2𝐿𝛽(𝑑𝜃(𝑚, 𝛽)⋅𝑑3,𝛽)2}, and a coupling energy 𝐶(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺{ 6𝐸𝐼2,𝛽 𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑3,𝛽 𝐿𝛽 + 𝜃(𝑚, 𝛽)⋅𝑑2,𝛽 )2 + 6𝐸𝐼3,𝛽 𝐿𝛽(𝑑𝑣(𝑚, 𝛽)⋅𝑑2,𝛽 𝐿𝛽 − 𝜃(𝑚, 𝛽)⋅𝑑3,𝛽 )2}. (31) 3.1. Integral representation of the energy In order to facilitate the passage to the continuum and make contact with homogenization theory (Cioranescu and Donato,1999), we rewrite the energy (28) in integral form. To this end, we introduce a number of auxiliary structures. Let 𝑃(𝑚, 𝛽)be the Voronoi cell of bar (𝑚, 𝛽)in an infinite metamaterial, i.e., 𝑃(𝑚, 𝛽) = {𝑥∈R𝑛∶ dist(𝑥, ℎ(𝑚, 𝛽)) ≤ dist(𝑥, ℎ(𝑚′, 𝛽′)),∀(𝑚′, 𝛽′) ∈ Z𝑛× {1,…, 𝑀}},(32) i.e., 𝑃(𝑚, 𝛽)is the locus of points that are closer to be bar (𝑚, 𝛽)than to any other bar in the metamaterial. The measure 𝑉𝛽∶= |𝑃(𝑚, 𝛽)|(33) does not depend on 𝑚by periodicity and translation invariance. Figs. 4 and 5depict examples of Voronoi cells for a two-dimensional honeycomb structure and a three-dimensional octet lattice, respectively, by way of illustration. We note that the bar domains 𝑃(𝑚, 𝛽)are polygons, for 𝑛= 2, and polyhedra, for 𝑛= 3, consisting of simplices incident on the spanning segment of the bar, cf. Figs. 4 and 5. We denote by the collection of all such simplices, which together constitute a triangulation of R𝑛. We note that the vertex set of includes — but its larger than — the set of joints of the metamaterial. Consider now a metastructure maximally contained in 𝛺. We denote by 𝛺the restriction of to 𝛺. Precisely, 𝛺is the union of the Voronoi cells 𝑃(𝑚, 𝛽)of all bars (𝑚, 𝛽) ∈ 𝐽𝛺. We denote by [𝛺]the domain spanned by 𝛺. By this construction, 𝛺is a triangulation of [𝛺]. However, [𝛺]≠𝛺in general, though the discrepancy is rendered small upon scaling, in a sense to be made precise. Let 𝑢(𝑙, 𝛼)be a displacement of the metastructure and 𝑢(𝑥)a piecewise linear function, not renamed, supported on the triangulation 𝛺, Mechanics of Materials 196 (2024) 105073 6 M.P. Ariza et al. such that 𝑢(𝑥(𝑙, 𝛼)) = 𝑢(𝑙, 𝛼),∀ (𝑙, 𝛼) ∈ 𝐼𝛺.(34) Then, 𝐷𝑢(𝑥)is a piecewise constant function also supported on the triangulation 𝛺. We remark that the function 𝑢(𝑥)is not uniquely defined by the discrete variables 𝑢(𝑙, 𝛼), as additional degrees of freedom arise at corners of the Voronoi cells which do not correspond to joints. The value of 𝑢at these degrees of freedom can be fixed, for example, by minimizing the mean square oscillation of 𝐷𝑢 inside each cell of the metastructure, thus rendering the definition unique. In the case of simplicial metamaterials, as for example triangular metamaterials in two dimensions, 𝐷𝑢 can be made constant on each cell of the metastructure. Using the identity 𝐿𝛽𝐷𝑣(𝑥)𝑑1,𝛽 =𝑑𝑣(𝑚, 𝛽), 𝑥 ∈ℎ(𝑚, 𝛽),(35) and the piecewise-constant property of 𝐷𝑣(𝑥), we obtain 𝐴(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺 1 𝑉𝛽∫𝑃(𝑚,𝛽) 𝐸𝐴𝛽 2𝐿𝛽(𝐿𝛽𝐷𝑣(𝑥)𝑑1,𝛽 ⋅𝑑1,𝛽 )2𝑑𝑥. (36) Introducing the piecewise-constant axial moduli A(𝑥) = ∑ 𝑚∈Z𝑛 𝑀 ∑ 𝛽=1 A𝛽𝜒𝑃(𝑚,𝛽)(𝑥), 𝑥 ∈R𝑛,(37) with 𝜒𝑃(𝑚,𝛽)(𝑥)as in (1) and coefficients A𝛽=𝐸𝐴𝛽𝐿𝛽 𝑉𝛽 (𝑑1,𝛽 ⊗ 𝑑1,𝛽)⊗(𝑑1,𝛽 ⊗ 𝑑1,𝛽),(38) the axial energy (36) further reduces to the integral form 𝐴(𝑢;𝛺) = ∫[𝛺] 1 2A(𝑥)𝐷𝑣(𝑥)⋅𝐷𝑣(𝑥)𝑑𝑥. (39) We note that the inhomogeneous moduli A(𝑥)thus defined are periodic with the periodicity of the metamaterial. Proceeding likewise, we arrive at the integral form of the bending energy 𝐵(𝑢;𝛺) = ∫[𝛺] 1 2B(𝑥)𝐷𝜃(𝑥)⋅𝐷𝜃(𝑥)𝑑𝑥 (40) where B(𝑥) = ∑ 𝑚∈Z𝑛 𝑀 ∑ 𝛽=1 B𝛽𝜒𝑃(𝑚,𝛽)(𝑥), 𝑥 ∈R𝑛(41) are piecewise-constant bending moduli with coefficients B𝛽=𝐺𝐼1,𝛽 𝐿𝛽 𝑉𝛽 (𝑑1,𝛽 ⊗ 𝑑1,𝛽)⊗(𝑑1,𝛽 ⊗ 𝑑1,𝛽) +𝐸𝐼2,𝛽 𝐿𝛽 𝑉𝛽 (𝑑2,𝛽 ⊗ 𝑑1,𝛽)⊗(𝑑2,𝛽 ⊗ 𝑑1,𝛽) +𝐸𝐼3,𝛽 𝐿𝛽 𝑉𝛽 (𝑑3,𝛽 ⊗ 𝑑1,𝛽)⊗(𝑑3,𝛽 ⊗ 𝑑1,𝛽). (42) Finally, we turn to the coupling energy. As before, let 𝑢(𝑙, 𝛼) = (𝑣(𝑙, 𝛼); 𝜃(𝑙, 𝛼))be a displacement of the metastructure. In contrast to the fully conforming interpolation used in the foregoing, we now let 𝑣(𝑥)be a piecewise linear function and 𝜃(𝑥)a piecewise constant function, not renamed, both supported on the triangulation 𝛺, and such that 𝑣(𝑥(𝑙, 𝛼)) = 𝑣(𝑙, 𝛼),∀ (𝑙, 𝛼) ∈ 𝐼𝛺,(43a) 𝜃(𝑥) =  𝜃(𝑚, 𝛽), 𝑥 ∈𝑃(𝑚, 𝛽),∀ (𝑚, 𝛽) ∈ 𝐽𝛺.(43b) As before, 𝐷𝑣(𝑥)is a piecewise constant function supported on 𝛺. Using again the identity (35) and the piecewise-constant property of 𝐷𝑣(𝑥), we obtain 𝐶(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺{(44) 1 𝑉𝛽∫𝑃(𝑚,𝛽) 6𝐸𝐼2,𝛽 𝐿𝛽(𝐷𝑣(𝑥)𝑑1,𝛽 ⋅𝑑3,𝛽 +𝜃(𝑥)⋅𝑑2,𝛽 )2𝑑𝑥+ 1 𝑉𝛽∫𝑃(𝑚,𝛽) 6𝐸𝐼3,𝛽 𝐿𝛽(𝐷𝑣(𝑥)𝑑1,𝛽 ⋅𝑑2,𝛽 −𝜃(𝑥)⋅𝑑3,𝛽 )2𝑑𝑥 }. Equivalently, we have 𝐶(𝑢;𝛺) = ∑ (𝑚,𝛽)∈𝐽𝛺{(45) 1 𝑉𝛽∫𝑃(𝑚,𝛽) 6𝐸𝐼2,𝛽 𝐿𝛽((𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))𝑑1,𝛽 ⋅𝑑3,𝛽 )2𝑑𝑥+ 1 𝑉𝛽∫𝑃(𝑚,𝛽) 6𝐸𝐼3,𝛽 𝐿𝛽((𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))𝑑1,𝛽 ⋅𝑑2,𝛽 )2𝑑𝑥 }, in terms of the Hodge-* operators (14) and (22). Introducing the piecewise-constant coupling moduli C(𝑥) = ∑ 𝑚∈Z𝑛 𝑀 ∑ 𝛽=1 C𝛽𝜒𝑃(𝑚,𝛽)(𝑥), 𝑥 ∈R𝑛,(46) with coefficients C𝛽=1 𝑉𝛽 12𝐸𝐼2,𝛽 𝐿𝛽 (𝑑1,𝛽 ⊗ 𝑑3,𝛽)⊗(𝑑1,𝛽 ⊗ 𝑑3,𝛽) +1 𝑉𝛽 12𝐸𝐼3,𝛽 𝐿𝛽 (𝑑1,𝛽 ⊗ 𝑑2,𝛽)⊗(𝑑1,𝛽 ⊗ 𝑑2,𝛽), (47) the coupling energy (45) reduces to the integral form 𝐶(𝑢;𝛺) = ∫[𝛺] 1 2C(𝑥)(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))⋅(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))𝑑𝑥, (48) which penalizes differences between the deflection-compatible rotations and the mean rotations of the bars. We again note that the inhomogeneous moduli C(𝑥)thus defined are periodic with the periodicity of the metamaterial. 3.2. Scaling Next, we consider a sequence of metastructures composed of the same metamaterial but spanning increasingly larger self-similar domains 𝛺∕𝜖, where 𝜖 > 0is a small scaling factor, Fig. 6. The corresponding metastructure is defined by the index set 𝐽𝛺∕𝜖∶= {(𝑚, 𝛽) ∈ Z𝑛× {1,…, 𝑀} ∶ ℎ(𝑚, 𝛽)⊂ 𝛺∕𝜖}.(49) Note that, in this scaling, the geometry of the metamaterial at the unitcell scale as well and the material properties and the geometry of the cross sections of the bars, remain unchanged along the sequence. In particular, the bar lengths 𝐿𝛽, cross-sectional areas 𝐴𝛽and moments of inertia 𝐼𝑖,𝛽 are not scaled according to 𝜖. By contrast, the number of bars in the metastructure grows as #𝐼𝛺∕𝜖∼𝑂(𝜖−𝑛)and the metastructure spans an infinite metamaterial in the infinite body limit 𝜖→0. With 𝜖→0, we consider the sequence of energy functions 𝐸𝜖(𝑢;𝛺) = 𝜖𝑛𝐸(𝑢𝜖;𝛺∕𝜖),(50) with 𝑥𝜖=𝜖−1 𝑥, 𝑣𝜖(𝑥𝜖) = 𝜖−1 𝑣(𝑥), 𝜃𝜖(𝑥𝜖) = 𝜃(𝑥).(51) The factor 𝜖𝑛in (50) accounts for the expected elasticity scaling of the limiting energy and is included to ensure a proper limit. From (39),(40) and (48), we have the identity 𝐸(𝑢;𝛺∕𝜖) = 𝐴(𝑢;𝛺∕𝜖) + 𝐵(𝑢;𝛺∕𝜖) + 𝐶(𝑢;𝛺∕𝜖) = Mechanics of Materials 196 (2024) 105073 7 M.P. Ariza et al. Fig. 6. Scaling of domain 𝛺with constant metamaterial. Fig. 7. Scaling of metamaterial with constant domain 𝛺. ∫[𝛺∕𝜖] 1 2A(𝑥)𝐷𝑣(𝑥)⋅𝐷𝑣(𝑥)𝑑𝑥 + ∫[𝛺∕𝜖] 1 2B(𝑥)𝐷𝜃(𝑥)⋅𝐷𝜃(𝑥)𝑑𝑥 + (52) ∫[𝛺∕𝜖] 1 2C(𝑥)(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))⋅(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))𝑑𝑥, for displacement fields 𝑢(𝑥)in the class of piecewise-linear or piecewiseconstant functions supported on defined in Section 3.1. The structure of (52) is noteworthy. Thus, the energy consists of three terms: an axial energy; a bending energy, including torsion; and a coupling energy that couples deflection-compatible rotations and free rotations. Energies of the form (52) were termed micropolar by Eringen (1966) and set forth an example of a theory of generalized continua. Such theories are often used to introduce an ad hoc micromechanical length into the material law that renders the energy non-local, in an effort to account for size effects. By contrast, in the present context the micropolar terms are physical and introduced into the energy by bending (see Fig. 7). Rescaling the domain back to 𝛺according to (50) and (51), we obtain 𝐸𝜖(𝑢;𝛺) = (53) ∫𝜖[𝛺∕𝜖] 1 2A(𝑥∕𝜖)𝐷𝑣(𝑥)⋅𝐷𝑣(𝑥)𝑑𝑥+ 𝜖2∫𝜖[𝛺∕𝜖] 1 2B(𝑥∕𝜖)𝐷𝜃(𝑥)⋅𝐷𝜃(𝑥)𝑑𝑥+ ∫𝜖[𝛺∕𝜖] 1 2C(𝑥∕𝜖)(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))⋅(𝐷𝑣(𝑥)− ∗ 𝜃(𝑥))𝑑𝑥, where we have dropped the prime for simplicity and the new displacement fields 𝑢(𝑥)are in the class of piecewise-linear or piecewiseconstant functions supported on the scaled triangulation 𝜖. 3.3. Continuum limit We wish to ascertain the asymptotic behavior of the forced energyminimizing configurations (𝑢𝜖)of the discrete energies (53) as 𝜖→0, or continuum limit. Specifically, we wish to determine whether a continuum limiting energy 𝐸0exists with the property that the sequence of minimum discrete potential energies converges to the minimum continuum potential energy for all admissible forcings. Then, we say that 𝐸0is the variational limit, or 𝛤-limit, of the sequence (𝐸𝜖)(Dal Maso, 1993). Thus, suppose that the metastructure is acted upon by macroscopic distributed actions defined by the functions 𝑓0(𝑥) ∶= (𝑞0(𝑥); 𝑚0(𝑥)),(54) where 𝑞0(𝑥)and 𝑚0(𝑥)are distributed forces and moments, or torques, per unit volume, respectively. The work of the applied actions is then ⟨𝑓0, 𝑢⟩∶= 𝑉 𝑁∑ 𝑙∈Z𝑛 𝑁 ∑ 𝛼=1 𝑓0(𝑥(𝑙, 𝛼)) ⋅𝑢(𝑙, 𝛼),(55) where 𝑉is the volume of the unit cell of the spanning Bravais lattice of the metamaterial and the factor 𝑉∕𝑁represents the volume per joint class. Its value follows from the property that all joint classes have the same number of joints per unit volume, or joint density. Consider now a sequence (𝐸𝜖)of scaled energies defined in (53). Assume now that the sequence of energies is uniformly stable, or equicoercive. Specifically, we assume that there is a constant 𝐶 > 0such that 𝐸𝜖(𝑢)≥𝐶‖𝑢‖2,(56) for some suitable energy norm ‖⋅‖. This equicoercivity condition expresses the property that 𝐸𝜖(𝑢)is bounded below for all 𝜖by a structurally stable quadratic energy, i.e., a positive quadratic energy devoid of zero-energy modes. This property effectively prevents the emergence of finely oscillating unstable modes that corrupt the solution, at no cost in energy, as 𝜖→0. The norm entering (56) needs to be chosen carefully, and a suitable choice may not be immediately apparent for some metastructures. Besides the usual fact that it should only control 𝑢modulo rigid-body motions, one needs to consider several conditions. First, the metamaterial needs to be appropriately truncated at the boundary, eliminating for example dangling bonds from the domain on which ‖⋅‖is computed. Second, the norm needs to be strong enough to ensure continuity of the energy up to boundary terms. Third, the ancillary degrees of freedom need to be eliminated as discussed after (34), and the resulting metamaterial needs to be sufficiently rigid to provide coercivity. For example, for a triangular metamaterial in two dimensions, the axial energy immediately provides control of the 𝐿2norm of the strain (𝐷𝑣+𝐷𝑣𝑇)∕2, Mechanics of Materials 196 (2024) 105073 8 M.P. Ariza et al. which leads to control of 𝑣and 𝐷𝑣 via Korn’s inequality and appropriate boundary data. Using this property, the coupling energy gives control of 𝜃, and finally the bending energy gives control of 𝐷𝜃∕𝜖. The situation for the honeycomb lattice is more complex as mechanisms exist for the axial energy, which correspond to zero-axial-energy deformations in addition to rigid-body motions. These zero-energy modes need to be stabilized by the bending and coupling energies in the long-wavelength limit relevant to the macroscopic material behavior. With these provisos, the equilibrium configuration 𝑢𝜖of the metastructure then follows by minimizing the potential energy 𝐹𝜖(𝑢) ∶= 𝐸𝜖(𝑢) − ⟨𝑓0, 𝑢⟩.(57) Denote by 𝑚(𝐸𝜖;𝑓0)the corresponding minimum potential energy. Then, the variational limit, or 𝛤-limit, 𝐸0of the sequence (𝐸𝜖)is characterized by the property that lim 𝜖→0𝑚(𝐸𝜖;𝑓0) = 𝑚(𝐸0;𝑓0),(58) for all admissible forcings 𝑓0(𝑥), where 𝑚(𝐸0;𝑓0)is the minimum of the limiting potential energy 𝐹0(𝑢0) ∶= 𝐸0(𝑢0) − ⟨𝑓0, 𝑢0⟩,(59) (cf., e.g., Dal Maso (1993, Theorem 13.5, Corollary 13.7)). Variational convergence, together with the equicoercivity condition (56), in turn ensures the weak convergence of the sequence of minimizers (𝑢𝜖)to the minimizer 𝑢0of 𝐸0for all admissible forcings 𝑓0. In view of the structure of the energy (28), the corresponding variational limit as 𝜖→0extends variational discrete-to-continuum problems (Braides and Gelli,2006;Alicandro and Cicalese,2004) to account for bending. Alternatively, appealing to the density of piecewise linear and constant displacements over the increasingly finer triangulations 𝜖, the energy (53) can be extended (by infinity) to continuum displacement fields without altering the continuum limit. The resulting limiting problem then falls within the class of variational homogenization problems (cf., e.g., Cioranescu and Donato (1999)). In addition, in view of the 𝜖2factor in (53) we expect the bending energy to be negligible in the limit, and the limiting continuum energy to be of the form 𝐸0(𝑢0;𝛺) = ∫𝛺 𝑊0(𝐷𝑣0(𝑥), 𝜃0(𝑥))𝑑𝑥, (60) for some effective energy density 𝑊0(𝐷𝑣0, 𝜃0), to be determined, and for a certain class of sufficiently regular and finite domains 𝛺. In particular, the effective energy density 𝑊0(𝐷𝑣0, 𝜃0)is expected to be independent of 𝛺within the admissible class. 3.4. Discrete fourier transform analysis Since, by assumption, the limiting energy density is independent of the domain 𝛺, we can conveniently identify 𝑊0(𝐷𝑣0, 𝜃0)by considering metamaterials of infinite extent. In order to make calculations explicit, it proves convenient to revert to the discrete form (24) or the energy, with 𝐽=Z𝑛× {1,…, 𝑀}. We can exploit the periodicity of the metamaterial by testing the energy with discrete harmonic displacements 𝑢(𝑙, 𝛼) = 1 (2𝜋)𝑛∫𝐵 𝑢(𝑘, 𝛼) e𝑖𝑘⋅𝑥(𝑙,𝛼)𝑑𝑘, (61) where 𝐵is the Brillouin zone of the metastructure, 𝑥(𝑙, 𝛼)are the joint coordinates (5) and 𝑢(𝑘, 𝛼)is the discrete Fourier transform of the discrete displacement field 𝑢(𝑙, 𝛼), cf. Appendix A. A straightforward calculation using Parseval’s theorem (A.5) then gives 𝐴(𝑢) = 𝑀 ∑ 𝛽=1 1 (2𝜋)𝑛∫𝐵 1 𝑉{𝐸𝐴𝛽 2𝐿𝛽|||𝑑 𝑣(𝑘, 𝛽)⋅𝑑1,𝛽 ||| 2}𝑑𝑘, (62) for the axial energy, 𝐵(𝑢) = 𝑀 ∑ 𝛽=1 1 (2𝜋)𝑛∫𝐵 1 𝑉{ 𝐺𝐼1,𝛽 2𝐿𝛽|||𝑑 𝜃(𝑘, 𝛽)⋅𝑑1,𝛽 ||| 2+ 𝐸𝐼2,𝛽 2𝐿𝛽|||𝑑 𝜃(𝑘, 𝛽)⋅𝑑2,𝛽 ||| 2+ 𝐸𝐼3,𝛽 2𝐿𝛽|||𝑑 𝜃(𝑘, 𝛽)⋅𝑑3,𝛽 ||| 2}𝑑𝑘, (63) for the bending energy, and 𝐶(𝑢) = 𝑀 ∑ 𝛽=1 1 (2𝜋)𝑛∫𝐵 1 𝑉{ 6𝐸𝐼2,𝛽 𝐿𝛽||| 𝑑 𝑣(𝑘, 𝛽)⋅𝑑3,𝛽 𝐿𝛽 +  𝜃(𝑘, 𝛽)⋅𝑑2,𝛽 ||| 2 + 6𝐸𝐼3,𝛽 𝐿𝛽||| 𝑑 𝜃(𝑘, 𝛽)⋅𝑑2,𝛽 𝐿𝛽 −  𝜃(𝑘, 𝛽)⋅𝑑3,𝛽 ||| 2}𝑑𝑘, (64) for the coupling energy, where 𝑉is the volume of the unit cell of the spanning Bravais lattice of the metamaterial and we write 𝑑 𝑣(𝑘, 𝛽)=e𝑖 2𝑘⋅𝑑𝑥𝛽𝑣(𝑘, 𝛽+)−e−𝑖 2𝑘⋅𝑑𝑥𝛽𝑣(𝑘, 𝛽−),(65a) 𝑑 𝜃(𝑘, 𝛽)=e𝑖 2𝑘⋅𝑑𝑥𝛽 𝜃(𝑘, 𝛽+)−e−𝑖 2𝑘⋅𝑑𝑥𝛽 𝜃(𝑘, 𝛽−),(65b)   𝜃(𝑘, 𝛽) = 1 2(e𝑖 2𝑘⋅𝑑𝑥𝛽 𝜃(𝑘, 𝛽+)+e−𝑖 2𝑘⋅𝑑𝑥𝛽 𝜃(𝑘, 𝛽−)).(65c) Collecting all joint displacements into an array, not renamed, 𝑢(𝑘) ∶= (𝑢(𝑘, 1),…, 𝑢(𝑘, 𝑁)),(66) the energy then takes the compact form 𝐸(𝑢) = 1 (2𝜋)𝑛∫𝐵 1 2D(𝑘)𝑢(𝑘)⋅𝑢∗(𝑘)𝑑𝑘, (67) in terms of the dynamical matrix D(𝑘)of the metamaterial. We note that D(𝑘)is Hermitian, positive definite and depends smoothly on the wavevector 𝑘. We additionally note that D(𝑘)follows explicitly from (62),(63),(64) and (65), by a simple rearrangement of terms, as a function of the beam properties: 𝐸𝐴𝛽,𝐺𝐼1,𝛽 ,𝐸𝐼2,𝛽 and 𝐸𝐼3,𝛽 ; and the metamaterial geometry: 𝐿𝛽,𝑑𝑥𝛽,𝑑1,𝛽 ,𝑑2,𝛽 and 𝑑3,𝛽 . In the Fourier representation, an application of Parseval’s theorem (A.5) to the work function (55) further gives ⟨𝑓0, 𝑢⟩∶= 1 (2𝜋)𝑛∫𝐵  𝑓0(𝑘)(1 𝑁 𝑁 ∑ 𝛼=1 𝑢∗(𝑘, 𝛼))𝑑𝑘, (68) where  𝑓0(𝑘)is the ordinary Fourier transform of 𝑓0(𝑥). The scaling (50) and (51) defines the sequence of scaled energies 𝐸𝜖(𝑢) = 1 (2𝜋)𝑛∫𝐵∕𝜖 1 2D𝜖(𝑘)𝑢(𝑘)⋅𝑢∗(𝑘)𝑑𝑘, (69) where the scaled dynamical matrix D𝜖(𝑘)is defined by the property D𝜖(𝑘)𝜁⋅𝜁∗=𝜖−2 D(𝜖𝑘) (𝜉;𝜖𝜂)⋅(𝜉∗;𝜖𝜂∗), ∀𝜁∶= (𝜉;𝜂)∈[C𝑛]𝑁× [C𝑛(𝑛−1)∕2]𝑁, 𝑘 ∈𝐵∕𝜖, (70) Likewise, as 𝜖→0the work function (68) can be expressed independently of 𝜖as ⟨𝑓0, 𝑢⟩∶= 1 (2𝜋)𝑛∫R𝑛  𝑓0(𝑘)(1 𝑁 𝑁 ∑ 𝛼=1 𝑢∗(𝑘, 𝛼))𝑑𝑘. (71) provided that 𝑢(𝑘)is extended to zero outside 𝐵∕𝜖and  𝑓0(𝑘)is kept fixed throughout the sequence. Suppose now that the metamaterial is loaded by macroscopic forces 𝑓0(𝑥). Then, from the Fourier representations (68) and (69) the corresponding equilibrium equations are D𝜖(𝑘)𝑢𝜖(𝑘) =  𝑓(𝑘), 𝑘 ∈𝐵∕𝜖, (72) Mechanics of Materials 196 (2024) 105073 9 M.P. Ariza et al. where again we collect all the joint forces into an array  𝑓(𝑘) ∶= (1 𝑁 𝑓0(𝑘),𝑁 …,1 𝑁 𝑓0(𝑘))∶= 𝐿 𝑓0(𝑘).(73) We note that the operator 𝐿thus defined localizes the continuum forces to each of the joint classes. In this representation, an appropriate form of the equicoercivity condition (56) is that D𝜖(𝑘)𝜁⋅𝜁∗≥𝐶(|𝑘|2|𝜉|2+|𝜂|2),(74) for all 𝜁∶= (𝜉;𝜂) ∈ C𝑛×C𝑛(𝑛−1)∕2 and some 𝐶 > 0. We note that this condition effectively implies a corresponding energy norm ‖⋅‖, Eq. (56). Under these conditions, the matrix D𝜖(𝑘)is non-singular and the equilibrium problem (72) can be solved pointwise for 𝑘∈𝐵∕𝜖, with the result 𝑢𝜖(𝑘) = D−1 𝜖(𝑘) 𝑓(𝑘), 𝑘 ∈𝐵∕𝜖, (75) which explicitly characterizes the sequence of energy-minimizing displacements (𝑢𝜖)over their corresponding Brillouin zones 𝐵∕𝜖. We note that we can equivalently extend the functions 𝑢𝜖(𝑘)to all of R𝑛by zero. The corresponding inverse Fourier transforms 𝑢𝜖(𝑥) are then said to be band-limited and are referred to as the Whittaker- Shannon interpolation of the discrete displacements 𝑢𝜖(𝑙, 𝛼)(Mallat, 2009;Español et al.,2013). The corresponding minimum potential energy is 𝑚(𝐸𝜖;𝑓0) ∶= − 1 (2𝜋)𝑛∫𝐵∕𝜖 1 2D−1 𝜖(𝑘) 𝑓(𝑘)⋅ 𝑓∗(𝑘)𝑑𝑘. (76) Passing to the limit as in (58) with the aid of (74) and dominated convergence, we obtain lim 𝜖→0𝑚(𝐸𝜖;𝑓0) ∶= 𝑚(𝐸0;𝑓0)=−𝐸∗ 0(𝑓0) = −1 (2𝜋)𝑛∫R𝑛 1 2D−1 0(𝑘) 𝑓0(𝑘)⋅ 𝑓∗ 0(𝑘)𝑑𝑘, (77) where D0(𝑘) = (lim 𝜖→0𝐿𝑇D−1 𝜖(𝑘)𝐿)−1, 𝑘 ∈R𝑛,(78) is the effective dynamical matrix in the continuum limit and 𝐸∗ 0(𝑓0) is the dual continuum energy. A final Legendre transform gives the effective continuum energy as 𝐸0(𝑢0) = 1 (2𝜋)𝑛∫R𝑛 1 2D0(𝑘)𝑢0(𝑘)⋅𝑢∗ 0(𝑘)𝑑𝑘, (79) for all continuum displacements 𝑢0(𝑥). The limiting equilibrium problem for the continuum displacements 𝑢0(𝑥)consists of minimizing the homogenized potential energy (59) for given forcing 𝑓0(𝑥). In the Fourier representation, the solution is 𝑢0(𝑘) = D−1 0(𝑘) 𝑓0(𝑘), 𝑘 ∈R𝑛,(80) which combined with (75) gives 𝑢𝜖(𝑘) = D−1 𝜖(𝑘)𝐿D0(𝑘)𝑢0(𝑘), 𝑘 ∈𝐵∕𝜖. (81) This relation shows how the continuum displacements 𝑢0(𝑥)are approximated by discrete displacements 𝑢𝜖(𝑙, 𝛼)along the sequence that are the result of energy relaxation at the unit-cell level. We additionally note that, for metamaterials with one class of joints only, 𝑁= 1, the matrix 𝐿reduces to the identity, whereupon (78) in turn reduces to Ariza and Ortiz (2005) D0(𝑘) = lim 𝜖→0 D𝜖(𝑘),(82) i.e., D0(𝑘)is the pointwise limit of D𝜖(𝑘)as 𝜖→0, or long wavelength limit. 3.5. Structure of the continuum limit We proceed to elucidate the structure of the limiting continuum energy (79) as a quadratic form. A lengthy but straightforward calculation, consigned to Appendix B for the sake of continuity, shows that D0(𝜆𝑘)𝜁0⋅𝜁∗ 0=D0(𝑘) (𝜆𝜉0;𝜂0)⋅(𝜆𝜉∗ 0;𝜂∗ 0).(83) for every 𝜆 > 0and 𝜁0∶= (𝜉0;𝜂0) ∈ C𝑛×C𝑛(𝑛−1)∕2,(84) i.e., rescaling 𝑘by 𝜆in the continuum dynamical matrix D0(𝑘)is equivalent to rescaling the deflection amplitude 𝜉0by 𝜆, while simultaneously keeping the rotation amplitude 𝜂0invariant. From this property, it follows that the energy density admits the representation 1 2D0(𝑘)𝜁0⋅𝜁∗ 0=𝑊0(𝑖𝑘 𝜉0, 𝜂0),(85) where 𝑊0is a quadratic form, and the energy (79) becomes 𝐸0(𝑢0) = 1 (2𝜋)𝑛∫R𝑛 𝑊0(𝑖𝑘 𝑣0(𝑘), 𝜃0(𝑘))𝑑𝑘. (86) An inverse Fourier transform then gives the energy as 𝐸0(𝑢0) = ∫R𝑛 𝑊0(𝐷𝑣0(𝑥), 𝜃0(𝑥))𝑑𝑥, (87) or (60) for a domain 𝛺, as surmised. Further restrictions on 𝑊0(𝐷𝑣0(𝑥), 𝜃0(𝑥))result from material-frame invariance, Eq. (17), which requires that 𝑊0(𝛽, 𝜃)=𝑊0(𝛽+𝑤, 𝜃+ ∗ 𝑤),(88) for all deformation gradients 𝛽∈R𝑛×𝑛, rotation angles 𝜃∈R𝑛(𝑛−1)∕2 and rotation matrices 𝑤∈𝑠𝑜(𝑛). Choosing 𝑤= −skw𝛽, Eq. (88) specializes to 𝑊0(𝛽, 𝜃)=𝑊0(𝜀, 𝜃− ∗ skw𝛽)(89) for some other function 𝑊0, not renamed, of the strains 𝜀= sym𝛽 and the difference between the free rotations 𝜃and the deflectioncompatible rotations ∗ skw𝛽. Inserting (89) into (87), we obtain the representation 𝐸0(𝑢0) = ∫𝛺 𝑊0(𝜀0(𝑥), 𝜃0(𝑥)− ∗ skw𝐷𝑣0(𝑥))𝑑𝑥, (90) where 𝜀0(𝑥) = sym𝐷𝑣0(𝑥)(91) are the local strains. We thus conclude that, to lowest order, the effective continuum energy (90) of linear metastructures is a special case of linear micropolar elasticity, in the sense of Eringen (1966), in which the energy density is independent of the curvature, or bending strain, 𝐷𝜃0(𝑥). Because of this special form, the energy (87) is local, i.e., it exhibits the local elasticity scaling 𝐸0(𝑢′ 0;𝛺′) = 𝜆𝑛𝐸0(𝑢0;𝛺),(92) under the transformation 𝑥′=𝜆𝑥, 𝑣′ 0=𝜆𝑣0, 𝜃′ 0=𝜃0.(93) In particular, to lowest-order the continuum energy does not account for size effects. We additionally note that the variational convergence properties of the energy carry over to potential energies involving forcing terms. We also expect the energy-minimizing configurations of the continuum problem to restrict properly to the boundary in the sense of Mechanics of Materials 196 (2024) 105073 16 M.P. Ariza et al. Tadmor, E.B., Phillips, R., Ortiz, M., 1996b. Mixed atomistic and continuum models of deformation in solids. Langmuir 12 (19), 4529–4534. Timoshenko, S.P., Young, D.H., 1965. Theory of Structures, second ed. McGraw-Hill Book Company, New York. Ulloa, J., Ariza, M.P., Andrade, J.E., Ortiz, M., 2024. Fracture and size effect in mechanical metamaterials. J. Mech. Phys. Solids (Submitted for publication). Weinberg, K., 2023. Data-driven finite element computation of microstructured materials. PAMM 23 (4), e202300285. Weiner, J.H., 2002. Statistical mechanics of elasticity, second ed. Dover Publications, Mineola, N.Y.. Wu, L., Wang, Y., Chuang, K., Wu, F., Wang, Q., Lin, W., Jiang, H., 2021. A brief review of dynamic mechanical metamaterials for mechanical energy manipulation. Mater. Today 44, 168–193. Xia, X., Spadaccini, C., Greer, J., 2022. Responsive materials architected in space and time. Nat. Rev. Mater. 7 (9), 683–701. Zhang, P., Yu, P., Zhang, R., Chen, X., Tan, H., 2023. Grid octet truss lattice materials for energy absorption. Int. J. Mech. Sci. 259, 108616. Zheng, X., Lee, H., Weisgraber, T., Shusteff, M., DeOtte, J., Duoss, E., Kuntz, J., Biener, M., Ge, Q., Jackson, J., Kucheyev, S., Fang, N., Spadaccini, C., 2014. Ultralight, ultrastiff mechanical metamaterials. Science 344, 1373–1377.