A study of deformation localization in nonlinear elastic lattices
Raj Kumar Pal, Federico Bonetto, Luca Dieci, Massimo Ruzzene
aa r X i v : . [ c ond - m a t . s o f t ] J a n A study of deformation localization in nonlinear elasticlattices
Raj Kumar Pal a , Federico Bonetto b , Luca Dieci b and Massimo Ruzzene a,c, ∗ a School of Aerospace Engineering, Georgia Institute of Technology, Atlanta GA 30332, USA b School of Mathematics, Georgia Institute of Technology, Atlanta GA 30332, USA c School of Mechanical Engineering, Georgia Institute of Technology, Atlanta GA 30332, USA ∗ Corresponding author. E-mail: [email protected]
January 31, 2018
Abstract
The paper investigates localized deformation patterns resulting from the onset ofinstabilities in lattice structures. The study is motivated by previous observationson discrete hexagonal lattices, where the onset of non-uniform, quasi-static deforma-tion patterns was associated with the loss of convexity of the interaction potential,and where a variety of localized deformations were found depending on loading con-figuration, lattice parameters and boundary conditions. These observations are hereconducted on other lattice structures, with the goal of identifying models of reducedcomplexity that are able to provide insight into the key parameters that govern the on-set of instability-induced localization. To this end, we first consider a two-dimensionalsquare lattice consisting of point masses connected by in-plane axial springs and verticalground springs. Results illustrate that depending on the choice of spring constants andtheir relative values, the lattice exhibits in-plane or out-of plane instabilities leadingto folding and unfolding. This model is further simplified by considering the one-dimensional case of a spring-mass chain sitting on an elastic foundation, which may beconsidered as a discretized description of an elastic beam supported by an elastic sub-strate. A bifurcation analysis of this lattice identifies the stable and unstable branchesand illustrates its hysteretic and loading path-dependent behaviors. Finally, the latticeis further reduced to a minimal four mass model which undergoes a folding/unfoldingprocess qualitatively similar to the same process in the central part of a longer chain,helping our understanding of localization in more complex systems. In contrast to thewidespread assumption that localization is induced by defects or imperfections in astructure, this work illustrates that such phenomena can arise in perfect lattices as aconsequence of the mode-shapes at the bifurcation points. Introduction
Localized deformations resulting form instabilities arise naturally in a wide range of phys-ical systems, and across multiple length scales. Manifestations include plastic twinning inmetals [1], localized buckling of epithelial cells in biological media [2], and Chevron foldsin rocks [3], among others. The study of instabilities, both at the microstructural scalein materials and at the macroscopic structural level, is an area of renewed interest, anda timely research topic. Several studies have focused on exploiting the formation of pat-terns resulting from instabilities in engineered materials and structures to enable a varietyof functionalities that are useful to applications ranging from flexible electronics to archi-tected adaptive materials [4]. Instabilities and the ensuing pattern formations, or topologicalchanges, for example, govern phase transitions in materials such as shape memory alloys [5]and attempts have been made at replicating similar principles in structural components.Extensive research has also been devoted to the study of instabilities in thin films on softsubstrates [6, 7], periodic composites [8] and lattices [9, 10, 11]. Relevant examples includethe investigation of global pattern formation in periodic composites and lattices [12], andthe onset of herringbone patterns in compressed thin films [13]. More recently, Bertoldiand Boyce demonstrated how instabilities in soft polymers induce changes in the frequencyband structure of periodic phononic systems [14], while Pal et al. [15, 16] investigated thestatic and dynamic properties of hexagonal lattices and demonstrated that instabilities canlead to the surface confinement of elastic waves. Engineered defect distributions that inducedesired deformation patterns in thin shells have been investigated in [17], where the onsetof instabilities is shown to be associated with the breaking of discrete lattice translationalsymmetry and can be predicted by a phonon stability analysis on a unit cell.Of particular interest to the present study is the investigation of localized deformationsassociated with the onset of instabilities. Methodologies for the prediction and design oflocalized patterns are the objectives of numerous studies and are considered open challengestowards the understanding of failure as well as the engineering of desired interfaces that actas tunable elastic waveguides. The investigation of localization resulting from post-bucklingin lattices and periodic media is presented for example in [15, 18]. Although the generalprinciples leading to localization in continuous media, which manifests as discontinuousstrain distributions, is typically assessed by examining the loss of ellipticity in the Hessianof the strain energy [19], its relation to the microstructure and the effect of the macroscopicgeometry including boundary conditions remain elusive. Indeed, in contrast to a phononstability analysis on a single unit cell for identifying global instabilities, no such recipeexists for localization. Several studies [20, 21] have demonstrated, both numerically andexperimentally, how buckling at the microstructural level evolves into localized deformations.In this context, notable are the works of Papka and Kyriakides [22, 23] who investigated thecrushing of honeycomb cellular lattices under a variety of loading conditions. More recently,d’Avila et al. [24] have demonstrated that the onset of localization in a periodic compositedepends on the effective tangent stiffness of the composite and occurs only if this stiffnessin the loading direction is negative.The current study is motivated by the observation of localized deformation patterns fol-2owing the onset of instabilities in discrete, hexagonal lattices [15]. Such deformations occuras a result of the presence of nonlinearities corresponding to large displacements, but arenot associated with the existence of defects and imperfections, which differs from the typ-ical assumptions made in most studies on localization. Some of the localized deformationpatterns observed in [15] are here presented as background to the investigation presentedherein. Specifically, the paper focuses on progressively simpler lattices with the objectiveof identifying configurations that are characterized by localized deformations and that pos-sibly lend themselves to analytical treatment and may provide useful insight into the mostrelevant parameters that govern the behavior of interest. An overview of the considered lat-tice configurations is shown in Fig. 1, which presents the original in-plane hexagonal lattice(Fig. 1(a)) and an elastically supported square lattice capable of both in-plane and out-of-plane motion (Fig. 1(b)). The dimensionality of the latter problem is further reduced byextracting the single lattice strip shown in Fig. 1(c), which then leads to the elementarycase of the elastically supported one-dimensional (1D) spring mass chain shown in Fig. 1(d).The hexagonal lattice consists of masses connected by longitudinal springs denoted as con-necting lines in the Fig. 1(a), and includes angular springs, denoted by the red arcs, thatprovide a restoring moment proportional to the relative rotations of neighboring springs.In all other figures in Fig. 1, the thin lines describe longitudinal springs, that provide arestoring force aligned with the spring itself and proportional to the relative displacementsof the two nodes it connects. The onset of localized deformation is initially based on obser-vation of the deformed equilibrium configurations resulting from the application of strains,which are evaluated through a Newton-Raphson procedure. For the low dimensional cases ofFig. 1(c), 1(d), a linear bifurcation analysis and the application of a shooting method supportthe interpretation of the observed solutions, the prediction of stability characteristics, theevolution of the equilibrium deformed configurations for increasing levels of applied strainas well as the evaluation of loading path-dependent or hysteresis behaviors. A rich set ofequilibrium solutions which evolve from global patterns, to strongly localized ones that areassociated with folding and unfolding characteristics. These are reminiscent of experimen-tal observations made on thin elastic membranes on soft supports as presented for examplein [7]. This study illustrates how deformed localization occurs as a result of nonlinearitiesassociated with large deformations. It also shows that a systems as simple as a 4 mass 1Dchain can help understand this phenomenon. This suggests that simple configurations asinvestigated herein may provide the basis for the formulations of a general framework forthe investigation of instability-induced localization, along with guidelines for the design ofassemblies with engineered localization patterns.The paper is organized as follows. Following this introduction, Sec. 2 provides an overviewof localization patterns observed in the hexagonal lattice of Fig. 1(a). Next, Sec. 3 presentstwo distinct kinds of instability that arise in square lattices interacting with ground springs:a localized deformation and a globally uniform pattern. Section 4 investigates this localizeddeformation in an analogous one-dimensional lattice using a combination of numerical sim-ulations and bifurcation analysis. Finally, the conclusions of this work are summarized inSec. 5. 3 a)(b)(c) (d)
Figure 1: Summary of discrete lattices considered for observation and analysis of instability-induced localization: hexagonal lattice investigated in [15] (a), square lattice with groundsprings (left) and detail of unit cell (right) (b), one-dimensional square lattice strip (c), and1D spring mass chain (d). 4
Background: instability-induced localization in hexag-onal lattices
Instabilities in lattices leading to localized deformations or global pattern formation dependon the interplay of several factors, including lattice geometry, material parameters, boundaryand loading conditions. The effect of some of these factors is here briefly illustrated asbackground and motivation for the investigations to follow. We consider the quasi-staticbehavior of the hexagonal lattice of Fig. 1(a), which consists of point masses connected bylinear longitudinal springs k a , and includes angular springs k t that oppose the change inangle between neighboring springs [15, 16]. The lattice is constrained to deform in the x, y plane and undergoes large displacements resulting from the imposed set of displacements.We determine the equilibrium configuration of a finite assembly of 32 ×
23 unit cells byminimizing the potential energy of the lattice at each incremental step of the prescribeddisplacement using a Newton-Raphson procedure. This procedure is described in detailin [15]. The behavior of the lattice is controlled by a single nondimensional parameter η = k t /k a L , which quantifies the relative values of the axial and torsional springs, with L being the distance between the two sub-lattice sites in a unit cell. The case of η = 6 × − here considered corresponds to a lattice characterized by a non-convex energy and thatexhibits complex deformation patterns under uniform loading conditions due to instabilities.We first illustrate the lattice deformation corresponding to an imposed vertical displace-ment leading to a total normal strain for the finite lattice of δ y = − .
15. As the imposeddisplacement is progressively increased, the deformation evolves from being initially affine(linear in x and y ) to subsequently achieving global patterns as in Fig. 2(a). These pat-terns, associated with a long wavelength instability, break the translation invariance of thedisplacement field in the x -direction and their onset coincides with with the loss of rank-oneconvexity in the potential energy of a single unit cell [15]. Next, we discuss the case of a bi-axial loading corresponding to imposed displacements along both the horizontal and verticaldirection, leading to normal strains for the lattice domain respectively equal to δ x = − . δ y = 0 .
2. The resulting deformed configuration shown in Fig. 2(b) is characterized bya localized deformation along a “zig-zag” interface. This pattern of localized deformationis attributed to the competition between two factors: localized deformation along the unitvector direction and the imposed affine boundary conditions.Finally, we illustrate how the localization pattern can vary with boundary conditions.Figure 2(c) displays the deformation pattern when periodic boundary condition are imposedalong with a prescribed horizontal normal strain δ x = − .
12. In contrast to the case inFig. 2(b), localization occurs along one of the lattice vector directions, at an angle 2 π/ x -axis. Indeed the change in localization pattern compared to the previous caseis attributed to relaxing the affine displacement constraint on the top and bottom surfacenodes. Thus the deformation patterns that result in these nonlinear lattices depend on The periodic constraint relating a top node x t with its corresponding bottom node x b takes the form x t − x b = (0 , w ). The width w is constant across the lattice, but it is not constrained and allowed to changefrom its initial value. As a result, the lattice expands in the y direction due to Poisson effect a) (b)(c) Figure 2: Deformed configurations having different localization patterns for different bound-ary conditions. (a) Biaxial loading with displacement prescribed on all boundary nodes; (b)compression prescribed on horizontal direction and periodic conditions along vertical bound-ary nodes; (c) compression prescribed on horizontal direction, free to expand and periodicconditions in vertical direction.their geometry, lattice parameters or material properties, loading and boundary conditions.Elucidating the interplay between these factors is a daunting challenge, and is the subjectof ongoing investigations. To address some of the questions that the observations presentedabove pose, we here focus on the effects of lattice parameters under a single, uni-axialloading condition. To this end, we elect to investigate the lattices in Fig. 1(b)-1(d), whichare characterized by a simpler geometry described by orthogonal lattice vectors, but whichhowever have a potentially richer behavior related to the ability of deform both in-plane andout-of-plane. The interaction of these deformations and their relation to the onset of variousinstabilities will be investigated in lattices of increasing simplicity in an attempt to obtaingeneral insight into the unique mechanism observed herein.6 a) (b)
Figure 3: Two kinds of instabilities in square lattices subjected to compression. The defor-mation is in-plane and globally uniform at high stiffness parameter γ = 1 . γ = 0 . Let us consider a square lattice of N × M nodes, that is ( N − × ( M −
1) unit cells, whichincludes point masses at the nodes and three types of springs connecting them. In contrastto the previous hexagonal lattice case where the masses were restricted to move in the xy -plane, here each point mass has 3 translational degrees of freedom and can move in threespatial dimensions. The diagonal springs and the springs connecting nearest neighbors bothhave stiffness k . There is a ground spring of stiffness k g , which resists motion in the out-of-plane z direction. The quasi-static behavior of a lattice of a given dimension is governed bya single nondimensional stiffness parameter, γ = k g /k . In all the subsequent calculationsin this work, we investigate the behavior of a lattice under uniaxial compression. Both thein-plane displacement components are prescribed at all the boundary nodes corresponding tothis uniaxial strain, while the out-of-plane displacement component is free at all the nodes. We first investigate the quasi-static behavior of this lattice under uniaxial compression along x -direction with the strain denoted by δ x . To this end, we conduct numerical simulationsto determine the deformed configuration by increasing the strain incrementally from theundeformed configuration. At each strain level, the equilibrium configuration is obtainedby minimizing the potential energy E of the lattice subject to the appropriate boundary7onditions. The problem may be expressed asminimize X E ( X ) = k M − X m =1 N − X n =1 (cid:2) ( k x m,n − x m +1 ,n k − a ) + ( k x m,n − x m,n +1 k − a ) + (cid:16) k x m,n − x m +1 ,n +1 k − √ a (cid:17) + (cid:16) k x m +1 ,n − x m,n +1 k − √ a (cid:17) (cid:3) + k g M X m =1 N X N =1 z p subject to x q = (1 + δ x ) X q , y q = 0 q = 1 , . . . , N b ,z q = 0 q = 1 , . . . N c . (1)Here N b is the total number of boundary nodes, N c is the nodes on the boundary alongthe y -direction, while ( X p , Y p , Z p ) and x p = ( x p , y p , z p ) denote, respectively, the undeformedand deformed positions of node p and X = ( x , , x , , . . . , x M,N ). Finally, a is the distancebetween nearest neighbor masses in the lattice.Figure 3 illustrates the typical deformation field for two values of stiffness parameter γ . The two chosen γ values are representative of the behavior of lattices with high andlow γ . Figure 3(a) displays the in-plane deformation field when γ = 1 .
0. The out-of-planedisplacement is zero and globally uniform patterns form after the onset of an instability.The displacement field is affine for small strains δ near the undeformed configuration. Thisrange of δ may be determined numerically using a procedure similar to the stability analysispresented later in Sec. 4.2. With increasing strain, this affine solution loses stability leadingto the evolution of globally uniform patterns with a wavelength of two unit cells alongthe loading direction and uniform along the transverse in-plane direction. We remark herethat this deformation mode arises as a bifurcation since the affine solution satisfies theequilibrium equations. Figure 3(b) displays the deformation field for a square lattice with alower stiffness parameter value γ = 0 .
1. Here the deformation happens along the out-of-planedirection too and leads to a localization which consists in one layer of the lattice folding overits adjacent layer. We observe that the deformation field is essentially two dimensional withno displacement in the transverse in-plane direction ( y ). Thus we see that the patterns thatform after the onset of a bifurcation depend on the material property, stiffness parameter γ in this case, and it can range from globally uniform patterns to localized folding behavior. Let us examine the sequence of deformations that result in a localized folding pattern in alattice with γ = 0 .
1. Since there is no displacement in the transverse in-plane direction, weconsider a single strip of the square lattice with N × M = 2 × a) (b) (c)(d) (e) (f) Figure 4: Snapshots of the folding transition process from the fully expanded and undeformedconfiguration (a) to the folded configuration (f), γ = 0 .
1. The process involves two distinctbifurcations, with the vertical deformation first increasing globally and then getting localizedat the center as the lattice is compressed. 9igure 5: Schematic of a 1 D lattice model, with point masses at the nodes interacting withnearest neighbor horizontal springs and vertical ground springs. Each mass has two degreesof freedom and can move in the xz -plane.displacement is maximum at the center and decreases away from it. As the strain increasesfurther (Fig. 4(c)), the out-of-plane displacement increases for the center nodes, while itdecreases for nodes away from the center. The out-of-plane displacement becomes localizedat the two center nodes of the lattice and it decreases rapidly to zero away from these nodes.As the strain increases, the center square folds over (Fig. 4(e)-(f)), thereby resulting in alocalized zone. Note that this localized deformation field will have zero energy at a strainlevel δ = − / ( M − ≈ − .
29 when all the axial springs are unstretched and the out-of-planedisplacement is zero. The existence of multiple zero energy states (undeformed and foldedconfigurations) implies that the potential energy is not quasi-convex. It also points to theexistence of negative stiffness regions and the possibility of this lattice to exhibit hysteresis. Anegative effective stiffness results in a snap-through type behavior as the lattice transitionsfrom the unfolded to the folded configuration or vice versa. We examine this rich set ofbehaviors arising due to this localized solution in the next section using a one-dimensionalmodel which turns out to be quite adequate for our purposes.
To understand the mechanism of the snap-through behavior during localization, let us furthersimplify the square lattice and consider an equivalent 1 D (one-dimensional) model obtainedby taking one section of the square lattice deprived of the diagonal springs. Figure 5 displaysa schematic of the considered lattice model. There are N nodes along a strip indexed from1 to N . The potential energy E of the lattice may be written as E ( X ) = N X i =1 k g z i + N − X i =1 k k x i − x i +1 k − a ) . (2)where x i = ( x i , z i ) and X = ( x , x , . . . , x N ). Here we can assume that the undeformed axialsprings have unit length ( a = 1). Minimizing this energy subject to the boundary conditionsyields the equilibrium configuration. We illustrate the typical displacement response historyas the folding transition happens in this model, see Fig. 6.10igure 6: Sequence of deformed configurations as the lattice transitions from the folded tothe unfolded configuration.We first conduct numerical simulations on a lattice of N = 10 masses by subjectingit to both a compressive strain from the unfolded configuration and a tensile strain fromthe folded configuration. For both cases, we denote the strain by δ , measuring it fromthe reference unfolded configuration. The numerical simulations are again conducted usinga Newton-Raphson solver. This numerical solver based on incrementally applying strainfaces convergence issues in the presence of instabilities and bifurcations. To overcome theselimitations, we develop alternative techniques in Sec. 4.3. This 1 D model also undergoes a folding based localization and we will analyze the defor-mation process leading to this localization in detail. In all our subsequent calculations, weconsider a chain of 10 point masses (nodes) connected by axial springs of stiffness k . Similarto the square lattice, each mass is also connected to the ground by a spring of stiffness k g and the quasi-static behavior is then governed by the non-dimensional stiffness parameter γ = k g /k . The chain starts from x = 0 and the undeformed lattice spans from x = 0to x = 9. This is referred to as the unfolded configuration. The first mass is kept fixedwhile the last mass is subjected to a compressive displacement. The two end masses arealso restricted from moving vertically ( z = z = 0). As the prescribed compressive dis-placement increases, the lattice folds at the center and now spans from x = 0 to x = 7.This is referred to as the folded configuration. Thus we have that δ = 0 for the unfoldedconfiguration while δ = − / (a) -4.5 -3 -1.5 0 1.5 3 4.5-0.3-0.1500.150.3 (b) -0.2 -0.15 -0.1 -0.05 0-0.5-0.2500.250.5 (c) -0.2 -0.15 -0.1 -0.05 0-0.3-0.2-0.100.10.20.3 (d) Figure 7: Vertical displacement of the nodes as the lattice is pulled from the folded configu-ration. Two bifurcation points are observed at low stiffness parameter γ = 0 .
05 (a, c), whilea single bifurcation and a jump in displacement is obtained at high γ = 0 . x = 0 to x = 7 in the axial direction. As the lattice is pulled from both ends,it reaches a final unfolded and undeformed configuration ranging from 0 to 9. The horizontaland vertical axes show, respectively, the compressive strain δ measured from the referenceunfolded configuration and the vertical displacement of all the nodes in between. Figure 7(c)displays the displacements for stiffness parameter γ = 0 .
1. At small displacements from theinitial configuration, the axial displacement field is affine in the lattice and the verticaldisplacement is zero. At around δ = − . δ = − . δ = − . x = 9. We thus see two distinct mode-shapesarising from the onset of two bifurcation points as the lattice deforms from the folded to theunfolded configuration. We address the mechanism of transformation from one mode-shapeto another using a reduced 4-mass model later in Sec. 4.3.2.Let us now examine the displacement when a lattice with a higher stiffness parametervalue γ = 0 . z to z ) as it deforms from the folded to the unfoldedconfiguration. In contrast to the low stiffness parameter case, here the vertical displacementremains zero for larger tensile strains until around δ = − .
1, when the vertical displacementjumps sharply to a non-zero value. The vertical displacement of the adjacent nodes hasopposite sign and resembles the second mode-shape discussed above. As the tensile strainincreases, the vertical displacement decreases and reaches zero at around δ = − . To investigate the local stability of the equilibria computed above, we need to look at theHessian H of the energy E ( X ), in (2). For a chain having N masses it can be expressed as H = I N ⊗ (cid:18) γ (cid:19) + N − X p =1 K h,p , (3)where I N is the N × N identity matrix and ⊗ denotes the tensor product. Here K h,p is thecontribution due to the p -th axial spring having a deformed length r p , with an angle θ p with13espect to the horizontal axis, see Fig. 5. It can be written as K h = D p ⊗ ( R p K loc,p R Tp ) . (4)where the D p is a N × N matrix what all 0 entries but for the ( p, p + 1) principal minor, D p = p p + 10 · · · · · · · · · · · · p · · · − · · · p + 1 0 · · · − · · · · · · · · · · · · · · · , while R p is the rotation of angle θ p and K loc,p is related to the stiffness of the axial springin the local coordinate system aligned with the spring axis: R p = (cid:18) cos θ p − sin θ p sin θ p cos θ p (cid:19) , K loc,p = − r p . (5)It is easy to observe that if all the springs are horizontal, that if θ p = 0 or π for every p , thenthe x and the z part of the Hessian H decouple. More precisely, let P be the permutationmatrix that send X to e X = ( x , x , . . . , x N , z , z , . . . , z N ) = ( X, Z ) then we get f H = P H P T = (cid:18) H x H z (cid:19) . Moreover one can verify that the x part of the Hessian H x is always positive definite. Thusto analyze the stability of an equilibrium with no displacement in the z direction we justneed to look at the eigenvalues of H z .Let us consider a lattice compressed from its initial unfolded configuration to a strain δ .The affine deformation induced by such a strain is given by x p = (1 + δ ) p and z p = 0 forall p . This is clearly an equilibrium configuration. Since the first and last mass are fixed,the stability of this configuration is controlled by the eigenvalues of the ( N − × ( N − H uz = N − α ζ · · · ζ α ζ · · · N − · · · ζ α ζN − · · · ζ α , (6)where α = γ + 2 δ δ , ζ = − δ δ . (7)14 Figure 8: Tensile strain in the axial spring with stiffness parameter γ at the onset of unfoldingfor three different chain lengths. A chain of six masses can essentially capture the onset ofbifurcation from the folded configuration of a much large chain with reasonable accuracy.By looking for eigenvectors X p of the form X p,q = sin( ω p ( q − q = 2 , . . . , N −
1, oneimmediately gets ω p = πp/ ( N −
1) with associated eigenvalue λ p = γ + 2 δ δ (1 − cos( ω p ))and p = 1 , . . . , N −
2. Thus at the critical strain δ ∗ u = − γγ + 2 (1 − cos( ω N − )) λ N − becomes negative and the affine configuration looses stability. The associated modeshape is given by X N − ,q = sin (cid:18) ( N − q − πN − (cid:19) that shows a maximum value at the center which explains the high deformation, which leadsto subsequent folding at the center. The strain value and mode shape predicted by the abovestability analysis is in excellent agreement with the numerical solution.Let us now consider what happens when a folded lattice is stretched. The folded referenceconfiguration is given by z p = 0 for all p while x p = p − p ≤ N/ x p = p − p > N/
2. Thus the strain of this configuration is δ = − / ( N − z p = 0 for all p we need all springs butthe central folded one to have the same strain ∆ while the central one has strain − ∆. Wethus get x p = ( p − p ≤ N/ x p = ( p − − (1 − ∆) if p > N/
2. The15otal strain δ of this configuration is thus δ = ∆ − / ( N − δ , for which this configuration looses stability.The Hessian H fz in this configuration is similar to (6), that is H fz = α ζ · · · ζ α ζ · · · · · · ζ N − α N − ζ N − · · · ζ N − α N − , (8)where α p = γ + 2∆1 + ∆ ζ p = − ∆1 + ∆for all p but for α N = α N +1 = γ + ∆1 + ∆ − ∆1 − ∆ ζ N = ∆1 − ∆ . Let Q be the reflection matrix such that ( Q X ) p = X N − p . Clearly H fz commutes with Q sothat we know all eigenvectors X of H fz must satisfy Q X = ± X . This means that they areeither symmetric or anti-symmetric with respect to the reflection p → N − p . Taking thisinto consideration, we find that the search for the eigenvalues of H fz can be reduced to thatof the eigenvalues of the two ( N/ − × ( N/ −
1) matrices given by H ± = α ζζ α ζ
0. . .0 α ζζ β ± , (9)where α = γ + 2∆1 + ∆ , ζ = − ∆1 + ∆ (10)and β + = γ + ∆1 + ∆ β − = γ + ∆1 + ∆ − − ∆ . Here β + refers to the symmetric eigenvalues while β − to the anti-symmetric ones.As before it is natural to look for eigenvectors X p of the form X p,q = sin( ω p ( q − q = 2 , . . . , N/ −
1. This leads to the condition( β ± − α ) sin( ωM ) = ζ sin( ω ( M + 1))where M = N/ − ωM ) = sin( ω ( M + 1))16o that we get the solutions ω p = πp/ (2 M + 1) with eigenvalue λ p = α + 2 γ cos( ω p ). It iseasy to see that these eigenvalues are always positive for ∆ >
0. In the anti-symmetric casewe get sin( ω ( M + 1)) = 3 + ∆1 − ∆ sin( ωM ) . (11)Observe that since 3 + ∆1 − ∆ > M + 1 M then (11) admits also exactly one complex solution ω r = iν r . The associated eigenvectoris thus X rq = sinh ν r ( q −
1) with eigenvalue λ r = α + 2 γ cosh( ν r ). A complete basis ofeigenvectors is then obtained from the real solution of (11). It is not hard to see that if ω isreal than the associated eigenvalue is positive for ∆ > ν r for finite r . If we take the limit M → ∞ and assume that ν r has a limit, we obtain, after some algebra,cosh( ν r ) = 12 (cid:18) − ∆ + 1 − ∆3 + ∆ (cid:19) = 5 + 2∆ + ∆ − − ∆ so that λ r = ( γ + 4)∆ + 2( γ + 2)∆ − γ ∆ + 2∆ − ∗ = − ( γ + 2) ± p ( γ + 2) − γ + 4 . (12)To evaluate how close the above asymptotic value of ∆ ∗ is to the real value for finite N we computed ∆ ∗ numerically. Since the eigenvalues of H − in (9) are distinct, only one canbecome zero at the onset of the instability and thus the stability of this configuration canbe evaluated by simply examining the sign of the determinant of the matrix H .Figure 8 displays the critical tensile strain ∆ ∗ at the onset of the unfolding bifurcationtransition for a range of stiffness parameter values. Three distinct chain lengths are consid-ered, with N = 4 , γ . Furthermore, this critical value in a finitechain converges rapidly to the corresponding value in an infinite chain. Indeed, the criticalstrain in a chain with 6 masses is quite close to that in the 32 mass chain. Having clarified the folding transition and identified the onset of bifurcation points usinga stability analysis, let us now turn attention to analyzing the transition path between thetwo zero energy configurations. We use two numerical techniques for this purpose. We firstdevelop and use a shooting method based procedure to investigate this folding transition indetail in Sec. 4.3.1. We use this method to illustrate the stability regions as γ increases and17he associated hysteretic behavior. For a moderate number of masses, 10 in our present study,we found this approach hereafter outlined useful in understanding the bifurcation structure.We are aware that this approach may have limitations for large number of masses. Wethen use a classical arc length solver based continuation scheme in Sec. 4.3.2 on a furthersimplified four mass chain which confirms the folding transition. Here we cast the governing equations for the equilibria into the evolution of a dynamicalsystem. Note that the governing equation of an interior node p depends only on its nodallocation and the location of its nearest neighbor nodes ( p − p + 1). The governingequations may be used to solve for the nodal locations of p + 1 node if the other two nodallocations are provided. Accordingly, we may write a relation of the form f ( x p +1 ) = f ( x p , x p − ) (13)where x p = ( x p , z p ). To determine explicit expressions for x p +1 above, as before we let x p +1 − x p = r p cos θ p and z p +1 − z p = r p sin θ p . Figure 5 displays a schematic of a part of thechain with the relevant variables.With this setup, the equilibrium equations for the mass p , obtained by resolving theforces due to the springs in the horizontal and vertical directions, may be written as( r p −
1) cos θ p = ( r p − −
1) cos θ p − , (14a)( r p −
1) sin θ p = ( r p − −
1) sin θ p − + γz p (14b)These equations can be unambiguously solved to obtain r p and θ p as a function of x p − and x p . An explicit formula for x p +1 is then easily obtained using x p +1 = x p + r p cos θ p , z p +1 = z p + r p sin θ p . (15)As discussed in Sec. 4.2, the onset of instability for both the folded and unfolded con-figurations happens along an anti-symmetric mode shape with respect to the reflection Q ,see discussion after (8). We thus expect that the deformation process from the unfoldedto the folded configuration will pass only through anti-symmetric configurations, that isconfigurations for which x − x − p = x p and z − p = − z p . This is well verified in Fig. 7.To find an equilibrium configuration for the chain, we first fix x and as a consequence x . Observe that, due to the symmetry of the configuration, this is equivalent to fixing θ and r that from now on we will call θ and r . Once x and x are fixed, we can use (13) tocompute x up to x and as a consequence we get a full configuration for the chain.At this point there is no need for z to be 0. For a range of value of r we search forthe value of θ for which z = 0. Finally we relate back x to the strain δ . We call thisprocedure the shooting method to find numerically an equilibrium configuration. Observethat high iterations of f in (13) potentially present stability issues. This is the reason whywe decided to compute only anti-symmetric configurations starting from the center of the18 (a) -0.2 -0.15 -0.1 -0.05 000.10.20.30.40.5 (b) -0.2 -0.15 -0.1 -0.05 000.10.20.30.40.5 (c) -0.2 -0.15 -0.1 -0.05 000.10.20.30.40.5 (d) Figure 9: (a) Vertical displacement of the center mass with the axial position of the last massduring the folding process for different stiffness parameter γ values. The Newton-Raphosnmethod solution (black circles) match the shooting method solution (curves) in both loadingand unloading for γ = 0 . γ = 0 . f .To gain further insight into the difference between the high and low stiffness parameterresponse (Fig. 7) and to understand the displacement jump in Fig. 7(d), let us examine thesolution as the center spring folds. Figure 9(a) displays the vertical displacement of the centermass z with compressive strain δ for 4 distinct stiffness parameter values. These curves areobtained from the shooting method solution described above. For low γ , there is a singlesolution for each strain value in the range δ ∈ [ − / , γ increases, there are three possible solutions in a range of δ values. This multiplicityof solutions is evident as the curve bends backward yielding two stable and one unstablesolution for each δ value. By comparison, in Fig. 9(b) we show with black circles the stablesolution computed with the previously described Newton-Raphson approach, as the latticeis stretched from the folded to the unfolded configuration. Obviously, there is in excellentagreement with the shooting method solution (solid curve) over the entire range of thefolding to unfolding transition. Furthermore, the Newton-Raphson method solution of latticecompression from the unfolded to the folded configuration also matches with this solution.Let us now examine the response of a lattice having a higher stiffness parameter ( γ = 0 . θ values and our Newton-Raphson method solution jumps to a stablebranch, which corresponds to a lower θ value. As the displacement increases, the solutionfollows this branch until θ = 0 at around δ = − .
049 after which the vertical displacementof all the nodes is zero. Finally, let us observe the behavior of this lattice as it is compressedto deform from the unfolded to the folded configuration. Figure 9(d) displays the Newton-Raphson method solution for this case. The Newton-Raphson method solution matcheswith the shooting method solution until the onset of instability, when the tangent to the redcurve in the figure becomes vertical. Beyond this point, our Newton-Raphson solver fails toconverge to a stable equilibrium solution. The black circle at around δ = − .
138 indicatesthe stable solution which a physical system is likely to reach as the prescribed compressivestrain increases, which corresponds to a folded configuration and then follows this horizontalbranch (coinciding with the previous unfolding case). Thus we observe how the solutioncan be path dependent for high stiffness parameter γ values when there are multiple stablesolutions.Having provided evidence for the existence of multiple stable solutions, let us see howthe lattice energy and force vary as the lattice deforms from one equilibrium configurationto another. Figure 10(a) displays the lattice energy in the vertical axis as it deforms fromthe folded to the unfolded configuration. The strain δ on the horizontal axis ranges from − / ≃ − .
22 to 0. Three curves corresponding to distinct values of stiffness parameter γ areshown. In each case, the lattice energy is zero at both the folded and unfolded configurations,resulting in a double well potential in this range of deformations. At low γ values, thereis a single energy value for each strain δ , while at high γ , there are multiple energy valuesconsistent with the observations above. For γ = 0 .
2, we observe that the energy curve has20 (a) -0.2 -0.15 -0.1 -0.05 0-0.1-0.0500.050.10.150.20.25 (b)
Figure 10: (a) Lattice energy and (b) force acting on the ends of the lattice as it deformsfor different stiffness parameter values. As the stiffness increases, the force is not unique ina range of displacements and allows for hysteresis during the folding process.three segments and folds back at around δ = − .
138 and − .
1. The middle segment joiningthe two end segments corresponds to the unstable branch.Let us now illustrate how the corresponding horizontal force F acting on the chain variesduring the folding process. It is defined as the derivative of the potential energy with respectto the chain length, given by F = 1 N − ∂E∂δ and it is computed by taking the numerical derivative of the energy illustrated in Fig. 10(a).Figure 10(b) illustrates the force for the same three stiffness parameter values γ . In all cases,the force is zero in the folded and unfolded configurations at δ = − / δ = 0. For lowstiffness parameter values ( γ = 0 .
05 and γ = 0 . δ . The force displacement response has three segments: two linear andone segment in the middle joining them. The linear segments correspond to affine horizontaldeformation about the folded and unfolded configurations. The middle segment correspondsto the rotation of the center spring as θ varies from 0 to π and the effective stiffness of thechain is negative in this regime as the force decreases with increasing prescribed displacementor tensile strain. Furthermore, there are multiple values of force in a range of strains δ forthe γ = 0 . γ , the lattice follows a different path around the unstable branch. For example,in Fig. 10(b), the force jumps down from 0 .
24 to − .
088 while unfolding and it jumps upfrom − .
042 to 0 .
17 while folding. 21 emark 1.
We finally briefly remark on the alternative possibility of inducing hysteresis byexploiting the negative stiffness zone. Instead of prescribing displacement, if we prescribe ahorizontal force at the ends of the chain, then the lattice solution will jump at the end of thefirst segment for all the three stiffness values, for example at a force level F = − . from δ = − . while folding and at a force level F = 0 . from δ = − . while unfolding.Thus our results illustrate the potential for hysteresis and localization guided programmablesmart materials in lattice based media. Here we consider the simplest case of N = 4, that is a chain of 4 point masses anchoredat its extremes. As in previous sections, the unfolded and folded configurations of the fourmasses are (0 , , (1 , , (2 , , (3 ,
0) and (0 , , (1 , , (0 , , (1 , . We can (and will) parametrize the x -coordinate of the last mass as 3 l , with 1 / ≤ l ≤ l by the linear relation δ = l − ∈ [ − / , . Calling x i = ( x i , z i ) the position of the i -th point mass, the energy of the system, ex-pressed in nondimensional form by normalizing it by a reference energy k a / E = γ X i =1 z i + 12 X i =1 ( k x i +1 − x i k − (16)where the first sum is the potential energy of the vertical restoring force and the last sum isthe potential energy of the springs. In this last sum x and x represent the fixed extremeof the chain, namely x = (0 , x = (3 l, N = 4, we assume that the equilibrium position satisfies x = 3 l − x and z = − z . Let us call r the length of the segment from the middle point to the location ofthe second point mass x , and call θ the angle formed between the middle point and thehorizontal. This way, the 4 point masses get located at the points: x = (0 , , x = (3 l/ − r cos θ, r sin θ ) , x = (3 l/ r cos θ, − r sin θ ) , x = (3 l, . Clearly we have only the variables r and θ to consider and the potential to minimizebecomes E ( r, θ ) = γr sin θ + r r − lr cos θ + 9 l − ! + 12 (2 r −
22o that by taking derivatives we get the necessary conditions E r ≡ γr sin θ + (2 r − l cos θ ) − (cid:18) r − lr cos θ + 9 l (cid:19) − ! + 2(2 r −
1) =0 E θ ≡ γr sin θ cos θ + 3 lr sin θ − (cid:18) r − lr cos θ + 9 l (cid:19) − ! =0 . (17)Now, observe that there are two solutions valid for all γ : Unfolded r = l/ θ = 0; Folded r = − l and θ = π .Unlike our previous numerical experiments, here we implement a classic pseudo-arc-lengthcontinuation algorithm starting with the trivial solution branch ( l/ ,
0) for l ranging from l = 1 to l = 1 /
3, and monitoring branch points originating from this branch. Furtherfollowing this bifurcating branch, we will eventually connect the branch ( l/ ,
0) to the othersolution branch ( − l , π ). Stability monitoring is done by looking at the Hessian eigenvalues.In the graphs below, we show what happens for several different values of γ . The mostinsightful visualizations are obtained displaying the values of cos( θ ), with the two trivialbranches corresponding to the values 1 and −
1; when not confusing, we also display theheight value, as in the figures of the previous sections. We will show in bold-face the stableportions of the equilibria branches. In spite of different values of γ , the picture we obtain inthis simple case of N = 4 is quite similar to what we observed in the previous sections for N = 10.For γ small ( γ = 0 . − / , γ grows, the stable connection between the two trivial branches develop a fold andthere are two stable equilibrium solutions for a range of values of the strain: one of thetrivial branches and (portion of) the connecting branch; see the case of γ = 1 .
05 in Figure12(a)-(b).At the same time, as we increase γ further, there is a range of values of the strain whereboth trivial branches are stable, as well as part of the connecting branch of equilibria; see γ = 2 in Figure 13(a). Eventually, for γ sufficiently large ( γ = 7 in Figure 13(b)), theconnecting branch is made up entirely of unstable equilibria and the branch points on thetrivial branches occur outside the range of values of the unfolded/folded configurations. Remark 2.
As an alternative to the continuation technique we used, we also double checkedour results by a more algebraic technique, obtaining the same results. Dividing the second Namely, values where there is another solution branch passing through them. Strain, δ -1-0.500.51 c o s ( ϑ ) (a) -0.8 -0.6 -0.4 -0.2 0 Strain, δ -0.100.10.20.30.40.50.6 H e i g h t o f ce n t e r m a ss , z (b) Figure 11: γ = 0 .
1. Stable equilibria connection between unfolded and folded configurations. -0.8 -0.6 -0.4 -0.2 0
Strain, δ -1-0.500.51 c o s ( ϑ ) -0.2593-0.2662 (a) -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 Strain, δ H e i g h t o f ce n t e r m a ss , z (b) Figure 12: Partially stable equilibria connection between unfolded and folded configurations.Hysteretic behavior. 24
Strain, δ -1-0.500.51 c o s ( ϑ ) -0.4 -0.1264 (a) -0.8 -0.6 -0.4 -0.2 0 0.2 Strain, δ -1-0.500.51 c o s ( ϑ ) -0.7 0.1215 (b) Figure 13: (a) γ = 2. Both folded and unfolded branches are stable for a range of strainvalues. (b) γ = 7. As γ becomes large, the branch of equilibria connecting folded andunfolded states is fully unstable; the states themselves are stable for all strain values. equation in (17) by r sin θ , and with some algebra, (17) rewrite as (2 γr cos θ + 3 l ) (cid:18) r − lr cos θ + 9 l (cid:19) =9 l (2 r − l cos θ ) γr cos θ − l ( γr sin θ + 2 r −
1) =0 . (18) From (18) , the second equation becomes γr cos θ − l (( γ + 2) r −
1) = 0 that is cos θ = 3 l (( γ + 2) r − γr (19) that replaced in the first of (18) gives (3 l (( γ + 2) r −
1) + 3 lr ) (cid:0) γr − l (( γ + 2) r −
1) + 18 γl r (cid:1) − γl r = 0 or l − l (5 γ +16) r + 18 l ( γ +3)(2 γ +7) r − l ( γ +3) ( γ +4) r − γ ( γ +3) r + 4 γ ( γ +3) r = 0 . From the roots of this quintic, which we do by finding them as eigenvalues of the associatedcompanion matrix (and keeping only the real ones for which the relative cos θ given by (19) is in [ − , ), we obtain the possible equilibrium solutions. Conclusions
To understand, predict and control localization patterns in lattices, we considered a squarelattice on an elastic substrate. We exemplified how the type of instability changes fromin-plane to out of plane with increasing stiffness parameter. The latter instability evolvesto a localized deformation with increasing strain and it is investigated in detail using aone-dimensional lattice. Using a shooting method based approach, we identified the stablebranches of deformation as the lattice deforms from one stable configuration to another witha localized deformation. For lattices with low stiffness parameter, there is only one stablesolution as the localized pattern forms, while if this spring stiffness exceeds a critical value,there are two stable solutions in a range of displacements. The lattice energy and forceresponse demonstrate the potential for inducing hysteresis due to the existence of the twostable solutions. Finally, we considered the simplest model possible, 4 point masses with thetwo endpoints fixed. In this case, we used a classical continuation technique and were ableto confirm the transition from unfolded to folded quasi-static equilibrium configurations.In summary, we have illustrated that, in contrast to many other existing works where lo-calization is unpredictable and sensitive to the presence of defects, localized deformation canbe predicted precisely based on geometry and symmetry considerations. Our analysis frame-work may open avenues for controlled design of interfaces for waveguiding and programmablesmart materials applications.
References [1] MH Yoo. Slip, twinning, and fracture in hexagonal close-packed metals.
MetallurgicalTransactions A , 12(3):409–418, 1981.[2] N Murisic, V Hakim, IG Kevrekidis, SY Shvartsman, and B Audoly. From discrete tocontinuum models of three-dimensional deformations in epithelial sheets.
Biophysicaljournal , 109(1):154–163, 2015.[3] JG Ramsay. Development of chevron folds.
Geological Society of America Bulletin ,85(11):1741–1754, 1974.[4] DM Kochmann and K Bertoldi. Exploiting microstructural instabilities in solids andstructures: From metamaterials to structural transitions.
Applied Mechanics Reviews ,69(5):050801, 2017.[5] K Bhattacharya.
Microstructure of martensite: why it forms and how it gives rise tothe shape-memory effect , volume 2. Oxford University Press, 2003.[6] Fan Xu, Michel Potier-Ferry, Salim Belouettar, and Heng Hu. Multiple bifurcationsin wrinkling analysis of thin films on compliant substrates.
International Journal ofNon-Linear Mechanics , 76:203–222, 2015.267] Luka Pocivavsek, Robert Dellsy, Andrew Kern, Sebasti´an Johnson, Binhua Lin,Ka Yee C Lee, and Enrique Cerda. Stress and fold localization in thin elastic membranes.
Science , 320(5878):912–916, 2008.[8] G Geymonat, S M¨uller, and N Triantafyllidis. Homogenization of nonlinearly elasticmaterials, microscopic bifurcation and macroscopic loss of rank-one convexity.
Archivefor rational mechanics and analysis , 122(3):231–290, 1993.[9] CL Kane and TC Lubensky. Topological boundary modes in isostatic lattices.
NaturePhysics , 10(1):39–45, 2014.[10] MS Anderson and FW Williams. Vibration and buckling of general periodic latticestructures. 1984.[11] Jayson Paulose, Anne S Meeussen, and Vincenzo Vitelli. Selective buckling via statesof self-stress in topological metamaterials.
Proceedings of the National Academy ofSciences , 112(25):7639–7644, 2015.[12] N Triantafyllidis and MW Schraad. Onset of failure in aluminum honeycombs undergeneral in-plane loading.
Journal of the Mechanics and Physics of Solids , 46(6):1089–1124, 1998.[13] JW Hutchinson and X Chen. Herringbone buckling patterns of compressed thin filmson compliant substrates.
J. Appl. Mech , 71:597–603, 2004.[14] K Bertoldi and MC Boyce. Mechanically triggered transformations of phononic bandgaps in periodic elastomeric structures.
Physical Review B , 77(5):052105, 2008.[15] RK Pal, M Ruzzene, and JJ Rimoli. A continuum model for nonlinear lattices underlarge deformations.
International Journal of Solids and Structures , 96:300–319, 2016.[16] RK Pal, JJ Rimoli, and M Ruzzene. Effect of large deformation pre-loads on the waveproperties of hexagonal lattices.
Smart Materials and Structures , 25(5):054010, 2016.[17] A Lee, FL Jim´enez, J Marthelot, JW Hutchinson, and PM Reis. The geometric roleof precisely engineered imperfections on the critical buckling load of spherical elasticshells.
Journal of Applied Mechanics , 83(11):111005, 2016.[18] C Combescure, P Henry, and RS Elliott. Post-bifurcation and stability of a finitelystrained hexagonal honeycomb subjected to equi-biaxial in-plane loading.
InternationalJournal of Solids and Structures , 88:296–318, 2016.[19] JW Rudnicki and JR Rice. Conditions for the localization of deformation in pressure-sensitive dilatant materials.
Journal of the Mechanics and Physics of Solids , 23(6):371–394, 1975. 2720] LJ Gibson and MF Ashby.
Cellular solids: structure and properties . Cambridge univer-sity press, 1999.[21] SD Papka and S Kyriakides. In-plane compressive response and crushing of honeycomb.
Journal of the Mechanics and Physics of Solids , 42(10):1499–1532, 1994.[22] SD Papka and S Kyriakides. Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb.
Acta materialia , 46(8):2765–2776, 1998.[23] SD Papka and S Kyriakides. Biaxial crushing of honeycombs:—part 1: Experiments.
International Journal of Solids and Structures , 36(29):4367–4396, 1999.[24] MPS d’Avila, N Triantafyllidis, and G Wen. Localization of deformation and loss ofmacroscopic ellipticity in microstructured solids.