Reduced Order Modeling of Dynamic Mechanical Metamaterials for Analysis of Infinite and Finite Systems

Dynamic mechanical metamaterials (MMs) feature sub-wavelength micro-structures that interact with stress waves, exhibiting exotic functionalities. Numerous exciting potentials have been proposed for MMs, including wave attenuation [1–3], negative refraction [4–7], cloaking [8–11], and insulators [12,13]. The systematic design of metamaterials requires a comprehensive understanding of the dynamic behaviors of MM itself as well as manufacturing design constraints. The first natural step is to find the eigenfrequency band structure and mode shapes of the design. The former encodes major characteristic information of a MM unit cell such as resonance frequencies, wave velocities, and band gaps, while the latter is needed to determine scattering in finite structures and to evaluate interactions between modes fully. Common methodologies of band structure calculation include plane wave expansion [14–17], transfer matrix method [18–20], and finite element method (FEM) [21]. In almost all cases, the design and analysis of these dynamic systems are plagued by geometric complexity and computational burden. The major challenges include (1) unclear design-performance relationship and (2) expensive computational cost for calculating the dynamic response. Therefore, for design optimization purposes, a reduced order modeling (ROM) method is clearly more suitable as it allows for simple and fast computation limited to frequency range or modalities of interest.

In this paper, we present a ROM approach for fast computation of MM problems, which can be used for different study setups, including eigenfrequency band calculation and computation of time-dependent dynamic responses of finite structures. The metamaterials considered here are assumed to comprise 2D micro-structural designs with beam-like elements. The materials are assumed to have no loss or gain mechanisms, but the inclusion of linear viscoelastic response can be considered as a natural future expansion. Detailed numerical and experimental studies on similar micro-structures can be found in previous works [22,23].

Extensive reduced order modeling techniques have been developed for vibration problems, e.g., dynamic condensation [24], improved reduced system [25], and system equivalent reduction expansion process (SEREP) [26]. These reduction methods in general employ certain transformation matrices that map the full set of degrees-of-freedom (DOF) to a reduced set of DOF. For wave propagation problems, especially metamaterial problems, the existing model order reduction methods are limited and less applicable because the system matrices and the eigenfunctions are dependent on the wavevector. The wavevector dependence leads to the frequency and mode variation in the band structures and is a key element in metamaterial dynamics. Therefore, novel reduction schemes that can preserve the wavevector dependence and band accuracy are needed for metamaterial problems. To this end, Hussein [27] introduced reduced Bloch mode expansion (RBME) for fast computation of band structures. The RBME method employs selected Bloch eigenfunctions to reduce the dimensionality. A similar method, Bloch mode synthesis (BMS) [28,29], is an extended sub-structuring technique that describes the structural DOF by normal and constraint modes. Both the RBME and BMS methods utilize selected eigenmodes to construct transformation matrices that reduce the size of the full matrices. These transformation-based methods effectively reduce the number of equations, but the resulting matrices are no longer representing the physical quantities (stiffness and inertia), therefore less suitable for geometric or material design problems. Additionally, these methods could not be applied to time/frequency domain computations of finite arrays. Nevertheless, these methods have been shown useful for topology optimization [30] in terms of reducing the computational cost. An alternative scheme is to develop discrete models comprised of masses and springs. The discretized mass-spring representation has been widely accepted in the literature as it offers analytical formulations that simplify the computational effort while retaining essential physics. It has proven to be beneficial for various design aspects such as feasibility analysis [31], reliability assessment [32], and design space mapping [33]. For higher-order systems operating at high-frequency ranges, an excellent example of modeling discrete weakly coupled MMs has been introduced by Matlack et al. [13], where the model reduction is performed using the Schrieffer–Wolff transformation so that the modes in the frequency range of interest are decoupled. However, this method could only be applied to narrow-band dynamics.

While the mass-spring representation can significantly reduce the computational effort, certain vibration modes may exhibit mixed coupling between DOF. To accurately capture such dynamics, the elastic spring elements cannot be simple 2DOF elements. In our previous work [23], the elastic spring elements are physically represented as beams. The reduced stiffness and mass matrices can then be derived using simple strength of materials analysis. Such an approach provides analytical matrices that operate on physical DOF and is naturally suitable for tuning the response via control of physical dimensions and material choices [33]. They also make interpreting the modal physics straightforward. For example, such a beam-based discrete model allows for accurate identification of the level repulsion [34] and coupling between the DOF. However, the selection of DOF has been mostly a heuristic step that may affect the results. In addition, approximating the structural components as standard beam elements does not generally match the actual response of the beam-like elements as accurately as needed.

The present work introduces a systematic implementation of the structural-element-based ROM approach that overcomes these limitations. A generalized ROM procedure, parameterized in terms of effective structural stiffness parameters and discrete DOF inertia, is developed and is shown to be applicable to a large family of 3D printable MM designs. The conceptual idea of the proposed ROM method takes advantage of the fact that the 3D printable MMs operating as low frequency (long wavelength) locally resonant systems are often comprised of slender plate- or beam-like elements. In addition, in most cases, only the low frequency dynamics of the MM are of particular interest for practical applications, which reside in the subspace spanned by the lowest few eigenmodes, representable using a few carefully selected DOF. Modeling the system with these “master” DOF can reduce the computational cost while maintaining high fidelity of the underlying physics. A structural assembly system with symbolic matrices is used to represent the repeating unit cell (RUC). By optimizing the energy fitness compared with numerical results, one can find the effective stiffness and inertia parameters of this structural assembly. Such a ROM unit cell can accurately predict the eigenfrequency band structures with minimal computational effort. This approach improves upon existing model order reduction methods in handling problems in metamaterials by maintaining eigen-solution accuracy within the Brillouin zone and providing parameterized matrices. It incorporates the propagating nature of waves, rather than just the modal response of finite structures. In MM systems, small variations in geometry can drastically change the overall response. Using an analytical model that characterizes the MM with a small number of parameters is therefore advantageous for understanding the influence of each component and fine-tuning the design. These resulting ROMs can also be extended for modeling finite-sized arrays, and the reduction in DOF can accelerate the computation process significantly, especially in time-dependent problems.

This paper is organized as follows. The general procedure of ROM development is first introduced in Sec. 2. Then, two examples are given in Sec. 3, showing the accurately reproduced band structures. Finally, further uses of the proposed ROM approach are discussed in Sec. 4, including time domain and impact modeling of finite structures and tractable study of exceptional points and level of repulsion. The conclusions and future outlook of this work are discussed in Sec. 5, where it is emphasized that the proposed approach will lead to efficient modeling and design discovery of mechanical metamaterials, with accurate construction of cellular discrete models.

The full procedure of ROM construction is as follows:

1. Assign a set (n_p) of primary nodes in the continuum unit cell, from whose displacement and rotation values the full continuum mode shapes may be approximated;

2. Perform eigenfrequency simulations at a few (n_k) selected wavevectors to determine the frequency ω^c, down-sampled continuum mode shapes Φ^p, and the modal energies $T¯c$⁠, $V¯c$ for the lowest n_m modes. The superscript ^c indicates the global quantities measured from the continuum calculations, and the superscript ^p denotes the down-sampled mode shapes;

3. Construct the symbolic stiffness K^f(β) and mass M^f(μ) matrices based on the unit cell geometry (layout, connectivity, symmetry) and positions of the n_p primary nodes as well as those of the additional n_d dependent ones (to be eliminated based on Floquet periodicity);

4. Apply Floquet boundary condition to the symbolic matrices so that the equations of motion only involve the DOF at the n_p primary nodes, yielding the symbolic matrices K^p(β, k) and M^p(μ);

5. Find the effective mass matrix M^p (i.e., parameters μ) by optimizing the kinetic energy fitness (matching the ROM results with the continuum $T¯c$⁠), for the lowest n_m modes, using the measured ω^c and Φ^p;

6. Identify the slave DOF (whose contribution to kinetic energy is negligible) and perform static condensation so that the number of studied DOF is reduced (from 3n_p to n_r). This step establishes the numerical-valued matrix M^r (a sub-matrix of M^p), the symbolic matrix K^r(β, k), and reduces the measured modes from Φ^p to Φ^r;

7. Find the effective stiffness parameters β by optimizing the potential energy fitness (matching the ROM modal matrix $V¯r$ with the already determined $T¯r$⁠), for the lowest n_m modes, using ω^c and Φ^r and the established ROM mass matrix M^r;

8. Use the established matrices K^r and M^r to compute the band structure, or adjust them for other types of problems.

In the next section of this paper, the detailed ROM construction steps are demonstrated through examples.

To demonstrate the ROM parameterization process, two MM unit cells are selected. The square unit cell in Fig. 1(a) features an H-shaped resonator mass, while the hexagonal cell in Fig. 1(b) has two split resonators. The two RUCs are modeled in FEM using the same material (typical alumina), with Young’s modulus E = 300 GPa, Poisson’s ratio ν = 0.22, and density ρ = 3900 kg/m³. Both designs have the lattice constant a = 10 mm. One should first assign a set of n_p nodes in the continuum model, where the deformation will be sampled and used for mode shape projections. These sampling points should, in principle, be located at the mass centers of structural components, or at the intersections between beam elements. Further detailed analysis and principles of DOF selection are given in Refs. [35,36]. The chosen sampling nodes in the two examples are denoted by the dots in Fig. 1. Notice that there is no node assigned at the top or right edge of the cell frames, due to the known Floquet periodicity Eq. (1) for infinite arrays. The mode shapes Φ^p will be represented by the deformation vector u^p containing the displacement u_x, u_y in x − y plane and the rotation θ_z at the chosen nodes. The rotational component is derived from the curl of the continuum displacement field $θz=(∂uy∂x−∂ux∂y)/2$⁠.

Then, one performs finite element simulations at a few selected k points in the irreducible Brillouin zone (IBZ). Preferably, these k points should be far away from each other. In the shown examples here, we select n_k = 4 points at k = [0, 0], [π/a, 0], [0, π/a], [π/a, π/a]. At each wavevector point, one collects the eigenfrequency results for the lowest n_m modes. The number n_m is determined such that: (1) the n_mth frequency $fnm$ covers the frequency range of interest and (2) the locations associated with dominant deformations for the lowest n_m modes are included in the pre-selected nodes. We chose n_m = 6 for the square cell and n_m = 8 for the hexagonal one. In this step, the collected information includes the diagonal frequency matrix $ωc∈Rnm×nm$⁠, the diagonal kinetic energy matrix $T¯c∈Rnm×nm$⁠, the diagonal potential energy matrix $V¯c∈Rnm×nm$⁠, and the mode shape matrix $Φp∈C3np×nm$⁠.

Taking advantage of the geometrical symmetries in the continuum model, the number of unknown parameters in the ROM property matrices can be reduced. For example, beam 5-6 and beam 7-6 are symmetric with respect to node 6 in Fig. 2(a). Then, node 5 will have the same mass and rotational inertia as node 7. The two beams will share the same local stiffness matrix as well (in the un-rotated local coordinate system). Under this idealized description, the structural symmetries lead to degeneracy in the eigenvalues. However, any physical or numerical realization of such systems will have the tendency to become non-degenerate due to any small asymmetry. The benefits of ROM for understanding the physics of modal degeneracy will be discussed in Sec. 4.1 and Appendix  B.

With the effective stiffness β and inertia μ parameters determined, the ROM procedure is completed. The optimized modal energy fitness ensures the fidelity of the ROM. It is observed that the optimized matching of modal relation leads to the accurate reproduction of the eigenfrequencies as well as the mode shapes at the pre-calculated n_k wavevector locations. Beyond these locations, the eigen-analysis results are also well extrapolated because of the symbolic implementation of the structural stiffness and Bloch–Floquet periodicity. One can solve for the eigenfrequency and mode shapes through the analytical formulation Eq. (2) for any given wavevector k. Such computation will be extremely fast due to the compactness of the matrices. Figure 5 shows the eigenfrequency band structures for the two studied examples, plotted in the dimensionless wavenumber space Q_x,y = k_x,ya, where a = 10 mm is the lattice constant. The surfaces are generated based on ROM, while the dots represent FEM results. It can be seen that the ROM provides close approximations of the band structures. Notice that the ROM construction only requires the simulations at four different wavevector locations [Q_x, Q_y] = [0, 0], [0, π], [π, 0], and [π, π]. The number and locations of these input simulations are, however, not fixed to the given ones.

It should be noted that the minimizing the matching error e_(j)(β) → 0 in Eq. (20) is a necessary yet insufficient requirement for the ROM system $Kr(β),Mr$ to produce the exact eigen-solutions $Φ(j)r,ω(j)c$ obtained from FEM simulations. The left multiplication of $Φ(j)r$ in Eq. (20) reduces the number of equations since the mode shape matrix is rectangular. However, the number of unknowns in β is limited, and Eq. (20) is collected for multiple wavenumber points. Therefore, such an optimization scheme creates an over-determined problem for seeking the limited set of ROM parameters representative of the unit cell properties. With known eigen-states, the ROM produces the same modal energy matrices as the higher-order FEM system. Then, it is observed that the resulting ROM leads to eigenvectors that closely agree with the FEM results. An alternative way to identify β is to create a multi-objective optimization problem in which one also optimizes the ROM mode shape accuracy while minimizing the modal energy error given by Eq. (20). In practice and in the examples here, the secondary optimization objective of maintaining mode shape accuracy is omitted and only used as a sanity check, leading to significant computational cost savings but insignificant or no loss to accuracy. The latter outcome is understood to be a consequence of the symbolic development of dynamic matrices based on the internal cell topology (beam connectivity).

This approach advances the well-established model order reduction methods such as SEREP [26] in attacking MM problems in the sense that (1) the proposed ROM maintains the eigen-solution accuracy for any wavevector and (2) it provides analytical and parameterized matrices instead of numerical ones. It expands such methods by including the propagating nature of waves instead of modal response of finite structures. In MM systems, the micro-structural features play vital roles in the dynamic properties. A small variation in the geometry could lead to a drastic change in the overall response. Therefore, an analytical model with parameterized structural elements is particularly advantageous for understanding the influence of each component and fine-tuning the design. Several applications of the method are discussed in the next section.

The proposed ROM approach has a wide application spectrum, as the matrices are parameterized by the physical properties (structural stiffness and inertia) and the modeled DOF are physical deformations instead of generalized coordinates. Therefore, the developed ROMs preserve the necessary physical ingredients for further analysis. The dependence of ROM parameters on the geometric dimensions and material properties may be curve-fitted for design purposes with continuous functions, see Ref. [33]. The fidelity of these models is inherently guaranteed by the optimized energy relations and the accurate production of band structures and mode shapes. While the ROM is capable of generating the band structure accurately and efficiently, the band computation is not the ultimate goal of the ROM, rather it is the basis and a starting point.

One immediate application is optimizing unit cell designs for desired eigenfrequencies. The ROM characterizes the continuum unit cell with a finite number of stiffness and inertia parameters and provides the analytical formulation of the eigenfrequency bands. It is then intuitive to tune the structural parameters and associated geometry for desired eigenfrequency performance (wave speeds and band gaps) based on the analytical model. The detailed steps are omitted here. A simplified example can be found in the previous work [33]. Other application examples are discussed in the following.

The micro-geometry and periodicity of MMs add an extra layer of complexity to the analysis of band topology and scattering response. The high dimensionality of traditional models presents challenges in understanding physical phenomena and interpreting results. In MM and phononic band structures, many apparent crossing points may exist between eigenfrequency branches. It is important to classify these crossings as either degeneracy points (real crossings) or level repulsions (avoided crossings) [34]. Despite a small quantitative difference between the two types of crossings, this discrepancy can lead to misunderstandings of modal natures and scattering responses. Level repulsion indicates mixing of modes, caused by the coupling between DOF, resulting in unexpected energy transfer in scattering analysis. Conversely, real crossings indicate fully decoupled modes [23].

A demonstration of band identification can be seen in Figs. 6–7, where the 1D band structure along the Γ − X direction is analyzed for the previously shown square cell. Figure 6 shows good overall agreement between the ROM and FEM results. However, zooming into a region with an apparent crossing (indicated by the circle in Fig. 6) reveals a discrepancy. The FEM results with default meshing, shown as the dashed curves in Fig. 7(a), clearly indicate that the two relevant branches appear repulsed with each other. On the other hand, the ROM results, shown by the dashed curves in Fig. 7(b), indicate a real crossing between the two branches.

To confirm that the real crossing indicated by ROM is a correct observation, we calculate the geometric phase along a prescribed wavenumber path in the complex domain, see Appendix  B for details. The computed geometric phase is zero, suggesting a real crossing point indeed. It is noticed that the discrete model possesses the same symmetry group as the idealized continuum one, while the FEM model may not have such symmetric properties due to the mesh imperfection. To further investigate the source of such a discrepancy, a manual perturbation to the parameterized ROM quantities is performed to break the two-fold symmetry of the ROM system. Then, repulsed branches are found in the symmetry-broken ROM band structure, as shown by the solid curves in Fig. 7(b), similar to the FEM results with default meshing. In this case, an exact geometric phase of π is accumulated after two loops in the prescribed wavenumber path, indicating the existence of exceptional point [34,39] in the complex domain, which is a known companion of level repulsion. Such an analysis using geometric phase calculations is rather easily implemented with the ROM formulation, but it can be extremely challenging for FEM because of the need for evaluating high-dimensional eigenvectors of large-sized non-Hermitian matrices. As the computational complexity of eigen-problems is of the order $O(n3)$⁠, the cost of ROM is significantly reduced (for both the band computation and the topological invariant computation).

The ROM analysis with perturbation in symmetry reveals the source of discrepancy, i.e., the root of this apparent repulsion is not associated with the cell response, but rather due to asymmetry in the numerical mesh. We then validate this conclusion by adjusting the FEM mesh with enforced two-fold symmetry. The resulting corrected FEM band, shown by the solid curves in Fig. 7(a), exhibits a real crossing point, which is the correct topology for an idealized symmetric model. This identification task requires repeated simulations and mesh refinements with the FEM approach. On the other hand, the analytical representation of the ROM formulas allows for easy distinction [23,39] between real crossings and level repulsions. This analysis shows that the studied crossing point is a symmetry-protected degeneracy (rather than an accidental one), which is unstable and prone to forming repulsed branches due to imperfections in the material or geometry of an actual specimen. The same is true for continuum FE models, which may lack symmetry due to their meshes.

Correct identification of these band topologies is crucial for understanding the dynamic behaviors in the scattering response, especially with the presence of the potential symmetry breaking that may lead to further misunderstandings. The ROM approach is more efficient for band sorting purposes (and can easily be leveraged in the calculation of topological invariants based on path integrals), and it provides benefits in understanding the difference between an ideal model and a realistic one, as well as in distinguishing the “normal” and “accidental” degeneracy.

The developed ROM matrices also allow for computation of the equi-frequency contours k_x(ω, k_y) or k_y(ω, k_x), i.e., to find the wavevector solutions k_x or k_y at prescribed k_y or k_x and frequency values. For complex k components, the solved eigen-modes are the propagating and evanescent waves that constitute the basis solution for the oblique scattering problem. Such problems are rarely studied for metamaterials but are of prime importance for analyzing the scattering dynamics [40]. Solving this type of problem in FEM is currently impractical (but not impossible) because the real-valued frequency ω can not be easily enforced in the traditional eigenfrequency study for given complex k components. It is also generally not straightforward to assign a wavenumber component as an eigenvalue to be found, though for phononic media an elegant mixed eigenvalue approach is presented in Ref. [40]. For metamaterials with complex internal features, the proposed ROM approach would be an ideal alternative. Using the ROM matrices, this problem can be simply solved by finding the global minima of the determinant det[K^r(k_x, k_y) − ω²M^r] in the complex (k_x, k_y) space for prescribed ω values. Examples of such contour calculations are deferred to focused studies on their use.

It must be noted that the modeling of finite arrays requires an extra step in terms of characterizing the boundary elements. The structural components of the main body remain the same as the ones obtained from the unit cell ROM. However, the outermost elements need to be quantified separately due to the different boundary conditions assigned to them, especially when the MM array is in contact with a different homogeneous medium, or exposed to traction and/or displacement boundary conditions. The detailed approach is subject of current research and is initially carried out by optimizing the static response or impedance of edge cells with the help of a few FEM simulations [42,44].

A general reduced order modeling technique for periodic mechanical metamaterials is introduced. The method uses a limited number of simulations at selected wavevector locations to establish the reduced system matrices, which are parameterized based on structural connectivity. Effective parameters are extracted by matching modal energies. This approach expands upon previous model order reduction techniques by considering the wave’s propagating nature, leading to accurate eigen-solutions for any wavevector. The parameterized and analytical matrices generated by the ROM method offer valuable insights into the micro-structural influence on overall structural behavior and provide significant assistance in fine-tuning of design. Additionally, the ROM approach leads to fast computation of the dynamic responses of finite-sized arrays, with a significant reduction in computational effort compared to FEM.

To summarize, the highlights of this work are: (1) the proposed ROM method can be easily applied to any periodic metamaterial that has beam-like components; (2) the reduced order matrices allow fast and accurate computation of band structures and the dynamic responses in frequency and time domains; and (3) the ROM method can further benefit design optimization due to its computational efficiency.

The essential limitation of the proposed work is that the ROM method can only be applied to MM micro-structures comprised of beam-like elements (and potentially plate-like elements). Other types of micro-structures, such as layered media, or unit cells with solid inclusions in a solid matrix, are not the suitable modeling target for the proposed ROM (instead, one can use RBME [27]). In addition, the proposed ROM approach is only applicable to systems with stiffness coupling between nearest neighbors, i.e., the long-range interaction [45] is not considered.

The presented method could contribute to the modeling and design of finite and periodic mechanical metamaterials by reducing the computational effort. The micro-structure and periodicity of MMs lead to exciting dynamic properties and present theoretical questions in the physics of micro-structured media. The ROM method, equipped with the vibration and strength of material domain knowledge, can offer concise descriptions of the micro-structural dynamics and is a solid analytical tool for studying metamaterial dynamics. In addition, the presented method leads to significant improvements in computational efficiency and is a promising candidate for further design optimization of graded MM arrays with data-driven techniques.

An immediate topic of research is the handling of edge cells due to their different dynamic response, while the interior cells appear to be very well represented by the ROM based on infinitely periodic media. Future work to be implemented is to adjust the element stiffness matrix formulation Eq. (A2) for compatibility with 3D (and potentially composite) beam and plate elements. In addition, the optimization approaches for matching global quantities between the FEM and ROM can be adjusted, so that lossy elements (viscoelastic material) are allowed in the system. These future potentials would extend the modeling capability to 3D designs and enlarge the feasible design space.

In terms of theoretical advancements, it is shown that the ROM representation allows for the efficient identification of band topology, as the geometric phase can be computed with the small-sized matrices and minimal effort. As a future direction, we suggest further truncating the ROM system to a second-order one, to approximate the local topology near those (apparent) crossing regions of the band structure. Then, only the two relevant modes and their associated DOF are involved. The truncated ROM matrices would lead to a simple yet elegant representation for analytical investigation of the band topology, e.g., see Ref. [39]. It is desired to use the derived mode shapes near such locations as the basis for high-sensitivity detection.

The authors wish to thank US Army Research Laboratory for continued support throughout this effort. This research was supported by DEVCOM ARL through Cooperative Agreements W911NF-17-2-0173 and W911NF-20-2-0147.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

The eigenfrequency trajectories of the parametric evolution in wavenumber Q are shown in Fig. 11. We first consider the case of evolution near the real crossing zone with the symmetric ROM setup, as shown in Fig. 11(a). Both two initial states, denoted by $■$ and ♦, are completely recovered after one cycle. The geometric phase picked up by either of the trajectories is exactly zero. No mode-switching behavior can be found. When the paths cross each other in f space, ensuring that there is no discontinuity in the mode shape will select the right choice of path.

For the case of asymmetric ROM with level repulsion, the trajectories are shown in Fig. 11(b). The state evolution continuously follows the complex f trajectory. Due to the existence of an exceptional point [34] and the Riemann sheet structure of the eigenfrequency surfaces, the eigen-mode switches after one cycle instead of recovering back to the original mode, i.e., $■→⧫$♦ and vice versa. After two cycles, the state goes back to the initial mode with an accumulated extra geometric phase of π: {▪, ♦}→{– ▪, –♦}. To fully restore the initial starting mode, four cycles will be needed. Similar observations have been made in a micro-cavity experiment in which an EP is encircled [49]. We refer to Refs. [46,49] for detailed theory and Ref. [50] for an experimental study of dynamically encircling an exceptional point.