Tunable 3D printed composite metamaterials with negative stiffness

The paper proposes a class of tunable metamaterials that use inclined beams to achieve instability in a rigid system. Three different beam tilt angles, 25°, 45°, and 60°, are evaluated in the form of unit cells using quasi-static compression tests and numerical simulations. Snap-through behavirous are characterised by structural stiffness and buckling load. Periodic and gradient structures are assembled and analysed by arranging the unit cells in rows and columns. Size effect analyses and parametric studies are carried out on various unit-cell arrangements and different beam angles. The proposed metamaterials are manufactured through fused filament fabrication 3D printing technology with a composite material, onyx. The results from experiments, finite element analysis, and analytical models are compared and evaluated. The structural stiffness and buckling load are shown to be positively related to the inclination angle of the tilted beams. The number of rows of unit cells governs the nonlinear mechanical response (number of snap-throughs) of multiple-layered structures. By increasing the number of rows and columns of unit cells, which are less prone to manufacturing defects, the reliability and repeatability of the structural properties of periodic/gradient structures could be improved. A design plot is also provided to predict and tune the snap-through behaviour of multiple-layered structures via beam angles and unit-cell arrangements.


Introduction
Mechanical metamaterials exhibit exotic or extraordinary properties that cannot be achieved in conventional materials [1].These mechanical properties are mainly decided by their intrinsic architecture design [2], surpassing those of their constituent materials.Metamaterials, that features zero/negative poisson's ratio [3,4], lightweight and high stiffness, * Author to whom any correspondence should be addressed.
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.high energy absorption (EA) [5], and shape-morphing [6] have already been realised.Recently, structural instability has been harnessed in metamaterial design, inspired by nature that often reacts by collapsing to internal and external effects.For example, snap-through behaviour (a type of instability) leads to the leaf closure of Venus flytrap [7].Thus, the ability to use instability in metamaterial design to attain novel properties has become a trending topic.
Snap-through is an instability that induces a transfer of a body from one configuration to another when load exceeds a critical value [8].Elastic beam is a well-embraced design element to achieve snap-through.It allows for large deformation and is liable to buckling, which equip metamaterials with exceptional properties, such as negative stiffness, enhanced EA, multi-stability, damping, and programmability.Negative stiffness is normally exhibited as a snap-through instability within an elastic system [9], requiring a reduction in applied force to generate an increase in displacement [10].Theoretically, elastic beams undergo a buckling instability when compressed axially and retain their initial configurations when the load is removed.This instability provides nonlinear but reversible design components for metamaterials [1].Correa et al [10,11] investigated the compressive behaviour of negative stiffness beams that were configurated in a honeycomb shape.By comparing these novel structures with conventional honeycomb structures, they found that negative stiffness honeycombs could fully recover after compression (also exhibits higher positive stiffness).Based on the early design, a three-dimensional conformal structure was proposed and demonstrated to be capable of impact protection from multiple loading directions [12].Such reversible and repeatable instability of beams provide a new strategy for enhancing EA.
In addition to reversible structures, studies on creating multi-stable configurations have emerged in recent years.Restrepo et al [13] introduced a novel material by periodically arranging negative stiffness unit cells, which consisted of curved beams and exhibit multiple stable states.Such materials possess excellent EA capacity, characterized by a long and serrated loading and unloading plateaus.Moreover, metamaterials with more than one stable configuration was reported to trap the input energy in the form of elastic deformation of tilted beams [14].The energy-trapping mechanism relied on the reversible change among different stable geometries of the prescribed beams [15].Certain types of applied force were required to reverse the beam structure back to its original shape.The energy-trapping attribute of multi-stable materials could effectively isolate vibration, offering a more costeffective solution to controlling vibration of any machinery or equipment from affecting their surroundings.With more variations, Zhang et al [16] developed a set of multi-stable meta-structures exhibiting both level and inclined stable configurations.The multi-stability was demonstrated to be used in design of morph-changing structures, vibration isolator, damper, and energy absorbers [17][18][19][20][21].
Instead of triggering instability from compression, Rafsanjani et al [2] proposed metamaterials (figure 1(a)) that underwent large deformations and switched between two configurations under tension.The phenomenon was realised through cosine-shaped beams.It was also pointed out that mechanical responses could be tuned by changing the amplitude of beams.Besides using single material to fabricate the buckling beams, studies on multi-material beam designs were also conducted.A structure with sandwich beams was 3Dprinted with rigid material [22], polyamide (sandwich plates), and flexible material, thermoplastic polyurethanes (sandwich core).Figure 1(b) shows the response of sandwich structures under cyclic compression load.This hybrid compliant structure was demonstrated to exhibit negative stiffness and high level of recoverability, which has potential for shock isolation, vibration control, and deployable structure.
Table 1 summarises the materials that were used in previous literatures to fabricate inclined beams.Most studies have employed either soft elastomers or stiff metal.However, the use of extremely soft or stiff materials is not practical.A certain level of stiffness and strength [23], robustness against environment [24], and reliability regardless of manufacturing defects are required in different engineering applications [17][18][19].Additionally, the materials mentioned previously were homogenous and suitable for small-scaled manufacturing, there is a lack of understanding of using compliant metamaterials using inhomogeneous materials and also practical for larger-scaled fabrication.Therefore, a type of carbon-fibre filled polymer is used in this study, which offers high strength, toughness, and chemical resistance etc.To process reinforced polymer, the fused filament fabrication (FFF) technique is utilised.Compared to conventional manufacturing methods, such as injection moulding and casting, FFF is toolless, zero/nearzero waste, and less challenging in iterating unit cells [25].However, the inherent defects in the printed parts from FFF could limit their mechanical performance, such as voids and weak bonds between fibre and matrix [26,27].
Inspired by materials in nature (e.g.bamboos, human bones) that could adapt to the changing environment, gradient metamaterials were demonstrated to obtain a continuousspatial variation of properties [33][34][35][36] by gradually varying their internal structure.For instance, the graded auxetic structures were found to yield higher resistance towards blast loading [37], higher EA [38], and experience gradual failure to sustain larger deformation under compression [39] comparing to the uniform designs.Tailored mechanical properties could be achieved, thus making them suitable for a wide variety of engineering applications [40,41].However, most studies focused on uniform structures [2,[10][11][12][13][14] with negative stiffness, while graded compliant structures have been limitedly explored to date.In 2021, Chen et al [30] proposed a gradient honeycomb structure with negative stiffness, revealing its superior energy control performance over the uniform honeycomb structures.
Here, we introduce a series of metamaterials with the ability of tuning peak load and EA capacity.Unit cells feature straight beams with different tilted angles.These unit cells exhibit snap-through behaviours with different mechanical properties.By arranging unit cells, both uniform and graded layered structures are proposed.Samples are manufactured with carbon-fibre reinforced nylon via 3D printing.Quasistatic compression tests, numerical simulations, and theoretical calculations are used to investigate the snap-through responses.X-ray µCT is conducted to examine the printing quality of unit cells.This study aims to provides insights into buckle-induced instability within rigid structures and to achieve tunable snap-through behaviours for future designs of protective structures, sensors, and actuators.

Specimen design and fabrication
The design of unit cells is shown in figure 2, controlled by three geometric parameters of beam (figure 2(a)).These three parameters capture two features, i.e. slenderness ratio and tilted angle.To investigate the effect of the tilted angle of the beams on the mechanical behaviours, three different unit-cells are proposed with different inclined angles, 25 • , 45 • , and 60 • (figure 2(b)).To control the variables, the same slenderness ratio (0.16) and depth (7.5 mm) are applied to all unit cells.Inclined angles and slenderness ratios were chosen according to the design map generated by Shan et al [14], which ensures occurrence of snap-through behaviour with different mechanical properties.
Samples (figure 2) are fabricated with a composite material, onyx, using MarkTwo printer, (Markforged Inc., USA).In this study, onyx is provided by Markforged, composed of nylon (80 vol%) and short carbon fibres (20 vol%).Carbon fibres are 200 µm in length, and pre-orientated in the printing direction [42].A 0.4 mm diameter extruder is used for filament deposition.Printing parameters are the same for all specimens, including 0.1 mm layer height, single-layer outline, and solid infill with a −45 • / + 45 • pattern.
Structures consisted of multiple unit cells are also designed and manufactured to further evaluate the tilted beams effect.Here, three types of structures are proposed, including twoby-one vertical structures (see figure 2(d)), two-by-two structures (see figure 2(e)), and one three-by-three structure (see figure 2(f)).Two-by-one vertical structures are designed by stacking two same unit cells together in two different ways (oriented and mirrored).Two-by-two structures are simply designed by combining two of the two-by-one structures in horizontal direction.Additionally, a three-by-three gradient structures, which is consisted of all three unit-cells, is proposed to investigate the interaction among different unit-cells.

Quasi-static compression tests
Uniaxial compression tests were conducted on three different unit cells as well as their correspondingly periodic structures.Experiments were performed using Shimadzu universal testing machine (AGS-X series), with a 50 kN load cell.A loading rate of 2 mm min −1 was applied to all tests.While unit cells were compressed to a maximum crosshead displacement of 7 mm, double and triple layered structures were compressed to 14 mm and 21 mm to maintain an equal level of compression for each design.Three specimens were tested for each design to ensure repeatability.Force-displacement curves and real-time videos were captured from the experiment to evaluate the responses of different structures.

Numerical simulation and parametric studies
A finite element (FE) model was developed using Abaqus/Explicit 2020 (Dassault Systems SIMULIA Corp.) to analyse the response of the unit cells and structures under quasi-static compression.Considering that the thickness of the specimens is comparatively smaller compared to their height and width, all specimens were simplified as 2D structures.Besides, symmetric analysis was used throughout the modelling because of the symmetry feature of designed unit cells.The symmetry was taken into account by applying the y-axis symmetry boundary condition of along the symmetrical axis of the unit cell in x-y plane (see coordinates in figure 3(c)).All the specimens were considered isotropic in the x-y plane in spite of the anisotropy of the constituent material.The reasons mainly lie in the highly oriented short fibres and +45 • / − 45 • alternating infill patterns.Current studies have demonstrated that short fibres are highly aligned in the nozzle moving direction for filament fused fabrication with carbon fibre reinforced polymers [42][43][44][45][46].For mark 2 printer, +45 • / − 45 • is pre-set as default infill trajectory and could not be customised by users.Given that the printing path is alternating between the +45 • infill and −45 • infill for reasonably large number of layers (75 layers in total), the in-plane anisotropy of each layer could be sufficiently compensated.Hence, the property of the as-fabricated specimens could be considered isotropic in terms of the orientations of the short fibres.Nevertheless, the inclined beams were the regions of interest comparing to the rest parts of the specimens, which were considered rigid and unreformable due to the geometric design.According to the printing paths preview, the beams were fabricated only with the contour wall with the short carbon fibres highly oriented along their longitudinal direction.Therefore, the beams within the specimens were technically symmetric and isotropic in the x-y plane.
Given that the metallic compression plates of Shimadzu machine are much stiffer than the 3D-printed onyx samples, two rigid plates were modelled for simplification.A normal contact behaviour was defined between the specimen and rigid plates by hard contact formulation, while penalty friction formulation with a friction coefficient of 0.3 was used to describe the tangential behaviours [42].It was assumed that structures were under plane stress condition with a thickness of 7.5 mm in z-direction (see coordinates in figure 3(c)).Two reference points were generated, with rigid body constraints applied between the reference points and rigid plates.Besides, the bottom plate was defined fixed in all directions, while displacement in the negative y-direction was applied to the top plate to implement compression.
The constitutive behaviour of onyx was simplified to be elastic-perfect plastic.Material properties used in the model were obtained from tensile tests on 3D-printed onyx specimens under ASTM D638 standard, which provided a Young's modulus of 1800 MPa and a yield strength of 61 MPa.Density and poisson's ratio are 1.2 g cm −3 and 0.3, respectively, adopted from onyx datasheet provided by Markforged material supplier [47].
Structures were discretised with 4-node bilinear plane stress quadrilateral (CPS4R) elements.To optimise the computational efficiency and accuracy, a mesh size convergence study was conducted firstly.This involved both top and bottom plates and the tilted beams.Four different mesh sizes (0.15, 0.2, 0.3, 0.4 mm) were evaluated for the tilted beams, and four different mesh sizes (0.2, 0.4, 0.6, 0.8 mm) were assessed for the plates.Results (figure 3) show that 0.2 mm and 0.4 mm are the optimal mesh size to capture the response of beams and plates, respectively, to guarantee desired accuracy and reduce the computational cost at the same time.
Firstly, the models of unit cells were validated by the experiment, and then simulations involving periodic structures were performed.Besides, parametric studies were carried out using the validated numerical models.Here, we proposed different three-by-three structures, considering all the combinations of different arrangements and orientations of the unit cells.

Analytical calculations of critical buckling load (Pcr)
To further understand structural behaviours of proposed unit cells, an established analytical model was used under certain assumptions to calculate the theoretical buckling load P cr .
Given the dimensions of tilted beams (5 mm length, 7.5 mm width, 0.8 mm thickness), a wide and thin plate buckling theory is used [48].Critical buckling load, which is defined as the minimum force to buckle a body, is mainly decided by the elastic modulus and cross-sectional (bending) stiffness of plate.In this case, the following assumptions were made to simplify the structural response of tilted plate: • All the tilted plates deform elastically before buckling occurs, • The cross-section of a tilted plate is initially constant along its length and remains rigid in its own plane during deformation, • Expect for the tilted plates, elements of the unit cells experience no deformation or buckling.
According to rees (2009) [48], the critical load that induces tilted plate buckling is expressed as: where E s is elastic modulus of onyx material, I is second moment of area of the inclined plate cross-section, and l e is the equivalent length.For a plate with rectangular cross-section, I is calculated as: where b is the width of plate that is tangential to the deformation plane, and t is the thickness of the plate.
Since the equivalent length of plate is highly dependent on the connection type of the plate end and its adjacent structures, three simplified types of joints were considered, including fixed-fixed ends, fixed-pinned ends, and pinned-pinned ends.Here, l e is calculated according to [48]: where ν is poisson's ratio of onyx material.

X-ray computed micro-tomography (µCT) analysis
To examine the 3D printing quality of the printed samples, x-ray µCT was performed on three different unit-cells with Bruker SkyScan 1275 (Bruker SkyScan, Kontich, Belgium).
The purpose of multi-model imaging was to obtain a 3D reconstruction of the internal microstructure for the specimens in a non-destructive way.The scanning process was performed at a voltage of 40 kV and current of 100 µA with a voxel size of 21 µm.Then, NRecon software was used to collect 2D image projections and reconstruct into a 3D representative of examined specimen.Details of all cross-sections (2D-slicing) were analysed and visualised through CTVox software.Here, only qualitative analysis was conducted, while porosity and pore information were not reported.

FE model validation and analytical calculation of unit cells
Experimental results for unit cells under uniaxial compression are presented in figure 4(a) and compared to the numerical simulations for the purpose of validating the FE models.The force-displacement curves yield similar trends and features between experiments and numerical simulations, including the elastic stiffness, peak force, and negative stiffness.A good agreement is achieved until the unit cells experience peak forces, while small discrepancies appear thereafter.The results indicate that all unit-cells experienced instability at the early stage of the compression tests.Elastic deformation could be observed at the initial phase for all cases.Once the vertical loads reach a certain magnitude, the beams start to buckle.This leads to reduction in the force, which is a negative stiffness phase, due to stress redistribution during buckling.The force continues to decline until the tilted beams contact themselves or collide with the rest of the structures.As densification initiates, the force increases again with the compression load.force-displacement curves.The following equations were used to calculate the effective axial stress (σ * yx ) and effective axial strain (ε * yx ) in y direction (compression direction): where R y is the total reaction force of unit cells and A y is the cross-sectional area of the whole unit-cell that is perpendicular to y-direction.L y Denotes the height of unit cells (length of unit cells along y-direction), while δ y is the total displacement experienced by unit cells along y-direction.
As indicated in figure 4(b), unit cell 1 yields the highest structural stiffness and critical force.Unit cell 2 exhibits the second highest structural stiffness and experiences a lower peak force than that by unit cell 1.In other words, structural stiffness and peak load are positively related to tilt angle (α) of the beam.To quantitively compare the experimental, numerical, and analytical results, table 2 lists the values for peak force P cr for the three unit-cells.
For unit cell 1, P cr captured from the experiment and simulation lies in between the theoretical predictions from the fixed-end and fixed-pinned end assumptions.For unit cell 2 and unit cell 3, fixed-pinned ends and pinned ends assumptions could predict the critical buckling load well, respectively.From analytical point of view, the different boundary conditions for different unit cells could be explained using the force decomposition in figure 5(a).When the inclined beam is subjected to a vertical force (F), F could be decomposed to an axial compression force along the beam (F 1 ) and a shear force perpendicular to the beam (F 2 ).F 1 contributes to the overall compression of the beam along its length, while F 2 causes bending and rotation due to the induced bending moment and exerts a reaction force on the support of the beam end.According to length relationship for a right triangle, the shear force could be calculated as: With the increase of the inclined angle α, F 2 decreases.Hence, the reaction force of the beam support increases from unit cell 1 to unit cell 2 and unit cell 3 (α decreases from 60 • to 45 • and 25 • ).Assuming the same printing quality for all three unit-cells, the inclined beams in unit cell 1 would experience the smallest rotation comparing to unit cell 2 and unit cell 3. Therefore, the corresponding boundary conditions from the analytical model were demonstrated to predict different unit cells.
To investigate the influence of beam designs on the boundary conditions of beam ends, nodes that are close to both ends of the beams were tracked in numerical simulations to observe the movement of the beam ends.As shown in figures 5(b) and (c), all ends undergo small rotations within the elastic regions.Unit cell 3 experiences the most rotation angles for both the Note: all results are based on the actual dimensions of printed specimens.bottom and top ends, followed by unit cell 2 and unit cell 1.This explains the pined-end boundary condition of unit cell 3. Additionally, there is no rotation for bottom end of unit cell 1 before 0.001 s, indicating one-end rigid connection type.However, discrepancies still exist between the experiment and analytical calculations based on corresponding assumptions.This could be attributed to the variation in cross-section along the beam, which results from sample fabrication and leads to inaccurate prediction from analytical model.Additionally, the 3D-printed joints do not satisfy the theoretical boundary conditions in that they are neither perfectly fixed nor perfectly pinned.

Experimental and numerical results on multiple-layered structures
Responses of multiple-layered structures, which are formed by more than one unit-cell, are investigated through experiment and simulations.Two-layered structures and three-layered structures are considered in this work.

Responses of two-by-one and two-by-two structures.
Experimental and numerical results are compared among different two-layer vertical structures, including two-by-one structures and two-by-two.Figures 6(a) and (b) exhibits the force-displacement curves for the two-by-one structures (oriented and mirrored).Similar mechanical properties were observed between the two-by-one structures and unit cells.Specifically, both highest structural stiffness and maximum force are exhibited in the two-by-one structures inspired from unit cell 1.With a decrease in the tilted beam angle, structural stiffness and peak load of structures experience a decreasing trend.
Overall, simulations and experimental results are found to follow the same trends but are not completely aligned.Simulation results indicate two local peak forces for all the designs.Experimental results show similar features except for the oriented pairs of unit cell 1.The absence of the second peak in the experiments could be attributed to the sliding between the specimens and compression plate, which could be observed in the recording of testing and was not noticed in the simulations by setting up proper boundary conditions.To compare the responses of different structures, configurations of specimens during beam buckling are also shown in figures 6(a) and (b).It can be noticed that deformation/buckling of the unit cell 1 two-by-one structures happens asymmetrically, while the tilted beams of the vertical structures from unit cell 2 or unit cell 3 behave almost symmetrically along the y-axis.The asymmetric buckling is attributed to the defects of 3D printing, which is explored later in this study.
Structures consisted of unit cell 2 and unit cell 3 achieve two peak loads as predicted by the numerical models.However, the responses of structures from unit cell 2 are closer to the simulation results compared to that of unit cell 3.It is worth noting that the structures generated from different unitcells yield two peak loads at different displacements (1 mm and 5 mm for structures inspired from unit cell 1; 1 mm and 8 mm for structures inspired from unit cell 2; 1.5 mm and 6.5 mm for structures inspired from unit cell 3).The displacement gaps between the two peak loads are roughly 4 mm, 7 mm, and 5 mm, respectively.It is also interesting to note that the second densification starts at a displacement of 8 mm approximately for both the structures inspired from unit cell 1 and unit cell 3.However, the structures generated based on unit cell 2 initiate densification later, around 11 mm to 12 mm.Discrepancies appear between the numerical and experimental results with respect to the magnitudes of reaction forces and elastic stiffness.The critical experimental force thresholds for all the designs are smaller than the values predicted by the FE model.Moreover, the specimens experience negative stiffness over a smaller span of forces and displacements in the experiments.These variations could be elucidated from several aspects.First, the viscoelastic behaviour of nylon (main composition of onyx) was not considered in the numerical models.Accordingly, the stress relaxation of nylon could not be captured [11].Moreover, the FE model does not take into account of fracture, which may have happened due to the stress concentration at the intersections of the angled beams and the top/bottom plates.
Similar to the two-by-one vertical structures, the two-bytwo structures respond to compression with two snap-throughs (phase AB and phase CD in figure 6(c)).As two-by-two structures are designed by simply combining two of the two-byone structures in x-direction, the peak loads experienced by the two-by-two structures are twice that of the two-by-one structures.
Figure 6(c) presents the force-displacement curves obtained from the experiments and simulations, with stress distribution corresponding to two snap-through phases (taking the structures generated from unit cell 2 as an example).As the compression load increases, the beams on the top row buckle first.With snap-through happening, the stress redistributes within the whole structure and the force starts decreasing (phase AB).Once the bended beams contact the middle plate, the same instability phenomenon happens to the bottom ones (phase CD).It should be noted here that the mirrored designs yielded similar snap-through behaviours as the oriented ones.Ideally, if the structure is perfectly symmetric along both the x-axis and y-axis, both top and bottom beams would buckle at the same time.Figure 7(a) depicts the force-displacement curve obtained from the FE model using the mesh that is symmetric along both x-axis and y-axis.This symmetric mesh was created by meshing only a unit cell (a quarter of the mirrored structure), using the mirroring transformation to produce the mesh for the other three unit-cells, and stitching them together.Accordingly, this symmetric mesh corresponds to a specimen that is assembled by gluing four identical printed unit-cells together (printing path within the structures shown in figure 7(c)).However, in reality, the specimen was printed following the printing path in figure 7(d), which is not symmetric along neither axis.Therefore, the symmetric mesh could not represent the as-fabricated specimens accurately.As shown in figure 7(a), the force-displacement curve from the symmetric mesh is very different from what recorded in the experiment.In experiment, two snap-throughs and local peak forces were captured, indicating two buckling events sequentially.However, the simulation result only shows one snap-through and it corresponds to all beams buckling at the same time.Furthermore, the as-fabricated specimens could not be perfectly symmetric also due to the printing defects from FFF technology.As the beams were slender, they tended to be sensitive towards defect-induced asymmetry.In order to take asymmetry into consideration in FE models, the mirrored parts were intentionally not meshed symmetrically.Instead of buckling simultaneously, the top and bottom beams of the mirrored structure experienced snap-through at different displacements (figure 7(b)).The similar behaviours between oriented and mirrored structures demonstrated the effectiveness of the slender beam designs to achieve snap-through instability.The slenderness of the beam enabled the structure to respond sensitive to the asymmetry induced by the manufacturing defects or misaligned loading conditions, thereby achieving two snap-throughs from top beams and bottom beams.
Two peak loads are achieved at the displacements of 1 mm and 8 mm approximately, consistent with two-by-one structures.Theoretically, two peak forces predicted was 706 N according to table 2 (unit cell 2).However, the critical load P cr obtained from the experiments is slightly smaller than the analytical prediction.This could be explained by the partially inelastic deformation of beams and manufacturing defects.In addition, the second peak force is slightly lower than the first one.This is resulted from the small elastic deformation of the bottom beams before phase CD, which can be seen from the stress distribution in figure 6(c) (green spots on the bottom beams in phase AB).It should also be noted that the experimental force-displacement curves share similar features as those of the numerical simulations, yet there exist some discrepancies among the three testing results.This indicates a certain level of 3D printing quality issue that can lead to imperfections and defects.Even though all the specimens for each design were manufactured from the same batch, the printing qualities differed among them.
As all the two-layer structures are designed symmetrically, they are expected to deform/buckle symmetrically regarding the compression loading direction.However, the responses appear to be asymmetrical to a certain level, which also lead to the discrepancies between the experimental and numerical results.Given that the compression plates from Shidmazu Universal testing machine are perfectly parallel, the quality of the 3D printed samples could be a predominant reason.Here, we examine the 3D-printing quality of the unit cells and assume vertical structures attain similar features to certain degree.
Firstly, the thicknesses of the tilted beams are measured to compare with its design value (table 3) using a calliper.As we can see from the table, tilted beams as manufactured are thicker than that it was designed and not consistent regardless of the tilt angle.This could be attributed to the limited resolution of the FFF printer, which is 0.4 mm because of the diameter of the filament nozzle.It is also interesting to note that all the beams on the left are thicker than the ones on the right when comparing the thicknesses of the beams within each unit cell.Reason behind this phenomenon is not clear could be the non-uniform processing parameters over the print area during manufacturing [49].Nevertheless, this can be one of the reasons why some specimens are not deformed symmetrically during the experiment.
To further examine the printing quality, micro-CT scanning is used.Figure 8 shows the binary images of the middle crosssection of all three unit-cells, in which black colour refers to pores and white colour represents solid material.As per manufacturing slicing input, all specimens are printed with 100% solid infill in a −45 • / + 45 • pattern (figure 8(a)).Therefore, no pores (black area) inside the specimens are expected if the 3D printing quality is perfect.However, the binary images (figure 8(b)) indicate that there exists a certain level of porosity for all unit cells.Pores can be classified into two major types in terms of their locations, i.e. pores among the −45 • / + 45 • infill, and pores in between walls and infill.
According to the printing toolpaths, pores among infill are only located in the top and bottom plates.As beams are much slenderer compared to the top and bottom plates, pores among infills are not within the scope of study here.The filament gaps in between infill and walls, which construct tilted beams, are more of a concern to the structural behaviour under compressive load.Pore size varies at four different connections between beam ends and plates, especially for unit cell 1.To be specific, gaps on the left joints are observed to be larger than the ones on the right in figure 8(b).This phenomenon explains the asymmetrical deformation of vertical structures inspired by unit cell 1.Furthermore, this observation again explains the discrepancy between the theoretical and experimental load (P cr ) in section 3.1.

Response of three-by-three structures.
As the structures from unit cell 2 exhibit the best printing quality and repeatability, a three-by-three uniform periodic structure (three unit-cell 2 in both columns and rows) is designed and manufactured.Figure 9(a) compares the force-displacement curves obtained from the compression tests and numerical simulations.A certain level of repeatability is maintained in testings, and a good agreement could be found between the experiments and models, which establishes the stable and good manufacturing quality of the samples.Three local peak forces are achieved, and this triples the capacity of a single unit cell 2. Deformation is also shown in figure 9(a), indicating the buckling sequence from the top layer to the bottom layer.
Experiments and simulations were also conducted on threeby-three gradient structures (figure 9(b)).From the experimental and simulation responses, it is observed that the beams in the bottom row (unit cell 3) buckle first, which leads to the first local peak load.Followed by the buckling of the beams in the middle row (unit cell 2), the tilted beams in the top row (unit cell 1) buckle at the end.The forcedisplacement curve obtained from the simulation agrees well with the experimental result prior to the buckling of the middle row.However, the tilted beams inclined at 60 • start to buckle before the middle layer densifies.In addition, notable sliding could be observed between the top compression plate and the surface of the structure.Hence, the last local peak force, which is supposed to be achieved by the undeformed 60 • beams, could not reach the same value as estimated in the simulation.
A notable observation from the responses of both the uniform and graded structures (figures 9(a) and (b)) is that the experimental results is more consistent as compared to the two-layered structures (figures 6(a) and (b)).In other words, the mechanical response of three-by-three structures is less sensitive to the manufacturing imperfections [50].By increasing the number of layers and columns, more constraints in the horizontal and vertical directions are added to unit cells within the structure.Therefore, the structural response becomes more controllable and symmetric.
To investigate the responses of three-by-three graded structures with different arrangements of unit cells, a parametric study is performed using the validated FE model.Figure 10(a) summarises the force-displacement curves for all the structural combinations, regarding the same orientation of the different unit cells.Figure 10(b) shows the response of all the structural combinations with mixing oriented and mirrored unit cells.
All the structures yield the same response under quasi-static compression load, irrespective of the arrangement and orientation of the unit cells.Specifically, the tilted beams unit cell 3 buckle first and lead to the first local peak load.This is followed by the buckling of beams in unit cell 2, whilst the tilted beams in unit cell 1 buckle last.Three local peak forces are obtained from low to high at the displacements of 1 mm, 7 mm, and 14 mm, respectively.As the FE model incorporated a quasi-static compression process with low velocity, all components in the structure experienced the effect of applied force at the same time.According to the experimental and numerical results for unit cells in section 3.1, three different unit cells yield different critical buckling loads.Hence, with the increasing applied load, unit cell 3 (α = 25 • ) buckles first, followed by unit cell 2 (α = 45 • ), and unit cell 1 (α = 60 • ) at the end.
From the results for multiple-layered structures, it could also be inferred that the number of rows of unit cells along the y-direction governs the nonlinear mechanical response, i.e. the number of snap-throughs.

Parametric studies of multiple-layered structures
To extend the present study for achieving various design goals, the size effect of periodic structures was conducted, and this was followed by exploration of more angle (α) values besides those considered (25 • , 45 • , and 60 • ).Here, the peak load (P cr ) and EA were evaluated, corresponding to structural strength and EA capability.
Size effect analysis was conducted on periodic structures, from one-by-one to five-by-five, generated from unit cells 1, 2, and 3. To compare the EA of different periodic structures, quasi-static compression tests were simulated for structures with different periodic size and different unit cells.The corresponding force-displacement curves were obtained, and the areas under the curve provided the energy absorbed by tilted the beams within periodic structures.Figures 11(a)-(c) depicts EA versus compression displacement for the periodic structures of different sizes.For all the calculations to determine EA, the integration of force-displacement curve was performed up to the initiation point of densification.
To better understand the mechanisms underlying the EA variation pattern with the size of structure, specific energy absorption (SEA) and maximum stress are compared among different periodic structures with the same unit cell (figures 11(d)-(f)).Given that the tilted beams are the only deformed elements during compression, the SEA is calculated by considering only the mass of the beams.The convergence of SEA and the maximum stress could be captured from the  two-by-two structures for all three cases.By using the mechanical properties from two-by-two structures, EA and P cr experienced by periodic structures with any size can be predicted.
Based on size effect analysis, several simulations were conducted on three-by-three periodic structures with α ranging from 25 • to 60 • (including 25      , 55 • , 60 • ).Peak load and EA per unit cell were calculated and shown in figure 12(b).
By using this plot, for any multiple-layered structure, as schematically shown in figure 12(a), the mechanical response can be predicted and tailored to a certain level.The multiplelayered structure is consisted of number of N x unit cells in a row and number of N y unit cells in a column.All the unit cells share the same inclined beam angle within each row, while different unit cells could be chosen for different rows.The maximum angle of tilted beam within this structure is α max .For such structure, the number of snap-throughs equals to N y .The maximum peak force experienced by the structure is the product of P cr (for α max ) and N x .The total EA is a sum of the EA of all unit cells.

Conclusions
In this paper, inclined beams were used to make an otherwise rigid system unstable.In the form of unit cells that were put under a nearly static compression load, the responses of different tilt angles were studied.By putting unit cells in horizontal rows and vertical columns and stacking them on top of each other, the peak load and EA of periodic and graded metamaterials were studied.From this, we can learn the following: • For single unit-cells, structural stiffness and buckling load positively correlate to the tilt angle α of tilted beams.On this account, snap-through behaviour could be tuned by varying the tilt angle (α) given a fixed slenderness ratio of beams.• For multiple-layer structures, the number of rows of unit cells govern the nonlinear mechanical response, i.e. the number of snap-throughs.• For uniform-layered structures, displacement gaps between the sequential snap-throughs can be tailored by using different unit cells (4 mm for unit cell 1, 7 mm for unit cell 2, and 5 mm for unit cell 3).• Three-by-three gradient structures yield nearly the same force-displacement curves regardless of the arrangement and orientation of the unit cells.Three local peak forces were obtained from low to high, attributed to three snapthroughs caused by unit cell 3, unit cell 2, and unit cell 1 (in sequence).• Though defects exist in FFF-printed specimens, the reliability and repeatability of structural properties of periodic structures can be improved by increasing the number of rows and columns of unit cells.As more restraints in horizontal and vertical directions are added, structural response appears to be less prone to manufacturing defects.• Snap-through behaviour of the periodic/graded structure can be customised by the different arrangements and different inclination angles of unit cells with.Maximum peak load capacity and EA capability can be predicted according to the design map.

Figure 2 .
Figure 2. (a) Schematic design, showing the control parameter of tilted beams; (b) design of different unit cells, with same slenderness ratio and different tilted angles; (c) pictures of 3D printed specimens (from left to right, corresponding to unit cell 1-3).Schematic 2D drawings and corresponding 3D-printed specimens for: (d) two types of two-by-one vertical structures, oriented and mirrored (taking unit cell 1 as examples); (e) two types of two-by-two structures, oriented and mirrored (taking unit cell 1 as examples); (f) schematic design of a three-by-three gradient structure.

Figure 3 .
Figure 3. (a) Mesh size convergence studies for tilted beams; (b) mesh size convergence studies for top and bottom plates; (c) example of a meshed 2D unit cell 3 (α = 25 • ), showing mesh details for plates and beams.In magnified meshes, the rhombus mesh of the beam has a side length of 0.2 mm, and the square mesh of the bottom plate has a side length of 0.4 mm.

Figure 4 (
c) shows the deformation of unit cell 1 during the test and predicted by the simulation, which demonstrates the accurate predictions from the FE model.To compare the performance of different designs, the effective stress-strain curves (figure 4(b)) are derived from the

Figure 4 .
Figure 4. (a) Comparison of the force-displacement curves obtained from the uniaxial compression tests (the shade shows the standard deviations from three repeating testings) and numerical simulations of unit-cells.(b) Comparison of the effective stress-strain curves among unit cells, derived from the experimental results.(c) Deformation of unit cell 1 (α = 60 • ) at the displacement of 0 mm (initial state), 0.5 mm (initiation of buckling), 1.0 mm, 2.0 mm, 3.0 mm, and 4.0 mm (initiation of densification).

Figure 5 .
Figure 5. (a) Decomposition of applied force at the end of the inclined beam; simulation results regarding rotation of the beam ends within the elastic range for unit cells (b) bottom end (c) top end.

Figure 6 .
Figure 6.Experimental and numerical results for: (a) the two-by-one structures of oriented pairs of unit-cells (specimen pictures shows deformation of structures after the first peak load); (b) two-by-one structures mirrored pairs of unit-cells; and (c) two-by-two structures (generated from unit cell 2), with stress distribution corresponding to the two snap-through phases (phase AB and phase CD).

Figure 7 .
Figure 7. (a) Force-displacement curves obtained from experiment and FE model using a symmetric mesh generated from a quarter of the structure, with the deformation showing top and bottom beams buckling simultaneously.(b) Force-displacement curves comparison between experiment and FE model using an asymmetric mesh.(c) Printing path within the specimen as described using a symmetric mesh, corresponding to simulation in (a).(d) Printing path for the as-fabricated specimen, corresponding to experiment and simulation in (b).

Figure 8 .
Figure 8.(a) Visualization of unit cell manufacturing procedure, taking unit cell 2 (α = 45 • ) as an example.Blue colour represents the wall, white colour indicates the infill, and black colour is the printing platform.(b) Micro-CT binary images showing the middle-height cross-section for unit cells (from left to right: unit cell 1, unit cell 2, and unit cell 3).Pores are indicated by the black colour, while solid parts are represented by the white colour.Two types of pores are defined using magenta colour (pores in between walls and infills) and bright green colour (pores among infills).

Figure 10 .
Figure 10.Force-displacement curves from FE models for three-by-three gradient structures with different combinations of unit cell 1, unit cell 2, and unit cell 3: (a) considering purely oriented or mirrored ('xxx' denotes the unit-cell number in the order of bottom-to-top stacking sequence); (b) considering both oriented and mirrored cases (xxx_x denotes the unit-cell number in the order of bottom-to-top stacking sequence, while the number after '_' is 'the number of the unit cell that orients differently from the other two').Note: simulation graphs only show left-half of structures as the FE models were built by applying a symmetry boundary condition along y-axis.

Figure 12 .
Figure 12.(a) Schematic drawing for a multiple-layered structure (periodic/graded), with number of Nx unit cells in a row and Ny in a column.The maximum angle of tilted beam within the structure is αmax.(b) Design map showing the peak force and absorbed energy versus the tilted beam angle throughout the snap-through phenomenon for a single unit cell.

Table 1 .
Summary of material type and property for flexible beam fabrication in literatures.

Table 2 .
Critical buckling load (Pcr) comparison among the experiment, numerical simulation, and analytical calculations.

Table 3 .
Thickness measurements of the tilted beams compared to design values.Note: viewpoint (left/right and top/bottom) is set when looking down at the printer platform.