Surrogate-model-based design optimization of multiple locally resonant metamaterials for low-frequency vibration reduction

Abstract

This study presents a surrogate-model-based optimal design framework to enhance low-frequency vibration attenuation of multiple locally resonant metamaterials (M-LRMs). The geometric parameters of the M-LRM unit cell are defined as design variables, enabling systematic tuning of resonator configurations. For each band gap, a performance index is defined to reflect the downward shift of the lower band-edge frequency and the increase in band-gap width, and these indices are integrated into a single objective function using a weighted-sum scalarization method. The relationship between the design variables and the performance indices is approximated by a surrogate model, and global optimization is performed in the surrogate design space using a hybrid metaheuristic algorithm (HMA) to identify an optimal unit cell configuration. Dispersion analysis confirms that, compared to the initial design, wider band gaps are formed with reduced lower band-edge frequencies. Frequency response function (FRF) analysis of a finite panel incorporating the optimized M-LRM demonstrates that the vibration response is attenuated within the band-gap frequency ranges predicted by dispersion analysis, validating the proposed framework for finite structures. The proposed approach streamlines the design of M-LRMs for low-frequency, broadband vibration reduction and offers extensibility to diverse unit cell configurations and lightweight panel applications.

Keywords

multiple locally resonant metamaterials band gap vibration reduction surrogate model design optimization

Introduction

Modern industrial structures—including those in the automotive, home appliance, and aerospace sectors—have consistently pursued weight reduction of constituent panels to reduce cost and improve productivity.^1–4 However, lightweight panels are inherently vulnerable to low-frequency vibrations in real operating environments, which can degrade perceived product quality and cause durability issues.^5–8 Conventional mitigation strategies, such as structural redesign and the application of damping materials have been extensively employed to address these problems.^9–11 Nevertheless, these approaches often involve trade-offs in the form of increased mass, added manufacturing complexity, and limited effectiveness over narrow frequency ranges.^12,13 Accordingly, there is a growing demand for an alternative technology capable of suppressing low-frequency vibrations while preserving structural lightness.

In this context, vibration metamaterials have emerged as a promising solution. Vibration metamaterials are artificially engineered structures designed to exhibit dynamic properties not found in natural materials, enabling the formation of band gaps that inhibit elastic-wave propagation within targeted frequency ranges. Among them, locally resonant metamaterials (LRMs) consist of periodically arranged unit cells, in which band gaps emerge in the vicinity of the resonators’ natural frequencies as a result of effective dynamic properties, such as negative effective mass.^14–18 Owing to their reliance on local resonance mechanisms, LRMs can generate low-frequency band gaps with relatively small added mass.^19–22 Moreover, the center frequency and bandwidth of these band gaps can be systematically tailored through geometric design, resonator configuration, and material selection. LRMs can be implemented by attaching them to, or partially embedding them within, existing panel structures without substantial modification of the overall geometry or thickness, making them particularly practical for engineering applications.

Numerous studies have investigated the realization of the band gaps in specific frequency ranges by applying LRMs to panel structures.^23–27 For example, Jung et al.²³ applied LRMs to automotive panels and demonstrated effective noise-and-vibration reduction; while Pires et al.²⁴ employed LRMs to mitigate flow-induced noise and vibration in duct systems. Hyun et al.²⁵ proposed LRMs exploiting defect modes to suppress vibration energy and simultaneously enable energy harvesting in the low-frequency range. In the aerospace domain, Manushyna et al.²⁶ verified the vibration reduction capability of LRMs by applying them to an upper launch-vehicle structure. Additionally, Jung et al.²⁷ demonstrated the feasibility of LRM design optimization for forming band gaps below 500 Hz through topology optimization based on simulated annealing. Despite these advances, most previous studies have focused on single-resonator (mass–spring) configurations, which inherently tend to generate only a single band gap. Although such a band gap can provide strong attenuation, its effectiveness is restricted by the inherently narrow bandwidth.

To address this limitation, multiple locally resonant metamaterials (M-LRMs) have been proposed, in which multiple resonators with distinct dynamic characteristics are integrated within a single unit cell to generate multiple band gaps. Yu and Kook²⁸ introduced an M-LRM in cooperating two resonators within a unit cell for panel vibration reduction and confirmed, through numerical analysis and experiments, the formation of multiple band gaps below 1 kHz. Their resonator design allowed the mass component to be readily inserted into the stiffness component, thereby reducing manufacturing complexity and improving experimental reproducibility. Miao et al.²⁹ presented a metamaterial based on composite cylindrical resonators and employed a genetic algorithm optimization to achieve multiple band gaps below 500 Hz. Nevertheless, practical applications of M-LRMs to structures remain limited, and existing studies have largely been restricted to concept-level demonstrations. In particular, systematic performance optimization tailored to specific low-frequency target ranges has not yet been sufficiently investigated.

In this context, the present study proposes a surrogate-model-based optimal design framework for an M-LRM unit cell aimed at achieving broadband band-gap formation in the low-frequency range. The surrogate model is constructed using a limited number of representative samples, thereby substantially reducing the computational cost associated with repeated dispersion analyses during optimization. To promote low-frequency, wide band gaps, a performance index is defined for each band gap to simultaneously reflect the minimization of the lower band-edge frequency and the maximization of the corresponding band-gap width, and the indices are integrated into a single objective function through weighted-sum scalarization. The design variables are parameterized such that, for each resonator, the mass component is represented by its height, while the stiffness component is represented by the control points of a quadratic Bézier curve that defines its geometric profile. This parameterization enables smooth and physically meaningful geometric variations and allows the resonator geometry to be automatically updated throughout the optimization process. A surrogate model is constructed to approximate the relationship between the design variables and the performance indices, and global optimization is performed using a Hybrid Metaheuristic Algorithm (HMA) coupled with the weighted-sum formulation to efficiently explore the design space. This procedure yields an approximate Pareto-optimal solution, from which a representative optimal design satisfying the target performance criteria is selected. To assess the practical effectiveness of the optimized design, frequency response function (FRF) analysis is conducted on a finite panel with the optimized unit cell periodically arranged. The results demonstrate pronounced vibration attenuation over frequency ranges closely aligned with the predicted band gaps, thereby confirming both the validity and applicability of the proposed design methodology.

The remainder of this paper is organized as follows. Section 2 describes the configuration of the proposed M-LRM unit cell, in which multiple resonators are integrated for panel applications, and examines the associated band-gap characteristics through dispersion analysis based on finite element analysis. Section 3 formulates the surrogate-based optimization problem, including the definition of the objective function, geometric design variables, and relevant constraints. Section 4 presents the design of experiments procedure and details the construction and validation of the surrogate model. Section 5 discusses the optimal design results obtained using the weighted-sum approach and evaluates the vibration attenuation performance of the optimized design through frequency response function analysis of a finite panel. Finally, Section 6 summarizes the main findings and conclusions of the study.

Multiple locally resonant metamaterial (M-LRM)

A multiple locally resonant metamaterial (M-LRM) is an engineered periodic system composed of unit cells that incorporate multiple local resonators with distinct dynamic characteristics. Unlike conventional locally resonant metamaterials based on a single resonator, the presence of multiple resonators within a unit cell enables the formation of several resonance-induced band gaps distributed over different frequency ranges. From the perspective of vibration analysis, these band gaps originate from the localized resonance of individual resonators and effectively inhibit wave propagation in the corresponding frequency intervals, thereby extending the achievable attenuation bandwidth.

In practical implementations, an M-LRM is typically realized as a periodic arrangement of unit cells attached to or integrated with a host structure, where each resonator can be modeled as an equivalent mass–spring subsystem. The collective dynamic interaction between the host structure and the multiple resonators gives rise to effective dynamic properties that cannot be achieved with single-resonator configurations, providing enhanced flexibility in tailoring low-frequency, broadband vibration attenuation characteristics.

Unit cell configuration

The geometry and material configuration of the proposed M-LRM unit cell are illustrated in Figure 1. The material properties assigned to each component of the unit cell are summarized in Table 1, and all materials are assumed to be isotropic. The present unit cell configuration is developed by partially extending and modifying the design originally proposed by Yu and Kook.²⁸ A support structure fabricated from PLA plastic is formed on the top surface of an aluminum host plate, on which two local resonators, denoted as Resonator-A and Resonator-B, are mounted.

Figure 1.

Schematic of the M-LRM unit cell: (a) three-dimensional geometry and material configuration of the host structure and resonators, (b) principal in-plane geometric dimensions on the top surface (unit: mm), and (c) principal in-plane geometric dimensions of the cross-section (unit: mm).

Table 1.

Material properties of the M-LRM unit cell.

Material	Density (kg/m³)	Young’s modulus (GPa)	Poisson’s ratio
PLA plastic	1400	4.16	0.4
Steel	7900	200	0.4
Aluminum	2738	68	0.38

The overall three-dimensional geometry and material arrangement of the unit cell are shown in Figure 1(a). Each resonator is composed of a mass component, modeled as a steel block, and a stiffness component, modeled as a PLA frame, such that the resonance frequency is governed by the combined mass–stiffness characteristics. By employing steel for the mass component instead of PLA, both the density and Young’s modulus increase for an identical geometry and volume. Consequently, the deformation of the mass component becomes negligible compared to that of the stiffness component, and the lumped-mass assumption is well justified. This material choice shifts the resonance frequency toward the low-frequency range and enhances the concentrated-mass effect, which is advantageous for forming band gaps at low frequencies compared to configurations using PLA for the mass component.

To account for the dynamic interaction between the two resonators, their geometries and dimensions are intentionally designed to be different. The principal geometric dimensions of the top surface and the cross-sectional view are presented in Figure 1(b) and (c), respectively. This asymmetric design leads to a separation of the resonance frequencies of Resonator-A and Resonator-B, thereby inducing the formation of multiple band gaps within the target frequency range. Accordingly, by appropriately tuning the mass and stiffness parameters of the resonators, the M-LRM unit cell can be designed to generate band gaps at prescribed frequencies. When such an M-LRM configuration is applied to a panel structure, effective vibration attenuation in the low-frequency range can be achieved.

Identified band gap of M-LRM unit cell

The band gap of a locally resonant metamaterial is identified from the dispersion relation $ω (k)$ , which describes the relationship between the angular frequency $ω$ and the wave vector $k = (k_{x}, k_{y})$ .^23,27,30,31 The M-LRM is modeled as an infinite periodic array, and the displacement field satisfies the Bloch–Floquet theorem.^23,27,31 Accordingly, the Bloch–Floquet periodic boundary conditions are imposed on the unit cell boundaries, as illustrated in Figure 2. These boundary conditions are applied to the opposing unit cell boundaries along the $x$ and $y$ directions in accordance with prescribed wave vector $k$ , such that the dynamic response of a single unit cell represents that of the entire periodic structure. In the present study, material damping was not included in the baseline dispersion analysis of the unit cell, as the primary objective was to clearly identify the fundamental band-gap formation mechanism induced by local resonance at the conceptual design stage. Then, the resulting governing equation is expressed as the following generalized eigenvalue problem:

[K (k) - ω^{2} M] u = 0,

(1)

where $K (k)$ is the wave vector dependent stiffness matrix from the application of Bloch–Floquet periodic boundary conditions and $M$ is the mass matrix of the unit cell, and $u$ is the discretized displacement eigenvector.³¹ The mass and stiffness matrices in equation (1) are assembled from the finite-element discretization of the Unit Cell under Bloch–Floquet boundary conditions. The analysis was implemented in COMSOL Multiphysics. The formulation follows a well-established finite-element framework widely used in band-structure analysis of periodic systems. Further details can be found in the related literature.^23,32 The displacement field follows the Bloch form

u (r, t) = \tilde{u} (r) e^{- i k \cdot r} \cdot e^{i ω t},

(2)

where $r = (x, y)$ denotes the position vector and $\tilde{u}$ represents the periodic mode shape of the unit cell. For a given k, solving equation (1) yields the corresponding eigenfrequencies $ω$ that satisfy equation (2). In principle, the dispersion relation must be evaluated by sweeping $k$ over the entire Brillouin zone. However, when the lattice symmetry is present, the analysis can be restricted to the irreducible Brillouin zone (IBZ) without loss of generality.^23,27 The shape of the IBZ is determined by the plane crystallographic group defined by rotational, mirror, and glide symmetries, in addition to the geometry of the underlying two-dimensional Bravais lattice.³³ Although the unit cell is based on a square Bravais lattice, the unit cell motif—specifically the geometry and spatial arrangement of the two resonators—does not fully satisfy rotational and mirror symmetries. As a result, the conventional triangular IBZ path commonly adopted for a square lattice ( $Γ \to M \to X \to Γ$ ) is not applicable in this case. Therefore, the square-lattice IBZ path shown in Figure 2 ( $X \to Γ \to M \to X \to T \to Γ$ ) is employed, and the dispersion curves are obtained by sweeping $k$ along this path. In the present study, the IBZ path was discretized using 61 wave-vector points, and equation (1) was solved at each point to obtain the corresponding eigenfrequencies. The resulting dispersion characteristics are presented in Figure 3.

Figure 2.

Bloch–Floquet periodic boundary conditions applied to the unit cell and the corresponding irreducible Brillouin zone (IBZ) path for a square Bravais lattice ( $X \to Γ \to M \to X \to T \to Γ$ ).

Figure 3.

Representative results along the IBZ path: (a) dispersion curves of the M-LRM unit cell, (b) mode shapes at the lower and upper band-edge frequencies of the first band gap ( $f_{1}^{low} = 461$ Hz, $f_{1}^{high} = 481$ Hz), (c) mode shapes at the lower and upper band-edge frequencies of the second band gap ( $f_{2}^{low} = 618$ Hz, $f_{2}^{high} = 709$ Hz).

Figure 3(a) presents the dispersion curve computed along the IBZ path, while Figure 3(b) and (c) show representative mode shapes at the band gap edges. The set of $(k, f)$ points corresponding to real eigenfrequencies on the dispersion curves defines the frequency ranges that allow wave propagation. In contrast, continuous frequency intervals in which no dispersion branch exists for all $k$ are identified as band gaps. The computed results indicate that the first band gap spans approximately 461–481 Hz ( $Δ f \approx 20$ Hz), while the second band gap extends from approximately 618 to 709 Hz ( $Δ f \approx 91$ Hz). Furthermore, inspection of the mode shapes at the band gap edges reveals distinct vibration characteristics. At the lower band-edge frequencies ( $f_{1}^{low}, f_{2}^{low}$ ), the responses are dominated by local resonant modes primarily associated with Resonator-B and Resonator-A, respectively. In contrast, the upper band-edge frequencies ( $f_{1}^{high}, f_{2}^{high}$ ) exhibit hybridized modes involving strong coupling between the resonators and the host plate. Such resonance-induced modal interactions constitute the primary physical mechanism responsible for the formation of multiple band gaps in the proposed M-LRM unit cell.

More specifically, the dispersion branches near the two band gaps in Figure 3(a) further clarify the stop-band formation mechanism. The relatively flat branches near the lower band edges indicate reduced group velocity, suggesting that the vibration energy is predominantly localized in the resonator subsystem rather than being efficiently transmitted through the host plate. By contrast, the upper band edges are associated with modal coupling between the resonators and the host plate, indicating that the deformation is distributed more broadly over the unit cell. In addition, the branches around the second band gap remain separated over a wider frequency interval than those around the first band gap, consistent with the larger width of the second band gap. These observations further support the multiple-band-gap behavior of the proposed M-LRM unit cell.

Frequency response analysis of a finite M-LRM panel

As discussed in Section 2.2, the M-LRM is idealized as an infinite periodic system, whereas practical implementations necessarily involve finite-sized structures. To validate the band gaps predicted by the dispersion analysis, the frequency response function (FRF) was evaluated on a finite panel composed of a 5 × 5 array of unit cells, as illustrated in Figure 4.¹⁸ For computational efficiency, symmetry boundary conditions were imposed along the left and bottom edges (solid lines), allowing the analysis domain to be reduced to one quarter of the full panel. Fixed boundary conditions were applied to the right and top edges (dashed lines). A point load with an amplitude of 1 N was applied at point A, and the displacement response was measured at point B to evaluate the FRF. Consistent with the baseline dispersion analysis in Section 2.2, material damping was not included in the present FRF analysis. In addition, within a linear frequency-response framework, variations in excitation amplitude affect the response magnitude proportionally, but do not directly alter the band-edge frequencies or band-gap widths. Accordingly, the present finite-panel analysis aims to verify whether the band-gap frequency ranges predicted by the dispersion analysis are consistently manifested as vibration-attenuation ranges linear, small-amplitude conditions. It should be noted, however, under high-amplitude vibration conditions, geometric and/or material nonlinearities of the resonators may modify the band-gap characteristics; such nonlinear effects are beyond the scope of the present study.

Figure 4.

Geometry of the finite panel with a 5 × 5 array of unit cells, including boundary and loading conditions for the frequency response analysis.

The FRF results shown in Figure 5 demonstrate that the displacement response is significantly attenuated over frequency ranges that closely coincide with the first and second band gaps predicted by the dispersion analysis in Figure 3. This agreement confirms that the band-gap-based vibration reduction mechanism remains effective when the M-LRM is implemented in a finite panel configuration.

Figure 5.

Frequency response function measured at point B. The gray shaded regions denote the first and second band-gap frequency intervals identified from the dispersion analysis.

Formulation of the optimization problem

This section formulates the optimization problem for the M-LRM with the aim of achieving broadband band-gap formation in the low-frequency range. Based on the unit-cell configuration introduced in Section 2, the objective function and geometric design variables are defined, and the overall optimization problem is systematically formulated.

Definition of the band-gap performance index

To enhance vibration attenuation performance in the low-frequency range, a band-gap performance index is defined based on a normalized band-gap bandwidth that simultaneously accounts for both the lower band-edge frequency and the width of each band gap. Let $f_{n}^{low}$ denote the lower band-edge frequency of the n-th band gap and let $Δ f_{n}$ denote its corresponding band-gap width; the band-gap performance index ${BW}_{n}$ is defined as

{BW}_{n} = \frac{Δ f_{n}}{f_{n}^{low}}, n = 1, 2 .

(3)

Maximizing ${BW}_{n}$ directly reflects the design objective of achieving low-frequency, broadband vibration attenuation. Specifically, ${BW}_{n}$ can be increased by decreasing $f_{n}^{low}$ or increasing $Δ f_{n}$ , thereby shifting the band gap toward lower frequencies and broadening the bandwidth. Consequently, the maximization of ${BW}_{n}$ naturally balances these two competing effects and provides a unified measure for evaluating band-gap performance in the low-frequency regime. Since the unit cell considered in this study incorporates two local resonators, two band gaps that are relevant within the target low-frequency range ( $n = 1, 2$ ) are selected as optimization targets.

Definition of design variables, Bézier parameterization, and sensitivity analysis

The design variables are defined as geometric parameters associated with the mass and stiffness components of the resonators. In this subsection, the parameterization of the resonator mass and stiffness components is first introduced, followed by the definition of the corresponding design variables and a sensitivity analysis performed at the initial design point.

Figure 6 illustrates the geometric design regions assigned to the mass components of Resonator-A and Resonator-B. The dashed regions indicate the allowable design domains for the resonator mass components. The design variables are defined as the heights of the mass components, denoted by $M^{A}$ and $M^{B}$ , respectively. The bottom contact surface of each mass component is constrained to remain geometrically coincident with the top contact surface of the PLA support structure on the $x$ – $y$ plane. As a result, only geometric variation in the positive z-direction is permitted.

Figure 6.

Design regions for the resonator mass components of Resonator-A and Resonator-B.

Accordingly, the design-variable vector associated with the mass components is defined as

x_{M} = [M^{A}, M^{B}]^{T},

(4)

which is subject to the following bound constraints:

x_{M}^{L} {\leq x}_{M} \leq x_{M}^{U} .

(5)

Accordingly, the optimization adjusts the design variable vector $x_{M}$ only within the predefined allowable range $[x_{M}^{L}, x_{M}^{U}]$ .

The cross-sectional profile of the resonator stiffness component is modeled using a quadratic Bézier curve. The fundamental characteristics of this curve are briefly summarized. As illustrated in Figure 7, a quadratic Bézier curve is defined by three control points $(P_{1}, P_{2}, and P_{3})$ and passes exactly through the two endpoint control points $P_{1}$ and $P_{3}$ . The intermediate control point $P_{2}$ governs the tangent directions at the endpoints as well as the overall curvature of the curve. Owing to the convex-hull property, the Bézier curve entirely lies within or on the boundary of the triangle formed by the three control points, ensuring stable and robust shape control against perturbations of the control points. These properties are exploited to define the upper and lower boundaries of the stiffness component, and a selected subset of the control points is parameterized as design variables. The detailed parameterization scheme and the associated boundary constraints are presented in Figure 8 and equations (6)–(11).

Figure 7.

Formation principle of a quadratic Bézier curve.

Figure 8.

Geometric design regions and Bézier curve parameterization of the stiffness components for Resonator-A and Resonator-B.

For each resonator $X \in (A, B)$ , the upper boundary of the stiffness component is defined by a quadratic Bézier curve specified by three control points. The control points are denoted as

p_{i}^{X} = {(x_{i}^{X}, y_{i}^{X})}^{T}, X \in (A, B), i = 1, 2, 3 .

(6)

The control points are parameterized as

\begin{matrix} p_{1}^{A} = {[0, 0]}^{T}, p_{2}^{A} = p_{1}^{A} + {[- 1.5, 0]}^{T}, \\ p_{3}^{A} = p_{2}^{A} + {[- S_{x}^{A}, S_{y}^{A}]}^{T} . \end{matrix}

(7)

\begin{matrix} p_{1}^{B} = {[0, 0]}^{T}, p_{2}^{B} = p_{1}^{B} + {[2.5, 0]}^{T}, \\ p_{3}^{B} = p_{2}^{B} + {[S_{x}^{B}, S_{y}^{B}]}^{T} . \end{matrix}

(8)

Here, $p_{1}^{X}$ is fixed at the origin; $p_{2}^{X}$ is obtained by translating $p_{1}^{X}$ by a prescribed distance along the $x$ -direction. The third control point $p_{3}^{X}$ is displaced from $p_{2}^{X}$ by the design-variable components $S_{x}^{X}$ and $S_{y}^{X}$ in the x- and y-directions, respectively. This parameterization ensures that variation in $p_{3}^{X}$ primarily governs the curvature and overall shape of the Bézier curve, while preserving geometric continuity and robustness.

The lower boundary of the stiffness component is defined by a set of copied control points obtained by translating each upper-boundary control point uniformly in the negative $y$ -direction by the thickness design-variable component $S_{h}^{X}$ , that is,

p_{i, copy}^{X} = p_{i}^{X} - {[0, S_{h}^{X}]}^{T}, i = 1, 2, 3 .

(9)

Accordingly, the upper and lower boundaries remain mutually parallel, and the thickness of the stiffness components is maintained at a constant value $S_{h}^{X}$ over the entire design region. The endpoints of the two curves are connected by straight line segments to form a closed cross-sectional profile. Based on this parameterization, the stiffness-related design variable vector for resonator $X$ is defined as

x_{S}^{X} = {[S_{x}^{X}, S_{y}^{X}, S_{h}^{X}]}^{T}, X \in (A, B),

(10)

which is subject to the following bound constraints:

{(x_{S}^{X})}^{L} \leq x_{S}^{X} \leq {(x_{S}^{X})}^{U} .

(11)

Accordingly, the optimization process adjusts $x_{S}^{X}$ only within the predefined admissible range $[(x_{S}^{X})^{L}, (x_{S}^{X})^{U}]$ . With this compact set of variables $(S_{x}^{X}, S_{y}^{X}, S_{h}^{X})$ , both the curvature and thickness of the stiffness component can be controlled in a stable and robust manner, thereby enabling efficient exploration of the design space during optimization.

The sensitivities were evaluated by the central difference method while fixing all other variables. The initial values and admissible bounds of the design variables are listed in Table 2. To enable a consistent comparison among variables with different admissible intervals, the perturbation size for each variable was set to 1%, 3%, 5%, 10%, and 15% of its full admissible range. The sensitivities were then nondimensionalized using the admissible range of each variable and the initial value of the corresponding performance index, such that the sign indicates the direction of influence, and the magnitude represents the relative importance near the initial design point.

Table 2.

Design variables and their admissible bounds (lower and upper, unit: mm).

Design variables	Initial value (mm)	Lower limit (mm)	Upper limit (mm)
$S_{x}^{A}$	1.5	0.5	2.5
$S_{y}^{A}$	0	−2	2
$S_{h}^{A}$	2	1.5	2.5
$M^{A}$	6	4	8
$S_{x}^{B}$	2.5	1.5	3.5
$S_{y}^{B}$	0	−2	2
$S_{h}^{B}$	2	1.5	2.5
$M^{B}$	7	5	9

Figure 9 summarizes the dimensionless sensitivities of the design variables. A positive sensitivity indicates that increasing the corresponding design variable improves the band-gap performance index, whereas a negative sensitivity indicates deterioration. For ${BW}_{1}$ , $S_{h}^{A}$ , $S_{x}^{B}$ , and $M^{B}$ exhibit positive contributions, while $S_{x}^{A}$ , $M^{A}$ , and $S_{h}^{B}$ show negative effects. The sensitivities of $S_{y}^{A}$ and $S_{y}^{B}$ are relatively smaller. For ${BW}_{2}$ , $S_{x}^{A}$ , $M^{A}$ , and $S_{y}^{B}$ generally have positive effects, whereas $S_{h}^{A}$ and $M^{B}$ are consistently negative. The sensitivities $S_{y}^{A}, S_{x}^{B}$ , and $S_{h}^{B}$ remain close to zero over most of the perturbation range. Although the sensitivity magnitudes vary slightly with perturbation size, the overall signs and relative ranking remain largely unchanged, indicating stable local trends near the initial design point. Notably, most design variables exhibit opposite sensitivity signs for ${BW}_{1}$ and ${BW}_{2}$ , revealing an inherent trade-off between the two performance indices. This observation supports the need for simultaneous multi-objective optimization.

Figure 9.

Dimensionless sensitivities of the band-gap performance indices with respect to the design variables: (a) ${BW}_{1}$ and (b) ${BW}_{2}$ .

Formulation of the optimization problem

Based on the performance indices defined in Section 3.1 and the design variables specified in Section 3.2, the optimization problem for simultaneously enhancing the performance of two band gaps is formulated as follows.

\begin{matrix} Find : x = (S_{x}^{A}, S_{y}^{A}, S_{h}^{A}, M^{A}, S_{x}^{B}, S_{y}^{B}, S_{h}^{B}, M^{B})^{T} \in Ω, \\ Ω = [x^{L}, x^{U}] \\ \max_{x \in Ω} : J_{i} (x) = α_{i} {BW}_{1} (x) + (1 - α_{i}) {BW}_{2} (x), i = 1, \dots 9 \\ Subject to : α_{i} \in {0.1, \dots, 0.9} \end{matrix}

(12)

Here, equation (12) defines the complete optimization problem, including the design-variable vector $x,$ the admissible design region $Ω$ , the weighted-sum objective function, and the associated weighting scheme. The admissible design region $Ω$ is specified by the lower and upper bounds $x^{L}$ and $x^{U}$ listed in Table 2. The objective function is constructed via a weighted-sum scalarization of the two band-gap performance indices ${BW}_{1}$ and ${BW}_{2}$ , and, for each weighting factor $α_{i}$ , the resulting single-objective function is maximized independently. The discrete set of weighting factors employed in equation (12) is used to explore the trade-off between the two band-gap performance indices. Accordingly, the weighted-sum approach is adopted as a practical scalarization method to identify representative compromise solutions. Since the optimization is performed for a discrete set of weighting factors, the resulting solutions should be interpreted as representative Pareto-optimal designs rather than complete coverage of all Pareto-optimal regions. The optimization search is initialized and constrained within the admissible region $Ω$ defined in Table 2.

Through this unified formulation, the optimization simultaneously promotes a reduction in the lower band-edge frequencies and an expansion of the corresponding band-gap widths, thereby yielding M-LRM designs that maximize band-gap performance in the targeted low-frequency range.

Surrogate modeling framework

As formulated in Section 3, the optimization problem requires repeated evaluation of the unit cell responses at multiple wave-vector points $k$ along the IBZ path for each design vector $x$ . Since a single dispersion analysis is computationally expensive, the total computational cost increases rapidly during the optimization process, particularly in the design-space exploration stage where a large number of evaluations are required.

To alleviate this computational burden, a surrogate model is constructed to provide a low-cost and efficient approximation of the input–output relationship between $x$ and the band-gap performance indices $({BW}_{1} (x), {BW}_{2} (x))$ . The surrogate model is constructed using a set of representative samples generated via a design of experiments (DOE) strategy, enabling it to capture the global trends of the underlying design space. When integrated into the optimization framework, the surrogate model allows rapid evaluation of objective-function values for a large number of candidate designs, thereby substantially reducing the number of expensive full-order analyses. As a result, the overall computational cost is significantly decreased, while enabling a denser and more effective exploration of the design space within a fixed computational budget.

This section describes the construction and validation of the Kriging surrogate model employed in the optimization framework. The sampling strategy and the generation of representative training samples using design of experiments (DOE) are first introduced. The formulation and construction of the Kriging-based surrogate model based on the sampled data are then described. Finally, the predictive performance of the surrogate model is assessed using quantitative error metrics, thereby confirming its suitability for application in the subsequent global optimization process.

All implementations of the surrogate model were carried out in the PIAnO software environment.³⁴

Design of Experiments (DOE)

To construct an accurate surrogate model with a limited number of high-fidelity analyses, an appropriate sampling of the design space is required. In the present study, the underlying numerical analyses are deterministic, as the unit-cell responses are obtained from finite-element-based dispersion calculations. As a result, random measurement noise is negligible, and the primary source of error in surrogate modeling arises from approximation bias rather than statistical variance.³⁵

Under these conditions, sampling strategies developed for the Design and Analysis of Computer Experiments (DACE) are well suited, as they aim to provide uniform coverage of the admissible design space. Accordingly, this study employs Optimal Latin Hypercube Sampling (OLHS) to generate the training samples for surrogate model construction.^35,36 The sampled design points are distributed within the prescribed bounds of the design variables listed in Table 2, ensuring adequate coverage of the design space for subsequent surrogate modeling. The resulting DOE samples serve as the basis for constructing the Kriging surrogate model described in the following subsection.

Construction of a Kriging-based surrogate model

To efficiently approximate the relationship between the design variables and the band-gap performance indices, a Kriging-based surrogate model is employed in this study. Kriging, which is rooted in Gaussian process regression, has been widely used in engineering optimization owing to its strong capability to represent nonlinear input–output relationships and its high predictive accuracy with a limited number of training samples.^36–41 A detailed theoretical background and mathematical formulation of the Kriging model can be found in the referenced literature and are therefore not repeated here for brevity.

In the present framework, the Kriging surrogate model is constructed using a training dataset consisting of DOE-generated samples and the corresponding $({BW}_{1}, {BW}_{2})$ obtained by post-processing the results of high-fidelity finite-element dispersion analyses performed for each sample. The surrogate model approximates the mapping from $x$ to ${BW}_{1}$ and ${BW}_{2}$ , thereby enabling rapid evaluation of the objective function without resorting to repeated full-order analyses.

Ordinary Kriging with a constant trend is adopted, and the model hyperparameters are identified using standard likelihood-based estimation procedures. Once constructed, the Kriging surrogate provides smooth and reliable predictions of band-gap performance over the admissible design space. This predictive capability allows the surrogate model to be effectively integrated into the subsequent global optimization process, in which a large number of candidate designs must be evaluated efficiently. The predictive performance of the constructed Kriging surrogate model is quantitatively assessed in the following subsection to verify its suitability for use in the optimization framework.

Evaluation of the predictive performance of the Kriging surrogate

The predictive performance of the Kriging surrogate model is assessed using the coefficient of determination $R^{2}$ and the root mean square error (RMSE). These two metrics are widely used to evaluate both the fidelity and generalization capability of surrogate models and are defined as:

R^{2} = 1 - \frac{\sum_{i = 1}^{m} {(y (x_{i}) - \hat{y} (x_{i}))}^{2}}{\sum_{i = 1}^{m} {(y (x_{i}) - \bar{y} (x_{i}))}^{2}},

(13)

RMSE = \sqrt{\frac{1}{m} \sum_{i = 1}^{m} [y (x_{i}) - \hat{y} (x_{i} {)]}^{2}} .

(14)

Here, $m$ is the number of test samples; $x_{i}$ is the $i$ -th input vector; $y (x_{i})$ is the reference response obtained from finite-element-based dispersion analysis (i.e. eigenvalue analysis based on equation (1)); $\hat{y} (x_{i})$ is the corresponding surrogate prediction; and $y (x_{i})$ is the mean value of the test responses. In general, a value of $R^{2}$ closer to 1 indicates higher predictive capability of the surrogate model, whereas a smaller RMSE reflects higher predictive accuracy with reduced residual errors. Acceptable thresholds for these metrics depend on the specific problem characteristics and prior experience. In the present study, the surrogate model is regarded as sufficiently accurate when $R^{2}$ exceeds approximately 0.85 and RMSE is on the order of $10^{- 2}$ .³⁶

The validation procedure is as follows. First, an initial set of DOE samples is generated within the admissible design region $Ω$ using OLHS. Among these samples, 10% are randomly selected as test points (validation set), while the remaining 90% are used to construct the Kriging surrogate model. The predictive performance metrics $R_{Test}^{2}$ and ${RMSE}_{Test}$ are then computed at the test points to evaluate generalization performance. If the prescribed accuracy criteria are not sufficiently satisfied, additional samples are adaptively added within $Ω$ while maintaining the test-point ratio, and the surrogate model is reconstructed and re-evaluated. This procedure is repeated until the predictive performance meets the prescribed tolerance.

Table 3 summarizes the predictive performance of the final Kriging surrogate model constructed with a total of 150 DOE samples. As shown, the surrogate model achieves $R_{Test}^{2}$ = 0.8602 for ${BW}_{1}$ and $R_{Test}^{2}$ = 0.9030 for ${BW}_{2}$ , while the corresponding RMSE values remain below 0.04. These results indicate good agreement between the surrogate predictions and the high-fidelity finite-element dispersion analyses, confirming that the surrogate model provides sufficient accuracy for use in the subsequent optimization process.

Table 3.

Predictive performance of the Kriging surrogate model for ${BW}_{1}$ and ${BW}_{2}$ with DOE size 150.

$DOE : 150$	$R_{Test}^{2}$ (-)	${RMSE}_{Test}$ (-)
${BW}_{1}$ (-)	0.8602	0.0321
${BW}_{2}$ (-)	0.9030	0.0314

It should be noted that $R^{2}$ and $RMSE$ provide useful global measures of surrogate accuracy over the design space; however, they do not fully capture local predictive behavior in highly sensitive regions such as band edges. Therefore, in the present study, these metrics are not used in isolation but are complemented by additional validation procedures. In particular, surrogate predictions are compared with direct FEA results at selected Pareto-optimal designs, and a region of interest (ROI)-based adaptive refinement is employed to reduce local prediction errors. This combined approach ensures reliable predictive performance for optimization.

Global optimization using a surrogate-assisted metaheuristic algorithm

The Kriging-based surrogate model constructed in Section 4 is employed as a computationally efficient evaluation engine, and the optimization problem formulated in Section 3.3 is solved using a global optimization strategy. Specifically, a Hybrid Metaheuristic Algorithm (HMA) is adopted to explore the surrogate-based objective landscape efficiently. The HMA integrates differential evolution (DE) and cuckoo search (CS) within a bi-population framework, allowing design variables to be handled directly in the continuous domain without additional encoding. This real-valued formulation simplifies implementation and enables the simultaneous treatment of multiple continuous design variables.⁴² In the present study, rank-iMDDE (an improved constrained DE) is combined with modified CS (mCS), yielding a balanced search behavior between exploration and exploitation. This hybrid configuration enhances population diversity and numerical stability, while improving convergence speed and solution quality, as demonstrated in prior studies.^43,44 By coupling the surrogate model with the HMA, the optimization can efficiently identify high-quality candidate designs with substantially reduced computational cost compared to direct high-fidelity optimization. The use of Kriging and HMA in the present study reflects a practical, problem-oriented methodological choice, rather than a claim of universal superiority over alternative approaches.

Validation at a selected Pareto-optimal design

A selected Pareto-optimal design is chosen from the Pareto solution set obtained using the HMA, and the predictions of the Kriging surrogate model at this design are compared with the FEA-based reference results to evaluate the relative errors of the band-gap performance indices. The selected Pareto-optimal design corresponds to the weighting factor $α_{5} = 0.5$ , which represents a balanced trade-off between the two band-gap performance indices, ${BW}_{1} (x)$ and ${BW}_{2} (x)$ .

Table 4 summarizes the band-gap performance indices and their constituent components for the initial design and the selected Pareto-optimal design. Compared with the initial design, both band gaps in the selected Pareto-optimal design exhibit lower band-edge frequencies ( $f_{1}^{low} : 461 \to 344 Hz$ , $f_{2}^{low} : 619 \to 519 Hz$ ) and substantially increased bandwidths ( $Δ f_{1} : 20 \to 50 Hz$ , $Δ f_{2} : 91 \to 168 Hz$ ). Consequently, both ${BW}_{1}$ and ${BW}_{2}$ increase significantly, indicating that the optimization objective of achieving wider band gaps in the low-frequency range has been successfully attained.

Table 4.

Comparison of band-gap performance indices and their components between the initial design and the optimal design ( $α_{5} = 0.5$ ).

$α_{5} = 0.5$	Initial design	Optimal design ( $α_{5}$ )
${BW}_{1}$ (-)	0.0432	0.1442
$f_{1}^{low}$ (Hz)	461	344
$Δ f_{1}$ (Hz)	20	50
${BW}_{2}$ (-)	0.1465	0.3231
$f_{2}^{low}$ (Hz)	619	519
$Δ f_{2}$ (Hz)	91	168

The relative error $δ_{n}$ between the surrogate-model predictions and the FEA-based reference values at the selected Pareto-optimal design is defined as follows:

δ_{n} = | \frac{{BW}_{n, Kriging} (x) - {BW}_{n, Original} (x)}{{BW}_{n, Original} (x)} | \times 100 %, n \in {1, 2} .

(15)

Where ${BW}_{n, Kriging} (x)$ denotes the Kriging surrogate prediction and ${BW}_{n, Original} (x)$ denotes the reference values obtained from the FEA-based dispersion analysis. In the present study, when $δ_{n}$ converges to within 5%, the corresponding prediction error is regarded as being within an acceptable range.

As reported in Table 5, the selected Pareto-optimal design $(α_{5} = 0.5$ ) yields relative errors of $δ_{1} = 18.31 %$ and $δ_{2} = 5.08 %$ . Although the global accuracy criteria based on the prediction metrics (RMSE and R²) presented in Section 4.3 are satisfied, relatively large local prediction errors are observed at this specific design point. This behavior can be attributed to the pronounced nonlinearity and non-stationary characteristics of the response in the vicinity of band-edge regions. Near the band edges, strong modal interactions—such as mode veering and branch crossing—lead to sharp variations in the slope and curvature of the dispersion characteristics. As a result, the response becomes highly sensitive to small perturbations in the design variables. Consequently, even when the surrogate model exhibits good global predictive performance, localized prediction errors $δ_{n}$ may remain relatively large in these highly sensitive regions.

Table 5.

Relative error $δ_{n}$ between the Kriging surrogate predictions and the FEA-based reference values at the selected Pareto-optimal design ( $α_{5} = 0.5$ ).

$α_{5} = 0.5$	${BW}_{n, Kriging} (x) [-]$	${BW}_{n, Original} (x) [-]$	$δ_{n} [%]$
${BW}_{1}$	0.0853	0.0721	18.31
${BW}_{2}$	0.1697	0.1615	5.08

Motivated by this observation, a region of interest (ROI) is defined based on the empirical distribution of the Pareto-optimal solutions, and additional infill samples are introduced within this region to further refine the surrogate model.^39–41 Specifically, for each discrete weighting factor $α_{i}$ , the optimal design-variable vector $x^{(i)}$ that maximizes the weighted-sum objective function $J_{i} (x)$ is defined as follows:

x^{(i)} = \arg \max_{x \in Ω} J_{i} (x), i = 1, \dots, 9

(16)

For each design variable $x_{j}$ ( $j = 1, \dots, 8$ ), the minimum and maximum values among the Pareto-optimal designs are subsequently identified to delineate the ROI for adaptive sampling.

x_{j}^{L, ROI} = \min_{i = 1, \dots, 9} x_{j}^{(i)}, x_{j}^{U, ROI} = \max_{i = 1, \dots, 9} x_{j}^{(i)}, j = 1, \dots, 8 .

(17)

Based on these bounds, the ROI design domain is defined as

Ω_{ROI} = [x \in Ω | x_{j}^{L, ROI} \leq x_{j} \leq x_{j}^{U, ROI}, j = 1, \dots, 8] .

(18)

Here, $Ω_{ROI}$ represents the minimal admissible region that encloses all Pareto-optimal solutions. Subsequent OLHS infill sampling and surrogate reconstruction are therefore performed exclusively within $Ω_{ROI}$ .

As illustrated in Figure 10, the vertical range bars indicate the admissible intervals each design variable, while the filled circular markers denote the corresponding lower and upper bounds. For all design variables, the ROI bounds become noticeably narrower than the initial bounds, resulting in a more compressed admissible design space. By concentrating additional samples in this reduced region, the effective sample density within $Ω_{ROI}$ increases, which enables the surrogate model to better capture localized variations associated with band-edge behavior and strong local nonlinearity. Consequently, the prediction bias observed near sensitive band-edge regions is alleviated, leading to improved accuracy and reliability of the surrogate model. The overall procedure up to this point is summarized in the flowchart shown in Figure 11.

Figure 10.

Comparison of the initial bounds $Ω$ and the ROI bounds $Ω_{ROI}$ for each design variable.

Figure 11.

Workflow of surrogate-assisted global optimization with prediction-metric verification and ROI-based adaptive OLHS infill sampling.

Table 6 summarizes the prediction performance of the Kriging surrogate after increasing the DOE size to 189. Starting from the initial 150 samples reported in Table 3, an additional 30 infill samples were generated within the ROI. Moreover, nine Pareto-optimal designs obtained from the global optimization results were included.

Table 6.

Prediction performance of the Kriging surrogate for ${BW}_{1}$ and ${BW}_{2}$ after ROI-based infill sampling (DOE size = 189).

$DOE : 189$	$R_{Test}^{2}$ (-)	${RMSE}_{Test}$ (-)
${BW}_{1}$	0.8796	0.0278
${BW}_{2}$	0.9034	0.0278

For model validation, the test set consists of the same 15 validation samples separated from the initial DOE, supplemented by 20% (six samples) of the ROI infill points, yielding a total of 21 test samples. The resulting prediction metrics are $R_{Test}^{2} = 0.8796$ , ${RMSE}_{Test} = 0.0278$ for ${BW}_{1}$ , and $R_{Test}^{2} = 0.9034$ , ${RMSE}_{Test} = 0.0278$ for ${BW}_{2}$ , These results indicate improved robustness and stability of the surrogate model compared with the initial DOE-based construction. In the subsequent re-optimization, to distinguish the process from the initial optimization performed over the full design bounds in Table 2, the initial design vector is set to the arithmetic mid-point of the ROI bounds, that is,

x_{j}^{(0)} = \frac{1}{2} (x_{j}^{L, ROI} + x_{j}^{U, ROI}) .

(19)

This setting reduces bias toward the boundaries of $Ω_{ROI}$ and improves the stability of the initial exploration phase.

Table 7 summarizes $δ_{n}$ for the Kriging surrogate predictions and the FEA-based reference results at $α_{5} = 0.5$ after re-optimization. Following ROI-based infill sampling and surrogate reconstruction, $δ_{1}$ and $δ_{2}$ are reduced to 3.04% and 2.27%, respectively, and both are below 5%, indicating convergence within the acceptable error tolerance. These results suggest that focused infill sampling within the ROI effectively mitigates localized prediction errors and enhances the reliability of the surrogate model in the vicinity of the Pareto-optimal solutions.

Table 7.

Relative error $δ_{n}$ between the Kriging surrogate prediction and the FEA reference at the selected Pareto-optimal design ( $α_{5} = 0.5$ ) after ROI-based re-optimization.

$α_{5} = 0.5$	${BW}_{n, Kriging} (x) [-]$	${BW}_{n, Original} (x) [-]$	$δ_{n} [%]$
${BW}_{1}$	0.1290	0.1252	3.04
${BW}_{2}$	0.3642	0.3561	2.27

The weighted-sum based multi-objective optimization results

The weighted-sum objective functions $J_{i} (x)$ are employed to obtain optimal solutions corresponding to different weight combinations $α_{i}$ , thereby approximating the Pareto front in the ${BW}_{1}$ – ${BW}_{2}$ objective space. Figure 12 illustrates the resulting Pareto front together with the locations of the optimal solutions obtained for each $J_{i} (x)$ . Table 8 summarizes the corresponding component metrics of each optimum solution and their variations relative to the initial design. Here, the term initial design refers to the baseline configuration presented in Table 2 and is distinguished from the initial values set during the subsequent ROI-based resampling and re-optimization process. This convention is adopted consistently throughout this section, as the purpose of the comparative analysis is to quantify the performance improvement from the original baseline to the final optimized designs.

Figure 12.

Pareto front on the ${BW}_{1}$ – ${BW}_{2}$ plane and distributions of the optima of $J_{i} (x)$ (including the initial design).

Table 8.

Comparison of component metrics at the optima of $J_{i} (x)$ (including changes relative to the initial design).

Design case	$f_{1}^{low}$ (Hz)	$f_{2}^{low}$ (Hz)	$Δ f_{1}$ (Hz)	$Δ f_{2}$ (Hz)	$Δ f_{1} + Δ f_{2}$ (Hz)
Initial	461	619	20	91	111
$J_{1} (x)$	394 (−14.53%)	469 (−24.23%)	8 (−60%)	189 (107.69%)	197 (77.48%)
$J_{2} (x)$	395 (−14.32%)	478 (−22.78%)	13 (−35%)	191(109.89%)	204 (83.78%)
$J_{3} (x)$	373 (−19.09%)	491 (−20.68%)	43 (115%)	182 (100%)	225 (102.7%)
$J_{4} (x)$	364 (−21.04%)	494 (−20.19%)	44 (120%)	178 (95.60%)	222 (100%)
$J_{5} (x)$	360 (−21.91%)	500 (−19.22%)	45 (125%)	178 (95.60%)	223 (100.9%)
$J_{6} (x)$	359 (−22.13%)	517 (−16.48%)	48 (140%)	178 (95.60%)	226 (103.6%)
$J_{7} (x)$	362 (−21.48%)	541 (−12.6%)	52 (160%)	175 (92.31%)	227 (104.5%)
$J_{8} (x)$	357 (−22.56%)	570 (−7.92%)	57 (185%)	156 (71.43%)	213 (91.89%)
$J_{9} (x)$	360 (−21.91%)	612 (−1.13%)	61 (205%)	112 (23.08%)	173 (55.86%)

As shown in Figure 12 and Table 8, a larger value of the weighted objective $J_{i} (x)$ does not necessarily indicate a globally superior design. For instance, although $J_{1} (x),$ which is strongly biased toward ${BW}_{2}$ , attains the largest weighted-sum value (0.365), the width of the first band gap decreases substantially from 20 to 8 Hz (–60%), leading to an imbalanced performance between the two band gaps. Moreover, while the lower band-edge frequency of the second band gap is significantly reduced ( $f_{2}^{low}$ : 619 → 469 Hz, −24.23%), the downshift of $f_{1}^{low}$ remains limited.

The extrema of each performance metric are summarized as follows. The minimum $f_{1}^{low}$ is achieved by $J_{8} (x)$ at 357 Hz (−22.56%), whereas the minimum $f_{2}^{low}$ occurs for $J_{1} (x)$ at 469 Hz (−24.23%). The maximum bandwidths $Δ f_{1}$ and $Δ f_{2}$ are obtained by $J_{9} (x)$ (61 Hz, 205%), and $J_{2} (x)$ (191 Hz, 109.89%), respectively, while the maximum total bandwidth $Δ f_{1} + Δ f_{2}$ is achieved by $J_{7} (x)$ (227 Hz, 104.5%). Considering the overall balance between low-frequency downshift and bandwidth expansion, $J_{5} (x)$ is identified as the most reasonable compromise solution. Although its total bandwidth ( $Δ f_{1} + Δ f_{2}$ $= 223 Hz$ ) is only 4 Hz smaller than the maximum value obtained by $J_{7} (x)$ , both lower band-edge frequencies shift downward simultaneously with a balanced reduction ( $f_{1}^{low} : 461 \to 360$ Hz, −21.91%; $f_{2}^{low} : 619 \to 500$ Hz, −19.22%). This balanced improvement demonstrates that the weighted-sum formulation can effectively capture trade-offs between competing objectives and identify designs that enhance low-frequency, broadband band-gap performance in a physically meaningful manner.

To quantify the computational efficiency of the proposed surrogate-assisted framework, the computational cost is summarized in Table 9. In the present study, a single full-order FEA-based dispersion analysis required approximately 640 s for one design candidate. Since the final Kriging surrogate model was constructed using 189 DOE samples, generating the corresponding high-fidelity training dataset required approximately 120,960 s (about 33.6 h). This initial cost is primarily associated with the construction of the surrogate model. As shown in Table 9, each surrogate-based HMA optimization required only about 236–333 s, with approximately $2 \times 10^{4}$ – $3 \times 10^{4}$ objective-function evaluations. If the same number of evaluations were carried out using the full-order FEA model, the computational cost would be prohibitively high (on the order of several thousand hours). Once the surrogate model was established, these repeated evaluations were performed efficiently at negligible cost compared to full-order simulations. These results clearly demonstrate the significant reduction in computational time achieved by the proposed framework and highlight its practical advantage for repeated design exploration under multiple weighting conditions in multi-objective optimization

Table 9.

Elapsed time, number of iterations, and function calls for surrogate-assisted HMA optimization of each weighted-sum objective function based on the final Kriging model (DOE size = 189).

HMA	Elapsed time (s)	Number of iterations (-)	Number of function calls (-)
$J_{1} (x)$	236	47	20,108
$J_{2} (x)$	272	54	23,209
$J_{3} (x)$	264	53	22,766
$J_{4} (x)$	301	61	26,310
$J_{5} (x)$	333	69	29,854
$J_{6} (x)$	286	58	24,981
$J_{7} (x)$	298	57	24,538
$J_{8} (x)$	288	53	22,766
$J_{9} (x)$	273	54	23,209

Figure 13(a) compares the dispersion relations of the initial design and the selected Pareto-optimal design $J_{5} (x)$ . For $J_{5} (x)$ , two distinct band gaps are clearly formed, with both lower band-edge frequencies shifted toward lower values relative to the initial design, accompanied by a noticeable increase in band-gap widths. Figure 13(b) presents the corresponding unit cell configurations of the initial and optimized designs, respectively, while Table 10 summarizes the numerical values of the design variables before and after optimization, together with their absolute changes. As reported in Table 10, the geometric modifications introduced by the optimization can be summarized as follows: (i) increases in the resonator mass heights $M^{A}$ and $M^{B}$ ; (ii) decreases in the resonator stiffness thicknesses $S_{h}^{A}$ and $S_{h}^{B}$ ; and (iii) increases in the in-plane offsets $S_{x}^{A}$ , $S_{y}^{A}$ , $S_{x}^{B}$ , and $S_{y}^{B}$ . This combined adjustment of the mass-related and stiffness-related parameters increases the effective mass while decreasing the effective stiffness of the resonators. Consequently, the local resonance frequencies are shifted downward, leading to lower band-edge frequencies and widened band gaps. These results provide clear physical insight into how the optimized geometric configuration of the M-LRM unit cell contributes to enhanced low-frequency and broadband vibration attenuation.

Figure 13.

Comparison of dispersion characteristics and unit cell configurations between the initial design and the selected Pareto-optimal design $J_{5} (x)$ : (a) dispersion curves and (b) unit cell configuration of the selected optimal design $J_{5} (x)$ (shaded regions highlight geometric modifications relative to the initial configuration).

Table 10.

Design variables of the selected Pareto-optimal design $J_{5} (x)$ : initial values, optimized values, and absolute changes with respect to the initial design.

Design variables	Initial value (mm)	Optimal value (mm)	Change from initial (mm)
$S_{x}^{A}$	1.5	1.9	0.4
$S_{y}^{A}$	0	1.386	1.386
$S_{h}^{A}$	2	1.818	−0.182
$M^{A}$	6	6.696	0.696
$S_{x}^{B}$	2.5	3.433	0.933
$S_{y}^{B}$	0	0.933	0.933
$S_{h}^{B}$	2	1.728	−0.272
$M^{B}$	7	7.311	0.311

In addition, as described in Section 3.2, geometric constraints were incorporated to prevent unrealistic or difficult-to-manufacture shapes, and the optimized geometry was therefore obtained with manufacturability considered at the conceptual design stage. However, a detailed manufacturability assessment based on specific fabrication processes, as well as experimental verification of structural integrity, are beyond the scope of the present study and will be addressed in future work.

Frequency response analysis of the selected optimal design

To validate the band gaps predicted by the dispersion analysis for the selected Pareto-optimal design, frequency response functions (FRFs) were computed for a finite panel consisting of a 5 × 5 array of unit cells. All modeling assumptions, boundary conditions, and loading configurations were identical to those described in Section 2.3.

As shown in Figure 14, the FRF of the optimized design exhibits pronounced attenuation of the displacement response over frequency ranges that closely coincide with the band gaps identified in the dispersion curves of Figure 13. In particular, within both the first (360–405 Hz) and second (500–678 Hz) band-gap intervals, the response magnitude of the optimal design is consistently lower than that of the initial design. These results demonstrate good agreement between the dispersion-based band-gap predictions and the vibration response of the finite structure. The observed reduction in displacement response confirms that the band-gap-induced vibration suppression mechanism remains effective in the finite panel configuration, thereby validating the practical applicability of the optimized M-LRM design.

Figure 14.

Displacement FRF for the initial design and the selected Pareto-optimal design $J_{5} (x)$ . Gray shaded regions denote the first and second band-gap frequency ranges identified from the dispersion analysis of $J_{5} (x)$ .

Conclusion

This study presented a surrogate-model-based optimal design framework for shifting the band gaps of multiple locally resonant metamaterials (M-LRMs) toward lower frequencies while simultaneously enlarging their band-gap widths. Band-gap characteristics were systematically identified through dispersion analysis of the unit cell, and optimal designs were obtained for different trade-offs using a weighted-sum formulation combined with a Kriging surrogate model and a Hybrid Metaheuristic Algorithm (HMA).

The performance of the optimized designs was quantitatively validated. For a selected Pareto-optimal design, both band gaps exhibited clear downward shifts in their lower band-edge frequencies together with substantial increases in bandwidth compared with the initial configuration. Furthermore, frequency response function (FRF) analysis of a finite panel demonstrated pronounced vibration attenuation over frequency ranges consistent with the predicted band gaps, confirming the effectiveness of the dispersion-based design in a practical finite-structure setting.

Overall, the proposed M-LRM design framework provides an effective and practical strategy for achieving low-frequency, broadband vibration reduction in lightweight panel structures. The methodology is general and can be extended to other metamaterial configurations and structural applications requiring efficient vibration control.

Footnotes

Handling Editor: Mahdi Mofidian

ORCID iD

Junghwan Kook

Ethical considerations

This article does not contain any studies with human participants or animals performed by any of the authors.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Author contributions

Kwangjin Kim: Conceptualization, Methodology, Investigation, Validation, Writing—Original Draft. Junghwan Kook: Conceptualization, Resources, Writing—Review & Editing, Supervision.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (RS-2025-25415734). This work was partly supported by the KETEP grant funded by the Ministry of Trade, Industry & Energy (MOTIE), Republic of Korea (No. 20241K00000010; Robot Utilization Production Support Center for SMR).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

All data generated or analyzed during this study are included in this article.

References

Lin

Jiang

, et al. Use of high strength steel sheet for lightweight and crashworthy car body. Mater Des 2003; 24: 177–182. https://doi.org/10.1016/s0261-3069(03)00021-9

Zhang

. Advanced lightweight materials for automobiles: a review. Mater Des 2022; 221: 110994. https://doi.org/10.1016/j.matdes.2022.110994

Chen

Feng

, et al. Structural optimization design of BIW using NSGA-III and entropy weighted TOPSIS methods. Adv Mech Eng 2023; 15: 16878132231220351. https://doi.org/10.1177/16878132231220351

Candela

Sandrini

Gadola

, et al. Lightweighting in the automotive industry as a measure for energy efficiency: review of the main materials and methods. Heliyon 2024; 10: e29728. https://doi.org/10.1016/j.heliyon.2024.e29728

Kumar

Walsh

Krylov

. Structural–acoustic behaviour of automotive-type panels with dome-shaped indentations. Appl Acoust 2013; 74: 897–908. https://doi.org/10.1016/j.apacoust.2012.11.017

Gur

Pan

Wagner

. Sound package development for lightweight vehicle design using statistical energy analysis (SEA). Report no. 0148-7191. SAE Technical Paper, 2015.

Koukounian

Mechefske

. Computational modelling and experimental verification of the vibro-acoustic behavior of aircraft fuselage sections. Appl Acoust 2018; 132: 8–18. https://doi.org/10.1016/j.apacoust.2017.11.004

Kim

H-G

Nerse

Wang

. Topography optimization of an enclosure panel for low-frequency noise and vibration reduction using the equivalent radiated power approach. Mater Des 2019; 183: 108125. https://doi.org/10.1016/j.matdes.2019.108125

Rao

. Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes. J Sound Vib 2003; 262: 457–474. https://doi.org/10.1016/s0022-460x(03)00106-8

10.

Zhou

Shao

, et al. Research and applications of viscoelastic vibration damping materials: a review. Compos Struct 2016; 136: 460–480. https://doi.org/10.1016/j.compstruct.2015.10.014

11.

Yoon

. Structural topology optimization for frequency response problem using model reduction schemes. Comput Methods Appl Mech Eng 2010; 199: 1744–1763. https://doi.org/10.1016/j.cma.2010.02.002

12.

Bein

Bös

Mayer

, et al. Advanced materials and technologies for reducing noise, vibration and harshness (NVH) in automobiles. In: Rowe J (ed.) Advanced materials in automotive engineering. Elsevier, 2012, pp.254–298.

13.

Shin

Bolton

. The identification of minimum-weight sound packages. Noise Control Eng J 2018; 66: 523–540. https://doi.org/10.3397/1/376644

14.

Liu

Zhang

Mao

, et al. Locally resonant sonic materials. Science 2000; 289: 1734–1736. https://doi.org/10.1126/science.289.5485.1734

15.

Liu

Chan

Sheng

. Analytic model of phononic crystals with local resonances. Phys Rev B 2005; 71: 014103. https://doi.org/10.1103/physrevb.71.014103

16.

Gao

Zhang

Deng

, et al. Acoustic metamaterials for noise reduction: a review. Adv Mater Technol 2022; 7: 2100698. https://doi.org/10.1002/admt.202100698

17.

Nakayama

. Acoustic metamaterials based on polymer sheets: from material design to applications as sound insulators and vibration dampers. Polym J 2024; 56: 71–77. https://doi.org/10.1038/s41428-023-00842-0

18.

Xiao

Wen

Wang

, et al. Theoretical and experimental study of locally resonant and Bragg band gaps in flexural beams carrying periodic arrays of beam-like resonators. J Vib Acoust 2013; 135: 041006. https://doi.org/10.1115/1.4024214

19.

Chen

Wang

. Band gaps in the low-frequency range based on the two-dimensional phononic crystal plates composed of rubber matrix with periodic steel stubs. Physica B Condens Matter 2013; 416: 12–16. https://doi.org/10.1016/j.physb.2013.02.011

20.

Chen

J-S

Chen

T-Y

Chang

Y-C

. Metastructures with double-spiral resonators for low-frequency flexural wave attenuation. J Appl Phys 2021; 130: 130. https://doi.org/10.1063/5.0053579

21.

Cai

Hui Wu

, et al. Realizing polarization band gaps and fluid-like elasticity by thin-plate elastic metamaterials. Compos Struct 2021; 262: 113351. https://doi.org/10.1016/j.compstruct.2020.113351

22.

Jung

Goo

Wang

. Investigation of flexural wave band gaps in a locally resonant metamaterial with plate-like resonators. Wave Motion 2020; 93: 102492. https://doi.org/10.1016/j.wavemoti.2019.102492

23.

Jung

Kim

H-G

Goo

, et al. Realisation of a locally resonant metamaterial on the automobile panel structure to reduce noise radiation. Mech Syst Signal Process 2019; 122: 206–231. https://doi.org/10.1016/j.ymssp.2018.11.050

24.

Pires

Sangiuliano

Denayer

, et al. The use of locally resonant metamaterials to reduce flow-induced noise and vibration. J Sound Vib 2022; 535: 117106. https://doi.org/10.1016/j.jsv.2022.117106

25.

Hyun

Jung

Park

, et al. Simultaneous low-frequency vibration isolation and energy harvesting via attachable metamaterials. Nano Converg 2024; 11: 38. https://doi.org/10.1186/s40580-024-00445-2

26.

Manushyna

Hülsebrock

Kuisl

, et al. Application of vibroacoustic metamaterials for structural vibration reduction in space structures. Mech Res Commun 2023; 129: 104090. https://doi.org/10.1016/j.mechrescom.2023.104090

27.

Jung

Goo

Kook

. Design of a local resonator using topology optimization to tailor bandgaps in plate structures. Mater Des 2020; 191: 108627. https://doi.org/10.1016/j.matdes.2020.108627

28.

Kook

. Design and validation of multiple locally resonant metamaterial. J Korean Soc Nondestruct Test 2024; 44: 1–7. https://doi.org/10.7779/jksnt.2024.44.1.1

29.

Miao

Yin

Yang

, et al. Design of multi-bandgap metamaterial plate based on composite cylindrical resonators. Mater Des 2025; 250: 113570. https://doi.org/10.1016/j.matdes.2024.113570

30.

Claeys

Vergote

Sas

, et al. On the potential of tuned resonators to obtain low-frequency vibrational stop bands in periodic panels. J Sound Vib 2013; 332: 1418–1436. https://doi.org/10.1016/j.jsv.2012.09.047

31.

Nunez-Solano

Kook

Puyana-Romero

. Experimental study of coupled acoustic-mechanical bandgaps. J Mech Sci Technol 2023; 37: 4973–4980. https://doi.org/10.1007/s12206-023-0904-9

32.

Elabbasi

. Modeling phononic band gap materials and structures, https://www.comsol.com/blogs/modeling-phononic-band-gap-materials-and-structures (2016, accessed 14 October 2025).

33.

Maurin

Claeys

Deckers

, et al. Probability that a band-gap extremum is located on the irreducible Brillouin-zone contour for the 17 different plane crystallographic lattices. Int J Solids Struct 2018; 135: 26–36. https://doi.org/10.1016/j.ijsolstr.2017.11.006

34.

PIDOTECH. PIAnO signature, https://pidotech.com/ (2022, accessed 3 February 2025).

35.

Santner

Williams

Notz

. The design and analysis of computer experiments. Springer, 2003.

36.

Peng

Zhang

Shi

, et al. Low-frequency sound insulation optimisation design of membrane-type acoustic metamaterials based on Kriging surrogate model. Mater Des 2023; 225: 111491. https://doi.org/10.1016/j.matdes.2022.111491

37.

Wang

Zhou

, et al. Kriging surrogate model for optimizing outlet temperature distribution in low-emission combustors without dilution holes. AIP Adv 2024; 14: 055310. https://doi.org/10.1063/5.0198258

38.

Santos

Costa

CBB

Caballero

, et al. Multi-objective simulation–optimization via kriging surrogate models applied to natural gas liquefaction process design. Energy 2023; 262: 125271. https://doi.org/10.1016/j.energy.2022.125271

39.

Fuhg

Fau

Nackenhorst

. State-of-the-art and comparative review of adaptive sampling methods for kriging. Arch Comput Methods Eng 2021; 28: 2689–2747. https://doi.org/10.1007/s11831-020-09474-6

40.

Pan

. Adaptive optimization methodology based on Kriging modeling and a trust region method. Chin J Aeronaut 2019; 32: 281–295. https://doi.org/10.1016/j.cja.2018.11.012

41.

Wang

Duan

Gong

, et al. An evaluation of adaptive surrogate modeling based optimization with two benchmark problems. Environ Model Softw 2014; 60: 167–179. https://doi.org/10.1016/j.envsoft.2014.05.026

42.

Zhang

Ding

Jia

. A hybrid optimization algorithm based on cuckoo search and differential evolution for solving constrained engineering problems. Eng Appl Artif Intell 2019; 85: 254–268. https://doi.org/10.1016/j.engappai.2019.06.017

43.

Park

. An efficient hybrid metaheuristic algorithm for solving constrained global optimization problems. PhD Thesis, Hanyang University, 2016.

44.

Gong

Cai

Liang

. Engineering optimization by means of an improved constrained differential evolution. Comput Methods Appl Mech Eng 2014; 268: 884–904. https://doi.org/10.1016/j.cma.2013.10.019