Design of Locally Resonant Acoustic Metamaterials with Specified Band Gaps Using Multi-Material Topology Optimization

Locally Resonant Acoustic Metamaterials (LRAMs) have significant application potential because they can form subwavelength band gaps. However, most current research does not involve obtaining LRAMs with specified band gaps, even though such LRAMs are significant for practical applications. To address this, we propose a parameterized level-set-based topology optimization method that can use multiple materials to design LRAMs that meet specified frequency constraints. In this method, a simplified band-gap calculation approach based on the homogenization framework is introduced, establishing a restricted subsystem and an unrestricted subsystem to determine band gaps without relying on the Brillouin zone. These subsystems are specifically tailored to model the phenomena involved in band gaps in LRAMs, facilitating the opening of band gaps during optimization. In the multi-material representation model used in this method, each material, except for the matrix material, is depicted using a similar combinatorial formulation of level-set functions. This model reduces direct conversion between materials other than the matrix material, thereby enhancing the band-gap optimization of LRAMs. Two problems are investigated to test the method’s ability to use multiple materials to solve band-gap optimization problems with specified frequency constraints. The first involves maximizing the band-gap width while ensuring it encompasses a specified frequency range, and the second focuses on obtaining light LRAMs with a specified band gap. LRAMs with specified band gaps obtained in three-material or four-material numerical examples demonstrate the effectiveness of the proposed method. The method shows great promise for designing metamaterials to attenuate specified frequency spectra as required, such as mechanical vibrations or environmental noise.


Introduction
LRAMs have some interesting phenomena, one of which is that they can form subwavelength band gaps.This gives LRAMs significant potential for applications in isolating low-frequency vibrations or noise using a small-sized structure.Liu et al.'s work demonstrated that composites consisting of LRAMs exhibit subwavelength band gaps [1].Later, an analytical model was introduced where effective mass densities might become negative near local resonances [2].Acoustic metamaterials, with their inherent deep subwavelength nature, have triggered many exciting investigations, as reviewed in [3,4].
In this work, a topology optimization method is used to obtain multi-material LRAMs with the objective of isolating specified frequency ranges of vibration or noise.Using this proposed method, LRAM band gaps can be designed as required within a frequency range.The frequency limitations when the unit cell size and materials are specified can be found in Roca et al.'s work [5].Unlike traditional methods, topology optimization methods can provide innovative structural configurations under different objectives and constraints [6].
As reviewed in Li et al.'s paper [15], since topology optimization techniques were introduced in acoustic metamaterial designs, many innovative works have been achieved.The target of this paper is to study the band-gap optimization problems of LRAMs, a topic that has attracted great interest from researchers.In 2003, Sigmund et al. utilized a topology optimization technique for the first time to obtain wider band gaps in phononic crystals [16].T. Matsuki et al. [17] used the optimization technique to design LRAMs.Yang et al. [18] maximized the first band gap of LRAMs using the concept of effective mass density.Roca et al. [5] introduced an approach to maximize LRAM band gaps, where a multiscale homogenization framework combined with model order-reduction techniques was introduced to build the optimization model.Zhang et al. [19] designed band gaps in acoustic metamaterials based on a material-field series expansion framework.To obtain low-frequency broad band gaps, Sun et al. [20] improved the scatterer filling scheme in two hierarchical honeycomb metamaterials.Li et al. [21] sought anisotropic hierarchical honeycomb acoustic metamaterials with multiple broad band gaps using a topology optimization method.
Since an LRAM typically consists of a matrix, coating, and scatterer, the material parameters of these components can be significantly different [1,2,18].For instance, the stiffness of the coating is lower by several orders of magnitude compared to that of the scatterer.This makes it difficult to use a single material to obtain LRAMs with band gaps relying on the local resonance mechanism, whereas the use of multiple materials can make this easier.This paper introduces a method that can utilize multiple materials to obtain LRAMs with band gaps that satisfy specified frequency constraints.Clear boundaries among the individual constituent materials are beneficial for LRAM band-gap optimization with multiple materials.Therefore, this work utilizes a parameterized level-set method (PLSM).Almost two decades ago, the LSM was introduced to address structural optimization problems [22,23].Later, by modifying the original method, the PLSM was presented for shape and topology optimization [24,25].Wang and Wang [26] presented a model for multi-material problems, in which n = 2 m material phases are denoted by m level-set functions.Wang et al. [27] introduced a model for multi-material problems, where a combined formulation of level-set functions is used to represent each phase.Using multiple materials, Vogiatzis et al. [28] adopted a reconciled LSM to address the optimization problems of negative Poisson's ratio metamaterials.Liu and Ma [29] introduced another model that employed p variables to denote p materials along with the void.Xia and Shi [30] employed an ESO technique to tackle hole nucleation problems in an LSM for multi-material optimizations.Using the alternating active-phase algorithm and the LSM, Sha et al. [31] effectively addressed optimization problems involving multiple materials.Besides LSMs, there are certainly other types of multi-material approaches.Bendsøe and Sigmund [32] proposed a strategy using the SIMP method to interpolate two materials and a void phase.Hvejsel and Lung [33] expanded on the aforementioned SIMP method by introducing a generalized interpolation scheme that can be applied to any number of materials.Gao and Zhang [34] discussed density-based multi-material interpolation techniques through numerical tests and theoretical analysis.Tavakoli and Mohseni [35] suggested an alternating active-phase algorithm, where the problem was partitioned as a series of two-phase sub-problems.Zuo and Saitou [36] presented the ordered SIMP approach, aimed at improving computational efficiency.Huang and Xie [37] introduced an ESO method for optimization problems involving both single and multiple materials.
Despite considerable progress, the topology optimization method for LRAMs is still not mature.Most current research does not involve obtaining LRAMs with specified band gaps.However, such LRAMs have significant application potential.This is because LRAMs can form subwavelength band gaps, which enable them to isolate vibrations or noise in specified low-frequency ranges.Utilizing multiple materials in LRAM optimization may help achieve results with a sufficiently wide band gap while simultaneously meeting other objectives, such as mass reduction.Therefore, this work proposes a multi-material topology optimization approach for designing LRAMs with specified band gaps.The band-gap calculation in this method uses a homogenization framework [38,39], which establishes both a restricted subsystem and an unrestricted subsystem to determine the band gaps.These subsystems are specifically tailored for LRAMs, facilitating the opening of band gaps based on the local resonance mechanism during optimization.Sensitivity analysis only needs to consider these two subsystems.Consequently, sensitivity information related to band gaps can be calculated without relying on the Brillouin zone.This implies that the proposed method simplifies the optimization calculation.In the representation model of this multi-material optimization method, each material, except for the matrix material, is represented by a similar combinatorial formulation of level-set functions.This model reduces the direct conversion between materials other than the matrix material, thereby enhancing the band-gap optimization of LRAMs.The effectiveness of this method is tested by applying it to solve two types of LRAM optimization problems.The first test involves maximizing the band-gap width while ensuring it encompasses a specified frequency range, and the second test focuses on obtaining light LRAMs with a specified band gap.Numerical examples successfully yield LRAM designs with the required band gaps using three or four materials, demonstrating the feasibility and effectiveness of this approach.In our future work, we will manufacture the obtained LRAMs and analyze their properties experimentally, similar to studies of other types of materials [40,41].
This paper is structured as follows.Section 2 introduces a simplified calculation method for LRAM band gaps.Section 3 proposes the mathematical model for multi-material LRAMs.Section 4 utilizes numerical examples to prove the feasibility of the presented method.The conclusions are provided in Section 5.

A Simplified Calculation Method for LRAM Band Gaps Based on a Homogenization Framework
Generally, the elastic wave field in LRAMs can be represented as a periodic function.Following the Bloch theorem, band-gap analysis can be simplified to consider only a single LRAM cell [42][43][44][45][46].The square unit cell under study and its corresponding irreducible Brillioun zone are presented in Figure 1.According to existing research [16,47,48], the computation of band diagrams is reduced by only considering the boundary of the irreducible Brillouin zone.Therefore, in conventional optimization methods, calculating the band gap of LRAMs typically involves traversing the Γ-X-M-Γ path of the Brillouin zone.The computational cost is significantly impacted by the number of selected wave vectors [16,49,50].This paper introduces a simplified calculation method for LRAM band gaps based on a homogenization framework.The eigenfrequencies obtained from two subsystems are used to determine the LRAM band gaps.

Dispersion Curve Analysis of LRAMs
Generally, one can identify whether there are band gaps through dispersion curve analysis of LRAMs.The propagation characteristics of elastic waves in isotropic elastic materials have been extensively researched.In this paper, the analysis method used in Li et al.'s work [49] is adopted.The governing equation is where γ is the location vector, λ and µ represent the Lame's coefficients, u denotes the displacement vector, and ρ is the mass density.For two-dimensional LRAMs, the elastic wave propagation mode includes two coupled in-plane modes and one out-of-plane mode.
In this paper, only in-plane waves are considered.
In an infinite periodic structure, according to Bloch's theorem, the displacement vector u(γ, k) is expressed as where u k (γ) is a periodic function of the amplitude displacement and k = k x , k y represents the wave vector.By substituting Equation (2) into the governing equations for in-plane waves, one obtains a generalized eigenvalue equation that can be solved using finite-element methods [51,52]: where K and M are the stiffness and mass matrices, respectively; v denotes the eigenvectors; and ω represents the angular frequency.

The Establishment of Subsystems Based on the Homogenization Framework
The homogenization framework utilizing the Craig-Bampton reduction technique [38] and the framework based on the extended Hill-Mandel principle [39] is tailored for modeling the phenomena involved in LRAMs.Both frameworks can be used to calculate LRAM band gaps.When applying the framework presented by Roca et al. [39] in practice, boundary conditions need to be considered [5].Therefore, by incorporating techniques from this framework with the Craig-Bampton method [53], the band-gap calculation method can be derived in a simple and direct way.The method presented in Roca et al.'s work [39] is used to establish the relationship between the microstructure and the macroscale and to capture the local resonance phenomenon.The Craig-Bampton technique is employed to generate compact, reduced models of microstructures, aiding in the computation of local resonance effects.Based on these two works, an unrestricted subsystem and a restricted subsystem are established to help calculate the band gaps of LRAMs.This method incorporates the characteristics of the two previously mentioned homogenization frameworks and can be easily implemented using finite-element methods.
In the homogenization scheme, to maintain the separation of scales, it is essential to ensure that the microstructure size is far smaller than the wavelength.This work is based on two main hypotheses: (1) kinematic admissibility between scales; and (2) the principle of virtual work, which states that the total internal virtual work remains constant for a dynamic system regardless of the imposed kinematics.Through homogenization, the solution to the microscale problem can be used to derive the macroscale constitutive response.Since this paper focuses on establishing the optimization computational model, more detailed information about the theoretical aspects can be found in [38,39].
The selected representative volume element (RVE) is the unit cell of the LRAM.This paper only studies two-dimensional problems without considering plane rotation.The response of the RVE to the applied force can be evaluated by solving where M represents the mass matrix, K denotes the stiffness matrix, u and ü represent the displacement and the acceleration, and F stands for the applied force.The responses can be split into two systems using the decomposition method.One is a quasistatic system, and the other is an internal dynamic system.The overall structural response can be represented as the sum of its quasistatic response and its internal dynamic response.The internal dynamic system can be built as an unrestricted subsystem or a restricted subsystem, through which one can determine the LRAM band gap [39].The contribution of the restricted subsystem to the macroscale consists of the rigid body displacement and the internal dynamic response.The internal dynamic response can be projected onto the space spanned by the eigenmodes of the microstructure, with the prescribed nodes fixed.This entire process is illustrated in Figure 2. The unrestricted subsystem is a generalized eigenvalue problem, whose contribution to the macroscale is discussed in Section 2.3.2. where For simplicity, the typical two-dimensional square unit cell is considered, as depicted in Figure 2. The system is divided into tied nodes and retained nodes to facilitate the application of periodic boundary conditions.The tied nodes consist of the right and top edges and the upper-right vertex.In the calculation, the retained nodes are subdivided into prescribed nodes and free nodes to facilitate problem-solving [38].Since the responses of tied nodes can be obtained through other nodes based on periodic boundary conditions, the LRAM cell can be partitioned into prescribed nodes, denoted by '1', and remaining nodes, denoted by '2', to derive the response where U µ and Üµ represent the displacement and acceleration associated with each node.
In calculating the quasistatic response, the mass contribution of the system is omitted.The following relationship is derived using the condensation method: The target of this paper is to study LRAM band-gap optimization problems, so we do not pay much attention to the quasistatic response.
In the restricted dynamic subsystem, the prescribed nodes do not participate in the internal vibrations.Based on the second equation of system (7), one can obtain Use ( •) to indicate that the periodic boundary condition is imposed; thus, Through Equation ( 10), the following subsystem can be obtained: where M * µ2 and K * µ2 represent the mass and stiffness matrices obtained after imposing the restrictions, respectively.Û * µ2 is the fluctuation field vector, and D represents a matrix coupling the accelerations of the micro-and macroscales.
The macroscale inertial force derived is as follows: Derivation details for Equations ( 11) and ( 12) are given in Appendix A.

A Simplified Calculation Method for LRAM Band Gaps Using Subsystems
In this section, the contribution of the microfluctuation field to the macroscale is analyzed.Subsequently, a simplified calculation method for the band gap of LRAMs is established.
In this paper, wave solutions are considered.The vector of nodal values Ûµ is where U represents the amplitude function of Ûµ , x denotes the spatial coordinate, ω represents the angular frequency, n κ signifies the wave propagation direction, and κ stands for the matching wavenumber.

The Map of the Restricted Subsystem to the Macroscale
This subsection analyzes the contribution of the microfluctuation field from the restricted subsystem to the macroscale.For the interior free vibration system, the eigenvalue problem is given as where ψ * µ2 is the natural vibration mode and ω * µ is the corresponding eigenfrequency.Solving the generalized eigenvalue problem of interior dynamics yields eigenmodes.These eigenmodes are then used as a reduced basis in subsequent calculations to capture the local resonance phenomenon and derive the macroscale inertial force where R represents the effective density tensor, the derivation of which can be found in Appendix A. Q is a coupling matrix.The symbol Ω * µ represents a diagonal matrix that holds the natural frequencies ω * µ of the restricted subsystem.The derivation process is given in Appendix B.1.
Using R e to represent the effective pseudo-density tensor, one can obtain where

The Map of the Unrestricted Subsystem to the Macroscale
Consider the generalized eigenvalue problem formed by the unrestricted subsystem, where the design domain is not partitioned into prescribed and retained nodes where N µ represents the shape function.( •) denotes that the periodic boundary condition has been imposed.The macroscale inertial force β can be obtained through the generalized eigenvalue problem of the unrestricted subsystem where ψN µ represents the natural vibration mode, and ω N µ is the corresponding eigenfrequency.The derivation process is given in Appendix B.2.

Band-Gap Analysis through the Two Subsystems
The average value of elements on the main diagonal of the effective pseudo-density matrix R e is the effective mass density.According to Equation (17), the values of elements on the main diagonal of R approach negative infinity as ω approaches ω * from the right, and R e obviously follows the same trend.
From Equation ( 19), it can be found that when ω approaches ω µ , With so the effective, frequency-dependent pseudo-density R e is 0. With the above work, consider the frequency range [ω * µ , ω µ ] determined by the natural frequencies of the restricted and unrestricted subsystems.Within this frequency range, R e initially is unbounded and negative-definite, eventually evolving to 0 at the end of the range.Therefore, [ω * µ , ω µ ] can be used to determine the band gap of LRAMs [39]. Figure 3 presents the dispersion curves and effective mass density (EMD) curves obtained from different LRAMs.The band gaps shown in the two types of curves are compared to demonstrate that the simplified band-gap calculation method is feasible.The lattice size of the LRAM cell is 0.02 × 0.02 m.The design domains presented in the figures contain 3 × 3 cells.The materials utilized include lead (Pb), rubber (NR), and epoxy (EP), with their corresponding material parameters detailed in Table 1.In the dispersion curves in Figure 3, the yellow zone denotes the band gap determined by the bands.In the EMD curves in Figure 3, the yellow zone representing the band gap is determined by the eigenfrequencies of the restricted and unrestricted subsystems.As observed in the dispersion curves and EMD diagrams, the yellow zones representing the first band gap only show a difference at the lower edge, with a value of less than half a percent (less than 1 Hz).Thus, this simplified calculation method can be utilized in the optimization process.Since the simplifying assumptions are more valid as the system approaches a quasistatic situation, this paper only studies the optimization problem of the first band gap.

A Multi-Material Optimization Method for LRAMs
A level-set-based method is introduced to address LRAM band-gap optimization problems using multiple materials in this work.The compactly supported radial basis function (CSRBF) is employed in the development of the PLSM [54].

The Uniform Multi-Material Description Model
Inspired by the "color" level-set model [26] and the uniform multiphase materials interpolation (UMMI) method [34], a level-set-based uniform multiphase description model is introduced.In this model, each material, except for the matrix material, is denoted by a similar combination form that utilizes level-set functions, as shown below where χ k represents the kth material and p is the total number of material phases.Using p = 3 as an illustration, two variables are employed to represent three materials within the design domain where H is the smeared Heaviside function where ∆ specifies the breadth of the numerical approximation of H.
To be more illustrative, take Figure 4 as an example to illustrate the description model.Three materials are represented by two level-set functions in this illustration.As demonstrated in Equation ( 23), when representing materials apart from the matrix material, the respective level-set function takes on a value of 1, whereas the remaining level-set functions are set to 0. This implies that the material type can only change to another type when one level-set function changes from 1 to 0, another changes from 0 to 1, and the retained ones remain unchanged.Otherwise, the material type will change to the matrix material.The mutual exclusivity of material descriptions makes direct conversion between materials other than the matrix material less smooth in the presented model.According to our numerical experience, this is beneficial for LRAM band-gap optimization.One possible reason is that the coating's stiffness is significantly lower by several orders of magnitude compared to that of the scatterer.This means that direct conversion between the coating and the scatterer may lead to a relatively large change in the objective value.

The CSRBF-Based Parameterized Level-Set Method
According to existing research [22,23], variables in LSMs are updated using the Hamilton-Jacobi partial differential equation (PDE) given below: where t represents pseudo-time, V q = V • − ∇ϕ |∇ϕ| denotes the normal velocity toward the outside, and ∇(•) represents the gradient.
By decomposing the original Hamilton-Jacobi PDE into a system of coupled ordinary differential equations (ODEs) using CSRBFs, the method known as the parameterized level-set method is proposed [54,55].The level-set function ϕ(x, t) is interpolated by n CSRBFs, which can be given as where the design variable α(t)= [α 1 (t),α 2 (t), . . ., α n (t)] consists of expansion coefficients that vary with pseudo-time t. g(x)= [g 1 (x),g 2 (x), . . ., g n (x)] is the vector of CSRBFs.This work adopts CSRBFs with C2 continuity.At the specified x, g i (x) with r as the support radius is given as follows: Then, Equation ( 26) transforms into the following governing equation where The optimization proceeds by updating α(t) until a convergence result is achieved.

Mathematical Model for LRAM Band-Gap Optimization Problems with Multiple Materials
Based on the work described above, multi-material band-gap optimization models for LRAMs with different objectives and constraints are built.In these models, the natural frequencies of the two subsystems established in Section 2 are utilized to determine the LRAM band gap.

Maximizing the Band-Gap Width While Ensuring It Encompasses a Specified Frequency Range
Obtaining a wider band gap at a lower frequency is significant for applications.This is typically the objective of LRAM band-gap optimization.However, in certain cases, the band gap of the LRAM obtained should encompass a specified frequency range to isolate the vibration or noise as required.The optimization model for maximizing the band-gap width while encompassing a specified frequency range is established as follows: (31) where N denotes the finite-element quantity, m represents the number of material phases, ω 0 represents the center frequency, and ω r denotes the deviation of the band-gap boundaries from the center frequency.These parameters are utilized to specify the required frequency range.

Obtaining Light LRAMs with a Specified Band Gap
To precisely control the propagation of vibrations or sound waves, optimization problems aimed at obtaining LRAMs with a specified band gap are also studied.One significant application of these LRAMs is in filter design.Meanwhile, achieving a lighter structure is a crucial objective considered in the design process.Here, the optimization model for obtaining light LRAMs with a specified band gap is established.(32) where the same parameters have the same meanings as those in the mathematical model in Equation (31).ρ i is the mass density, and V i is the volume fraction.ω o1 and ω o2 are the target frequencies of ω * µ and ω µ , respectively.According to our numerical experience, it may be difficult to achieve convergence when setting the lower and upper limits of the band gap at specified values.So, in implementation, the lower and upper frequency constraints will be set as certain frequency ranges.ω r1 is used to represent the deviation of ω * µ from ω o1 , and ω r2 is used to represent the deviation of ω µ from ω o2 .Take the specified target frequency range of 500-700 Hz as an example, the deviation of the obtained frequencies from the specified constraint frequencies is permitted within a range of ten percent.ω o1 is 500 Hz and ω o2 is 700 Hz. ω r1 and ω r2 are obtained as 0.1 × (0.5 × (700 − 500)) = 10 Hz.So, ω * µ ranges from 490 to 510 Hz, and ω µ ranges from 690 to 710 Hz.

Sensitivity Analysis
As the method presented is gradient-based, conducting sensitivity analysis is necessary in the optimization process.The obtained sensitivity information is used to determine the search direction.
The sensitivity information for the objective described in Equation ( 31) with respect to χ can be calculated as follows: Normalizing the eigenvectors to the global mass matrix, one can obtain the sensitivities of eigenfrequencies ω by The sensitivity information for the objective described in Equation ( 32) with respect to χ can be calculated as follows: where χ i denotes the ith material.v 0 is a vector whose elements denote the volume fraction of solid finite elements.The derivative of H(ϕ) is the Dirac function, given by Considering Equations ( 27) and ( 36), the sensitivity of f is given as Further details about the sensitivity of K and M are presented in Appendix C. Note that the sensitivity analysis shown above is limited to the problem with a single eigenvalue.When the system has multiple eigenvalues, the sensitivity information should be calculated using the method presented in [56][57][58].
This work adopts the method of moving asymptotes (MMA) [59] to update the design variables.

Results and Discussion
This section presents several 2D three-material and four-material numerical cases to confirm the feasibility and utility of the approach.The impact of the initial guesses on the optimized results is also reported.The square unit cell studied has square symmetry [16,60,61], exhibiting mirror symmetry along the horizontal, vertical, and two diagonal lines.The lattice size of the studied LRAM cell is 0.02 × 0.02 m.The results presented in the topology figures contain 3 × 3 cells.All calculations were conducted using MATLAB R2020b.The materials used include lead, aluminum, rubber, and epoxy, with their material parameters detailed in Table 1.The MMA parameters used are a 0 = 1, a = 1, c = 1000, d = 1, and move = 0.1.The meaning of these parameters can be found in Svanberg's paper [62].In these numerical examples, only the band gap between the third and fourth bands is studied.The initial guesses and corresponding band diagrams are presented in Figure 5, unless otherwise specified.Obviously, there are no band gaps in the initial designs.It should be pointed out here that the final band diagrams presented in all numerical cases are calculated using the conventional method based on the boundary of the irreducible Brillouin zone.In all band diagrams in this section, black vertical dotted lines are used to divide different boundary ranges of the irreducible Brillouin zone, and blue and red horizontal dotted lines are used to indicate the band gap zones.In the optimized results, the black region is lead, rubber is depicted in red, and the green region corresponds to epoxy.The scatterer is made of lead, which has a high mass density, while the coating is made of rubber, a soft material.This is in line with the band-gap formation mechanism in LRAMs, as described in the existing literature [1,2,18,63,64].
To investigate the impact of different mesh resolutions on the final results, mesh resolutions of 80 × 80 and 160 × 160 are also employed for the optimization.It can be observed that only slight differences appear at the lower limits.The objective values are 3.02, 3.14, and 3.05, respectively.In the three resulting topologies, lead occupies 0.63, 0.66, and 0.67 of the total volume, and rubber occupies 0.24, 0.23, and 0.22 of the total volume.The remaining region is occupied by epoxy.These results indicate that the mesh resolution impacts the final optimized topology.However, with the same objectives and constraints, the resulting topologies will be similar, and their band gaps will meet the specified frequency constraints while being wide enough.Since the numerical examples aim to verify the feasibility of the presented approach, a resolution of 40 × 40 is used in the remaining examples.
Four materials are also used to obtain the required LRAM under the same objective and constraints.The obtained topology and band diagram are presented in Figure 8. Figure 9 illustrates the evolutionary history of the four-material case.In the obtained topology, black, yellow, red, and green denote lead, aluminum, rubber, and epoxy, respectively.The band gap of the optimized topology is 388.9-1559.5 Hz, which is also very close to 388.3-1559.5 Hz, as determined by the simplified calculation method.The objective value is 3.01.The final volume fractions of lead, aluminum, rubber, and epoxy are 0.6225, 0.005, 0.24, and 0.1325, respectively.Due to changes in material distribution and the band gap that appeared, the objective value experienced acute changes in the initial iteration steps of the optimization.Subsequently, the optimization converges, meeting the frequency constraints.As shown in Figure 9, the aluminum contained in the initial guess almost disappears (decrease to 0.005) in the obtained topology.According to the band-gap formation mechanism [1,2,18,63,64], LRAMs can be composed of scatterers, coatings, and matrices.The scatterer is composed of a material with a high mass density, while the coating is composed of a soft material.Once the coating is determined, the band-gap width is influenced by the difference in mass density between the matrix and scatterer.Furthermore, the higher the mass density of the scatterer and the softer the coatings, the lower the initial frequency of the band gap.The mass density of lead is much higher than that of aluminum.Therefore, it is easier to achieve a low-frequency band gap with a better objective value using lead as the scatterer instead of aluminum.This is a possible reason why the volume fraction of aluminum decreases and is substituted by lead.

Numerical Examples Where the Obtained Band Gaps Should Encompass 1000-2000 Hz
In this subsection, three-material and four-material examples where the obtained band gaps should encompass the frequency range of 1000-2000 Hz are presented.Figure 10 presents the topologies obtained and their corresponding band diagrams.In the optimized results, black, yellow, red, and green represent lead, aluminum, rubber, and epoxy, respectively.
The final band gaps of the three-material and four-material cases are 552.3-2964.3Hz and 597.6-3011.5 Hz, respectively.The band gaps determined by the simplified method are 552.1-2964.3Hz and 597.4-3011.5 Hz, with only slight differences at the lower limits.In the three-material case, lead, rubber, and epoxy account for 0.65, 0.15, and 0.2 of the final design in volume, respectively.In the four-material case, lead, aluminum, rubber, and epoxy occupy 0.48, 0.04, 0.28, and 0.2 of the final design in volume, respectively.The objective values of the two obtained results are 4.37 and 4.04, and the masses are 3.17 and 2.52, respectively.In the four-material case, the final objective value decreased by 7.6% compared to the three-material case, while the mass of the structure decreased by 20.4%.Sections 4.1.1 and 4.1.2show that the method successfully obtains LRAMs with sufficiently wide band gaps encompassing the specified frequency ranges.This demonstrates the method's effectiveness in using multiple materials to solve band-gap optimization problems in LRAMs.It can be seen from the numerical examples of three and four materials that using more types of materials in optimization does not necessarily result in a wider band gap.This is because LRAM band gaps are influenced by many factors, including material properties, topology shapes, symmetries, filling fractions, and more.However, when the lattice size is specified, material properties become the predominant factor in determining LRAM band gaps.In the four-material examples, since a certain volume of lead is substituted with aluminum, it can be difficult to obtain a wider band gap compared to the three-material examples under certain circumstances due to differences in the parameters of the two materials.But using more kinds of materials to design LRAMs can obtain results that meet more goals and constraints.For example, in Section 4.1.2,the addition of aluminum in the four-material case results in a lighter structure compared to the three-material case, with a much smaller decrease in the band-gap width relative to the mass.In Table 2, a comparison is presented to show that a lighter design can be obtained while similar performance is achieved using more materials.In this section, numerical examples of obtaining light LRAMs with a specified band gap are presented.Following the rule mentioned in Section 3.3.2, the band gaps of the optimized results are set as 500 ± 10-700 ± 10 Hz, 700 ± 10-900 ± 10 Hz, 800 ± 40-1600 ± 40 Hz, and 700 ± 90-2500 ± 90 Hz, respectively, regardless of whether they are three-material or four-material cases.

Three-Material Cases
Lead, rubber, and epoxy are utilized to obtain the required LRAMs in these cases.The optimized results are shown in Figure 11.In cases (a) to (d), lead occupies 0.17, 0.1, 0.14, and 0.34 of the total volume, while rubber occupies 0.52, 0.42, 0.33, and 0.16, respectively.The remaining region is occupied by epoxy.The masses of these obtained structures are 1.22, 0.92, 1.07, and 1.89, respectively.The band gaps of the optimized topologies are 497.3-705.2Hz, 701-904.2Hz, 794.4-1566.1 Hz, and 711.4-2468.5 Hz, respectively.The band gaps determined by the proposed simplified method are also provided, which are 496.1-705.2Hz, 699.9-904.2Hz, 795.5-1566.1 Hz, and 711.1-2468.5 Hz, respectively.Figure 11 shows that LRAMs with the specified band gaps are successfully obtained.The proposed method has the potential to design devices with specified functions, such as filters.
(a)  The results depicted in Figure 12 demonstrate that this approach can also employ four materials to obtain LRAMs with a specified band gap.In this optimization problem, whether using three or four materials, the final band gaps and those determined by the proposed simplified method show differences of less than half a percent (less than 2 Hz), also only at the lower limits of the band gaps.This demonstrates that the simplified band-gap calculation method can be applied to the optimization process.
With more kinds of materials providing more design possibilities, multi-material topology optimization methods may have the ability to achieve results that meet more objectives and constraints.Since structural mass is an important factor that impacts the structure's application, it should be taken into consideration during the design stage.Table 2 presents a comparison of the changes in the band gaps and masses of the resulting topologies using three or four materials described in Section 4.2.
Table 2 shows that compared to using three materials, lighter designs with the required band gaps can be obtained using four materials.While this conclusion may not be general, it demonstrates the possibility of utilizing a greater variety of materials to design lighter LRAMs without significantly altering the band-gap width.Therefore, this multi-material approach is capable of obtaining LRAMs with the required band gap while also reducing their mass.

The Impact of Initial Guesses on the Final Optimized Results
To investigate the impact of the initial guess on the final result, numerical examples with different initial guesses are provided in Table 3.The objective of these examples is to obtain LRAMs with the widest band gap while ensuring that the band gap encompasses the specified frequency range of 400-600 Hz.In these cases, if the initial guess is not asymmetrical, the method described in Dong et al.'s paper [65] is used to calculate the band gap.
Table 3 shows that different initial guesses can result in different final topologies.This indicates that numerous locally optimal results exist in multi-material LRAM band-gap optimizations.However, if all the materials are sufficiently discretely distributed in the initial guesses, topologies with band gaps that fulfill the objectives and constraints can be obtained.As a result, the proposed method can significantly assist in designing LRAMs with the required band gaps.
As observed from the numerical examples in this section, the proposed method can use multiple materials to obtain LRAMs with specified band gaps.This method can be employed to design metamaterials that attenuate vibration and noise as required.However, challenges remain that affect its practical application, such as selecting materials for specific applications and incorporating the impact of manufacturing uncertainty into the optimization process.

Conclusions
In this paper, based on the homogenization framework, a level-set-based multimaterial topology optimization method is proposed for designing LRAMs with specified frequency constraints.The method has the following features: (a) A simplified calculation method for LRAM band gaps is established using the homogenization framework.This method utilizes a restricted subsystem and an unrestricted subsystem to determine the LRAM band gap.
(b) In the presented uniform level-set-based multi-material description model, each material, except for the matrix material, is denoted by a similar combination formulation of level-set functions.This reduces direct conversion between materials other than the matrix material.
(c) The presented multi-material optimization method can obtain LRAMs with band gaps that fulfill specified frequency constraints.The two problems of obtaining LRAMs with the maximum band-gap width while ensuring that this band gap encompasses a specified frequency range, and obtaining light LRAMs with a specified band gap, are tested using multiple materials.
The results of numerical cases using three or four materials confirm the feasibility and effectiveness of the presented method in obtaining LRAMs with the required band gaps.Compared with similar optimization strategies presented by Roca et al. [5], the method proposed here removes limitations such as considering the matrix fixed frame as infinitely stiff and the coating material as massless.Moreover, the presented method does not restrict the analysis to a single dimension and can obtain LRAMs with the required band gaps using multiple materials.This method can obtain LRAMs with wide low-frequency band gaps based on the local resonance mechanism.Thus, it shows great promise for attenuating specified frequency spectra as required, such as mechanical vibrations or environmental noise.
The proposed method can use multiple materials to solve LRAM optimization problems, potentially achieving results that meet more objectives and constraints.For instance, this work shows that this approach is capable of using a greater variety of materials to obtain LRAMs with the required band gaps while simultaneously reducing their mass.However, it should be mentioned that due to the limitations of simplifying assumptions, when the system approaches the quasistatic case, the assumed simplifications hold more strongly.As the frequency increases, the simplified band-gap calculation method may become less accurate as the separation of scales approaches its limits.However, according to our numerical experience, in most first-band-gap optimization problems, the proposed method can obtain the required LRAMs that meet objectives and constraints.Additionally, higher-frequency band gaps can be easily obtained by optimizing LRAMs with a smaller lattice size.
Multi-material topology optimization is often criticized.One reason is that the obtained results face challenges in properly addressing interfaces between different materials in the manufacturing process.The fabricated metamaterials should ensure that the different materials (where material parameters can be significantly different) are properly connected, with minimal defects and a high degree of precision.Three-dimensional printing technology may help facilitate practical applications.Additionally, the impact of material interfaces and manufacturing uncertainty is not considered in this work.These aspects will be explored in our future research.Applying the method to solve three-dimensional and large-scale problems, as well as designing other types of metamaterials, will also be studied in future work.
As mentioned before, the displacements of the prescribed nodes are set as the macroscale displacements u.Together with Equations (A13) and ( A18 represents the corresponding eigenfrequency. Using q * µ2 to denote the amplitude of the eigenmodes, the solution Û * µ2 can be obtained as Û * µ2 =Ψ * µ2 q * µ2 .(A36) Substituting Equation (A36) into Equation (11) and then multiplying both sides of Equation ( 11) by Ψ * T µ2 , one can obtain and by taking Equation (A34) into consideration, one can obtain where q * u2 represents the modal amplitude vector.ω is used to denote the excitation macroscale frequency, so q * µ2 = − ω 2 q * µ2 .(A40)

Appendix C. Detailed Derivation of the Sensitivity of the Stiffness and Mass Matrices
The global stiffness matrix or global mass matrix can be assembled by the element stiffness matrix or element mass matrix, respectively.The design domain contains p kinds of materials, described as One can compute the element stiffness matrix K e by where K m e , m = 1, 2, • • • , p represents the mth material's contribution to K e .Similarly, the element mass matrix M e can be calculated by where M m e , m = 1, 2, • • • , p represents the mth material's contribution to M e .Using p = 3 as an example, three materials are represented by two level-set functions.Using Equation (24), K e and M e can be represented as and Using H m to denote H(ϕ m ), the sensitivity of K e with respect to ϕ 1 can be calculated as ∂K e ∂ϕ 1 = ((1 The sensitivity of K e with respect to ϕ 2 is Similarly, the sensitivities of M e with respect to ϕ 1 and ϕ 2 are and respectively. In Equations (A58) to (A61), δ i , i = 1, 2 is the derivative of H(ϕ i ), i = 1, 2, defined in Equation (36).
Using Equation ( 27), the sensitivity of K e with respect to α m is The sensitivity of M e with respect to α m is ∂M e ∂α m = g(x) −1 ∂M e ∂ϕ m .(A63)

Figure 1 .
Figure 1.The square unit cell and its corresponding irreducible Brillouin zone (blue region).

Figure 2 .
Figure 2. Illustration of the total response of the structure and the representation of the periodic boundaries in a typical 2D square unit cell.Based on this, the subsystems are established.Use β to represent the local macroscale homogenized force

Figure 3 .
Figure 3. Different LRAM unit cells with corresponding dispersion curves and effective mass density curves.

Figure 4 .
Figure 4.The description model involves representing three materials using two level-set functions: ϕ 1 and ϕ 2 .

Figure 5 .
Figure 5.Initial guesses and corresponding band diagrams of numerical cases in Sections 4.1 and 4.2 (in the band diagrams, the blue and red lines overlap due to the absence of band gaps): (a) three materials; (b) four materials.

4. 1 .
Numerical Examples of Maximizing the Band-Gap Width While Ensuring It Encompasses a Specified Frequency Range In this subsection, numerical examples of maximizing the band-gap width encompassing a specified frequency range are presented.In the three-material cases, lead, rubber, and epoxy are utilized.In the four-material examples, lead, aluminum, rubber, and epoxy are employed.4.1.1.Numerical Examples Where the Obtained Band Gaps Should Encompass 400-600 Hz First, three materials are used to obtain LRAMs with the required band gaps.The design domain is partitioned into linear four-node elements with a resolution of 40 × 40. Figure 6 presents the obtained topology and corresponding band diagram.

Figure 6 .
Figure 6.Optimized solution and corresponding band diagram (the band gap should encompass 400-600 Hz) obtained with a mesh resolution of 40 × 40.

Figure 8 .
Figure 8. Optimized solution and corresponding band diagram (the band gap should encompass 400-600 Hz) obtained using four materials.

Figure 9 .
Figure 9. Iteration histories of the topology, objective, and volume fractions of the four-material case (the band gap should encompass 400-600 Hz).

Figure 10 .
Figure 10.Optimized solutions and corresponding band diagrams (the band gap should encompass 1000-2000 Hz) obtained using three or four materials: (a) three-material case; (b) four-material case.
), one can obtain β =I T ( M11 I 1 ü + M12 I 2 ü M21 I 1 ü + M22 I 2 ü The Map of the Subsystems to the Macroscale Appendix B.1.The Map of the Restricted Subsystem to the Macroscale Consider the generalized eigenvalue problem defined by Equation (14), the solution of which fulfills the mass normalization condition Nth natural vibration mode and ω * (N) µ

Table 1 .
Parameters of materials utilized to obtain LRAMs.

Table 2 .
The changes in the band gap and mass from three-material to four-material cases.

Table 3 .
The results obtained with different initial guesses.