## Abstract

A finite-difference time-domain (FDTD) simulation coupled with an immersed-boundary method is used to investigate sound attenuation through both two-dimensional (2D) and three-dimensional (3D) cylinder arrays. The focus is on sound attenuation behaviors near Bragg’s bandgap frequencies for periodic structures. Both 2D and 3D simulations show that the finite cylinder arrays produce significant sound attenuation near the bandgap frequencies, with more attenuation effects in the 2D cylinder arrays because of the uniformity of sound source and neglected structure diffraction in the third dimension. When extended to 3D simulation, which can accommodate physically realistic conditions, sound attenuation near Bragg’s frequencies is reduced in comparison with 2D results. The 3D simulation also reaches a better agreement when comparing with the measurement data from the literature. Results and discussions on arrangement of cylinder arrays to achieve better sound attenuation effects are also presented.

## 1 Introduction

Sound propagation over sonic crystals has attracted much interest and been studied during the past few decades. In the late 1980s, several studies showed that a transparent material could become opaque for any light wave vector if a strong modulation of the refractive index in the three dimensions of the space is attained [1]. In the 1990s, it was demonstrated that sculptures built by the periodic arrays of cylinders inhibited the sound transmission for certain frequency ranges related to this modulation just as photonic crystals do with light [2]. Thus, the sonic crystals are defined as a periodic arrangement of structures made of sound hard scatters, and sound can be attenuated in a certain range of frequency [3].

For an infinitely periodic structure, there is a range of frequency known as bandgap or stop-band. The sound attenuation level at the bandgap range usually increases. The phenomenon follows Bragg’s diffraction law and those bandgap frequencies are theoretically predictable. The Bragg’s scattering effect is usually used along with the resonators to maximize the sound blockage effect and widen the bandgap range [4,5]. Therefore, periodic structures can be designed so that the interested acoustic signals can be enhanced or blocked. The application can be realized by systems such as Helmholtz resonators, periodic lattices insulators, or hybrid systems. The pattern of each unit of the periodic structure can be very complicated, especially with the recent development of the acoustic metamaterial and acoustic metasurface [4,6,7]. However, to demonstrate the influence of system-level arrangement instead of single unit pattern to the bandgap and Bragg’s scattering effect, we will use the simplest circular shape cylinders in our following study. Our previous two-dimensional (2D) study already included different ground conditions [8]. When the problem is extended to physically realistic three dimensions (3D), acoustic reflection and diffraction from the surrounding geometries and environmental conditions can be even more complicated. Therefore, in order to design the cylinder arrays to effectively attenuate sound waves, computational simulations are necessary and very helpful.

In addition, it is important to verify, in both 2D and 3D, how sound attenuation behaves near the theoretically predicted bandgap frequencies for practical, finite periodic structures. Jiang et al. [9] studied bandgaps of 2D periodic tube arrays in various angles. They reported that the second bandgap varies with different array geometries, while the first bandgap still matches Bragg’s law. Montiel et al. [5] did an analytical and numerical study comparing a 5 × 51 cylindrical with C-shape scatters. The slit-chamber feature allows the scatters to perform as resonators, and therefore, the system bandgap frequencies do not simply follow Bragg’s scattering law. Recently, Yu et al. [4] presented a design based on the finite element method (FEM) that acoustic metasurface can be constructed in chamber shape. A series of sub-chamber perform as a quarter-wave tube, which provides not only noise insulation but also ventilation and heat conduction. All these studies show that numerical simulation is an effective way to design and analyze acoustic bandgaps, which also motivated this study.

So far, plane wave expansion analysis [10–13], multiple-scattering theory (MST) method [14,15], boundary element methods, and FEMs [4,16–20] are widely used in solving sound propagation around objects. Finite-difference time-domain (FDTD) methods, as a type of time-domain techniques, can simulate the time evolution of acoustic pressure propagation. When acoustic waves are propagating around complicate geometries, the FDTD simulation can capture the development and interactions among incident, reflection, and diffraction waves, which is difficult to observe by using frequency-domain techniques. Particularly, when the geometries of the propagation problem are realistic and complicated, many times these problems cannot be simplified to be two-dimensional. Therefore, we need to use a 3D FDTD method to solve the problem. With the advent of high-performance computers, FDTD methods have evolved to be a powerful and effective way for simulating sound propagation around complex geometries, different media, and moving objects. By using parallel computation techniques, the time consumption for simulation can also be significantly reduced [21]. The method has already been implemented in a 2D case [8].

The purpose of this study is three folds. First, we validate our implementation of 3D linearized Euler equation (LEE) solver by comparing simulation results with the measurement results from an existing literature [3]. Second, to demonstrate the differences between 2D and 3D simulation and emphasize the importance of 3D simulation, we conduct several studies based on cylinder arrays and discuss the source and distance effects. Third, the same cylinder array simulations are used to study the Bragg’s diffraction law in 3D and discuss the sound blockage effects due to different cylinder array configurations.

In this study, the linearized Euler equations are employed for sound propagation in air, and the Zwikker–Kosten equations [22] are used for sound propagation inside cylinder arrays that are assumed to be made of porous materials. Different arrangements of cylinder arrays are studied for sound attenuation in both 2D and 3D.

The paper is organized as follows: In Sec. 2, governing equations for acoustic propagation and the immersed-boundary method are described. Then, numerical simulation examples including 2D and 3D simulations are studied along with the comparison with measurement data. Finally, we discuss the results and offer the conclusion.

## 2 Numerical Method and Model Description

*f*

_{u}and

*f*

_{p}are fictitious body forces to enforce the velocity and pressure to accommodate the governing equations in air and inside a rigid or porous object, which are given by [2]

*u*_{av},

*p*

_{av}, and

*α*

_{av}are the time-averaged velocity, pressure, and specific volume, respectively;

**,**

*u**p*, and

*α*are the corresponding acoustic fluctuations, with

*α*given as

*γ*is the ratio of specific heat. The coefficients, Ω,

*c*

_{s}, and

*σ*in Eqs. (3) and (4), are porous medium porosity, structure factor, and flow resistivity, respectively. A second-order finite-difference numerical scheme both in time and space [27] is applied for the simulations in this study. It should be noted that a good agreement has been reached between our simulated results and MST in our previous work and details can be found in Ref. [8].

## 3 Results and Discussion

### 3.1 Sound Propagation Simulation in Two-Dimensional Cylinder Arrays.

To compare differences between 2D and 3D simulations and investigate the bandgap phenomenon, regular lattice cylinder array arrangements are used in the study. Figure 1 shows the geometrical setup of the 4 row- and 3 column- (4 × 3) simulation, with each of the cylinder diameter at 2 m. The computational domain is 10 m in the *x*-direction and 15 m in the *y*-direction. A uniform grid mesh is used in the simulation, with the grid size determined to ensure that there are at least 20 grid points within one shortest wavelength for the frequency range in the simulation. Meanwhile, the Courant–Friedrichs–Lewy number is set to be 0.3 to ensure stable and efficient simulation. Details of the scheme and its verification can be found in Ref. [27]. The simulation time is 65 ms, which is long enough to allow all the reflection and diffraction waves to fully pass all the receivers.

*x*

_{0},

*y*

_{0}) and any location (

*x*,

*y*) in the domain. The round circles in Fig. 1 represent the cylinder array. The horizontal distance,

*H*

_{y}, the vertical distance,

*H*

_{z}, and the square lattice distance,

*L*, are set to be 5 m, 3.5 m and 1 m, respectively. The diameter of each cylinder is 0.5 m. Porosity Ω, structure factor

*c*

_{s}, and flow resistivity

*σ*in Eqs. (3) and (4) are chosen to be 0.3, 3, and 4 × 10

^{8}Pa × s × m

^{−2}, respectively, which follows our previous study and represents the acoustic rigid material of the cylinders. In order to study the relation between bandgap and cylinder array layout, four different configurations, 3 × 3, 4 × 3, 5 × 3, and 4 × 4, are simulated. Four receivers are put behind the cylinder arrays in the centerline along the

*y*-direction at 9 m, 10 m, 11 m, and 12 m to monitor acoustic pressure at these locations. The perfectly matched-layer (PML) boundary conditions [28] are specified around the four boundaries of the simulation domain. The PML thickness is equally set to be 1 m to ensure a non-reflection condition. In Fig. 1, the PML on the left side represented by dot-dashed line means the left area will stay as air until the whole wave-front passes the line at

*y*= 1 m. Once the wave-front completely passes the dot-dashed line, the left area will be changed to the PML, so that the waves reflected back from the cylinders will not enter the computational domain.

Figure 2 shows pressure contours of four different configurations of cylinder arrays at time 25 ms when the wave-front almost passes the cylinders. The complicated wave structures due to multiple reflections and diffractions are observed. Figure 2(d) clearly shows diffraction effects at the wave-front. Acoustic pressures near the wave-front become weaker when the number of cylinders increases, by looking at the wave-front patterns of the contour plots in Fig. 2. The increase of cylinders in the propagation direction has more significant effects (comparing Fig. 2(a) for the 3 × 3 case with Fig. 2(d) for the 4 × 4 case) than the increase of cylinders in the perpendicular direction to propagation (Figs. 2(b) and 2(c) for the 4 × 3 and 5 × 3 cases, respectively). These phenomena also show in the 3D cases to be discussed later.

_{test case}and PSD

_{free space}are the acoustic pressure power spectrum densities for the cylinder array case and the corresponding free-space case, respectively. Figure 3 gives plots of sound attenuations of the four different receiver locations in each of the test cases. Comparing the SA readings in Fig. 3, we can find that the first SA peak is always around 170 Hz, the first dash-dot line in the figures, which can be explained by Bragg’s diffraction law. The theoretical Bragg bandgap frequency

*f*is given as

*n*is an integer number,

*c*is the speed of sound (340 m/s), and

*L*is the sonic crystals lattice constant (1 m). The resultant first bandgap value is 170 Hz, which matches the numerical results when

*n*= 1. The second and third peaks can also be found at 340 Hz (

*n*= 2) and 510 Hz (

*n*= 3) in most of the test cases, although they are not as obvious as the first peak in some cases.

The SA curves in Fig. 3 also show that if more cylinder rows are added on top and bottom of 3 × 3 and 4 × 3 arrays (from 3 × 3 to 5 × 3), more SA fluctuations are observed in Fig. 3(c), especially in the mid and high-frequency ranges. Meanwhile, the first peak SA value is decreasing. When more cylinder array rows are added, stronger reflection and diffraction effects are observed at locations close to the cylinder array, resulting in reduced noise blockage, particularly at the low-frequency range.

If one more column of cylinders is added on right-hand side of 3 × 3 and 4 × 3 arrays (from 3 × 3 to 3 × 4 and from 4 × 3 to 4 × 4), we still observe the first peak at the same frequency. However, the magnitude of the peak decreases more. Moreover, Fig. 3(d) shows that the best low-frequency sound attenuation location moves away from the cylinder array, from 10 m to 12 m, which suggests that noise at locations close to the cylinder array may not be properly blocked, especially at the low-frequency range.

### 3.2 Sound Propagation Simulation in Three-Dimensional Cylinder Array.

To compare the differences between 2D and 3D results, most of the test setup in 2D test cases are carried over to the 3D simulation. Figure 4 presents the simulation setup in 3D simulation from top, front, and isometric views for a 4 × 3 test case. The computational domain in the *z*-direction is 4 m, from −1 m to 3 m, where *z* = 0 is the ground. The same type of Gaussian point source, as specified in Eq. (6), is placed in the middle of the domain at (5 m, 0 m, 1 m) with $r=(x\u2212x0)2+(y\u2212y0)2+(z\u2212z0)2$. The cylinder array is 2 m high in the *z*-direction. The horizontal distance, *H*_{y}, the vertical distance, *H*_{z}, and the square lattice distance, *L*, are set to be 5 m, 3.5 m, and 1 m, respectively, to ensure a consistent setup with the 2D simulation. Also, the material of the cylinders remains as acoustic rigid material. Four-cylinder array combinations, 3 × 3, 4 × 3, 5 × 3, and 4 × 4, are tested. Four receiver locations are equally distributed along the *y*-direction in the middle of the domain from (5 m, 9 m, 1 m) to (5 m, 12 m, 1 m). PMLs with 1-m thickness are specified around the simulation domain in all directions. The left side PML represented by dot-dashed line is only active when the wave-front completely passes this line.

Although the simulation setup is very similar between the 2D and 3D cases, two major differences are noted: one is the 2D point source (cylindrical wave) versus the 3D point source (spherical wave); the other is the third dimension (the *z*-dimension) propagation. Since the sound attenuation is calculated based on the free-space solution in each of the corresponding 2D and 3D cases, the relative influence of the 2D and 3D sources is removed. In addition, because the top and bottom boundary conditions are specified as PMLs, which mimic a none-reflection free-space condition, the domain size in the *z*-direction should not influence the final conclusion, assuming the PMLs are working correctly. We have carefully tested the PML boundary-condition implementation and made sure that the boundary conditions in the third dimension are accurate. As for the size of the cylinder in the *z*-dimension, it can possibly influence the blockage, as the effective geometry for blockage is reduced if the cylinder length (the size in the *z*-dimension) is reduced. This is actually the reason why we select a relatively smaller size of the cylinder length of 2 m versus the cylinder diameter also at 2 m to show the effect. Any longer length would be in between the 3D blockage effect shown here and the 2D blockage effect. Any shorter length would change the geometry that can no longer be considered as a cylinder because it would be more like a disc.

Figure 5(a) shows part of the pressure contours of the 3D simulation at the sections of *x* = 5 m and *z* = 1 m at the simulation moment *t* = 25 ms. Comparing it with Fig. 2(a), the min–max contour levels of the 3D simulation are much lower than that of 2D simulation. This is because acoustic pressures in 3D decay much faster than that of 2D waves, e.g., at the rate of 1/*r* in 3D versus $1/\u221ar$ in 2D. Hence, the lower magnitude of a 3D contour level is expected. The 3D point source versus 2D point source also changes wave interaction patterns. However, the trend that the increase of cylinders causes weaker wave-front is the same as in the 2D cases. In addition, the increase of cylinder in the propagation direction has more significant effects (comparing Fig. 5(a) for the 3 × 3 case with Fig. 5(d) for the 4 × 4 case) than the increase of cylinders in the perpendicular direction to propagation (Figs. 5(b) and 5(c) for the 4 × 3 and 5 × 3 cases, respectively). These behaviors are the same as those in the 2D cases.

However, the SA levels between the 2D and 3D cases change. Figure 6 gives the SA of the different cylinder array arrangements. It should be noted that the SA here is calculated with respect to a 3D free-space simulation result. Therefore, the source power influence is removed from the calculation. Comparing the SA between 2D and 3D in Figs. 3 and 6, we can find the following two facts. (1) The SA curves are very similar in shape, while the biggest difference is the magnitude. For example, in 3 × 3 configuration, the SA of first peak in Fig. 3(a) is 100 dB. However, SA of first peak in Fig. 6(a) can only reach 40 dB. The counterpart of a point source in the 3D cases is a line source in the 2D cases, and the finite-length 3D cylinders in the 3D cases were represented by infinitely long 2D cylinders. These two factors introduce stronger diffraction effects in the third direction (the *z*-direction), which can significantly reduce the blockage performance of the cylinder arrays. (2) The first two peak frequencies always appear at the bandgap frequencies around 170 Hz and 340 Hz in Fig. 6 which implies that Bragg’s diffraction law is still valid in 3D wave propagation problems.

Similarly, when one more row is added on each side of the 3 × 3 array to become a 5 × 3 array, the simulation results in Fig. 6(c) show that SA of the second bandgap has obvious changes, and a similar phenomenon can also be found at the third or the fourth bandgap in higher frequency ranges. When cylinder row numbers grow, more obvious bandgaps can be obtained. Jiang et al. [9] did a similar study but with a 15-deg angle between the incident wave and the cylinder array, and their experiment showed very similar results.

Furthermore, when one column is added based on the 4 × 3 array to become a 4 × 4 array, SA at the first bandgap in Fig. 6(d) reduces significantly which agrees with the 2D study results in Fig. 3(d). Meanwhile, the best attenuated frequency and location both shift toward the high-frequency range, suggesting a thicker cylinder array can enhance the sound blockage performance in the mid and high-frequency ranges and push the best blockage location away from the cylinder array.

### 3.3 Sound Propagation Simulation Results in Comparison With the Measurement Data for the Bandgap Behaviors.

*y*-direction, 4.25 cm in the

*x*-direction, and 25 cm in the

*z*-direction. This setup is the same as the test section used in Gupta’s experiment [3]. The sound source is specified at the left side of the domain from

*x*= 0 cm to

*x*= 3.8 cm and from

*z*= 6.25 cm to

*z*= 18.75 cm to mimic the FWB-MX970RS speaker used in the test. A linearly combined sinusoidal plane wave in the form of Eq. (9) is used in the study.

*f*is the frequency of the plane wave. The plane wave propagates through five cylinders. The diameter of the cylinder is 3 cm, and the height is 25 cm. The space between each cylinder is 8.5 cm. A receiver is placed near the end of the domain at (2.125 cm, 124 cm, 12.5 cm), 10 cm away from the last cylinder center. The PMLs are used at top, bottom, and right boundaries of the domain to ensure a plane wave like condition and mimic acoustic absorbing foam used in the experiment. The 2D simulation is just one slice of the 3D simulation in

*x–y*plane, with the setup in Fig. 7(a). A uniform mesh size, 0.5 mm, is used in computation. The total simulation time duration is 25 ms.

Figure 8 gives 2D and 3D simulation results of acoustic pressure contours at *t* = 7.5 ms. In the 3D cases shown in Figs. 8(b) and 8(c), it is clear that sound pressure is strengthened around the first cylinder and then gradually weakened in the wave propagation direction. In the 2D case in Fig. 8(a), the level of sound pressure does not seem to be reduced. However, it should be noted that the pressure fluctuation range is reduced. The propagation pattern actually shows that the pressure fluctuation magnitude is reduced downstream, fluctuating within 0–50 (the raw pressure level, not really dB) in the time history. When this data set is mapped to the frequency domain in terms of the power spectral density, the magnitude of the sound pressure is actually reduced. Later in Fig. 9 for the SA level, it actually shows that the 2D SA levels are higher than the 3D SA levels, consistent with the results shown previously.

Figure 9 shows the SA comparisons among 2D and 3D simulation results and measurement data from Gupta et al. [3]. All of the three curves in Fig. 9 demonstrate that the bandgap frequency ranges are between 1500 Hz and 2500 Hz, and 3500 Hz and 4500 Hz, which match the bandgap frequency ranges of the measurement data. However, the 3D simulation results agree with the measurement much better than those of 2D simulation. The plane wave source and the infinitely long cylinders in the 2D simulation are unable to represent physically realistic sound wave propagation in the measurement which is 3D. Therefore, 3D simulation results agree with measurement data better. In addition, the good agreement of sound attenuation frequency ranges around the theoretical bandgap frequencies provides additional evidence that validates the simulation results.

It should be noted that the purpose of the comparison is actually not to match the quantitative values with the measurement data at every frequency, because it is not practical for the simulation to specify exactly the same source configuration as the measurement, even in the 3D simulation. Rather, it is to compare the general trend, particularly to observe if Bragg’s bandgap indicated by the measurement is also captured in the simulation. In addition, whether 2D or 3D, or both, can capture the Bragg’s bandgap and which one can do better in comparison with the measure. In this sense, the 3D simulation shows a better agreement with the measurement results.

## 4 Conclusion

The sound blockage effects of cylinder arrays are studied for both two- and three-dimensional simulations. The validity of the finite-difference time-domain method incorporated with the immersed-boundary method is shown in an agreement between the simulated bandgap frequencies and the theoretical values.

In the two-dimensional simulations, when the number of cylinders in the vertical direction is fixed, the usage of more cylinders in the wave propagation direction results in better sound blockage effects. It is not the case when the receiver is close to the arrays because the sound pressure level is enhanced instead of attenuated due to reflections. When the number of cylinders in the wave propagation direction is fixed, by adding more rows in the vertical direction, more obvious bandgaps can be observed, which indicates that better sound attenuation can be achieved in the mid- and high-frequency ranges. Three-dimensional studies also show that more cylinders in the vertical direction may push the best sound attenuation location at the low bandgap frequency away from cylinder array. Therefore, selecting an appropriate number of cylinder barriers is an important factor for sound attenuation.

Additional simulation of cylinder arrays for sound attenuation is performed to compare with the measurement [3]. In comparison with the 2D simulation results, the 3D simulation results agree better with the experimental data, showing two frequency ranges of sound attenuation around the theoretical bandgap frequencies, again validating the simulation model.

## Acknowledgment

This research was partly supported by the U.S. Army under a cooperative agreement W911NF-14-2-0077.