Reprogramming Static Deformation Patterns in Mechanical Metamaterials

This paper discusses an x-braced metamaterial lattice with the unusual property of exhibiting bandgaps in their deformation decay spectrum, and, hence, the capacity for reprogramming deformation patterns. The design of polarizing non-local lattice arising from the scenario of repeated zero eigenvalues of a system transfer matrix is also introduced. We develop a single mode fundamental solution for lattices with multiple degrees of freedom per node in the form of static Raleigh waves. These waves can be blocked at the material boundary when the solution is constructed with the polarization vectors of the bandgap. This single mode solution is used as a basis to build analytical displacement solutions for any applied essential and natural boundary condition. Subsequently, we address the bandgap design, leading to a comprehensive approach for predicting deformation pattern behavior within the interior of an x-braced plane lattice. Overall, we show that the stiffness parameter and unit-cell aspect ratio of the x-braced lattice can be tuned to completely block or filter static boundary deformations, and to reverse the dependence of deformation or strain energy decay parameter on the Raleigh wavenumber, a behavior known as the reverse Saint Venant’s edge effect (RSV). These findings could guide future research in engineering smart materials and structures with interesting functionalities, such as load pattern recognition, strain energy redistribution, and stress alleviation.


Introduction
Metamaterials is a term that has become synonymous with materials with unusual properties that enable exotic or anomalous functionalities, first advocated by Veselago in 1967 [1], in the context of physics of materials with negative permittivity and negative permeability, resulting in a negative refractive index. A validation of this fascinating theory by illustrative works, such as Pendry and Smith [2][3][4], has excited research on applications of photonic metamaterials [5][6][7][8][9], such as wave guiding and imaging resolution enhancement. The theory of negative refractive index has seen similar consequences in sound and wave mechanics when artificial materials are developed with a negative bulk modulus and negative mass density. The ability to control sound waves is therefore realized and the existence of bangaps in material dispersion curves can be utilized for sound and seismic wave insulation and reflection [10][11][12][13].
Lubenski, Bertoldi, and others have also studied topological metamaterials [14][15][16], lattice materials that are on the verge of mechanical instability due to their topological index and so exhibit zero-frequency modes or floppy modes. More interesting is how the release of internal stress in these lattices can lead to localization of those zero modes at a boundary surface without propagation. Such surface qualities allow for applications, such as adhesion and grip enhancement in rubber tires, as well as compliancy and energy absorption in rigid load bearing elements [17,18]. Other studies [19] Materials 2018, 11 on isostatic lattice network have been able to achieve any arbitrary global deformation by periodic actuation of substructure members. A recent study by Karpov [20] on the reverse Saint-Venant edge effect in periodic highly non-local lattice networks has also introduced the concept of deformation decay spectrum, where the presence of bandgaps in such a spectrum leads to blockage of an applied Raleigh mode deformation and the ability to reverse its decay rate. The existence of bandgaps and their effect on static Raleigh wave mode propagation in periodic lattices could be related to acoustic metamaterials and how their bandgap characteristics define dynamical wave propagation. The study [20], even though restricted to 1DoF lattices, is the first to provide an avenue to discuss how arbitrary modes of static deformation could be programmed in a periodic lattice network. Our current interest is to enable programming of a 2D, nonlocal, x-braced lattice to manifest interesting functionalities, such as boundary deformation localization, deformation pattern recognition, strain energy redistribution, and stress and strain alleviation inside the lattice. The use of periodic material systems has showed a great importance in structural and materials engineering, where design and analysis are undertaken at the repetitive unit cell level to provide elegant approaches to cost-effective solutions. Geometrical periodicity has helped in analyzing many interesting micro-structural properties exhibited by metamaterials. Some of the earlier quasistatic analysis of patterned and repetitive structures involved the use of the discrete field analysis as a solution method for the governing system of finite difference equations [21,22]. A compact matrix form of these equations with a discrete convolution operator was also suggested [23], and discrete Fourier transform (DFT) methods were used to write computationally efficient semianalytical solutions for arbitrary force loads in terms of lattice Green's function operators [24]. The Green's function approach also allowed for various domain reduction techniques in molecular mechanics of materials, including single and multilayered graphene and carbon nanotubes [25][26][27][28].
Studies of static deformation in periodic lattices have been done in recent works [29][30][31][32] using several methods, such as the transfer matrix approach for 1D (beam like) structures. Using the transfer matrix approach, these authors [32] realized exponential decay of self-equilibrated end-load components, and rather simple polynomial behavior of tensile and bending displacement fields, even in complex architectural trusses and microstructured beams. Exponential decay of self-equilibrated sets of forces and couple moments at the beam ends was discussed as a manifest of the Saint-Venant principle in a discrete elastic media, and the corresponding decay rates were related to the transfer matrix eigenvalues.
Fully analytical solutions for discrete 2D materials and structures are interesting due to their immense technological significance realized in the recent decades. It was recently shown [20] that with a combination of the Fourier and transfer matrix methods we can produce a fundamental solution as a static surface mode or harmonic, which propagates in the material volume without any shape transformation. Only the amplitude of these modes decays exponentially at a known rate, λ < 1, depending on the Fourier parameter, q ∈ (−π, π), termed the static Raleigh wave solution: Here, the index, n, increases toward the material interior; the index, m, varies along the side edges of a discrete two-dimensional material or structure; and d nm is a vector of all displacement components in any group of repetitive structural nodes numbered (n, m). C(q) is an arbitrary complex amplitude of the Raleigh wave, so that the real and imaginary parts of the solution (Equation (1)) separately represent two possible physical states of static deformation; the decay rate, λ(q), and polarization vector, h, are the eigenvalue and normalized half-eigenvector of a transfer matrix written in terms of partial Fourier images of lattice force constants, as explained in Section 2.
A monotonous increase of λ(q) with q in Equation (1) is a manifest of the Saint-Venant principle in discrete 2D material systems, because q is a basic measure of unevenness of the Raleigh mode (Equation (1)), whose mean square deformation gradient over index, m, is proportional to q 2 . Thus, any acute mode with a higher q should generally decay faster in the material volume than a smooth mode with a smaller q. However, anomalous behavior of the λ(q) dependence is possible, as was recently shown [20], where a group of vertical asymptotes (asymptotic bandgaps), and therefore negative slope intervals, can be introduced. This imply that certain small-scale unevenness may in fact prevail over coarse ones and propagate farer in the material interior. The ability of an engineering structural material to support anomalous propagation of fine-scale details of deformation in the material volume was named the reverse Saint-Venant's effect (RSV) [20], and the material itself-the RSV metamaterial. Another practical significance of the asymptotic bandgap in static Fourier spectra of RSV metamaterials is the ability to completely detain certain types of deformation on the material surface, including some rather smooth ones.
A systematic understanding of static load processing and modification in RSV metamaterials require simple analytic and numerical tools for reconstruction and testing of various deformation patterns in the material volume for any type of practical natural and essential boundary conditions. In this paper, we outline a simple semianalytical approach to reconstruct analytical solutions for any essential (displacement) or natural (forcing) boundary condition applied at one edge of any nonlocal lattice, where the opposite edge is free and indefinitely remote, and periodic boundary conditions are applied at two other opposite edges. A bandgap design approach is also presented to provide a designer the flexibility to program the 2D x-braced lattice load bearing and strain energy storage capabilities, and to engineer lattice materials capable of blockage and polarization of particular Raleigh waves (Equation (1)) on the materials surface.

Displacement Transfer Matrix and Polarization Vectors
We start with the governing equation of equilibrium of an arbitrary nonlocal elastic medium with only pair-wise elastic interactions in the form of non-buckling bars, springs, or linearized interatomic bonds [20,24,[26][27][28]: Vectors, d nm and f nm , are comprised of all displacement and external force components in an arbitrary repetitive group of lattice nodes numbered (n, m). Matrices, k, represent the configuration and intensity of the elastic interactions between nodes of the current group and all its neighbor groups (n', m'). The repetitive group is selected as large enough for the n'-summation to run from n − 1 to n + 1 only, while the m'-summation (along the lattice edge) may run for an arbitrary range. Structural periodicity implies dependence of k only on the differences n-n' and m-m', rather than separate dependences on the current and running indices, which would be the case for non-periodic lattices. Assume that we know an essential boundary condition, d 0m , or a natural boundary condition, f 0m , on the lattice edge and we aim to determine the displacement solution, d nm , where n > 0. This solution will describe the state of free static deformation in the lattice interior, arising in response to a particular boundary condition.
For all n > 0, we may write a homogenous governing equation: Taking the summation over n in Equation (3) from n − 1 to n + 1, we may write: The Fourier domain form of Equation (4) is obtained after performing the discrete Fourier transform (DFT) over its terms to get: where K n (q) = ∑ m k nm e −iqm and d n (q) = C(q)h(q)λ n (q), according to Equation (1). Hence, the displacement transfer matrix, H(q), is formulated from Equation (5) as: The static Raleigh wave mode solution in Equation (1), as mentioned earlier, is dependent on the eigensystem of matrix, H(q), since the size of this matrix is 2R × 2R where R is the number of degrees of freedom.
At a repetitive group of lattice nodes, there will be 2R eigenvalues, λ(q). However, these eigenvalues will come in reciprocal pairs as a consequence of the symplectic [33][34][35] nature of the displacement transfer matrix, H(q). If λ(q) is an eigenvalue of H(q), then 1/λ(q) is also an eigenvalue, and they could be either a real or complex conjugate pair. For a convergent physical solution (Equation (1)), we select only R eigenvalues, such that |λ(q)| ≤ 1. The 2R-component eigenvectors will have , where the bottom half-vector is the top half-vector multiplied by λ(q) [20].
Considering an x-braced lattice with two degrees of freedom per node, see Figure 1, the 4 × 4 transfer matrix, H(q), is consistent with the structure: where k represents the relative stiffness of the diagonal bars over vertical or horizontal bars. The vertical and horizontal bars have the same stiffness. Solving the eigenvalue problem, (H(q) − Iλ(q)) h(q) λ(q)h(q) = 0, the eigenvector is derived Equations (A1)-(A4) in a compact form as: Here, C(q) could be any real or complex number for a given value, q, and this number is chosen to normalize the half-eigenvectors, so that |h(q)| = 1 at all q. It can be seen in Equation (8) that the bottom half-of the eigenvector is equal to the top half-vector multiplied by the eigenvalue, λ(q). The polarization vector for a complex eigenvalue can be derived by substituting the complex conjugate eigenvalue pair, λ(q) = µ ± ωi, into the top half-vector of Equation (8) and regrouping of the real and imaginary parts to get: The polarization vector in the instance of a real eigenvalue, λ(q), is derived directly from the polarization vector, h(q), of the complex eigenvalue by eliminating the imaginary part, ω, in Equation (9) to obtain the form: The static Raleigh wave mode solution in Equation (1), as mentioned earlier, is dependent on the eigensystem of matrix, ( ), since the size of this matrix is 2R × 2R where R is the number of degrees of freedom. At a repetitive group of lattice nodes, there will be 2R eigenvalues, ( ) . However, these eigenvalues will come in reciprocal pairs as a consequence of the symplectic [33][34][35] Here, ( ) could be any real or complex number for a given value, , and this number is chosen to normalize the half-eigenvectors, so that | ( )| = 1 at all . It can be seen in Equation (8)

2DoF Static Raleigh Wave Solution
Since the Raleigh mode solution is being built for a periodic spatial domain along the index, m, of a lattice, the static Raleigh wave solution (Equation (1)) must be constructed such that it is a real-valued cyclic solution by obeying the following symmetry condition about the mid-plane of a periodic domain: Constructing d nm , we first substitute λ(q) and h(q) into Equation (1) and obtain a solution containing real and imaginary parts. Basically, that part of the solution satisfying (11) is considered the real-cyclic solution. Should the eigenvalues be complex-valued, λ(q) = µ ± ωi, it would be convenient to utilize the polar coordinate form, λ(q) = ρe iθ , where ρ = |λ(q)| and θ = Arg(λ(q)). Since eigenvalues are complex conjugate, after substituting them into d nm , we obtain four (4) possible solutions (two real parts and two imaginary parts), although not satisfying the symmetry condition Equation (11). However, we can reduce this to the solution in Equation (12) that satisfies the condition in Equation (11) by summing the two real parts and subtracting the two imaginary parts. Below, we state the real-valued cyclic solution forms Equations (A5)-(A12) for the complex and real eigenvalues as: Complex Eigenvalue: Real Eigenvalue: In Equations (12) and (13), the coefficients (a, b, c, d) reconstruct the polarization vector, h(q), of the Fourier parameter, q, that constructs real-cyclic Raleigh mode solution in Equation (1).

Arbitrary Essential Boundary Condition Solution
Our discussion till now has been concerned with static Raleigh wave solution forms, which are essentially harmonic (cosine, sine), but in practice, applied boundary conditions could be in several forms, such as a gaussian, impact, triangular, etc., and therefore this section is dedicated to finding a solution form that would satisfy an arbitrary essential boundary condition, d 0m . Considering Equation (1), a general solution is obtained by the summation of all possible modes of the decomposed essential boundary condition as: where h 1 (q) h 2 (q) is a matrix of the column-vector components. In Equation (14), C 1 (q) and C 2 (q) are some Fourier coefficients to be determined from boundary conditions. Considering an arbitrary boundary deformation, d 0m , a semi-analytical procedure is developed for finding C 1 and C 2 by multiplying both sides of Equation (14) with the conjugate transpose of a normalized polarization vector, h * 1 (q ), and performing the DFT on both sides of the same equation. Rearranging, Since, Equation (16) can be simplified as: Repeating the above procedure with h * 2 (q ), we obtain a similar expression as in Equation (18) and solve for C 1 (q) and C 2 (q) as: where d 0 (q) = ∑ M−1 m=0 d 0m e −iqm represents the DFT of the arbitrary essential boundary condition, and substituting C 1 (q) and C 2 (q) into Equation (14), a general solution is obtained that can fully represent static deformation in the lattice interior.

Natural Boundary Condition Solution
Having formulated a general solution form for analyzing arbitrary essential boundary conditions, it would be appropriate to deal with scenarios of a natural (forced) boundary condition, which has a greater practical appeal. Considering the Fourier form of the equilibrium governing Equation (4), at n = 0, we can write: where f 0 (q) = ∑ M−1 m=0 f 0m e −iqm represents the arbitrary natural boundary condition and Equation (20) represents the effect of neglecting all the nodal set to the left of the boundary where the natural boundary condition is applied and ignoring boundary stiffness interaction with the same by having the term, 1 2 K 0 (q), in Equation (20). Decomposing d 0 (q) and d 1 (q) from Equation (20) into their Fourier components and solving for C 1 (q) and C 2 (q), we get: Substituting C 1 (q) and C 2 (q) into Equation (14) and adding G(n) = n K −1 (0) −1f 0 (q), a linear polynomial term is constructed to account for uniform deformation [32], after considering the possible canonical modes of static deformation of the Jordan block of the displacement transfer matrix, H(q). Hence, the general solution for any natural boundary condition is stated as:

Raleigh Wave Mode Bandgap Design
Boundary deformation blockage or localization relates to the feature of asymptotic bandgaps in the deformation decay spectrum [20] for a periodic lattice structure, where q corresponds with a zero eigenvalue (λ = 0) and, subsequently, the RSV effect, since η(q) = − log λ(q) starts to decrease in value as we increase q: Growth in the fineness of the static Raleigh wave mode corresponding to a slower deformation decay. The deformation decay spectrum is a map of the distribution of the decay parameter, η(q), over q; examples are shown in Section 8. The aim of this section is to develop a relationship between k and q for finding polarization vectors, h(q), that correspond to a bandgap in the lattice deformation decay spectrum. Knowing that H(q) is a square matrix, it is valid that its determinant is equal to the product of its eigenvalues, detH(q) = ∏ n i=1 λ i (q), according to Vieta's rule and so applying this property, a condition for attaining a zero eigenvalue, λ(q) = 0, could be stated as when detH(q) = 0. Applying this condition, a zero-eigenvalue relationship, Equations (A13) and (A14), is derived for a 2DoF x-braced lattice as: Figure 2 shows a plot of k against q from Equation (23). This plot prescribes the stiffness parameter, k, of a 2DoF x-braced lattice that can be designed as a deformation blocker or filter, and the dark arrows on the plot describe the direction of polarization vectors, h(q), that would be arrested at the boundary surface of the prescribed x-braced lattice. determinant is equal to the product of its eigenvalues, det ( ) = ∏ ( ) , according to Vieta's rule and so applying this property, a condition for attaining a zero eigenvalue, ( ) =0, could be stated as when det ( ) = 0. Applying this condition, a zero-eigenvalue relationship, Equations (A13) and (A14), is derived for a 2DoF x-braced lattice as:   (23). This plot prescribes the stiffness parameter, k, of a 2DoF x-braced lattice that can be designed as a deformation blocker or filter, and the dark arrows on the plot describe the direction of polarization vectors, h(q), that would be arrested at the boundary surface of the prescribed x-braced lattice.
It is also possible to create a bandgap phase diagram for design purposes by introducing an aspect ratio, , as a system parameter Equation (A15). In such a scenario, Equation (23) has the form: 2 1 + + cos + 2 cos + cos = 0 = breadth heigth (24) A plot of against k from Equation (24), as shown in Figure 3, represents a design map for generating bandgaps when modelling a 2DoF x-braced periodic lattice. It also shows the range of k permissible in the design for a specific . A typical x-brace lattice shown in Figure 1, with = 1 , from Figure 3 has k ranging from 0-1.41, which was also seen in Figure 2. Figure 3 also shows that when is between 0 and 0.5, there is no restriction on the parameter, k, for which a bandgap would exist. It is also possible to create a bandgap phase diagram for design purposes by introducing an aspect ratio, α, as a system parameter Equation (A15). In such a scenario, Equation (23) has the form: A plot of α against k from Equation (24), as shown in Figure 3, represents a design map for generating bandgaps when modelling a 2DoF x-braced periodic lattice. It also shows the range of k permissible in the design for a specific α. A typical x-brace lattice shown in Figure 1, with α = 1, from Figure 3 has k ranging from 0-1.41, which was also seen in Figure 2. Figure 3 also shows that when α is between 0 and 0.5, there is no restriction on the parameter, k, for which a bandgap would exist.  In the above section, we dealt with boundary displacement blockage because of bandgap existence in the deformation decay spectrum of a lattice structure, where in the analysis of 2DoF systems, we required only a single = 0 to reprogram a lattice to possess this displacement blocking

Polarizing Structures: Case of Repeated Zero Eigenvalues ( &  → )
In the above section, we dealt with boundary displacement blockage because of bandgap existence in the deformation decay spectrum of a lattice structure, where in the analysis of 2DoF systems, we required only a single = 0 to reprogram a lattice to possess this displacement blocking feature. However, a case of repeated zero eigenvalues, where  1 = 0 and  2 = 0, presents a unique class of structures that could be termed as polarizers, with the mechanical property of reprogramming an arbitrary vector of a Raleigh deformation mode at = 0 into a desired polarization vector, ( ), at = 1 that would be completely blocked by the lattice structure at that point ( = 1). Implementing a numerical searching procedure, it is possible to obtain repeated zero eigenvalues, but for a 2DoF lattice structure, the system transfer matrix, ( ), generates only a single independent eigenvector, { ( ) } , instead of two. Such a structure is deemed to have two (Equation (2)) modes of static deformation according to the Jordan canonical form of ( ) [32], written as: The two possible solution modes, (1) and (2) , resulting from the above Jordan block would have the forms: (1) = 0 (26) (2) = 0 + 0 −1 (27) where is the only independent eigenvector and is the subsequently obtained generalized eigenvector. In the case of a 2DoF system with repeated zero eigenvalues, when the polarization vector, ( ), of the lattice structure is applied, we observe that the mode, (1) , controls the Raleigh mode static deformation, and displacements are blocked at = 0. Applying an arbitrary vector in Equation (1), we note that (2) , a combination of an exponential and polynomial modes, control static deformation in the lattice and displacements are blocked at = 1 instead, as shown below: At = 1, In the above section, we dealt with boundary displacement blockage because of bandgap existence in the deformation decay spectrum of a lattice structure, where in the analysis of 2DoF systems, we required only a single λ = 0 to reprogram a lattice to possess this displacement blocking feature. However, a case of repeated zero eigenvalues, where λ 1 = 0 and λ 2 = 0, presents a unique class of structures that could be termed as polarizers, with the mechanical property of reprogramming an arbitrary vector of a Raleigh deformation mode at n = 0 into a desired polarization vector, h(q), at n = 1 that would be completely blocked by the lattice structure at that point (n = 1). Implementing a numerical searching procedure, it is possible to obtain repeated zero eigenvalues, but for a 2DoF lattice structure, the system transfer matrix, H(q), generates only a single independent eigenvector, h(q)

λ(q)h(q)
, instead of two. Such a structure is deemed to have two (Equation (2)) modes of static deformation according to the Jordan canonical form of H(q) [32], written as: The two possible solution modes, j (1) n and j (2) n , resulting from the above Jordan block would have the forms: j where h is the only independent eigenvector and g is the subsequently obtained generalized eigenvector. In the case of a 2DoF system with repeated zero eigenvalues, when the polarization vector, h(q), of the lattice structure is applied, we observe that the mode, j (1) n , controls the Raleigh mode static deformation, and displacements are blocked at n = 0. Applying an arbitrary vector in Equation (1), we note that j (2) n , a combination of an exponential and polynomial modes, control static deformation in the lattice and displacements are blocked at n = 1 instead, as shown below: At n = 1, j At n = 2, j From the above analysis, we can state that a lattice structure yielding repeated zero eigenvalue polarizes an arbitrary Raleigh mode at the position, n = 1, and this wave completely disappears at n = 2.

Illustrative Examples
We begin our illustrations with several examples of the 2DoF x-braced lattice (Figure 1), with a relative stiffness parameter, k = 0.93, and a deformation decay spectrum of           Deformation configuration (scaled) of the 2DoF x-braced lattice (k = 0.93, q = 1 5 π, m = 0 → 10 , n = 0 → 4 ) given by the analytical solutions (Equations (32) and (33)) in (a,b), respectively. Now, we would like to program an x-braced lattice having a deformation decay spectrum shown in Figure 4 to observe a block in Raleigh mode propagation and the RSV effect. Since the bandgap (λ 2 = 0) in Figure 4 exists at q = 8π 11 ≈ 0.73, the RSV effect must be seen in the subsequent Fourier parameters and so we consider the Fourier parameters, q = 9π 11 and 10π 11 . The Raleigh mode solutions for all three cases are constructed as: : d nm = 0.0015 n 0.6609 cos 8 11 πm 0.7504 sin 8 11 πm : d nm = −0.0756 n 0.4495 cos 9 11 πm 0.8901 sin 9 11 πm : d nm = 0.1039 n 0.9823 cos 10 11 πm 0.1524 sin 10 11 πm (36) The Raleigh mode deformation configurations of Equations (34)-(36) for a lattice vertical dimension, M = 22, are as shown in Figure 7. This figure shows a complete blockage of the Raleigh wave mode at = 8π 11 , as expected, and comparing the finer modes, q = 10π 11 , with the coarser mode, q = 9π 11 , we observe the RSV edge effect: The coarser mode decays faster than a finer one. In Figure 8, the deformation configurations for the Raleigh mode of q = 8π 11 applied at n = 0 are analyzed by varying the stiffness parameter, k. The ability to reprogram the x-braced lattice for blockage is realized as shown in Figure 8 by tuning k to the value corresponding to the band gap, which is k 3 = 0.93. We also observe how the rate of Raleigh mode decay in the lattice interior is programmed by tuning the stiffness parameter, where at k 1 = 0.15, we observe a slow Raleigh mode decay, and at k 2 = 0.6, a much faster decay occurs. Now, we would like to program an x-braced lattice having a deformation decay spectrum shown in Figure 4 to observe a block in Raleigh mode propagation and the RSV effect. Since the bandgap ( = 0) in Figure 4 exists at = ≈ 0.73, the RSV effect must be seen in the subsequent Fourier parameters and so we consider the Fourier parameters, = and . The Raleigh mode solutions for all three cases are constructed as:   , we observe the RSV edge effect: The coarser mode decays faster than a finer one. In Figure 8, the deformation configurations for the Raleigh mode of = applied at = 0 are analyzed by varying the stiffness parameter, . The ability to reprogram the x-braced lattice for blockage is realized as shown in Figure 8 by tuning to the value corresponding to the band gap, which is = 0.93. We also observe how the rate of Raleigh mode decay in the lattice interior is programmed by tuning the stiffness parameter, where at = 0.15, we observe a slow Raleigh mode decay, and at = 0.6, a much faster decay occurs. The distribution of strain energy inside a lattice, when mapped, could show interesting features and a spectrum of such a distribution could help when programming a lattice to harness the functionalities, such as strain energy redistribution and resilience in design. The total strain energy per vertical layer of a unit cell thickness is calculated by summing the strain energy stored at each associate substructure, Equations (A16) and (A17), over the lattice index, m, at a given horizontal index, n, in the x-braced lattice (Figure 1) is plotted against the value, n, for the example in Figure 7. Figure 9 shows that strain energy stored along the lattice index, n, follows similar trends as the decay of deformation along the index, n. Namely, the plot of energy for the mode, q = 8π 11 , exhibits the fastest decay followed by q = 9π 11 and then q = 10π 11 , due to the RSV effect. To present a comprehensive picture to show that the pattern of strain energy decay in a lattice is analogous to its deformation decay, we also plot the strain energy along index, n ( 1 → 4 ), relative to n = 1 for all possible static Raleigh modes in Figure 10. This figure shows a monotonous decay of the strain energy, as the parameter, q, of the Raleigh mode increases from 0 to 7π 11 . When q is 8π 11 , there is a much faster decay due to a bandgap (where λ = 0) and after this point, as we increase q, the normalized strain energy at each lattice index, n, starts to increase, which manifests the RSV behavior.
the stiffness parameter, . The ability to reprogram the x-braced lattice for blockage is realized as shown in Figure 8 by tuning to the value corresponding to the band gap, which is = 0.93. We also observe how the rate of Raleigh mode decay in the lattice interior is programmed by tuning the stiffness parameter, where at = 0.15, we observe a slow Raleigh mode decay, and at = 0.6, a much faster decay occurs. Figure 9. Strain energy along the index, n, in an x-braced lattice with k = 0.93: The fastest total strain energy decay occurs at q = 8 11 π, followed by q = 9 11 π, and then q = 10 11 π due to the reverse Saint-Venant's effect.
Materials 2018, 11, x FOR PEER REVIEW 13 of 19 Figure 9. Strain energy along the index, n, in an x-braced lattice with k = 0.93: The fastest total strain energy decay occurs at = , followed by = , and then = due to the reverse Saint-Venant's effect.
The distribution of strain energy inside a lattice, when mapped, could show interesting features and a spectrum of such a distribution could help when programming a lattice to harness the functionalities, such as strain energy redistribution and resilience in design. The total strain energy per vertical layer of a unit cell thickness is calculated by summing the strain energy stored at each associate substructure, Equations (A16) and (A17), over the lattice index, m, at a given horizontal index, n, in the x-braced lattice (Figure 1) is plotted against the value, n, for the example in Figure 7. Figure 9 shows that strain energy stored along the lattice index, n, follows similar trends as the decay of deformation along the index, n. Namely, the plot of energy for the mode, = , exhibits the fastest decay followed by = and then = , due to the RSV effect. To present a comprehensive picture to show that the pattern of strain energy decay in a lattice is analogous to its deformation decay, we also plot the strain energy along index, n (1 → 4), relative to = 1 for all possible static Raleigh modes in Figure 10. This figure shows a monotonous decay of the strain energy, as the parameter, q, of the Raleigh mode increases from 0 to . When q is , there is a much faster decay due to a bandgap (where  = 0) and after this point, as we increase q, the normalized strain energy at each lattice index, n, starts to increase, which manifests the RSV behavior. Figure 10. Relative strain energy against the Fourier parameter, q. A plot of strain energy at n relative to strain energy at n =1 against the Fourier parameter, q.
At this instance, we illustrate the behavior of a polarizing lattice structure by analyzing two xbraced lattices of equal spatial dimensions; the first lattice with k = 0.4714, and the second lattice with k = 1.0834. The deformation decay spectrums of these two lattices are shown in Figure 11a,b, respectively.
From Figure 11, it can be referenced that at = , the first x-braced lattice (k = 0.4714) has eigenvalues,  >0 and  <0, and the second x-braced lattice (k = 1.0834) produces eigenvalues,  ≈ 0 and  ≈0 . We verify the polarization behavior of the second x-braced lattice, with repeating eigenvalues that are approximately zero, by applying an arbitrary Raleigh mode deformation in Equation (1) or Equations (12)(13). So, instead of using the required polarization vector , with = 0.7677 and = 0.6408, for constructing the solution to the zero eigenvalue, we apply a random At this instance, we illustrate the behavior of a polarizing lattice structure by analyzing two x-braced lattices of equal spatial dimensions; the first lattice with k = 0.4714, and the second lattice with k = 1.0834. The deformation decay spectrums of these two lattices are shown in Figure 11a,b, respectively. that displacements propagate in the first lattice because one of its eigenvalues is not zero. However, in the second lattice, we see that although displacements are not blocked at = 0, they propagate one unit cell distance further and then get completely blocked at = 1. Therefore, the second xbraced lattice can polarize an arbitrary Raleigh wave into a polarized Raleigh wave. The polarization behavior can be checked by simply deriving the required polarization vector from the calculated displacements at = 1 and comparing it with the analytical solution by Equation (12) or (13).  (14), (20), and (21), and in Figure 13, we show the lattice deformation configuration and the normalized strain energy distribution contour map. However, the x-braced lattice in this case experiences a very slow decay in deformation and strain energy compared to its boundary as deformation at a node would be composed of all possible static Raleigh modes. As we move along the lattice, n, it can be realized that the strain energy at the boundary is much more concentrated, and moving away from the boundary, it begins to take a Gaussian form till the energy distribution becomes rather uniform in the material interior.  From Figure 11, it can be referenced that at q = 7 9 π, the first x-braced lattice (k = 0.4714) has eigenvalues, λ 1 > 0 and λ 2 < 0, and the second x-braced lattice (k = 1.0834) produces eigenvalues, λ 1 ≈ 0 and λ 2 ≈ 0. We verify the polarization behavior of the second x-braced lattice, with repeating eigenvalues that are approximately zero, by applying an arbitrary Raleigh mode deformation in Equation (1) or Equations (12)(13). So, instead of using the required polarization vector h, with b = 0.7677 and c = 0.6408, for constructing the solution to the zero eigenvalue, we apply a random polarization, with b = 0.5139 and c = 0.8579, to both lattices and calculate their nodal displacements. The deformed shape of the two x-braced lattices can be seen in Figure 12, showing that displacements propagate in the first lattice because one of its eigenvalues is not zero. However, in the second lattice, we see that although displacements are not blocked at n = 0, they propagate one unit cell distance further and then get completely blocked at n = 1. Therefore, the second x-braced lattice can polarize an arbitrary Raleigh wave into a polarized Raleigh wave. The polarization behavior can be checked by simply deriving the required polarization vector from the calculated displacements at n = 1 and comparing it with the analytical solution by Equation (12) or (13). polarization, with = 0.5139 and = 0.8579 , to both lattices and calculate their nodal displacements. The deformed shape of the two x-braced lattices can be seen in Figure 12, showing that displacements propagate in the first lattice because one of its eigenvalues is not zero. However, in the second lattice, we see that although displacements are not blocked at = 0, they propagate one unit cell distance further and then get completely blocked at = 1. Therefore, the second xbraced lattice can polarize an arbitrary Raleigh wave into a polarized Raleigh wave. The polarization behavior can be checked by simply deriving the required polarization vector from the calculated displacements at = 1 and comparing it with the analytical solution by Equation (12) or (13).  (14), (20), and (21), and in Figure 13, we show the lattice deformation configuration and the normalized strain energy distribution contour map. However, the x-braced lattice in this case experiences a very slow decay in deformation and strain energy compared to its boundary as deformation at a node would be composed of all possible static Raleigh modes. As we move along the lattice, n, it can be realized that the strain energy at the boundary is much more concentrated, and moving away from the boundary, it begins to take a Gaussian form till the energy distribution becomes rather uniform in the material interior.   (14), (20), and (21), and in Figure ??, we show the lattice deformation configuration and the normalized strain energy distribution contour map. However, the x-braced lattice in this case experiences a very slow decay in deformation and strain energy compared to its boundary as deformation at a node would be composed of all possible static Raleigh modes. As we move along the lattice, n, it can be realized that the strain energy at the boundary is much more concentrated, and moving away from the boundary, it begins to take a Gaussian form till the energy distribution becomes rather uniform in the material interior. The examples shown in this paper can be modelled using a commercial finite element analysis software package such as ANSYS, where the top and bottom edges of the lattice are constrained to have only horizontal displacements, rendering a seamless cyclic model [24]. Comparing deformations from the above examples to their model solutions from ANSYS, the difference was of the order, 10 −6 , which shows the high accuracy of our analytical results.

Conclusions
In this paper, we have discussed in detail the 2DoF general Raleigh wave mode solutions for analyzing essential boundary and natural boundary conditions. The concept of bandgap design for deformation blockage and achieving the RSV effect was also introduced. Such bandgap analysis has been shown to be readily applicable to any fundamental Raleigh wave solution due to the solution dependence on a single zero eigenvalue ( = 0) corresponding with the Fourier parameter, q. The bandgap relationships presented can serve as tools for programming the unit-cell aspect ratio and stiffness parameter of an x-braced lattice to block a specific static Raleigh mode or filter out irrelevant modes when an applied boundary condition is a combination of several modes. The case of repeated zero eigenvalues has also been shown to present a unique class of nonlocal lattice that can serve as polarizers to induce blockage at = 1 of an arbitrarily polarized Raleigh wave. Solutions to non-Raleigh mode boundary conditions were shown to depend on all possible Fourier modes with controllable decay parameters, providing an opportunity to program the overall strain energy distribution in the material sample. An equivalent continuum theory [36], analysis of strain energy spectral density, and information entropy of deformation [37] for nonlocal mechanical metamaterials can be an interesting separate study for the future.
An understanding of the methodologies presented in this study can be key in driving future research on RSV metamaterials, where uniqueness in deformation and strain energy distribution patterns are harnessed to design smart materials and structures with interesting functionalities and properties, such as load pattern recognition, high resilience, and stress alleviation. The examples shown in this paper can be modelled using a commercial finite element analysis software package such as ANSYS, where the top and bottom edges of the lattice are constrained to have only horizontal displacements, rendering a seamless cyclic model [24]. Comparing deformations from the above examples to their model solutions from ANSYS, the difference was of the order, 10 −6 , which shows the high accuracy of our analytical results.

Conclusions
In this paper, we have discussed in detail the 2DoF general Raleigh wave mode solutions for analyzing essential boundary and natural boundary conditions. The concept of bandgap design for deformation blockage and achieving the RSV effect was also introduced. Such bandgap analysis has been shown to be readily applicable to any fundamental Raleigh wave solution due to the solution dependence on a single zero eigenvalue (λ = 0) corresponding with the Fourier parameter, q. The bandgap relationships presented can serve as tools for programming the unit-cell aspect ratio and stiffness parameter of an x-braced lattice to block a specific static Raleigh mode or filter out irrelevant modes when an applied boundary condition is a combination of several modes. The case of repeated zero eigenvalues has also been shown to present a unique class of nonlocal lattice that can serve as polarizers to induce blockage at n = 1 of an arbitrarily polarized Raleigh wave. Solutions to non-Raleigh mode boundary conditions were shown to depend on all possible Fourier modes with controllable decay parameters, providing an opportunity to program the overall strain energy distribution in the material sample. An equivalent continuum theory [36], analysis of strain energy spectral density, and information entropy of deformation [37] for nonlocal mechanical metamaterials can be an interesting separate study for the future.
An understanding of the methodologies presented in this study can be key in driving future research on RSV metamaterials, where uniqueness in deformation and strain energy distribution patterns are harnessed to design smart materials and structures with interesting functionalities and properties, such as load pattern recognition, high resilience, and stress alleviation.
Author Contributions: The authors contributed equally to the paper content. Writing-original draft, L.D.; writing-review & editing, E.K.
Funding: This work is supported by the U.S. National Science Foundation via Grant #1634577.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
Finding eigenvectors of the transfer matrix, H(q), of 2DoF x-braced lattice: To find eigenvectors, h(q) , of H (q, we use the row reduction echelon method Writing down equations, z is 0, which makes the eigenvector zero, but since an eigenvector cannot be zero, we take z as any real or complex number. However, taking z = 1 and solving for (x, y, w, z) we get: The above expression can be rewritten as: Constructing real-valued cyclic Raleigh wave solutions: