Identification of mechanically representative samples for aperiodic honeycombs ✩

Honeycombs are a class of metamaterials widely used in engineering. Traditionally they have been comprised of periodic arrays of hexagonal, triangular and square cells however there have been many studies into novel variations of these lattices. Recent works have focused on the mechanical properties of metamaterials based on lattices with aperiodic order. This promising new field of study offers the potential for a wider range of available mechanical properties. With this potential comes challenges in producing, testing and simulating these structures. Analysis cannot be carried out on a unit cell as aperiodic patterns lack translational symmetry. Furthermore, it is unknown how much of the structure is required for a representative sample. This study uses a statistical approach to investigate how patch size influences the ability to accurately estimate the mechanical properties of aperiodic honeycombs. By exploiting a numerical framework requiring minimal computational resources, 1600 simulations were carried out on randomly sampled patches. This was supported by mechanically testing a targeted set of 40 additively manufactured honeycombs. It was found that in most cases increases in patch size resulted in consistent reductions in variation of properties with the samples varying by less than 5.3% from the mean when considering Young’s modulus


Introduction
Mechanical metamaterials are engineered materials that are designed to exhibit desirable properties that differ from the base materials from which they are made [1].Amongst the range of modern metamaterials are a class of cellular solids known as engineering honeycombs.Comprised of two dimensional lattices that have been stacked into the third dimension, honeycombs have found many applications in engineering due to high specific stiffness [2] and energy absorption [3].Effective mechanical properties of these cellular solids are attributed to the geometry of the unit cell along with the properties of the base material [4].Traditionally they have been comprised of periodic arrays of cells such as the hexagon, square and triangle [5] however there have been many studies into novel variations of these lattices [6].These shapes tessellate to tile the plane periodically meaning they display translational symmetry as well as a range rotational symmetries limited to the 2, 3, 4 and 6-fold associated with standard crystallographic structures [7].Fuelled by the development of additive manufacturing, an abundance of novel variations of these honeycombs have been developed using a number of strategies.The study of these cellular structures are well established and effective mechanical properties are often determined numerically, analytically and experimentally by considering a unit cell [4,8].
In 1976 Roger Penrose [9] discovered a system of two shapes that could tile the plane aperiodically meaning the resulting pattern was ordered yet lacked translational symmetry.Although lacking translational symmetry these patterns display rotational symmetries other than the 2-, 3-, and 6-fold symmetries displayed by periodic lattices.Then in 1984, Dan Schechtman [10] discovered quasicrystals, aperiodic arrays of atoms that contradicted the long held understanding that crystals only existed as periodic arrays of atoms.Along with surprising elastic, conductive and thermal behaviours [11], these quasicrystals were also found to be more isotropic than conventional crystals [12].In recent years, this has led to studies into the behaviour of mechanical metamaterials based on aperiodic patterns.Moat et al. [13] investigated the mechanical properties of honeycombs based on the Penrose P3 tilings experimentally and demonstrated that where dimensional stability under plastic deformation is important, the P3 honeycombs showed significant differences in deformation modes when compared with periodic honeycombs.Kim et al. [14] investigated the strength, toughness and stiffness of composites based on the P3 and compared https://doi.org/10.1016/j.mtcomm.2024.108505Received 3 January 2024; Received in revised form 12 February 2024; Accepted 27 February 2024 them with equivalent periodic lattices.It was found that the P3-based composites exhibited the highest degree of isotropy in elastic response.Somera et al. [15] classified the mechanical behaviour of a range of aperiodic lattice structures as either bending dominated or stretch dominated based on numerical simulations.Imediegwu et al. [16] used asymptotic expansion homogenisation (AEH) to investigate the mechanical properties of seven aperiodic lattice structures and showed that all offered near isotropic properties when considering elastic modulus and Poisson's ratio.Clarke et al. [17] introduced an isotropic zero Poisson's ratio metamaterial based on the 'hat' monotile while Jung et al. [18] investigated composites constructed with the same lattice and found they exhibited greater strength, toughness and stiffness when compared to periodic lattices.Wang et al. [19] introduced a novel lattice structure based on the 'hat' monotile and observed high strength, toughness and damage tolerance through finite element analysis and compression tests on additively manufactured samples.Rieger and Danescu [20] explored the macroscopic elasticity of structures based on the 'hat' tiling whilst considering the influence of the size of the domain on the mechanical behaviour of the structure.
One formal requirement of an aperiodic tiling is that it can tile an infinite plane.When used as a basis for metamaterials, this is desirable as a continuous material can be produced without maximum size restrictions.In order to maintain aperiodicity, tile configuration varies throughout the pattern, however this variation creates challenges when attempting to define a size of patch or domain that is representative of the continuous pattern.In this context, the size of a patch is a measure of the number of tiles in the patch, which can be formally defined as the ratio of the area of the tile to the area of the patch.When using these lattices to construct cellular materials, where tiles become cells, this may become problematic.If an engineer were to design a component based on the mechanical properties of this material, they first must ascertain whether they can make that assumption based on the relative scales of the component and cellular structure.Unlike for periodic lattices, a unit cell cannot be defined for lattices based on aperiodic tilings.Instead, individual tiles form groups of tiles that repeat but not periodically.Furthermore any different finite patch of cells may yield different mechanical properties due to number and orientation of different cells.These issues become apparent when attempting to investigate the mechanical properties of aperiodic cellular structures.When selecting a patch of an aperiodic structure for mechanical testing or numerical analysis, it must be ascertained whether any results obtained are true for a large patch of the structure.Limitations to the sample size that can manufactured and tested, may require an estimation of error when analysing results for smaller samples.
This work explores how patch size influences the ability to accurately estimate the mechanical properties of continuous aperiodic honeycombs from finite patches.By randomly sampling many patches of varying sizes from different locations within a large patch, a statistical approach is applied based on the premise that when a representative patch is achieved, variance in properties will converge to zero.This is made possible by using an adaption of AEH [21], allowing for the collection of elastic modulus and Poisson's ratio data from hundreds of simulations for each lattice.In addition, mechanical tests on additively manufactured samples are used to validate these results and show the extent to which yield strength and compressive strength vary for a reduced set of targeted samples.Furthermore, the numerical and experimental methods of extracting Young's modulus and Poisson's ratio do not imitate each other with respect to loading, and any variance related to the method may be significantly different.With the emergence of studies into this promising class of material, patch size studies are an important step in the design of experiments and simulations.Previous studies into the mechanical properties of structures based on aperiodic patterns have not identified the extent to which their results are representative, and so cannot be assumed to apply to all structures based on those patterns.The work presented here explores, for the first time, the criteria which need to be met to ensure that results from simulation and testing are representative for a range of patterns.

Generation and selection of samples
Methods for generating periodic tessellations cannot be applied to aperiodic tessellations.This is because no unit cell can be defined for an aperiodic tiling and there are no local matching rules that ensure assembly with long range order.Two main approaches are used to generate aperiodic tilings, the ''cut-and-project'' and the ''substitution'' methods.For simplicity, the second method is used in this work.The substitution method forces aperiodicity by applying substitution rules to each prototile.Aperiodic tilings can be generated by starting with an initial prototile, and subdividing it into a set of prototiles obeying the substitution rules [7].Sets of prototiles are then subdivided again using the same rules until the required patch resolution is obtained.Substitution rules for the Penrose-P2, Penrose-P3, Robinson Triangle, Tübingen Triangle, Ammann-Beenker, Pinwheel and Square-Triangle tilings are described in detail in previous works [16].As an example and for completeness, substitution rules for the 12-fold Socolar tiling can be seen in Fig. 1.The prototile types are shown in rows (A, B, C, D) and prototile generations are shown in columns (I, II, III).The process of substituting a prototile for a group of smaller prototiles is referred to as ''deflation'' in this paper.Patterns were selected based on the availability of substitution rules whilst also ensuring the patterns were geometrically diverse with respect to symmetry and cell shape.The patterns have a range of rotational symmetries, the Socolar, Ammann-Beenker and Square-Triangle have 12-, 8-and 6fold symmetries respectively.The Penrose-P2, Penrose-P3, Robinson Triangle and Tübingen Triangle all display 5-fold symmetry yet have differing cell geometries.Although the Robinson Triangle and Tübingen Triangle patterns are both constructed with golden triangles and golden gnomons, they exist in different relative scales leading to distinct patterns.The Pinwheel is considered to have infinite rotational symmetry due to the rotation of deflated prototiles by an irrational angle.
A bespoke Python program, developed by the authors, was used to generate all samples.The geometries for all cells constituting each lattice are coded into the program with attributes such as a list of vertices along with methods for substituting tiles, rotating vertices and shrinking cells.Cell objects are initiated with four arguments, the xcoordinate of the cell centre point, the y-coordinate of the cell centre point, the length of a cell wall and the angle of rotation.After selecting the type of lattice, cell wall thickness, target density (  ), size of sample and rotation, an initial cell much larger than the sample is created and deflated until the resulting area of the largest cell is smaller than a predefined maximum area whilst any cell completely outside the sample is discarded.The cells are shrunk so that cell walls are of the required thickness, in this case defined by the capabilities of the 3D printing equipment discussed later.The initial cell wall length is then varied within an optimisation loop so that the total area of the shrunken cells (  ) and sample cross-sectional area (  ) is such that (  −   )∕  ≈   .This initial cell wall length is then used for all subsequent samples.A file containing vertices of cells prior to shrinking is exported for FEA simulations or 3D models of experimental samples are produced using FreeCAD Python modules.Top and bottom plates are added when 3D models are generated so that the samples can be mechanically tested.
In order to investigate how properties vary for different patches depending on the size of the patch, AEH was used for the reduction in computational resources compared with other techniques, and its adaptability to aperiodic lattices [21].50 randomly selected samples were generated from a large patch for sample sizes of 20, 30, 40 and 50 mm.These sample sizes were chosen to aid practical comparisons and effectively vary the amount of the pattern contained within the sample.One form of quantifying the size of a patch is by calculating the cell area (  ) to sample area (  ) ratio (  ∕  ).Values for all samples and sizes can be found in Table 1 where the area of the largest cell has been used for calculations.Samples can be scaled to the same size since the physical size is not of interest; it is the resolution of the pattern that is being investigated.Samples were generated for each aperiodic tiling by starting with an initial cell orders of magnitude larger than the sample, the cell having a wall length optimised to give a relative density of 0.45 with a cell wall thickness of 0.5 mm on the level of deflation where largest cell ceases to be larger than the predefined maximum area (the correct quantity of deflations have been carried out).The resulting pattern is always the same, as the same rules are applied an equal amount of times on the same initial cell.It is therefore not necessary to generate this large patch by acknowledging the fact that any cells completely outside a sample window can be ignored after each deflation without changing the cells that fall withing the sample area.This is true as long as deflated cells do not overlap the original cell as with the Socolar.Here, the sample window is temporarily enlarged by a fraction of cell wall length for the current deflation.On translating the initial cell by randomly generated vectors and fixing the sample window at (0, 0), samples are taken from a significantly larger patch without the need to generate such a large patch.Fig. 2 gives a graphical depiction of this, the large patch has been deflated seven times, in reality many more deflations are carried out.
The above approach is not reasonably practicable for experimental samples due to the time required to print, prepare and test such quantities.Instead, five different samples of each tiling were selected from a finite patch, with the central patch labelled as ''sample 0''.All samples were 50  50 mm and were taken from the same location for each tiling.Samples were visually inspected to ensure significant variation in pattern, although aperiodic tilings lack translational symmetry, sections of tiling do repeat throughout a continuous patch [7].Fig. 3 shows where these samples were extracted from the original 150  150 mm patch.

Modelling
Asymptotic expansion homogenisation (AEH) [22] is a widely used technique to obtain the effective or average properties of materials that can be assumed to be built up of periodic arrays of micro-architecture.Effective properties for lattices can be approximated with a reduction of the problem size leading to significantly less computational resources required in comparison with other FEA methods.Analysis of large quantities of samples becomes viable using this technique.Moreover, AEH offers significant advantages for analysis of aperiodic cellular structures as other methods rely on the existence of a unit cell [21].The analysis of aperiodic structures is made possible by the addition of a homogenised region (-region) surrounding the sample allowing for the application of periodic boundary conditions shown in Fig. 4(a).The macroscopic properties from an initial iteration are assigned to the added -region and the macroscopic properties calculated again.This process is repeated until convergence is achieved.
Numerical simulations of 1,600 honeycombs were carried out by element-based assignment of material properties using FEniCS [23,24], an open source Python package for finite element analysis.Here the structures are treated as composites consisting of material and void, the void being assigned as weak isotropic material, in this case E = 100 Pa.Files for each sample are produced using the pattern generating software as detailed in Section 2.1.Each file containing the vertex coordinates of all cells for one sample are read by the simulation software.A uniform triangular mesh is created and material properties are then assigned to each mesh element based on proximity to the cell facets, with fractional properties assigned to elements that fall partially within the structure.The Young's modulus of the base material was taken to be 3292 MPa based on orientation dependant values supplied by Ultimaker.The material assignment process is described in detail within recent works by Imediegwu et al. [21] along with a mesh convergence study and convergence study of -region size.An example of material assignment can be seen in Fig. 4(b).The number of mesh elements was set to ensure struts had seven elements across their width.The honeycombs were simplified to 2-dimensional lattices with displacements obtained using Hooke's law for plain strain conditions.An effective elasticity matrix for each honeycomb was obtained from simulations and the effective Young's modulus and Poisson's ratio were extracted from the resulting compliance matrix.By using this method analysis can be automated and all resulting samples have the same individual cell size, FEA meshes are equivalent and each cell wall has the same amount of elements across it.Essentially, the only parameter that is varied is the individual cell area to patch area ratio.

Production of samples
Physical samples were produced as set out in Clarke et al. [5] using Fused Deposition Modelling (FDM).FDM is an extrusion based additive manufacturing (AM) system where thermoplastic polymers are extruded from a nozzle.The nozzle lays down material following a toolpath and the model is built up layer by layer.FDM can produce  cost effective models quickly, however the method often introduces anisotropy into models [25].Moreover, any printed specimen may exhibit different mechanical properties depending on toolpath and print parameters such as layer height [26].Datasheet values for mechanical properties do not always yield equivalent results when used for FEA models and compared with experiments [27].For this reason care must be taken to maintain consistency and ensure that any comparisons are valid.
Forty physical samples were produced with dimensions 50  50  50 mm in accordance with ASTM D1621 [28].An example of a Penrose-P3 honeycomb can be seen in Fig. 5(a).Cura 4.12 [29] was used to slice the model for printing and a cell wall thickness of 0.5 mm was chosen as it yielded toolpaths that were consistent for all lattices and build quality was satisfactory for mechanical testing.The samples were then printed on an Ultimaker S3 FDM printer with a nozzle diameter of 0.4 mm.All samples were produced in Ultimaker white polylactic acid (PLA) and were built up of 0.2 mm layers to achieve the required sample thickness of 50 mm.

Mechanical testing
Mechanical testing was carried out on an Instron 6800 mechanical testing machine with a 50 kN load cell.All samples were compressed to a maximum strain of 0.05 at a constant strain rate of 0.01/min using crosshead displacement measurement.Strain was measured using digital image correlation (DIC) for greater accuracy and to negate the need to account for machine compliance.DIC is a full-field optical strain measurement method where images acquired during tests are compared with an image of an undeformed sample and local displacement vectors are measured by tracking subsets of pixels [30].In order for the algorithm to achieve a unique solution, these subsets must contain random features, achieved by applying a random speckle pattern to the samples with black spray paint.An undeformed image of each pattern can be seen in Fig. 6.Images and data were collected at a rate of 1 Hz using LaVision DaVis 10.2.1.81611[31] software.Clarke et al. [5] used live virtual linear strain gauges where subsets of pixels at each end of the strain gauge are tracked in real time removing the need to process entire images.In that work, normal strain was measured with a vertical gauge and Possion's ratio was successfully calculated with the addition of a horizontal gauge.The periodic nature of these honeycombs enabled such a strategy due to the repeating rows of cells.Aperiodic honeycombs lack repeating rows of cells and so it becomes unclear where these horizontal gauges should be placed.For this reason, images from the elastic region were processed with subsets sizes of 15 and step sizes of 4 enabling horizontal strain to be measured using a box strain gauge spanning all complete cells horizontally and approximately 15 mm of central cells vertically.The vertical gauge is placed between the bottom and top surfaces to measure the normal engineering strain.Fig. 5(b) shows the location of both strain gauges.The Poisson's ratio is calculated by dividing the horizontal strain from the box gauge by the strain obtained from the vertical strain gauge.Poisson's ratio is quoted as the mean of all values over the elastic region.The effective compressive strength is defined as the first local maximum in the stress strain curve.The effective Young's modulus and the effective yield strength were calculated using a bespoke Python algorithm: the Young's modulus by finding the slope of the linear region; and the yield strength by finding where the chord modulus ceases to equal the tangent modulus.

Statistical variation of cell orientation
Statistical data was collected for all modelled samples in order to investigate the extent to which patches vary with respect to geometric arrangement of cells.This variance can be demonstrated by considering the number of occurrences of each shape at each orientation.Fig. 7 shows this using the Robinson triangle lattice as an example.For the purposes of generating tilings, the Robinson Triangle tiling has four tiles, the golden triangle, the golden gnomon and their respective reflections.The reflected cells are equivalent with respect to structural integrity so the resulting honeycombs can be considered to have two different cells at ten orientations each.On calculating the mean occurrence of a tile at each orientation for all sets of samples, the maximum percentage deviation from the mean for that sample set and tile orientation was calculated.Each point in Fig. 7 represents the maximum percentage deviation from the mean observed for a cell orientation, the sets being the 50 randomly generated samples generated for each patch size.It can be seen that the range of values decreases along with the median as patch size is increased.There is significant scatter in maximum deviation from the mean obtained for tile orientations in 20 mm samples.Five golden triangle tile orientations had a maximum deviation of 100% from the mean due to the absence of this orientation in some patches.Nine of the ten golden gnomon orientations had maximum deviations of between 61% and 65% (to 2 s.f).This distinct contrast between the shapes can be attributed to the higher frequency of the gnomons within a large patch where the ratio of gnomons to triangles tends to the golden ratio.This analysis is not exhaustive, there are other metrics such as occurrences of edge connectivities or occurrences of groups of tiles however it is sufficient to show how the geometry of the patches vary with patch size.Plots for the P2, P3, Ammann-Beenker, Tübingen triangle, Square-Triangle and Socolar can be found in Supplementary data.It is however not appropriate to apply this strategy to the Pinwheel.For all other lattices, there is only a finite number of tile orientations that are possible whereas for the pinwheel, an infinite patch will contain an infinite amount of tile rotations.It is therefore more appropriate to consider the number of tile rotations in a patch and noting this will tend to infinity with increasing patch size.Fig. 8 shows the number of tile orientations found in each randomly sampled patch of Pinwheel tiling.Each point indicates that at least one sample contained the corresponding quantity of tile orientations.Larger points indicate more samples contained the corresponding quantity.It can be seen that the spread in results does not vary significantly with patch size, however the maximum and minimum number of orientations increase substantially.More importantly, the total number of tile orientations observed for all 20 mm samples was 50 whereas the maximum number observed in any individual sample was 20.This suggests that not only can a sample have a different quantity of tile orientation, the tile orientations constituting a sample may be completely different from another.At patch sizes of 50 mm the maximum number observed in one sample was 35 and the total number of tile orientations observed for all samples was 62.

Variations in elastic properties obtained numerically
In Figs. 9 and 10, statistical data for the distribution of mechanical properties obtained through simulations of randomly sampled patches is presented in box plots.In addition, swarm plots of the same data set are superimposed onto the box plots with each point representing one sample.All plots in each column have the same axis for comparison.The first column of plots shows the percentage Young's modulus () deviated from mean (100( −   )∕  ) and the second column shows the percentage Poisson's ratio () deviated from the mean (100( −   )∕  ).For the plots in the third column, an adaptation of the isotropy factor presented by Zhu et al. [32] is used.In this case only orthogonal values for Poisson's ratio ( 0 and  90 ) are considered.The isotropy factor is defined as where  0 is the Poisson's ratio at 0 • and  90 is the Poisson's ratio at 90 • for the same sample.Values range between 0 and 1 with 0 being an isotropic material.
For patch sizes of 20 mm, the Young's modulus of P3 samples vary between −16.33% and 13.09% of the mean, this reduces to between −6.51% and 3.87% for 40 mm samples.There is however an increase for 50 mm samples with samples varying between −6.45% and 6.73%.This is also true when considering the first and third quartiles where 50% of samples ranged between −1.08% and 1.77% of the mean for 40 mm samples and between −2.22% and 1.64% for 50 mm samples.Increase in patch size does however provide consistent reductions in deviation when considering Poisson's ratio with values ranging between −7.4% and 9.35% of the mean for 20 mm samples and −2.97% and 2.89% for 50 mm samples.This reduction is also apparent in the first and third quartile values with 50% of samples falling between −0.78% and 0.81% when the patch size is increased to 50 mm.Larger samples also show a reduction in maximum isotropy factor with values falling from 0.09 to 0.03 for 20 mm and 50 mm samples respectively.50% of samples have an isotropy factor of less than 0.02 for 30, 40 and 50 mm patch sizes.The P2 honeycombs show much less variation in Young's modulus than the P3 honeycombs throughout the size ranges, with samples ranging between −5.52% and 4.74% from the mean for 20 mm samples and −3.8% and 2.34% for 50 mm samples.The P2 also shows consistent reduction in spread of data when increasing patch size with respect to maximum/minimum and interquartile range.This trend is not observed for Poisson's ratio as there is an increase in spread when patch size is increased from 40 mm to 50 mm due to two low outliers.There is however a slight reduction in interquartile range with 50% of samples falling between −1.5% and 1.64% of the mean for 40 mm samples and between −1.5% and 1.61% for 50 mm samples.There is also a steady reduction of maximum isotropy factor with increasing patch size, values fall from 0.09 for 20 mm patches to 0.02 for 50 mm patches.

When considering Young's modulus and Poisson's ratio for Robinson
Triangle honeycombs, consistent reductions in deviation from the mean are observed with increasing patch size.Like the P3 honeycombs, there is a significant reduction in spread of Young's modulus when increasing patch size from 20 mm to 50 mm with values of −13.79% to 13.76% and −4.72% to 3.3% respectively.The same is observed for Poisson's ratio with deviations from the mean varying between −17.22% and 11.43% for 20 mm patches and between −6.34% and 5.09% for 50 mm patches.There is however a slight increase in interquartile range between 40 and 50 mm patch sizes.Maximum isotropy factors fall from 0.12 for 20 mm patches to 0.03 for 50 mm patches with values increasing slightly between 30 and 40 mm patch sizes due to high outliers, there is however a reduction in interquartile range.
Box plots for the Tübingen Triangle honeycombs show a steady reduction in deviations from the mean with regards to Young's modulus.Values reduce from between −13.44% and 9.39% to between −3.21% and 2.09% for 20 and 50 mm patch sizes respectively.Values for deviation from the mean of Poisson's ratio vary between −17.65% and 12.27% for 20 mm samples, this reduces to between −3.98% and 4.17%.There is however a slight increase in spread between 30 mm and 40 mm.Maximum isotropy factors reduce from 0.11 for 20 mm patches to 0.02 for 50 mm patches.A slight increase in interquatile range can be observed between 30 and 40 mm samples.
The Ammann-Beenker honeycombs show the expected reduction in maximum/minimum range with patch size for all three properties.Deviation from the mean of Young's modulus varies from between −13.44% and 14.52% for 20 mm patches to between −5.3% and 4.84% for 50 mm patches.The reduction in spread of Poisson's ratio is not as profound as with Young's modulus with no significant improvement until patch size is increased to 50 mm.Spread in data reduces from between −10.07%and 7.72% to between −4.38% and 3.53% for patch sizes 50 mm.Maximum isotropy factor falls from 0.08 for 20 mm samples to 0.05 for 50 mm samples where 50% of samples fall between 0.01 and 0.02.
The Pinwheel honeycombs show the most substantial variation in properties of all the patterns.At 20 mm patch size, deviation of Young's modulus from the mean varies between −26.58% and 34.23%.Although increasing patch size reduces this variation, Young's modulus still varies between −13.57% and 12.8% for 50 mm patches.Even higher deviations from the mean can be observed for Poisson's ratio.Variations between −37.61% and 38.64% can be seen for 20 mm patches and this is reduced to between −18.66% and 17.07% for 50 mm patches.Isotropy factors are as high as 0.21 for 20 mm samples.This maximum only reduces to 0.12 for 50 mm samples due to a high outlier, however 75% of samples have an isotropy factor of below 0.04 for this size.
The Square-Triangle honeycombs show consistent reductions in maximum deviation from the mean with increasing patch size for all properties.At 20 mm patch size, Young's modulus varies between −9.07% and 8.8% of the mean.When patch size is increased to 50 mm, the Square-Triangle shows the least maximum deviation from the mean of all the patterns with variations of between −2.45% and 1.83% observed.Poisson's ratio deviates from the mean by between −7.28% and 8.72% for 20 mm patches and between −3.68% and 2.7% for 50 mm patches.The Square-Triangle displays consistently low isotropy factors with a maximum value of 0.05 for 20 mm patches and 0.01 for 50 mm patches.At 50 mm, 75% of samples had an isotropy factor of less than 0.005.
Like the Square-Triangle honeycomb, the Socolar honeycomb shows limited deviation from the mean for both Young's modulus and Poisson's ratio.Young's modulus varies between −8.42% and 5.55% for 20 mm patches and between −3.13% and 2.73% for 50 mm patches.Consistent reductions in variations are observed with increasing patch size.Poisson's ratio deviates from the mean between −6.74% and 7.28% for 20 mm samples and displays the lowest maximum/minimum range of between −2.71% and 2.69% when compared with other patterns.Although the range in data remains almost constant for 30 and 40 mm samples, consistent reductions in interquartile range can be observed with increasing patch size.Maximum isotropy factor reduces from 0.07 to 0.02 for 20 and 50 mm samples respectively with 75% of samples having an isotropy factor of less than 0.01 for the latter.Increasing patch size provides consistent reduction in isotropy factor.
As expected, patch size greatly influences the spread of properties obtained through randomly sampling patches, however the effect varies for each pattern.For some patterns, increases in patch size provide consistent reductions in spread, whereas for others, increases in spread are observed.The minimum isotropy factor for all patterns and all patch sizes is approximately 0. This suggests that at least one orthogonally isotropic sample was found for each group of samples.

Young's modulus, yield strength and compressive strength
Fig. 11 shows stress/strain curves obtained from mechanically testing 50  50  50 mm PLA honeycombs.Results from testing five  There is relatively good agreement in mean Young's modulus between experiments and modelling for the Robinson triangle honeycomb.A value of 569.7 MPa was obtained experimentally and 526.6 MPa was obtained from simulations.A similar spread of Young's modulus data was also observed for experiments and simulations.Yield strength and compressive strength varied very little with the exception of sample 3. Yield strength varied between 5.05 and 5.64 MPa, and compressive strength varied between 9.51 and 9.8 MPa for samples 0, 1, 2 and 4.These values drop to 4.55 MPa and 8.52 MPa for sample 3. A different shape plot can also be observed for this sample with a flatter appearance after yield as opposed to other samples showing a steep reduction in stress after the first local maximum.
As with the Robinson triangle honeycombs, a good agreement in mean Young's modulus between experiments and modelling for the Tübingen Triangle honeycombs is observed.The mean values obtained from experiments and simulations are 637.6MPa and 614.7 MPa respectively.There is however a significantly higher spread of Young's modulus data for experimental results.This variance can also be seen in results for yield strength and compressive strength.Yield strength deviated between −21% and 10.3% of the mean and compressive strength deviated between −7.88% and 6.9%.The appearance of plots also vary between samples with sample 4 showing a steep fall in stress after the first local maximum, sample 1 showing a shallower reduction and the remaining samples falling in between.
The Young's modulus of experimental samples varies very little for the Ammann-Beenker honeycombs.A maximum of 355.35 MPa and a minimum of 340.8 MPa was measured equating to a deviation of between −1.81% and 2.71% from the mean.The Mean value of 347.1 MPa obtained from experiments does however vary significantly from the value of 228.4 MPa obtained from simulations.Samples vary very little with respect to yield strength and compressive strength.Yield strength varied between 3.12 and 3.34 MPa relating to deviations from the mean of −4.06% and 2.71%.There was even less variation in compressive strength with values of between 6.53 and 6.73 MPa measured, equating to deviations from the mean of between −3.01% and 2.16%.Plots differ slighty after the first local maximum with some samples showing a steeper decline in stress with strain.Experimental Young's modulus results for the PinWheel honeycombs mostly fall within the range predicted by modelling, with only 2 samples falling marginally outside.This reflects the large variations in properties obtained from simulations rather than low variations in experimental results or a similar mean obtained from both methods.Experimental results fall between 432.6 and 492.2 MPa and mean values are 460.6 and 417.19 MPa for experiments and simulations respectively.Moderate variations in Yield strength and compressive strength are observed, samples yielded between 3.56 and 3.88 MPa (between −4.6% and 3.96% of the mean) and values of between 6.81 and 7.62 Mpa (between −5.55% and 5.69% of the mean) were obtained for compressive strength.Samples also show very different plots after yield.
The Square-Triangle honeycombs show the least variation in Young's modulus with regards to modelling results however this is not apparent in experimental results.Values vary between 590.95 and 746.16MPa, a deviation from the mean of between −11.58% and 11.63 MPa%.Three samples do however lie within 0.5% of the mean.Mean values from experiments and simulations are in close agreement with values of 668.4 and 645.9 MPa found respectively.Variations in yield strength are between 5.74 and 6.71 MPa equating to a deviation from the mean of between −8.13% and 7.39%.Much lower variation in compressive strength is observed with values of between 10.81 and 11.4 MPa obtained.This equates to a deviation from the mean of between −2.44% and 2.89%.All plots display a very similar shape.
The Socolar honeycombs show more spread in results for experiments when compared to simulations.Young's modulus varies between 409.76 and 477.64 MPa, representing deviations from the mean of between −6.57% and 8.90%.Mean values differ greatly with 438.6 MPa measured experimentally and 343.9 MPa obtained from simulations.Yield strength and compressive strength show similar variation in values as Young's modulus.Samples yielded between 3.77 and 4.17 MPa, a deviation of between 4.07% and 6.11%.Compressive strength varied between 7.24 and 8.05 MPa representing deviations from the mean of between −5.85% and 4.68%.Plots for all samples show a similar overall shape.

Poisson's ratio
Fig. 12 shows plots of Poisson's ratio versus normal strain obtained from mechanical testing five samples of each pattern.The shaded area represents the spread of Poisson's ratio data obtained from simulations.In most cases, Poisson's ratio increases with normal strain.For simplicity, the average value over the elastic region is quoted when making comparisons.Broadly, modelling and experiments show remarkable agreement.
For the P3 honeycombs, Poisson's ratio varies between 0.63 and 0.73, a deviation from the mean of between −4.8% and 10.27%.This is greater than that from modelling however maximum values were equal and mean values were 0.66 an 0.71 for experiments and simulations respectively.
Experimental Poisson's ratio values for the P2 honeycombs lie within the range predicted by modelling with exception of Sample 1. Values vary between 0.24 and 0.29 with a mean value 0.28 compared with a value of 0.30 obtained from simulations.
For the Robinson-Triangle honeycombs, experimental results show more variation than results from simulations with Poisson's ratio varying between 0.29 and 0.37.Mean values of 0.32 and 0.38 were obtained from experiments and simulations respectively.Poison's ratio results for the Tübingen Triangle are very similar to the those of the Robinson-Triangle.Values vary between 0.29 and 0.36 with a mean value of 0.33 compared with a value of 0.37 obtained from simulations.
Experimental results for the Ammann-Beenker honeycombs show more variation in Poisson's ratio than those obtained from simulations.Values vary between 0.68 and 0.81 for experiments and between 0.69 and 0.75 for simulations.Mean values from both methods are in close agreement with values of 0.74 and 0.72 obtained from experiments and simulations respectively.
The Pinwheel is the only pattern where all experimental results fall within the range predicted by modelling.Again, this reflects the large variations in Poisson's ratio obtained from simulations.Values of between 0.27 and 0.33 were acquired experimentally whereas values from modelling fell between 0.25 and 0.35.Mean Poisson's ratio for both methods were equal (to two decimal places) at 0.30.
Simulations of Square-Triangle honeycombs showed very little variation in Poisson's ratio, this is also true for experimental results.Values varied between 0.32 and 0.35 for experiments, a deviation of between −4.19% and 4.80% from the mean.There is however a substantial difference in mean Poisson's ratio where values of 0.33 and 0.41 were obtained from experiments and simulations respectively.There was no overlap in experimental and numerical results.
For the Socolar honeycombs, more variation in experimental results can be observed than for modelling results.Poisson's ratio varied between 0.54 and 0.66 with the mean value of 0.58 falling just outside the range of 0.62 to 0.65 predicted by modelling.

Discussion
In general, increasing patch size provides a reduction in the variation of mechanical properties obtained from simulating large quantities of samples.In most cases deviations from the mean are reduced to below ±5%.This reduction, in part, can be attributed to statistical uniformity of the samples where the quantity and orientation of features contained within the patch approaches that of a continuous patch.Consideration must also be given to the method of how these results are obtained.Due to the nature of the lattices, there are always partial cells at the sample edges.These are formed by cropping cells so they fit within the sample boundary.Although every effort has been made to ensure equivalence of all samples, inevitable differences in how partial cells contribute to the effective properties are introduced.They may have a positive or negative effect on Young's modulus depending on the size and orientation of the partial cell wall.These effects are hard to quantify however it can be assumed that they reduce as cell area (  ) to sample area (  ) ratio is reduced.As   ∕  is increased, the proportion of the structure comprised of partial cell walls is more likely to increase.As an example, the mean percentages of partial cell walls for the P3 are 14.8%, 10.2%, 7.4% and 6.0% for 20, 30, 40 and 50 mm patches respectively.This trend was found for all the simulated lattices.The effect of partial cells on the range of mechanical properties is difficult to demonstrate with aperiodic lattices as isolating this effect is challenging, changing   ∕  not only will change the percentage of partial cells, but also the cells that contribute to the structure.It can however be demonstrated with periodic lattices.It is not appropriate to randomly sample a periodic pattern as any variation in sample can be achieved by translating the sample window by distances of less than the period of the pattern.For these purposes it is however simpler to demonstrate this with isotropy factors.Fig. 13 shows isotropy factors for triangle and hexagon lattices at varying sizes.These lattices have been chosen as analysis on unit cells show them to be isotropic [4].As 90 • is not a multiple of the angle of rotational symmetry for these lattices, partial cells are different on each orthogonal boundary.The samples have 2 fold rotational symmetry and this is reflected in the fact that they are not orthogonally isotropic.Isotropy can however be increased by increasing patch size.It is worth noting that isotropy factors for the hexagon are in the region of the maximum values found for the Square-Triangle.
Increasing patch size, in general, provides consistent reduction in the variation of mechanical properties with variations reaching acceptable values for all but the P3 honeycombs and Pinwheel honeycombs.The P3 honeycomb does not show the same reduction in Young's modulus deviation from the mean as with other patterns.It is unclear whether this is because the patch must be much larger to incorporate more cells and therefore a more representative configuration of cells.P3 tilings also contain continuous one dimensional features known as worms or ribbons [33].If these features have a significant effect on mechanical properties, the presence and orientation of them may be more consequential than providing a statistically representative patch.Young's modulus and Poisson's ratio values obtained for the Pinwheel honeycombs vary significantly.The maximum and minimum values for 50 mm patches are approximately equal to the worst case 20 mm samples for other patterns.This is also true for maximum isotropy factors.The substitution rules for the Pinwheel tiling give rise to copies and reflections of the original tile rotated by an irrational angle.This property leads to the possibility of an infinite amount of distinct rotations within an infinite patch.This may lead to the hypothesis that such variation in orientation might give rise to higher isotropy, however 50 mm patches with these parameters can contain cells from only four deflations leading to a limited number of rotations.As stated earlier patches could contain cells at completely different angles than other patches.Substitution rules maintain the edges of triangles produced by previous deflations leading to continuous features spanning samples in many possible orientations and configurations.Mean Young's modulus values for experiments are in good agreement with those obtained from simulations with the exceptions of the P3, the Ammann-Beenker and the Socolar honeycombs.In all three cases, the Young's modulus measured experimentally is higher than that predicted by modelling.On inspection of Fig. 12 it can be seen that these samples have relatively high Poisson's ratio in comparison with other patterns.As discussed by Clarke et al. [5], simulations do not directly imitate the constraints imposed on experimental samples.The necessity for rigid top and bottom plates on experimental samples applies horizontal constraints that are not present in simulations.It is reasonable to expect that these surface plates have more of a constraining effect on samples with high Poisson's ratio.This effect is demonstrated using a simplified model, with results displayed in Fig. 14.A square sample of homogeneous, isotropic material in plain strain conditions was modelled using Fenics.Dirichlet boundary conditions were applied to the top and bottom nodes simulating the constraints imposed by the plates.By varying Poison's ratio and holding force constant, displacements were used to estimate correction factors ( 0 ∕ * ) for calculating the equivalent Young's modulus of an unconstrained model ( 0 ) from the effective Young's modulus ( * ) of a constrained model with a given Poisson's ratio.Using Ammann-Beenker sample 4 as an example, and taking the maximum Poissons ratio measured in the elastic region (0.85), a correction factor of 0.83 is obtained.This brings the measured Young's modulus down from 340.8 MPa to 282.9 MPa.A more refined method is required to accurately predict correction factors for honeycombs.

Fig. 2 .
Fig. 2. A graphical depiction of how patches of aperiodic lattices are randomly sampled from a large patch.Seven levels of deflation are shown in the large patch, in reality many more are carried out.The Robinson Triangle tiling is shown as an example.

Fig. 3 .
Fig. 3. Graphic showing where 50  50 mm samples are taken from a 150  150 mm patch.The Robinson Triangle tiling is shown as an example.

Fig. 4 .
Fig. 4. Example of a simulated Ammann-Beenker sample: (a) The sample of honeycomb (-region) with the extended homogenised region (-region) and periodic boundary conditions; (b) An area of uniform mesh with discontinuous material assignment.

Fig. 5 .
Fig. 5. Manufactured samples: (a) A Penrose-P3 honeycomb; (b) An image taken from DIC software.A vertical strain gauge and horizontal box strain gauge are shown.

Fig. 6 .
Fig. 6.Images of PLA honeycombs taken from DIC software.A central patch of each pattern is shown.All samples are 50 x 50 x 50 mm.

Fig. 7 .
Fig. 7. Plot to show the maximum percentage deviation from the mean for occurrences of each cell at each orientation for 50 randomly sampled 20, 30, 40 and 50 mm patches (2 tiles, 10 rotations).Each point represents the maximum deviation from the mean for occurrences of one tile at one orientation.The Robinson triangle is shown as an example, data for other lattices can be found in supplementary data.

Fig. 8 .
Fig. 8. Plot to show the number of tile orientations found in each randomly sampled patch of Pinwheel tiling.Each point indicates that at least one sample contained the corresponding quantity of tile orientations.Larger points indicate more samples contained the corresponding quantity.

Fig. 11 .
Fig. 11.Plots showing stress/strain curves obtained from mechanical testing of PLA honeycombs.Tests on 5 differing patches are shown for each pattern.The spread of Young's modulus data obtained from simulations is represented by a shaded area.

Fig. 12 .
Fig. 12. Plots showing Poisson's ratio versus normal strain obtained from mechanical testing of PLA honeycombs.Tests on 5 differing patches are shown for each pattern.The spread of Poisson's ratio data obtained from simulations is represented by a shaded area.

Fig. 14 .
Fig. 14.Plot showing factors to correct for changes in effective Young's modulus ( * ) due to constrained samples with varying Poisson's ratio where  0 is a sample with 0 Poisson's ratio or no constraint.

Table 1
Cell area (  ) to sample area (  ) ratios for simulated samples.