A computational method for determining the linear elastic properties of 2D aperiodic lattice structures

Abstract

This paper develops a framework for determining the linear elastic properties of non-periodic lattice structures. An element-based material assignment methodology is implemented that facilitates the generation and analyses of arbitrary patterns on a structured mesh. An adapted numerical homogenization strategy features the inclusion of a homogenized region in the neighbourhood of the domain boundary that validates the implementation of periodic boundary conditions for an arbitrary finite patch of a periodic or non-periodic lattice structure. To demonstrate the method, the linear elastic properties of an aperiodic lattice pattern based on the Penrose (P3) pattern is evaluated. Such a structure exhibits order without translational symmetry and consequently lacks a repeating unit cell. The isotropic performance of the aperiodic lattice structure is investigated and compared to that of the well-known square periodic lattice. The framework opens the door to the investigation and analyses of other novel cellular structures which are not based on a repeating unit cell. Additive manufacturing facilitates the physical realization of such lattice structures, presenting them as viable alternatives to conventional periodic structures in the aerospace and bio-engineering industries.

Keywords

Penrose tiling aperiodic order lattice homogenization non-periodic

Introduction

Technological advances in additive manufacturing (3D printing) have relaxed the manufacturing constraints on complex geometries, motivating exploration of the space of material design for next-generation lightweight structures.¹ In response, material scientists and engineers are investigating novel ‘designed-for-purpose’ materials to broaden the scope of achievable material performance.^2–4 For example, materials are being designed to replicate bone structures that exhibit remarkable robustness to arbitrary load paths, with capability to efficiently sustain multiple loading modes.⁵ Such materials derive their mechanical properties from their meso-architectures as well as the chemical composition of the base material they are made from,^6,7 and one area of research is in the use of cellular lattice structures.^7–9 For example, cellular lattices have been integrated into the design of bespoke orthopaedic implants.^10,11

Cellular structures are typically defined according to periodic repetition of unit cells.⁹ The structural behaviour of periodic lattice structures is well understood,^12,13 and by investigating such structures over multiple scales, homogenization techniques have been used to determine their effective (homogenized) mechanical properties.^14,15 These effective properties determine the macroscopic behaviour of the periodic lattice,¹⁶ and explicit analytical expressions have been derived for common periodic lattices, such as the square and hexagonal lattices, which describe their elastic behaviour when used as thin-walled honeycomb structures.^12,17 However, a fundamental limitation associated with the use of periodic lattice structures is the requirement for over-designing to accommodate the existence of undesirable mechanical anisotropy.

In natural materials, such as bone, there exist pseudo-random structures which cannot be represented in repeatable fashion by periodic cellular structures.¹⁸ To avoid the limitations of periodicity in designed structures, and taking advantage of the capabilities of additive manufacturing, alternative non-periodic cellular geometries have been explored, such as conformal lattices where meso-architectures of the lattice are locally orientated according to principal stress directions,^19,20 or other geometries resulting from structural topological optimization.^21,22 But simulating and analyzing the resultant behaviour of these geometries is difficult due to their underlying non-periodic geometry. The in-plane mechanical properties of cellular structures are often determined by the process of homogenization which involves determining effective elastic properties of a homogeneous equivalent phase. Most approaches to executing homogenization can be classified as analytical^23,24 numerical^15,16 or experimental.^25,26 Many numerical strategies use asymptotic expansion homogenization to determine effective mechanical properties such as the effective stiffness tensor, ultimate strength or yield strength. For periodic structures, this involves defining a representative volume element (often a repeating unit cell) and applying periodic boundary constraints on opposite cell boundaries on a structured mesh boundary.^27,28 Strain deformation analyses are then executed to yield the stress distributions from which effective properties are determined.²⁹ Even where the structures are non-periodic or are composites containing randomized inclusions of particles, efficient strategies have been developed to determine their properties by imposing periodicity in the neighbourhood of the boundaries representative volume element.^30–32 Despite the reliability associated with the numerical homogenization technique described, the efficiency of the process relies significantly on the existence of a repeating unit cell which lies at the heart of lattice periodicity. Non-periodic lattice structures do not feature a repeating unit cell, and as a result a direct application of these techniques is inappropriate.

This paper presents a general method for simulating and analysing the linear elastic properties of non-periodic two-dimensional lattice structures. It builds on work by Imediegwu et al.^33,34 where a numerical framework was presented in which lattice unit cells are defined on a structured mesh by assigning material properties spatially within the domain of a discretized unit cell according to a unique vector of radii parameters. The approach boycotted the challenges of manual geometry meshing for complex lattice structures and presented a means to evaluate parametric sweeps of geometric variables with automated mesh generation. We extend this methodology for non-periodic lattice structures by implementing an element-based material assignment on an arbitrarily chosen patch of a non-periodic pattern. The approach satisfies the requirement for identical mesh nodes on opposite boundaries of the numerical domain for the correct implementation of periodic boundary conditions. To satisfy the requirement for geometric periodicity in the neighbourhood of the boundary of the numerical domain, a homogenized region is introduced around the patch of interest.³⁵ Asymptotic expansion homogenization is iteratively applied to determine effective properties of the extended patch, with property reassignment to the homogenized region, until convergence of the mechanical properties are attained.

As a demonstration of the method, an analysis of aperiodically-ordered lattice structures is presented. Aperiodic crystallographic structures are ordered but not symmetric under any translation.^36,37 They form a relatively new field of inquiry, modelling the structure of quasicrystals, which were first discovered in the 1980s.^38,39 As a consequence, the concept of aperiodicity is often taken congruous to the idea of quasi-periodicity.^40,41 Due to their unusual structure, metallic quasicrystals show remarkable micromechanics, more akin to ceramics than metals and they can be considerably more isotropic than periodic crystals – giving nearly uniform mechanical properties in all directions.⁴² As such, mimicking the form of quasicrystals at the mesoscopic scale offers an interesting possibility for cellular structures, and their linear elastic properties are worthy of investigation, but difficult to simulate and analyze using methods that take advantage of periodically repeating geometry. Randomized techniques for stochastic metamaterials are also unsuitable since aperiodic tessellations possess long-range-order despite lacking local tessellation rules.³⁶ Recent investigations of mechanical properties of quasi-periodic metamaterials demonstrate that quasi-periodic structures offer intriguing possibilities for mechanical applications.⁴³ Beli et al.⁴⁰ demonstrated that quasicrystals with higher order symmetries allow for high equivalent stiffness characteristics. They also show that these intriguing aperiodic structures exhibit more uniform strain energy distributions in comparison to other periodic configurations. Somera et al.⁴¹ demonstrated that the elementrary pattern shape plays a key role to the overall mechanical behaviour of these aperiodic lattices. The mechanical properties of an aperiodically-ordered lattice structure was recently investigated based on a lattice spring model.⁵ However, we deem the patch for the aperiodic lattice structure simulated by Kim et al.⁵ as non-representative - sufficient scale separation between the size of the cells and the computational domain is a strict requirement for numerical homogenization. To ensure sufficient cells of the aperiodic lattice are captured within the patch domain of our work, we randomly vary the centre of the numerical domain relative to the centre of the aperiodic pattern to generate a stochastic sampling of material properties for a range of pattern resolutions. We ensure that scale separation is sufficiently large to enforce invariance in mechanical properties for the resolution of choice through mesh and pattern resolution sensitivity studies. The derived framework supports the determination of directional Young’s modulus via rotational analyses, providing insight into direction-dependent and pattern resolution-dependent mechanical behaviour of a lattice structure based on the Penrose (P3) tiling. The result is the development of a robust numerical toolkit for generating and investigating aperiodically-ordered lattice structures. The family of aperiodic lattice structures are yet to be characterized and there are potentially an infinite number of aperiodic tilings. As a consequence, it is essential that the technique introduced in this work is rapid – suited to a wide range of aperiodic lattice geometries at arbitrary orientations. The framework is also applicable to periodic lattice structures, offering a means to validate its implementation against classical unit cell homogenization.

The next section introduces an aperiodic pattern and its generation. Details concerning the numerical strategies implemented follow leading up to the validation of the methodology. The section culminates in computational implementation of the method. Section 3 shows the performance of the exemplar lattice structure against Hashin-Shtrikman bounds. It also shows rotational analyses results in comparison to the well-known square lattice structure. Finally, Section 4 provides a summary of the work presented and presents a hint of future research directions.

Methods

The proposed framework is best introduced in the context of an example. For this purpose, we focus on non-periodic lattice structures based on the P3 aperiodic pattern, first introduced by Penrose.⁴⁴ The next sections present the pattern along with methods for generating lattice structures based on the pattern. We then introduce the framework for computationally determining the mechanical properties of these aperiodic lattices and validate the framework via consideration of the well-understood square-based lattice structure.

The Penrose (P3) pattern

The P3 pattern is an example of an aperiodic tiling. An aperiodic tiling is defined as one built up of cells (prototiles) in an ordered manner but without periodicity. Aperiodic order is an emerging mathematical field of inquiry that governs the crystalline structure of quasi-crystals. As a consequence, aperiodic structures are often associated with the term ‘quasi-periodicity’.^41,43 The P3 pattern is composed of two rhombi prototiles, with acute angles $72^{°}$ and $36^{°}$ (denoted as prototile type A and type B respectively). Though it is possible to attain a periodic assembly of these two rhombi prototiles, the P3 pattern features a non-periodic assembly of both prototiles. Indeed, such an array is strictly aperiodic, lacking translational symmetry yet capable of completely tiling the entire $R^{2}$ -plane. By ‘aperiodicity’, we refer to a pattern which is ordered but lacks translational symmetry. This differs significantly from a ‘non-periodic’ assembly such as a composite structure with generalized voronoi-shaped particles,³² particulate-reinforced materials³⁰ or cork-based agglomerates,³¹ for example. It is impossible to assemble (as possible with periodic repeating patterns) the P3 pattern using the prototiles by implementing local rules. In general, this property of aperiodic tilings is consistent with the crystallographic concept of long-range-order such that local tessellations correlate intrinsically with detached tessellations in the aperiodic system.³⁶ As a consequence, aperiodic patterns are often created by two main methods.⁴⁵ A ‘cut and project’ method generates aperiodic patterns by taking hyperslices through high-dimensional periodic patterns, whilst ensuring the hyperslice does not cut through any vertices of the high-dimensional periodic structure. For example, the P3 pattern is a two-dimensional slice of a five-dimensional hypercube.³⁷ A second approach, called the substitution method, is implemented in this work. In the substitution method, an aperiodic pattern is generated by substituting each prototile by scaled-down prototiles based on unique scaling ratios as illustrated in Figure 1. For conciseness, we refer to this scaling-down process as a ‘deflation’ going forward. An aperiodic pattern generator algorithmically deflates prototiles over several deflation levels to generate the aperiodic pattern of choice. We refer to the first prototile generated at algorithm commencement as the prime prototile, which has dimensions much larger than that of the patch of interest in order to fully cover the patch with tiling after deflation is complete. A stopping criteria for the deflation process is often predefined to meet pattern resolution requirements.

Figure 1.

Penrose-P3 deflation.

Regarding the P3 patterns of this work, pattern generation commences with a prime rhombus prototile of type A, spanning a ‘large extent’ of the $R^{2}$ plane relative to the domain of a comparably ‘much smaller’ finite patch of interest. However, the choice of a prime prototile of type A or B is inconsequential to the P3 pattern generation provided sufficient deflation processes are executed. Section 2.9 investigates the mechanical property fluctuations associated with the arbitrary choice of a finite patch of the generated P3 pattern at different pattern resolutions.

Pattern generation

This section presents the methodology for the pattern generator algorithm. The algorithm has also been generalized to the creation of periodic patterns used for comparison. Consider a square patch of side, $L$ within $R^{2}$ with centre point coordinates $(x, y) = (\frac{L}{2}, \frac{L}{2})$ such that the origin $O (0, 0)$ lies at the bottom-left corner of the patch. The algorithm commences generation by creating the prime prototile, which may be of arbitrary centre point coordinates. For simplicity, we choose the centre point coinciding with the patch centre, $(x^{*}, y^{*}) = (x, y)$ and orientation, $θ^{*} = 0$ as shown in Figure 1(A,I), that is prototile type A, generation I. The side length of the prime prototile, $l^{*}$ is chosen such that

l^{*} > > L

(1)

The algorithm assigns the prime prototile to a repository of created prototiles. For ease of replication and given $(x^{*}, y^{*})$ , $l^{*}$ and $θ^{*}$ , we present the specific deflation rules for the P3 tiling shown in Appendix Tables A1 and A2 where $l_{n}$ and $(c_{x}, c_{y})$ are the side length and centre point coordinates for all current generation level prototiles, $ψ = \frac{π}{10}$ for the P3 pattern and $l_{n + 1}$ is the scaled side length for next level prototiles given by

l_{n + 1} = \frac{l_{n}}{φ}

(2)

The scaling parameter, $φ = 2 \cos \frac{π}{5} \approx 1.618$ is often referred to as the Golden ratio. These rules govern the creation (positioning and orientation) of all type A and B prototiles of the next generation level. The repository is then updated by substituting each prototile of the previous generation with the newly created prototiles. This substitution process continues until predefined pattern resolution requirements are satisfied (see Section 2.9). Figure 1 shows a geometrical parallel to the equations of Tables A1 or A2 as type A and type B prototiles of generation I are deflated over two generations II and III. In this work, we define the resolution of any P3 pattern, $r_{p}$ as:

r_{p} = \frac{a_{A}}{L^{2}}

(3)

where $a_{A}$ is the area of prototile type A at deflation termination and $L^{2}$ is the area of the patch. To create the aperiodic lattice patch, prototiles in the repository at deflation termination which lie outside the patch domain are eliminated. Finally, governed by the Sutherland-Hodgman algorithm,⁴⁶ prototiles lying across the patch boundary are clipped, leaving only relevant portions within the patch domain. Figure 7(a) shows an arbitrary clipped section of a P3 pattern at deflation termination (Figure 2).

Figure 2.

Arbitrary patch of a Penrose-P3 pattern, clipped by a square boundary.

Element-based material assignment

To generate the corresponding P3 lattice structure for onward finite element analyses, the domain enclosed by the aperiodic patch is discretized into a structured mesh and an element-based material assignment ensues. This methodology is motivated by the following finite element-based considerations:

A geometry-based representation requires a physical boundary for the application of boundary conditions and/or forces for the finite element analyses, significantly altering the effective mechanical properties evaluated.

The element-based approach affords better control of mesh discretization which facilitates the application of periodic boundary conditions.

The element-based approach facilitates automatic mesh generation for arbitrary lattice patterns of this work, presenting a technique to rapidly simulate the geometry and mechanical properties of a large family of novel aperiodic structures at arbitrary orientations. It is worthy of note to emphasize that the material assignment methodology has its limitations, particularly for low-density lattices. Mesh resolution must be increased with decreasing minimum thickness to sustain the reliability of the strategy. In this work, mesh discretization is a function of the minimum thickness such that at least 6 elements always lie across the minimum thickness for any given pattern.

The details of the element-based material assignment methodology is given as follows:

Consider the discretized section of the patch shown in Figure 3(a) with enlarged nodes for emphasis. For illustrative purposes, only a single cell is shown. The generated mesh is a structured triangular mesh. A structured mesh facilitates the implementation of periodic boundary conditions in the neighbourhood of the patch boundaries. Generated by the pattern generator of Section 2.2, green points in Figure 3(a) denote prototile vertex points of an illustrative type I rhombus prototile, with its sides shown in red dashed lines. Each mesh node must be mapped to a binary space based on its coordinate location relative to the line segments joining any pair of prototile vertex points. For each line segment, the black nodes (in the viscinity of that line segment) satisfy the following two conditions:

They lie within a distance $\frac{t}{2}$ to the line segment, where $t$ is the predefined thickness of struts of the lattice structure.

The orthogonal projection of their coordinate positions on the line segment belongs to the set [0,1].

Mathematically, the coordinate positioning of the black nodes are such that

\begin{matrix} Ψ (x, \frac{t}{2}) \leq 0 \end{matrix}

(4)

and

0 < Φ (x) < 1

(5)

where

Ψ (x, \frac{t}{2}) = \frac{| (b - a) \times (a - c) |}{| b - a |} - \frac{t}{2}

and

Φ (x) = \frac{〈 (b - a), (c - a) 〉}{〈 (b - a), (b - a) 〉}

For any line segment, the line end vertex points have coordinates $a$ and $b$ while $c$ is the coordinate of any particular node. The $operation 〈, 〉$ denotes the vector inner product operation. The white nodes do not satisfy the conditions of equations (4) and (5). Finally, the discretized elements are assigned material properties based on each element’s node fraction, $n_{f}$ given as:

n_{f} = \frac{N_{m}}{3}

(6)

where $N_{m}$ is the number of each triangular element’s black nodes. Node fraction values for each element would vary between 0 and 1 as shown in Figure 7(c). Local Young’s modulus on each element is determined by

E_{e} = n_{f} E_{s} + (1 - n_{f}) E_{v}

(7)

where $E_{s}$ and $E_{v}$ represents the Young’s modulus of the isotropic base material and void (represented as a very weak isotropic material) respectively. As a consequence, elements with intermediate node fraction values, lying across the solid-void boundary, would have intermediate Young’s modulus values. Figure 3(b) illustrates the generation of the aperiodic lattice structure on a coarse structured mesh discretization. The structure-void interphase smoothens with mesh refinement as will be demonstrated for aperiodic patches of this work in subsequent sections. This approach to modelling lattice structures is inspired by the computational efficiency associated with modelling complex geometries and the ease of executing parametric sweeps through pattern design variables such as the prime prototile orientation, $θ^{*}$ , without the need for explicit geometry re-meshing.

Figure 3.

Illustration of element-based material assignment for aperiodic lattice structure generation: (a) node classification and (b) element classification based on node fraction, $n_{f}$ .

Numerical analyses assumptions

The description of the numerical analyses of the framework commences with a note on the background assumptions governing its implementation. Considerations have been limited to two-dimensional linear elastic theory so that the general relationship between the Cauchy stress tensor, $σ$ and strain tensor, $ε$ at material point is given by

σ = λ_{s} Tr (ε) I + 2 G_{s} ε

(8)

where $Tr$ is the trace function of a matrix and $I$ is the identity matrix. The bulk modulus, $K_{s}$ and the Lamé constants $λ_{s}$ and $G_{s}$ are given by

\begin{matrix} K_{s} & = \frac{E_{s}}{3 (1 - 2 ν_{s})} \\ λ_{s} & = \frac{E_{s} ν_{s}}{(1 + ν_{s}) (1 - 2 ν_{s})} \\ G_{s} & = \frac{E_{s}}{2 (1 + ν_{s})} \end{matrix}

(9)

respectively, where $ν_{s}$ is Poisson’s Ratio of the base isotropic material.

The relationship between strain and displacement, $u$ are given by

ε_{ij} = \frac{1}{2} (u_{i, j} + u_{j, i})

(10)

Equations (8) to (10) define the linear elastic isotropic constitutive equations governing our problem formulations. Static equilibrium mechanics is also assumed so that

σ_{ij, j} + F_{i} = 0

(11)

where $F_{i}$ represents all external forces acting on the patch domain, $Ω$ . All numerical analyses is modelled according to plane strain hypothesis.

Mesh convergence studies

Mesh studies are critical to the successful implementation of the method presented in this work. A course mesh will present challenges with effective representation of the pattern and will increase the region of intermediate densities between the material and void regions. However, high mesh resolution, though giving converged mechanical properties, increases computational cost exponentially. An efficient meshing strategy that supports accurate numerical simulations but minimized computational cost is critical. Mesh convergence studies of this work involved discretizing the SRVE domain for a fixed P3 and SQR lattice patterns over increasing number of elements along each dimension. The number of elements were incremented from 100 to 800 elements in steps of 20 elements. Figure 4 shows the convergence plots for the components of the elasticity matrices of the square and P3 lattices of this work. All numerical analyses have been executed at 600 elements on each dimension to ensure convergence of the mechanical properties.

Figure 4.

Mesh convergence, normalized elasticity components, $E_{xx}$ versus mesh refinement, $n_{e}$ for both square and P3 lattices.

Homogenization with boundary phase inclusion

The effective mechanical properties of aperiodic structures can be determined by an adaptation to asymptotic expansion homogenization. Aperiodic lattice structures are not comprised of a repeating unit cell and hence it appears counter-intuitive to suggest that their mechanical properties can be determined by a similar procedure. In this work, asymptotic expansion homogenization is implemented to derive properties of aperiodic structures. This is supported by a theoretically-consistent adaptation that permits the use of periodic boundary conditions. The theory of homogenization is well documented in literature^47,48 and will not be reintroduced in this work. However, the adaptation of the homogenization technique as applied to aperiodic structures is presented as follows:

Consider the aperiodic lattice structure generated by the element-based material assignment of Section 2.3. For ease of reference, the patch, which is a statistical representative volume element (SRVE), is referred to as the $α$ -region going forward as shown in Figure 5(a). Observe that the mesh resolution minimizes the elements with intermediate properties so that lattice struts are well defined at the scale of the patch. Also observe that the pattern resolution constitutes a representative patch of the P3 aperiodic pattern – a fact which is proven in Section 2.9. Periodic boundary conditions are applied to opposite edges of the SRVE – an association of the degrees of freedom of nodes on opposite edges. The homogenization technique as described by Hassani and Hinton¹⁵ is then executed to yield an initial estimate of $α$ -region’s elasticity matrix, $E^{α}$ . We refer to this initial evaluation of effective elastic property of the $α$ -region as iteration zero. The orientation of the SRVE as shown in Figure 5 is preserved for all simulation results of this work and will be omitted hereafter for conciseness.

Figure 5.

Simulation patches for aperiodic lattice structure: (a) $α$ -region of $υ_{f} = 0.5$ for iteration 0 and (b) extended ( $α$ -region and $β$ -region) domain for all other iterations.

Next, a second region is introduced that bounds the $α$ -region. The inclusion of this region, referred to as the $β$ -region inclusion, constitutes the homogenized phase. The effective properties determined in iteration 0 are assigned to the elements of the $β$ -region. The effective properties of the extended domain (both $α$ and $β$ -regions) are then recomputed and reassigned to the $β$ -region in the next iteration. A minimum thickness of the $β$ -region is determined in Section 2.8. This property recomputation of the extended domain and reassignment to the $β$ -region is executed over several iterations until convergence of the trace of the elasticity matrix, $Tr (E)$ . Observe that homogenized properties in the neighbourhood of the edges of the extended domain permit the correct application of periodic boundary conditions for all homogenization processes on the extended domain. Figures 6 shows the convergence plots for the relevant components of the SRVE’s elasticity matrix. Figures 7 and 8 show the displacement and stress magnitude plots at convergence of material properties. The computed elasticity matrix after homogenization is given as

E^{α + β} = [\begin{matrix} E_{11} & E_{12} & E_{16} \\ E_{22} & E_{26} \\ sym & E_{66} \end{matrix}] = [\begin{matrix} 423.8 & 252.5 & 0 \\ 425.2 & 0 \\ sym & 87.2 \end{matrix}] MPa

(12)

for $E_{s} = 2000 MPa$ , $E_{υ} = 10^{- 4} MPa$ and $ν_{s} = 0.3$ . Observe that the elasticity matrix obtained is symmetric, consistent with physical requirements. Extension-shear coupling is also non-existent as components $E_{16} = E_{26} = 0$ , even for an arbitrary finite patch of the P3 lattice structure. Mesh refinement is converged for all simulations of this work (Figure 9).

Figure 6.

Elasticity matrix convergence plots: (a) $E_{11}$ and $E_{22}$ component convergence, (b) $E_{33}$ component convergence and (c) $E_{12}$ component convergence.

Figure 7.

Displacement magnitude, $u$ (deformation scale = 0.1:1): (a) $u$ due to $ε_{11} = (1, 0, 0)^{⊤}$ , (b) $u$ due to $ε_{22} = (0, 1, 0)^{⊤}$ and (c) $u$ due to $ε_{12} = (0, 0, 1)^{⊤}$ .

Figure 8.

Stress magnitude, $σ$ (deformation scale = 0.1:1): (a) $σ$ due to $ε_{11} = (1, 0, 0)^{⊤}$ , (b) $σ$ due to $ε_{22} = (0, 1, 0)^{⊤}$ and (c) $σ$ due to $ε_{12} = (0, 0, 1)^{⊤}$ .

Figure 9.

Square periodic structure for $v_{f} = 0.5$ : (a) representative volume element (RVE) and (b) extended domain (PATCH).

Validation

To validate the framework implemented in this work, the mechanical properties of an arbitrary SRVE of a known periodic structure is determined. The results obtained are compared with the classical homogenization technique on a representative unit cell of the same square periodic structure. Material properties as documented by Gibson and Ashby¹² also serve as a benchmark as applicable to thin-walled honeycombs. Figure 10 shows good correlation of the classical homogenization technique with the representative volume element (RVE) and the adapted homogenization technique on the extended domain (PATCH) for the relevant components of the elasticity matrix and for $0 \leq υ_{f} \leq 1$ . The Ashby model for the square periodic honeycombs (ASHBY) agrees with both homogenization techniques for thin-walled honeycombs when volume fraction, $υ_{f} < 0.3$ as documented.

Figure 10.

Validation of the adapted AEH technique.

Minimum $β$ -region thickness, $t_{β}$

The homogenized phase is essential to the methodology because it permits the use of periodic boundary conditions in the neighbourhood of the boundaries of the extended domain. However, its inclusion enlarges the computational domain. Too thin, and the $β$ -region inclusion will introduce errors in the evaluated material properties. Too thick, and its inclusion significantly degenerates computational efficiency for a fixed mesh refinement. To determine a minimum thickness of the $β$ -region for the chosen pattern resolution of this work, it was steadily increased until convergence was attained in the computed mechanical properties. Figure 11 demonstrates that $t_{β} = 0.11 L$ meets the requirements for accuracy and efficiency for the $β$ -region thickness, regardless of the lattice structure simulated. We have ensured pattern equivalence by enforcing equal pattern resolution and volume fraction, denoted $r_{p}$ and $v_{f}$ respectively.

Figure 11.

$β$ -region thickness convergence plots: (a) Penrose-P3 structure (P3) and (b) square-periodic structure (SQR).

Effect of variations in finite SRVE centre and pattern resolution

To investigate the sensitivity of evaluated mechanical properties to pattern resolution, we execute a random sampling of the SRVE centre relative to the prime prototile centre point. For each pattern resolution as shown in Figure 12(a), the Young’s modulus in the x-direction was evaluated for 1000 randomly positioned SRVEs. Each boxplots depicts the variation in property evaluated. The plots present conventional statistical information associated to box plotting and the three pattern resolutions are shown in Figure 12(c). It can be observed that a resolution of $r_{p} = 0.001$ yielded less than $2 %$ variation within the data’s interquartile range. On the other hand, the coarsest P3 lattice of the study, with resolution $r_{p} = 0.0071$ yielded just over $15 %$ variation within the data’s interquartile range. For the same random sampling, Poisson’s ratio variations were similar to those of the Young’s modulus, though much less scattered about the median Poisson’s ratio ( $< 10$ %) within the data’s interquartile range (see Figure 12(b)). The circles in Figure 11(b) represent outlier Poisson’s ratio values for 1 and 4 SRVE samples at resolutions, $r_{p} = 0.0027$ and $0.0071$ respectively. The three resolutions depicted in Figure 12 correspond to the last three sequential deflation levels of the P3 lattice structures of this work. All results correspond to mechanical properties of the P3 lattice structure at pattern resolution, $r_{p} = 0.001$ .

Figure 12.

Boxplots illustrating the effect of variations in the finite SRVE centre for different pattern resolutions of the P3 lattice structure: (a) effective Young’s modulus in x-direction, (b) effective Poisson’s ratio and (c) pattern resolutions.

Computational implementation

All numerical formulations of this work have been implemented with FEniCS,^49,50 an open source finite element solver based on the Python language. SRVE domain for typical simulations of this work featured mesh refinement comprising $4 \times 600^{2}$ Continuous Garlekin (Lagrange) triangular elements generated by the built-in mesh generator within FeniCS. This automated mesh generation eliminates the requirement for explicitly creating mesh files for each pattern, especially for analyses involving parametric sweeps of design variables. The pattern generation and element-based material assignment on an Intel(R) Core(TIM) i7-5600U CPU at 2.60 GHz system with two processors is evaluated in approximately 4 min for the pattern resolutions of our work. On this system, each iteration of the adapted homogenization procedure is resolved in approximately 40 s, facilitated by parallelization of the three strain deformation analyses. Converged properties are obtained typically within less than 20 iterations. Solutions yielded by the direct solvers are visualized with Paraview – an open source mesh data visualization application.

Results and discussion

Performance against Hashin-Shtrikman bounds

Having established validation for the algorithm used to derive the elasticity matrix of aperiodic lattice structures, some mechanical properties of the Penrose P3 were investigated. The plots in Figure 12 show the variation of bulk and shear moduli of the P3 and SQR lattice patterns. The Hashin-Shtrikman bounds show the bulk and shear moduli limits for any composite made up of a mechanical mixture of two or more isotropic and homogeneous elastic phases.⁵¹ In this work, our patterns have been modelled as composites predominantly composed of two materials, a base isotropic material and void (represented as a very weak isotropic material). A very small volume fraction ( $< 0.07$ ) consists of intermediate material properties as a consequence of the element-based material assignment methodology which creates the pattern. However, these intermediate materials also prevent singularities at the material/void interphases within our numerical framework. Observe that the P3 and SQR lattice structures have similar effective bulk moduli, particularly in the range where $υ_{f} > 0.5$ , for uniform thickness variation. However, the P3 lattice structure proffers higher shear moduli for all volume fractions, though more significantly in the range $υ_{f} > 0.5$ .

Rotational analyses

To measure the degree of anisotropy associated with patterns in this work, rotational analyses were conducted. By varying the orientation of the prime rhombus in increments of $1^{°}$ for $0 < θ^{*} < 360$ , lattice structures were generated for the P3 and SQR patterns and the directional Young’s moduli and Poisson’s Ratios were determined.⁴ Such a parametric sweep of the orientation variable was facilitated by the element-based material assignment methodology, boycotting the challenges of geometry re-meshing. Figures 13 and 14 clearly depict the high dependence of the properties of the SQR lattice structure to orientation. The SQR lattice is much stiffer when struts are aligned in the directions of applied strain but much weaker otherwise. This also accounts for its comparatively lower shear modulus in Figure 15(b). The SQR lattice pattern’s Poisson’s ratio also changes significantly with orientation. This is not the case for the P3 lattice structure. It exhibits near-isotropic mechanical properties in Young’s modulus and Poisson’s ratio, demonstrating remarkable insensitivity to orientation variations. Even though the P3 lattice pattern offers comparatively lower Young’s modulus for a given volume fraction, its near-isotropic nature confirms speculation about a significant correlation between isotropic material behaviour and the existence of higher rotational symmetries associated with aperiodic patterns. To further investigate the degree of isotropy in comparison to the square lattice, we leverage the elastic anisotropy measure for 2D crystals as presented by Li et al.⁵² where given the elasticity matrix of any given 2D lattice system, a universal elastic anisotropy index, $A^{SU}$ can be derived given by:

\begin{matrix} A^{SU} = & ([\frac{1}{4} (C_{11} + C_{22} + 2 C_{12}) (S_{11} + S_{22} + 2 S_{12}) - {1]}^{2} \\ + 2 \frac{1}{16} (C_{11} + C_{22} - 2 C_{12} + 4 C_{66}) \\ \times (S_{11} + S_{22} - 2 S_{12} + S_{66}) - {1]}^{2})^{\frac{1}{2}} \end{matrix}

(13)

where the $C_{xx}$ terms are the components of the stiffness matrix and the $S_{xx}$ terms are the components of the compliance matrix for the lattice system. A perfectly isotropic lattice system would have $A^{SU} = 0$ , with increasing $A^{SU}$ values representing increasing anisotropy.

Figure 13.

Directional Young’s modulus, $\frac{E}{E_{s}}$ .

Figure 14.

Directional Poisson’s ratio, $ν$ .

Figure 15.

Effective bulk and shear moduli against Hashin-Shtrikman bounds: (a) normalized bulk modulus, $\frac{K}{K_{s}}$ and (b) normalized shear modulus, $\frac{G}{G_{s}}$ .

The universal elastic anisotropy index for the square and P3 lattice systems of this work are $A_{s}^{SU} = 1.2952$ and $A_{p}^{SU} = 0.000231$ respectively for $υ_{f} = 0.50$ . This confirms ‘near-isotropic’ elastic properties for the P3 in comparison to the square lattice.

Conclusion

Summary

A robust finite element framework has been developed and validated for determining the effective properties of non-periodic lattice structures. The framework is demonstrated by evaluating the mechanical properties of an aperiodically-ordered lattice structure based on the Penrose (P3) tiling, facilitated by an element-based material assignment methodology. An algorithm for designing P3 patterns was created following specific substitution rules, from which lattice structures were developed. The framework was validated by its application to the square periodic lattice structure. It also agrees with the Gibson-Ashby model for square periodic thin-walled honeycombs. The near-isotropic nature of the P3 lattice structure to pattern orientation was established, even for an arbitrary finite SRVE. The methodology implemented is robust and supports the rapid investigation of a wide range of novel aperiodic structures at arbitrary orientations. The framework is suited to periodic cellular structures as well, forming the basis for validation by comparison to unit cell homogenization.

Future directions

This study primarily documents a technique for determining the properties of aperiodic cellular structures, an example of which is the P3 aperiodic lattice structure. There exists a myriad of other aperiodic patterns from which lattice structures can be created and investigated. Infact, many aperiodic lattice structures generalize to three dimensions. This will be the focus of future work – an exploration of the material space populated by aperiodically-ordered lattice structures, leading up to the creation of a library of aperiodic lattice structures with their mechanical properties. We are aware of hexagonal and triangular lattice structures which also exhibit isotropic material behaviour.¹² However, the topology of aperiodic lattice structures integrate the element of aesthetics to possible structural design. Moreover, preliminary experimental investigation suggest astonishing non-linear response of aperiodic lattice structures to large deformation compared to periodic lattice structures. Consequently, computational validation of the non-linear elastic behaviour of aperiodic lattice structures is a worthy research direction. Ultimately, the integration of aperiodic lattice structures towards optimal structural design remains a key drive in the research of these emerging alternatives to periodic cellular structures within the aerospace and bio-engineering industries.

Footnotes

Appendix: Penrose P3 substitution rules

Table A2.

Type B substitution rules.

Type	Centre x-coordinate	Centre y-coordinate	Orientation
A	$c_{x} - \cos θ (\frac{l_{n}}{2} \cos ψ) + \sin θ (\frac{l_{n}}{2} \sin ψ)$	$c_{y} - \cos θ (\frac{l_{n}}{2} \sin ψ) - \sin θ (\frac{l_{n}}{2} \cos ψ)$	$θ - \frac{3 π}{5}$
A	$c_{x} + \cos θ (\frac{l_{n}}{2} \cos ψ) + \sin θ (\frac{l_{n}}{2} \sin ψ)$	$c_{y} - \cos θ (\frac{l_{n}}{2} \sin ψ) + \sin θ (\frac{l_{n}}{2} \cos ψ)$	$θ + \frac{3 π}{5}$
B	$c_{x} + \cos θ (l_{n + 1} \sin ψ \cos ψ) - \sin θ (\sin ψ (l_{n} - l_{n + 1} \sin ψ))$	$c_{y} + \cos θ (\sin ψ (l_{n} - l_{n + 1} \sin ψ)) + \sin θ (l_{n + 1} \sin ψ \cos ψ)$	$θ - \frac{3 π}{5}$
B	$c_{x} - \cos θ (l_{n + 1} \sin ψ \cos ψ) - \sin θ (\sin ψ (l_{n} - l_{n + 1} \sin ψ))$	$c_{y} + \cos θ (\sin ψ (l_{n} - l_{n + 1} \sin ψ)) - \sin θ (l_{n + 1} \sin ψ \cos ψ)$	$θ + \frac{3 π}{5}$

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the EPSRC for funding this research (Grant Ref: EP/V047108/1).

ORCID iD

Chikwesiri Imediegwu

References

Wang

Zhang

Daynes

, et al. Design of graded lattice structure with optimized mesostructures for additive manufacturing. Mater Des 2018; 142: 114–123.

Panetta

Zhou

Malomo

, et al. Elastic textures for additive fabrication. ACM Trans Graph 2015; 34(4): 1–135.

Schumacher

Bickel

Rys

, et al. Microstructures to control elasticity in 3d printing. ACM Trans Graph 2015; 34(4): 1–13.

Schumacher

Marschner

Gross

, et al. Mechanical characterization of structured sheet materials. Compos B 2018; 37(4): 1–2.

Kim

Libonati

, et al. Designing tough isotropic structural composite using computation, 3D printing and testing. ACM Trans Graph 2019; 167: 736–745.

Benedetti

Plessis

Ritchie

, et al. Architected cellular materials: a review on their mechanical properties towards fatigue-tolerant design and fabrication. Mater Sci Eng R 2021; 144: 100606.

Maconachie

Leary

Lozanovski

, et al. Slm lattice structures: properties, performance, applications and challenges. Mater Des 2019; 183: 108137.

Zhu

Chen

Skouras

, et al. Two-scale topology optimization with microstructures. ACM Trans Graph 2017; 36(5): 1.

Lumpe

Stankovic

. Exploring the property space of periodic cellular structures based on crystal networks. Proc Natl Acad Sci USA 2021; 118(7): e2003504118.

10.

Sigmund

. On the optimality of bone microstructure. In: IUTAM symposium on synthesis in bio solid mechanics, Copenhagen, Denmark, 24–27 May 1998, 1999, vol. 69, pp.221–234. Cham: Springer.

11.

Colabella

Cisilino

Fachinotti

, et al. Multiscale design of elastic solids with biomimetic cancellous bone cellular microstructures. Struct Multidiscip Optim 2019; 60: 639–661.

12.

Gibson

Ashby

. Cellular solids: structure and properties, chapter 4, 2nd ed. Cambridge: Cambridge University Press, 1997, pp.93–174.

13.

Chen

Skouras

Zhu

, et al. Computational discovery of extremal microstructure families. Sci Adv 2018; 4(1): eaao7005.

14.

Arabnejad

Pasini

. Mechanical properties of lattice materials via asymptotic homogenization and comparison with alternative homogenization methods. Int J Mech Sci 2013; 77: 249–262.

15.

Hassani

Hinton

. A review of homogenization and topology optimization i -homogenization theory for media with periodic structure. Comput Struct 1998; 69: 707–717.

16.

Hassani

Hinton

. A review of homogenization and topology opimization ii -analytical and numerical solution of homogenization equations. Comput and Struct 1998; 69: 719–738.

17.

Jiang

Yang

. Advanced honeycomb designs for improving mechanical properties: a review. Compos B 2021; 227: 109393.

18.

Sigmund

. Topology optimization of multi-scale structures: a review. Struct Multidisc Optim 2021; 63(3): 1455–1480.

19.

Liao

Dai

, et al. Anisotropic design and optimization of conformal gradient lattice structures. Comput Aided Des 2020; 119: 102787.

20.

Wang

Gao

. Design and optimization of conforming lattice structures. IEEE Trans Vis Comput Graph 2021; 27(1): 43–56.

21.

Bendsøe

Sigmund

. Topology optimization: theory, methods and applications, chapter 1, 2nd ed. Berlin, Germany: Springer-Verlag, 2003, pp.1–68.

22.

Zhang

Zhou

. Feature-driven method for structural optimization. vol. 1, chapter 6, 1st ed. Amsterdam, the Netherlands: Elsevier, 2020, pp.157–240.

23.

Wang

McDowell

. In-plane stiffness and yield strength of periodic metal honeycombs. J Eng Mater Technol 2004; 126(2): 137–156.

24.

Masters

Evans

. Models for the elastic deformation of honeycombs. Compos Struct 1996; 35: 403–422.

25.

Vigliotti

Pasini

. Stiffness and strength of tridimensional periodic lattices. Comput Methods Appl Mech Eng 2012; 229-232: 27–43.

26.

Moat

Muyupa

Imediegwu

, et al. Compressive behaviour of cellular structures with aperiodic order. Results Mater 2022; 15: 100293.

27.

Rodrigues

Coelho

Fernandes

. Multiscale modelling of bone tissue –remodelling and application to scaffold design. Comput Methods Appl Sci 2011; 20: 15–33.

28.

Yera

Rossi

Méndez

, et al. Topology design of 2D and 3D elastic material microarchitectures with crystalsymmetries displaying isotropic properties close to their theoritical limits. Appl Mater Today 2020; 18: 100456

29.

Abad

Arabnejad

Pasini

. Fatigue design of lattice materials via computational mechanics: application to lattices with smooth transitions in cell geometry. Int J Fatigue 2013; 47(42): 126–136.

30.

Gentieu

Catapano

Jumel

, et al. Computational modelling of particulate-reinforced materials up to high volume fractions: linear elastic homogenisation. J Mater Des Appl 2019; 233: 1101–1116.

31.

Delucia

Catapano

Montemurro

, et al. Determination of the effective thermoelastic properties of cork-based agglomerates. J Reinf Plast Compos 2019; 38: 760–776.

32.

Delucia

Catapano

Montemurro

, et al. A stochastic approach for predicting the temperature-dependent elastic properties of cork-based composites. Mech Mater 2020; 145: 103399.

33.

Imediegwu

Murphy

Hewson

, et al. Multiscale structural optimization towards three-dimensional printable structures. Struct Multidiscip Optim 2019; 60: 513–525.

34.

Imediegwu

Murphy

Hewson

, et al. Multiscale thermal and thermostructural optimisation of three-dimensional lattice structural. Struct Multidiscip Optim 2022; 65: 13.

35.

Colera

Kim

. Asymptotic expansion homogenization analysis using two-phase representative volume element for non-periodic composite materials. Multiscale Sci Eng 2019; 1: 130–140.

36.

Baake

Grimm

. Aperiodic order: a mathematical invitation, chapter 1, 1st ed. University Printing House, Cambridge: Cambridge University Press, 2013, pp.1–58.

37.

Baake

Grimm

. Aperiodic order: crystallography and almost periodicity, chapter 1, 2nd ed. University Printing House, Cambridge: Cambridge University Press, 2017, pp.1–35.

38.

Shechtman

Blech

Gratias

, et al. Metallic phase with long-range orientational order and no translational symmetry. Phys Rev Lett 1984; 53: 1951–1953.

39.

Levine

Steinhardt

. Quasicrystals: a new class of ordered structures. Phys Rev Lett 1984; 53(26): 2477–2480.

40.

Beli

Rosa

MIN

De Marqui

, et al. Mechanics and dynamics of two-dimensional quasicrystalline composites. Extreme Mech Lett 2021; 44: 101220.

41.

Somera

Poncelet

Auffray

, et al. Quasi-periodic lattices: pattern matters too. Scr Mater 2022; 209: 114378.

42.

Dubois

. Useful quasicrystals, chapter 1–2, 1st ed. Singapore, Singapore: World Scientific Publishing Co. Pte. Ltd., 2005, pp.1–128.

43.

Wang

Sigmund

. Quasiperiodic mechanical metamaterials with extreme isotropic stiffness. Extreme Mech Lett 2020; 34: 100596.

44.

Penrose

. Pentaplexity: a class of non-periodic tilings of the plane. Math Intell 1979; 2: 32–37.

45.

Senechal

. Quasicrystals and geometry, chapter 1–4, 1st ed. Cambridge: Cambridge University Press, 1995, pp.1–121.

46.

Jungang

. An advanced Sutherland-Hodgman algorithm, Journal of Xi'an University of Posts and Telecommunications 2013; 18(3): 80–82.

47.

Pinho-da Cruz

Oliveira

Teixeira-Dias

. Asymtotic homogenisation in linear elasticity. part I: mathematical formulation and finite element modeling. Comput Mater Sci 2009; 45: 1073–1080.

48.

Oliveira

Pinho-da Cruz

Teixeira-Dias

. Asymptotic homogenisation in linear elasticity. part II: finite element procedures and multiscale applications. Comput Mater Sci 2009; 45: 1081–1096.

49.

Alnaes

Blechta

Hake

, et al. The Fenics project version 1.5. Arch Numer Softw 2015; 3(100): 9–23.

50.

Logg

Mardal

Wells

. Automated solution of differential equations by the finite element method. chapter 1. Cham: Springer, 2012, pp.1–69.

51.

Hashin

Shtrikman

. A variational approach to the elastic behavior of multiphase materials. J Mech Phys Solids 1963; 11(42): 127–140.

52.

Shao

Gao

, et al. Elastic anisotropy measure for two-dimensional crystals. Extreme Mech Lett 2020; 34: 100615.

A computational method for determining the linear elastic properties of 2D aperiodic lattice structures

Abstract

Keywords

Introduction

Methods

The Penrose (P3) pattern

Pattern generation

Element-based material assignment

Numerical analyses assumptions

Mesh convergence studies

Homogenization with boundary phase inclusion

Validation

Minimum β -region thickness, t β

Effect of variations in finite SRVE centre and pattern resolution

Computational implementation

Results and discussion

Performance against Hashin-Shtrikman bounds

Rotational analyses

Conclusion

Summary

Future directions

Footnotes

Appendix: Penrose P3 substitution rules

Declaration of conflicting interests

Funding

ORCID iD

References

Minimum $β$ -region thickness, $t_{β}$