Equally probable positive and negative Poisson's ratios in disordered planar systems

Abstract

Auxetic materials, characterised by a negative Poisson's ratio, have properties that are different from most conventional materials. These are a result of the constraints on the kinematics of the material's basic structural components, and have important technological implications. Models of these materials have been studied extensively, but theoretical descriptions have remained largely limited to materials with an ordered microstructure. Here we investigate whether negative Poisson's ratios can arise spontaneously in disordered systems. To this end, we develop a quantitative description of the structure in systems of connected basic elements, which enables us to analyse the local and global responses to small external tensile forces. We find that the Poisson's ratios in these disordered systems are equally likely to be positive or negative on both the element and system scales. Separating the strain into translational, rotational and expansive components, we find that the translational strains of neighbouring basic structural elements are positively correlated, while their rotations are negatively correlated. There is no correlation in this type of system between the local auxeticity and local structural characteristics. Our results suggest that auxeticity is more common in disordered structures than the ubiquity of positive Poisson's ratios in macroscopic materials would suggest.

---

1 Introduction

Materials with a negative Poisson's ratio (PR), called auxetics, have positively correlated horizontal and vertical strains, and thus expand (contract) in both directions when stretched (compressed).¹ Various such materials can be found in nature, including polymers, foams, minerals, and even nuclei of stem cells.^2,3 A range of auxetic structures has also been made artificially.^4,5 Auxetic behaviour in disordered materials is not uncommon; for example, a crumpled ball of paper expands in all directions when stretched between two fingers. Yet, most investigations of the relation between internal structural characteristics and large scale auxeticity have focused on ordered systems.^6–13 As such, few studies exist of auxeticity in disordered systems, including behaviour of three-dimensional folded sheets^14,15 and perturbing slightly ordered structures by defects.^16,17 Thus, auxeticity of disordered structures is still far from fully understood.

We address this problem by modelling auxeticity in isostatic structures that consist of minimally connected constituent units that can freely fold, expand and contract.^18,19 Isostatic systems are statistically determinate—a particularly convenient property for our purposes, as we can relate the stress to the local microstructure by calculating the inter-element forces directly, without requiring elasticity theory. This feature also makes these structures ideal for modelling the jamming transition in systems in which the elements are macroscopic grains or colloids. We use arrangements of isostatic systems to analyse random disordered isostatic systems, and investigate the microscopic drivers of auxeticity on the local and global scales.

To be isostatic, a mechanically stable structure has only to satisfy a minimal connectivity criterion: that the mean number of force-carrying inter-element contacts per element is equal to a specific value z_c.²⁰ Generically, z_c = d(d + 1) and d + 1 for d-dimensional systems of frictionless and frictional non-spherical elements, respectively.²¹ Analysis of stress transmission in these materials explains the ubiquity of non-uniform stress states exhibited by particulate media.^22–25 A first-principles continuum stress theory for two-dimensional isostatic granular media has been developed, based on a parameterisation of the inter-element forces into ‘loop forces’, where loops are the elementary voids, or cells, enclosed by individual elements.^26,27

Since z_c = 3 in two dimensions, we choose to focus on planar systems of triangles, connected to nearest neighbours at the vertices by conceptual frictionless hinges, as illustrated in Fig. 1. For later use, we assign a direction to the edges of each triangle, t, making it a vector, r⃑_ct, with c being the cell that the edge borders. Thus defined, these vectors circulate clockwise around the triangles.²⁶ For each r⃑_ct, we define a dual vector R⃑_ct, extending from the centroid of triangle t to the centroid of cell c. These two vectors form the diagonals of a quadrilateral, known as a quadron.²⁸ The quadrons are fundamental volume elements that tessellate the system space perfectly. Moreover, the structure of each quadron can be quantified unambiguously by a local structural tensor, r⃑_ct⊗R⃑_ct, allowing us to quantify the local disordered structure anywhere in the system.^26,28,29 This is convenient for relating local structural characteristics to local PR.

Fig. 1 Close-up view of individual triangles in our system. The edges of triangle t are made into vectors, r⃑_ct, that circulate the triangle clockwise. These vectors make a disordered graph of randomly shaped triangles, connected to their nearest neighbours at the vertices. The vector R⃑_ct extends between the centroid of triangle t and that of its neighbouring cell, c. The vectors r⃑_ct and R⃑_ct are the diagonals of a quadrilateral volume element, which is called quadron.²⁹ Inter-triangle forces (red) are calculated as the difference between the loop forces l⃑_c (blue) of neighbouring cells.

In mechanical equilibrium, the forces of triangles on their neighbours (inter-triangle forces or ITFs), can be parametrised using loop forces, with one loop force l⃑ per cell enclosed by triangles. As such, the force that triangle t exerts on its neighbour t′ is the difference between the loop forces of the two cells, c and c′, straddling the common joint of these triangles²⁶ (see Fig. 1):

f⃑_tt′ = l⃑_c′ − l⃑_c (1)

By construction the loop forces automatically satisfy force balance on each triangle, leaving only the torque balance conditions. This cuts by two thirds the number of equations to solve, considerably reducing the computational effort.

By calculating the fabric tensor, the quantitative description of the local structure of our isostatic systems can be coarse-grained to the overall system.^26,27,30 Furthermore, the fabric tensor is related to local rotational strains,³¹ and can be used to identify and isolate rotational strains in auxetic materials.^18,19 In this paper, we extend this latter development and apply it to random systems.

2 Theory

We first determine the number of equations required to calculate the inter-triangular forces in a disordered planar system consisting of N (≫1) triangles. The contours of the triangles form a graph of N_v connected vertices, N_e edges, and N_f faces. Identifying the N_b boundary triangles in this structure, it is convenient for the purpose of our analysis to enclose this graph within a frame to which the boundary triangles are connected by one of their vertices. For such a graph, Euler's topological relation for two-dimensional graphs in the plane is,³²

N_v − N_e + N_f = 1, (2)

where N_v, N_e and N_f are the numbers of the graph's vertices, edges and faces, respectively.

In our graph, the vertices are the analogues of the contact points in granular assemblies. The edges are the triangle edges and the additional N_b edges formed by the frame. The faces consist of the triangle faces and the internal elementary voids (called cells or loops here), which include the additional N_b voids formed by the frame. Putting these numbers together, we have: N_f = N_c + N + N_b faces, (the N_c internal cells, the N loops around each triangle and the N_b boundary cells), N_e = 3N + N_b, and N_v = (3N + N_b)/2. Substituting these into Euler's relation (2), we have

Since one loop force is associated with each cell, there are N_c + N_b loop forces altogether, which yield, using (3), 2 + N + N_b unknowns. Furthermore, each ITF is defined as the difference between two loop forces, meaning we can fix one of the loop forces at will.²⁶ This reflects the nature of the loop forces field as a potential field, of which the local ITF is the gradient.²⁶ Having exactly N torque balance conditions to determine the loop forces, we therefore need to fix N_b/2 boundary forces on the vertices connected to the bounding frame. This leaves an exact set of equations to solve for the loop forces. We can then solve for the internal N/2 loop forces, from which we can determine the ITFs.

Specifically, we can define the following N × N matrix A, loop forces vector X, and boundary forces vector B.³³ As before, the first N rows correspond to the torque balance on all N triangles, the next N_b rows correspond to the boundary forces on N_b/2 vertices, and the final 2 rows normalise the 2 force components of the first loop force:

andwe have

Ax = B (7)

We can solve for the loop forces by Gaussian elimination.

3 Model and methods

We generated numerically 15 000 disordered planar assemblies of triangles as model systems. First, we distributed 200 points randomly within a rectangle of 2000 × 2000 in arbitrary units. We then generated a Voronoi tessellation around the points of the system, obtaining a planar cellular structure, with exactly three edges emanating from each node. Connecting the midpoints of these edges around each node produces a disordered assembly containing about 300 triangles in each system. With exactly 3 = d + 1 contacts per triangle, we made this system isostatic by imposing one force on every other boundary triangle vertex. We then used (7) to calculate the loop forces (see, e.g.Fig. 2a) and, from those, the ITFs by (1). Typically, the solutions contain regions in which the loop forces have smoothly-changing directional orientation, separated by sharp gradients, as illustrated on the right side of Fig. 2a. The sharp gradients correspond to localised force chains.^26,27,34,35 These force chains can be observed experimentally in granular systems.^36–38

Fig. 2 The force field response to two antipodal tensile forces, applied at the centre of opposite boundaries (in blue, not to scale). Shown in red are: (a) the loop forces, l⃑_c; (b) the corresponding ITFs.

Starting from the initial stress state, we calculated the mean force amplitude at a vertex for each system, f̄, and applied two equal and opposite small tensile forces, of amplitude about δf = 1 × 10⁻³f̄, antipodally on the left and right boundaries, along the x-axis near the middle of the system (shown as blue forces in Fig. 2). We then calculated the changes in the loop forces and the ITFs.

The changes in the ITFs give rise to vertex displacements, resulting in triangle rotations, translations and expansions (or contractions), such that the triangle edges remain straight. It should be noted that the static determinacy constrains all the degrees of freedom and, therefore, the vertex displacements are only the result of the triangles’ elastic compliance. In the absence of inertial effects, each vertex's displacement is proportional to the force on it, and we expect vertices to displace towards the external tensile force it is closest to. This divides the system naturally into two halves with an imaginary half-line at x = 0. The vertex's displacement is then in the direction of the force pulling toward the closest boundary, scaled by a proportionality factor ζ that is 0 on the half-line and 1 on the corresponding boundary. A typical example of a displacement field is shown in Fig. 3a.

Fig. 3 (a) Vector [small rho, Greek, vector]_ti points from the centroid of triangle t to its vertex i, and this vertex is displaced by δ[small rho, Greek, vector]_ti. (b) Displacement of vertices in response to 100 small increments of the boundary forces.

The displacements data were then recorded for later analysis, increased the magnitudes of the boundary tensile forces, calculate the new ITFs, and repeated the process. The accumulated displacements resulting from 100 successive repetitions of the process are shown in Fig. 3b.

In the following, we refer to the gradient of the displacement as strain, for short, although this term is normally reserved for its symmetric part. The strain can be defined for individual triangles and can be separated into a translational, an expansive and a rotational component:

ε_t = ε^trans_t + ε^exp_t + ε^rot_t. (8)

In our system of non-rigid triangles, all these components are relevant,¹⁹ while in systems consisting of rigid triangles, only the translational and rotational components contribute to the total strain.²⁶

The translational strain captures the displacements of the triangles’ centres of mass and it is defined as

where R⃑_t is the vector from the origin to the centroid of triangle t, R⃑₀ is the position of the entire system's centre of mass, and δR⃑ is the mean displacement of the vertices of t. The expansive component of the strain depends on the outward or inward displacements of the vertices of the triangle, projected on the vector ρ_tiwhere i = 1, 2, 3 are the three vertices of the triangle. The expansive strain is then

Each triangle hosts exactly three quadrons, as shown in Fig. 1, each described by a structure tensor C_ct = r⃑_ct⊗R⃑_ct.²⁸ The triangle's structure tensor is the sum of the structure tensors of its three quadrons:

The antisymmetric part of this tensor gives the (generically, non-convex) area associated with the triangle's three quadrons, A_t,where C^T_t is the transpose of C_t and Θ is the two-dimensional π/2 rotation matrix (the Levi-Civita operator). The symmetric part of C_t,describes the local rotation of the triangle relative to the global mean.^26,27 It has also been shown³¹ that the rotational strain satisfiesIn terms of the triangle variables, the rotational component of the displacement is the displacement perpendicular to the direction of the expansionTherefore, the rotation angle of vertex i isand the rotational strain of the entire triangle is

4 Results

Defining the PR of triangle t as

ν_t = ε_t,yy/ε_t,xx, (19)

provides a local definition of the PR and makes it possible to study the spatial distribution of this property across the system.

In Fig. 4a, we plot a histogram of the local PRs for data collected from over 40 000 triangles from 100 different system realisations. It is broadly symmetric, with the mean and standard deviation being −0.2 and 1.8, respectively, demonstrating that there is hardly any sign preference for the local PRs.

Fig. 4 (a) Histogram of the PRs of individual triangles, ν_t, calculated from over 40 000 triangles from 100 simulated systems. (b) Histogram of the PR of global systems, calculated from over 15 000 systems of 300 triangles each. (c) Mean (blue; left axis) and standard deviation (red; right axis) of system of 150 to 1050 triangles each. The mean is fitted well by PR = −1.8 × 10⁻⁴N + 0.8 and the standard deviation by [small sigma, Greek, macron]_PR = 3.8 × 10⁻³N + 5.7, with N the number of triangles.

We also computed the global PR of each system, defined as

Here, L_α is the system extension in the α direction, defined as the difference between the mean positions of the vertices at opposite boundaries, and ΔL_α is the mean change in that distance. In Fig. 4b we show the histogram of ν_sys across 15 000 systems of 300 triangles each. This histogram is also symmetric, with its mean and standard deviation 0.9 and 7.0, respectively. Interestingly, this histogram's standard deviation is much larger than for the individual triangles, which is a result of its tails falling slower than exponential, although faster than algebraically.

To examine the dependence of these results on system size, we repeated the simulations for systems of sizes ranging from N = 150 to N = 1050 triangles each. We find that the distributions’ standard deviations are an order of magnitude larger than their means and are broadly symmetric around 0, as can be observed in Fig. 4. Zooming in on the means and standard deviations of the global PR distributions over the 7-fold size variation, we observe that the mean decreases slightly with N, from 0.82 to 0.66. This decrease is not statistically significant and is also consistent with no change within 95% confidence. In contrast, the standard deviation increases significantly with size by about 53%, from 6.31 to 9.70, indicating, with 95% confidence, that the width of the PR value distribution increases proportionally to N.

To gain insight into the relationship between the Poisson's ratio and the local structure, we investigated the local correlations between structural, elastic and strain properties, with the structural characteristics quantified using the above quadron description. To this end, we first express the stress on triangle t in terms of the loop forces of the cells surrounding it.²⁶ The increase in the boundary forces gives rise to the following change in the stress

Given a triangle's stress and strain, we also calculated its compliance matrix, S_t, using σ_t = S_tε_t.

In particular, we correlated the traces of the different strain components of individual triangles with the trace of that triangle's total strain, for 100 systems and 3 × 10⁴ triangles. Irrelevant outliers, within the 2% highest and lowest strains, were disregarded. We find that the expansive component is the dominant contributor to the symmetric part of the strain, with a correlation of γ = 0.93 between Tr{ε^exp_t} and Tr{ε^tot_t} (Fig. 5a). The second largest contributor is the translational strain, with a correlation of γ = 0.24 between Tr{ε^trans_t} and Tr{ε^tot_t} (Fig. 5b). In contrast, there is hardly any correlation between the rotational component and the total strain, γ = 0.03 (Fig. 5c). As expected, the importance of expansion in the overall strain of individual triangles causes a negative correlation of γ = −0.10 between the expansive strains of neighbouring triangles (Fig. 5d). The translational strains of neighbours are positively correlated, with γ = 0.67, which is induced by our method of determining vortex displacements. Interestingly, while the rotational strain does not correlate well between nearest neighbours, with γ = −0.06, the rotation angle of neighbours is anti-correlated, with γ = −0.35, which is an aspect of the ‘anti-ferromagnetic’-like rotational nature of such systems.^{11,31,39–41}

Fig. 5 The relative contribution of strain components to the total strain of individual triangles. We show the correlation between the traces of a triangle's: (a) expansive and total strain; (b) translational and total strain; (c) rotational and total strain. In (d) we show correlation between the trace of a triangle's expansive strain and the average of the trace of the expansive strain components of the triangle's neighbours.

We also investigated the correlations between the one-triangle PR and local structural characteristics, including the triangle: size, shape, quadron sizes, displacement, and any of the strain contributors. We found no significant correlation between any of these quantities.

5 Conclusions and discussion

We studied the poorly understood issue of auxetic behaviour in disordered structures, modelled here by planar statically determinate, or isostatic, systems of triangular elements. Specifically, we addressed the question of whether a negative Poisson ratio (PR) can emerge spontaneously without design or whether a positive PR is more common, as its ubiquity in macroscopic materials might suggest.

We simulated 15 000 disordered planar systems, each of about 300 randomly-shaped connected compliant triangles, stretched them incrementally, and computed the PR for each triangle and for the entire assembly. Our main finding is that positive and negative PRs are almost equally likely to occur in such disordered systems. This conclusion holds on both the triangle scale (Fig. 4a) and the entire system (Fig. 4b). This suggests that large-scale systems with this layout would be almost PR-neutral, displaying a behaviour that is neither conventionally elastic nor fully auxetic. This result may seem counter-intuitive in view of the prevalence of positive PRs in macroscopic structures. However, it reinforces that these are not conventional elastic systems and that PR of such structures should be carefully analysed.

The mean contributions of the translational, expansive and rotational strains to the total strain of a triangle were found to be 23.2 ± 0.05%, 77.3 ± 0.05% and 0.0 ± 0.05%, respectively. The corresponding correlations between those and a triangle's total strain were 0.24 ± 0.005, 0.93 ± 0.005, and 0.03 ± 0.005, respectively. We also found a significant anti-correlation of −0.10 between the expansive strains of neighbouring triangles, as well as a considerable anti-correlation of −0.35 between their rotations. The latter is directly related to antiferromagnetic-like rotational dynamics in such systems, discussed in granular systems.^31,39–41 We emphasise that, by imposing static determinacy, we constrained movement in our systems and, had the triangles been rigid, they would have no degrees of freedom to deploy and not displace. Thus, the displacements we observe arise from local expansions and contractions, which fits well with the correlations we observe.

We also searched for correlations between the sign of the PR and local microstructural descriptors. Such correlations, if found, could lead to informed design of auxetic disordered structures. However, we found none in this study. It is possible that such correlations have been washed out by the strong anti-correlation between nearest-neighbour expansive strains. To study this issue in more depth, a similar study should be carried out on systems of fully rigid triangles⁹ with fewer boundary constraints. We conjecture that such systems would show stronger correlations between the local PR and the translational and rotational strains. The latter is directly related to the local structural rotational disorder, defined in ref. 26 and 27 as the symmetric part of the local tensor Q_t(r⃑) and we therefore expect that such systems would show local structure–PR correlations.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the Engineering and Physical Sciences Research Council (1506089) (C. M. V.), and the UK Medical Research Council and Wellcome Trust (K. J. C.).

References

1. R. S. Lakes, Science, 1987, 235, 1038–1040.
2. S. Pagliara, K. Franze, C. R. McClain, G. W. Wylde, C. L. Fisher, R. J. M. Franklin, A. J. Kabla, U. F. Keyser and K. J. Chalut, Nat. Mater., 2014, 13, 638–644.
3. A. C. Hodgson, C. M. Verstreken, C. L. Fisher, U. F. Keyser, S. Pagliara and K. J. Chalut, Lab Chip, 2017, 17, 805–813.
4. A. Alderson and K. L. Alderson, Proc. Inst. Mech. Eng., Part G, 2007, 221, 565–575.
5. J. N. Grima, D. Attard, R. Gatt and R. N. Cassar, Adv. Eng. Mater., 2009, 11, 533–535.
6. K. E. Evans, Endeavour, 1991, 15, 170–174.
7. D. R. Reid, N. Pashine, J. M. Wozniak, H. M. Jaeger, A. J. Liu, S. R. Nagel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E1384–E1390.
8. J. N. Grima and K. E. Evans, J. Mater. Sci. Lett., 2000, 19, 1563–1565.
9. J. N. Grima, A. Alderson and K. E. Evans, Phys. Status Solidi B, 2005, 242, 561–575.
10. J. N. Grima, P. S. Farrugia, C. Caruana, R. Gatt and D. Attard, J. Mater. Sci., 2008, 43, 5962.
11. J. N. Grima and K. E. Evans, J. Mater. Sci., 2006, 41, 3193–3196.
12. D. Attard, E. Manicaro and J. N. Grima, Phys. Status Solidi B, 2009, 246, 2033–2044.
13. R. Gatt, L. Mizzi, J. I. Azzopardi, K. M. Azzopardi, D. Attard, A. Casha, J. Briffa and J. N. Grima, Sci. Rep., 2015, 5, 8395.
14. A. S. Balankin, D. S. Ochoa, E. P. León, R. C. M. de Oca, A. H. Rangel and M. A. M. Cruz, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 125421.
15. O. Bouaziz, J. Masse, S. Allain, L. Orgéas and P. Latil, Mater. Sci. Eng., A, 2013, 570, 1–7.
16. J. N. Grima, S. Winczewski, L. Mizzi, M. C. Grech, R. Cauchi, R. Gatt, D. Attard, K. W. Wojciechowski and J. Rybicki, Adv. Mater., 2015, 27, 1455–1459.
17. T. Mukhopadhyay and S. Adhikari, Mech. Mater., 2016, 95, 204–222.
18. R. Blumenfeld, Mol. Simul., 2005, 31, 867–871.
19. R. Blumenfeld and S. F. Edwards, J. Supercond. Novel Magn., 2012, 25, 565–571.
20. J. C. Maxwell, Philos. Mag., 1864, 27, 294–299.
21. R. C. Ball and D. V. Grinev, Phys. A, 2001, 292, 167–174.
22. J. P. Wittmer, P. Claudin, M. E. Cates and J.-P. Bouchaud, Nature, 1996, 382, 336–338.
23. J. P. Wittmer, M. E. Cates and P. Claudin, J. Phys. I, 1997, 7, 39–80.
24. M. E. Cates, J. P. Wittmer, J.-P. Bouchaud and P. Claudin, Phys. Rev. Lett., 1998, 81, 1841–1844.
25. R. Blumenfeld, S. F. Edwards and S. M. Walley, Granular Systems, The Oxford Handbook of Soft Condensed Matter, 2016, ch. 5, pp. 167–222.
26. R. C. Ball and R. Blumenfeld, Phys. Rev. Lett., 2002, 88, 115505.
27. R. Blumenfeld, Phys. Rev. Lett., 2004, 93, 108301.
28. R. Blumenfeld and S. F. Edwards, Phys. Rev. Lett., 2003, 90, 114303.
29. R. Blumenfeld and S. F. Edwards, Eur. Phys. J. E: Soft Matter Biol. Phys., 2006, 19, 23–30.
30. M. Gerritsen, G. Kreiss and R. Blumenfeld, Phys. Rev. Lett., 2008, 101, 098001.
31. R. C. Ball and R. Blumenfeld, Philos. Trans. R. Soc. London, Ser. A, 2003, 361, 731.
32. H. S. M. Coxeter, Regular Polytopes, Dover Publications, New York, 3rd edn, 1973, p. 368.
33. G. Camisasca, MSc thesis, University of Cambridge, 2013.
34. J.-P. Bouchaud, P. Claudin, D. Levine and M. Otto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 4, 451–457.
35. R. Blumenfeld, New J. Phys., 2007, 9, 160.
36. R. P. Seelig and J. Wulff, Trans. Am. Inst. Min., Metall. Pet. Eng., 1946, 166, 492–504.
37. T. Wakabayashi, Proc. 7th Jpn. Nat. Cong. Appl. Mech., 1957, 153.
38. C. h. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, S. Majumdar, O. Narayan and T. A. Witten, Science, 1995, 269, 513–515.
39. R. Blumenfeld, S. F. Edwards and M. Schwartz, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 32, 333–338.
40. R. M. Baram, H. J. Herrmann and N. Rivier, Phys. Rev. Lett., 2004, 92, 044301.
41. R. Blumenfeld and S. F. Edwards, J. Phys. Chem. B, 2009, 113, 3981–3987.

---