Exploring Adaptive Behavior of Non-linear Hexagonal Frameworks

Non-linear structural responses offer a rich design space that can be integrated into the development of novel (meta)-materials to enable transformative capabilities. By exploiting robust, repeatable, non-linear elasticity the (meta)-material's performance can be tailored to significantly exceed what has been demonstrated through traditional linear design paradigms. Embracing geometric non-linearity offers the designer the potential for truly bespoke large-range elastic behavior. Herein we explore the behavior of bistable elements that can be arranged into a space-filling triangular lattice, constructed in a hierarchical manner. We develop an analysis of the system that can be readily extended to more complex hierarchies, such as the hexagonal arrangement considered herein, whilst retaining physical insight. Specifically, we present the geometric restrictions of lattice elements that govern the existence of stable stress-free states and subsequently characterize the transitions between such states by searching for the minimum energetic transition pathway. Through our approach, we explore how hierarchical multi-stable systems can be analyzed, thus contributing to the development of truly bespoke adaptive (meta-)materials.


INTRODUCTION
Traditional structural design practice typically utilizes stiff components whose deformations are purposefully minimized in order to exploit linear elastic behavior. Where large deformations are required, such as in morphing or deployable systems, mechanical joints are typically used to achieve the desired transformative behavior. With the advent of new manufacturing capabilities and advanced analysis tools, a purely linear design philosophy is increasingly being challenged. One example of this paradigm shift is the development of compliant mechanisms that utilize reversible but large structural deformations to achieve adaptability, oftentimes facilitated by snapping across regions of static instability. Advantages of such compliant mechanisms are decreased part count, reduced design complexity and the introduction of preferable manufacturing techniques (Howell, 2001). Most importantly, in compliant systems the structural topology can be tailored to obtain desirable transformative behavior (see e.g., Lu and Kota, 2003;Jutte and Kota, 2008;Krishnan et al., 2013).
Harnessing structural instabilities is not without its own challenges, see e.g., Champneys et al. (2019) who outline recent advances in the field. One prominent example is the buckling of an axially compressed cylindrical shell, which, beyond the classical periodic buckling waveforms distributed over the entire domain of the cylinder, has a strong proclivity to localize the buckling mode over only a small portion of the domain, thereby forming localized dimples (Groh and Pirrera, 2019a).
High-speed photography experiments of axially loaded cylinders suggest that a single dimple forms to initiate buckling, and this single dimple multiplies circumferentially to complete a single row of buckles (Eßlinger, 1970). This pattern formation, known as homoclinic snaking or cellular buckling, is also observed in numerical simulations (e.g., Hunt et al., 1999Hunt et al., , 2000aGroh and Pirrera, 2019a).
Owing to the rotational invariance of the cylinder, each unique combination of dimples can occur anywhere around the circumference, and the spatial multiplicity of these localizations implies a large set of possible trajectories to instability, with each trajectory affine to a particular imperfection signature (Groh and Pirrera, 2019b). The resulting complexity of multiple and oftentimes entangled post-buckling paths of this symmetric problem creates unique challenges for numerical methods, for the interpretation of results, and for obtaining physical insight into problem. However, the sequence of instabilities can be controlled more reliably by breaking the system's inherent symmetry through engineered imperfections (Lee et al., 2016) or through tailored material properties (White et al., 2015).
Similar difficulties are indeed also faced in the (meta)material community that is developing novel material systems constructed from multi-scale arrangements of mechanically multi-stable snapping components (Mullin et al., 2007;Florijn et al., 2014;Nadkarni et al., 2014;Rafsanjani et al., 2015Rafsanjani et al., , 2019Shan et al., 2015). As the periodicity of the unit cells is scaled-up, e.g., by placing multiple unit cells in rows and columns to create a 2D array, the complexity in the number of possible equilibrium states scales geometrically . This is because under load many different deformation sequences can ensue; cells may snap individually, together as a row or column, or together as an entire array. The specific sequence of buckling is highly sensitive to the initial configuration of the system and thus presents close similarities to the cylinder buckling problem. The sequence of instabilities can thus be choreographed by introducing a bias into the system, either by varying the properties of individual unit cells or by adding initial imperfections (Che et al., 2016). These complexities often make obtaining physical insight into the behavior more difficult and the reliance on numerical methods can lead to computationally expensive design and optimization processes. Clearly, it is desirable to harness structural nonlinearities in a manner that allows physical intuition to be retained and that permits an algorithmic approach to design. One approach to meeting this challenge is to exploit constrained non-linear behavior, which can be readily understood, to develop hierarchical systems. For example, the helical lattices investigated by Pirrera et al. (2013) may be thought of as well-understood 1D non-linear springs whose behavior is tuned by varying stiffness, pre-stress and characteristic lengths (a macroscopic example is shown in Figure 1). An extension to the analytical model of Pirrera et al. (2013) for lattices with strips of nominal width, thereby including their transverse curvature and in-plane strains, demonstrates a robust consistency with finite element simulation and experimental characterizations (McHale et al., 2020b). Furthermore, the helical lattice can be tuned to exploit external fields, such as thermal effects  FIGURE 1 | Prototype of a helical lattice structure that behaves as a one-dimensional non-linear spring. McHale et al., 2020a), to modulate the effective stiffness and stability landscape. The potential of helical lattices to obtain truly bespoke elastic responses was considered by Dixon et al. (2019). By exploiting the helical lattice as a structural base-unit for nonlinear systems of coupled lattices, an algorithmic tuning process that can match, to an arbitrary precision, any energetic function up to a constant was developed (Dixon et al., 2019).
Emergent behavior is not restricted to the helical lattice alone and there exists a range of non-linear units from which emergent behavior is observed as system complexity is increased, such as cylindrical ribbons (Lachenal et al., 2012;Aza et al., 2019;Carey et al., 2019) and snapping arches . This paper seeks to extend the design concepts explored in Dixon et al. (2019) beyond their current one-dimensional nature through an exploration of hierarchical space-filling lattices. In doing so, we take one step closer toward the design of truly bespoke (meta-)materials.
At smaller length-scales, the architecture of non-linear (meta-)materials can be informed by structural behavior, thereby reducing reliance on time-consuming and complex investigations based on material science-as is the case for tuning the properties of shape memory alloys (Bhattacharya and James, 2005). Instead, using a structural approach, instabilities can be integrated into the design and exploited to obtain a unique richness in system response (Reis, 2015;Bertoldi, 2017). By employing these structural response mechanisms at sufficiently small length scales, a desirable macroscopic continuum response can be achieved (O'Donnell et al., 2016;O'Donnell et al., 2019). Herein, we propose a new approach for creating novel material behavior from a hierarchy of non-linear elements; namely, rather than utilizing complex topologies to achieve the desired non-linear behavior, we can subsume this complexity into well understood non-linear "springs" acting as base-units of a hierarchical design. By combining robust base-units in simple lattices, in this case triangular lattices, we can achieve desirable morphing characteristics.
The paper proceeds as follows. First, we identify the base bistable element from which our lattices derive. Then we explore the stable stress-free configurations of a triangular framework, the so-called first hierarchy, and describe its transitional behavior. Finally, we extend the lattice to the second hierarchy, namely a hexagonal framework, and identify the geometric restrictions that dictate the existence of stable stress-free states. We then explore the least-energy transitions between these states.

ONE-DIMENSIONAL BISTABLE ELEMENTS
We seek to investigate the behavior of frameworks composed of one-dimensional bistable elements which are non-linear springs characterized by a non-dimensional length l ∈ R + and nondimensional energy (l) 0. We assume that the two stable equilibria, which have zero energy, occur at lengths l = l 0 ∈ (0.5, 1) and l = 1. We also assume that the unstable equilibrium at l p ∈ (l 0 , 1) has energy (l p ) = 1, which is thus the energy barrier B between the stable equilibria.
The energy profile of such elements can be tuned by tailoring their characteristic properties and can even match a target energy profile with an arbitrary level of accuracy (Dixon et al., 2019). However for computational ease we choose l p = l 0 +1 2 and select to be a fourth-order polynomial with coefficients c i that satisfy The resulting energy is illustrated in Figure 2.
(2) We will identify "+" with 1 and "−" with −1 so signs can be added and multiplied, see section 5 below.

UNCOUPLED SYSTEMS. I. MINIMA, SADDLES AND MAXIMA
Consider planar hinged structures, constructed from identical bistable elements, for which the element lengths can be prescribed independently of each other, restricted only by satisfaction of the triangle inequality, see Figure 3. (The twoelement system in Figure 3A requires also specification of the fixed-support separation, L sep .) We refer to such systems as uncoupled. Their energy is simply the sum of the energy of each element,

Triangles. I. Minima
In particular, the triangular structure in Figure 3B has 8 stable stress-free configurations, namely when its sides are of length either l 0 or 1. We label the (interior) angles of these configurations as follows. Consider a triangle with sides s 1 , s 2 , s 3 ∈ {l 0 , 1}. We denote the angle subtended by the sides s 1 , s 2 by a symbol chosen from {φ 0− , φ 00 , φ 0+ , φ ±− , φ ±+ } as follows: The first subscript is 0 if s 1 = s 2 and is ± otherwise.
From the cosine rule, Moreover, by considering the non-equilateral configurations we obtain, φ 0+ + 2φ ±− = π; (5b) and a simple calculation shows that the angles are ordered as.

Triangles. II. Saddles and Maximum
Restricting the lengths of the elements to the interval [l 0 , 1], we observe that local minima of the system energy correspond to configurations where all elements are at a local minimum of , i.e., of length either l 0 or 1. The system's only maximum occurs when all elements are at a local maximum of , i.e., of length l p . Saddle points, i.e., extrema that are neither local maxima or minima (equivalently, where the Hessian is indefinite), correspond to configurations where some elements are extended to the local maximum l p with the remaining elements at either minima.
Thus, in Figure 4, the energy minima are located at the vertices, the maximum at the center and saddles at the centers of the edges/faces. In particular, in Figure 4B, saddles with an energy of 1 occur along the midpoint of each edge and saddles with an energy of 2 at the center of each face.

Generalization to Larger Systems
Considering Figure 5, it is clear that the triangular framework can be extended to larger frameworks through addition of pairs of elements to any free edge (Graver, 2001). Provided the lengths of the new elements satisfy the triangle inequality as part of the newly formed triangle, their size can be altered independently. Thus we may describe the energy of such an independent system of N elements through Equation (3). In this case, analogous to the three element system, the energy minima are located precisely at vertices of the hypercube [l 0 , 1] N . Saddles are located along all 1D edges (these have energy 1) and, more generally, on each face (a hyperplane), with energies determined by the number of unstable elements activated. A single interior maximum is located at the center of the hypercube.

UNCOUPLED SYSTEMS. II. TRANSITIONS BETWEEN MINIMA
When transitioning between stable states it is often desirable to select the path that minimizes the actuation energy. The barrier to actuation can be considered analogous to the Maxwell energy (see e.g., Hunt et al., 2000b;Wadee and Edmunds, 2005;Champneys et al., 2019). For a system to transition between two minima, it must overcome an energetic barrier, the smallest possible barrier is characterized by one or a series of adjacent saddle points. As individual elements have an activation energy of (l p ) = 1, any transition path must encounter an energetic barrier that is at least this magnitude. Thus, for systems of independent elements, the energetically optimal paths between minima are paths traversing the edges of the hypercube [l 0 , 1] N . Intuitively, such a path makes physical sense, by actuating each element in isolation the individual element's energy barrier is never exceeded.

Two-Element System
Consider the simple system described in Figure 3A. To ensure the existence of four stable equilibra, we set L sep = 1. Figure 6A illustrates the energy landscape and shows the expected four stable equilibria and central maximum. The normalized negative gradient, − ∇ N ||∇ N || , in Figure 6B, highlights the stability of the system and matches the expected behavior, Figure 4A.
A systematic exploration of the domain can be conducted by investigating the transition from a fully contracted to a fully extended configuration, i.e., moving from a configuration where all element lengths are l 0 to one where all lengths are 1. This is always possible for an uncoupled system subject to suitable boundary conditions. This corresponds to two opposite corners of the square domain in Figure 4A. Consider the parametrized path between the two stable configurations given by, where n ∈ {1, 2, . . . , N}, α ∈ [1, ∞) and t ∈ [0, 1]. As α → ∞ this curve approaches an edge traversal path. The vertices to be visited along the path can be selected by permuting the order of the powers of t. Owing to the system's symmetry, this is sufficient to explore the entire domain, even though each path sweeps only a portion of the domain. We can define the energy barrier of a path for a given α as, However for visualization it is more convenient to use the variablet K = (1 − t) 1 K , where K ∈ N is a scaling parameter that can be increased as required to account for the number of elements in the system. Figure 7A illustrates the paths of Equation (7) for a range of α. Figure 7B shows the system's energy during traversal along the path, notably we observe the emergence of two distinct peaks and the reduction in peak magnitude as α increases. Finally, Figure 7C, demonstrates that as the path approaches an edge traversal the energy barrier tends toward 1, i.e., the energetic barrier of actuating a single element.

Three-Element System
Similarly, for the three-element system in Figure 3A, paths for a range of α are illustrated in Figure 8A. Figure 8B shows the system's energy during the traversal, notably we observe the emergence of three distinct peaks and the reduction in peak magnitude as α is increased. Finally, Figure 8C, demonstrates that as the path approaches an edge traversal the energy barrier tends toward 1.

HEXAGONAL FRAMEWORK. I. STABLE STRESS-FREE STATES
We now proceed to analyse a hexagonal framework to illustrate the emergent behavior observed for non-independent systems. Consider six triangles from section 3.1 that form a planar hexagon as illustrated in Figure 9A. In this section we characterize the stable stress-free states of such hexagons.  We track the positions of the nodes of a hexagon through complex numbers: The central node is at 0 and the other six nodes, numbered anticlockwise, are at z i ∈ C, i = 0, . . . , 5, subject to the constraints: where addition of indices is mod 6. Set so |w i | is the proportion between the lengths of the adjacent radial edges and arg w i is the interior angle between adjacent radial edges. It follows that a hexagon can be identified with w i ∈ C, i = 0, . . . , 5, which satisfy, j i=0 w i ∈ {l 0 , 1, 1 l 0 }, j = 1, . . . , 5, With every hexagon we associate a sign diagram, which is a (finite) sequence of signs, s 0 s 1 s 2 s 3 s 4 s 5 , (11a) modulo 6-cycle permutations, where as defined in Equation (2). Now consider the sign diagram Equation (11a), modulo 6cycle permutations. Note that 1. From Equation (10b), the signs "+" and "−" should occur in equal number. 2. Considered cyclically "+" and "−" should alternate with "0" allowed to intervene, since from Equation (10a) neither two "+"s nor two "−"s can occur sequentially.
From Item 1, there are four classes of sign configurations. These, labeled by the numbers of each sign, are tabulated in Table 1. We now consider each of these configurations.
(for C) Thus, this configuration can be realized only in a regular hexagon. There are two possibilities: Either all edges of the hexagon are of side l 0 or all edges of the hexagon are of side 1.

Configuration 3-0-3
From Item 2 above, the only possible sequence of signs is considered cyclically. From section 3.1, each interior angle could be either φ ±− or φ ±+ . We can choose p ∈ {0, . . . , 6} of these to be φ ±− provided Since, from Equation (6), φ ±− < π 3 < φ ±+ , p / ∈ {0, 6}. A numerical investigation shows that the only solution is p = 2 and l 0 = −1+ √ 5 2 ≈ 0.6180, which is the inverse golden ratio. A number of hexagons are possible, one of which is illustrated in Figure 10A. However, since we are interested in generic configurations, i.e., those that are feasible for any l 0 ∈ (0.5, 1), henceforth we shall omit details of hexagons that are exist only for some values of l 0 .

Configuration 1-4-1
Modulo multiplication by −1 there are three sub-configurations considered cyclically. To see this note that we can choose to traverse the circuit such that (i) "−" precedes "+, " and (ii) the number of "0"s between "−" and "+" is no more than the number of "0"s between "+" and "−." The last sub-configuration is its own negative; the negatives of the first two are considered cyclically, respectively. As for interior angles, corresponding to each "−" or "+" we may choose either φ ±− or φ ±+ . Assume, with no loss of generality, that we traverse the interior angles so that the "−" precedes the "+." Then each "0" between "−" and "+" can be either φ 00 or φ 0+ ; and each "0" between "+" and "−" can be either φ 00 or φ 0− .
The generic configurations are listed in Table 2.
This allows the total energy to be written as 11 C (l) = 11 (l) + G (l).
The analysis in section 5 identified the global minima with zero energy. We note that the system may, in general, possess additional local minima, but these must have an increased energetic state.
Using the approach identified in section 4, a parametric sweep of the 11-dimensional hypercube [l 0 , 1] 11 can be conducted. Figure 11 shows the system's energy along the path and the emergence of 11 distinct peaks and convergence of the barrier energy to 1 as the path approaches an edge traversal.
We can explore the energy function G in a similar manner, see Figure 12. Unlike the independent case, the system's response is dependent on the relative position of the secondary minimum l 0 . There is a significant region of the domain where the energetic barrier does not exceed 1. Geometrically, this indicates that in order for the other independent elements to transition, the dependent length must be extended and contracted in tandem, i.e., there is no path that allows the length of l 12 to remain constant. Additionally, this energetic barrier indicates that the actuation path does not require the dependent length l 12 to be excessively contracted or extended, but that it may undergo multiple transitions across its own local maximum. Additionally, we observe that as l 0 approaches 0.5 the energy barrier associated with the geometric constraint significantly begins to increase for edge traversal paths. This indicates that in order to accommodate such large geometric reconfiguration the dependent element must be compressed, or extended, significantly. The behavior of the combined system follows from these observations indicating that a path with an energy barrier of 1 cannot be achieved, but that for many geometries an edge traversal remains an effective approach to minimize the energetic cost.

CONCLUSIONS
Structural non-linearity offers an expanded design space that can be exploited for the development of adaptive (meta-)materials. The traditional approach to developing transformative (meta-)materials has often relied upon optimization through empirical material science investigations, or the introduction of complex topological structures. Here, instead, we embed non-linearity, in the form of bistability, at the elemental level, recasting the problem as the analysis as simple frameworks of nonlinear springs. In our analysis we investigate the energetic description of these frameworks. In doing so we identify geometric conditions governing the existence of global energetic minima of the system. Many frameworks can be constructed such that independent variation of the elemental lengths is possible, subject only to satisfaction of the triangle inequality. For these systems, all minima, saddles, and the (unique) maximum, can be obtained directly by considering a hypercube defined by the relative position of the elemental local minima. We use a systematic path search to determine the minimum energetic barrier to transition between self-equilibrated states, analogous to the Maxwell energy. These results support the expected conclusion that actuation of the elements in a sequential manner minimizes the energetic barrier, corresponding to an edge traversal of the hypercube.
For more constrained systems, such as the hexagonal framework we present herein, not all elemental lengths can be selected independently. In these systems, the ratio of the relative positions of elemental local minima determine the geometric feasibility of system-wide global minima. We present the complete set of constraints governing geometric feasibility and identify special elemental ratios that permits adaptive states that are not generically permissible.
In order to explore the transitional behavior of the hexagonal system we consider its energetic description as the combination of independent energy, and an additional "geometric" perturbation. The geometric perturbation arises from the dependency of one elemental length on the remaining 11 elements of the system. Our transitional analysis demonstrates that the dependent element introduces an energetic penalty that cannot be eliminated. That is, we can no longer devise a path that actuates each element in turn. The extent of the energetic penalty is dependent on the elemental ratios, however there are significant regions of the design space that possess nominally identical penalties.
Planar frameworks are of interest even in three-dimensional contexts, e.g., the cross-sectional form of a morphing aerofoil. Moreover, the energetic behavior of uncoupled frameworks is identical in two and three dimensions. For coupled hexagonal systems, as explored in sections 5 and 6, out-ofplane deformations of the nodes may be energetically preferable but, as with the present planar analysis, the system will seek geometrically-preferable configurations that can be identified through global energy minima. The benefits of exploiting non-planar responses of rectangular lattices to develop smart surfaces was explored by Zhang et al. (2017). Three-dimensional frameworks composed of non-linear elements will be considered in future work.
The approach we have presented provides an initial foray into the exploration of the behavior of hierarchical multistable frameworks using a systematic analytical approach. Our approach offers physical intuition into system behavior via two mechanisms. Firstly, the behavior of an independent system can be readily understood owing to its characterization via the hypercube that captures all stable stress-free states. For constrained systems, we can isolate the effects of the constraints, via a perturbation function. By replacing complex topological forms with elements of non-linear stiffness we have provided an alternative approach to tuning global response mechanisms of adaptive (meta-)materials. There is a balance to be struck between elemental non-linearity and topology complexity to yield the most efficient design. This paper begins an exploration of this concept.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
MO'D and IC formulated the initial research problem. MO'D, IC, and MT developed the analytical framework. MO'D, IC, and RG conducted the numerical analysis. All authors contributed to the analysis, discussion, and preparation of the manuscript.