Analytic analysis of auxetic metamaterials through analogy with rigid link systems

TL;DR

This paper introduces an analytical approach to study auxetic metamaterials by drawing analogies with rigid link systems, providing accurate predictions of elastic instability and post-buckling behavior to aid design.

Contribution

It presents a novel analytical methodology linking lattice mechanisms with mechanical behavior, improving understanding and design of auxetic metamaterials.

Findings

Excellent agreement with finite element simulations

Accurate prediction of elastic instability onset

Insight into lattice mechanisms and mechanical responses

Abstract

Recent progress in advanced additive manufacturing techniques has stimulated the growth of the field of mechanical metamaterials. One area particular interest in this subject is the creation of auxetic material properties through elastic instability. This paper focuses on a novel methodology in the analysis of auxetic metamaterials through analogy with rigid link lattice systems. Our analytic methodology gives extremely good agreement with finite element simulations for both the onset of elastic instability and post-buckling behaviour including Poisson's ratio. The insight into the relationships between mechanisms within lattices and their mechanical behaviour has the potential to guide the rational design of lattice based metamaterials.

Figures

Equations

$$U$$

$${\eta }_{i,j}=-{\eta }_{i+1,j}.$$

$${\theta }_{i,0}=-{\theta }_{i,1},$$

$${\theta }_{i,N_y+1}=-{\theta }_{i,Ny}.$$

$${\theta }_{i,j}={\theta }_{i+1,j},$$

$${\theta }_{i,j}=-{\theta }_{i+1,j}.$$

$${\eta }_{i,j}=0\forall i,j.$$

$$(2{\kappa }_1+\tau -Fl){\theta }_{n,m}-{\kappa }_1({\theta }_{n,m+1}+{\theta }_{n,m-1})+({\delta }_{m,1}+{\delta }_{m,Ny})2{\kappa }_1=0$$

$$\mathbf{A}\Theta =0$$

$$F_{\text{min}}=\frac{2{\kappa }_1+\tau -2{\kappa }_1\cos \left(\frac{\pi }{N_y+1}\right)}{l},$$

$${\Theta }_i=A\sin \left(\frac{i\pi }{n+1}\right).$$

$${\theta }_{i,j}=\pm \theta$$

$${\eta }_{i,j}=\gamma (N_x,\tau ,{\kappa }_2){\theta }_{i,j}.$$

$$U=\Omega {\theta }^2-F_{i}l{N}_x{N}_y(1-\cos \theta )$$

$$F_{\text{min}}=\frac{2\Omega }{N_x{N}_{y}l}.$$

$${\kappa }_1={\kappa }_2=\frac{E{I}_{\kappa }}{l_c}$$

$$\tau =\frac{2E{I}_{\tau }}{l_c},$$

$$U\approx ({\alpha }_i-F{\beta }_i)B^2+{\zeta }_{i}F{B}^4,$$

$$B=\{\begin{array}{ll}0& \text{mbox}forF<F_{\text{mbox}min}\\ \sqrt{\frac{{\alpha }_i-F{\beta }_i}{2{\zeta }_{i}F}}& \text{mbox}forF>F_{\text{mbox}min}\end{array}$$

Full text

1Aalto Science Institute, Aalto University, FI-02150 Espoo, Finland

2Department of Applied Physics, Aalto University, FI-02150 Espoo, Finland

3Department of Mathematics, Massachusetts Institute of Technology, USA

4School of Mathematics and Statistics, University of St Andrews, St Andrews KY16 9SS, Scotland, UK

5Department of Physics and Astronomy, James Madison University, Harrisonburg, VA 22807, USA

1 Introduction

With the close control of a materials internal geometry, comes the potential to design a particular response to load into the architecture of a material. Such structures, where geometry rather than material properties govern the macroscopic response of a solid, are termed “mechanical metamaterials” [1, 2]. Perhaps the most well-known class of these kind of materials are those with an auxetic reponse to external load [3, 4]. Such structures, when compressed (stretched), contract (expand) in the direction perpendicular to the applied load.

Auxetic behaviour has been observed in layered ceramics [5], re-entrant foams [3], origami [6] and kirigami [7, 8] geometries, structures with rotating elements [10], dimpled sheets [9] and other carefully designed architectures. Some of these geometries are examples of an emerging paradigm whereby the use of elastic instability is viewed as a route to new functionality rather than a mode of failure [11, 12]. This concept has been notably used in the creation of materials endowed with negative poisson ratio, where spontaneous symmetry breaking associated with buckling creates an ordered collapse and an auxetic response beyond some threshold value of loading [13, 14, 15]. Alongside auxetic response, metamaterials exhibiting frustrated mechanics [16] and hysteretic [17] behaviours have been designed and fabricated using similar buckling induced motifs. More broadly, the utilisation of elastic instabilities in soft materials resulting in novel mechanical behaviour has been bought under the umbrella term “designer matter” [18]. This emerging field is notable for embracing truly interdisciplinary work, encompassing fields including engineering, physics, architecture and increasingly applied mathematics.

In this paper, we establish a fundamental connection between the buckling-induced auxetic behaviour of a lattice made up of a soft material and the mechanisms (“floppy modes”) designed into a lattice of rigid links. We import the knowledge of the mechanisms of the rigid link lattices and apply it to analytically establish the buckling and postbuckling behaviour of complex continuum lattice systems. In order to facilitate this analogy, we utilise a specific void shape allowing the straightforward calculation of the appropriate values for the torsional springs in the system. We note however, that the underlying mechanisms within the continuum lattice are relatively insensitive to the nature of the voids [19], thus we hypothesise that the rigid-link analysis proposed here elucidates the fundamental mechanisms present within a wide range of elastic instabilities in lattice based materials. In this paper, we choose to focus on the auxetic properties of lattice based materials, however, the methodology presented here can be applied to more complex lattice problems, including mechanical hysteresis behaviour and/or frustrated mechanics [16, 17]. Such behaviour can be studied through the introduction of more complex mechanisms which permit further degrees of freedom and/or constraints within the rigid-link lattice. It is also noted that while good agreement between applied load and resultant deformation is observed here, the model is unsuitable to model displacement controlled relationships (with the exception of the “far-from-thershold” regime) due the the rigid link nature of the links considered here.

This paper focuses on analytic solutions for the rigid link lattice system for the onset of elastic instability and its associated eigenmode, the postbuckling stiffness and associated properties such as the Poisson’s ratio of the system (found to be negative for certain regions of parameter space). We present close agreement for all of these properties with an analogous continuum system investigated through finite element methods. We propose that the use of rigid link methods can be a powerful tool in the analysis of elastic systems providing insight into the fundamental mechanisms/modes present within the lattice and reducing the computational power required for such investigations.

2 Lattice stability

The deformation of the rigid link lattice presented in Fig. (1, left) is parameterised by a set of angles $\{\theta_{i,j}\}$ denoting the rotation of the initially vertical elements, while another set of angles $\{\eta_{i,j}\}$ is used to represent the rotations of the initially horizontal elements. The indices $i$ and $j$ take integer values and denote the position of the element: element $(i,j)$ has its centre initially positioned at $(iL,jL)$ where $L$ is the length of the links. There are no horizontal element on the upper and lower boundaries ($j=0,N_{y}$). The energy of any deformation described by the parameters $\{\theta_{i,j},\eta_{i,j}\}$ can be calculated as,

[TABLE]

where $\kappa_{1}$, $\kappa_{2}$ and $\tau$ are the stiffnesses of the torsional springs in the system (springs with strength $\kappa_{1}$ ($\kappa_{2}$) penalise rotation of initially vertical (horizontal) elements relative to their neighbours while $\tau$ springs penalise relative rotation of the vertical and horizontal elements whose centres remain coincident). The term proportional to $\kappa_{1}$ contains the Kronecker delta function $\delta_{i,j}$ to account for our choice of boundary constraints. The last term in Eq. (1) is the external work done on the system by the force $F_{i}$ (at each $i$-th vertical element). Due to the connectivity of the lattice and fixed length of elements, we impose a set of constraints. Working to first order in $\{\theta_{i,j},\eta_{i,j}\}$, it can be shown that for a deformation to be compatible with the connectivity of the lattice, the distance between the locations of $\eta_{i,j}$ and $\eta_{i+1,j}$ takes a constant value for all $j$. Furthermore it is required that

[TABLE]

We set the boundary conditions of the system to be,

[TABLE]

2.1 Symmetry relations

Following some of the authors previous work [21], it can be assumed that the energy of the system is minimised when one of two symmetry relations are assumed:

[TABLE]

These modes will be referred to as translationally symmetric and mirror symmetric modes, respectively. In the following subsections, we establish the buckling behaviour of the system subject to these two possible symmetry relations.

2.1.1 Translational symmetry

Utilising the symmetry relationships presented above, the energy of a given deformation can be greatly simplified. In the case of translational symmetry ($\theta_{i,j}=\theta_{i+1,j}$) we see that energy minimisation with respect to $\eta_{i,j}$ enforces that

[TABLE]

Working to first order in $\theta_{i,j}$, direct calculation shows that the minimum energy configuration exists when

[TABLE]

is satisfied for all values of $n$ and $m$. This requirement can be rewritten in matrix form,

[TABLE]

where $\Theta=(\theta_{i,1},\theta_{i,2},...,\theta_{n,N_{y}})^{T}$. It is noted that in matrix form the first two terms of Eq. 8 create a tridiagonal symmetric Toeplitz matrix, the remaining term mean that the full expression for $\mathbf{A}$ deviates slightly from this form. Buckling of the system into a mode with translational symmetry will occur if the loading on the system is sufficient to create (atleast) one zero eigenvalue of $\mathbf{A}$. For suitably large values of $N_{y}$, neglecting the last term in Eq. (8) yields a good approximation to the system (for physically relevant parameters here, $N_{y}>5$ is sufficient). This approximation allows for analytic calculation of the eigenvalues/vectors. We find that buckling of the system into a mode with translational symmetry will occur if the loading on the system exceeds the threshold

[TABLE]

and the associated mode is found to be

[TABLE]

For small values of $N_{y}$, the eigenmode of this system can be obtained numerically.

2.1.2 Mirror symmetry

The second symmetry between neighbouring columns we consider here is that of mirror symmetry ($\theta_{i,j}=-\theta_{i+1,j}$). Due to the restriction that the mid-point of vertical and horizontal bars are coincident throughout the deformation process, it can be shown that,

[TABLE]

for some value of $\theta$ and that the sign of any rotation in the lattice is the opposite of its nearest neighbours. Furthermore through the minimisation of energy with respect to $\eta$, we can also derive the an expression for the rotation of initially vertical elements:

[TABLE]

The energy of the whole system for a given deformation can then be expressed as

[TABLE]

for some $\Omega$ which depends on parameters describing the lattice. Working to first order in $\theta$, we can thus establish that the minimum energy configuration corresponds to non-zero values of $\theta$ (buckled configurations), provided $F$ is above

[TABLE]

3 Rigid link as a continuum approximation

If we now consider a continuum lattice structure, as shown in figure 1 (right), the slender beams within the framework serve as hinge points when the lattice deforms beyond the buckling threshold. The resistance to bending of these slender beams is easily obtainable and, thus, we are able to calculate the effective values of $\kappa_{1}$, $\kappa_{2}$ and $\tau$ of the continuum lattice. From standard beam theory [20], it can be found that the appropriate stiffnesses of these springs are given by

[TABLE]

where $E$ is the Young’s Modulus of the material, $I_{\kappa}$ and $I_{\tau}$ are the second moment of area of the slender elements making up the $\kappa_{1}$, $\kappa_{2}$ and $\tau$ springs and $l_{c}$ is the length of the slender elements. Therefore, using the expressions for buckling of the rigid link lattice found in the previous section, we can predict the buckling load of the more complex continuum lattice system.

3.1 Linear stability

We note that the minimum value of the two buckling loads of the lattice (corresponding translational or mirror symmetry) will be the physically relevant mode. Thus, we are able to make predictions about the regions of the parameter space describing the continuum lattice and how they will yield modes that correspond to either of the two symmetries mentioned in the previous section. For a set values of $N_{x}=10,N_{y}=9$, $l_{c}=3.75$mm and $L=20$mm, the thicknesses of the slender beam elements are varied through the dimensionless parameters $a$ and $b$ (see figure 1 (right)). This in turn sets values of $\kappa$ and $\tau$ in the rigid link lattice through Eqs. (16 & 17).

Predictions of which symmetry mode is present within the deformed structure are then found through the rigid link methodology. The results are shown in a phase diagram, for the space $(a,b)$, where a boundary between the translationally symmetric mode (long wavelength) and mirror symmetric mode (short wavelength) is analytically drawn (see figure 2). Good agreement between these predictions and the modes computed in finite element work is observed. Furthermore, we can predict the buckling load of the continuum structure through substituting Eqs. (16 & 17) into Eqs. (10 & 15). The agreement in buckling load between the rigid link method and the finite element simulations is shown in figure 3, where increasingly better quantitative match is found for more slender beam elements. It is also noted that the linear buckling modes, as predicted in Eqs. (11, 12 & 13) agrees well with the modes found in the linear analysis through finite element studies, as shown in figure 4.

3.2 Postbuckling

The postbuckling of the rigid lattice can also be described analytically. Beyond the buckling threshold, due to the incompressibility of the links, the system remains in a configuration closely approximated by the modes predicted by the linear analysis. The magnitude of these modes can be predicted through energy considerations. Substituting Eq. (11) or Eqs. (12 & 13) into Eq. (1), for the translationally symmetric or mirror symmetric mode respectively, we find that in both cases the energy of the system can be expressed as,

[TABLE]

where $B$ characterises the magnitude deformation present within the system (in the case of a translationally symmetric mode, $B=A$ from Eq. (11), while for the mode with mirror symmetry, $B=\theta$ from Eq. (12)), $\alpha_{i},\beta_{i}$ and $\zeta_{i}$ are constants that take values dependent on the mode being investigated. It can then be shown that the minimum energy configuration is realised as

[TABLE]

Thus, considering Eqs. (16 and 17), we make predictions about the nature of the postbuckling behaviour of the continuum lattice presented in figure 1. These predictions, alongside the postbuckling behaviour found through finite element simulations are shown in figure 5, where it is noted that the rigid link analysis appears as a limit to which the continuum lattice converges in the limit of decreasing magnitudes of imperfections in the system. In figure 1, the value of $F_{\mbox{\tiny min}}$ used in Eq. (19) has been taken from the finite element simulations to decrease the error (these errors can be seen in figure 3). Nonetheless, the functional form of the graph shows good agreement between the two methods.

3.3 Auxetic metamaterials

The antisymmetric mode of this system is strongly associated with auxetic behaviour of the lattice [13, 14, 15], where the Poisson’s ratio of the structure is a function of the loading parameter [13]. We stress that in the limit of large deformations, in many cases, this dependence has a well defined limit (as investigated in [13]). By utilising the rigid link system, we are able to obtain estimates for this limit. For the geometry investigated here, it is found that, for large strains, localisation of deformations occur close to the boundaries. In figure 6 we plot the minimum value of the Poisson’s ratio observed in FEM simulations (observed immediately before the localisation of deformation) and the Poisson’s ratio predicted the rigid link analysis for various values of $a$ and $b$ and for the parameters of the system are given in the caption of figure 4.

4 Discussion

In this paper, we have elucidated the fundamental mechanisms behind the auxetic behaviour of a broad class lattice based metamaterials, and presented a new methodology in the analysis of such materials. We have applied symmetry arguments to simplify the problem before finding fully analytic solutions for the behaviour of the rigid link lattice. We have shown excellent agreement between the rigid link mechanism results and the mechanical response of buckling lattices made up of a soft isotropic material, encapsulating both the linear stability studies and post-buckling behaviour. While in this work we have focused on the auxetic properties of the lattices, the methodology is well suited to the analysis of materials exhibiting frustrated mechanics and hysteretic behaviour. We hypothesise that the methodology presented here can be a useful tool in the design of structures with the mechanical response programmed into the geometry of the material.

5 Appendix: FEM methods

Finite element studies were undertaken using COMSOL Multiphysics 5.2 [22]. Both the linear and stationary studies were preformed on the same mesh. The mesh density varied a depending on the aspect ratio of the members considered, however approximately 200 - 1500 mesh elements were used per unit cell. Mesh refinement studies were undertaken to check for convergence of results. In all simulations, a Young’s Modulus of 170MPa and Poisson’s ratio of 0.25 has been used.

We used both linear buckling analyses for figure 3 and stationary studies for figure 5. In both cases, the boundary conditions considered here were such that the upper and lower boundaries of the lattice were permitted to move freely in the $x$ direction, while the left and right boundaries were premitted to translate freely in any direction. The upper, lower, left and right boundaries were not permitted to rotate. In the linear buckling studies, an applied load was applied to the upper and lower boundaries and the linear buckling load was found. In the post buckling studies, the upper boundaries had a displacement in the $y$ direction imposed while the lower boundaries were fixed in the $y$ direction. The reaction forces on the upper and lower boundaries were then measured. For the stationary studies, small imperfections in the shape of the first eigenmode (found through the linear buckling studies) were then added to the mesh through the use of MeshPerturb 1.0 [23].

For the measurements of the Poisson’s ratio, the displacement of the central points of the nodes within the structure were recorded. These displacements were then analysed to find the Poisson’s ratio, the displacements of a single unit cell towards the centre of the lattice were taken to be representative of the structure.

6 Funding

DRK acknowledges funding support from Academy of Finland and Aalto Science Institute. CZ acknowledges funding from Aalto Science Institute. LT acknowledges funding from Aalto Science Institute thematic program “Challenges in large geometric structures and big data”.