The role of unit cell topology in modulating the compaction response of additively manufactured cellular materials using simulations and validation experiments

Additive manufacturing has enabled a transformational ability to create cellular structures (or foams) with tailored topology. Compared to their monolithic polymer counterparts, cellular structures are potentially suitable for systems requiring materials with high specific energy-absorbing capability to provide enhanced damping. In this work, we demonstrate the utility of controlling unit-cell topology with the intent of obtaining a desired stress–strain response and energy density. Using mesoscale simulations that resolve the unit-cell sub-structures, we validate the role of unit-cell topology in selectively activating a buckling mode and thereby modulating the characteristic stress–strain response. Simulations incorporate a linear viscoelastic constitutive model and a hyperelastic model for simulating large deformation of the polymer under both tension and compression. Simulated results for nine different cellular structures are compared with experimental data to gain insights into three different modes of buckling and the corresponding stress–strain response.

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Cellular structures (or foams) are increasingly being considered as materials for lightweighting in auto, aerospace, construction, and defense industries as a way to achieve better fuel and energy efficiencies in logistics flow.Additive manufacturing (AM) has enabled significant control over the geometric parameters of a unit cell, thereby facilitating an expanded space for the design of a cellular structure (or cellular material) resulting in a wider range of possible mechanical responses [1][2][3].Current methods of vat photopolymerization (or commonly stereolithography (SLA)) have enabled the fabrication of an expanded variety of complex structures and geometries not possible through other forms of manufacturing.The ability to design the mechanical response of cellular structures by varying the geometrical aspects of the unit cell has garnered growing attention.
Stochastic foams and periodic cellular structures have emerged as cellular materials of particular interest.Stochastic foams have a complex microstructure consisting of an interconnected network of cell walls, commonly forming closed pores, resulting in randomly packed cells as a consequence of the foaming process [4].They have excellent light-weighting ability and are widely used in aerospace structural applications, impact mitigation in vehicles, packaging, and applications requiring impermeability to gas and fluids [5,6].In contrast to stochastic foams, open-cell foams have an interconnected network of ligaments (or surfaces), forming open pores.Specifically, periodic cellular materials have a well-determined periodic geometry that can be controlled by a small number of design parameters [7,8].They have excellent energy absorption and damping characteristics.
Cellular materials, when subjected to large unconfined compression, exhibit a non-linear stress-strain response with characteristic elastic, collapse, and densification regimes.These materials possess the unique capability to store and dissipate energy, thereby making them useful for applications requiring cushioning or damping phenomena in diverse industrial applications.Large deformation analysis not only helps in determining the cellular solid's loadbearing capacity (load at the initiation of collapse or yield), toughness, and overall response to applied forces but also allows for the evaluation of its ability to undergo buckling and recovery after removal of the load, which is essential for designing safer and more reliable products [9][10][11][12].
Cellular materials can display multiple modes of deformation.Examples of such modes of deformation are bending, buckling, shearing, and folding.Prior experiments on cellular materials have commonly shown the formation of a localized compaction band accompanied by the collapse of a few cells within a narrow region of the specimen.This localized collapse was followed by a progressive collapse of the entire cellular structure, facilitated by a cellcompaction (collapse) wave that traverses the length of the specimen [13].It is possible to describe the non-linear stress-strain response using force balance laws of the classical mechanics [14,15] combined with viscoelastic constitutive response of the underlying material, e.g.polymer, however, such techniques are limited to simple structures undergoing modes of deformation that can be described analytically.With the advancements in AM, the deformation analysis of more sophisticated topologies (such as Gyroid) is more challenging to model using a closed-form solution alone.In this study, advanced finite element (FE) simulations were performed that enabled an understanding of the deformation behavior of AM cellular materials at the mesoscale, i.e. unit-cell level, which is crucial for predicting and subsequently optimizing their mechanical performance.
Several previous studies have investigated the effect of unit-cell topology on the mechanical behavior of cellular structures [16][17][18][19][20].While many studies have focused on a few topologies with selective variations in the unit cell such as thickness, strut angle, and density, limited progress has been made in understanding the complex mechanisms of deformation that lead to an expanded range of mechanical responses from a wider range of unit-cell topologies.In this work, we investigated how varying the structure of the unit cell activates different mechanisms of deformation, leading to strikingly different stress-strain responses of the cellular structures.
Polyurethanes offer versatility, toughness, heat resistance, and durability, making them ideal for engineering-grade materials [21,22].They outperform other polymers, such as polylactic acid, polyamides, thermoplastic elastomers, polyolefins, and acrylonitrile butadiene styrene (ABS), in terms of semi-flexible material performance [23][24][25].Selecting the right feedstock and identifying an appropriate cellular structure for a specific application demands an understanding of the process-structure-property relationship that can be used to optimize the performance of the final material for specific part requirements [26][27][28][29].
The primary focus of this study was to understand the effect of mesoscale topology on the compressive response and compare the energy-absorption capability of polymer-based cellular structures.To achieve this, we printed nine different cellular structures using vat photopolymerization with polyurethane-based resin.Our approach involved subjecting all structures to uniaxial (unconfined) compression tests to acquire their engineering stress-strain response.Mesoscale simulations provided insights into the mechanics of non-linear deformation.Simulations resolved the shape and modes of deformation at the unit-cell level to assist in identifying mechanisms of buckling that were correlated with the measured stress-strain response and imaging from experiments.Integrated simulations and experiments were then used to understand the effect of cell topology on energy absorption capacity and efficiency of the cellular structures under compression.

Material and methods
We aim to characterize the mechanical properties of AM cellular materials and the associated mesoscale response through a combination of approaches involving computational simulation, AM technique, and compression testing.Methods used in this investigation are detailed in the subsequent sub-sections.

Cellular structures
From an extensive array of possible topologies, we selected six strut-based cellular structures namely simple cubic (SC), Kelvin, Gibson Ashby (GA), octet, body-centered (BC), and diamond, and three triply periodic minimal surfaces (TPMS) namely gyroid, Schoen IWP, and Schwarz primitive.These topologies are widely used while adopting cellular materials in applications for cushioning, damping, and energy absorption [30][31][32][33][34][35][36].In this study, we fixed the material volume fraction (VF) to 20% and the overall part dimensions to 30 mm × 30 mm × 30 mm for all the structures.This sample size was chosen considering the feasibility of printing, testing, and post-processing the samples.A smaller sample size would compromise the quality of the printed internal cell structures due to the resolution of the printer while making the sample size too large would increase the printing time and post-processing time in addition to the high computation requirements to run the full-scale simulations.
Computer-aided design (CAD) models were developed for the aforementioned cellular structures using the parametric design platform in Grasshopper (a plugin for Rhinoceros 7) [37].We constructed unit cells for each topology; the unit cells were then replicated along three orthogonal coordinate axes to fulfill the requirements of the overall part dimensions.For all strut-based structures, the unit cell models were created by assigning thicknesses to the wireframe construction (with 2.5 mm cell size) resulting in the target VF of 20%.TPMS-based structures were generated from the isosurfaces of the nodes that satisfy the nodal approximation equations.The nodal approximation equations for Gyroid (Φ Gyroid ) [20], Schoen IWP (Φ IWP ) [38], and Schwarz Primitive (Φ Primitive ) [39] are respectively given in equations ( 1)-( 3): In the aforementioned equations, x, y, and z are the coordinates of the nodes on the isosurface, and 'L' denotes the unit cell length of the TPMS structure (refer to figure 1).Initially, the size 'L' of the TPMS structures was set to 2.5 mm, similar to the strut-based structures.However, this resulted in too fine of features for successful printing and post-processing of TPMS structures on the selected AM system (explanation follows).Hence, these TPMS were also printed with a larger unit cell size (L = 5 mm).Consequently, the number of unit cells was 6 × 6 × 6 for the three TPMS structures, compared to 12 × 12 × 12 for all other structures.
The unit cell solid models and the CAD of macro-architectures obtained by replicating the unit cells in the orthogonal directions are shown in figure 2. All CAD models were converted to stereolithography (STL) formats, which were used for printing.The STLs were additionally used to create FE models for FE simulations of compression.Detailed methodologies employed for conducting both the experiments and simulations are provided in the subsequent sections.

Additive manufacturing and testing
The cellular materials were printed using SLA, a vat photopolymerization method, on a Formlabs Form3+ (Formlabs, Somerville, MA, USA) using Elastic 50A resin, a photopolymer resin classified as polyurethane.Default printing parameters were adopted and controlled by Preform software version 3.24.2.The layer thickness was set to 100 µm and printed without using any external or internal support features.The printed solids were washed with isopropyl alcohol (IPA) (Fisher Scientific, Waltham, MA, USA) to remove all excess resin from the surfaces.They were submerged for three days with the IPA replaced daily.We then dried them in a convection oven (BINDER TM , Bohemia, NY, USA) at 30 • C for 30 days.Each structure was printed in sets of six, and from these, three samples without visible defects were selected.30 mm × 30 mm × 30 mm bulk specimens were also designed and printed as a control.The printed bulk specimens were washed with IPA for 10 minutes to remove any residual resin on the surface.After that, they were dried in a convection oven at 30 • C for 3 days.
A Keyence VHX-6000 digital microscope (Osaka, Japan) was used to characterize and verify the strut size and overall conditions of the post-processed cellular materials.The observations were conducted at 20X, 30X, and 50X magnification with full ring lighting.Compression testing (unconfined) was performed to evaluate the cellular materials by using an Instron 3343 (Instron Bluehill, Norwood, MA, USA) with a 1 kN load cell.The displacement rate was set as 0.2 mm s −1 and the maximum stress limit was set to 0.1 MPa.A total of 4 compression cycles were applied to each sample prepared in accordance with the protocols in ASTM D575-91 standard.The four consecutive cyclic compression tests were intended to determine whether or not there was any appreciable change in the stress-strain response after multiple-cycle testing.Figure 3 provides an overview of the key steps in the experimental process.

Finite element modeling
Advanced FE methods (FEMs) available in SIMULIA's Abaqus package (version 2020) were used to simulate and gain insights into the nonlinear material response of the cellular structures.FEM has been successfully used for simulating cellular materials [40][41][42][43].To achieve accurate results, high-fidelity simulations require a fine discretization (good quality mesh) that will resolve the sub-structures of the unit cell, an accurate definition of the material properties combined with an appropriate constitutive model, and relevant boundary conditions.In this study, we used a hexahedral element to create the high-quality mesh for all nine cellular structures, calibrated a hyperelastic constitutive model and viscoelastic material properties of the pristine material, and performed simulations under unconfined compression using Abaqus/Explicit that solves elastodynamic governing equations.Details of the workflow in developing FEM models are as follows: 2.3.1.Finite element mesh generation.Mesh discretization was accomplished using Cubit version 16.06.0.Three-dimensional 8-noded hexahedral meshes (Element C3D8 within Abaqus) were created from STL files using the overlay grid method called Sculpt [44] (a companion application in Cubit version 16.06.0).We chose an 8-noded hexahedral element (C3D8) as it is commonly preferred over a 4-noded tetrahedral element (C3D4) when modeling soft, nearly incompressible materials because the first-order tetrahedral element exhibits a tendency for volumetric locking, commonly leading to an artificial increase in stiffness [45].While higher order tetrahedral element (C3D10) can mitigate the tendency for volumetric locking, it tends to be computationally expensive in terms of memory requirement and time required  To assess the quality of the FE mesh, the scaled Jacobian (sJ) measure for the generated meshes was calculated.A lower value of sJ indicates a poor-quality mesh, and sJ < 0 signifies an invalid or inverted mesh.Notably, the sJ of the elements in the generated mesh exceeds 0.2, which is a satisfactory threshold for accomplishing FE simulations without excessive mesh distortion or degeneration errors that are commonly encountered during large deformation scenarios.

Calibration of the constitutive models.
A hyperelastic constitutive model was calibrated from uniaxial tension and compression experiments.The tensile tests were performed with ASTM D412 Die C specimens measuring 1.91 mm in thickness and a gauge length of 50 mm.The specimens were printed flat on the build plate (in the XY direction).Printed specimens were then washed with IPA for 10 minutes and subsequently dried in a convection oven at 30 • C for 3 days.Tensile tests were performed on three printed specimens using an Instron 3343 with a 1 kN load cell.Tests were performed to failure using a displacement rate of 0.2 mm s −1 .The nominal stress-strain data from three tensile tests were averaged.For the compression experiment, three cubes measuring 30 mm × 30 mm × 30 mm were subjected to compression using a platen velocity of 0.2 mm s −1 .Stress-strain data combined from the tensile and compression experiments was used to calibrate a hyperelastic model.Several different hyperelastic models were assessed, including Neo-Hookean, Yeoh [46], Ogden [47], and Mooney-Rivlin [48,49] for their ability to faithfully reproduce experimental data (refer to supplementary table S1 and figure A2).Amongst these models, the first-order polynomial (Mooney-Rivlin) hyperelastic model was chosen based on the high accuracy of fit and inherent stability.The strain energy potential (U) for the Mooney-Rivlin hyperelastic model is given by [50]: where C 10 , C 01 and D 1 are material constants that require calibration.Ī1 and Ī2 are the first and the second invariants of the deviatoric component of the left Cauchy-Green deformation tensor and J el is the Jacobian of deformation gradient tensor.The initial shear modulus is µ 0 = 2(C 10 + C 01 ) and the initial bulk modulus is K 0 = 2/D 1 .Further, the Poisson's ratio is related to initial values of µ 0 and K 0 by ν = (3K 0 /µ 0 − 2)/(6K 0 /µ 0 + 2).The values of material constants obtained from the calibration process were C 10 = 139 019.19 Pa, C 01 = 345 441.15 Pa at room temperature.In our simulations, we assigned a Poisson's ratio of 0.48, and consequently D 1 = 8.37 × 10 −08 1/Pa, K 0 ≈ 23.9 MPa and K 0 /µ 0 = 37.Note that Abaqus commonly assumes K 0 /µ 0 = 20, corresponding to a Poisson's ratio of 0.475 for some elastomers, which is closer to our approximation.Viscoelastic material characterization was accomplished using a stress relaxation test.Three full-size ASTM D412 Die C specimens were printed for the test.Following printing, the specimens were post-processed similarly to the tensile specimens.The relaxation test was carried out on an Instron 3343 with a 1 kN load cell, where the specimen was strained to 21% and held for 3.5 h.The two-term Prony series was calibrated from the measured relaxation test data.In the form of a Prony series, the time-dependent normalized shear modulus (g R ) is: Here, g i and τ i are normalized shear modulus and relaxation time for ith term in the series.The calibrated parameters are g 1 = 3.19 × 10 −2 , τ 1 = 20.4 s, g 2 = 3.8 × 10 −2 , and τ 2 = 3879.2s.The calibrated hyperelastic and viscoelastic models are plotted along with the experiment data in figure 5.

Boundary conditions and contact properties.
The FE discretized geometries of cellular structures with aforementioned material definitions were positioned in between two rigid platens with platen motion along the Z-axis.One of the platens (lower on the Z-axis) was fixed, i.e. zero displacements and rotations were imposed on all six degrees of freedom.The other platen (higher on the Z-axis) was assigned a constant velocity along the Z-axis and was constrained from moving in other directions, i.e. zero displacements and rotations were imposed on the remaining five degrees of freedom.
The compression experiments were conducted at a velocity of 0.2 mm s −1 , which corresponds to a nominal strain rate of 0.006 67 s −1 .At such a low strain rate, it is impractical to run explicit simulations due to the constraint on the maximum FE time-step, which in our case would require a simulation run-time of over a month.At such a low strain rate, the inertial effects can be neglected assuming quasi-static approximation and the FE simulation can be accelerated using either mass scaling or time scaling to achieve reasonable computation with negligible effect on the overall mechanical response.Since the material under consideration exhibits time-dependent behavior, we opted for semi-automatic mass scaling with the scale factor of 10 6 wherein the kinetic energy was negligible (<4%) compared to the internal energy, thus ensuring the validity of the quasi-static approximation.We used the 'general contact' algorithm [50] that essentially handles the contact between strut to strut, strut to surface, and between all other components of the cellular materials and with the platen.The default 'hard' contact pressure-overclosure relationship was used to define the normal behavior and 'penalty' friction formulation with a friction coefficient of 0.75 [51] was used to define the tangential behavior of the contact property.

The effect of post-processing on the quality of the printed structures
Optimizing the post-processing stage after printing is as crucial as the choice of resin.This step for soft polymer resins like polyurethane plays a critical role in determining the quality and mechanical properties of the final specimen [52].We established a consistent post-processing method for the cellular and bulk samples by controlling three key factors: washing time, drying temperature, and duration.These parameters were qualitatively evaluated based on the following criteria: tackiness (due to residual resin), appearance (including changes in color and quality), equilibrium weight (after residual IPA evaporation), and variability in the compressive stress-strain response (refer to supplementary table S2).Based on these parameters, a 3-day wash in IPA and 30 days at 30 • C was arrived at as the best process for our cellular structures.The washing condition for the bulk specimen was simplified to 10 minutes of wash in IPA and 3 days of drying at 30 • C. The cellular materials required a considerable time for washing and drying to remove the residual resin inside of the structures completely.Comparatively, the bulk specimen required much less time due to the lack of pore spaces.This post-processing approach allowed us to reduce structural defects and warping of the cellular structures, such as tilting and misalignment in certain areas while ensuring reliable results from compression tests.Following these steps, the printed cellular materials exhibited a linear shrinkage of approximately 6% compared to their original size, all while maintaining a VF of 20%.This shrinkage occurs during the post-curing process of remaining acrylate and methacrylate monomers.It is facilitated by strong covalent bonds between carbon atoms in different monomer units, providing an alternative long-range connection mechanism influenced by van der Waals forces [53,54].The isotropic nature of the shrinkage suggests likely mechanisms involving volume swelling and evaporation of residual liquid, however, understanding the exact causes is out of the scope of this paper.
Figure 6 displays the comparison of the thicknesses of the STLs and the printed samples.The dimensions of the STLs were determined from unit cell CAD models whereas the dimensions of the printed samples were approximated from optical images (refer to supplementary figure A3) and averaged.Only the dimensions at two locations of each printed sample were used for the comparison.The differences between the thickness of the STL and the printed samples were significant in the 12 × 12 × 12 (unit cell size 2.5 mm) structure and go as high as 48% in the IWP structure mainly due to very small thickness requirements (145 µm for IWP STL) considering the resolution limits of the printer.Furthermore, at this size, not only did the shape of the cell deviate from the model in certain areas, but the number of unit cells per cell layer also increased, while the inter-sheet spacing decreased [18].The reduction in inter-sheet spacing makes it challenging to completely remove the residual resin internally.Therefore, significant printing defects were observed in 12 × 12 × 12 TPMS (refer to supplementary figure A4).Hence, the 6 × 6 × 6 TPMS structures were also printed that have unit cell sizes and thicknesses twice as large as the 12 × 12 × 12 TPMS.The compression experiments were only conducted on the 6 × 6 × 6 TPMS structures.Figure 7 displays the high print quality and washed internal cavities and struts of the nine structures after post-processing.
AM printing of TPMS structures for the cell size of 2.5 mm resulted in chunks at the bottom of the printed structure.This critical print error affected characterizing the structure's representative compression behavior, and further optimization of the print parameters would be needed to resolve this error in future studies.Furthermore, we refrained from using supports for this printing process because isolating their effect from the overall stress-strain response or removing supports after printing was unfeasible.Cycling through the stress-strain response was crucial before obtaining a representative stress-strain response to preclude irreversible inelastic deformation [55].

Experimental results on compression and hysteresis response
Each printed cellular sample was compressed to 0.1 MPa and released to the original uncompressed platen position.The compression and release cycles were repeated four times on the same sample without a delay between each cycle.A digital caliper was used to measure the dimensions of the three samples before and after each compression test.The average dimensions were 28.15 ± 0.1 mm before testing and 27.94 ± 0.1 mm immediately after testing.The slight difference in dimensions was attributed to incomplete viscoelastic recovery.The thickness was measured again 30 minutes after completion of the compression test, and the cellular samples were able to recover their original dimension completely.
Figure 8 shows stress-strain curves for four consecutive compression cycles for one of the three samples of each topology.These stress-strain data account for large-deformation, viscoelastic behavior exhibited by the pristine cellular material.As the specimen undergoes subsequent cycles of compression, the stress-strain curve demonstrated repeatability for both the loading and unloading responses [56].The stress-strain curves from the 3rd cycle to the 4th cycle are nearly identical for all nine structures.Therefore, we used the stress-strain curves of the 4th cycle as a representative stress-strain response during the loading and unloading of the cellular structure.
The loading and the unloading parts of the stress-strain curve do not overlap, resulting in a significant area contained within the loop, termed hysteresis.Hysteresis is an indication of the amount of energy stored in the pristine material and sources of energy dissipation, such as friction.Comparison of stress-strain between different structures indicates both stress-strain response and hysteresis is affected by the mesoscale structure.
The stress-strain curves for the compression of the cellular structures of SC, Kelvin, and GA samples exhibit an initial linear elastic response followed by a non-linear regime at the yield point.What is common about these three structures is their manifestation of global buckling when first stressed, transitioning thereafter to local cell buckling behavior.This coincides with the yield transient point.This mechanism corresponds to overcoming the resistance to relative motion between the sub-structures or ligaments of the unit cell, allowing the ligaments to move against each other through non-affine transformation and create a change in shape leading to cell collapse.
After the yield point, a stress drop occurs with increasing strain, termed strain softening.The heterogeneous nature of deformation associated with cell collapse has a noteworthy characteristic.The collapse zone did not spread homogeneously throughout the entire material right away.Instead, the deformation intensified in the local collapse zone while the cells in the vicinity deformed in a non-affine manner, exhibiting the non-local influence of the collapse band.As a result, the deformed specimen ceases to retain its original shape.This non-affine deformation leading to the shape change process gradually extends to the entire material, leading to the transition of the stress-strain data from the plateau to densification.Simulated results presented in the subsequent section 3.3.1 provide further insights at the mesoscale, inclusive of the insights for the other cellular structures [57,58].

Finite element simulation results on compression of the cellular structures
The FE methods were successful in simulating the compaction of all nine cellular structures to a uni-axial compressive strain as high as 75%.Nominal stress-strain data obtained from simulations are compared with the experiments in figure 9.The stress-strain plot for each structure is accompanied by the corresponding snapshot (front view) from the simulation at 45% strain.The dark solid lines in the stress-strain plots are from FE simulations.The redshaded region represents the total span of measured stress-strain data obtained from multiple uniaxial compression experiments.The experimental data set considered here shows bounds in the stress-strain data obtained from testing following protocols discussed in the earlier section.The observed variability in the stress-strain data was attributed primarily to inconsistency in printing and to some minor extent material irregularity and measurement error.

Deformation of the cellular structures during compression.
Based on the compression of nine cellular structures studied in this work, we observed three distinct mechanical responses based on buckling response that is consistent in both experiments and simulations.We identified the following modes of deformation: (1) global structural buckling, (2) local cell buckling, and (3) bending-dominated collapse (no buckling).The schematic of these behaviors is shown in figure 10.
The term global structural buckling here refers to the non-uniform, macroscopic buckling of the entire specimen at once attributed to a few localized deformation bands.Similar to the buckling of the column under compression, in this case, the entire length of the specimen is affected by such buckling phenomenon and subsequently the specimen length would be used as the characteristic length to formulate a buckling load.
Here local cell buckling refers to the case wherein the structure as a whole does not buckle at once; however, there is local buckling of individual unit cells.From simulations, the first instance of buckling was observed to initiate at a layer of cell walls with the weakest strength (likely due to the variation in printed dimensions).Once a layer buckles, it was observed to influence the stability of unit-cell layers solely in its immediate vicinity, thus triggering a wave of unit-cell buckling or collapse that then propagates and traverses the gauge length of the structure leading to the layer-by-layer gradual collapse of the specimen, subsequently leading into the densification.
With regards to the third case of bending-dominated collapse, neither the cells-walls nor the overall structure experience any form of buckling as in the first two cases.The deformation is a result of non-buckling modes of deformation, e.g.bending.Each of these cases is discussed in detail henceforth.In the following subsections, we present further results from simulations of representative cellular structures corresponding to the individual cases of buckling.

3.3.1.1.
Case I: global structural buckling.SC, Kelvin, and GA structures demonstrate global structural buckling.Such macroscopic buckling (i.e.buckling with wavelengths that are of the order of the length of the specimen) has been observed in the square lattice structures [59] and some 2D periodic soft porous structures [60].The combined effect of strut thickness and the specific topology of the unit cell results in these structures demonstrating global buckling.In Figure 10.Schematic showing modes of collapse (buckling and/or bending) seen in the simulated cellular structures.1. Global structural buckling: the entire structure undergoes macroscopic buckling.2. Local cell buckling: every unit cell undergoes buckling, however, the structure as a whole does not buckle.3. Bending dominated collapse/no buckling: the overall deformation is characterized by the bending of the struts at joints within the individual cell, and there is no apparent buckling instability.
these structures, the strength of the load-supporting members that would otherwise activate the local cell buckling is typically higher than the stress that activates the global buckling, thus the structural level buckling occurs before the local buckling [59].In SC, Kelvin, and GA, the non-uniform global buckling was seen to initiate with the formation of a band that spanned over a few layers of the unit cell.The band demonstrated a highly localized gradient in the strain.Meanwhile, the remaining cells that were outside the band collectively assisted in the overall buckling of the specimen.Consequently, different cells in the specimen exhibited different deformation and resulting shapes contingent upon their respective positions within the structure.
To further illustrate the progression of structural level buckling, the compression of the SC structure is shown as a representative in figure 11.The figure shows the global structural buckling that was activated at approximately 9% strain.The band composed of a few layers of cells is indicated by the highlighted rectangle.This instability corresponds with the onset of a plateau demonstrated by the stress-strain curve.With the progression of non-uniform structure-level buckling, the stress-strain plateau persisted until approximately 60% strain, at which point the cell struts began to close in and impinge upon each other, thereby leading to the densification regime.
Regarding figures 9(a)-(c), for the Kelvin structure, the structural buckling was activated at 23% strain whereas for the GA structure, buckling was activated at a compressive strain of 4%.This relative difference in the activation strain for buckling indicated higher strength to activation of structural buckling.A lower strain to activation of structural buckling was associated with corresponding lower stress for the onset of plateau, e.g. the SC structure has a plateau onset stress and strain that are lower than those of the Kelvin structure.Furthermore, a lower onset strain leads to a larger strain span over which the plateau sustains, as is the case for the SC structure.Interestingly, the GA structure has the lowest plateau onset strain and stress compared to both SC and Kelvin structures.For the GA structure, the plateau stress is mildly increasing with strain, and the strain span over which the plateau sustains is also lower, compared to the SC and Kelvin structures.Therefore, the mode of buckling, the stress-strain behavior, and the onset of plateau stress are all dependent on the topology and can be modulated by altering the shape, size, and structure of the unit cell.However, such structures are prone to macroscopic-level buckling and might not be suitable for structural stability applications.

Case II: local cell buckling.
The distinct feature of this mode of deformation is the buckling was restricted over the length of the unit cell, hence the term local buckling.Subsequently, the unit-cell length would be used as the characteristic length to formulate a buckling load.Local buckling of the unit cells appeared to initiate in the specific regions (typically near the top and bottom platens) and traversed the structural dimensions thereon, accomplishing the collapse of the specimen.The length of the specimen did not seem to influence the stress for the onset of buckling.The prominent feature in this process of buckling-induced collapse was a layer of unit cells in the form of a plane, which spanned specimen dimensions in the transverse directions to the loading, with the normal to the plane aligned along the compression loading.Several instants of such a layer of unit-cell buckling were progressively activated.
We observed this buckling mechanism for the gyroid, Schoen IWP, primitive, and octet structures under compression.For example, figure 12 shows snapshots from a compression experiment and corresponding FE simulation of the gyroid structure.Snapshots are presented at different strains, as indicated by markers on the stress-strain curve.At low strains, the Gyroid structure demonstrated uniform deformation accompanied by linearity in the stressstrain response.The onset of buckling was identified at a strain of approximately 18%.Deformation of the specimen soon after 18% strain was non-affine, a contribution from the post-buckling non-linear response of the buckled row/s of cells and linear response of the intact cells.Subsequent increases in the strain progressively activated the buckling of other rows of cells leading to the gradual collapse of a fractional volume that extended to the remainder of the specimen.The progressive buckling of cells corresponded with the plateau in the stress-strain curve that also appeared to last until the repetitious mechanism of individual cell-buckling and collapse traversed the entire specimen.A previous study by Higuera et al [61] indicates that gyroid structure does not necessarily buckle layer-by-layer.For higher densities above 25%, the collapse occurred for all layers simultaneously.Further increasing the strain resulted in compressing the unit cell in the post-buckling regime, resulting in the gradual densification of the specimen and corresponding with the onset of the densification regime.
Figures 9(d)-(f) shows the comparison between stress-strain responses of gyroid, Schoen IWP, and primitive TPMS structures, all of which demonstrate this mode of buckling.The primitive structure buckled at a lower stress (plateau at 0.015 MPa) than the gyroid (plateau at 0.023 MPa) and the Schoen IWP (plateau at 0.025 MPa) as shown in figures 9(d)-(f).The primitive structure thus demonstrated lower compressive strength than the gyroid and Schoen IWP, consistent with observations made by Abueidda et al [62].Furthermore, Schoen IWP, gyroid, and primitive in that order seem to indicate qualitatively higher transverse strain for the same compressive strain, as seen in figures 9(d)-(f).We leave this aspect of research on the influence of unit-cell topology on Poisson's ratio for future study.
Despite the lattice structure, the octet structure demonstrated the Case II mode of buckling.The relative orientations of lattices were found to be favorable towards mitigating the sudden modes of instability, e.g.structural buckling seen in the SC.The Octet structure collapsed at a lower stress with a plateau at 0.011 MPa relative to any of the TPMS structures.The lower stress may be attributed to its thin struts, 343 µm in thickness, joined together in an intricate fashion with lower strength and buckling resistance.
Yet another interesting aspect of the Case II buckling mode is that every unit cell ultimately undergoes buckling in the collapse regime of the stress-strain, unlike the Case I mode of buckling in SC, Kelvin, and GA structures.Thus, the non-linear compression behavior of Case II structures is likely to be more influenced by the mechanical properties of the individual cell topology.The structures demonstrate Case II buckling because of their ability to undergo buckling at the unit cell, which appears to have a lower critical load for buckling than that for the structure.Moreover, the buckling-dominated plateau regime of the stress-strain is much more regularized.This characteristic can be ideal for constructing cellular materials requiring structural stability.

Case III: no buckling, only bending-dominated collapse.
In this category, the deformation is dominated by bending, and there is no buckling at the cellular or structural levels.Such compression behavior with no apparent instability was demonstrated by BC and diamond structures.For example, the BC cellular structure compaction is shown in figure 13.The sub-structures of the unit cell, which in this case are made of struts, are arranged such that the struts can accommodate global deformation simply by changes in the relative rotation angle between adjacent struts.This resulted in only the strut bending at the joints and without a requirement for contribution from buckling.The figure is accompanied by the stressstrain curve for the compaction with markers indicating strain for corresponding snapshots in figure 13(b).The figure shows an excellent quantitative agreement on the stress-strain and qualitatively on the extent of dilation between the simulation and experimental results.Due to the lack of a buckling mechanism in this case, the stress-strain curve demonstrated only two distinct regimes (linear deformation and densification) without the stress plateau.For BC, the transition from the bending-dominated deformation to the densification regime occurs at around 50% strain.
For all the aforementioned buckling cases, the evolution of the buckling and bending mechanisms at the scale of a unit cell is shown in figure 14.Regarding the SC structure shown in figure 14(a), which demonstrated Case I buckling, the extent of buckling-induced deformation of individual unit cells within the specimen was non-uniform, i.e. was dependent on the location of the unit cell within the structure.Cells in the compaction band experienced a relatively larger extent of deformation compared to the cells that were further away from the compaction band.We included a snapshot of a unit cell in this localized deformation band.For this specific unit cell, the direction of loading was along the axial direction of the vertical strut.This resulted in the progressive buckling of the vertical strut under compression.The horizontal strut does not appear to undergo a significant deformation until the densification of the specimen.
Regarding the gyroid structure shown in figure 14(b), which demonstrated Case II buckling, all cells buckled and collapsed similarly.Since TPMS unit-cell structures are made of planar cell walls rather than ligaments or beams, multiple planar modes of deformation likely coexist, e.g.wrinkling and folding.We used buckling as an overarching term for all the modes of deformation experienced by this structure.
Regarding the BC structure shown in figure 14(c), which demonstrated Case III bending mode of deformation, the struts rotated relative to each other; however, none of these cells demonstrated any buckling.The cell deformed via the relatively compliant rotational degree of freedom at the joint between ligaments.This mode of unit-cell deformation has been documented elsewhere [63,64].The low compression stiffness in the linear regime of the stress-strain curve can be traced to the rotational stiffness of these joints at the unit-cell level.Compared to the SC and Gyroid structures, the BC structure demonstrates a relatively compliant behavior at the macroscale with only two distinct regimes in its stress-strain response.The lack of buckling instability at the scale of the unit cell for the BC structure results in a stress-strain response without a stress plateau.

Energy absorption density of the cellular structures.
An energy-absorbing material, e.g.within a packaged layered composite, typically dissipates a portion of the total incoming kinetic energy during an impact while protecting sensitive or critical structures from experiencing excessive stress [65].The cellular structures constructed from polymeric materials can commonly undergo large deformation and absorb energy through cell bucking and bending.Such possible non-linear modes of deformation modulate the incoming impact stress, thereby limiting the transmitted stress that can be identified by the long plateau in the stress-strain curve.Due to the ability of compliant polymers to handle large deformation, these materials typically do not exhibit apparent fracture or damage.The energy absorption density, W(ε), which is the area under the stress-strain curve until a strain ε, and the efficiency parameter, η(ε), which is the ratio of the energy density and the peak stress, (σ p ) are often used as the indicators in characterizing cellular structures under compression and impact [41,65,66].We further elaborate on these expressions and use them to analyze the energy absorption capacity of cellular structures in this study: W(ε) and η(ε) are plotted for different levels of strains and stresses obtained from FE simulations, shown in figure 15.Such plots are commonly useful in identifying a particular cellular structure better suited for known stress and strain thresholds per application requirements.TPMS structures demonstrated higher energy absorption density compared to strutbased structures at moderate to high strain levels as shown in figure 15(a).Simulations on 12 × 12 × 12 sized TPMS structures were conducted to consistently compare energy absorption and efficiency measures.Schoen IWP 12 × 12 × 12 has the highest energy absorbing density.This was followed by the gyroid 12 × 12 × 12, Schoen IWP 6 × 6 × 6, gyroid 6 × 6 × 6, and primitive 12 × 12 × 12.The high energy absorption density of TPMS structures is likely associated with the fact that all cells undergo buckling and the onset of buckling was at a relatively higher plateau stress compared to other structures.Interestingly, the Kelvin structure demonstrated energy density greater than Primitive 6 × 6 × 6 and nearly as high as Primitive 12 × 12 × 12 even though it is a strut-based structure that underwent macroscopic structural level buckling.Among all the structures, BC has the lowest energy absorption density at strains before its densification.
The W(ε)-stress and η(ε)-stress plots are often used in identifying suitable cellular materials [67,68].Identifying whether the cellular structure can absorb a particular amount of energy within a particular stress or strain limits is often desirable to safeguard the protected object [68].Figure 15 shows the plot of these characteristic measurements versus stress in (c) and (d).At low stresses, BC and GA have higher energy densities and efficiency measures compared to other structures.IWP and gyroid structures have much higher energy absorption density but also demand elevated peak stresses during compression.Therefore, even with high energy densities, IWP and gyroid have similar peak efficiency parameters (IWP 0.43 and Gyroid 0.38) compared to SC (0.39), Octet (0.40), and Kelvin (0.41).
For identifying materials with a required energy absorption characteristic, the energy absorption density is typically assessed up to the strain when densification begins (i.e. the onset of densification strain [66], ε d ).One of the different methods to determine the onset of densification strain is from efficiency vs. strain plot (figure 15(b)).It is the strain at the last peak of the efficiency curve [41] where dη(ε) dε = 0.The onset of densification strains (ε d ) for different structures is provided in supplementary table S3.The energy absorption densities and the peak stresses until the onset of densification strains for different cellular structures are plotted in figure 16.Schoen IWP 12 × 12 × 12 structure has the highest energy absorption density followed by gyroid 12 × 12 × 12.These structures have energy densities more than twice their 6 × 6 × 6 counterparts.However, the peak stresses in these structures are also significantly higher (Schoen IWP 12 × 12 × 12: 0.059 Mpa and gyroid 12 × 12 × 12: 0.057 MPa).This suggests the possibility of modulating the energy absorption capacity and peak stress by altering the unit cell size.Among strut-based structures, diamond demonstrated the highest energy absorption density when accessed until the onset of densification even though it only underwent bending and no buckling.This is because of the monotonic increase in the stress with applied strain and the relatively high onset of densification strain (64.6%).BC has the lowest energy absorption density among all the structures in this study because of the low densification strain (50.8%) and compliant linear deformation regime.
This study accomplished AM printing of structures using Elastic 50A Resin, a polyurethane elastomer.Elastic 50A has attributes including high elongation, energy restitution, and tear resistance, is transparent, and has a pliable texture reminiscent of molded silicone parts.We optimized post-processing conditions involving washing and curing to ensure that there was not any leftover resin that would otherwise affect the mechanical performance.
This work further developed simulations for polymer-based cellular materials under large compression.Simulations of large compression through densification were accomplished by overcoming commonly encountered FE mesh distortion and degeneration errors by deliberating with the resulting Jacobian for elements used in mesh discretization.Mesh discretization for simulations involved using a narrow distribution of FE mesh sizes that sufficiently resolved the cell wall thickness to accurately capture the strain at which buckling instability was seen in the experiments.A mesh discretization involving a minimum of four elements across the cell wall combined with a linear shape function was found suitable for simulating all modes of buckling.In addition, a hyperelastic material model calibrated from tension and compression data on the elastomer was sufficient to accurately describe the pristine elastomer's mechanical response undergoing large deformation.
We accomplished insights into multiple modes of buckling and post-buckling response that resulted in specimen densification.SC, Kelvin, and GA structures demonstrated a structural mode of buckling.The structural length of the specimen had a part in the process of buckling.Essentially, a buckling mode that is typically seen for long slender columns was observed for these structures under compression, despite the lowest aspect ratio of the specimens.Consequently, the specimen length was suggested as the characteristic length for the onset of buckling for these structures.On the contrary, the Case II, the unit-cell mode of buckling was restricted to the length of a unit-cell, without the structural length directly partaking in the process of buckling.Lattice-based structures BC and diamond lacked either structural or unit-cell buckling modes.Consequently, the stress-strain curve of these structures exhibited linear deformation and densification without the presence of a stress plateau.
This study is not without limitations.The calibration of the constitutive models in the simulation was based on uniaxial tensile and compressive loading.For future studies, incorporating experimental data from biaxial, planar, and volumetric tests would give a more robust representation of the constitutive model for the pristine material.Simulations did not study the effect of varying VF, size of the unit-cell, the specimen size effect, and printing-induced variability in properties of the pristine material.The specimen dimensions considered in the three orthogonal directions were equal.If the specimens were to be slender, i.e. one dimension much larger than the other two dimensions, this may demonstrate structural buckling mode in the TPMS (Case II) and lattice (Case III) structures, as well.Furthermore, only uniform unit cell designs were used for all cellular structures.The work by Dai et al [69] further provides insights on incorporating gradients in the gyroid lattice structure.From their results, the gradient gyroid lattice exhibits a relatively improved energy absorption value than the uniform gyroid Finally, the integrated modeling and experimental approach in this study facilitated insights into different modes of buckling influenced by various unit-cell topologies.Mechanistic insights from this study can assist material engineers and scientists in identifying suitable polymer-based cellular structures for applications requiring energy absorption properties.

Conclusions
This study identified post-processing protocols for Elastic 50A Resin, a polyurethane elastomer.SLA was used to obtain cubes made of the same material VF but with different unit cell mesoscale structures, resulting in widely varying mechanical responses.The hysteresis response of nine cellular structures was observed through four compression cycles, and the shape of the fourth hysteresis loop was in qualitative agreement with the one obtained from the third cycle, thereby precluding an irreversible inelastic response from the stress-strain.
Mesoscale simulations that resolved the sub-structure associated with the unit cell were able to capture various buckling modes as observed in the experiments.Buckling modes were found to be influenced by the topology of the unit cell.SC, Kelvin, and GA structures demonstrated a structural mode of buckling.Octet and TPMS structures, namely Gyroid, IWP, and Primitive demonstrated unit-cell mode of buckling.Lattice-based structures BC and diamond lacked either structural or unit-cell buckling modes.Consequently, the stress-strain curve of BC and diamond structures exhibited linear deformation and densification without a stress plateau regime.
Validation of individual mesoscale simulations was accomplished for multiple AM cellular topologies by comparing simulated responses against experimental measurements of the stress-strain and corresponding buckling modes from the same specimen geometry as used in the simulation.Comparison of simulated stress-strain and buckling modes are in excellent agreement with the corresponding experimental data.Simulated results were then used to investigate the mechanical performance of the cellular materials.We compared energyabsorbing efficiency, subsequently leading us to identify the IWP structure as the highest energy density among all the cellular structures considered in this study.

Figure 1 .
Figure 1.Unit cell isosurfaces for TPMS structures (a) Gyroid, (b)Schoen IWP, and (c) Schwarz primitive.Each nodal position on individual surfaces satisfies the nodal approximation equation (Φ = 0) given by equations (1), (2) and (3), respectively.The length of the unit cell L is the length of the bounding box indicated in the figures.

Figure 2 .
Figure 2. 3D CAD models of cellular structures with part size 30 mm × 30 mm × 30 mm and volume fraction (VF) 20%.The insets show the magnified images of the corresponding unit cell for individual topologies.

Figure 5 .
Figure 5. (a) Calibration of the hyperelastic constitutive model from uniaxial tensile and compression experiments.(b) Calibration of the viscoelastic two-term Prony series from stress relaxation testing.

Figure 6 .
Figure 6.Comparison of cell-wall dimensions between the intended geometry from STL (blue bars) and the corresponding printed structure (gray bars).Average dimensions are reported for the cell walls of the printed structures.The error bars indicate the total span for all the measured values.

Figure 7 .
Figure 7. Representative optical images at 20x magnification showing uniformity quality of additively manufactured cellular structures after post-processing.

Figure 8 .
Figure 8. Experimental characterization of loading and unloading stress-strain curves over four cycles of compressive for all nine cellular structures.

Figure 9 .
Figure 9. Finite element simulations of nine cellular structures under uniaxial compression.Stress-strain response and front-view snapshot at 45% nominal strain are presented for individual structures from FE simulation.

Figure 11 .
Figure 11.Compression of the simple cubic (SC) cellular structure.Qualitative comparison between (a) experiment and (b) simulation snapshots at different strains.The simulated results are colored using the von Mises stress value.The corresponding simulated stress-strain result and the experimental data are plotted in (c).

Figure 12 .
Figure 12.Compression in Gyroid cellular structure.Qualitative comparison between (a) experiment and (b) simulation snapshots at different strains.The simulated results are colored using the von Mises stress value.The corresponding simulated stress-strain result and the experimental data are plotted in (c).

Figure 13 .
Figure 13.Compression in body-centered (BC) cellular structure.Qualitative comparison between (a) experiment and (b) simulation snapshots at different strains.The simulated results are colored using the von Mises stress value.The corresponding simulated stress-strain result and the experimental data are plotted in (c).

Figure 14 .
Figure 14.Evolution of buckling and bending mechanisms with strain for individual unit cells of (a) SC, (b) gyroid, and (c) BC cellular structures subjected to compression.The location of a unit cell within the corresponding cellular structure is highlighted using a red square in the leftmost column.

Figure 15 .
Figure 15.Energy density and efficiency plotted as a function of strain (a), (b) and peak stress (c), (d) for various cellular structures.Comparison is also made between two different unit-cell sizes for Gyroid, IWP, and Primitive.

Figure 16 .
Figure 16.Comparison of energy absorption density (W) and peak stress of simulated cellular structures.The energy density was calculated as the area under the uniaxial stress-strain curve over the range of strains from zero through the onset of densification (ε d ).