Numerical Modelling and Optimization of Two-Dimensional Phononic Band Gaps in Elastic Metamaterials with Square Inclusions

Abstract

A numerical simulation study on elastic wave propagation of a phononic composite structure consisting of epoxy and tungsten carbide is presented for low-frequency elastic wave attenuation applications. The calculated dispersion curves of the epoxy/tungsten carbide composite show that the propagation of elastic waves is prohibited inside the periodic structure over a frequency range. To achieve a wide bandgap, the elastic composite structure can be optimized by changing its dimensions and arrangement, including size, number, and rotation angle of square inclusions. The simulation results show that increasing the number of inclusions and the filling fraction of the unit cell significantly broaden the phononic bandgap compared to other geometric tunings. Additionally, a nonmonotonic relationship between the bandwidth and filling fraction of the composite was found, and this relationship results from spacing among inclusions and inclusion sizes causing different effects on Bragg scatterings and localized resonances of elastic waves. Moreover, the calculated transmission spectra of the epoxy/tungsten carbide composite structure verify its low-frequency bandgap behavior.

1. Introduction

Elastic metamaterials are architecturally engineered structures that exhibit unusual properties to control the propagation of elastic waves through their structures [1,2]. Elastic metamaterials can be composed of periodically arranged phononic crystals (PnCs) that control elastic waves at a specific frequency range where the incoming wave cannot propagate—named stop bands or band gaps. The observation of the band gap results from two possible mechanisms: Bragg scattering and localized resonances of elastic waves [3]. The observation of the Bragg scattering results from the interaction of elastic waves scattered by periodic unit cells. The localized resonance is elastic waves locally oscillated around the center of the inclusion structure, mainly depending on the inclusion structure itself. Therefore, by tuning the geometry, size, arrangement, and shape of the scatterers within the soft host, it could be possible to achieve band gaps with a wider frequency range by combining these two mechanisms [4]. This is useful for a variety of applications such as noise and vibration isolation [5,6,7,8,9], energy harvesting [10,11,12,13], and acoustic and elastic filters and waveguides [14,15,16,17].

PnCs are typically designed by placing a high-stiffness material inclusion in a low-stiffness (soft) matrix [18]. Previous studies have reported that multiple bandgaps can be realized using different shapes of inclusions within a soft matrix. For example, Cheng and Shi achieved low-frequency band gaps using two- and three-component sphere inclusions [19]. Wang et al. achieved complete multiple band gaps using two-dimensional phononic crystal slabs with periodic cross-like holes [20]. Mizukami et al. achieved attenuation of broadband vibration using continuous carbon fiber-reinforced plastics with steel disks [21]. Elastic band gaps of PnCs could also be achieved by other different shapes of inclusions, such as comb-like inclusions [22], spindle shaped [23], quasi-Sierpinski carpet [24], and square/triangle shaped [25,26]. Among these studies, the geometric effect of the inclusions on the bandwidth and central frequency of the band gap are discussed by changing the size, inclusion number, and filling fraction of the unit cell. However, the effect of irregular size and locations of inclusions within a unit cell are rarely discussed but will be included in this work.

Another critical factor in determining the bandwidth and central frequencies of the bandgaps is the material properties of the inclusions and matrix used in a phononic crystal [27], which includes mass density ratio, the shear modulus ratio and Poisson’s ratios. Several studies have shown that the central frequencies of the band gap can generally be increased by increasing Poisson’s ratio of the inclusion/host material [28,29,30]. In addition, the increase in density and Young’s modulus ratios increase the bandwidth of the band gaps [31]. Recently, tungsten carbide inclusion and epoxy matrix has been used in different elastic metamaterials for broadband elastic attenuation applications owing to their high Young’s modulus ratios [32,33]. For example, Lu et al. have demonstrated that a phononic crystal composed of tungsten carbide and epoxy can achieve wider broadband phononic band gaps [34]. Li et al. achieved wide elastic bandgaps using the optimized structure of tungsten carbide spheres embraced in an epoxy matrix [35].

To realize a broad elastic bandgap in the low-frequency range, this work presents a numerical study of a periodic elastic composite composed of tungsten carbide within the epoxy host. The effects of square inclusion size within the unit cell (filling fraction), the stiffness of square inclusion, number of square inclusions, rotation angles of square inclusions, graded inclusion size, and irregular inclusion locations on their bandgap are studied systematically. A relationship between the resulting band gaps and the geometric dimensions is provided. The mechanism of the broad phononic bandgap is analyzed by solving its eigenmodes near the phononic bandgap. Two types of mechanisms, including Bragg scattering and localized resonances resulting in elastic bandgaps, are discussed. Moreover, analysis of the transmission loss also demonstrates the broad bandgap metrics of the proposed elastic metamaterials. This study is potentially beneficial as a guideline to design an elastic metamaterial with the broadband bandgap.

2. Numerical Methods

Elastic wave propagation through use of a phononic crystal was studied using a two-dimensional finite element method implemented in Solid Mechanics COMSOL modules [36]. The governing equation is given by Equation (1):

$$\rho \frac{{\partial }^{2}u}{\partial t^2}=f_v+\nabla ·\sigma .$$

(1)

where ρ and u are the density and displacement of the elastic materials, the volumetric force f_v is the force per deformed volume, and σ is the stress tensor. In this work, the proposed phononic crystals consist of high-stiffness inclusions placed in a low-stiffness matrix. The low-stiffness matrix in this study is epoxy, whereas the high-stiffness inclusions are composed of tungsten carbide [34]. The material properties are summarized in

Table 1. The material properties of the matrix and inclusion used in the study.

Material/Property	Elastic Modulus E (GPa)	Density ρ (kg/m3)	Poisson’s Ratio ν
Epoxy	4.344	1180	0.3679
Tungsten Carbide	387.6	13,800	0.3459

.

The dispersion relations describing the wave propagation behavior in a periodic phononic crystal was obtained by solving the eigenvalue problem with the periodic boundary condition. The Floquet–Bloch periodic boundary conditions were applied to the four edges of the unit cell to ensure a periodically repeating nature, as indicated by Equation (2).

$$u(\overrightarrow{r})=u_k(r)e^{j\overrightarrow{k}·\overrightarrow{r}}.$$

(2)

where k is the periodicity vector with a length |k| = 2π/Λ, in which Λ is the period in that direction and r is the position vector. Modal analysis was performed to obtain the eigenfrequency corresponding to a given Bloch wave vector k = (k_x, k_y). The study was performed on a range of wave vectors covering the irreducible Brillouin zone (IBZ), which spans from Γ to X, X to M, and M to Γ, as shown in Figure 1a. The PnC is infinite and periodic in the x, y directions with a lattice constant a_o = 1 m. For the elastic wave transmission simulation, the proposed epoxy/tungsten carbide composite unit cells were centered in the simulation domain with a 1 m distance away from the absorption boundary. An incident force of 1 N was applied in front of the composite structure in the propagating direction, and edges of the simulation domain were applied with the perfectly matched layers [37], which absorb all outgoing wave energy. The transmission loss as a function of frequency is given by

$$TL(dB)=10log\left(P_{in}/P_{out}\right).$$

(3)

where P is the power derived from mechanical energy flux, and P_in and P_out are the power of the incident and transmitted elastic waves.

3. Results and Discussions

3.1. Elastic Composites with Square Inclusions and Epoxy Host

i.

Single-square inclusion/epoxy host with different Young’s moduli.

Firstly, the unit cell composed of a single-square inclusion within a soft host was studied by varying its filling fraction (FF). The filling fraction was defined as FF = (a_i)²/(a_o)², where a_i is the side length of the square inclusion and a_o is the unit cell length, as shown in Figure 1b. The length a_o was kept constant at 1 m, and the a_i was varied from 0.316 to 0.948 m, corresponding to different FFs from 0.1 to 0.9. The effect of the stiffness of the inclusion in the matrix was studied by changing its Young’s modulus ratio (YMR), defined as YMR = E_i/E_o, where E_i is Young’s modulus of the inclusion material and E_o is Young’s modulus of the host material. The material considered here for the matrix is epoxy with E_o = 4.3438 GPa, and Young’s modulus of the inclusion E_i was varied to yield YMRs of 9, 90, and 300. Figure 1c show the dispersion curves of phononic crystals with the single-square inclusion at FF = 0.5, where the shaded region represents the bandgap and points A and B are located at the maximum and minimum frequencies of the bandgap, respectively.

Figure 1d shows that in the case of YMR = 300, the bandwidth increases as the filling fraction increases, reaching a maximum of 3280 Hz at FF = 0.8, after which it starts to decrease with any further increase in FF. When YMR was 9, the elastic composite bandwidth showed a slight change with an increased filling factor value, which results from the lack of the significant contrast in the elastic stiffness of the inclusion to that of the matrix. In general, a higher YMR of the elastic composite has a broader bandgap, and its bandwidth increases and reaches a maximum point at a specific FF. Then, the bandwidth would decrease with increase in its FF. The relationship between bandwidth and filling fractions are strongly related owing to the corresponding changes in geometric dimensions in the periodic unit cell, which affects the resonance frequency of the elastic eigenmodes [26,27,28,29,30,31,32,33,34,35,36,37,38]. In Figure A1a,b of Appendix A, the lower bound frequency (LBF) and upper bound frequency (UBF) are plotted against the filling fraction of the square inclusion at YMR = 9, 90, and 300. The simulation result shows that the relationship between UBF and filling fraction followed a similar trend to that of the bandwidth in Figure 1d, which indicates that the variations in bandwidth are more dominated by UBF, and this difference between LBF and UBF would be discussed based on the eigenmode analysis in the part (iii) of Section 3.1.

ii.

Increasing number of square inclusions within a unit cell.

To target a wider bandgap, in the following sections, the YMR of the inclusion/matrix was fixed at 90, corresponding to the ratio of tungsten carbide to the epoxy. The phononic crystals discussed above consisted only of one inclusion in each unit cell. Here, the number of square inclusions within the unit cell was increased to four and nine to investigate the effect of the inclusion number within the unit cell on the phononic band structure. The unit cells were designed as shown in the inset of Figure 2a,b, where the filling fraction was calculated as FF = 4(a_i)²/(a_o)² in the case of four-square inclusions within the unit cell and FF = 9(a_i)²/(a_o)² in the case of nine-square inclusions within the unit cell. The relationship between bandgap bandwidth and their filling fractions of three different unit cells was studied by changing a_i while fixing the a_o at 1 m.

The dispersion curves of phononic crystals with the four- and nine-square inclusions are shown in Figure 2a,b, where the shaded region represents the bandgap. The simulation results show that increasing the number of inclusions not only increased the bandwidth but resulted in an upward shift in the frequency range. The bandwidth was increased from 1316 Hz for one inclusion to 2678 Hz and 3915 Hz for four and nine inclusions. Interestingly, the LBF and UBF show a nonmonotonic relationship with the filling fraction, as shown in Figure 2c,d. This relationship causes the bandwidth of the bandgap to reach a maximum at a filling faction of 0.6.

iii.

Eigenmode analysis.

To analyze what causes the LBF and UBF of the tungsten carbide epoxy phononic crystal to not follow the same trend with filling fractions, their deformation distributions of eigenmodes at the lowest and highest frequencies of the bandgap were calculated, as shown in Figure 3. The simulation results show that at the lower frequency near the bandgap, the localized elastic wave modes are excited and the inclusion/epoxy vibrations behave as mass-spring oscillators [39], as indicated in Figure 3b,c. The phononic crystal exhibits a bandgap at wavelengths less than gap width among inclusions and inclusion sizes [4]. The nonmonotonic shift of the LBF results from two different trend of inclusion sizes and the width among inclusions with increase in the filling fraction. Therefore, LBF has a minimal value when the filling fraction is near 0.4 for a different number of inclusions within a periodic unit cell, as shown in Figure 2c. Meanwhile, according to a previous study conducted by Cui et al. [40], local resonances mainly depend on the resonator natural frequency and total mass density ratio. Therefore, the YMR does not have a pronounced effect on the LBF, as shown in Figure A1a.

At the UBF near the bandgap, elastic Bragg scattering effects were observed between the square inclusion and matrix at a higher frequency, as shown in Figure 3d,e. The Bragg scattering phenomena were induced by destructive interferences when the wavelength of the elastic wave was comparable to the periodicity of the scatterers embedded in the host elastic medium [41,42]. The first band gap was opened at wavelengths larger than the wavelength allowed by the Bragg mechanism, which is also supported by Kushwaha’s study [43]. Therefore, the observed wide broad bandgap in this study resulted from the coupling between two types of mechanisms [4]. The introduction of additional inclusions at a constant FF, as discussed in the previous section, reduced the periodicity of the inclusions, gap width among inclusions, and inclusion sizes. Thus, the reduction in these dimensions led to higher UBF and LBF, as shown in Figure 2c–e.

3.2. Effect of Rotation Angle of Square Inclusions

The effect of the orientation of the square inclusion was investigated by varying the rotation angle of the square inclusions from θ = 0°, 20°, 45°, 70° to 90° for different inclusion numbers, as shown in Figure 4a. The bandwidths of the bandgap in the phononic crystal with different inclusion numbers were plotted against the rotation angle θ for different filling fractions from 0.1 to 0.49 in Figure 4b,d. Due to the square shape of the inclusion, the simulation results show that the bandwidth versus rotation angle is symmetrical, as compared to the results between 0° to 45° and 45° to 90°. For the FF < 0.49, the bandwidth of the phononic bandgap shows a negligible change with an increase in the rotation angle of inclusions. Except when the rotation angle is 45°, the bandwidths of the bandgaps of the phononic crystal are increased owing to the reduced periodicity of square inclusions.

When the FF is near 0.49, the bandwidth significantly decreases at a rotation angle of 45° due to extremely narrow gaps among inclusions. In Figure A2a,b of Appendix A, the LBF and UBF versus rotation angles of the four-square inclusion composite, for different filling fractions from 0.4 to 0.45, are presented. The simulation results show that the bandwidth of the bandgap is decreased with a decrease in the FF, mainly due to the reduction in the UBFs. The elastic vibration at the UBF is related to the Bragg scattering [43], as indicated by the eigenmode analysis in the previous section. When the spacings among square inclusions are reduced with an increasing FF, the eigenmode caused by the Bragg scattering shifts to a lower frequency (higher wavelength).

3.3. Effect of Size of Square Inclusions

The effect of using squares of different side lengths within a single unit cell was investigated. In the four-square inclusions, the inclusion sizes were controlled by the square side lengths a_i1 and a_i2, and the filling fraction was calculated as FF = 2(a_i1² + a_i2²)/a₀. In the case of nine squares, three square side lengths were defined as a_i1, a_i2, and a_i3. The filling fraction was FF = 3(a_i1² + a_i2² + a_i3²)/a₀. In each case, the side lengths were defined, such as a_i1 < a_i2 < a_i3, to achieve the gradient effect, as shown in Figure 5a,b.

In the case of four inclusions, phononic crystals with five different filling factions were studied, the square dimensions of which are listed in

Table 2. The dimensions of the phononic crystals with 4 gradient square inclusions.

#	FF	ai1	ai2
1	0.1	0.1	0.2
2	0.2	0.15	0.2783
3	0.3	0.2	0.3316
4	0.4	0.25	0.3708
5	0.5	0.3	0.4

. Similarly, in the case of nine-square inclusions, phononic crystals with four different filling factions are presented in

Table 3. The dimensions of the phononic crystals with 9 gradient square inclusions.

#	FF	ai1	ai2
1	0.105	0.05	0.1
2	0.2175	0.1	0.15
3	0.375	0.15	0.2
4	0.5775	0.2	0.25

. The lower bound frequency (LBF), upper bound frequency (UBF), and bandwidth versus filling fraction are shown in Figure A3a,d. The simulation result shows that for four- and nine-square inclusions the gradient dimensions of inclusions lead to a higher LBF and a lower UBF, which indicates a narrower bandgap as compared to the equally sized inclusions, as shown in Figure 5c,d. This can be explained by the effective size of the phononic crystal with square inclusions of graded sizes. For example, when the FF = 0.3, the effective size of graded square inclusions was 0.2658 m, which is less than the effective size of normal inclusions, 0.2738. The larger bandwidth is due to the larger effective inclusion size, as indicated by the study of bandwidth versus filling fraction in Section 3.1.

3.4. Effect of Arrangement of Square Inclusions

In addition to the effect of number, sizing, and orientation of inclusions on the band gaps generated by the phononic crystal, the effect of the arrangement of the inclusions was investigated. A phononic crystal with four-square inclusions and a filling fraction of 0.2 is considered in this section. The location of the inclusions was varied by changing the distance r from the center of the unit cell to the center of the inclusion, as shown in Figure 6a. Figure 6b,d show that the locations of the inclusions significantly affect bandwidth of the bandgap. When the r is greater or less than the 0.3535 m, corresponding to the r in Figure 6c, the bandwidth is reduced. This study indicates that the gap width among inclusions is critical for controlling the bandwidth. A similar conclusion can be found in Section 3.1 for phononic crystals with high filling fractions.

3.5. Transmission Loss in Elastic Metamaterial

Transmission loss analysis was performed on an elastic metamaterial with different numbers of unit cells to verify the effectiveness of the metamaterial in preventing the transmission of the elastic wave through the material within the bandgap frequencies. Figure 7a shows the transmission loss spectra with one-, three-, and five-unit cells with a filling fraction of 0.125 for each. The peak transmission loss occurred at the bandgap regions, where the metamaterial prohibits wave propagation. It was verified that increasing the number of unit cells results in higher transmission loss and more defined peaks, enabling more effective prohibition of elastic wave propagation through the material.

The transmission spectra were obtained for normal and oblique incidences of the incoming wave to verify the omnidirectional absorption of the metamaterial, shown in Figure 7b. It was shown that the metamaterial is more effective towards normal incident waves than oblique (45°) incident waves, where the transmission loss peak reduced from −140 to −100 dB. The corresponding deformation distributions of the metamaterials with normal and oblique incidences are at 1700 Hz, as shown in Figure 7c, where the peak transmission loss occurs. It became evident that the five-array metamaterial is active in terms of prohibiting any wave propagation through the material, resulting in zero transmittance in the specified frequency range. By tuning the geometric and material parameters within the PnC, the elastic metamaterial proposed could be utilized in a variety of passive vibration isolation applications such as in small electronic devices, sensitive lab equipment, or even in civil engineering structures. The simplicity of the design, relying solely on square inclusions, enables the ease of manufacturing of this metamaterial to be implemented in a suitable application.

4. Conclusions

The band gap generation in phononic crystals, consisting of square inclusions in the soft host, with various sizes, filling fractions, material parameters, and arrangements was investigated in this work; in general, it can be concluded that a wide bandgap of the phononic crystal can be formed by combining two mechanisms: high-frequency Bragg scattering and low-frequency local resonances. Bragg scattering results from multiple elastic wave scatterings of periodic inclusions and local resonances results from the elastic natural resonant frequencies, which is more dependent on the inclusion size and gap width among inclusions. A nonmonotonic relationship between the bandwidth and filling fraction of the composite results from different responses between geometric dimensions and these two mechanisms. Tuning of band gaps to a larger bandwidth can be followed by the following procedures:

• increasing the number of inclusions within the periodic unit cell to increase the bandwidth and cause a blue shift in the central bandgap frequency;
• increasing the filling fraction of the unit cell to increase the bandwidth, and maximum bandwidth is reached near FF = 0.6 when using a tungsten-carbide/epoxy phononic crystal;
• rotating the angle of inclusions causes a minor effect on the bandwidth of the bandgap, except when the spacing among inclusions is reduced to near or less than the vibration wavelength;
• inclusion with a graded size leads to lower bandwidth, as compared to that of equally sized inclusions in a phononic crystal;
• a narrow spacing between the inclusions drastically reduces the bandwidth and UBF, and this can be found in the study of the rotation and arrangement effect;
• simulate transmission spectra is a way to verify the near-zero transmittance of the proposed elastic metamaterials over the frequency range within the phononic bandgap.