Designing recipes for auxetic behaviour of 2-d lattices
DDesigning recipes for auxetic behaviour of 2-d lattices
Daniel J. Rayneau-Kirkhope , ∗ and Marcelo A. Dias , † Aalto Science Institute, School of Science, Aalto University, FI-02150 Espoo, Finland Department of Applied Physics, Aalto University, FI-02150 Espoo, Finland and Nordita, Royal Institute of Technology and Stockholm University,Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden (Dated: February 22, 2016)We present an analytical model to investigate the mechanics of 2-dimensional lattices composedof elastic beams of non-uniform cross-section. Our approach is based on reducing a lattice to asingle beam subject to the action of a set of linear and torsional springs, thus allowing the problemto be solved through a transfer matrix method. We show a non-trivial region of design space thatyields materials with auxetic properties for strains greater than that required to trigger elasticinstability. The critical loading required to make this transition from positive to negative Poisson’sratio is calculated. Furthermore, we present lattice parameters that provide direction-dependentdeformation modes offering great tailorability of the mechanical properties of finite size lattices.Not only is our analytical formulation in good agreement with the finite element simulation results,but it provides an insight into the role of the interplay between structure and elastic instability,and gives an efficient methodology to pursue questions of rational design in the field of mechanicalmetamaterials.
INTRODUCTION
Poisson’s ratio is a material constant defined as thenegative of the ratio between the transverse and longi-tudinal strains in the direction of an applied load. Ourdaily experiences frequently tell us that when a materialis compressed in a given direction, it will expand in the di-rection perpendicular to the applied load: materials con-forming to this ubiquitous behavior are characterized bypositive values of their Poisson’s ratio. In contrast, andperhaps against our common intuition, there exist exam-ples of both man-made and naturally occurring materialsthat contract (expand) in the direction perpendicular toan applied compressive (tensile) load. Such materials arecharacterised by a negative Poisson’s ratio and are oftenreferred to as auxetic materials.Since the first description of a material exhibiting anegative Poisson’s ratio [1], numerous examples and ap-plications have been reported: natural layered ceram-ics [2], fabric reinforcement[3], and low-stiffness auxeticyarns and fabrics [4] among many others. The surpris-ing behavior of auxetic materials continues to challengeour intuition and did the same to many influential physi-cists, for instance, in 1964 Richard Feynman noted that“it is reasonable that [Poission’s ratio] should be gener-ally positive, but it is not quite clear that it must be so”[5]. However, recent advancements in the field of aux-etic metamaterials have presented fresh new looks intomaterial sciences, structural design, and soft condensedmatter physics [6–10].It is often the case that auxetic properties of a ma-terial are a result of interplay between a structure’s in-ternal geometry and the constituent material properties.Such structures, where one or more mechanical propertyis dependent on the geometry of the sample’s substruc-ture (rather than purely the material composition), fall into the class of mechanical metamaterials [11]. Exam-ples of the exploitation of the internal geometry’s abil-ity to achieve novel global mechanical properties includepentamode metamaterials [12, 13], mechanical cloaks forflexural waves [14], seismic metamaterials [15] as well asauxetic mechanical metamaterials [6].Recent work centering on a square 2-dimensional lat-tice of circular holes has shown that through close controlof internal geometry, elastic instability can be utilised asa route to uni-directional, planar, and auxetic behaviour[6–8]. This example is of particular interest because itprovides a systematic way to tune macroscopic mechan-ical responses through a wide range of Poisson’s ratiosand stiffnesses [6, 7, 9]. Furthermore, the auxetic natureof the structure is “switched on” above a critical strain;before this strain the material has a positive Poisson’sratio and this allows for further tailorability in the ma-terial’s mechanical response [6]. In this 2-dimensionallattice, uniaxial compressive loads induce a short wave-length elastic deformation mode: the lattice structureundergoes a high degree of reorganization as alternatingmutually orthogonal ellipses are formed in place of thecircular voids and the structure collapses. Hence, thereconfiguration results in an auxetic response of the lat-tice. This elastic instability and the consequent breakingof internal lattice symmetry is shown in figure 1-(a). It isnoted that a square lattice made up of homogenous slen-der beams may also, in theory, present a similar bucklingmode where the lattice deforms with a short wavelengthmode, as depicted in figure 1-(d). It has been demon-strated however, that the short wavelength mode is neverpreferred because the stresses required for the onset ofthis instability are always higher than the stresses thattrigger the long wavelength modes depicted in figure 1-(c)[16, 17].In this article, we consider a square lattice made of a r X i v : . [ c ond - m a t . s o f t ] F e b component beams described by two of their thicknesses:the central section of the beam has thickness t , while theend sections of the beam have thickness t , part of a lat-tice made up of these elements is shown in figure 1-(b).We find that through manipulation of the beam thick-ness and the length of the central section considered, theshort wavelength deformation mode can be accessed ata lower critical load than the long wavelength mode. Incontrast to methods based on Bloch-wave analysis [18–21], we present an analytical model for this finite sizedlattice, including boundary effects, and show good agree-ment between our model and finite element simulations.This approach to calculating the buckling load of the lat-tice gives us an efficient method for the exploration of thedesign space. We capture a non-trivial region of param-eter space that would yield metamaterials with auxeticproperties for loading greater than that required to trig-ger elastic instability. Within this region of design space,we also calculate the critical load that would be requiredto switch the properties of the material from positive tonegative Poisson’s ratio. Furthermore, by considering alattice with different beam properties in the elements par-allel and perpendicular to the direction of applied load-ing, we explore the possibility of creating samples whichhave direction-dependent failure modes: that is, whenloaded in one direction, the material will fail with shortwavelength failure and thus would be described by a neg-ative Poisson’s ratio; however, when loaded in a direc-tion perpendicular to this, the sample would exhibit longwavelength failure and thus no auxetic properties.We now lay out the organisation of this article. Insection 2, we introduce the necessary elements that de-scribe the lattice’s unit cell for our model example. Ourapproach is based on the classical Euler-Bernoulli beam,which has been modified to incorporate non-uniform elas-tic beams. In section 3, we introduce the general method-ology, namely the transfer matrix approach. This methodallows us to describe the in-plane deflection of elasticbeams in the most general way possible by taking intoconsideration a series of attached linear and torsionalsprings as well as thickness variation. We analyticallyexpress the transfer matrices for these three forms of dis-continuities, unify this treatment in order to describe thegeneral transfer matrix across the finite direction of thelattice, and deal with the boundary conditions of the sys-tem. In section 4, we show how symmetry analysis makesthe homogenisation of the unit cell to the lattice struc-ture possible and discuss how the linear and torsionalsprings are determined from the reaction forces and mo-ments that the lattice inflict on a localised unit cell. Insection 5, we present our results which are comprisedof good agreement between the Finite Element Method(FEM) and our semi-analytic calculations. In section 6,we summarize our main results and provide the readerwith our view of the broad impact of our method. (a) (b)(c) (d) FIG. 1. [Colour online] (a) 2-dimensional periodic porous lattices;(b) Beam lattice of varying thickness; (c) Long wave length modefor a beam lattice of constant thickness; (d) Short wave lengthmode for a beam lattice of constant thickness.
BEAM THEORY AND THE UNIT CELL l l l t t L Λ κ κ κ κ τ τ τ τ F F (b)(a)
FIG. 2. [Colour online] (a) The modified square lattice underinvestigation here and the notation used. (b) A schematic of themodel used in our analytic approach, a single vertical beam in thelattice (coloured blue in (a)) is taken and the effect of the latticeon the beam is encapsulated by a series of linear and rotationalsprings.
Let us consider a 2-dimensional lattice of height L T which is infinite in the horizontal spatial dimension.Through assuming certain symmetry relationships be-tween one vertical element of length L T and its neigh-bours, we are able to establish the behaviour of the wholelattice through the analysis of one beam element. Theinfluence of the lattice on the beam we here consider isencapsulated in a distribution of elastic support providedto the beam itself – this is shown schematically in figure(2) where much of the notation used in this paper is intro-duced. The strength of this elastic support is dependenton symmetry relations that we shall consider in the fol-lowing sections. The approach in this work adopts andexpands the transfer matrix formulation of previous work[22], where such matrices were used to calculate the op-timal placement of linear springs along a Euler-Bernoullibeam.The basis for this model is the classical Euler-Bernoullibeam equation [23, 24]. Here, we consider the generalsituation where the slender beam under considerationmay have a varying cross-section, or thickness, causingits second moment of area I to be a function of the La-grange coordinates ˜ x , I = I (˜ x ). The beam of length L T is subjected to a compressive force p (˜ x ) and externalbody loads acting perpendicular to the long axis of thebeam q (˜ x ). It can be shown that the beam’s deflectionfrom its initially straight configuration, ˜ y (˜ x ) obeys thefollowing ODE:d d˜ x (cid:18) EI (˜ x ) d ˜ y (˜ x )d˜ x (cid:19) + p (˜ x ) d ˜ y (˜ x )d˜ x = q (˜ x ) , (1)where E is the Young’s modulus of the material. Weshall consider two sources of elastic support: linear andtorsional springs. The former enters the balance equationexplicitly as, q (˜ x ), and is written q (˜ x ) = − ˜ K (˜ x )˜ y (˜ x ) , (2) where ˜ K (˜ x ) is the stiffness of the elastic support. The lat-ter form of elastic support, provided by torsional springs,accounts for the applied moments on the beam given by m (˜ x ) = −T (˜ x ) ϕ (˜ x ) , (3)where T (˜ x ) is the stiffness of the field of torsional springsand ϕ (˜ x ) describes the rotation of the beam at ˜ x relativeto its initial configuration. We adopt the small angleapproximation, which reads ϕ ≈ d˜ y/ d˜ x .In this work, we introduce beam thickness discontinu-ities in a unit cell, where these are represented in Eq. (1)by a specific choice of the functional form of I (˜ x ). Theexplicit choice of I (˜ x ) to be considered here, for the tran-sition between two regions of constant thickness as shownin figure 2-(a), is given by I (˜ x ) = I + ( I − I ) Θ(˜ x − l ) , (4)where I and I are the second moments of area for thesections of length l and l , respectively, and Θ(˜ x − l ) isthe Heaviside step function (Θ = 0 for x < l and Θ = 1for x > l ).As a matter of convenience, we define the followingnon-dimensional quantities: y ≡ π ˜ y/L T , x ≡ π ˜ x/L T , f ≡ pL T / ( EI π ) and K ≡ ˜ K L T / ( EI π ). Eq. (1) thentakes a dimensionless form,d d x (cid:26)(cid:20) r I −
1) Θ (cid:18) x − πl L T (cid:19)(cid:21) d y d x (cid:27) + f ( x ) d y ( x )d x + K ( x ) y ( x ) = 0 , (5)where r I ≡ I /I METHODOLOGY
Through considerations of symmetry, the bucklingthreshold calculations of the lattice will be reduced tothat of a single vertical beam with a variable cross sectionand an appropriate set of linear and/or torsional springsplaced along its length. These springs are used to repre-sent the influence of the horizontal elements in the latticeon the vertical beam. The buckling of this single beamis analysed through a transfer matrix formulation. Thisapproach is then used to infer the elastic properties of thelattice. The transfer matrices calculated here can be de-rived from considerations of the continuity/discontinuityof the solution to Eq. (1) and its derivatives. We ded-icate the rest of this section to the derivation of thesetransfer matrices for linear springs, torsional springs and thickness discontinuity.
Linear springs:
Let us assume only point like springsupports. In general, each linear spring placed on thebeam may be thought of as having independent stiff-nesses. Therefore, in order to describe the system, aset of spring constants { κ j } , for j ∈ { , , · · · , N − } ,must be defined. For point-like linear springs, we writethe distribution K ( x ) in Eq. (5) in the following way: K ( x ) = N − (cid:88) j =1 κ j δ ( x − x j ) . (6)This discrete set of supports divides the beam into N segments, in between these discrete positions { x j } theEuler-Bernoulli Equation, Eq. (5), with K ( x ) = 0 governsthe deflection of the beam. This equation can be solved inthe regions x ∈ ( x j , x j +1 ), for any j . Hence, the generalsolution is given by: y ( x ) = A j sin (cid:104)(cid:112) f ( x − x j ) (cid:105) ++ B j cos (cid:104)(cid:112) f ( x − x j ) (cid:105) + C j ( x − x j ) + D j , (7)where we have defined the boundaries to be placed at x ≡ x N ≡ π . Considering a single spring placedat x j and integrating the Eq. (5) over a small interval around x j , it is found that,lim x → x + j y ( x ) = lim x → x − j y ( x ) , lim x → x + j y (cid:48) ( x ) = lim x → x − j y (cid:48) ( x ) , lim x → x + j y (cid:48)(cid:48) ( x ) = lim x → x − j y (cid:48)(cid:48) ( x ) , lim x → x + j y (cid:48)(cid:48)(cid:48) ( x ) − lim x → x − j y (cid:48)(cid:48)(cid:48) ( x ) + κ j y ( x j ) = 0 . (8)Defining v j ≡ ( A j , B j , C j , D j ) T , the continuity relationsshown in Eq. (8) on the piecewise solution of Eq. (7) canbe captured in the following form, v j = T lin j · v j − , (9)where the transfer matrix is defined as T lin j = κ j sin [ √ f ∆ x j ] f / + cos (cid:2) √ f ∆ x j (cid:3) κ j cos [ √ f ∆ x j ] f / − sin (cid:2) √ f ∆ x j (cid:3) κ j ∆ x j f / κ j f / sin (cid:2) √ f ∆ x j (cid:3) cos (cid:2) √ f ∆ x j (cid:3) − κ j sin [ √ f ∆ x j ] f − κ j cos [ √ f ∆ x j ] f − κ j ∆ x j f − κ j f x j , (10)where ∆ x j ≡ x j − x j − . Torsional springs:
With the addition of torsionalsprings placed along the beam at each position x j , weconsider an applied moment given by m j = − τ j y (cid:48) ( x j ),where τ j ≡ T j L T / ( πEI ). Therefore, as in the previoussection, we define a set { τ j } of torsional spring constantsfor j ∈ { , , · · · , N − } . These additional moments af-fect the boundary conditions of the Eq. (5) [23]. It canbe shown that the expressions relating function y ( x ) and its derivatives on either side of the rotational spring are:lim x → x + j y ( x ) = lim x → x − j y ( x ) , lim x → x + j y (cid:48) ( x ) = lim x → x − j y (cid:48) ( x ) , lim x → x + j y (cid:48)(cid:48) ( x ) − lim x → x − j y (cid:48)(cid:48) ( x ) + τ j y (cid:48) ( x j ) = 0 , lim x → x + j y (cid:48)(cid:48)(cid:48) ( x ) = lim x → x − j y (cid:48)(cid:48)(cid:48) ( x ) . (11)Hence, given that the solution shown in Eq. (7) is validon either side of the torsional spring, we obtain the anal-ogous transformation to Eq. (9), i.e. v j = T tor j · v j − ,thus finding that the transfer matrix for torsional springto be given by T tor j = cos (cid:2) √ f ∆ x j (cid:3) − sin (cid:2) √ f ∆ x j (cid:3) τ j cos [ √ f ∆ x j ] √ f + sin (cid:2) √ f ∆ x j (cid:3) cos (cid:2) √ f ∆ x j (cid:3) − τ j sin [ √ f ∆ x j ] √ f τ j f
00 0 1 0 − τ j cos [ √ f ∆ x j ] √ f τ j sin [ √ f ∆ x j ] √ f − τ j f + ∆ x j . (12) Thickness variation:
Using the above methodology,we now derive relationships for y ( x ) and its derivativesacross a discontinuity in beam thickness. Here we con- sider a single change in beam thickness at x j moving fromone second moment of area I ( x j − < x < x j ) to an-other second moment of area I ( x j < x < x j +1 ). Here,Eq. (5) can be solved, when K ( x ) = 0, for two differentintervals: (i) the region where I ( x ) = I , y ( x ) = A j − sin (cid:104)(cid:112) f ( x − x j − ) (cid:105) ++ B j − cos (cid:104)(cid:112) f ( x − x j − ) (cid:105) + C j − ( x − x j − ) + D j − , (13)and (ii) I ( x ) = I : y ( x ) = A j sin (cid:34)(cid:114) fr I ( x − x j ) (cid:35) ++ B j cos (cid:34)(cid:114) fr I ( x − x j ) (cid:35) + C j ( x − x j ) + D j . (14)Then, integrating Eq. (5) over a small interval around x j , we may write the following continuity equations:lim x → x + j y ( x ) = lim x → x − j y ( x ) , lim x → x + j y (cid:48) ( x ) = lim x → x − j y (cid:48) ( x ) , lim x → x + j r I y (cid:48)(cid:48) ( x ) = lim x → x − j y (cid:48)(cid:48) ( x ) , lim x → x + j r I y (cid:48)(cid:48)(cid:48) ( x ) = lim x → x − j y (cid:48)(cid:48)(cid:48) ( x ) . (15)The above system, Eqs. (15), and the solutions given inEqs. (13) and (14), allow us to write the transfer matrixfor the transition I → I . Therefore, we arrive at itsexplicit form T → j = √ r I cos (cid:2) √ f ∆ x j (cid:3) −√ r I sin (cid:2) √ f ∆ x j (cid:3) (cid:2) √ f ∆ x j (cid:3) cos (cid:2) √ f ∆ x j (cid:3) x j . (16) Notation for lattice:
We now introduce the notationthat is used for the transfer matrix across the entire finitedirection of the lattice. The effect of a given horizontalelement of the lattice at x j acting on the vertical beamunder consideration will be encapsulated by linear andtorsional springs. Thus the transfer matrices given inEqs. (10) and (12) with the appropriate spring stiffnesseswill be used. On either side of the horizontal elements,there will be a discontinuity in beam thickness in thedirection of increasing x : prior to the springs we will havea step from I → I and after the horizontal elementa second step from I → I . Thus, we introduce thenotation, T j ≡ T → j T lin j T tor j T → j +1 . (17)representing the node at x j and the closest thickness dis-continuities. It can then be seen that the entire lattice is now conveniently described by the product of the matri-ces at each node given by Eq.(17): R ≡ T → N (cid:89) j =1 T j T → N +1 , (18)where T → and T → N +1 on either side take into accountthe vertical boundaries of the lattice being in a regionof thickness t , as is depicted in figure 2. Using thisnotation, we can relate v N − to v through the expression v = R v N − . (19) Boundary conditions and buckling:
Here we estab-lish boundary conditions for the end points of the beamcorresponding to clamped-clamped ends, which equatesto y (0) = y ( π ) = 0 and y (cid:48) (0) = y (cid:48) ( π ) = 0. It is noted thatat the boundaries, the thickness of the beam is t andthus Eq. (7) with x j = 0 or x j = x N − are the relevantsolutions in the regions x ∈ ( x , x ) and x ∈ ( x N − , x N )respectively. In terms of Eq. (7) these boundary condi-tions dictate that B + D = 0 (20) A N − sin( (cid:112) f ∆ x N ) + B N − cos( (cid:112) f ∆ x N ) ++ C N − ∆ x N + D N − = 0 (21) f A + C = 0 (22) f A N − cos( (cid:112) f ∆ x N ) − f B N − sin( (cid:112) f ∆ x N ) ++ C N = 0 . (23)Using the expression given in Eq. (19), we see thatthese expressions can be rewritten in terms of v =( A , B , C , D ) T as M v = 0 , (24)where the elements of the matrix M are given by M , = M , = M , = 1 , (25) M , = (cid:112) f (26) M ,i = R ,i sin( (cid:112) f ∆ x N ) + R ,i cos( (cid:112) f ∆ x N ) ++ R ,i ∆ x N + R ,i for i ∈ [1 ,
4] (27) M ,i = R ,i (cid:112) f cos( (cid:112) f ∆ x N ) − R ,i (cid:112) f sin( (cid:112) f ∆ x N ) ++ R ,i for i ∈ [1 , , (28)and all other values are zero. The buckling load is thengiven by the minimum value of f such that det( M ) = 0. RECIPES TO BUILD LATTICES
In this section we consider how a single vertical beamelement can be used to model the complex elastic be-haviour of a 2-dimensional lattice. As shown schemati-cally in figure 2, the approach is based on the replacement x j x j +1 x j +2 x j +3 x j x j +1 x j +2 x j +3 y i y i +1 y i y i +1 (a) (b) FIG. 3. [Colour online] The symmetries assumed in the deforma-tion mode of the lattice. (a) shows part of the anti-symmetric de-formation mode where y i ( x ) = − y i +1 ( x ). The applied moments onthe ends of the horizontal beams, which are a result of the deforma-tion of the vertical beams, are of reversed handedness (representedby yellow/red circular arrows for clockwise/counterclockwise). (b)shows part of the symmetric deformation mode ( y i ( x ) = y i +1 ( x ))where end moments are of the same handedness. of the horizontal beam elements within the lattice witha set of torsional and linear springs. It is shown in figure3 that the behaviour of the horizontal beams is depen-dent on the relationship between the deflection of a givenbeam, y i ( x ), and the deflection of its neighbours y i − ( x )and y i +1 ( x ). Here we introduce the assumption that thedeformation of one vertical lattice member will be relatedto its neighbours’ through either an antisymmetric rela-tionship, y i ( x ) = − y i +1 ( x ) or a symmetric relationship, y i ( x ) = y i +1 ( x ). These symmetries are shown schemat-ically in figure 3. It is noted that this assumption, al-though reducing the space of possible deformation modesthe lattice can take, only neglects higher order modes andis justified by finite element simulations presented laterin this work.Under this assumption, the strength of the linear andtorsional springs can be calculated for both the anti-symmetric figure 3-(a), or symmetric modes, figure 3-(b).Eq. (5) can be solved for the deflection of the horizontalbeam of length L by taking the appropriate function I ( x )with zero compressive load. This allows us to establishthe relationship between applied end moment, and thegradient of the deflection of the beam, through compari-son with Eq. (3), we are then able to calculate τ . In thesymmetric regime ( y i ( x ) = y i +1 ( x )) deformations resultin no end-to-end length change of the beam and thus κ s = 0 (29)and the end moments applied to the horizontal beam areof the same handedness as shown in figure 3, consequently τ s = 2 r I L T π ( l + 2 r I l ) . (30)In the antisymmetric regime ( y i = − y i +1 ), any deflection y i ( x ) at the point of the horizontal beam will result in a (a)(b) FIG. 4. [Colour online] The deformation modes for (a) thesymmetric and (b) antisymmetric mode obtained for r L = 0 . t = 0 . l = 0 . N = 5 and r I = 0 . change in the end-to-end length of the horizontal beam.The cost of this deformation is encapsulated in the linearspring, whose spring constant, κ , can then be calculatedaccording to the geometry of the horizontal member, κ a = 12 L T t π t ( t l + 2 t l ) . (31)In this same regime, it is noted that the applied momentson the horizontal beams are of opposite handedness ateach end. Thus, τ can be shown to be τ a = 6 L T R I L π ( l + 8 r I l + 6 r I l l ( l + 2 l )) . (32)Taking these two regimes separately, the buckling loadfor the anti-symmetric and symmetric modes can be cal-culated, the minimum of which will correspond to theactive mode for a given set of parameters describing thelattice system. RESULTS
Numerical investigations were undertaken in order tovalidate the results of the single beam model proposedin this work. The numerical scheme used for this taskwas finite element method (FEM) performed using COM-SOL Multiphysics 4.4 [25]. Our simulation setup usesthe 2D Structural Mechanics module together with the (b)(a)
FIG. 5. [Colour online] First invariant of the stress tensor foundthrough FEM simulations show the regions of compression (blue orinward arrows) and tension (red or outward arrows) for the defor-mation of (a) the symmetric and (b) antisymmetric modes.
Solid Mechanics interface. In this interface we use alinear constitutive law, i.e.
Hookean elasticity, as wellas a geometrically linear model. The material proper-ties of the lattice are chosen to have Young’s modu-lus E = 170 GP a , material’s Poisson ratio ν = 0, andmass density ρ = 2329 kg/m . The studies are carriedout through a linear buckling analysis with a parametricsweep over the beam thicknesses and their relative length.Mesh refinement to confirm convergence has been under-taken.For a given lattice, our analytic method can be usedto predict the loading at which buckling will occur forboth the symmetric and antisymmetric modes. The min-imum of these two loadings will correspond to the activemode. Typical deformation patterns for the two sym-metries are shown in figure 4 where the outline of theundeformed lattice is shown in the background in black,the red beams indicates the results of a finite elementlinear buckling analysis, and the blue curve indicates thepredicted deformation of the neutral axis of the com-ponent beams from the single beam approximation pre-sented here. It is noted that antisymmetric modes alwayscorrespond to short wavelength buckling where y i ( x ) isclose to zero at the nodes, thus indicating a high en-ergy cost in the stretching of the horizontal elements ascompared to bending them. Increasingly close agreementbetween FEM and our single beam model is found withan increasing aspect ratio of the component beams. Withdecreasing r I , bending becomes more and more concen-trated in the central region of the beam in both the sym-metric and antisymmetric modes.Studies on the post-buckling behaviour of porous struc-tures have shown that the effect of the pore shape dic-tates whether symmetric or antisymmetric buckling takeplace [7]. Figure 5, resulting from FEM simulations, canbe used to further visualise the competing mechanismsfor the short and long wavelength elastic deformation:we contrast heat maps of the trace of the stress tensor(first invariant) for the two modes of deformation, (a)symmetric and (b) antisymmetric. In the antisymmetricdeformation mode, it is observed that the majority of the curvature is localised to the central region of the beamswith reduced thickness t (see figure 2). The trace of thestress shows increased magnitudes in these regions of in-creased curvature, where moments induced on the beamof thickness t can be seen though regions of compressionand tension. The thicker lattice elements (that make across centred on the nodes of the lattice) experience verylittle deformation and instead rotate almost as rigid bod-ies. On the other hand, the symmetric deformation moderesults in a shearing effect on the square voids within thelattice. This is bought about through deformation of thethicker beam element (thickness t ), hence the symmet-ric mode occurs when the beam thicknesses t and t arecomparable. Therefore, long wavelength resulting fromthe compression-tension pattern in figure 5-(a) can be in-terpreted as the system promoting deformation close tothe joints of the thicker elements (the cross at the nodes).This effect can be seeing in a similar system, where thechoice of lattice that results in long wavelength has star-shaped voids [7, 26]. It is also observed here that in thelong wavelength mode the trace of the stress is dominatedby the small wavelength sub-oscillation (the long wave-length failure without the effect of the horizontal beamswould yield blue on the inside of the curvature and redon the outside) signifying the importance of deformationon the smaller length-scale, even when considering longwavelength failure [17].For a given lattice ( L , r L , and t fixed), the dependenceof the critical loading of the lattice and the ratio r I can beobtained. This dependence is shown in figure for threelattices with t = 0 . L = 0 . r L = 0 . , . .
5. It is found that for a given r L there is a transitionfrom the antisymmetric to the symmetric mode beingthe active deformation mode with increasing r I . Goodagreement with FEM simulation is found for both f min and the transition from short to long deformation mode.Figure 7 shows a further exploration of the design spaceof this system: for a given t and L , the values of r L and r I can be varied and the active deformation mode es-tablished. Through setting the parameters that wouldresult in antisymmetric deformation modes (short wave-length), a region of parameter space is shown in figure7. It is noted that for this range of design parametersthe post-buckling regime would exhibit auxetic materialproperties. FEM simulations are shown in figure 7 oneither side of the phase transition, close agreement is ob-tained.Finally, we explore the possibility of creating latticeswith direction-dependent responses to a fixed verticalloading. Here allow the value of r I to vary between thehorizontal and vertical lattice elements, we denote thesetwo parameters r h and r v respectively. We restrict our in-vestigation to the case where r L is equal for all the beamcomponents. For a given pair of values, r (1) I and r (2) I ,two simulations can be performed: r v = r (1) I , r h = r (2) I r I f c Antisymmetric modeSymmetric modeFinite element r I f c Antisymmetric modeSymmetric modeFinite element r I f c Antisymmetric modeSymmetric modeFinite element
FIG. 6.
Failure loading for 3 particular lattices: symmetric andanti-symmetric deformation pattern shown against results of finiteelement simulations. Solid lines represents anti-symmetric modewhich would result in auxetic post-buckling behaviour. Dashed linerepresents symmetic failure mode. Results are shown for latticeswith t = 0 . L = 0 . N = 5 and r L = 0 . , . and r v = r (2) I , r h = r (1) I . Figure 8 shows the region( r (1) I , r (2) I ) of phase space containing three distinct re-gions: (a) long wavelength for the two possible orienta-tions; (b) and (c) one orientation with short wavelengthmode and the other long; (d) short wavelength deforma-tion being present in both possible orientations. It is r I r L Short wavelength failure, FEMLong wavelength failure, FEM
Symmetric modeAntisymmetric mode
FIG. 7.
Phase diagram showing short vs long wavelength failuresin the r I - r L parameter space. Other parameters describing theframe are set as t = 0 . N = 5, and L = 0 .
1. The plots infigure are slices through this parameter space with set r L . d ) ( b )( c ) ( a ) ( a )( b )( c )( d ) r (1) I r (2) I r h = r (1) I r h = r (2) I r v = r (2) I r v = r (1) I FIG. 8. [Colour online] The response of a system with the verti-cal and horizontal beams discontinuity in second moment of areadescribed by parameters r v and r h respectively. All other param-eters are equal in the two sets of beams. For each pair of values( r (1) I , r (2) I ), two calculations are performed: r v = r (1) I , r h = r (2) I and r v = r (2) I , r h = r (1) I . Region (a) depicts the area of parameterspace where long wavelength instability is the active mode in bothorientations, (d) shows where short wavelength mode will be activefor both, while (b) and (c) show where the two orientations willgive different modes. noted that the effects of the boundaries parallel to thedirection of loading are relatively short ranged [6], thuswe hypothesized that although the theory presented hereis for lattices that are infinite in the horizontal direction(see figure 2), this dual response will be present for finitelattices as well. It is therefore likely that a square sampleconstructed with parameters taken from region (b) or (c)has auxetic properties when compressive load is appliedparallel to one set of component beams, but it behavesas a material with positive Poisson’s ratio when a load isapplied perpendicular to these beams. CONCLUSION
We have presented an analytical model, based on asingle beam approximation, that captures key featuresof the elastic instability of a 2-dimensional square latticewith non uniform component beams. We have shown ex-cellent agreement with FEM simulations for the criticalloading and deformation mode for a wide range of latticeparameters. Furthermore, we have utilised the efficientmethodology presented here to explore a large parame-ter space and have uncovered a non-trivial phase spacethat differentiates long from short wavelength deforma-tion. This phase space describes when the 2-dimensionallattice will behave as an auxetic material and, within thisregion of design space, the critical loading required tomake the transition from positive to negative Poisson’sratio. The model has then been used to make predic-tions about materials with lattice whose parameters dif-fer from horizontal to vertical beams. This has uncoveredthe possibility of direction-dependent auxetic structuralproperties: we have shown that there exist regions ofparameter space in which a square lattice, subjected tocompressive loads parallel to one set of beams, would be-have as a material with positive Poisson’s ratio but whenperpendicular to these beams, the structure may insteadpresent an auxetic behaviour. The model can also beused to predict the tuning of the load of transition frompositive to negative Poisson’s ratio. The critical loadingcould also be varied depending on its directionality withrespect to the beams on the sample.
ACKNOWLEDGEMENTS
The authors would like to acknowledge Mikko Alavafor insightful discussions regarding this work. D.R-K.thanks the support of the Academy of Finland throughits Centres of Excellence Programme (2014-2019). ∗ marcelo.dias@aalto.fi † daniel.rayneau-kirkhope@aalto.fi[1] R. Lakes, Foam Structures with a Negative Poisson’s Ra-tio., Science (New York, N.Y.) 235 (4792) (1987) 1038–40.[2] F. Song, J. Zhou, X. Xu, Y. Xu, Y. Bai, Effect of a Neg-ative Poisson Ratio in the Tension of Ceramics, PhysicalReview Letters 100 (24) (2008) 245502.[3] Z. Ge, H. Hu, Innovative three-dimensional fabric struc-ture with negative Poisson’s ratio for composite reinforce-ment, Textile Research Journal 83 (5) (2012) 543–550.[4] J. R. Wright, M. K. Burns, E. James, M. R. Sloan, K. E.Evans, On the design and characterisation of low-stiffnessauxetic yarns and fabrics, Textile Research Journal 82 (7)(2012) 645–654. [5] R. P. Feynman, R. B. Leighton, M. Sands, Feynman lec-tures on physics, Vol. II, California Institute of Technol-ogy, 1964.[6] K. Bertoldi, P. M. Reis, S. Willshaw, T. Mullin, NegativePoisson’s ratio behavior induced by an elastic instability.,Advanced materials (Deerfield Beach, Fla.) 22 (3) (2010)361–6.[7] J. T. B. Overvelde, S. Shan, K. Bertoldi, Compactionthrough buckling in 2D periodic, soft and porous struc-tures: effect of pore shape., Advanced materials (Deer-field Beach, Fla.) 24 (17) (2012) 2337–42.[8] J. Shim, C. Perdigou, E. R. Chen, K. Bertoldi, P. M.Reis, Buckling-induced encapsulation of structured elas-tic shells under pressure., Proceedings of the NationalAcademy of Sciences of the United States of America109 (16) (2012) 5978–83.[9] B. Florijn, C. Coulais, M. van Hecke, Programmablemechanical metamaterials, Phys. Rev. Lett. 113 (2014)175503.[10] C. Coulais, J. T. B. Overvelde, L. a. Lubbers, K. Bertoldi,M. van Hecke, Discontinuous Buckling of Wide Beamsand Metabeams, Physical Review Letters 115 (4) (2015)044301.[11] M. Kadic, T. B¨uckmann, R. Schittny, M. Wegener, Meta-materials beyond electromagnetism., Reports on progressin physics. Physical Society (Great Britain) 76 (2013)126501.[12] G. W. Milton, A. V. Cherkaev, Which elasticisty tensorsare realizable?, J. Eng. Mater. Technol. 117 (483).[13] M. Kadic, T. Buckmann, N. Stenger, M. Thiel, M. We-gener, On the practicability of pentamode mechanicalmetamaterials, Appl. Phys. Lett. 100 (191901).[14] N. Stenger, M. Wilhelm, M. Wegener, Experimentson elastic cloaking in thin plates, Phys. Rev. Lett.108 (014301).[15] S. Brule, E. Javelaud, S. Enoch, S. Guenneau, Experi-ments on seismic metamaterials: Molding surface waves,Phys. Rev. Lett. 112 (133901).[16] N. Ohno, D. Okumura, T. Niikawa, Long-wave bucklingof elastic square honeycombs subject to in-plane biax-ial compression, International Journal of Mechanical Sci-ences 46 (11) (2004) 1697–1713.[17] B. Haghpanah, J. Papadopoulos, D. Mousanezhad,H. Nayeb-Hashemi, A. Vaziri, Buckling of regular, chi-ral and hierarchical honeycombs under a general macro-scopic stress state, Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences 470(2014) 20130856.[18] R. Hutchinson, N. Fleck, The structural performance ofthe periodic truss, Journal of the Mechanics and Physicsof Solids 54 (4) (2006) 756 – 782.[19] M. S. Elsayed, D. Pasini, Analysis of the elastostatic spe-cific stiffness of 2d stretching-dominated lattice materi-als, Mechanics of Materials 42 (7) (2010) 709 – 725.[20] A. Vigliotti, D. Pasini, Stiffness and strength of tridimen-sional periodic lattices, Computer Methods in AppliedMechanics and Engineering 229232 (2012) 27 – 43.[21] K. Bertoldi, Stability of periodic porous structures, in:Extremely Deformable Structures, Springer, 2015, pp.157–177.[22] D. Rayneau-Kirkhope, R. Farr, K. Ding, Y. Mao, Bifur-cations in the optimal elastic foundation for a bucklingcolumn, Physics Letters A 375 (2) (2010) 67 – 72.0