Elastic Wave Isolation and Highly Directive Wave Propagation in Architected Two-Dimensional Metamaterials: An Accelerated Algebraic Approach

Achieving elastic wave isolation and noise reduction is essential for ensuring the structural integrity and operational safety of systems across aerospace, advanced machinery, modern architecture, and railroad infrastructure [1], [2], [3]. In contemporary engineering, however, beyond operating in isolation, these structural components are increasingly integrated with embedded sensor networks, active control architectures, and energy harvesting systems. From this interconnected perspective, excessive vibrations at low frequencies threaten much more than basic mechanical stability. While ambient vibrations can be opportunistically harvested to power autonomous sensor nodes, unmitigated structural noise distorts sensitive diagnostic data and disrupts precise control feedback, ultimately compromising the performance and safety of the broader system. Given these systemic vulnerabilities, mitigating elastic waves at low frequencies while preserving the capacity for energy recovery poses a critical yet persistent challenge. Traditional attenuation strategies consistently struggle in this regime, as the long wavelengths inherent to these frequencies typically require bulky treatments that impose severe weight penalties and spatial constraints [4].

To overcome these limitations, periodic metamaterials have emerged as a promising solution, facilitating lightweight wave manipulation through rationally engineered microstructures while offering the distinct versatility to function independently or integrate synergistically within coupled multi-physics systems. The primary functional hallmark of these structures is the generation of band gaps, specific frequency ranges where the propagation of elastic waves is strictly prohibited [5]. The formation of these attenuation zones is fundamentally governed by two distinct physical mechanisms: Bragg scattering, which suppresses wave propagation via destructive interference across the periodic model, and local resonance, wherein localized oscillators trap and dissipate wave energy [6], [7]. Beyond broadband wave isolation, the versatility of metamaterials is highlighted by their strong spatial selectivity, enabling the precise steering of waves along predefined paths within specific frequency ranges known as pass bands [4]. This capability can be further integrated into multi-physics systems, where the dual ability to simultaneously block and route waves functions as a dynamic mechanical filter. Accordingly, by shielding sensitive hardware from harmful noise while intentionally guiding useful signals, the structure actively ensures accurate data collection and maintains the stability of sensor feedback loops.

Motivated by these operational functionalities and driven by their systemic filtering capabilities, extensive research over the past two decades has been dedicated to optimizing metamaterials, particularly for wave isolation at low frequencies. In pursuit of enhanced performance, researchers have investigated a myriad of sophisticated topological configurations, ranging from foundational re-entrant [8], chiral hybrid [7], and hierarchical frameworks [9] to intricate star shaped hexagonal lattices [10]. Subsequent microstructural refinements have further tailored isolation profiles at low frequencies through the strategic integration of localized features, such as nodal cantilevers [11], curved and semicircular ligaments [5], [12], and hammer inspired resonant unit cells [13]. Moving beyond purely geometric strategies, contemporary research demonstrates that coupling these periodic structures with external active stimuli [14], [15], inertial amplification mechanisms (IAM) [16], [17], and acoustic black hole (ABH) phenomena [18], [19] can significantly expand the boundaries of wave suppression at low frequencies. Further practical implementations of these periodic multi-physics couplings are highly relevant in complex railroad infrastructure, where soil interaction dynamics are effectively captured using the Winkler- Pasternak-foundation model [20]. Within this context, the conventional concrete floating slabs supporting the rail (represented as a uniform beam in Figure 1a) can be replaced with periodic metamaterials. Furthermore, rather than relying on simple uniform beams equipped with basic piezoelectric patches (see Figure 1b), these sophisticated structures can be electromechanically coupled with inductor resistor shunting circuits [21]. This specific configuration operates as a dual purpose mechanism; it functions as a highly tunable electrical shock absorber for robust broadband vibration suppression while simultaneously harnessing the captured wave energy for power generation, rather than merely dissipating it.

Translating these advanced concepts into functional metamaterials necessitates rigorous computational frameworks to accurately capture their wave propagation dynamics. Consequently, a robust analytical landscape has evolved to mathematically elucidate these dispersion phenomena, prominently featuring the plane wave expansion method (PWEM) [22], the dynamic stiffness method (DSM) [23], [24], variational differential quadrature method (VDQM) [25], and the finite element method (FEM) [8]. Among these paradigms, FEM has solidified its position as the preeminent numerical standard. By imposing periodic boundary conditions (PBCs) onto a representative unit cell, FEM elegantly reduces the infinite periodic array into a localized eigenvalue problem. This enables the systematic extraction of comprehensive wave propagation characteristics by sweeping the wavevector along the boundaries of the Irreducible Brillouin Zone (IBZ) within the First Brillouin Zone (FBZ).

Despite its widespread adoption, the practical implementation of FEM is frequently encumbered by severe computational bottlenecks. These limitations primarily stem from the strict necessity for exceedingly dense spatial discretization to guarantee numerical accuracy, particularly when evaluating metamaterials characterized by slender ligaments and localized geometric discontinuities. Such intensive meshing drastically inflates the number of spatial degrees of freedom (DOFs) within the unit cell, which directly magnifies the scale and computational cost of solving the associated eigenvalue problem. The severity of this computational burden is further compounded in complex topologies with multiple periodic interfaces, as achieving mathematical convergence in these structures mandates a highly refined sampling of the wavevector space. Consequently, the sheer volume of these expensive numerical evaluations creates a cumulative delay that cascades through the entire workflow, ultimately rendering extensive parametric sweeps and iterative design optimizations prohibitively slow.

To circumvent these prohibitive computational costs, extensive research has sought to accelerate band structure calculations without sacrificing numerical accuracy. Among the various strategies explored, two primary methodologi-cal paradigms have emerged to address this challenge: the implementation of advanced numerical solvers engineered for rapid convergence [26], [27], and the formulation of reduced-order modeling (ROM) techniques [28]. In recent years, ROM has garnered widespread adoption, emerging as the preferred framework due to its highly favorable balance between computational efficiency and model fidelity. Fundamentally, these dimensionality reduction schemes are delineated by their operational domains, typically classified into physical-space and modal-space methodologies. Physical-space techniques such as static condensation [29] and the improved reduced system (IRS) method [30] compress the global system matrices via a transformation matrix that strategically retains a subset of essential DOFs strictly within the physical coordinate system. Conversely, modal-space schemes achieve drastic matrix compression by projecting the internal DOFs into a truncated modal subspace, while preserving the interface DOFs in the physical domain to seamlessly enforce PBCs and facilitate structural assembly.

This hybrid projection strategy is fundamentally rooted in the Craig-Bampton (CB) method [31], a technique within the broader framework of component mode synthesis [32]. The classical CB procedure is executed by artificially restraining the interface DOFs to extract a set of fixed-interface normal modes representing the internal dynamics. Concurrently, a set of static constraint modes is derived to accurately capture how boundary displacements statically deform the interior domain of the unit cell. To achieve substantial model size reduction, the set of fixed-interface normal modes is partitioned such that only the dominant low-frequency modes are retained, while the higher-frequency residual modes are truncated. This truncation significantly reduces the dimensionality of the reduced system while preserving the characteristics governing the low-frequency response of the structure.

The manner in which these truncated residual modes are mathematically treated subsequently defines the specific reduction framework employed. In the standard Bloch mode synthesis (BMS) formulation [33], the residual modes are completely discarded following truncation, thereby yielding a highly compact reduced system. While this strict elimination strategy greatly enhances computational efficiency, it also implicitly neglects the potential contribution of the discarded higher-frequency modes to the global dynamic behavior. Conversely, the generalized Bloch mode synthesis (GBMS) framework [28] addresses the accuracy limitations associated with strict modal truncation by retaining the quasi-static influence of the discarded modes through a residual compliance matrix. By approximating the stiffness contribution of the truncated high-frequency modes without explicitly including them in the reduced modal basis, GBMS partially restores the dynamic fidelity of the reduced system while maintaining a relatively compact model size. Building upon these foundational formulations, subsequent research has further extended the applicability and robustness of BMS-based reduction strategies. Notable developments include the integration of algebraic condensation techniques into the GBMS workflow [34], the mathematical extension of the BMS formulation to address the inverse $k(\omega)$ eigenvalue problem [35], and the adaptation of these reduction frameworks to complex multi-physics systems. Examples include electromechanical metamaterials incorporating piezoelectric resonant shunt damping [21] and acoustic structures exhibiting non-classical damping characteristics [36]. Furthermore, synergistic formulations coupling GBMS with the FEM [37] have been successfully implemented, enabling a more systematic balance between computational efficiency and accurate wave propagation analysis.

Although model reduction yields a highly efficient eigenvalue problem, robustly enforcing PBCs introduces a subsequent mathematical challenge. Conventional formulations map independent DOFs to dependent DOFs via explicitly derived, topology-dependent transformations, a laborious process requiring complete re-derivation for every unique unit cell geometry. While penalty-based strategies such as the virtual spring method (VSM) [38] attempt to resolve this by decoupling the PBCs to prevent the recalculation of system matrices at each wavenumber, their reliance on empirically tuned penalty parameters frequently induces numerical instability and spectral distortion. To circumvent these limitations, the null space method (NSM) [39] offers an exact, parameter-free, and numerically stable algebraic framework. By defining PBCs through the null space of linear constraints, the NSM guarantees well-conditioned systems and isolates wavenumber dependencies, thereby keeping property matrices constant throughout the dispersion analysis. Consequently, coupling this exact enforcement with modal-space reduction yields a highly efficient computational framework, rendering the extensive parametric optimization of complex metamaterials practically viable.

In comparison with established paradigms like the DSM [24], which provides benchmark-quality results by formulating exact, frequency-dependent element matrices, the developed framework demonstrates exceptional computational efficacy. While the DSM inherently entails the computationally demanding resolution of non-linear transcendental eigenvalue problems, the proposed methodology achieves an accuracy closely rivaling that of the DSM, yet operates at significantly accelerated speeds via a linearized algebraic formulation. Specifically, extracting eigenfrequencies within the DSM framework typically requires iterative root-finding algorithms that are highly susceptible to missing closely spaced modes. This numerical vulnerability is entirely bypassed by the direct eigensolution utilized in the present approach. Furthermore, when evaluated against alternative variational schemes such as the VDQM [25], the architectural superiority of the proposed framework becomes even more pronounced. Although VDQM successfully discretizes models with substantially fewer DOFs than standard FEM, it fundamentally lacks mathematical universality. Due to this limitation, its governing equations demand re-derivation and explicit verification for every distinct element, including straight, curved, or spatially gradient configurations. Notably, this rigid geometric limitation is a drawback also shared with the DSM. Conversely, the developed model inherently preserves the universal adaptability of FEM, seamlessly accommodating all such geometric complexities to any desired degree of accuracy while sustaining robust computational performance.

Despite recent advancements in computational mechanics, engineering a monolithic, lightweight metamaterial capable of effective low-frequency wave isolation across diverse environments remains a formidable challenge [40]. In terms of geometric design, the implementation of topological configurations governed by continuous parametric curves, including arcs [41], sinusoidal profiles [42], and cubic Bezier curves (CBCs) [43], remains noticeably underexplored in the current literature. Furthermore, from a methodological perspective, previous investigations have predominantly centered on band gap generation, thereby providing only marginal analysis of wave directionality and propagation characteristics within the pass bands. However, comprehensively understanding these pass band dynamics is crucial from a systems engineering standpoint, as predictable wave routing can be coupled with embedded sensor arrays to optimize structural health monitoring and data acquisition processes. Despite its system-level importance, accurately mapping these dispersion characteristics necessitates dense numerical sampling across the FBZ rather than the restricted IBZ, which amplifies the computational burden. While ROM frameworks offer a theoretical pathway to mitigate this prohibitive cost, their deployment to accelerate such exhaustive, full-domain dispersion analyses remains virtually unprecedented. This computational bottleneck is further compounded by a specific methodological void: the explicit integration of the NSM for enforcing PBCs within FEM-based ROMs of periodic metamaterials has not yet been documented. Addressing these coupled topological and computational voids therefore presents a critical opportunity to pioneer novel contributions in the predictive design and analysis of multi-physics metamaterial systems.

To bridge these identified gaps, this study introduces foundational contributions across two distinct but highly complementary fronts: computational methodology and structural topology. First, an advanced semi-analytical framework is developed to enable the highly efficient, algebraic analysis of wave propagation. The underlying mathematical model is rigorously formulated utilizing Timoshenko beam finite elements, establishing a generalized linear eigenvalue problem. To ensure computational tractability, dimensionality reduction is subsequently executed via the BMS. Within this reduced modal space, the NSM is systematically integrated to algebraically enforce PBCs, flawlessly extrapolating the unit cell dynamics to an infinite periodic array. By mathematically decoupling the PBCs from the physical domain meshing, this algebraic treatment effectively bypasses the numerical ill-conditioning that frequently emerges with manual transformation (MT) methods. The fidelity and computational superiority of this proposed BMS+NSM framework are then rigorously benchmarked against three established paradigms: a full-order model (FOM) utilizing NSM, an FOM employing a MT, and a ROM coupling BMS with MT. This comprehensive verification encompasses comparative evaluations of band structures, mode shapes, transmittance spectra, and iso-frequency contours, supplemented by a quantitative error analysis evaluating computation time against modal truncation limits.

Second, serving both as a rigorous test bed for the developed framework and as a standalone topological innovation, a novel, Bezier-based curved re-entrant hybrid metamaterial (CRHM) is proposed. Designed to overcome the low-frequency attenuation limitations inherent to conventional periodic structures, the CRHM architecture leverages a fully parametric, monolithic design that maximizes mitigation capabilities while preserving structural simplicity and low global mass. Specifically, the geometric continuity of the Bezier profiles eliminates the impedance mismatches found in standard angular designs, thereby suppressing the generation of mechanical scattering. Crucially, this single-phase configuration intentionally circumvents the manufacturing complexities associated with multi-material integration, ensuring seamless compatibility with scalable fabrication paradigms like advanced additive manufacturing and precision laser cutting. Beyond manufacturability, the core advantage of the proposed CRHM lies in its dynamic role within a larger cuber-physical framework. By establishing reliable wave manipulation directly at the structural level, the material acts as a foundational mechanical filter. This reliable passive attenuation significantly unburdens downstream active control mechanisms and enhances the signal-to-noise ratio for integrated sensory networks. Ultimately, this synergistic interaction between the tailored mechanical structure and the system’s operational control and sensing components elevates the overall resilience of the Multiphysics system. To fully quantify these benefits, the validated BMS + NSM framework is leveraged to execute an exhaustive characterization of the CRHM topology, detailing both band gap formation and directional wave propagation. The study concludes with a detailed parametric investigation, elucidating the precise influence of the Bezier curve parameters and unit cell size on the overarching system-level isolation performance.

The remainder of this article is organized as follows. Section 2 establishes the geometric configuration and mechanical modeling of the proposed CRHM topology. Section 3 details the mathematical derivation of the semi-analytical framework, culminating in a rigorous quantitative assessment of its numerical fidelity and computational efficiency. In Section 4, the low-frequency elastic wave isolation capabilities of the CRHM are analyzed, accompanied by a parametric investigation elucidating the influence of fundamental geometric variables on the resulting band structure. Finally, concluding remarks and key findings are delineated in Section 5.

This section establishes the fundamental geometric and theoretical frameworks underpinning the proposed metamaterial. First, the topological architecture of the individual unit cell is formally defined, emphasizing the integration of parametrically controlled curved ligaments. Subsequently, the principles of Bloch's theorem are introduced, providing the mathematical basis for extrapolating the unit cell dynamics into array propagation behaviors.

To address the challenge of achieving low-frequency wave isolation within stringent lightweight design constraints, this section establishes the architectural foundation of the proposed metamaterial. Because traditional monolithic structures inherently struggle to mitigate low-frequency waves without introducing prohibitive mass, a novel hybrid topology is conceptualized. This architecture synergistically integrates curved elements within a re-entrant honeycomb (RH) framework. However, rather than utilizing traditional simple curved walls, this configuration systematically employs parametrically defined curved ligaments, as illustrated in Figure 2. This precise integration culminates in the proposed CBC and RH hybrid metamaterial (CRHM). This geometric substitution is strictly physically motivated: it is well established that curved microstructural elements induce localized modes that seamlessly govern both local resonance and elastic scattering phenomena. Through their dynamic interaction with propagating elastic waves, these tailored curvatures ultimately facilitate the formation of novel, distinct, and significantly wider band gaps [12], [44].

Translating these localized phenomena into the practical realization of wave isolation necessitates a clear understanding of the system’s qualitative operating mechanism. During dynamic loading or physical testing, the overarching re-entrant framework acts as the primary stabilization method, providing the necessary global stiffness to maintain structural integrity under applied static loads. Concurrently, the engineered curved ligaments function as distributed dynamic absorbers; their tailored flexibility facilitates an elastic rebound action that actively traps and dissipates wave energy through local resonance. Furthermore, the precise positioning mechanism of these constituent ligaments within the periodic array ensures that these localized scattering effects are guided and distributed across the entire domain. Ultimately, this synergistic interplay between the stabilizing rigid frame and the rebounding resonant ligaments allows the metamaterials to reliably filter targeted wave frequencies while preserving its overall mechanical stability.

To mathematically parameterize the geometry of these resonant ligaments, the CRHM architecture employs a continuous design approach utilizing Bezier curves. A generic Bezier curve of degree $n$ is completely defined by a sequence of $n+1$ discrete spatial control points (denoted as ${P}_0, {P}_1, \dots, {P}_n$) and is mathematically constructed via a weighted summation of Bernstein basis polynomials. Serving as the fundamental blending functions for this curve generation, the $i$-th Bernstein polynomial of degree $n$, denoted as $B_i^n(\xi)$, is explicitly expressed as [43]:

where, $\xi \in [ 0, 1]$ represents the normalized parametric coordinate spanning the length of the ligament. Because these polynomials act as continuous weighting functions, they dictate the precise localized influence of each discrete control point on the overarching geometry. By taking the linear combination of these basis functions with their corresponding spatial control points, the generalized vector equation defining the continuous trajectory of the Bezier curve, ${P}(\xi)$, is formulated as [43]:

Within this parametric formulation, the degree $n$ dictates the total number of governing control points, thereby establishing the geometric flexibility of the resulting curve. To achieve a high level of shape controllability while preserving mathematical tractability and computational efficiency, the present study adopts a CBC formulation ($n=3$). This specific configuration provides an optimal balance between topological adaptability and structural regularity, affording the precise parametric control over the ligament profiles. Consequently, by substituting $n=3$ into Eqs. (1) and (2), the explicit representation of the spatial trajectory governing the curved ligaments is derived as:

where the starting and ending control points, ${P}_0$ and ${P}_3$, are constrained to lie on the $x$-axis and are separated by a fixed horizontal distance $w$. The first interior control point, ${P}_1$, is defined in terms of its radial distance $r$ from the curve’s midpoint $M$ and the angle $\beta$ that it forms with the $x$-axis. To preserve geometric symmetry, the second interior control point, ${P}_2$, is selected such that the resulting CBC remains symmetric with respect to the midpoint $M$. The explicit coordinate formulations for all four control points are summarized in Table 1. This geometric definition establishes a symmetric and fully parameterized representation of the curved ligaments, which serves as the foundation for the subsequent mechanical tailoring of the unit cell.

Building upon this CBC formulation, the localized internal geometry and consequently the effective mass and stiffness distributions can be precisely tuned through controlled manipulation of the governing control points. Such inherent tunability constitutes a key advantage of the proposed framework: by adjusting only the geometric parameters $r$, $\beta$, and $w$, a broad family of curve profiles can be generated, enabling systematic optimization of wave isolation performance without requiring modifications to the underlying topological layout. Moreover, the two‑segment structural form of the CBC enables seamless geometric compatibility with the RH components of the hybrid model, allowing the Bezier‑based ligaments to fully redefine the framework of the unit cell and establish a metamaterial architecture specifically tailored for robust low‑frequency wave isolation.

Building on this foundation, the investigation of wave interference within the infinite periodic array requires embedding the response of a single representative CRHM unit cell into the broader context of lattice periodicity. This is achieved by enforcing PBCs, which ensure that the displacement and traction fields on opposite boundaries remain kinematically and dynamically consistent with those of an infinite array. Imposing these constraints transforms the finite cell into a mathematically equivalent segment of the unbounded structure, enabling a rigorous characterization of its wave propagation behavior. Because this formulation inherently relies on the phase‑shift relationships prescribed by Bloch’s theorem, the following subsection establishes the foundations of periodicity concept and its role in governing wave motion in periodic media.

The theoretical foundation for analyzing the isolation performance of the periodic structures is firmly rooted in the framework of Bloch’s theorem. To apply this theorem, the geometry of the infinite periodic array must first be established. Specifically, the position of an arbitrary point $\mathbf{p}$ in any given unit cell $(n_1, n_2)$ is defined relative to a point $\mathbf{o}$ in a reference unit cell (e.g., the cyan cell $(n_1 = 0, n_2 = 0)$ in Figure 2d) via the following relation [8]:

In this expression, $\mathbf{r}_{\mathbf{o}}$ and $\mathbf{r}_{\mathbf{p}}$ designate the respective spatial position vectors, while $\mathbf{a}_1$ and $\mathbf{a}_2$ serve as the primitive translation vectors of the direct lattice. To properly evaluate wave dispersion, this physical-domain representation must be mapped onto the reciprocal (wavenumber) domain. Consequently, a set of reciprocal basis vectors $\mathbf{b}_j$ is introduced, governed by their mathematical relationship with the direct lattice vectors as follows:

where, $\delta_{ij}$ denotes the Kronecker delta and the indices $i, j \in \{1, 2\}$ correspond to the spatial dimensions. Within this established framework, the solution for the wave field $U(\mathbf{r}, t)$ assumes the form of a Bloch wave, expressed as:

Here, $\mathbf{r}$ and $t$ designate the spatial position vector and time, $\mathbf{i }= \sqrt{-1}$ is the imaginary unit, and $\mathbf{k} = [k_1, k_2]^\mathrm{T}$ defines the wave vector comprising the corresponding wavenumbers. Crucially, the amplitude function $\tilde{U}(\mathbf{r}, \mathbf{k})$ inherits the exact spatial periodicity of the underlying unit cell, satisfying:

where, $\mathbf{R} = n_1\mathbf{a}_1 + n_2\mathbf{a}_2$ defines the spatial translation vector of the lattice. Substituting this periodicity condition into the aforementioned Bloch wave formulation yields a spatially periodic expression for the total wave field:

Based on this fundamental relationship, the dynamic response $U_{\mathbf{p}}^{(n_1,n_2)}$ at point $\mathbf{p}$ can be directly derived from the corresponding reference response $U_{\mathbf{o}}^{(0,0)}$ as:

with $U_{\mathbf{p}}^{(n_1,n_2)} = U(\mathbf{r} + \mathbf{R}, t)$ and $U_{\mathbf{o}}^{(0,0)} = U(\mathbf{r}, t)$. Because the expression in Eq. (9) rigorously captures the geometric translation and corresponding phase shift between unit cells [8], it directly satisfies the required PBCs.Building upon this theoretical basis, this study proposes a novel computational framework to evaluate wave dispersion with superior precision and efficiency. The complete formulation of this proposed methodology is detailed in the following section.

This section develops a semi-analytical framework to characterize the wave propagation behavior of periodic structures, specifically applied to the proposed CRHM unit cell. This framework synergistically integrates three distinct computational methods into a unified workflow. First, the governing equations of motion for the discrete unit cell are established utilizing the FEM. Second, a BMS reduction technique is implemented to significantly enhance computational efficiency without compromising numerical precision. Finally, the methodology incorporates the NSM to algebraically enforce the PBCs, thereby extending the localized unit cell dynamics to govern the infinite periodic array.

To accurately capture the dynamic behavior of the proposed CRHM unit cell, particularly the complex deformations along its Bezier-profiled curved walls, the spatial discretization is performed using Timoshenko beam theory. Because this theory accounts for both transverse shear deformation and rotary inertia, it provides a robust representation of the wave-structure interactions within the architected geometry. Employing an energy-based approach to derive the FEM equations of motion, the continuous kinematic field of a representative beam element must first be defined. Specifically, the axial displacement $u(x, t)$, transverse flexural displacement $w(x, t)$, and cross-sectional rotation $\psi(x, t)$ (as depicted in Figure 3) are discretized via appropriate shape functions and expanded as follows:

In the preceding expansions, $\boldsymbol{\alpha}_u$, $\boldsymbol{\alpha}_w$, and $\boldsymbol{\boldsymbol{\alpha}}_\psi$ represent the vectors of admissible spatial shape functions, while $\mathbf{q}_u$, $\mathbf{q}_w$, and $\mathbf{q}_\psi$ denote the corresponding time-dependent nodal coordinate vectors. To rigorously capture the kinematic behavior dictated by the Timoshenko beam theory, specifically accounting for transverse shear deformation and rotary inertia in the curved members, appropriately formulated Hermite shape functions are employed [8]. By assembling the individual displacement components, the complete generalized nodal displacement vector for the element can be defined as $\mathbf{q}_e$. Consequently, to derive the element equations of motion via energy methods, the elemental kinetic energy ($T_e$) and potential strain energy ($U_e$) are formulated as:

and the potential energy as:

In these expressions, the material parameters are characterized by Young's modulus $E$, mass density $\rho$, and shear modulus $G = \frac{E}{2(1+\nu)}$, where $\nu$ represents Poisson's ratio. Geometrically, the cross-sectional area and the second area moment of inertia are defined as A = bt and I = $\frac{\text{b t}^3}{12}$, respectively, with b and t denoting the width and thickness of the beam. Based on these fundamental properties, the governing dynamics at the local level are derived by constructing the element Lagrangian, $\mathcal{L} = T_e - U_e$. Substituting this formulation into the standard Euler-Lagrange equation, $\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{q}}_e}\right) - \frac{\partial \mathcal{L}}{\partial \mathbf{q}_e} = \mathbf{0}$, yields the discretized equations of motion for the individual beam element:

Here, the elemental displacement vector is explicitly structured as $\mathbf{q}_e = [\mathbf{q}_u^\mathrm{T}, \mathbf{q}_w^\mathrm{T}, \mathbf{q}_\psi^\mathrm{T}]^\mathrm{T}$, while $\mathbf{M}_e$ and $\mathbf{K}_e$ denote the local mass and stiffness matrices, respectively. Subsequently, to establish a comprehensive and accurate model of the entire unit cell, these local matrices are mapped into a global framework utilizing standard FEM assembly protocols [8]. Mechanistically, this integration is achieved by executing a congruence transformation on the individual element matrices, followed by a systematic superposition of all elemental contributions. This rigorous mathematical aggregation unifies the discrete components into a cohesive continuous domain, culminating in the overarching global equations of motion for the unit cell as [8]:

where, $\mathbf{q}(t)$ is an $N$-dimensional state vector encompassing all DOFs of the unit cell, while $\mathbf{K}$ and $\mathbf{M} \in \mathbb{R}^{N \times N}$ represent the assembled global stiffness and mass matrices, respectively. By postulating a time-harmonic response of the form $\mathbf{q}(t) = \bar{\mathbf{q}}\mathrm{e}^{\mathrm{i}\omega t}$, the differential equations of motion are mapped into the frequency domain, yielding the standard generalized eigenvalue problem as:

In this expression, $\omega^2$ denotes the eigenvalue (representing the squared angular frequency) and $\bar{\mathbf{q}}$ is the corresponding complex eigenvector dictating the mode shape of the unit cell.

Despite the exactness of this FOM, direct analysis of the CRHM or any highly intricate unit cell topology rapidly becomes computationally prohibitive. This computational bottleneck stems primarily from the massive number of DOFs required to accurately resolve spatial dynamics. Specifically, capturing short acoustic wavelengths without numerical dispersion necessitates a highly refined mesh (typically requiring six to ten elements per wavelength), which inevitably generates vast, densely populated system matrices. Furthermore, direct FOM evaluations often suffer from severe numerical ill-conditioning and round-off errors at very low frequencies. This instability occurs because the stiffness matrix $\mathbf{K}$ overwhelmingly dominates the inertial term $\omega^2\mathbf{M}$ as $\omega \rightarrow 0$. To circumvent these formidable computational constraints while preserving the fidelity of the wave analyses, adopting a highly efficient dimensional reduction technique is imperative. Consequently, this study employs BMS, an advanced projection-based ROM framework specifically tailored for the dynamic analysis of periodic structures [33]. The theoretical formulation and implementation of this ROM strategy are comprehensively detailed in the subsequent section.

Conceived as an advanced adaptation of CMS [31], BMS leverages the CB method [32] to systematically condense the internal dynamics of the unit cell without sacrificing global accuracy. The fundamental operation of this reduction technique involves partitioning the global spatial vector, $\bar{\mathbf{q}}$, into mutually exclusive subsets: the interior DOFs, denoted as $\bar{\mathbf{q}}_I$, and the boundary (or interface) DOFs, denoted as $\bar{\mathbf{q}}_B$. Accordingly, the state vector is rearranged into a block-column format as $\bar{\mathbf{q}} = [\bar{\mathbf{q}}_I^\mathrm{T}, \bar{\mathbf{q}}_B^\mathrm{T}]^\mathrm{T}$. Applying this topological segregation, the generalized eigenvalue problem previously established in Eq. (15) can be recast into the following block-partitioned matrix form:

In this expanded formulation, the block subscripts $II$ and $BB$ denote the sub-matrices strictly associated with the interior and interface domains, respectively, while the off-diagonal sub-matrices $IB$ and $BI$ encapsulate the cross-coupling interactions between these two regions.To isolate and extract the fundamental dynamic characteristics of the substructure's interior, a clamped kinematic condition is temporarily imposed on the interface nodes, effectively setting $\bar{\mathbf{q}}_B = \mathbf{0}$. Under this rigid constraint, the top row partition of Eq. (16) structurally uncouples, yielding a reduced, fixed-interface eigenvalue problem that governs the internal DOFs as follows:

Here, the modal matrix $\mathbf{\Phi}_I \in \mathbb{R}^{N_{II} \times N_I}$ is populated by the column vectors of the fixed-interface eigenvectors, $\bar{\mathbf{q}}_I^k$. The corresponding eigenvalues, $\omega_k^2$, constitute the diagonal entries of the spectral matrix, $\lambda_I$. Within the BMS framework, the modal index $k$ spans from $1$ to $N_I$, where $N_I$ represents the truncated number of retained modes, and $N_{II}$ denotes the total number of full internal DOFs. These $N_I$ retained modes, typically identified as the dominant or master modes, are conventionally selected from the lower-frequency spectrum governed by a predefined frequency cut-off criterion [45] or through more advanced mode selection algorithms [46].

To complement this internal dynamic behavior, the static response of the interior DOFs to arbitrary kinematic excitations at the interface must also be accounted for. This is achieved via static condensation, which defines a set of static constraint modes, $\mathbf{\Psi}_{IB}$, formulated as:

Ultimately, a comprehensive CB reduction matrix is constructed by synthesizing the retained dynamic internal modes (extracted from Eq. (17)) with the static constraint modes (derived in Eq. (18)). Utilizing this projection basis, the full-order displacement vector, $\bar{\mathbf{q}}$, can be efficiently and accurately approximated as:

where, $\mathbf{T}_{\text{CB}}$ signifies the CB transformation matrix, which serves to map the reduced generalized coordinate vector, $\hat{\mathbf{q}}$, back to the full-order physical domain, $\bar{\mathbf{q}}$. This condensed state vector, $\hat{\mathbf{q}}$, constitutes a hybrid coordinate system; it comprises the modal participation factors, $\mathbf{\eta}_I$, which govern the truncated internal dynamics, alongside the explicitly retained physical interface coordinates, $\bar{\mathbf{q}}_B$. The sub-matrix $\mathbf{I}_{BB}$ represents an identity matrix commensurate with the dimensions of the interface DOFs. Substituting the kinematic approximation from Eq. (19) into Eq. (16) and applying a Galerkin projection through pre-multiplication by $\mathbf{T}_{\text{CB}}^\mathrm{T}$ yields the unit cell's reduced eigenvalue problem:

In this condensed formulation, the reduced mass and stiffness matrices are respectively defined via the congruence transformations $\hat{\mathbf{M}} = \mathbf{T}_{\text{CB}}^\mathrm{T} \mathbf{M} \mathbf{T}_{\text{CB}}$ and $\hat{\mathbf{K}} = \mathbf{T}_{\text{CB}}^\mathrm{T} \mathbf{K} \mathbf{T}_{\text{CB}}$. The resulting BMS-reduced system matrices possess a dimension of $(N_I + N_B) \times (N_I + N_B)$, where $N_B$ dictates the number of interface DOFs. Given that $N_I \ll N_{II}$, this dimension is drastically smaller than that of the original FOM, which scales as $(N_{II} + N_B) \times (N_{II} + N_B)$. This profound reduction in the system's DOFs translates to a remarkable enhancement in computational efficiency. Consequently, it unlocks the feasibility of conducting high-resolution wave analyses, extensive parametric studies, and iterative topological optimizations that would otherwise remain computationally prohibitive within a purely FOM framework.

Fundamentally, the reduced eigenvalue problem in Eq. (20) captures the localized elastodynamic response of a solitary unit cell. To seamlessly transition from this isolated domain to an infinite periodic medium, the strict imposition of PBCs governed by Bloch's theorem is imperative. A primary theoretical advantage of the BMS technique lies in its targeted condensation strategy: it systematically compresses the internal DOFs while deliberately leaving the physical interface coordinates ($\bar{\mathbf{q}}_B$) entirely intact. This explicit preservation supplies the essential physical nodes required for the direct enforcement of spatial periodicity. Thus, the BMS approach acts as a powerful precursor, generating a computationally lean yet dynamically authentic model optimized for boundary condition integration.

Consequently, this paper systematically contrasts two distinct strategies for implementing these periodic constraints. The initial method requires manually constructing a kinematic transformation matrix derived from the explicit geometric mapping of the unit cell boundaries. This traditional route, however, is fundamentally hindered by its reliance on invariant structural topologies, its vulnerability to manual derivation inaccuracies, and a heightened risk of numerical instability. Conversely, the second technique employs the NSM [47], providing a sophisticated algebraic formulation to extract a transformation that automatically and exactly fulfills the PBC requirements. This generalized, algorithm-driven methodology ensures computational robustness and completely circumvents manual errors, proving exceptionally beneficial for geometrically complex or arbitrary designs. The subsequent sections delineate the comprehensive mathematical execution of both comparative frameworks.

This subsection details the analytical procedure for manually constructing the kinematic transformation matrix necessary to enforce PBCs. The fundamental objective of this approach is to systematically partition the physical interface DOFs into independent (master) and dependent (slave) subsets, subsequently expressing the dependent coordinates as a direct function of the independent set via periodicity relations [8].

To concretize this methodology, the CRHM unit cell illustrated in Figure 2a is utilized as a representative spatial domain. The enforcement protocol commences by reformulating the reduced generalized coordinate vector, $\hat{\mathbf{q}}$, extracted from Eq. (20), and mapping it according to the periodic constraints established in Eq. (9). This projection is mathematically given by:

where, $\bar{\mathbf{q}}_1$, $\bar{\mathbf{q}}_2$, and $\bar{\mathbf{q}}_3$ designate the partitioned nodal displacement vectors corresponding to the unit cell interfaces. The operator $\mathbf{P}_k$ constitutes the manually assembled, wavevector-dependent transformation matrix. The vector $\tilde{\mathbf{q}}$ encapsulates the minimized set of strictly independent generalized coordinates, constructed alongside identity ($\mathbf{I}$) and null ($\mathbf{0}$) sub-matrices of appropriately matched dimensions. To project the dynamic system onto this reduced, independent basis, the explicit kinematic mapping from Eq. (21) is substituted directly into the reduced eigenvalue problem (Eq. (20)). Following this substitution, a Galerkin projection is executed by pre-multiplying the system by the Hermitian transpose of the transformation matrix, $\mathbf{P}_k^\mathrm{T}$, yielding the final governing equation as:

In this condensed formulation, $\tilde{\mathbf{K}}(k_1, k_2) = \mathbf{P}_k^\mathrm{H} \hat{\mathbf{K}} \mathbf{P}_k$ and $\tilde{\mathbf{M}}(k_1, k_2) = \mathbf{P}_k^\mathrm{H} \hat{\mathbf{M}} \mathbf{P}_k$ represent the final, wavenumber-dependent stiffness and mass matrices, respectively. Through this rigorous projection, the computational dimensions of the global matrices are drastically reduced to $(N_I + N_{B,i}) \times (N_I + N_{B,i})$, where $N_I$ denotes the number of retained dominant internal modes and $N_{B,i}$ signifies the quantity of strictly independent interface DOFs.

While the aforementioned manual transformation successfully enforces periodicity, its fundamental reliance on explicit kinematic mapping renders it highly sensitive to the specific geometric topology of the unit cell, limiting its versatility for complex or evolving architectures. To circumvent these topological restrictions, the subsequent section introduces a generalized, purely algebraic framework by utilizing null space projections that imposes exact PBCs.

To systematically resolve the topological dependencies and implementation constraints inherent in manual kinematic mappings, this study employs the NSM [39] to rigorously extrapolate the finite-domain dynamics of the isolated unit cell to an infinite periodic array. The defining computational advantage of the NSM lies in its ability to completely decouple the enforcement of PBCs from the specific finite element shape functions utilized to formulate the underlying spatial domain.

From an operational standpoint, the overarching objective of this algebraic technique is to re-parameterize the unit cell’s generalized coordinate vector, $\hat{\mathbf{q}}$ (derived in Eq. (20)), projecting it onto a purely linearly independent subspace through an exact, mathematically robust transformation matrix. This exact transformation matrix is constructed directly from the fundamental basis vectors that span the null space of the global boundary constraint matrix. By isolating this null space, the algorithm automatically resolves any complex, linearly dependent constraints, entirely circumventing the need for the manual master-slave partitioning. To contextualize this generalized framework using the BDRHM unit cell architecture, the enforcement protocol initiates by recasting the spatially distributed periodic constraints, as previously explicitly mapped in Eq. (21) into a unified, homogeneous linear algebraic system. This initial constraint reformulation is mathematically given by:

In this formulation, $\mathbf{G}_k$ defines the global, wavevector-parameterized boundary constraint matrix. As visually underscored by the vertical ellipses in Eq. (23), a profound advantage of this matrix architecture is its modularity: any supplementary periodicity constraints (necessitated by intricate internal topologies or multi-node interfaces) can be effortlessly appended as sequential rows. This automated row-stacking fundamentally eradicates the susceptibility to implementation errors that plague the derivation of topology-dependent manual transformation matrices.

The analytical objective now shifts to identifying a generalized coordinate vector, $\hat{\mathbf{q}}$, that concurrently satisfies both these algebraic periodic constraints and the reduced dynamics governed by the BMS eigenvalue problem. To achieve this, the NSM is deployed to extract a definitive set of basis vectors spanning the solution space of the homogeneous system in Eq. (23), mathematically representing the set of all vectors that the operator $\mathbf{G}_k$ maps to zero.

Dictated by the rank-nullity theorem, the null space of $\mathbf{G}_k$, denoted as $\mathcal{N}(\mathbf{G}_k)$, is constructed from a fundamental system of linearly independent solutions: $\mathcal{N}(\mathbf{G}_k) = \text{span}\{\boldsymbol{\chi}_1, \boldsymbol{\chi}_2, \dots, \boldsymbol{\chi}_{(N_I+N_B-h)}\}$. Here, $h \equiv \text{rank}(\mathbf{G}_k)$ quantifies the number of independent constraints, $N_I + N_B$ represents the total DOFs retained in the reduced BMS model, and $\boldsymbol{\chi}_m$ (for $m = 1 \dots N_I+N_B-h$) are the resultant NSM basis vectors. By systematically assembling these basis vectors as the constituent columns of a transformation matrix $\mathbf{Z}$, such that $\mathbf{Z} \in \mathbb{C}^{(N_I+N_B) \times (N_I+N_B-h)} = [\boldsymbol{\chi}_1, \boldsymbol{\chi}_2, \dots, \boldsymbol{\chi}_{N_I+N_B-h}]$, the constrained generalized vector $\hat{\mathbf{q}}$ can be exactly parameterized as:

where, $\boldsymbol{\gamma}$ denotes the unknown column vector of independent modal participation coefficients. By substituting Eq. (24) back into the reduced BMS eigenvalue problem (Eq. (20)) and executing a rigorous Galerkin projection via pre-multiplication by the Hermitian transpose of the basis matrix, $\mathbf{Z}^{\mathrm{T}}$, the unit cell's dynamics are condensed into a final, algebraically exact eigenvalue problem:

Within this governing formulation, $\bar{\mathbf{K}}(k_1, k_2) = \mathbf{Z}^{\mathrm{T}}\hat{\mathbf{K}}\mathbf{Z}$ and $\bar{\mathbf{M}}(k_1, k_2) = \mathbf{Z}^{\mathrm{T}}\hat{\mathbf{M}}\mathbf{Z}$ denote the ultra-compact, wavevector-dependent stiffness and mass matrices of the unit cell, respectively. These operators possess a minimized dimensionality of $(N_I + N_B - h) \times (N_I + N_B - h)$. Because these matrices are synthesized using a transformation basis ($\mathbf{Z}$) that intrinsically satisfies all PBCs, the finite domain of the unit cell is seamlessly and accurately extrapolated to capture the wave propagation behavior of the infinite periodic array.

A critical comparative analysis between this NSM formulation and the traditional manual mapping approach reveals several profound computational advantages: First, the framework is highly automated; the extraction of the null space is executed via standard numerical linear algebra routines (e.g., utilizing Singular Value Decomposition via MATLAB's null function), thereby eliminating manual master-slave partitioning. Second, the spatial phase shifts dictated by the wavevectors are natively embedded within the algebraic structure of the constraint matrix $\mathbf{G}_k$, seamlessly unifying geometric constraints with wave propagation physics. Third, the NSM actively condenses the eigenvalue problem by mathematically eliminating the $h$ dependent constraints, yielding highly efficient matrices of size $(N_I + N_B - h)$, in sharp contrast to the bulkier $(N_I + N_B)$ matrices retained in manual enforcement. Finally, as visually established, advancing to structurally intricate unit cells requires only the vertical expansion of the constraint matrix $\mathbf{G}_k$, without fundamentally altering the algorithmic logic. Consequently, this BMS+NSM semi-analytical paradigm establishes an algebraically exact, computationally lean, and exceptionally robust methodology perfectly suited for the dynamic characterization of complex periodic architectures.

To extract the dynamic wave propagation characteristics from either of the condensed eigenvalue formulations (Eqs. (20) and (25)), non-trivial solutions are isolated by enforcing a zero determinant on the characteristic dynamic stiffness matrix. This fundamental algebraic condition generates the characteristic equation, yielding the eigenvalues as the square of the angular frequency, $\omega^2$. By systematically sweeping the wavevector components $(k_1, k_2)$ across the defined boundaries of the reciprocal space, a distinct spectrum of eigenfrequencies is computed for each discrete coordinate pair. In general, two primary analytical representations are leveraged to interpret these spectral results: the one-dimensional band structure, synthesized by evaluating the wavevector strictly along the high-symmetry contour of the IBZ, and the comprehensive dispersion surfaces, generated by sweeping the wavevector across the entire domain of the FBZ. The boundary-restricted band structure is instrumental for a rigorous characterization of the phononic band gaps, explicitly pinpointing the bounding opening and closing frequencies, quantifying the absolute gap widths, and elucidating the underlying modal mechanisms governing wave isolation. Conversely, the full-domain dispersion surfaces provide critical phenomenological insights into the anisotropic and directional propagation behavior of elastic waves.

This section is dedicated to a rigorous evaluation of the computational efficiency and predictive fidelity of the proposed semi-analytical framework. To facilitate this assessment, an in-house code is developed within the MATLAB environment to simulate the dynamic behavior of the CRHM unit cell, utilizing the material and geometric specifications summarized in Table 2. The numerical implementation encompasses both the high-fidelity FOM and the reduced-order BMS architecture. To systematically extract the phononic band structures and dispersion surfaces, two distinct mathematical schemes are employed to enforce the PBCs across both modeling paradigms: the MT in Subsection 3.3 and the NSM in Subsection 3.4. For notational consistency in the ensuing discussion, models relying on the manual transformation are designated as FOM + MT and BMS+MT, whereas those adopting the null space formulation are labeled FOM+NSM and BMS+NSM. Finally, to quantitatively benchmark the spectral accuracy of the ROM representations against the FOM baseline, a relative percentage frequency error, $e_{\boldsymbol{\omega}}$, is introduced:

In this expression, $\boldsymbol{\omega}_f$ and $\boldsymbol{\omega}_r$ denote the eigenfrequencies computed by the FOM and approximated by the ROM, respectively. To provide a precise, localized comparison of these methodologies, this error metric will be specifically evaluated at a targeted high-symmetry point within the IBZ. Beyond this assessment of accuracy, the computational efficiency must also be quantified; therefore, excluding the preliminary FEM assembly time, the total computational expense associated with the band structure analysis is formulated as:

where $t_f$ and $t_r$ designate the total execution times for the FOM and BMS models, respectively, evaluated over $n_{\mathbf{k}}$ discrete wavevector points. Here, $t^{r}_o$ represents the initial, computational overhead incurred during the reduction phase, whereas $t_{\mathbf{k}}^{f}$ and $t_{\mathbf{k}}^{r}$ denote the processing time required to resolve the eigenvalue problem at a single wavevector coordinate for the FOM and BMS paradigms, respectively.

To ensure rigorous cross-verification of the four configurations (FOM + MT, BMS + MT, FOM + NSM, and BMS + NSM), the validation protocol is anchored by a high-fidelity FEM model developed in COMSOL Multiphysics 6.2, which serves as a reliable baseline with strictly preserved parameters. Structured into four progressive phases, the protocol initiates with a comparative benchmark of the proposed formulations to quantify their relative computational efficiency and spectral accuracy. This foundational stage involves extracting the $1\text{D}$ band structure of the model, alongside a comprehensive band gap analysis of a state-of-the-art curved metamaterial benchmark. To ensure these spectral predictions translate to accurate dynamic behaviors, the second phase focuses on modal reliability, rigorously cross-verifying the extraction fidelity of the unit cell mode shapes against the COMSOL-FEM baseline. Returning to the dispersion analysis and expanding beyond the initial $1\text{D}$ regime, the third phase subsequently maps $2\text{D}$ dispersion surfaces through iso-frequency contours. The validation ultimately scales from the unit cell to a finite array to analyze frequency transmission loss, providing a definitive, physical verification of the model's predictive capabilities.

Initiating the first phase of the validation protocol, the analysis commences with the discretized representation of the CRHM unit cell, as illustrated in Figure 4a. Within this geometric framework, the direct lattice vectors are indicated by hex-colored arrows, while the FEM nodes are systematically partitioned into interior and interface domains, denoted by gray and cyan circular markers, respectively. To execute the $1\text{D}$ band structure extraction, the wavevector $\mathbf{k} = (k_1, k_2)$ is constrained to the perimeter of the IBZ, highlighted by the cyan-shaded region in Figure 4b. This spectral evaluation trajectory sweeps through the high-symmetry nodes in the sequence $\text{O} \rightarrow \text{A} \rightarrow \text{B} \rightarrow \text{C} \rightarrow \text{O}$ (traced by orange arrows), strictly guided by the reciprocal basis vectors (hex-colored arrows). For comprehensive reference, the mathematical formulations governing the direct and reciprocal lattices, alongside the precise coordinates of these IBZ symmetry points, are detailed in Table 3.

Before computing these dispersion properties along the defined IBZ trajectory, a preliminary mesh sensitivity analysis was conducted (omitted here for the sake of brevity) to systematically mitigate the spatial discretization errors inherent to the FEM. Based on this analysis, the FOM is discretized using $402$ Timoshenko beam elements with a uniform characteristic length of $0.3\text{ mm}$. This configuration yields a global system comprising $403$ nodes and $1209$ DOFs, guaranteeing highly converged baseline results. All subsequent numerical simulations of this system are executed on a workstation equipped with a 13th Gen Intel^Ⓡ Core $^{\text {TM }}$ i7-13650 CPU (2.40 GHz) and 16 GB of RAM.

Having established this robust FOM baseline, the next step is to rigorously evaluate the trade-off between computational efficiency and spectral accuracy introduced by the BMS reduction. Therefore, a convergence study of the BMS framework was performed prior to the complete band structure extraction, utilizing the exact FOM+MT solution as the definitive benchmark.

Figure 5 delineates this convergence behavior. As illustrated in Figure 5a, the maximum relative frequency error drops precipitously as the number of retained internal modes increases, stabilizing at a numerical error floor of approximately $10^{-5}$ after the inclusion of roughly $50$ to $60$ modes. Figure 5b re-contextualizes this convergence against computational expenditure, demonstrating that this high-fidelity threshold is achieved at a minuscule relative computing time ($t_r/t_f < 0.01$). Crucially, the error and performance profiles for both the BMS+MT (dashed blue squares) and BMS+NSM (red circles) formulations perfectly overlap, confirming that the residual error is strictly an artifact of the CB modal truncation, independent of the PBC enforcement strategy. While $50$–$60$ modes prove sufficient, a highly conservative basis of $N_I = 120$ internal modal coordinates was deliberately selected for the subsequent validations to guarantee uncompromising precision.

Utilizing this established modal basis, the BMS framework facilitates a profound dimensional reduction. By coupling the $N_I = 120$ internal modes with the $N_B = 9$ physical interface DOFs essential for PBC enforcement, the governing system is mathematically condensed from $1209$ down to a mere $129$ DOFs which is an $89\%$ reduction in global dimensionality. Exploiting this highly efficient reduced system, the low-frequency ($0$–$2\text{ kHz}$) band structure is evaluated by sweeping across $n_{\mathbf{k}} = 150$ discrete wavevectors along the IBZ contour (see Figure 4b). To provide an independent, high-fidelity baseline for cross-verification, this spectral sweep is concurrently executed in COMSOL Multiphysics using the Structural Mechanics module. It is crucial to distinguish that while the proposed numerical frameworks (FOM and BMS) are implemented in MATLAB employing $1\text{D}$ Timoshenko beam elements, the COMSOL validation explicitly leverages full $2\text{D}$ solid modeling. Within the COMSOL environment, Floquet PBCs are applied to all interface nodes, and the unit cell domain is discretized with a super-refined triangular mesh, yielding a solved system comprising $19762$ DOFs. As depicted in Figure 6a, which evaluates these methods within the $2\text{ kHz}$ spectrum, the resulting dispersion branches generated by the proposed $1\text{D}$ frameworks, delineated by solid maroon lines (BMS + NSM), dashed cyan lines (BMS + MT), green triangles (FOM + MT), and solid black lines (FOM+NSM), exhibit perfect coincidence. This exact superposition among all four MATLAB methodologies, alongside their excellent visual agreement with the $2\text{D}$ COMSOL baseline (solid orange lines), strongly confirms the absolute predictive fidelity of the ROMs.

To rigorously substantiate these visual observations, a quantitative assessment of the models' performance at vertex point A within the band structure is detailed in Table 4. This tabulation catalogs selected eigenfrequencies alongside their associated relative errors and computational overhead, benchmarked against both the FOM reference and the commercial COMSOL baseline. In the low-frequency regime, the reduced-order approximations demonstrate pristine agreement with the full-order baseline, yielding a mathematically null error margin (0.000%). While this deviation naturally scales within higher frequency bands (culminating in a modest 2.9% discrepancy by the 64th mode), the overall fidelity remains exceptionally high. Crucially, this internal approximation error is entirely acceptable when contextualized against the high-fidelity 2D COMSOL model, which exhibits a maximum relative deviation of 3.6% from the 1D FOM at the same 64th mode. Furthermore, the BMS + MT and BMS + NSM formulations produce indistinguishable error profiles across the entire evaluated spectrum. Such unwavering stability highlights the structural reliability of the presented methodology, validating the framework as a highly efficient, automated mechanism for enforcing PBCs that mirrors the accuracy of heavy, commercial FEA packages.

Regarding computational efficiency, the Table further demonstrates the overwhelming performance superiority of the BMS architecture. This efficiency is fundamentally rooted in a drastic reduction of systemic DOFs. While the COMSOL benchmark requires a staggering $21{,}762$ DOFs, the framework condenses the system to $1{,}209$ DOFs (FOM), which the BMS further compresses to a mere $129$ DOFs, a size reduction of over two orders of magnitude compared to the commercial FEM. This translates directly to processing speed. Evaluating the complete sequence of $150$ $\mathbf{k}$-points necessitates a total time ($t_f$) of $68.3\text{ s}$ and $69.1\text{ s}$ for the FOM + MT and FOM+NSM configurations, and $155.72\text{ s}$ for the COMSOL baseline. In stark contrast, following a nominal reduction time ($t_o^r$) of just $0.39\text{ s}$, the ROMs demand a mere $0.0031\text{ s}$ and $0.0035\text{ s}$ per wavevector ($t_{\mathbf{k}}^r$). Consequently, the total analysis time ($t_r$) plummets to $4.38\text{ s}$ (BMS+MT) and $4.45\text{ s}$ (BMS + NSM), representing an approximate 35-fold acceleration over COMSOL and a 15-fold speedup over the FOM. Finally, the results verify that applying the NSM to the FOM introduces exactly zero algebraic error ($0.000\%$), a flawless PBC enforcement that is inherently preserved within the reduced systems. Ultimately, for applications requiring even more stringent error bounds, the framework's accuracy can be systematically enhanced by modestly expanding the set of retained internal modes ($N_{\text{I}}$), granting researchers precise control over the trade-off between computational fidelity and processing expense.

To further validate the predictive fidelity and overall robustness of the proposed BMS+NSM reduction scheme, its performance is benchmarked against, state-of-the-art metamaterial designs established in recent literature. Specifically, this comparative analysis evaluates the model’s accuracy by contrasting the band gap characteristics extracted from our computed band structure directly against the reference data depicted in Figure 2c of [12]. This validation employs a re-entrant curved unit cell configured with the specific angular parameters $\theta = -10^\circ$, $\psi_l = 30^\circ$, and $\psi_r = 30^\circ$. For a rigorous comparative analysis, the material and geometric parameters are adopted directly from the reference study: elastic modulus $E = 210$ GPa, mass density $\rho = 25 \times 10^3$ kg/m$^3$, and Poisson’s ratio $\nu = 0.25$. The constituent beams have a length $L = 0.125$ m and a slenderness ratio $\beta = \text{t}/L = 1/15$, with the cross-sectional area defined as A = bt and the second moment of area as $I = \text{bt}^3/12$ (where b and t denote the beam width and thickness, respectively).

Consistent with the reference model, the unit cell comprises three interconnected arms, each discretized using $10$ two-node Timoshenko beam elements. Sharing a central connection node, this spatial discretization yields a total of $31$ nodes. Since each node possesses $3$ DOFs, the resulting FOM comprises a total of $93$ DOFs. The band structure is subsequently computed by sweeping the wavevector across the boundaries of the IBZ. To maintain consistency with established literature, the eigenfrequencies ($\omega$) extracted from the eigenvalue analysis are normalized against the fundamental flexural frequency of an equivalent simply-supported beam, defined analytically as $\omega_0 = (\pi^2/L^2)\sqrt{EI/\rho A}$. Consequently, the dimensionless frequency is expressed as $\Omega = \omega/\omega_0$.

To determine the optimal truncation of the modal basis, a convergence analysis of the internal modes was performed to guarantee spectral fidelity. Although the resulting plots are omitted to prevent graphical redundancy, the convergence behavior closely mirrors earlier observations, confirming that a numerical error floor of approximately $10^{-5}$ is robustly achieved upon retaining merely $N_\mathrm{I}=10$ modes. Proceeding with this optimally condensed basis, Table 5 summarizes the locations and widths of the first four band gaps, providing a direct comparative assessment between the reference data and the predictions derived from both the reduced (BMS + NSM, BMS + MT) and full-order (FOM + NSM, FOM + MT) models.

As evidenced by these quantitative results, the proposed formulations demonstrate exceptionally tight agreement with the published findings. Notably, the limits of the first and second band gaps are captured with absolute precision ($0.000\%$ relative error) across all evaluated frameworks. Even for the higher-frequency third and fourth band gaps, the discrepancies remain strictly negligible, capped at a maximum relative error of merely $0.438\%$ for the full-order evaluations and $0.440\%$ for the highly condensed BMS models. This comprehensive benchmarking provides rigorous validation, confirming the formulation's capacity to accurately resolve complex band gap phenomena.

Crucially, this high degree of spectral fidelity is achieved concurrently with substantial gains in efficiency. By truncating the internal DOFs to $N_\mathrm{I}=10$, the global system size is drastically reduced from $93$ to a mere $19$ DOFs. Consequently, following a negligible reduction phase ($t_o^r \approx 0.012$ s), the BMS models significantly accelerate the total analysis time for $150$ $\mathbf{k}$-points, dropping from roughly $8.6$ s in the FOMs down to $1.93$ s (BMS + NSM) and $2.24$ s (BMS + MT). Beyond these raw performance metrics, the proposed framework inherently accommodates complex, multi-node boundary conditions through an automated and modular implementation of constraint equations. This structural advantage effectively mitigates the susceptibility to user errors traditionally associated with topology-dependent manual transformations. In summary, the synergistic combination of accelerated calculation times, uncompromised spectral fidelity, and algorithmic adaptability positions the BMS+NSM methodology as a highly potent numerical platform. It effectively streamlines the complex process of engineering and fine-tuning metamaterial structures for optimal wave isolation.

Building upon the established spectral accuracy and computational efficiency of the proposed framework, the analysis proceeds to validate the spatial characteristics of the unit cell by examining its mode shapes. Specifically, this verification focuses on the mode shapes at the high-symmetry point $\mathrm{A}_5$, $\mathrm{A}_7$, and $\mathrm{A}_9$ (refer to Figure 6a), which correspond to the wavevector coordinates $(k_1, k_2) = (\pi, -\pi)$ within the IBZ.

While the reduced eigenvectors $\mathbf{\tilde{q}}$ and $\mathbf{\breve{q}}$ (derived in Eqs. (22) and (25)) accurately capture the unit cell's governing dynamics, they reside within a generalized modal subspace and do not directly map to the physical field. Consequently, recovering the physical mode shapes for the BMS-reduced models necessitates a two-stage back-transformation. For the BMS+NSM formulation, the reduced generalized coordinates $\mathbf{\breve{q}}$ are first projected back into the full generalized coordinate space $\mathbf{\hat{q}}$ via the algebraically constructed null space basis $\mathbf{Z}_{\mathbf{k}}$, such that $\mathbf{\hat{q}} = \mathbf{Z}_{\mathbf{k}} \mathbf{\breve{q}}$. The absolute physical displacements $\mathbf{\bar{q}}$ are then fully reconstructed by pre-multiplying this intermediate vector with the CB transformation matrix, yielding $\mathbf{\bar{q}} = \mathbf{T}_{\mathrm{CB}} \mathbf{\hat{q}}$. The physical recovery process for the BMS+MT approach follows an identical two-step sequence, differing only in the initial expansion matrix. Here, the manually derived periodicity basis $\mathbf{P}_{\mathbf{k}}$ is applied to the reduced coordinates $\mathbf{\tilde{q}}$ to obtain $\mathbf{\hat{q}} = \mathbf{P}_{\mathbf{k}} \mathbf{\tilde{q}}$, followed by the same CB transformation $\mathbf{\bar{q}} = \mathbf{T}_{\mathrm{CB}} \mathbf{\hat{q}}$. In contrast, because the FOMs bypass the CB reduction entirely, their physical mode shapes are recovered via a single, direct back-transformation step: $\mathbf{\bar{q}} = \mathbf{P}_{\mathbf{k}} \mathbf{\tilde{q}}$ for the FOM + MT model, and $\mathbf{\bar{q}} = \mathbf{Z}_{\mathbf{k}} \mathbf{\tilde{q}}$ for the FOM + NSM model.

Figure 7 compares the CRHM mode shapes obtained from the developed BMS + NSM framework against BMS + MT, FOM + MT and COMSOL simulations, utilizing the material properties from Table 2 and the BMS truncation parameters outlined in Table 4. To facilitate a consistent visual assessment, the magnitude of the modal displacement fields, calculated as $\sqrt{U^2 + V^2}$ (where $U$ and $V$ denote the nodal deformation components along the Cartesian $x$- and $y$-axes, respectively), is normalized to a unified range of $[ 0, 1]$. Arranged sequentially from left to right, the columns depict the mode shapes obtained via COMSOL, BMS + NSM, BMS + MT, and FOM + MT, which are distinctly rendered using the jet, parula, hot, and sky colormaps, respectively. Similar to the previous subsection, a 2D solid model is used in COMSOL, while 1D Timoshenko beam modeling is adopted across the other frameworks. It should be noted that due to the similarity between FOM + NSM and FOM + MT, as well as space limitations, the FOM + NSM results are omitted. To provide a clear geometric reference, the undeformed configuration of the unit cell is superimposed in gray across all plots. The visual comparison reveals an exceptional degree of spatial correlation, confirming that the automated NSM models reproduce the fields predicted by the MT approach with remarkable fidelity. As elucidating the localized deformation patterns of the unit cell is critical for understanding the physical mechanisms governing the opening and closing of band gaps, these validated mode shapes will serve as the foundation for deeper phenomenological investigations in subsequent sections.

Although previous sections successfully validated the band structure of the CRHM along the boundaries of the IBZ, a comprehensive understanding of its wave isolation capabilities requires exploring propagation behavior across the entire reciprocal space. Relying solely on these 1D boundary profiles can obscure highly anisotropic effects and directional wave phenomena; thus, mapping the full spatial dispersion is a critical step for accurately characterizing the metamaterial’s performance. To address this and provide a final verification, the proposed BMS + NSM framework is employed to construct the complete 3D and 2D dispersion surfaces of the unit cell.

Given the established consistency among the various techniques, we streamline this spatial analysis by restricting the comparison strictly to the BMS + NSM framework and the reference FOM + MT model. {For this evaluation}, the unit cell is defined by the properties detailed in Table 2, with the BMS + NSM formulation utilizing 120 interior and 9 interface DOF. This specific DOF allocation directly mirrors the configuration used in preceding sections, which was rigorously validated by the study presented in Table 4.

To generate these surfaces and extract the corresponding iso-frequency contours, the FBZ is uniformly discretized using a dense $100 \times 100$ grid in the $(k_1, k_2)$ wavevector plane, across which the governing eigenvalue problem is solved. Figure 10 visualizes the results of this analysis for the first and fourth modes, plotted against the normalized wavenumbers $k_1 D / \pi$ and $k_2 D / \pi$. The left column (Figure 10a and Figure 10c)) displays the complete 3D dispersion surfaces, where shaded gray horizontal planes denote specific target frequencies. Geometrically, the intersections of the 3D surfaces with these constant-frequency planes yield the 2D iso-frequency contours projected in the right column (Figure 10b and Figure 10d). Within these contour maps, the proposed BMS+NSM predictions (solid maroon lines) are superimposed onto the FOM+MT reference solutions (dashed cyan lines). The two sets of curves exhibit virtually indistinguishable overlap across all examined frequency slices and dispersion branches. This exceptional agreement verifies the capacity of the BMS + NSM approach to capture complex dispersion topologies with high fidelity throughout the entire 2D reciprocal domain.

Beyond predictive accuracy, the BSM + NSM framework demonstrates a pronounced computational superiority. Specifically, under identical wavevector resolution and contour-extraction protocols, this method evaluates the iso-frequency profiles in merely $32.8$ s, compared to the $240.6$ s required by the full model, representing an approximate $9.2$-fold acceleration in runtime. In practice, such reductions in computational overhead are highly advantageous when these forward simulations must be evaluated repeatedly, such as during topological optimization or inverse parameter identification.

Beyond these computational advantages, accurately mapping these iso-frequency contours is not merely a validation exercise; rather, these profiles serve as the foundational input for characterizing wave propagation behavior within periodic structures, including the CRHM. Consequently, the subsequent sections will leverage these 2D dispersion properties to conduct an in-depth analysis of energy transport metrics, specifically focusing on group velocity, spatial directivity, and their combined influence on the structure's overall attenuation efficacy.

The preceding subsections analyzed wave propagation within an infinite periodic CRHM array using BMS-based frameworks, specifically comparing the MT and NSM algebraic schemes for PBCs. However, it is crucial to recognize that practical engineering applications are inherently finite in their dimensions. Unlike idealized infinite models, real-world structures are subject to boundary effects, spatial constraints, and truncated periodicity that can significantly alter their wave-scattering characteristics. Therefore, modeling spatially bounded configurations becomes indispensable for accurately capturing the true dynamical complexities and operational boundaries of these metamaterial systems. Motivated by the necessity to reconcile idealized theoretical predictions with physical applicability, the present subsection investigates the dynamic response of a truncated $3 \times 8$ CRHM model (depicted in Figure 8). Fundamentally, this finite-domain analysis is driven by a twofold objective: to explicitly corroborate the derived band gap regimes and to systematically quantify the actual wave isolation capacity inherent to a physically realizable construct [48].

To facilitate this assessment, the steady-state harmonic response of the finite structure is extracted utilizing a dedicated FEM framework programmed within the MATLAB environment. This computational procedure hinges on formulating and solving the global dynamic stiffness equation in the frequency domain, which is governed by the following mathematical relation [12]:

$\left(\mathbf{K}-\omega^2 \mathbf{M}\right) \mathbf{U}(\omega)=\mathbf{F}(\omega)$

where, $\textbf{K}$ and $\textbf{M}$ represent the global stiffness and mass matrices, $\textbf{U}(\boldsymbol{\omega})$ denotes the nodal displacement vector at a given angular frequency $\boldsymbol{\omega}$, and $\textbf{F}(\boldsymbol{\omega})$ is the applied harmonic force vector. To efficiently compute $\textbf{U}(\boldsymbol{\omega})$ at each discrete frequency step, MATLAB's backslash operator ($\backslash$) is utilized, ensuring superior numerical stability and computational speed compared to direct matrix inversion. To initiate the elastic waves, a $1$ N harmonic point load (indicated by the magenta arrow in Figure 8) is applied at the input node on the left boundary of the array (yellow point), and the resulting response is measured at the output node on the opposite end (green point). At each increment of the $0\text{--}2000$ Hz frequency sweep, the wave transmission level, $T$, is rigorously quantified on a logarithmic scale according to the following relation:

where $U_r$ and $U_b$ represent the steady-state displacement amplitudes at the excitation and measurement locations, respectively. Evaluating this metric across the designated frequency spectrum yields the transmittance diagram. This spectral representation serves a dual analytical purpose: firstly, it provides a direct quantitative measure of elastic wave isolation, wherein profoundly negative decibel values denote robust wave suppression. Secondly, it delineates contiguous attenuation regions that physically manifest the theoretical band gaps, thereby confirming the reliability and physical realizability of the preceding infinite-lattice dispersion computations.

While the aforementioned transmittance computations effectively isolate the attenuation arising from structural periodicity, accurately predicting the true magnitude of wave suppression in a physical prototype necessitates the inclusion of internal energy losses. Therefore, to replicate realistic physical conditions within the bounds of linear elasticity, the computational framework is augmented to account for inherent material damping. This dissipative behavior is mathematically integrated by prescribing a hysteretic loss factor, $\eta$, which transforms the baseline material properties into complex dynamic moduli, defined as $E_d = E(1 + i\eta)$ and $G_d = G(1 + i\eta)$. Consequently, the resulting global stiffness formulation assumes a complex-valued nature, thereby guaranteeing that the computed transmission profiles properly reflect spatial energy dissipation mechanisms. Following the extraction of the damped global nodal displacement vector, the continuous internal kinematic field of any constituent ligament can be rigorously recovered via local interpolation, leveraging the cubic Hermite shape functions of the beam elements [12].

While the in-house MATLAB formulation provides a comprehensive toolset for extracting both global nodal behaviors and localized internal kinematics, establishing its computational fidelity requires independent validation. Therefore, to provide rigorous cross-verification and benchmark the custom MATLAB framework, a secondary computational model of the identical array geometry is developed. This multi-platform validation strategy serves to eliminate potential algorithmic biases inherent to any single numerical environment. Specifically, the array's geometry is computationally generated within Abaqus 2023 and subsequently exported as a standard \texttt{.step} file into COMSOL Multiphysics 6.2. Within its environment, steady-state frequency-domain simulations are executed utilizing the structural mechanics module equipped with Timoshenko beam elements. Finally, edge probes are strategically positioned at the input and output boundaries of the COMSOL model to precisely extract the kinematic response vectors, $U_r$ and $U_b$.

Figure 6b presents the resulting transmittance diagram for the finite CRHM array, comparatively evaluating the two aforementioned FEM frameworks: the in-house MATLAB code (solid maroon lines) and the commercial COMSOL model (dashed cyan lines). To accurately capture realistic energy dissipation within both models, a uniform structural loss factor of $\eta = 0.01$ is prescribed (as detailed in Table 2), ensuring that the predicted transmission profiles faithfully reflect physical attenuation limits. Furthermore, to maintain consistency with the preceding Timoshenko beam-based formulations and the $3$ mm mesh size verified in Subsections 4.1 to 4.3, the full array model requires a mesh density of $180$ elements per unit cell to guarantee spatial convergence. Consequently, this discretization leads to a global system of $1254$ elements and $3726$ DOFs, wherein all constituent unit cells share uniform geometric parameters and material properties.

A comparative assessment of the resulting spectra reveals exceptional agreement between the two FEM methodologies. The models yield near-identical predictions regarding the spectral locations of the band gaps and the overall magnitude of the damped wave isolation, thereby firmly establishing the accuracy and reliability of the proposed framework. It is worth noting that marginal deviations in the transmissibility profile, most prominently localized around sharp resonant peaks and deep attenuation troughs, are entirely anticipated. These minor discrepancies are directly attributable to inherent variances in the underlying numerical architectures, encompassing minor differences in element formulations, disparate linear solver algorithms, and underlying numerical tolerances. Beyond contributing to these minor spectral variations, these inherent architectural disparities profoundly impact computational efficiency. When evaluating $400$ discrete frequency samples of transmittance $T$ across the $0$ to $2$ kHz bandwidth (utilizing a $5$ Hz resolution), the COMSOL model resolves the system in a mere $217$ s. Conversely, achieving an equivalent level of accuracy demands a higher execution time of $374$ s in MATLAB.

To further substantiate the analyses and validate the proposed model's proficiency in delineating pass band from band gap regimes, Figure 9 presents a comparative visualization of the normalized forced steady-state kinematic fields. These dynamic responses, computed via the in-house MATLAB-FEM and the commercial COMSOL-FEM, are evaluated at $868$ Hz (pass band) and $500$ Hz (band gap), corresponding to the frequency markers in Figure 6b.

To ensure rigorous graphical consistency, the deformed configurations map the real component of the complex steady-state field. Color intensities denote the normalized global displacement magnitude, defined mathematically as $\sqrt{U^2 + V^2}$ within the X-Y plane, and are scaled relative to the peak displacement of each respective simulation. The undeformed periodic array is superimposed in gray to provide a static geometric baseline. Identical normalization and scaling protocols are enforced across both computational platforms to facilitate a direct, unbiased comparative assessment of wave propagation versus attenuation mechanisms.

The reconstructed spatial fields demonstrate remarkable concordance between the MATLAB-FEM predictions (Figure 9a and Figure 9c) and the COMSOL-FEM benchmarks (Figure 9b and Figure 9d)). Specifically, at the $868$ Hz pass band frequency, both frameworks successfully capture the spatially extended, continuous wave propagation indicative of unhindered dynamic transmission, resulting in virtually indistinguishable deformation profiles across the finite array. In sharp contrast, when evaluated at the $500$ Hz band gap frequency, the incident waves are rendered distinctly evanescent. Under these conditions, the structural vibration is drastically suppressed, remaining deeply localized within the immediate vicinity of the excitation source. Both methodologies capture this rapid spatial decay with striking fidelity, yielding spatially coincident zones of localized strain energy and equivalent attenuation profiles. By seamlessly bridging these localized deformation mechanics with the global transmission behavior, these finite-array steady-state visualizations ultimately serve as a robust, direct physical validation of the infinite-domain BSM+NSM unit cell predictions. It should be noted that the provided results is strictly confined to small-amplitude wave motions and infinitesimal deformations. Consequently, amplitude-dependent dissipation and complex nonlinear effects such as geometric and material nonlinearities, large rotations, and contact or friction phenomena fall outside the present scope and are not considered in this study.

Having rigorously established the validity of the proposed computational framework, the focus now shifts toward quantifying the wave isolation characteristics of the CRHM. To this end, the dispersion band structure is derived utilizing the combined BMS+NSM approach. The emergent band gaps are systematically characterized by identifying their exact bounding (opening and closing) frequencies and elucidating the underlying physical mechanisms. Particular emphasis is placed on isolating the lowest and widest band gaps, as these spectral regions fundamentally dictate the low-frequency and broadband attenuation capabilities of the CRHM. This foundational dispersion analysis is subsequently enriched by a comprehensive evaluation of phase velocity, group velocity, and spatial wave directivity. Finally, the investigation bridges these theoretical infinite-domain insights with practical applications by evaluating the steady-state transmission responses of finite CRHM arrays.

Leveraging the validated BMS+NSM theoretical framework, the attenuation performance of the CRHM is systematically evaluated by computing its fundamental band structure. To visually benchmark these spectral characteristics, Figure 11 contrasts the baseline RH with the proposed CRHM, utilizing the geometric and material constants specified in Table 2. The complete 3D dispersion surfaces over the FBZ are presented for the RH in Figure 11a and the CRHM in Figure 11d, alongside their corresponding 2D band structures evaluated along the IBZ path in Figure 11b and Figure 11e, respectively.

Analyzing these 3D dispersion surfaces reveals the physical origin of the CRHM's enhanced filtering capabilities. While the conventional RH unit cell exhibits numerous overlapping and intersecting modal surfaces, the CRHM topology induces distinct frequency separations across the entire reciprocal domain. In the projected 2D band structure (Figure 11e), these separations manifest as eight complete and absolute band gaps, indicated by the shaded lime green regions. Because the practical efficacy of a metamaterial is heavily dictated by the spectral locations and spans of its lowest and widest band gaps [4], these regions warrant rigorous examination. Notably, the CRHM opens its lowest frequency band gap between $163$~Hz and $179$~Hz, alongside its widest band gap spanning from $332$~Hz to $706$~Hz. In stark contrast, the baseline RH model fails to open any absolute band gaps within the entire $0$ to $2$~kHz operational frequency regime.

To further quantify this broadband attenuation capacity, the cumulative band gap coverage, defined as the ratio of the forbidden frequency bandwidth to the entire analyzed spectral range, serves as a critical performance metric [4]. For the CRHM, the absolute sum of all eight band gap widths below 2~kHz amounts to 1377~Hz, yielding an exceptionally high spectral coverage of 68.85% (detailed in Table 6). This superior performance is physically driven by the strategic integration of CBCs. By replacing traditional straight walls and sharp vertices, these CBCs fundamentally alter the unit cell's internal mass-stiffness distribution, thereby inducing rich localized scattering effects that block wave propagation. Ultimately, these dispersion characteristics conclusively demonstrate the proposed design's exceptional proficiency in low-frequency wave isolation.

Beyond their isolated theoretical deployment as wave attenuators, it is crucial to evaluate the integration efficiency that translates these structural traits into system-level performance. Traditional metamaterials, such as conventional RH, frequently suffer from unpredictable high-frequency wave scattering and stress concentrations at their geometric junctions, resulting in poor performance in the sub-2~kHz regime (as demonstrated in Figure 11a and Figure 11b). Compounding this physical limitation, when deployed in active engineering environments, these scattering phenomena manifest as structural noise that corrupts sensory feedback. To compensate for this, control systems are forced to rely on intensive software-side signal filtering, which inevitably introduces latency and phase lags that threaten the stability of the entire control loop.

The proposed CRHM directly resolves these systemic issues; its continuous parametric curves allow it to function as an intrinsic mechanical filter. Unlike conventional structures that permit low-frequency noise to propagate and cause control spillover, the CRHM physically suppresses these disruptive waves within its optimized band gaps. This effectively pre-filters the incoming energy, maintaining a high signal-to-noise ratio prior to any digital processing. By inherently isolating the system from such multi-frequency disturbances, the CRHM significantly reduces the computational burden on active controllers and guarantees robust stability.

Finally, while the CRHM clearly outperforms traditional static structures, it also offers distinct advantages over contemporary advanced isolation strategies. Modern approaches frequently rely on highly complex mechanisms to achieve broadband isolation, such as the spatiotemporal modulation [49], or the intricate geometries of multi-body ABH systems [50]. Although theoretically powerful, these methods introduce significant operational complexities and manufacturing hurdles that complicate real-world deployment. In contrast, the CRHM operates as a streamlined, single-phase mechanical system that completely bypasses the need for active modulation or intricate multi-component assembly. This elegant simplicity facilitates seamless integration into diverse engineering applications, simultaneously ensuring robust operational safety and critical structural lightness without the need for auxiliary subsystems.

Investigating the modal kinematics of the unit cell provides critical insight into the physical principles governing band gap formation [40]. Specifically, the initiation and termination of a band gap are fundamentally linked to the distinct vibrational morphologies exhibited at its bounding frequencies. To elucidate these transition mechanisms, this subsection focuses on the lowest and widest band gaps. Figure 12 illustrates the normalized profiles extracted at the precise frequency limits of these zones for the CRHM, denoted by $\text{L}_{\text{o}}$, $\text{L}_{\text{c}}$, $\text{W}_{\text{o}}$, and $\text{W}_{\text{c}}$ in Figure 11c. For comparative reference, the undeformed unit cell configuration is superimposed as a solid outline.

As observed in Figure 12a, the onset of lowest band gap is driven by an offset resonance within the unit cell. This behavior effectively traps and attenuates the wave energy, thereby inhibiting propagation across the periodic array. Conversely, at the upper frequency bound (Figure 12b), the kinematic energy becomes delocalized and distributed throughout the entire unit cell. This globalized deformation precludes the localized resonance necessary for energy confinement, thus marking the closure of the lowest band gap.

Furthermore, Figure 12c demonstrates that at the lower limit of the widest band gap, wave energy is primarily confined within the CBC-based ligaments. This localization manifests as a flexural resonance mode, effectively trapping the energy and instigating the opening of the broad band gap. Finally, at the upper bounding frequency of this gap (Figure 12d), the response shifts dramatically: the mode shape becomes governed by the motion of the straight vertical cell walls, leaving the CBCs nearly quiescent. This decoupling indicates that the resonators cease to trap energy, thereby restoring wave transmissibility and resulting in the closure of the widest band gap.

Building upon the preliminary validation of the BMS+NSM unit cell framework in Subsection 4.4, this section assesses the wave isolation capabilities of the finite CRHM array. As depicted in Figure 8, the evaluated domain consists of a spatially uniform arrangement of unit cells, constructed using the material and geometric parameters detailed in Table 2. Furthermore, to accurately capture physical energy dissipation during wave propagation, a structural loss factor of 1% ($\eta$ = 0.01) is prescribed throughout the model.

With this numerical framework established, Figure 2b and Figure 2d provide a direct comparison to measure the improved attenuation of the CRHM topology against the RH. As shown in Figure 11b, the baseline RH shows weak wave isolation; its transmittance ($T$) mostly fluctuates near $0$~dB across the $0$--$2000$~Hz spectrum, with only a moderate drop approaching $-80$~dB around $300$~Hz. In contrast, the response of the CRHM (Figure 11d) shows multiple deep attenuation troughs that exactly match the shaded band gaps predicted by the dispersion analysis. Notably, the widest band gap (spanning approximately $330$–$700$~Hz) causes a large transmission loss, with $T$ falling below $-300$~dB. Other higher-frequency band gaps further suppress wave propagation, reaching transmittance values near $-220$~dB between $750$–$1000$~Hz, and approaching $-200$~dB in the $1500$–$1850$~Hz range. This direct comparative assessment conclusively validates the predictive fidelity of the unit cell framework while highlighting the vastly superior isolation capabilities engineered into the CRHM topology.

Having demonstrated the exceptional wave isolation capacities of the CRHM configuration, the analysis now transitions to evaluating its propagation characteristics within the pass band regimes. This shift in focus is essential, as understanding wave transport in these transmissive regions is a critical prerequisite for advanced applications involving spatial energy routing and waveguiding [12]. To systematically characterize how the unit cell geometry dictates these spatial wave kinematics, the direction-dependent dynamic response must be mapped across the 2D reciprocal space. Fulfilling this objective, Figure 13 presents the iso-frequency contours derived from the complete dispersion surfaces, evaluated at discrete frequencies marking the onset and closure of the lowest and widest band gaps. By analyzing the topological features of these spectral contours, the spatial trajectory of propagating waves can be directly inferred. Specifically, whereas perfectly circular profiles reflect standard, isotropic wave scattering, pronounced non-circular topologies reveal intense dynamic anisotropy within the lattice. Ultimately, it is this sharp directional dependence that provides the underlying physical mechanism for structural wave steering and precise directional filtering [12].

The spectral contour maps exhibit stark morphological variations that dictate the CRHM’s spatial kinematics. At the lower boundary of the lowest band gap (Figure 13a), the dispersion surfaces yield elongated, wavy contours aligned predominantly along the diagonal, signifying a strong directional preference in wave scattering. In stark contrast, the upper boundary of this band gap (Figure 13b) is characterized by densely packed, open diagonal lines. Notably, because the direction of energy flow (group velocity) is mathematically orthogonal to these iso-frequency lines, such flat contour topologies dictate that wave energy is highly channeled into non-diffractive beams, a classic hallmark of wave self-collimation [8]. Moving to the lower boundary of the widest band gap (Figure 13c), the response shifts into complex, lobed contour topologies, implying severe spatial dispersion and multi-directional energy routing. Conversely, the upper boundary of this wide gap (Figure 13d) exhibits distinctly localized, rhomboidal contours. The relatively flat facets of these diamond-like geometries suggest specific axes of zero curvature. This acts as a strong indicator of partial wave collimation [12], ultimately confirming the CRHM’s robust capacity for advanced, frequency-selective wave steering within its pass band regimes.

While the morphological analysis of the iso-frequency contours provides profound qualitative insight into the CRHM’s spatial kinematics, a rigorous quantitative assessment of the CRHM’s directional anisotropy is imperative. Hence, this subsection systematically investigates the in-plane wave propagation mechanics within the transmissive pass bands. To translate the observed geometric topologies into actionable physical properties, three critical kinematic metrics are evaluated: phase velocity, group velocity, and energy directivity. Each of these fundamental parameters is mathematically extracted directly from the previously established dispersion surfaces.

The phase velocity ($C_p$) dictates the absolute speed at which a monochromatic wave's phase front propagates along a prescribed spatial trajectory. It is analytically expressed as [40]:

where, $\omega$ signifies the angular eigenfrequency, $|\mathbf{k}| = \sqrt{k_1^2 + k_2^2}$ is the magnitude of the wave vector, and $\hat{\mathbf{k}} = \mathbf{k}/|\mathbf{k}|$ defines the unit directional vector. The angular dependency of the phase velocity is optimally visualized via polar representations, wherein the radial magnitude at any azimuthal angle corresponds to the phase celerity in that specific direction. Within these plots, perfectly circular profiles denote isotropic spatial propagation, whereas distorted or lobed geometries are definitive physical manifestations of dynamic anisotropy, highlighting preferential phase propagation pathways.

While phase velocity characterizes the wave fronts, the group velocity ($C_g$) is fundamentally more critical for structural applications, as it governs the vectoral speed and direction of energy transport (wave packets) through the metamaterial array. Mathematically defined by the gradient of the dispersion surface in the reciprocal space, it is evaluated as \citep{40}:

Crucially, these group velocity vectors are strictly normal to the iso-frequency contours. This orthogonality directly translates the geometric topologies discussed previously into physical energy routing: a dense clustering of $C_g$ vectors aligned along specific spatial axes signifies intense energy focusing, commonly referred to as wave beaming or collimation [8]. Conversely, angular sectors devoid of these vectors constitute propagation blind zones, wherein wave transmission is severely attenuated.

To accurately map this energy transport from the reciprocal wavevector domain into the physical spatial domain, directivity plots are constructed to visualize the permissible propagation angles at discrete frequencies. The physical energy propagation angle, $\varphi$, is calculated from the localized slope of the underlying iso-frequency contour as [48]:

where, $L_1$ and $L_2$ represent the physical dimensions of the unit cell, and $df(k_1)/dk_1$ dictates the tangential slope of the contour curve at the evaluated coordinate. Plotting these extracted angles in a polar format yields a highly intuitive, quantitative map of the metamaterials’s angular transmission spectrum. Congruent with the group velocity distribution, angular voids in these polar plots identify wave-inhibited blind zones. Ultimately, this rigorous directivity evaluation is an indispensable tool for the targeted design of metamaterial devices, such as structural waveguides and acoustic cloaks which intrinsically demand high-fidelity, direction-selective wave manipulation.

Figure 14, Figure 15, Figure 16, and Figure 17 systematically map the phase velocity (first row), group velocity (second row), and energy directivity (third row) at three distinct frequencies (columns) to analyze the wave kinematics at the edges of the CRHM’s lowest and widest band gaps.

Figure 14 details the specific dynamic transition at the lower edge of the lowest band gap. Analyzing the columns from left to right reveals a severe escalation in dynamic anisotropy as the operational frequency approaches the band gap threshold. Initially, the phase and group velocity contours (Figure 14a and Figure 14d) present as continuous, slightly elongated ovals, indicating wave transmission with a mild vertical bias along the $90^\circ$ and $270^\circ$ axes. However, as the frequency shifts (middle and right columns), profound spatial dispersion takes effect. The phase velocity contours (Figure 14b and Figure 14c) develop deep horizontal concavities at $0^\circ$ and $180^\circ$, while the group velocity profiles (Figure 14e and Figure 14f) constrict sharply, forming complex, highly localized transmission lobes.

This dramatic shift in energy routing is most intuitively visualized in the directivity analysis (third row). The omnidirectional, continuous propagation loop observed in Figure 14g rapidly decomposes into discrete directivity rays (Figure 14h and Figure 14i) as the band gap approaches. The dense clustering of these rays along specific spatial axes signifies intense wave beaming. Conversely, the widening angular gaps between these ray clusters definitively mark the onset of propagation blind zones angular sectors where physical energy transmission is strictly forbidden as the model transitions into the band gap.

Conversely, Figure 15 captures the dynamic transition occurring at the upper edge of the lowest band gap, revealing a drastically different wave propagation regime dominated by extreme spatial filtering. Unlike the lower edge, which transitions from continuous to fractured propagation, the upper edge emerges from the band gap in a state of severe wave collimation. The initial phase velocity profile (Figure 15a) consists of completely disconnected polar arcs, mathematically proving that horizontal wave propagation is forbidden. As the frequency evolves, this profile closes into a tightly pinched, figure-eight geometry (Figure 15b and Figure 15c).

The physical energy transport, dictated by the group velocity (Figure 15d, Figure 15e, and Figure 15f), is exclusively localized into sharp, vertically oriented multi-lobed structures, demonstrating that the structure severely penalizes any off-axis energy routing. This near-perfect wave collimation is unambiguously visualized in the directivity arrays (Figure 15g and Figure 15i). Transmission is initially confined to extremely narrow, sharply focused bidirectional cones along the $90^\circ$ and $270^\circ$ axes, surrounded by massive, contiguous blind zones that cover the vast majority of the angular spectrum. Although these transmission wedges slightly broaden as the frequency steps away from the band gap threshold, the CRHM effectively acts as a highly rigid, 1D spatial waveguide within this frequency envelope.

Subsequently, Figure 16 illustrates the dynamic evolution of wave kinematics approaching the lower edge of the widest band gap. This regime presents a distinct anisotropic topology that differs fundamentally from the extreme collimation observed at the lowest band gap’s upper edge. Analyzing the phase velocity (first row), the continuous, eccentric ellipses found deeper in the passband (Figure 16b and Figure 16c) experience a severe topological breakdown nearest the band gap threshold (Figure 16a), rupturing into disconnected, strongly pinched segments. This breakdown in phase propagation is mirrored by the group velocity (Figure 16d, Figure 16e, and Figure 16f). The relatively smooth energy flow profiles distort rapidly into complex, vertically biased petal-like lobes, signaling a highly localized and erratic routing of physical energy.

The definitive impact of this spatial dispersion is captured in the wave directivity maps (third row). While frequencies deeper in the passband Figure 16h and Figure 16i) exhibit dense, nearly omnidirectional transmission, the state immediately bordering the band gap (Figure 16g) fragments dramatically. Instead of generating the razor-thin beams, the directivity here coalesces into broad, diagonal transmission sectors. Crucially, this fragmentation introduces abrupt angular gaps, distinct blind zones that dominate along the primary horizontal ($0^\circ$ and $180^\circ$) and vertical ($90^\circ$ and $270^\circ$) axes. This confirms that the lower edge of the widest band gap imposes strong, but not maximally one-directional, wave collimation, operating as a broad, multi-sector spatial filter immediately prior to total wave attenuation.

Finally, Figure 17 delineates the wave kinematics emerging at the upper edge of the widest band gap, revealing a transition into a state of severe, four-fold orthogonal collimation. The phase velocity maps (first row) capture a profound topological breakdown; the continuous, deeply concave clover-like contour (Figure 17a) violently fractures into four disconnected, cardinal polar arcs (Figure 17b and Figure 17c). This indicates that the array rapidly loses the capacity to support diagonal phase propagation near the band gap threshold.

This behavior governs the physical energy routing, as demonstrated by the group velocity profiles (second row). The continuous energy boundaries seen further in the pass band (Figure 17e and Figure 17f) condense into a sharp, sharply constricted four-lobed cross (Figure 17d). This highly localized topology confirms that energy transport is strictly funneled along the fundamental structural vectors. The practical consequence of this dispersion is visualized in the wave directivity (third row). The relatively omnidirectional transmission fields (Figure 17h and Figure 17i) collapse into an intense, cross-shaped ray clustering (Figure 17g). By highly concentrating wave transmission along the $0^\circ, 90^\circ, 180^\circ,$ and $270^\circ$ axes, the structure simultaneously generates massive diagonal blind zones. Consequently, the CRHM effectively operates as a rigorous orthogonal wave-splitter at the upper boundary of the widest band gap, entirely forbidding diagonal energy propagation. Ultimately, this comprehensive analysis demonstrates that the CRHM exhibits highly tunable, frequency-dependent spatial dispersion at its band gap boundaries, passively toggling between extreme unidirectional beaming, broad multi-sector channeling, and strict orthogonal wave-splitting. This engineered anisotropy confirms the metamaterial’s exceptional capability to function as a versatile spatial filter, precisely steering wave energy along targeted vectors while imposing absolute directional blind zones.

The preceding analyses establish that integrating CBC elements into the baseline RH framework to yield the CRHM architecture fundamentally amplifies the structure's overall wave isolation capabilities. This profound enhancement in dynamic performance is explicitly evidenced by the generation of novel, low-frequency band gaps and the pronounced, frequency-dependent spatial collimation detailed in the prior section. Motivated by these compelling dynamic phenomena, this subsection systematically explores the geometric parameter space to elucidate how specific design variables dictate the formation and evolution of low-frequency band gaps. By mapping these structural sensitivities, the subsequent analysis aims to isolate optimal configurations that maximize subwavelength wave attenuation. To achieve this objective, a comprehensive parametric study is conducted to quantify the influence of both the intrinsic CBC geometric parameters and the overarching unit cell size on the emergent band gap properties.

A defining advantage of the CRHM's Bezier-based architecture is its inherent parametric control, which affords exceptional design versatility. Consequently, its wave attenuation capabilities can be precisely calibrated through nuanced geometric tuning, circumventing the need for extensive structural redesigns. To systematically exploit this morphological adaptability, a comprehensive parametric study is conducted to investigate the sensitivity of the band gap properties to primary Bezier curve variables, namely the control point distance ($r$) and the inclination angle ($\beta$). Accordingly, the band diagrams of the CRHM are derived using the proposed BMS+NSM framework (established in Eq. (25)), evaluating a discrete parameter space defined by $\beta \in \{30^\circ, 45^\circ,60^\circ\}$ and $r \in \{20, 30, 40, 50\}$ mm. The findings of this parametric investigation are comprehensively presented in Figure 18, charting the variation of band gap characteristics with respect to $r$, wherein each column discretely isolates the spectral response for a distinct angle ($\beta$). The upper panels (Figure 18a, Figure 18b, and Figure 18c) delineate the evolutionary trajectories of individual band gap ranges up to $2000$~Hz, while the corresponding lower panels (Figure 18d, Figure 18e, and Figure 18f) quantify the total band gap coverage.

A visual inspection of the dispersion diagrams reveals two prominent phenomenological trends driven by the geometric parameter $r$. First, as $r$ increases from $20$ mm to $50$ mm, the lowest band gap undergoes a conspicuous downshift, progressively lowering toward the lower-frequency spectrum. Second, larger values of $r$ induce a spectral fragmentation effect, manifesting as a higher multiplicity of discrete, often narrower, band gap patches distributed across the frequency domain. This demonstrates that augmenting the $r$ parameter effectively enhances the low-frequency attenuation capabilities of the CRHM.

Examining the quantitative coverage charts elucidates the substantial impact of both $r$ and $\beta$ on the overall attenuation bandwidth. For a relatively shallow angle ($\beta = 30^\circ$), the total band gap coverage exhibits a monotonic enhancement with respect to $r$, nearly doubling from $33.8\%$ at $r=20$ mm to $65.4\%$ at $r=50$ mm. While intermediate localized fluctuations are observed for steeper angles such as the minor coverage reductions at $r=40$ mm for both $\beta=45^\circ$ and $\beta=60^\circ$, the maximum for any given angle is consistently achieved at the upper geometric limit of $r=50$ mm.

Furthermore, a comparative analysis across the columns underscores the critical role of the angle $\beta$. Transitioning from $\beta = 30^\circ$ to $\beta = 60^\circ$ generally promotes a broader total attenuation bandwidth. This enhancement is particularly pronounced at lower $r$ values. For instance, at $r=20$ mm, the coverage expands significantly from $33.8\%$ ($\beta=30^\circ$) to $53.3\%$ ($\beta=60^\circ$).

However, a critical examination of these dispersion characteristics reveals a fundamental design trade-off between total attenuation capacity and spectral continuity. While leveraging higher geometrical parameters, specifically a steep angle ($\beta=60^\circ$) combined with a maximum radius ($r=50$~mm) successfully maximizes the band gap coverage (peaking at $69.0\%$) and effectively drives the onset frequency into the ultra-low regime, this comes at the cost of severe spectral fragmentation. As observed in Figure 18c, these highly-tuned configurations fracture the band gap into multiple narrow, discrete domains separated by intermediary pass bands. Conversely, lower-to-intermediate parameter pairings (e.g., $\beta=30^\circ$ at $r=30$ mm) yield a lower percentage but offer a highly contiguous, uninterrupted broadband attenuating region. Consequently, the selection of CRHM parameters must be dictated by the specific application: targeting lower $\beta$ values when a single, robust, continuous band gap is required, or maximizing $\beta$ and $r$ when the objective is to optimize the total cumulative percentage of band gaps across the spectrum.

While the CBC profile governs the nuanced wave isolation features of the CRHM, the global spatial footprint, defined by the unit cell dimension $D$, dictates its overarching operational frequency regime. Therefore, scaling the unit cell fundamentally alters the spatial periodicity of the lattice, thereby systematically shifting the Bragg scattering boundaries and local resonance thresholds. To quantify this size-dependency and evaluate the spectral scalability of the proposed CRHM, a parametric sweep is performed over the unit cell size, considering $D \in \{42, 44, 46, 48\}$ mm, while all other parameters remain fixed to the values listed in Table 2.

The resulting dispersion characteristics are computed utilizing the established BMS+NSM framework (Eq. (25)), with the findings comprehensively presented in Figure 19. Specifically, Figure 19a delineates the evolutionary trajectories and bandwidths of the individual complete band gap intervals up to $2000$~Hz. As anticipated by fundamental scaling laws, an incremental increase in $D$ induces a noticeable downshift in these distributions, shifting the band gap zones toward lower frequency spectra.

Alongside these frequency shifts, Figure 19b quantifies the total band gap coverage as a percentage of the evaluated $2000$ Hz domain. Despite the low-frequency migration of the individual band limits, the CRHM sustains a highly robust wave attenuation profile across all evaluated sizes. In fact, the overall band gap coverage remains consistently high (exceeding $63\%$), culminating in a peak coverage of $69.3\%$ for the $D = 48$ mm configuration. Ultimately, these results confirm that the proposed CRHM not only affords predictable frequency tunability through simple geometric scaling but also intrinsically preserves its broad wave-suppression efficiency regardless of spatial constraints.

The computational burden of evaluating elastic wave isolation phenomena in advanced architected materials remains a significant bottleneck in the rapid discovery and optimization of high-performance metamaterials. To directly address this challenge, the primary contribution of this research is the development of a novel, ultra-fast, and reliable semi-analytical framework tailored for the dynamic analysis of arbitrary metamaterials. Moving beyond the topological constraints of classical finite element methodologies, this work integrates BMS with the null space method to enforce PBCs algebraically. By synergizing this formulation with a CB component mode synthesis, internal DOF are systematically condensed, granting precise, user-defined control over truncation accuracy. This establishes a highly automated, rapidly convergent workflow capable of coupling complex metamaterial topologies at a fraction of the computational cost associated with FOMs. The fidelity and versatility of this generalized framework have been unequivocally validated through rigorous cross-verification. Predictions of dispersion band structures, unit cell mode shapes, and finite-array transmittance spectra demonstrate excellent agreement across diverse theoretical schemes (e.g., FOM+MT, FOM+NSM, and BMS+MT), alongside independent commercial COMSOL and in-house MATLAB FEM solvers. Furthermore, the developed formulation empowers deep spatial and spectral characterization over the entire FBZ. The extraction of phase velocities and the calculation of group velocities derived from 3D dispersion surface gradients (corroborated by spatial directivity plots) demonstrate the framework’s robust capacity to analyze and predict complex wave routing and spatial manipulation.

Serving as a rigorous test bed to demonstrate the framework’s predictive power and constituting the second major contribution of this work, a novel single-phase curved hybrid metamaterial (CRHM) is proposed. Designed to address the persistent engineering challenge of low-frequency, broadband wave isolation, the CRHM leverages Bezier-based parametric walls to precisely engineer its internal mass-stiffness distribution, thereby activating intense local resonances and localized scattering phenomena. As validated by comprehensive benchmarking (Table 7), the CRHM definitively surpasses existing curved topologies, establishing new performance frontiers by simultaneously achieving the lowest onset frequency, the widest attenuation bandwidth, and maximized spectral coverage. Beyond its utility as an independent attenuator, this architecture exhibits profound advantages when integrated into complex multi-physics systems. By employing parametric geometries, the CRHM inherently avoids the unpredictable high frequency structural noise that affects conventional designs, functioning instead as an intrinsic mechanical pre-filter. In active systems, physically suppressing these disturbances ensures an optimal signal to noise ratio before any digital processing even begins. As a result, this inherent isolation significantly reduces the computational load on active controllers, ensuring highly robust system stability. Crucially, these exceptional dynamic capabilities are realized within a geometrically streamlined, lightweight, and fabrication aware topology. By minimizing the reliance on complex profiles, the design is optimally primed for scalable additive manufacturing, such as multi-jet 3D printing, ensuring its seamless transition into safety critical real world applications. Finally, the successful conceptualization and validation of such advanced physical designs are entirely underpinned by the ultra-fast algebraic formulation developed in this study. Ultimately, this mathematical methodology yields a highly automated, rapidly convergent, and robust paradigm for dispersion analysis, providing a powerful and adaptable engine for the future discovery and optimization of increasingly intricate architected materials.

Synthesizing these computational and structural advancements, the primary findings are outlined below:

• The developed BMS+NSM framework facilitates the automated, high-fidelity evaluation of wave isolation characteristics, achieving a superior trade-off between computational efficiency and numerical precision. The reliability and absolute accuracy of this formulation are unequivocally validated through excellent agreement with alternative schemes (FOM+MT, FOM+NSM, and BMS+MT) across both dispersion band structure calculations and unit cell mode shape predictions.

• The strategic integration of parametric Bezier profiles into the RH architecture profoundly transforms its response. For the proposed CRHM, this innovation naturally triggers the opening of eight distinct, complete band gaps, originating at an exceptionally low onset frequency of $163$ Hz and yielding a remarkable spectral coverage of $69\%$ in the frequency range below $2$ kHz. Crucially, this demonstrates an extraordinary, inherent capacity for deep subwavelength wave isolation achieved strictly through topological formulation.

• Parametric evolution of the cubic Bezier profiles reveals a fundamental trade-off between attenuation capacity and spectral continuity. Aggressive geometric scaling maximizes total band gap coverage and lowers the onset frequency, but fractures the spectrum into discrete, narrow band gaps. Conversely, moderate geometries sacrifice peak coverage to preserve a contiguous, uninterrupted broadband attenuation region. This predictable tunability allows the architecture to be tailored for either absolute isolation bandwidth or continuous filtering.

• The CRHM exhibits exceptional spatial anisotropy, fundamentally transforming how energy propagates through the array. Iso-frequency contours and directivity evaluations verify that the precise manipulation of phase and group velocities enables highly directional energy transmission, allowing engineers to strictly control and route waves along predetermined directions.