## Abstract

In this article, we present a design methodology for resonant structures exhibiting particular dynamic responses by combining an eigenfrequency matching approach and a harmonic analysis-informed eigenmode identification strategy. This systematic design methodology, based on topology optimization, introduces a novel computationally efficient approach for 3D dynamic problems requiring antiresonances at specific target frequencies subject to specific harmonic loads. The optimization’s objective function minimizes the error between target antiresonance frequencies and the actual structure’s antiresonance eigenfrequencies, while the harmonic analysis-informed identification strategy compares harmonic displacement responses against eigenvectors using a modal assurance criterion, therefore ensuring an accurate recognition and selection of appropriate antiresonance eigenmodes used during the optimization process. At the same time, this method effectively prevents well-known problems in topology optimization of eigenfrequencies such as localized eigenmodes in low-density regions, eigenmodes switching order, and repeated eigenfrequencies. Additionally, our proposed localized eigenmode identification approach completely removes the spurious eigenmodes from the optimization problem by analyzing the eigenvectors’ response in low-density regions compared to high-density regions. The topology optimization problem is formulated with a density-based parametrization and solved with a gradient-based sequential linear programming method, including material interpolation models and topological filters. Two case studies demonstrate that the proposed design methodology successfully generates antiresonances at the desired target frequency subject to different harmonic loads, design domain dimensions, mesh discretization, or material properties.

## 1 Introduction

Devising effective means to control wave propagation is of paramount importance across length scales, from designing electronic devices to vibration isolation of critical structures. In search of effective means, an approach inspired by the concept of locally resonant metamaterials uses sub-wavelength resonators to create frequency bandgaps [1], i.e., a frequency range where waves do not propagate. A subclass of these metamaterials composed of surface-mounted or surface-embedded resonators are called locally resonant metasurfaces (LRMs). These LRMs are fundamentally different from phononic crystals whose working mechanisms, explained by Bragg scattering [2], depend upon the periodicity of their unit cells. Instead, LRMs rely on interactions between the propagating wave and the local resonator’s response, exhibiting scattering and nonconventional dispersion properties that lead to frequency bandgaps [3]. When these LRMs are used to control elastic guided wave propagation, e.g., Lamb or Rayleigh waves, they are called elastodynamic LRMs (ELRMs).

Since the introduction of ELRMs, their local resonator designs have been based on simple geometries, i.e., rods, holes, spheres, beams, or simply mass-spring systems. These designs often rely on the parametric tuning of geometrical features empirically [4–9] or experimentally [10–14] until the desired bandgap is achieved. Optimization and machine learning methods have been used to tailor dispersion curves [15–21] or effective material properties [22–25] mainly for 2D phononic crystals assuming periodicity of unit cells; however, a rational design methodology to tailor the local resonator’s response to particular wave excitations is missing. Our investigation [26–28] demonstrates that the local resonator’s *antiresonances* are the key factor in generating frequency bandgaps for ELRMs. Thus, the motivation of our work is how to design resonant structures by tailoring their antiresonances. This article establishes a generalized systematic design methodology to conceive resonant structures based on antiresonance matching. One of the methodology’s applications is the design of local resonators for ELRMs; however, it is applicable to any design problem requiring antiresonance manipulation.

The topology optimization (TO) method [29], which is of particular interest to this study, is a computational tool to design complex layouts in numerous applications [30–32], and it has been proven to be a robust design tool for dynamic problems [33,34]. Previous works on TO have attempted to manipulate antiresonances with harmonic-based formulations such as the minimum dynamic compliance [35–37], solving mostly 2D problems with only a few design variables and surrogate formulations to reduce computational load. 3D problems have significantly more design variables, resulting in complex and expensive computational tasks [38]. Moreover, the dynamic compliance could be undefined around an antiresonance frequency [36,39], making the minimum dynamic compliance problem not suitable for antiresonance manipulation. Thus, a harmonic-based TO is not feasible. Instead, we propose an approach that relies on manipulating eigenfrequencies, which are significantly less expensive to compute, in conjunction with a harmonic-informed selection of antiresonance eigenmodes to be used within the optimization process.

The eigenfrequency problem is well known in TO. It has been commonly used from three different approaches [38]: maximize the structure’s eigenfrequencies [40], generate a gap in-between eigenfrequencies [41], and match the structure’s eigenfrequencies with target frequencies [42]. The last approach, which is most relevant to this article, has been mostly used to match resonance eigenfrequencies without considering how dynamic forces interact with the structure, an important consideration for the antiresonance matching problem. This is because matching eigenfrequencies for resonances or antiresonances are fundamentally different problems. Resonances occur when a dynamic force excites the structure’s natural frequencies of resonance, namely, *resonance eigenfrequencies*. On the other hand, antiresonances occur when the dynamic force is out of phase with respect to the structure’s response, resulting in zero or near-to-zero displacement(s) at a particular location(s); the dynamic load cancels out with an equal but opposite reaction force, equivalent to constraining the displacement(s) at that point(s) [43]. The antiresonance frequencies of the dynamically loaded structure correlate with *specific* eigenfrequencies of the displacement-constrained structure [44]. Here, we refer to those displacement-constrained eigenfrequencies as the *antiresonance eigenfrequencies*. Contrary to resonance eigenmodes, several antiresonance eigenmodes are not excited by particular dynamic loads, making the selection of these eigenmodes a critical component of the antiresonance eigenfrequency matching approach. Thus, the commonly used eigenfrequency matching approaches do not work for the antiresonance matching problem as it needs additional information to appropriately select antiresonance eigenmodes considering specific dynamic loads. This article presents a new harmonic-informed eigenfrequency matching approach that properly identifies the particular antiresonance eigenmodes to be optimized.

This article presents a twofold contribution: (i) a computationally feasible design methodology for 3D structures exhibiting antiresonances at desired target frequencies using a density-based TO combining eigenfrequency matching and harmonic analyses, and (ii) a new analysis procedure to identify specific antiresonance eigenmodes using a harmonic-informed identification scheme, effectively preventing several problems reported in the literature for the eigenfrequency optimization problem. The remainder of this article is organized into three sections: Sec. 2 introduces the design methodology, Sec. 3 presents two case studies and their corresponding analysis, and Sec. 4 summarizes our findings and discusses future work.

## 2 Design Methodology

This section presents a design methodology seeking to generalize a topology optimization problem in which a structure’s dynamic response is tailored by manipulating antiresonance eigenfrequencies with a harmonic-informed selection of eigenmodes.

### 2.1 Optimization Problem.

*f*

_{T}and the structure’s antiresonance eigenfrequency

*f*

_{A}selected by a harmonic-informed analysis routine. Thus, the optimization problem formulation is expressed as follows:

*ρ*are the pseudo-densities, namely, design variables. Equation (2) presents the optimization constraints, where

*ρ*

_{e}is the pseudo-density of element

*e*, with

*N*

_{e}the total number of finite elements. The minimum pseudo-density, set to

*ρ*

_{min}= 1 × 10

^{−6}, prevents numerical problems.

*ρ*

_{e}

*V*

_{e}is the effective volume of element

*e*, and the maximum and minimum volume constraints are, respectively,

*V*

_{max}and

*V*

_{min}. [

*K*] and [

*M*] are, respectively, the global stiffness and mass matrices, and Φ

_{A}is the eigenvector corresponding to the eigenvalue

*λ*

_{A}= (2

*πf*

_{A})

^{2}.

### 2.2 Optimization Solver.

*LB*

_{j}=

*ρ*

_{j}− 0.01 and

*UB*

_{j}=

*ρ*

_{j}+ 0.01. The objective function in Eq. (1) is linearized using first-order Taylor series as follows:

### 2.3 Harmonic-Informed Selection of Eigenmodes.

Identifying appropriate antiresonance eigenfrequencies to be matched with the target frequency is critical to ensure feasible solutions. By analyzing the structure’s frequency response function (FRF) at each iteration, suitable antiresonances are identified. The eigenmodes corresponding to these antiresonances are selected as appropriate modes to use during the optimization process. The harmonic analysis-informed identification approach consists of the following steps:

- Compute antiresonance eigenvalues {
*λ*} and eigenvectors [Φ]:$[K][\Phi ]={\lambda}[M][\Phi ]$ Delete rigid body eigenmodes below a threshold of 1 Hz

Remove localized eigenmodes using the routine detailed in Appendix A.

- Compute the structure’s harmonic response to the harmonic force
*F*_{Harmonic}$[M]{x\xa8}+[C]{x\u02d9}+[K]{x}={FHarmonic}$ - Extract the averaged FRF at the points of interest:$averaged_FRF=mean(FRF)$
- Find local antiresonance peaks using the matlab function findpeaks:$[Amplitude,idx,Width,Prominence]=findpeaks(1/averaged_FRF)$
- Evaluate local peak frequency distance to target frequency
*g*_{q}:$Proximity=abs(Target_frequency\u2212Peak_frequency)$ - Evaluate metrics for each local antiresonance peak identified, i.e., antiresonance bandwidth, amplitude, prominence, and proximity to the target frequency
*g*_{q}. Define a combined antiresonance evaluation factor as follows:$factor=Amplitude+Prominence+Width\u2212Proximity$ Select the local antiresonance peak with the highest evaluation factor as the identified antiresonance frequency

*A*_{f}.Extract the harmonic displacement response

*x*_{H}at the identified frequency*A*_{f}.- Evaluate the modal assurance criterion (MAC) [47] between the harmonic displacement response
*x*_{H}and all the antiresonance eigenvectors $\Phi A$:$MAC=|\varphi ATxH|2(\varphi AT\varphi A)(\varphi ATxH)$ Select the eigenmode with the highest MAC coefficient as the identified antiresonance eigenmode; its corresponding eigenvalue

*λ*_{A}= (2*πf*_{A})^{2}is used during the optimization in Eq. (1).

The reader is referred to Appendix A for a step-by-step explanation.

### 2.4 Topology Optimization Parametrization.

*E*

_{e}is the interpolated material property for element

*e*,

*E*

_{0}is the original material property, and

*p*is the penalization factor. Note that

*E*may represent stiffness, mass, or damping. Typical penalization factors for stiffness (

*p*

_{K}) and mass (

*p*

_{M}) are, respectively,

*p*

_{K}= 3 and

*p*

_{M}= 1 [40,51–54]. However, this combination of penalization factors promotes the generation of spurious localized eigenmodes in low-pseudo-density regions as it is shown in Appendix B. Thus, this design methodology uses the same penalization factors, i.e.,

*p*

_{K}=

*p*

_{M}, regardless of the interpolation model.

### 2.5 Topology Optimization Problems in Eigenfrequency Approaches.

In density-based TO, it is necessary to address problems such as checkerboard solutions, mesh dependency, or nonfeasible solutions [31]. Moreover, the eigenfrequency matching approach presents other problems such as eigenmodes switching order or repeated eigenfrequencies [42], localized eigenmodes in low-density regions [40], and disconnected members [55]. The following solutions are proposed in this design methodology.

#### 2.5.1 Eigenmode Switching Order or Repeated Eigenfrequencies.

Ma et al. [42] introduced the mean-eigenvalue formulation to overcome these problems, although it only works for multiple eigenfrequencies. Kim and Kim [56] introduced a tracking mode scheme that allows the tracking of single eigenmodes, although the tracked eigenmodes must be prescribed in advance. In this design methodology, these problems are overcome by matching eigenvectors with a harmonic response at an antiresonance frequency, as it is explained in Sec. 2.3 and detailed in Appendix A.

#### 2.5.2 Localized Eigenmodes.

This problem was originally addressed by Neves et al. [57] by removing the low-density elements from the stiffness matrix; this approach modifies the design domain at each iteration, which is not feasible in density-based TO. Several approaches appear in the literature. Pedersen [40] modified the penalization stiffness-to-mass ratio (*p*_{K}/*p*_{M}) in low-density regions. Tcherniak [52] and Du and Olhoff [51] set the mass matrix penalization to zero or near to zero in low-density regions, therefore modifying the stiffness-to-mass ratio. Similarly, Qiao et al. [54], Rong et al. [58], and Zhang and Kang [59] used different penalization schemes for the stiffness and mass matrices combining SIMP and RAMP interpolation models. All these articles have reported improvements in the localized eigenmodes problem by modifying the stiffness-to-mass ratio, mostly for the maximizing eigenfrequencies problem [40]; however, for the antiresonance matching problem, we have not seen the same results, as explained in Appendix B. This design methodology considers a constant ratio *p*_{K}/*p*_{M} = 1, i.e., all interpolation models and penalization factors are applied equally for stiffness, mass, and damping. On another note, Li et al. [55] pointed out another kind of localized eigenmodes, not for low-density elements but for disconnected floating members. Li et al. proposed a recognition technique to remove these localized eigenmodes for a level-set TO. Here, we propose a new localized eigenmode recognition and elimination routine for density-based TO that removes both the low-density and the disconnected members localized eigenmodes, as mentioned in Sec. 2.3 and detailed in Appendix A.

#### 2.5.3 Disconnected Members.

*N*

_{s}is the set of elements within a filter radius

*R*

_{f}, such that

*N*

_{s}= {

*i*|

*R*

_{ie}≤

*R*

_{f}} with

*R*

_{ie}=

*c*

_{i}−

*c*

_{e}being the distance between the element centers

*e*and

*i*. The weights are

*w*

_{i}=

*R*

_{f}−

*R*

_{ie}.

*v*

_{i}and

*ρ*

_{i}are the volume and pseudo-density of element

*i*, respectively.

*β*controls the strength of the filter. In this design methodology, a continuation scheme is used, progressively increasing

*β*to remove intermediate-density material, therefore promoting structural definition. The

*η*parameter controls the volume during filtering. We preserve the volume after filtering by optimizing the

*η*value with a bisection method at each iteration.

### 2.6 Optimization Program.

An optimization program to design resonating structures based on the proposed design methodology is implemented. This program, composed of multiple modules, integrates the software matlab and abaqus through python communication scripts [61], as shown in the flowchart in Fig. 1. The optimization program starts with a finite element model data generation, including stiffness, mass, and damping matrices, mesh data with its nodal coordinates, element connectivity, loads, and boundary conditions. After this initialization, the optimization loop starts from the filters, followed by the finite element analysis (FEA) considering matrices modified by filters and material interpolation models. The harmonic-informed routine from Sec. 2.3 analyses the topology’s frequency response at each iteration to select the most suitable eigenmode for computing the objective function (Eq. (1)) and sensitivities (Eq. (4)).

The proposed optimization workflow has been generalized to design any kind of structure by matching resonance or antiresonance eigenfrequencies. It can accommodate different material interpolation models, maximum and minimum volume constraints, solid or void nondesign domains, numerical or analytical sensitivities, and multiple filtering schemes such as density, Heaviside, sensitivity, or move limits filters. Moreover, the integration with commercial fea software allows analysis of structures with complex geometries, multiphysics problems, and a variety of boundary conditions.

## 3 Case Studies

This section illustrates the design of resonant structures using the methodology proposed in Sec. 2. Two case studies present resonator designs for different harmonic loads, frequency ranges, design domain dimensions, mesh discretization, and material properties. Figure 2 shows design domains with their mesh discretization for both case studies. Note that a symmetry condition on the *xz*-plane is used to reduce computation time. Table 1 lists the initial parameters for both case studies.

Material properties, case study 1 | E = 3.2516 GPaρ = 1222.2 kg/m^{3}v = 0.33 |

Material properties, case study 2 | E = 69 GPaρ = 2730 kg/m^{3}v = 0.33 |

Target antiresonance frequency f_{T} | 30 kHz (case study 1) 50 kHz (case study 2) |

Starting point ρ | Homogenous pseudo-density at ρ_{e} = 0.5Solid nondesign 4 mm × 2 mm bottom base. Void nondesign all other first-layer elements |

Maximum allowed volume | $Vmax=90%$ |

Minimum allowed volume | $Vmin=10%$ |

Interpolation model | SIMP model |

Penalization factors | p_{K} = p_{M} = 1 |

Density filter radius R_{f} | 3 mm |

Heaviside filter β starting value | β = 20 |

β continuation scheme | β = β + 1; at every iteration |

Material properties, case study 1 | E = 3.2516 GPaρ = 1222.2 kg/m^{3}v = 0.33 |

Material properties, case study 2 | E = 69 GPaρ = 2730 kg/m^{3}v = 0.33 |

Target antiresonance frequency f_{T} | 30 kHz (case study 1) 50 kHz (case study 2) |

Starting point ρ | Homogenous pseudo-density at ρ_{e} = 0.5Solid nondesign 4 mm × 2 mm bottom base. Void nondesign all other first-layer elements |

Maximum allowed volume | $Vmax=90%$ |

Minimum allowed volume | $Vmin=10%$ |

Interpolation model | SIMP model |

Penalization factors | p_{K} = p_{M} = 1 |

Density filter radius R_{f} | 3 mm |

Heaviside filter β starting value | β = 20 |

β continuation scheme | β = β + 1; at every iteration |

The starting point considers all design variables to have the same pseudo-density, except for the bottom layer of elements, which considers a fully solid nondesign bottom base at which the harmonic load must be applied, surrounded by void nondesign elements. This is done to prevent numerical errors in the frequency response when harmonic loads are applied to void or near-to-void elements. The volume constraints are in place solely to prevent numerical instabilities in extreme cases, i.e., uncontrolled allocation of solid material or uncontrolled removal of material. However, the volume constraints do not play a major role in the proposed design methodology, as is usually done in topology optimization for lightweight structures. For a comparative study on how different volume constraints influence the solution, refer to Appendix C.

### 3.1 Case Study 1.

This case seeks to design a resonator that exhibits an antiresonance at 30 kHz subject to normal (out-of-plane) harmonic loads, as shown in Fig. 2(a). Figure 3 shows the convergence history, Fig. 4 shows the resultant optimized topology, and Fig. 5 shows the topology’s frequency response.

Figure 3 shows the objective function converging after iteration 52 with a value of 0.63. Several spikes are seen in the objective function before convergence. This phenomenon occurs when the harmonic-informed routine identifies antiresonance with a high evaluation factor but low proximity to the target frequency, as described in Sec. 2.3; therefore, the objective function increases rapidly. The selection of antiresonance closer to the target frequency is quickly recovered in the following iteration. This behavior is considered normal because the optimization explores the solution space for optimal solutions, especially during the first iterations, jumping to other potential solutions, and returning to the lowest value. Once the objective function has converged, the recognition of antiresonances is stable around the target frequency.

An antiresonance function has been computed by replacing the eigenfrequency *f*_{g} in Eq. (1) with the antiresonance *A*_{f} identified by the harmonic-informed routine of Sec. 2.3. The correlation between the objective function and the antiresonance function in Fig. 3 reveals whether the selected eigenfrequency accurately represents an antiresonance when the topology is subject to a specific harmonic loading condition. Figure 3 shows good agreement between the objective function and the antiresonance function, therefore validating our design methodology to tailor antiresonances using eigenfrequency matching with a harmonic-informed selection of eigenmodes.

Figure 4 presents both the optimized and postprocessed topologies, and Fig. 5 presents their corresponding frequency responses. The volume decreased from 50%, the starting point, to 15.27% at the final iteration, as shown in Fig. 4(a), where most of the solid material has been removed from the surroundings of a central solid body. After postprocessing, the resultant topology resembles a rounded Balloon-like shape. The FRFs presented in Fig. 5 evidence good agreement between the raw and postprocessed dynamic responses to harmonic loads, and most importantly, both FRFs evidence an antiresonance at 30 kHz; the design objective has been achieved successfully, therefore demonstrating our design methodology works for the case study.

### 3.2 Case Study 2.

This case seeks to design a resonator that exhibits an antiresonance at 50 kHz subject to shear (in-plane) harmonic loads, as shown in Fig. 2(b). Figure 6 shows the convergence history, Fig. 7 shows the optimized topology, and Fig. 8 shows the topology’s frequency response.

Figure 6 shows the objective function converging after iteration 50 with a value of 0.01. The antiresonance function, as explained in Sec. 2.1, was computed to evaluate the correlation between eigenfrequencies and antiresonances. Similar to case study 1, Fig. 6 reveals that the eigenfrequencies used to compute the objective function and the identified antiresonances from the FRFs have similar values, validating the harmonic-informed routine effectiveness to select eigenfrequencies by analyzing frequency responses, regardless of different harmonic loads, material properties, dimensions, mesh discretization, or target frequency. For this case study, the sudden peaks in the objective function shown in Fig. 3 are not seen. This objective function behavior shows a stable optimization process where the harmonic-informed selection of eigenmodes routine did not find any other antiresonances in higher frequency ranges during the optimization.

Figure 7 presents both the optimized and postprocessed topologies, and Fig. 8 presents their corresponding frequency responses. The volume decreased from 50%, the starting point, to 19.21% at the final iteration, as shown in Fig. 7(a), where most of the solid material has been removed from the surroundings of a central solid body. After postprocessing, the resultant topology resembles an asteroid-like shape. Figure 8 presents the FRFs for both topologies, and note that the raw topology’s FRF was computed from 10 to 150 kHz, while the postprocessed topology’s FRF was obtained up to 200 kHz to observe the closest resonance peak. In this case study, there are some differences between the raw and postprocessed topologies’ FRFs. Both topologies have antiresonances around the target frequency, but they do not align as perfectly as shown in Fig. 5. Those differences could be a consequence of using different types of finite elements in a structure subject mostly to shear stress; the raw topology uses hexahedral elements, while the postprocessed topology uses tetrahedral elements. Nonetheless, Fig. 8 shows acceptable agreement between the FRFs, exhibiting an antiresonance around 50 kHz, achieving the design objective, and validating our design methodology.

### 3.3 Locally Resonant Metasurfaces.

The design methodology proposed in Sec. 2 can be used to design any structure requiring specific antiresonances; one application is the design of topology-optimized resonators as demonstrated in Secs. 2.1 and 2.2. These resonators can be arranged to compose ELRMs, as it was introduced in Sec. 1. Figure 9 shows a conceptual design of an ELRMs composed of 56 topology-optimized resonators from Sec. 2.1 exemplifying a metasurface that reduces the propagation of surface waves at the resonator’s frequency of antiresonance when mounted on top of a surface.

## 4 Conclusions

A systematic design methodology to realize resonant structures exhibiting antiresonances at desired frequencies subject to specific harmonic loads was proposed as a structural optimization problem using density-based topology optimization. The key advancement to the well-known eigenfrequency matching approach [42] is the new harmonic-informed identification routine that ensures recognition of antiresonance eigenmodes while preventing common problems for this kind of eigenfrequency topology optimization. In other words, the topology optimization problem presented in this article is fundamentally different from other eigenfrequency matching approaches because it not only matches eigenfrequencies but also analyses frequency responses to select potential eigenmodes to become antiresonances using the identification routine presented in Sec. 2.3. Even though this approach is computationally more expensive than a simple eigenfrequency matching optimization, it is still faster to compute than a harmonic-based or full dynamic approach, making the design methodology feasible for solving large 3D dynamic problems. Moreover, well-known problems reported in the literature such as localized eigenmodes, repeated eigenfrequencies, and eigenmodes order switching [42,55,62] are simple to overcome using the harmonic-informed routine because appropriate eigenmodes are selected for use during the optimization. Additionally, the combination of a density filter with a Heaviside filter contributes to realizing well-defined structures whose dynamic response after postprocessing remains coherent with the original optimized topologies. To illustrate these contributions, we have demonstrated how the proposed design methodology provides a robust tool to realize resonant structures for different initial conditions, frequency ranges, material properties, or design domain definitions.

Interesting applications can be developed with the proposed design methodology, for example, an antiresonance at a given frequency is known to generate frequency bandgaps for locally resonant metasurfaces [26]; however, assessing a metasurface performance composed of topology-optimized resonators, such as the ones presented in this article, is yet to be studied. Beyond the design of acoustic or elastodynamic metasurfaces for wave propagation control at different length scales, the proposed design methodology has broader potential applications, including designing vibration control devices, sensors and actuators, energy harvesting devices, seismic shields, musical instruments, or any other application where the structure’s frequency response plays an important role. Moreover, the design methodology has been generalized to accommodate different design objectives, such as matching resonances and antiresonances simultaneously or introducing eigenfrequency gaps for broadband frequency applications.

## Acknowledgment

Any opinions, findings, conclusions, or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Computations for this research were performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer.

## Funding Data

• The National Science Foundation (Grant No. 1934527).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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

## Appendix A: Harmonic-Informed Selection of Antiresonance Eigenmodes

One of the most important components of the design methodology presented in this article is to appropriately identify and select antiresonance eigenmodes to compute the objective function and sensitivities. Section 2.3 introduces the harmonic-informed selection routine step by step; here, we present a detailed explanation of each step.

*Step 1*. Compute antiresonance eigenvalues {*λ*} and eigenvectors [Φ].

*λ*} and eigenvectors [Φ]:

*Step 2*. Delete rigid body eigenmodes below a threshold of 1 Hz.

*Step 3*. Remove localized eigenmodes.

An analysis routine that identifies and removes localized eigenmodes has been developed. Since the identified localized eigenmodes are fully removed from the optimization, the problem is effectively prevented. This routine removes localized eigenmodes following these steps:

Apply material interpolation models (SIMP or RAMP) and filters to the pseudo-density variable

*ρ*from Eqs. (6)–(9).Find solid or near-to-solid elements such that the pseudo-density is higher than 70%, i.e.,

*ρ*_{hard}=*ρ*_{e}> 0.7. Similarly, find void or near-to-void elements such that the pseudo-density is lower than 10%, i.e.,*ρ*_{soft}=*ρ*_{e}< 0.1.For each eigenvector, extract the displacement at every node for all

*ρ*_{hard}elements,*u*_{hard}, and at every node for all*ρ*_{soft}elements,*u*_{soft}. Do not consider nodes with boundary conditions applied whose displacement is already zero.Compute

*u*_{soft}/*u*_{hard}displacement ratios using the following metrics:- –
Percentile 90 ratio: PC90(

*u*_{soft}/*u*_{hard}) - –
Maximum value ratio: max(

*u*_{soft})/max(*u*_{hard}) - –
Average value ratio: mean(

*u*_{soft})/mean(*u*_{hard})

- –
Compute a compound ratio:

- –
$Ratiosoft/hard=(Percentile90+Maximumvalue+Averagevalue)$

- –
Identified localized modes whose compound ratio overpasses a threshold criterion established as Threshold = 1 × 10

^{3}, such that:- –
Localized mode identified if → Ratio

_{soft/solid}> Threshold

- –

*Step 4*. Compute the topology’s frequency response subject to harmonic loads.

*Step 5*. Analyze harmonic responses and identify an antiresonance frequency.

To properly identify antiresonances, the FRFs must be analyzed using the following steps:

Extract the FRFs on each of the nodes where the harmonic load has been applied.

- Compute the average FRF over all the aforementioned nodes:$averaged_FRF=mean(FRF)$
- Find local antiresonance peaks using the matlab function findpeaks:$[Amplitude,idx,Width,Prominence]=findpeaks(1/averaged_FRF)$
Evaluate local peak frequency absolute distance to target frequency

*f*:_{T}$Proximity=abs(Target_frequency\u2212Peak_frequency)$- Evaluate metrics for each local antiresonance peak identified, i.e., antiresonance bandwidth, amplitude, prominence, and proximity to the target frequency
*f*. Use the combined antiresonance evaluation factor:_{T}$factor=Amplitude+Prominence+Width\u2212Proximity$ Select the local antiresonance peak with the highest evaluation factor as the identified antiresonance frequency

*A*_{f}.

Note the following metric definitions [63]:

Amplitude: Signal value of a data sample that is either larger than its two neighboring samples or is equal to infinite.

Prominence: Measures how much the peak stands out due to its intrinsic height and its location relative to other peaks.

Width: Distance between the points where the signal intercepts a horizontal reference line positioned beneath the peak at a vertical distance equal to half the peak prominence.

Proximity: Measures the distance in Hertz (Hz) from the target frequency to the identified antiresonance frequency peak.

*Step 6*. Select an antiresonance eigenfrequency mode.

Once an antiresonance frequency *A*_{f} is identified, the corresponding topology’s harmonic displacement response is matched with an antiresonance mode shape (eigenvector) using the MAC [47]. The following steps are used to identify an appropriate eigenmode:

Extract the harmonic displacement response

*x*_{H}at the identified frequency*A*_{f}.- Compute the MAC between the harmonic displacement response
*x*_{H}and every antiresonance eigenvector $\Phi A$:$MAC=|\varphi ATxH|2(\varphi AT\varphi A)(\varphi ATxH)$ Select the eigenmode with the highest MAC coefficient as the identified antiresonance eigenmode; its corresponding eigenvalue

*λ*_{A}= (2*πf*_{q})^{2}is to be used during the optimization in Eq. (1).

## Appendix B: Impact of Penalization Factors on Localized Eigenmodes

Nonzero finite frequencies localized in the void or near-to-void regions are not preventable just by setting the stiffness and mass penalization factors to different values, e.g., *p*_{K} = 3, *p*_{M} = 1, or by modifying the stiffness-to-mass penalization ratio *p*_{K}/*p*_{M} in low-density regions. This problem is neither fully preventable by setting both penalization factors to the same value as we did in the case studies in Sec. 3, i.e., *p*_{K} = *p*_{M} = 1; however, the number of localized modes is greatly reduced, while the computation time is improved using equal penalization factors. To illustrate this statement, consider a simple demonstration. Figure 10 shows a design domain manually created, containing a solid nondesign stem surrounded by low-density elements. The stem comprises 2% of the available solid space; the remaining 98% of space has near-to-void density elements, a perfect scenario to observe localized modes in low-density regions.

The eigenproblem is solved using both the SIMP and the RAMP models, Eqs. (6) and (7), with both the same and different penalization factors for the stiffness (*p*_{K}) and mass (*p*_{M}) matrices. In the frequency range between 10 kHz and 50 kHz, the stem has three eigenmodes: two flexural modes and one extensional mode. Table 2 shows all the eigenfrequencies obtained with the SIMP model. Table 3 shows all the eigenfrequencies obtained with the RAMP model. We have used the localized eigenmode recognition routine from Appendix A to differentiate localized modes from real modes. Note that in the following tables, the localized eigenmodes are in italics, while the real modes are in black boldface. Finally, Table 4 shows the computation time required to solve each eigenproblem.

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | ||||
---|---|---|---|---|---|---|

18,511 | 17,876 | 12,511 | 15,524 | 17,787 | 18,983 | 19,947 |

19,099 | 19,153 | 12,537 | 15,629 | 17,963 | 19,000 | 19,988 |

23,272 | 21,770 | 12,710 | 15,650 | 17,984 | 19,018 | 20,055 |

24,233 | 22,670 | 12,756 | 15,709 | 18,093 | 19,035 | 20,061 |

26,876 | 25,644 | 13,070 | 15,902 | 18,123 | 19,069 | 20,155 |

35,773 | 35,104 | 13,154 | 15,942 | 18,130 | 19,231 | 20,454 |

36,377 | 35,369 | 13,488 | 16,525 | 18,240 | 19,295 | 20,488 |

37,952 | 37,443 | 14,004 | 16,633 | 18,254 | 19,301 | 20,578 |

40,050 | 38,567 | 14,418 | 16,731 | 18,293 | 19,354 | 20,587 |

40,636 | 39,061 | 14,588 | 16,847 | 18,429 | 19,364 | 20,608 |

43,469 | 42,500 | 14,595 | 16,896 | 18,551 | 19,368 | 20,670 |

44,489 | 43,492 | 14,611 | 17,087 | 18,569 | 19,502 | 20,816 |

45,850 | 45,199 | 14,620 | 17,308 | 18,609 | 19,529 | 20,925 |

46,572 | 46,117 | 15,128 | 17,567 | 18,635 | 19,547 | 21,056 |

48,305 | 47,864 | 15,188 | 17,613 | 18,685 | 19,662 | 21,059 |

49,712 | 15,469 | 17,631 | 18,885 | 19,663 | 21,145 |

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | ||||
---|---|---|---|---|---|---|

18,511 | 17,876 | 12,511 | 15,524 | 17,787 | 18,983 | 19,947 |

19,099 | 19,153 | 12,537 | 15,629 | 17,963 | 19,000 | 19,988 |

23,272 | 21,770 | 12,710 | 15,650 | 17,984 | 19,018 | 20,055 |

24,233 | 22,670 | 12,756 | 15,709 | 18,093 | 19,035 | 20,061 |

26,876 | 25,644 | 13,070 | 15,902 | 18,123 | 19,069 | 20,155 |

35,773 | 35,104 | 13,154 | 15,942 | 18,130 | 19,231 | 20,454 |

36,377 | 35,369 | 13,488 | 16,525 | 18,240 | 19,295 | 20,488 |

37,952 | 37,443 | 14,004 | 16,633 | 18,254 | 19,301 | 20,578 |

40,050 | 38,567 | 14,418 | 16,731 | 18,293 | 19,354 | 20,587 |

40,636 | 39,061 | 14,588 | 16,847 | 18,429 | 19,364 | 20,608 |

43,469 | 42,500 | 14,595 | 16,896 | 18,551 | 19,368 | 20,670 |

44,489 | 43,492 | 14,611 | 17,087 | 18,569 | 19,502 | 20,816 |

45,850 | 45,199 | 14,620 | 17,308 | 18,609 | 19,529 | 20,925 |

46,572 | 46,117 | 15,128 | 17,567 | 18,635 | 19,547 | 21,056 |

48,305 | 47,864 | 15,188 | 17,613 | 18,685 | 19,662 | 21,059 |

49,712 | 15,469 | 17,631 | 18,885 | 19,663 | 21,145 |

Note: The localized eigenmodes are in italics, while the real modes are in black boldface.

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | |
---|---|---|---|

18,679 | 18,203 | 16,455 | 36,994 |

18,791 | 19,048 | 17,131 | 37,906 |

23,271 | 23,271 | 17,159 | 38,805 |

24,232 | 24,232 | 18,100 | 39,216 |

26,875 | 26,875 | 19,004 | 40,077 |

35,774 | 35,774 | 25,296 | 42,393 |

36,377 | 36,377 | 25,722 | 42,557 |

37,952 | 37,952 | 26,836 | 43,743 |

39,967 | 39,156 | 28,320 | 44,560 |

40,050 | 40,050 | 30,738 | 44,836 |

43,469 | 43,469 | 31,459 | 45,650 |

44,489 | 44,489 | 32,421 | 46,911 |

45,850 | 45,850 | 32,932 | 47,159 |

46,572 | 46,572 | 34,157 | 47,828 |

48,305 | 48,305 | 35,985 | 49,636 |

36,755 |

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | |
---|---|---|---|

18,679 | 18,203 | 16,455 | 36,994 |

18,791 | 19,048 | 17,131 | 37,906 |

23,271 | 23,271 | 17,159 | 38,805 |

24,232 | 24,232 | 18,100 | 39,216 |

26,875 | 26,875 | 19,004 | 40,077 |

35,774 | 35,774 | 25,296 | 42,393 |

36,377 | 36,377 | 25,722 | 42,557 |

37,952 | 37,952 | 26,836 | 43,743 |

39,967 | 39,156 | 28,320 | 44,560 |

40,050 | 40,050 | 30,738 | 44,836 |

43,469 | 43,469 | 31,459 | 45,650 |

44,489 | 44,489 | 32,421 | 46,911 |

45,850 | 45,850 | 32,932 | 47,159 |

46,572 | 46,572 | 34,157 | 47,828 |

48,305 | 48,305 | 35,985 | 49,636 |

36,755 |

Note: The localized eigenmodes are in italics, while the real modes are in black boldface.

SIMP model | RAMP model | ||||
---|---|---|---|---|---|

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 |

68 | 71 | 523 | 70 | 69 | 90 |

SIMP model | RAMP model | ||||
---|---|---|---|---|---|

p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 | p_{K} = 1p_{M} = 1 | p_{K} = 3p_{M} = 3 | p_{K} = 3p_{M} = 1 |

68 | 71 | 523 | 70 | 69 | 90 |

These results demonstrate that using different penalization factors not only leads to significantly more localized modes but also to higher computation time, especially when the SIMP model is used. It also shows that when using the same penalization factors, fewer localized modes are obtained while improving computation time.

## Appendix C: Volume Constraints Study

Section 3 presents topologies where volume constraints are not active during the optimization process as they are not needed to converge to optimized solutions; the design objective is achieved regardless of volume constraints. Here, we present a comparative study on how different volume constraints influence the solution.

Consider the following demonstration to illustrate the volume constraint effect. For the same design problem of Sec. 2.1, a set of optimized topologies are obtained by varying the maximum volume constraint while keeping all the optimization parameters from Table 1 the same. The selected maximum volume constraints are 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, and 20%. Figure 11 presents the volume evolution as a function of the iteration number for each of the maximum volume constraints considered. The starting point is always the maximum volume fraction for all design variables, except for the bottom layer of elements, as explained in Sec. 3. Please note that the 10% volume constraint case was also considered, but this low percentage results in an overconstrained problem for which no feasible solution is found.

Note that in all the cases presented in Figure 11, the imposed volume constraint is never active, i.e., the optimization never attempts to push the volume to the constraint. Instead, the volume fraction is freely reduced by the optimization, converging to a certain volume fraction after 100 iterations. For the cases with the maximum volume constraint of 40%, 30%, and 20%, the convergence occurs around 17%. When the constraint is set to 50%, the convergence value is around a volume fraction of 27%. The convergence points for the other cases are not shown as they are not reached after 100 iterations. The convergence point is different because of disconnected members adding to the volume without decreasing the design objective; if the disconnected members are excluded, the converge point would be similar irrespective of the initially set volume constraint. This analysis indicates that the volume constraints are never active during the optimization process as they neither improve nor deteriorate the optimization results.

Figure 12 shows the topologies obtained at iteration 100 for each maximum volume constraint considered. Each subfigure shows the raw topology on the left-hand side and the binarized topology on the right-hand side. Note that the raw topologies show half of the structure as it includes a symmetry condition, while the binarized topologies present a postprocessed topology with its symmetry recovered.

The resultant topologies in Fig. 12 show that a volume constraint does not prevent disconnected members and that well-defined topologies are achieved regardless of the constraint value. Interestingly, when using the volume constraint of 90%, the resultant topology does not exhibit disconnected members. When the volume constraint is set to a lower value, the amount of available material is reduced; therefore, the optimization removes unnecessary material as seen for the cases of 40% and 30%; however, when the volume constraint is further reduced, the optimization finds solutions in a more constrained solution space, leading to different topologies and disconnected members, as shown in the case of volume constraint 20%. If the volume constraint is set to even lower values, i.e., 10% or 15%, the optimization becomes overconstrained and does not find feasible solutions. This observation suggests that disconnected members can be obtained with both high- and low-volume constraint values. Therefore, finding an appropriate constraint value requires a search that must be attempted heuristically for every design problem.

Disconnected members do not influence the solution as they do not contribute to the design objective. Consider the optimized topology shown in Fig. 13. The disconnected elements have been manually removed from the raw topology, and both frequency responses, with and without disconnected elements, have been computed. Please note that the typical postprocessing does not remove these disconnected elements directly from the raw topology. Instead, the raw topology is binarized to obtain fully solid members with properly defined boundaries, and then, disconnected parts can be discarded before meshing and simulating the final postprocessed topology.

Figure 14 compares the frequency responses for the original topology and the manually removed elements topology, demonstrating that both frequency responses are the same. Any differences come from postprocessing steps because small deviations from the raw topology might result in frequency response changes. This is especially critical for very small structural features that could change significantly after postprocessing. Despite some differences, the original and the postprocessed topologies achieved the design objective while maintaining high-quality metrics in their dynamic responses, as shown and discussed throughout the results section.

Differences between the raw and postprocessed frequency response can be minimized by adjusting the threshold value used during the binarization process. This threshold controls how the boundaries are defined; a higher threshold reduces the amount of solid material, and a lower threshold increases solid material. In the proposed design methodology, the threshold value is always set to 50% because this value does not promote solid material over void space, and vice versa. All the presented results were postprocessed using the same threshold value to maintain consistency; therefore, frequency response differences are expected.