## Abstract

Mechanical metamaterials offer exotic properties based on local control of cell geometry and their global configuration into structures and mechanisms. Historically, these have been made as continuous, monolithic structures with additive manufacturing, which affords high resolution and throughput, but is inherently limited by process and machine constraints. To address this issue, we present a construction system for mechanical metamaterials based on discrete assembly of a finite set of parts, which can be spatially composed for a range of properties such as rigidity, compliance, chirality, and auxetic behavior. This system achieves desired continuum properties through design of the parts such that global behavior is governed by local mechanisms. We describe the design methodology, production process, numerical modeling, and experimental characterization of metamaterial behaviors. This approach benefits from incremental assembly, which eliminates scale limitations, best-practice manufacturing for reliable, low-cost part production, and interchangeability through a consistent assembly process across part types.

## INTRODUCTION

The notion of rationally designing a material from the microscale to the macroscale has been a long-standing goal with broad engineering applications. By controlling local cell properties and their global spatial distribution and arrangement, metamaterials with exotic behavior can be achieved. The foundation for mechanical metamaterials comes from the study of cellular solids (*1*), where natural materials, such as wood and bone (*2*), or synthetic materials, such as stochastic foams, are understood as a network of closed or open cells (*3*). In the latter case, edges form a network of beams, and on the basis of the connectivity of these beams and their base material, macroscopic behaviors can be predicted analytically (*4*). It was from this insight that the field of architected materials formed, enabling design of periodic structures with tailorable properties such as improved stiffness over foams at similar density due to higher degrees of connectivity (*5*).

Advances in digital fabrication, specifically additive manufacturing, have enabled these complex designs to be realized. Seminal work demonstrated stiff, ultralight lattice materials (*6*), and has since been improved, resulting in mechanical metamaterials with superior stiffness and strength at ultralight densities (*7*) with multiscale hierarchy (*8*). Benefits of nanoscale features further expand the exotic property parameter space (*9*), and architectures featuring closed-cell plates have shown potential for approaching the theoretical limit for elastic material performance (*10*). Other designs seek to use compliance, which can be attained through internal geometric mechanisms (*11*) or through base materials capable of large strain (*12*). Internal architectures can be designed to transmit or respond to load in other nonstandard ways. Auxetic metamaterials exhibit zero or negative Poisson’s ratio (*13*). Internal, reentrant architectures produce contraction perpendicular to compressive loading, and expansion perpendicular to tensile loading, counter to traditional continuum material behavior (*14*). Chiral metamaterials exhibit handedness based on asymmetric unit cell geometry. These designs produce out-of-plane deformations, such as twist, in response to in-plane loading (*15*).

Nearly all of the aforementioned mechanical metamaterials are made with some forms of additive manufacturing, most of which are summarized in (*16*). These processes vary widely in terms of cost, precision, throughput, and material compatibility. The lower end of the cost spectrum, such as fused deposition modeling (FDM), also tends to have lower performance. Limits of thermoplastic extrusion include layer-based anisotropy (*17*) and errors resulting from build angles for complex three-dimensional (3D) geometry (*18*). Higher performance, and higher cost, processes such as selective laser melting (SLM) use materials such as stainless steel but require nontrivial setup for particulate containment and can suffer from layer-based anisotropy, thermal warping, and geometry irregularity (*19*). Some of the highest performance multiscale metal microlattice production techniques based on lithographic and plating processes are well studied and repeatable but are highly specialized and labor, time, and cost intensive. Polymerization, curing, plating, milling, and etching can require up to 24 hours from start to finish for sample preparation (*6*). Large-area projection microstereolithography (LAPμSL) is capable of producing lattices with micrometer-scale (10^{−6} m) features on centimeter-scale (10^{−2} m) parts (*8*) with significantly improved throughput, but extension to macroscale (>1 m) structures remains out of reach, due to practical limitations in scaling these processes and their associated machines.

The largest structure that can be printed with any given process is typically limited by the build volume of the machine. Therefore, substantial effort is focused on scaling up the machines. Meter-scale FDM platforms (*20*) and larger cementitious deposition machines (*21*) have been demonstrated, and coordinated mobile robots are proposed to achieve arbitrarily large work areas (*22*). However, there is a trade-off between precision, scale, and cost. Commercially available two-photon polymerization machines have resolution on the order of 1 μm (10^{−6} m), build size on the order of 100 mm (10^{−1} m), and cost on the order of $10^{6} per machine (*23*). Macroscale FDM machines boast build sizes of 10^{1} m (*24*) but are unlikely to have better than millimeter (10^{−3} m) resolution. Thus, roughly the same dynamic range (scale per resolution) is offered, but with costs approaching $10^{7} per machine, we see a possible superlinear cost-based scaling of achievable dynamic range. Building large, precise machines is expensive, and due to the inherent coupling of machine performance, size, and cost, there are significant challenges for realizing macroscale (>1 m) mechanical metamaterials with high quality and low cost.

An alternative approach to producing mechanical metamaterials seeks to decouple these aspects and, in doing so, overcome machine-based limitations. On the basis of reversible assembly of discrete, modular components, this method uses mechanical connections to build larger, functional metamaterials and structures out of smaller, mass-producible parts. The first demonstration of this approach used custom wound, centimeter-scale, carbon fiber–reinforced polymer (CFRP) components (*25*), resulting in an ultralight density lattice with improved elastic stiffness performance over the state-of-the-art metallic microlattice (*6*), due to the high modulus constituent material. Following this, larger-scale, octahedral voxel (volumetric pixel) building block units were made using commercial off-the-shelf high-modulus, unidirectional pultruded CFRP tubes connected with injection-molded glass fiber–reinforced polymer (GFRP) nodes, resulting in a macroscale (>1 m), high-performance, reconfigurable structure system (*26*). Following this, entire voxel units were made with injection molding of GFRP, yielding the first truly mass-producible discrete lattice material system with low cost, best-practice repeatability, and high performance (*27*). Discrete assembly offers scalability and functionality not achievable with traditional methods due to process and machine limitations.

Here, we present a construction system for mechanical metamaterials based on discrete assembly of a finite set of modular, mass-produced parts. We demonstrate experimentally the desired metamaterial property for each part type and, combined with numerical modeling results, display other unexpected and useful properties. A modular construction scheme enables a range of mechanical metamaterial properties to be achieved, including rigid, compliant, auxetic, and chiral, all of which are assembled with a consistent process across part types, thereby expanding the functionality and accessibility of this approach. The incremental nature of discrete assembly enables mechanical metamaterials to be produced efficiently and at low cost, beyond the scale of the 3D printer.

## RESULTS

### Subsystem characterization

First, we present the discrete material construction system and show how continuum behavior is achieved through design of the parts and their relative structural performance. Parts are designed to have their local beam properties govern global lattice behavior, resulting in an effective bulk material that behaves as if it were produced monolithically.

A lattice, or a mechanical metamaterial consisting of a periodic network of interconnected beams, can be described, and its performance predicted, analytically. We can describe lattices as stretch or bending dominated, based on how they resolve external forces as a function of their internal beam connectivity, which corresponds to Maxwell’s frame rigidity criteria extended to 3D (*5*). Stretch-dominated lattices, such as the octet, have higher connectivity (*Z* = 12) and higher stiffness to weight than bending-dominated lattices, such as the kelvin, which have lower connectivity (*Z* = 4) (*7*). In this work, we use the cuboctahedra lattice (referred to as cuboct) geometry, which is uniquely positioned between low and high connectivity (*Z* = 8) yet has been shown to have stretch-dominated behavior, both in microlattice implementation (*28*) and as discretely assembled vertex-connected octahedra (*27*).

In Fig. 1 (A to C), we show a new decomposition using face-connected cuboctahedra voxels, which produces the same lattice geometry but has additional benefits to be discussed here. Voxels are discretized into faces, which consist of beams and joints. There are two types of joints: Inner-voxel joints are the points at which six separate faces are joined to form a voxel, and intervoxel joints provide the vertex to vertex connections between neighboring voxels along a single face. A joint consists of nodes, which are the geometric features on the part providing the fastening area, and the fasteners, which are mechanical connectors. On the basis of the material and geometric properties of each subsystem, local properties can be controlled to ensure proper global, continuum behavior. In this case, our lattice should behave as an interconnected network of beams. Therefore, we wish to design joints to have significantly higher effective stiffness and strength than the beams they connect. In this way, the global effective stiffness and strength of the lattice are governed by those subsystems with the lowest relative value.

Following as-molded material characterization to calibrate analytical and numerical models (fig. S1), subsystems were then characterized in tests designed to isolate the critical performance aspects for proper system behavior. Rivets, intervoxel nodes, individual voxels (consisting of beams and inner-voxel joints), and multivoxel assemblies were tested. The specific goal is to quantify the degree to which voxel and multivoxel behavior is governed by stiffness and strength properties of the beams, rather than the joints. Experimental results are shown in Fig. 1D, with axial stiffness and critical load values noted.

Because each subsystem effectively acts across the same cross section (a single voxel), we can directly compare their yield strength using their observed failure loads. We see that the intervoxel node and fastener yield strengths are roughly two and four times the voxel yield strength, respectively. For axial stiffness, we treat single and multivoxel tests as effective springs in series. A single voxel then consists of five effective springs in series: top fasteners, top nodes, voxel, bottom nodes, and bottom fasteners. For springs in series, the equivalent axial stiffness is the reciprocal of the sum of the individual spring reciprocals*k _{i}* and small

*k*

_{1}, we see that

*k*

_{eq}equals

*k*

_{1}, indicating that the governing value is the lower spring stiffness. Using measured values for fasteners, nodes, and voxels, we see that the experimental value for the two-voxel assembly agrees with this analytical description and that both effective stiffness and strength are governed by voxel, and thus beam, properties. Additional details on the joint load paths and hysteresis effects are presented in the Supplementary Materials. Under cycling, the hysteresis rapidly decreased to a stable value, which for the stiffest lattice (cuboctahedral) was approximately twice the base material, corresponding to matching the hysteresis of a rigid rubber at a fraction of a percent of the density (

*53*,

*54*). This can be further reduced with preloaded joints (

*27*).

### Part types

Using this construction system, we present the discretely assembled mechanical metamaterials consisting of four part types: rigid, compliant, auxetic, and chiral, as shown in Fig. 2. Six face parts (Fig. 2A) are assembled to form voxels (Fig. 2B), which are then assembled to form multivoxel lattices (Fig. 2C). Details of the assembly procedure and throughput metrics can be found in the Supplementary Materials.

Rigid voxels resolve external loading through axial beam tension and compression, resulting in elastic, followed by plastic, buckling of struts. Lattices made with these parts show near-linear scaling of effective modulus, positive Poisson’s ratio, and yield strength determined by geometric and manufacturing process parameters. Compliant voxels are designed with corrugated flexure beams, a motif found in flexural motion systems (*29*), which resolve axial beam forces through elastic deformation of the planar flexures. Lattices made with these parts show consistent elastomeric behavior at even single voxel resolution and have a near-zero Poisson’s ratio. Auxetic voxels are designed as intersecting planes of reentrant mechanisms, which expand and contract laterally under uniaxial tension and compression, respectively. Lattices made with these parts show a negative Poisson’s ratio through a combined action of pure mechanism and flexural beam bending. Chiral voxels are designed with an asymmetric mechanism, which responds to in-plane loading by producing either clockwise (CW) or counterclockwise (CCW) rotation. When interconnected in three dimensions, this produces out-of-plane twist in response to uniaxial tension or compression. By combining CW and CCW parts, internal mechanism frustration can be avoided, enabling improved scalability over prior art. The four lattice types and their behaviors will be described in further detail in the following subsections.

### Rigid lattice behavior

The rigid lattice type exhibits relative modulus-density scaling, which matches previous results in literature but does so with a novel geometric decomposition. We present experimental and numerical results for the rigid lattice type in Fig. 3. The characteristic behavior of a unit cell voxel is shown in Fig. 3A. The geometry is isotropic along its primary axes, and it responds to loads through axial beam tension and compression. While individual voxels are dominated by underconstrained, mechanism behavior of the quadrilateral faces, when multiple voxels are joined, there is sufficient connectivity to provide rigidity through triangulation of neighboring voxel faces. As a result, effective modulus increases with increasing cell count, and this value eventually reaches an effective continuum value, as seen in Fig. 3D.

Having established that the global behavior is governed by the beam properties, now, we can correlate analytical models with experimental results for effective lattice behavior. Here, we will look at effective elastic modulus *E** and yield strength σ_{y}, the former corresponding to the linear portion of the stress strain curve under quasi-static loading, and the latter corresponding to the failure load divided by the specimen cross-section area. Stress-strain curves for lattice specimens of cube side voxel count *n* = 1 to 4 are shown in fig. S10, where an initial linear elastic regime is followed by a nonlinear elastic regime and plastic yield. Using load and displacement data, stress and strain values are calculated on the basis of lattice specimen size. The calculated moduli are shown with numerical results in Fig. 3D, in this case using the reduced order beam models as described in Materials and Methods. It can be seen that as voxel count *n* increases, *E** approaches a continuum value depending on the beam thickness, and thus relative density of the lattice. In the case of our built lattice, voxel cubes of side voxel count *n* = 1 to 4 have effective moduli relative to the continuum approximation (horizontal line, value for 10 × 10 × 10 determined numerically) of 9, 56, 73, and 89%, respectively. Discrepancy between experimental and numerical results is also calculated for specimens *n* = 1 to 4 to be 458, 10, 6, and 3%, respectively. This can be attributed to the ratio of internal to external beams increasing as voxel count increases (fig. S7). The internal beams, which are fully constrained and behave as a rigid network, asymptotically govern the effective global behavior.

These predicted effective lattice properties over the range of effective densities are plotted relative to constituent values in Fig. 3E. The slope of the curve connecting these points, plotted on a log/log chart, provides the power scaling value, which is used to analytically predict lattice behaviors at the macroscopic scale (*4*). Effective lattice modulus and density are related to constituent material modulus and density by *b* is 1 for stretch-dominated lattices and 2 for bending dominated. We find *b* = 1.01 for our rigid lattice. This scaling value had been shown previously for the monolithic (additively manufactured) cuboctahedron lattice (*28*) and for discretely assembled, vertex-connected octahedra (*27*), to which we now add our novel lattice decomposition. It should be noted that these effective values are from numerical simulations, not experiment, although we direct the reader to Figs. 3D and 4D for agreement between experimental and numerical results.

Next, we compare experimental yield stress results with analytical predictions of local beam failure based on relative density, as a function of beam thickness *t* and lattice pitch *P*. Here, we will use experimental data from the 4 × 4 × 4 specimen, as this is closest to demonstrating continuum behavior (effective modulus is 89% predicted continuum value). On the basis of the load at failure and lattice material and geometry, we can determine a given beam compressive failure load to be around 88 N. We determine the analytical critical beam load using either the Euler buckling formula or the Johnson parabola limit, depending on the compression member’s slenderness ratio (fig. S5). We determine our beam slenderness ratio to be 29.5, which is over the critical slenderness ratio of 19.7 (see the Supplementary Materials for complete calculation); thus, we use Euler buckling formula. Because the as-molded material properties vary, we determine the critical load to range from 70 to 108 N, with the mean value of 89 N very closely approximating the experimental value. Thus, there is good correlation between both stiffness and strength based on the design of our discrete lattice material.

### Compliant lattice behavior

The compliant lattice type exhibits quadratic scaling for effective stiffness, as well as consistency across voxel counts regarding continuum behavior and elastic limit values. We present experimental and numerical results for the compliant lattice type in Fig. 4. The characteristic behavior of a unit cell voxel is shown in Fig. 4A. While the load paths are topologically the same as the rigid voxel, as this is a function of lattice connectivity, the mechanism through which beams resolve these loads is different. Here, the planar-spring beams deform in combined axial and in-plane bending, as a controllable property of the compliant features we design. This produces several unique behaviors in this lattice type.

First, we can see from the experimental stress-strain curves that for similar strains, the compliant lattice shows linear elastic behavior up until the elastic limit (fig. S10B). The stress at which this transition occurs is consistent across voxel counts, from *n* = 1 to *n* = 4. Second, the effective modulus is also consistent across voxel counts. This is confirmed by simulations using reduced order beam models, as shown in Fig. 4D. Given the large range of linear to nonlinear and individual to continuum behavior seen in the rigid lattice, the compliant lattice is markedly different in its consistency. This behavior is attributable to the spring-like behavior of the beams, a similar observation to analytical models for stochastic foams (*30*). As cube specimen side length voxel count increases, so do the number of springs acting in parallel, which produces an effective spring stiffness *K*_{eff} = *K*_{1} + *K*_{2} + *K _{n}* …. However, as spring count increases, so does effective area, both proportional to side length squared. Thus, a single voxel has the same effective modulus as a 4 × 4 × 4 cube or an

*n*×

*n*×

*n*cube. This effect is reduced as beam-spring amplitude

*a*goes to zero, meaning it shows more asymptotic behavior similar to the rigid cuboct lattice.

Another property observed experimentally, and confirmed numerically, is a low, near-zero, Poisson’s ratio. Figure 4E shows the simulated effective Poisson’s ratios for the compliant and rigid voxel. At the largest compliant amplitude, we see a value of near zero. As the amplitude *a* of the compliant spring feature goes to zero, the Poisson’s ratio converges to around 0.15, which is the effective value for the entire parameter range of the rigid lattice.

Last, this lattice shows near quadratic stiffness scaling, in contrast to the near linear scaling shown by the rigid lattice, while having the same base lattice topology and connectivity as the rigid version (Fig. 3E)—meaning it has bending-dominated behavior with a stretch-dominated lattice geometry. The range of spring amplitudes as a function of lattice pitch *P* shown in Fig. 3E is *a* = 0.05, 0.1, 0.15, and 0.2, and these have scaling values of *b* = 1.72, 1.89, 1.93, and 1.95, respectively. This is attributable to the localized behavior of the spring-like beams. Whereas in the rigid lattice vertically oriented beams in compression are offset by horizontally oriented beams in tension, resulting in stretch-dominated behavior, here, global strain is a function of local spring-beam strain, which does not produce significant reactions at beam ends opposite an external load.

### Auxetic lattice behavior

The goal of the auxetic lattice type is to exhibit a controllable negative Poisson’s ratio. We present experimental and numerical results for the auxetic lattice type in Fig. 5. The characteristic behavior of a unit cell voxel is shown in Fig. 5A. Because of the internal architecture, which consists of interconnected, reentrant mechanisms seen elsewhere in literature (*14*), the cell responds to axial strain with a similarly signed transverse strain, resulting in a negative Poisson’s ratio ν, where *d* as a function of lattice pitch *P*, as shown in Fig. 5D.

Experimental results are shown in Fig. 5B. Lattice specimens are cubes of voxel width *n* = 1 to 4. Specimens were compressed to identical strain values (ϵ_{axial} = 0.2), and transverse strain was measured by visually tracking points using fiducials mounted to the nodes along transverse faces (*yz* plane) parallel to the camera. Experimental data can be found in fig. S10C. These points are slightly obscured due to reduced reentrant behavior at the edges of the lattice. In Fig. 5C, we show contour plot element translation in the *y* direction, which is out of plane and normal to the camera view. While this behavior is generally isotropic, it should be noted that the effect of the internal mechanisms is reduced at the corners/edges of the cube specimen, as shown in Fig. 5F. The median effective strain values are plotted in Fig. 5E over the range of parameters shown in Fig. 5D. The median was chosen to reduce the influence of the boundary conditions where ν ≈ 0. The experimental Poisson’s ratios, indicated as black squares, were measured using fiducial targets and motion tracking at the points indicated in Fig. 5B.

There are two main insights from this study. First is that the effective metamaterial behavior approaches a nominal continuum value as cube side length of voxel count *n* increases. For any reentrant distance, this behavior can be attributed to the increase of internal mechanism architecture relative to boundary conditions. Boundary conditions increase as a function of surface area proportional to *n*^{2}, while internal mechanism architecture increases as a function of specimen volume proportional to *n*^{3}. For lower values of *d*, the single voxel demonstrates lower values for Poisson’s ratio (increased auxetic behavior) compared to multivoxel specimens, but this is strongly influenced by boundary conditions and can be considered an outlier.

The second insight is that the effective Poisson’s ratio decreases (becomes more negative) as reentrant distance *d* is increased, for voxel specimens larger than *n* = 1. This can be understood by considering the continuous beams of the reentrant faces as a pseudo–rigid body model, where continuous flexural mechanisms are discretized as effectively rigid links connected by planar joints with torsional stiffness (i.e., a spring) (*31*). As *d* decreases, so does link length, causing less clearly defined boundaries between the rigid link and compliant spring joint (see the Supplementary Materials for further analysis). As a result, the rigid link behavior begins to dominate, causing higher overall effective stiffness and lower compliance, thus reducing the reentrant mechanism efficacy.

### Chiral lattice behavior

The chiral lattice type exhibits scalable twisting behavior, which is attributable to having two chiral part types and developing a construction logic to avoid internal frustrations. We present experimental and numerical results for the chiral lattice type in Fig. 6. The characteristic behavior of a unit cell voxel is shown in Fig. 6A. On the basis of the chirality orientation, the cell will respond to an axial strain with a macroscopic twisting in either the CW or CCW direction, in the plane normal to the direction of loading (i.e., loading in *z* direction results in twisting in *xy* plane). The effective chirality can be measured as degrees twist per unit strain.

Experimental results are shown next to their numerical simulations in Fig. 6E. Lattice specimens are designed as columns with 1:4 width to height ratio, similar to (*15*). The top half is CCW chiral lattice, and the bottom half is CW chiral lattice. This produces the largest net twist at the rigid interface between the two halves and allows fixed boundary conditions at the top and bottom. Chiral columns of 1 × 1 × 4, 2 × 2 × 8, and 3 × 3 × 12 were tested in compression to identical strain values (ϵ_{axial} = 0.05), and twist was measured by tracking a single point at the center of the lattice. Experimental results are shown in fig. S10D. The 1 × 1 × 4 column shows larger values for twist than does the 2 × 2 × 8 column. This is attributable to internal architecture, which is also the cause of the scalable twisting found over a range of beam sizes.

Experimental values for twist per strain are shown next to reduced order beam model simulation results in Fig. 6B, over a range of values for radius *r* of the face part as a function of lattice pitch *P*, with increasing column voxel width *n*. We observe an increased twist per axial strain for smaller values of *r*. This is attributable to the direct relationship between strain and twist as a function of the rotational mechanism. If we assume that a unit strain is translated into an arc length *s*, then the rotation angle θ increases as circle radius *r* goes to zero. However, given a nominal beam thickness *t*, there is a limit to how small *r* can become before the mechanism becomes ineffective. See the Supplementary Materials for further analysis.

There are several key takeaways from this. First, we see that performance does not decrease monotonically with increasing voxel count *n* but rather stabilizes to a continuum value. This is in contrast to comparable results in literature (*15*) and can be explained by looking more closely at the combination of CW and CCW part types. Done properly, internal frustrations—when CW and CCW faces are joined, they essentially cancel each other’s twist, resulting in zero twist per strain—can be avoided, as shown in (*32*) by using voids. In our case, we get improved twist performance by designing the internal architecture according to rules chosen to avoid frustration. This means that voxel types are directionally anisotropic, in contrast to the previous three lattice types, and further are spatially programmed to produce desired global effective behavior. Strategies for this spatial programming are shown in Fig. 6C. On the left, we show a beam with odd number voxel widths. Here, design rule #1 is to orient the net face chirality (represented as arrows) away from the column interior. The experimental lattices for *n* = 1 and *n* = 3 widths were built using rule #1. Design rule #2 was developed starting from *n* = 2, where the orientation of interior faces is ambiguous when following rule #1. Rule #2 introduces continuous, CW circumferential orientation of the interior chiral faces and was used in construction of the *n* = 2 experimental articles. Both rules are hierarchical, e.g., a rule #1 5 × 5 column contains a 3 × 3 and 1 × 1 column in its interior as shown in Fig. 6C. Simulations were performed for all column widths using both rules and show decreased twist response for rule #2, in agreement with experimental measurements. These rules were determined empirically and are not considered exhaustive but indicate the importance of rational design in this lattice type.

## DISCUSSION

Here, we presented a method for producing large-scale mechanical metamaterials through discrete assembly of modular, mass-produced parts. We showed that bulk, continuum behavior can be achieved through design of the parts and connections, ensuring that global behaviors are governed by local properties. We presented a finite set of part types, which exhibit a diverse range of behaviors. Rigid lattice types show linear stiffness-to-density scaling with predictable failure modes. Compliant lattice types show quadratic stiffness-to-density scaling, as well as unique bulk behavior at low cell count, such as near-zero Poisson’s ratio. Auxetic lattice types show controllable, isotropic negative Poisson’s ratio. Chiral lattice types show scalable transverse twist in response to axial strain, which is a result of two part types being used to prevent internal architectural frustration. All four part types showed good agreement with numerical results, and their behavior is predictable through analytical means. All lattice types are made the same way: Parts are injection-molded and assembled to make voxels, and voxels are similarly joined to build lattices. This is a low-cost, highly repeatable process that promises to enable mechanical metamaterials at macroscales (fig. S13).

There are several advantages resulting from discrete assembly, which make it stand out from existing fabrication methods currently available for producing metamaterials, which include increased functionality, repairability, reconfigurability, and scalability. While this work presented mechanical metamaterials, discretely assembled electromagnetic materials have been previously demonstrated. Passive and conductive parts have been assembled into heterogeneous, functioning 3D circuitry (*33*), and rigid, flexural, and actuated building block parts were used to assemble modular microrobots (*34*). These are millimeter-to-centimeter scale parts, and the extension of this approach to larger scales is expected to enable, mesoscale cellular robots. Because of the discrete nature of the construction, damaged or broken parts can be removed and replaced. This was demonstrated in prior work (*27*), where lattice specimens were tested to initial failure (plastic beam buckling and rupture) and then unloaded, the damaged voxel unit was removed and replaced, and then the specimen was tested again. Repaired specimens showed only 1.5% loss of effective stiffness and 5% loss of effective strength. Quasi-static reconfigurability was demonstrated through the assembly, disassembly, and reuse of macroscale (225-mm pitch) octahedral voxels (*26*). In that case, over 125 voxels were used to build a 5-m bridge capable of holding several hundred kilograms, then these were reconfigured into a boat, and then these were again reconfigured into a shelter. Scalability has been demonstrated in prior work, where over 4000 injection-molded octahedral voxel units were assembled into a 4.25-m wingspan ultralight lattice aerostructure (*35*). The parts were manually assembled, with a mass and volumetric throughput that was competitive with typical mesoscale additive processes such as SLM and FDM. The machine cost and process challenges associated with making such a lattice structure with either of those methods highlight the benefits of this approach. Scaling to part counts above 10^{3} will benefit greatly from assembly automation. Stationary gantry platforms have been fitted with end effectors for voxel transport and bolting operations (*36*), and mobile robots have been implemented to perform similar operations while locomoting on the lattice as they construct it (*37*). Stationary systems promise high throughput for a bounded work envelope, while mobile robots can be parallelized and require no global positioning due to local alignment features, offering benefits of autonomy and reliability. Automation will be critical for producing these metamaterials and structures in large quantities envisioned for commercial applications.

Injection molding as used here offers low cost and high repeatability, but it immediately limits which constituent materials can be used. Sheets of material could be used with subtractive processes such as milling, laser, or water-jet cutting to make voxel face parts, although redesign of the joints would be needed. Prior work has shown successful lattice production this way, using a snap fit connection, which needs a final adhesive or thermal bonding step to remove the final degrees of freedom at the joints (*38*–*41*). Natural materials such as wood can be used this way, and in the future, moldable bio-based resins with natural fibers are expected to be commercially available. Looking at scaling down our process, there are some practical limitations to both the part production and the assembly. Scaling down the parts by an order of magnitude (from 75-mm cell pitch to 7.5-mm cell pitch) should be possible based on current best practice microinjection molding and existence of commercially available microfasteners (see the Supplementary Materials for details). Scaling down further (submillimeter cell pitch) would require novel part production and joining methods, suggesting that this may be a regime where conventional additive processes are preferable. Rather than focus on absolute length scale, for our metamaterials, we are concerned with the ratio of cell size to smallest characteristic system size. Given the quasi-static loading in our case (*42*), we easily achieve subwavelength cell size while also demonstrating effective continuum properties as a function of local cellular architecture. Thus, the ability to compose macroscopic metamaterials blurs the boundaries between material and structure.

Last, we limited our study to a set of four distinct behaviors, shown as separate homogeneous lattices. Comparable demonstrations of these properties exist in prior art, but each has typically entailed dedicated development, whereas here we show a single scalable system capable of achieving this range with a consistent production process based on discrete assembly. Because of this, heterogeneous lattices can be made with this approach just as easily. Heterogeneous metamaterials have been shown to offer exponential combinatorial possibilities (*43*), as well as the ability to realize any arbitrary elasticity tensor (*44*). Furthermore, the design of new part geometries with blends of behavior is a promising next step for use in assembling spatially graded heterogeneous structures, which is one of the main benefits sought through additive processes (*45*) to achieve functionality seen in natural systems (*46*). By offering a simple yet diverse set of parts unified with a consistent assembly method, this work represents a significant step in lowering the barrier for entry to realizing the promise of mechanical metamaterials, especially for macroscale applications. Combined with hierarchical design tools and assembly automation, we foresee this research enabling emerging fields such as soft robotics, responsive aero- and hydrodynamic structures, and user-defined programmable materials, thereby further merging the digital and physical aspects of future engineering systems.

## MATERIALS AND METHODS

### Injection molding and assembly

Part production and assembly details are shown in fig. S1. Parts were injection-molded by Protolabs, a U.S.-based on-demand manufacturing service provider. To ensure low cost, parts were designed to be two-part moldable. While this is simple for most of the part, the inner-voxel tab and hole at 45° required a custom-designed opening, as shown in fig. S1C. Parts were assembled with ^{3}/_{32}-inch-diameter blind aluminum rivets, using a pneumatic rivet gun. The voxel assembly process is shown in fig. S1D. Voxel to voxel joints used the same process, as shown in fig. S1E. Metrics for assembly time and throughput are shown in table S1.

### Mechanical characterization

Small-scale tests to validate continuum behavior as shown in Fig. 1 were performed on an Instron 4411 testing machine using a 5-kN load cell. Lattice specimens for each type were tested in cubes of side length voxel count *n* = 1, 2, 3, and 4. Lattice tests were performed on an Instron 5985 testing machine using a 250-kN load cell. Specimens of a given lattice type were loaded to the same amount of relative strain at an extension rate of 10 mm/min. Both machines use Bluehill 2 software for data acquisition. Video was recorded using a Nikon D3400 camera. Video was analyzed using Tracker, an open source video analysis and modeling tool (https://physlets.org/tracker/).

### Numerical modeling

Fully meshed finite element analysis (FEA) simulations were used to check stress concentrations, but these typically incur higher computational costs (figs. S5 and S6), and therefore were limited to under 10 voxels. A static stress analysis solver based on NASTRAN was used in Autodesk Fusion 360’s built-in simulation environment. Larger lattice models were simulated using the Frame3DD library, a freely available numerical solver implementing Timoshenko beam elements (http://frame3dd.sourceforge.net/) along with a python interface, pyFrame3DD (https://github.com/WISDEM/pyFrame3DD). For analysis of asymptotic behavior of large lattices, Frame3DD was modified to incorporate sparse matrix math using CHOLMOD from the SuiteSparse library (https://github.com/DrTimothyAldenDavis/SuiteSparse). Python utilities were written to automate creating nodes, edges, faces, and voxels, as well as applying loadings and boundary conditions using spatial rules (e.g., fixing the bottom of a lattice and applying forcing to the top nodes). These simulations were validated against a commercial software with comparable sparse matrix solving capabilities (Oasys GSA v9.0).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/47/eabc9943/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported by NASA Cooperative Agreement #80NSSC19M0039, U.S. Army Research Lab Cooperative Agreement #W911NF1920117, and CBA Consortia funding.

**Author contributions:**B.J. designed and produced parts. C.C. developed code to perform sparse matrix numerical modeling. F.T. performed numerical modeling using best-practice commercial packages. A.P.R. performed mechanical testing of lattice and subsystems. M.O. led lattice test specimen assembly. N.G. provided system architecture guidance.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).