Three-Dimensional Phononic Band Gap
Calculations Using the FDTD Method and a
PC Cluster System
Po-Feng Hsieh, Tsung-Tsong Wu, and Jia-Hong Sun
Abstract—This paper aims at studying the band gap
phe-nomena of three-dimensional phononic crystals using the finite difference time domain (FDTD) method and a PC cluster system. In the paper, Bloch’s theorem is applied to the wave equation and to the boundary conditions of the periodic structure. We calculate the variations of displace-ments and take discrete Fourier transform to acquire the resonances of the structures. Then, the dispersion relations of the bulk acoustic wave can be obtained and the band gaps are predicted accordingly. On the other hand, because of larger data calculation in three-dimensional phononic crys-tals, we introduce the PC cluster system and parallel FDTD programs written with respect to the architecture of a PC cluster system. Finally, we discuss the numerical calculation of two-dimensional and three-dimensional phononic crystals consisting of steel/epoxy and draw conclusions regarding the band gap phenomena between these phononic crystals.
I. Introduction
T
herehas been research pertaining to artificial crys-tals starting with photonic cryscrys-tals since the 1980s. A photonic crystal is a periodic arrangement of dielectric materials that an electromagnetic (EM) wave of a spe-cific frequency range is forbidden to pass through—the so-called band gap phenomenon of the photonic crystal. Owing to the analogy of EM waves and acoustic waves, the band gap phenomenon also exists in phononic crystals, the periodic arrangement of elastic materials. The band gap phenomenon may potentially find many engineering applications, such as in acoustic wave filters, the acous-tic waveguide, and acousacous-tic resonators, which enhance the necessity to develop a fast and efficient wave propagation theory for the phononic crystal.Theoretical analyses on the band gaps of phononic crys-tals have been conducted during the last 15 years. For re-search in band gaps of bulk acoustic wave (BAW) modes of phononic structures, the mixed and transverse polar-ization modes have been carried out using the plane wave expansion (PWE) method [1]–[7]. The PWE method is also adopted to study level repulsions of bulk acoustic waves in composite materials [8]. The phononic
struc-Manuscript received March 30, 2005; accepted July 6, 2005. The authors gratefully acknowledge the financial support of this research by the National Science Council of ROC (NSC 93-2212-E-002-025). The authors are with the Institute of Applied Mechan-ics, National Taiwan University, Taipei, Taiwan (e-mail: [email protected]).
tures considered in the aforementioned studies include two-dimensional (2D) solid/air, solid/solid, and solid/fluid structures. It is worth noting that in conducting the acous-tic waves in the solid/fluid structure, a modified PWE method was used [4]–[7]. In addition to the band gap anal-yses of bulk modes in 2D phononic structures, there is some literature on the investigation of band gaps of sur-face acoustic waves (SAW) [9]–[11]. More detailed stud-ies of temperature effects of surface wave band gaps and electromechanical coupling coefficient in a piezoelectric phononic crystal are also reported [12], [13]. Instead of the PWE method, the multiple scattering theory (MST) has been applied to study the band gaps of bulk wave proper-ties in three-dimensional (3D) periodic acoustic compos-ites and the band structure of a phononic crystal consist-ing of complex and frequency-dependent Lam´e coefficients [14]–[17]. With regard to the numerical simulation of band gaps in phononic structures, the finite difference time do-main (FDTD) method was applied to interpret the experi-mental data of 2D systems consisting of cylinders of fluids (Hg, air, and oil) inserted periodically in a finite slab of aluminum host [18]. It was also used to calculate the pe-riodic solid-solid, solid-liquid, and solid-vacuum compos-ites [19]. The coupling characteristics of localized phonons in photonic crystal fibers and analyses of mode coupling in joined parallel phononic crystal waveguides using the FDTD method are reported in [20], [21].
In these studies, most phononic crystals were 2D structures; some 3D phononic crystals are reported. The 3D structures include the simple cubic (SC), body-centered cubic (BCC), and face-body-centered cubic (FCC) lattices consisting of water/mercury (liquid-liquid) [22], mercury/aluminum (liquid-solid) [23], and steel/polyester (solid-solid) [24], [25]. In these 3D phononic crystals, the band gaps were evaluated by both the MST and the PWE methods. However, the MST method has a limitation of calculating with overlap scatters, and the PWE method encounters convergence problems when the phononic crys-tal has a large elastic mismatch. The FDTD method is very suitable for dealing with different geometric structures and for handling the convergence problem but requires consid-erable calculation.
The purpose of the present study was to develop a par-allel FDTD method and a PC cluster system to handle and accelerate the calculation of band gap phenomenon of the phononic crystal in 3D cases. By using a heteroge-neous FDTD formulation, the material property and the
geometric shape can be changed easily and the formula-tion can be employed to study the band gap phenomenon in different geometric shapes, material properties, and de-sign patterns. We adopted the message passing interface (MPI) to write the parallel FDTD program in order to accelerate the whole computational speed in the PC clus-ter system. Owing to the high-speed computation of the FDTD method and the PC cluster system, calculations of 3D phononic crystal can be achieved quickly.
II. The Finite Difference Time Domain (FDTD) Method
Whether the FDTD method or the PWE method is used, the BAW theory of phononic crystal is based on the same physical concept. In this paper, the periodic bound-ary using Bloch’s theorem is used to describe the peri-odic arrangement of the phononic crystal. The principle of FDTD is mentioned in this section.
A. Equation of the Wave Propagation
According to the theory of elasticity, all the elastic ma-terials are considered in the linear elastic range. Equation of motion and constitutive law with the anisotropic and inhomogeneous material can be expressed as [26]:
ρ¨ui= τij,j+ ρ fi, (1)
τij = Cijklεkl, (2)
where ρ is the density of materials, ui is the displacement,
τij is the stress, fi is body force, Cijkl is elastic constant,
and εklis strain. Eq. (1) and (2) describe the behavior of an
infinitesimal element of an anisotropic material, and they can apply to the inhomogeneous structure. Therefore, the FDTD method based on (1) and (2) is suitable to analyze the phononic crystal by arranging the density and elastic constant periodically.
In this FDTD method, the displacement and stress vari-ables in materials are used as the bases to describe the equation of motion and constitutive law, and further to develop the difference equations. To cooperate with the staggered grids adopted in the following difference equa-tion, the anisotropy of elastic material is limited. An or-thotropic material with nine independent elastic constants is allowed in the program. The matrix form of the elasticity constant using an abbreviated notation is
CIJ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C11C12C13 0 0 0 C21C22C23 0 0 0 C31C32C33 0 0 0 0 0 0 C44 0 0 0 0 0 0 C55 0 0 0 0 0 0 C66 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (3)
where the subscripts 1, 2, 3, 4, 5, and 6 of material con-stant C represent xx, yy, zz, yz, xz, and xy, respectively.
From (1), (2), and (3), we can derive the 3D wave equa-tion in which the displacement and stress are regarded as the variables; thus the nine equations of wave equations in matrix form are
⎡ ⎣uu¨¨12 ¨ u3 ⎤ ⎦ = 1 ρ ⎡ ⎣ττ1121 ττ1222ττ1323 τ31 τ32τ33 ⎤ ⎦ ⎡ ⎣∂∂12 ∂3 ⎤ ⎦ + ⎡ ⎣ff12 f3 ⎤ ⎦ , (4) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ τ11 τ22 τ33 τ23 τ13 τ12 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C11C12C13 0 0 0 C21C22C23 0 0 0 C31C32C33 0 0 0 0 0 0 C44 0 0 0 0 0 0 C55 0 0 0 0 0 0 C66 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂1u1 ∂2u2 ∂3u3 ∂2u3+ ∂3u2 ∂1u3+ ∂3u1 ∂1u2+ ∂2u1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (5)
In (4) and (5), the subscript numbers 1, 2, and 3 of dis-placement u and stress τ denote the components along the x, y, and z directions in the Cartesian coordinate, respec-tively, and the symbol ∂irepresents the partial differential
in direction i.
Then we use the second-order center difference equa-tion and the first-order center difference equaequa-tion to ap-proximate the second-order differential equation and the first-order differential equation, respectively. Therefore we define the first-order center difference equation in space first as: D1[M (i, j, k, l)] = 1 ∆x1 M i + 1 2, j, k, l − M i−1 2, j, k, l D2[M (i, j, k, l)] = 1 ∆x2 M i, j +1 2, k, l − M i, j−1 2, k, l D3[M (i, j, k, l)] = 1 ∆x3 M i, j, k +1 2, l − M i, j, k−1 2, l . (6)
In (6) and the following formulas, the i, j, and k inside the bracket represent the position index of the x, y, and z directions, respectively, and the symbol l is the index of time step. The ∆xi, where i = 1, 2, and 3, represents
the distance of two nearby discrete grids in the x, y, and z directions. By replacing the second-order differential in time with the second-order center difference in time, we have ∂2 ∂t2[M (i, j, k.l)] = 1 (∆t)2 M i, j, k, l + 1 2 − 2M(i, j, k, l) + M i, j, k, l−1 2 . (7) After the discrete procedure, we have the difference equa-tion of the displacement and stress (8)–(16) (see next two pages).
In (8)–(16), the symbol Ui and Tij are the discrete
U1 i +1 2, j, k, l + 1 2 = 2U1 1 + 1 2, j, k, l − U1 i +1 2, j, k, l− 1 2 + (∆t) 2 ρ i +1 2, j, k D1 T11 i +1 2, j, k, l + D2 T12 i +1 2, j, k, l + D3 T13 i +1 2, j, k, l + f1 i +1 2, j, k, l + 1 2 (8) U2 i, j +1 2, k, l + 1 2 = 2U2 i, j +1 2, k, l − U2 i, j +1 2, k, l− 1 2 + (∆t) 2 ρ i, j +1 2, k D1 T21 i, j +1 2, k, l + D2 T22 i, j +1 2, k, l + D3 T23 i, j +1 2, k, l + f2 i, j +1 2, k, l + 1 2 (9) U3 i, j, k +1 2, l + 1 2 = 2U3 i, j, k +1 2, l − U3 i, j, k +1 2, l− 1 2 + (∆t) 2 ρ i, j, k +1 2 D1 T31 i, j, k +1 2, l + D2 T32 i, j, k +1 2, l + D3 T33 i, j, k +1 2, l + f3 i, j, k +1 2, l + 1 2 (10) T11(i, j, k, l + 1) = C11(i, j, k)· D1 U1 i, j, k, l + 1 2 + C12(i, j, k)· D2 U2 i, j, k, l +1 2 + C13(i, j, k)· D3 U3 i, j, k, l +1 2 (11) T22(i, j, k, l + 1) = C21(i, j, k)· D1 U1 i, j, k, l + 1 2 + C22(i, j, k)· D2 U2 i, j, k, l +1 2 + C23(i, j, k)· D3 U3 i, j, k, l +1 2 (12) T33(i, j, k, l + 1) = C31(i, j, k)· D1 U1 i, j, k, l + 1 2 + C32(i, j, k)· D2 U2 i, j, k, l +1 2 + C33(i, j, k)· D3 U3 i, j, k, l +1 2 (13)
T23 i, j +1 2, k + 1 2, l + 1 = C44 i, j +1 2, k + 1 2 · D3 U2 i, j +1 2, k + 1 2, l + 1 2 + D2 U3 i, j +1 2, k + 1 2, l + 1 2 (14) T13 i +1 2, j, k + 1 2, l + 1 = C55 i +1 2, j, k + 1 2 · D3 U1 i +1 2, j, k + 1 2, l + 1 2 + D1 U3 i +1 2, j, k + 1 2, l + 1 2 (15) T12 i +1 2, j + 1 2, k, l + 1 = C66 i +1 2, j + 1 2, k · D2 U1 i +1 2, j + 1 2, k, l + 1 2 + D1 U2 i +1 2, j + 1 2, k, l + 1 2 (16)
Fig. 1. The unit of the grids in the FDTD method.
i, j = 1 ∼ 3. The positions of the staggered grids are de-fined by these discrete wave equations. Fig. 1 shows the unit of the staggered grids in the calculation of the FDTD method, and the relationship between the displacement grids and the stress grids.
B. Boundary and Initial Conditions
To calculate the dispersion curves of phononic crystals, the periodic boundary condition of FDTD is developed. Because of the periodicity of phononic crystals, we applied Bloch’s theorem to fulfill the boundary condition. The pe-riodic boundary condition satisfying Bloch’s theorem is promoted by Tanaka et al. in [19] for 2D phononic crys-tal cases. In the present article, we deal with 3D phononic crystals based on the same principle. According to Bloch’s theorem, the displacement and stress components in the periodic structure can be expressed in a periodic function as follows:
Ui(x, t) = eik·xUi(x, t), (17)
Tij(x, t) = eik·xTij(x, t), (18)
where k = (k1, k2, k3) is a wave vector, and Ui and Tijare
periodic functions which satisfy the following relation: Ui(x + a, t) = Ui(x, t), (19)
Tij(x + a, t) = Tij(x, t), (20)
where a is a lattice translation vector with the components a1, a2, and a3 of each direction.
In [19], the equation of motion (1) and constitutive law (2) are translated to solve the periodic functions Ui and
Tij and to apply the periodic boundary conditions (19)
and (20). Alternatively, we develop the FDTD in Ui and
Tij, and therefore the periodic boundary conditions are
derived as:
Ui(x + a, t) = eik·aUi(x, t), (21)
Tij(x + a, t) = eik·aTij(x, t). (22)
Expanding (21) and (22), we have the discrete forms of the surface normal to the X, Y, and Z axes:
Normal to the X axis:
U1 −1 2, j, k, l = eik1a1U 1 a1− 1 2, j, k, l U2 0, j + 1 2, k, l = eik1a1U 2 a1, j + 1 2, k, l U3 0, j, k + 1 2, l = eik1a1U 3 a1, j, k + 1 2, l (23) T11(0, j, k, l) = eik1a1T11(a1, j, k, l) T12 −1 2, j + 1 2, k, l = eik1a1T 12 a1− 1 2, j + 1 2, k, l T13 −1 2, j, k + 1 2, l = eik1a1T 13 a1− 1 2, j, k + 1 2, l .
Normal to the Y axis: U1 i + 1 2, 0, k, l = eik2a2U 1 i +1 2, a2, k, l U2 i,−1 2, k, l = eik2a2U 2 i, a2− 1 2, k, l U3 i, 0, k + 1 2, l = eik2a2U 3 i, a2, k + 1 2, l (24) T12 i +1 2,− 1 2, k, l = eik2a2T 12 i +1 2, a2− 1 2, k, l T22(i, 0, k, l) = eik2a2T22(i, a2, k, l) T23 i,−1 2, k + 1 2, l = eik2a2T 23 i, a2− 1 2, k + 1 2, l . Normal to the Z axis:
U1 i +1 2, j, 0, l = eik3a3U 1 i +1 2, j, a3, l U2 i, j +1 2, 0, l = eik3a3U 2 i, j +1 2, a3, l U3 i, j,−1 2, l = eik3a3U 3 i, j, a3− 1 2, l (25) T13 i + 1 2, j,− 1 2, l = eik3a3T 13 i +1 2, j, a3− 1 2, l T23 i, j +1 2,− 1 2, l = eik3a3T 23 i, j +1 2, a3− 1 2, l T33(i, j, 0, l) = eik3a3T33(i, j, a3, l) .
After setting up the periodic boundaries, we need to define the initial condition for the calculation. Because the periodic structure needs to be vibrated in this method, we provide a small disturbance in a random position of the structure for the setting of the initial condition, and we have the initial disturbance which is set in any grid of the unit cell
Ui(x1, x2, x3) = δx1,xδx2,yδx3,z, (26) where δ is a delta function and (x1, x2, x3) is the random
point which excites a wide-band signal in the unit cell at the initial time step l = 0. Thus, all possible wave modes are transported inside the phononic crystal under consider-ation, and then the signal is expanded into a time Fourier series; this permits the finding of the information about possible types of waves.
Before the calculation is started, the specific wave num-ber in the first Brillouin zone needs to be defined for the boundary condition. During the calculation, the variation of the displacement must be recorded while time increases. After obtaining a sufficiently large number of displacement data, we Fourier-transformed the variation of the displace-ment into the frequency space. Then the resonant frequen-cies of the periodic structure can be obtained. The reso-nant points are the existing peaks in the frequency spectra; they identified the eigenfrequencies of the phononic crystal for a given wave vector.
Fig. 2. The hardware architecture of a PC cluster system.
III. Parallel Computing in the PC Cluster System
The traditional supercomputer no longer satisfies the demands of large scientific computing in that the hard-ware has met the bottleneck of design and the price of a supercomputer is too expensive to be afforded by general research institutions. So the PC cluster system is becom-ing the most popular high-efficiency computation machine in scientific and technologic fields nowadays. In this study, because of the high-speed computation and low cost of the PC cluster system, the complex research of 3D phononic crystals can be calculated and discussed by saving much of the CPU computation time.
As Fig. 2 shows, the PC cluster consists of several PCs with different operating systems and different hardware ar-chitectures. PCs are connected with each other by external network devices, such as network cards based on Ethernet. Thus the PCs in the system can use the messages passing around to communicate with each other. Therefore the whole system of PCs can be regarded as a complete par-allel machine.
The parallel model of the PC cluster is SPMD (single program multiple data). In the model of SPMD, the paral-lelism is controlled by the programs, not the CPU instruc-tions. The controller of the system sends the program to every computational node involved in the calculation, and commands every PC to do the task. In the parallel com-puting of the PC cluster system, we adopt MPICH [27], a high-performance portable implementation of the stan-dard MPI, to parallel the program. MPICH offers many kinds of parallel ways in the library so that the program-mer can easily pass the message from one of the computers to the others by simply calling the function offered from the library.
In this paper, we accelerate the calculation speed by paralleling the program. Among several kinds of parallel
Fig. 3. Parallel model—working queue.
methods for calculation, such as dividing the calculation time or space, here we propose the method of dividing the calculation by wave vectors. Unlike dividing the cal-culation by time or space, which requires numbers of data exchanging, by calculating different wave vectors around the irreducible part of the first Brillouin zone, there are fewer messages passing in PC cluster system so that the computation can achieve higher performance. There is no message passing, data exchanging, and synchronization be-tween the calculations with different wave vectors. There-fore, dividing the computation by the wave vector is ex-pected to obtain higher performance in PC cluster system. In this parallel method, as Fig. 3 shows, the working queue is created to store all the wave vectors of the whole computation in the server before the calculation. Then the client sends the request to the server for the jobs asked for in the calculation. If the working queue in the server is not empty, the new job will be sent to the client. After receiving the reply from the server, the client can do the calculation with the specific wave vector. Then the client can iterate the job requesting process until the working queue is empty. Therefore the whole computation is fin-ished completely.
Using this kind of parallel method, the whole computa-tion can save more calculacomputa-tion time by increasing numbers of CPUs. Because the numbers of CPUs in calculation are increased, jobs dispatched to each node are also reduced. And in this method there is no message passing between the clients and fewer data exchanging between the client and the server. Thus the total waiting time is spent mostly in the calculation and just a little in the waiting for mes-sage passing. Therefore this parallel method is very suit-able for use in the FDTD study of phononic crystals.
Fig. 4. (a) A 2D phononic crystal with square lattice arrangement; (b) The irreducible part of the first Brillouin zone.
IV. Examples of Phononic Band GAP Calculations (2D and 3D)
In this section, the above-mentioned BAW theory of a phononic crystal using the FDTD method is employed to discuss phononic band gaps with different crystal struc-tures. The results for 2D and 3D phononic crystals are presented, including 2D cases with square lattice arrange-ment, and 3D cases with SC lattice, BCC lattice, and FCC lattice arrangements. In all of the numerical calculations of phononic crystals, the background and filling materi-als chosen are epoxy and steel. These materimateri-als are well-known for the study of the property of both acoustic band gaps and total band gaps [21], [28]. Thus it is convenient to choose them for demonstrating the results of 2D and 3D phononic crystals in this study. The density and the elastic constants C11 and C44 of epoxy are assumed as
1180 kg/m3, 7.61 GPa, and 1.59 GPa, respectively, and those for steel are 7780 kg/m3, 264 GPa, and 81 GPa. To proceed with the calculation of the FDTD method, the phononic crystals are created by assigning the mate-rial constants of corresponding positions according to the shape of inclusions. In this study, although the phononic crystals are formed with isotropic background and filling materials, the anisotropic lattice structures demonstrate the anisotropic property of wave propagation, as shown below.
A. Two-Dimensional Phononic Crystals
In this subsection, we consider a 2D steel/epoxy phononic crystal with square lattice arrangement as shown in Fig. 4(a). The unit cell of the square lattice phononic crystal can be chosen as either the square with a full in-clusion inside or a square with four quarters of the inclu-sion. In this study, we adopt the former one for simplicity. The first Brillouin zone of the phononic crystal in square arrangement is shown in Fig. 4(b). Owing to the symme-try of the material and lattice, the irreducible part of the first Brillouin zone is an isosceles right triangle. Therefore we calculate the dispersion relations for a BAW along the boundary of the irreducible part of the first Brillouin zone. The symbols Γ, X, and M in Fig. 4(b) denote the vertexes
Fig. 5. Dispersion relation of a circular inclusion phononic crystal with square lattice arrangement.
of the first Brillouin zone of square-arrayed phononic crys-tals.
First, we discuss the band gap phenomena of a phononic crystal with circular inclusion. The radii of the cylinder are 9 and 15 grids while the lattice constant is 30 and 50 grids. The filling fraction of the 2D phononic crystal is defined as the ratio of the cross-section areas of the inclusion and the unit cell. Therefore, the filling fraction f is 0.283 for these cases. To verify the result of the FDTD method, the dispersion relation is also calculated by the PWE method. There are 1681 reciprocal lattice vectors (RLV) for the PWE method while each unit cell is divided into 30 by 30 and 50 by 50 grids for the FDTD method. The dispersion curves are shown in Fig. 5, where the lines with solid trian-gle symbols “” represent the results of the PWE method and the lines with hollow circle “◦” and hollow rhombus “” symbols represent those of the FDTD method. The dispersion relation shows the allowed elastic wave modes inside the phononic crystal. The vertical axis is the nor-malized frequency Ω = ωa/Ct and the horizontal axis is
the reduced wave vector defined as ka/π, where ω is the angular frequency, k the wave vector, a lattice constant, and Ctthe transverse wave velocity of the background
ma-terial.
As Fig. 5 shows, the results show good agreement between the PWE and the FDTD methods. The PWE method shows the total band gap from the normalized frequency 4.07 to 6.79, and the FDTD method shows the band gap from Ω = 4.08 to 6.67 and 4.05 to 6.67 for 30×30 and 50×50 grids per unit cell, respectively. The dispersion curves match each other in the lower frequency range and show slight differences in higher frequency range. The pos-sible reason is that the presentation of the inclusion shape is limited by the orthogonal grids system. Therefore, the
resonance frequencies of higher eigenmodes change slightly when the phononic crystal unit cell is defined by more grids. Basically, the result of the FDTD method with the unit cell divided into 30× 30 grids shows very high agree-ment with that for 50× 50 grids, but the calculation time is less than half that of the latter. It is also worth indicat-ing that the numbers of reciprocal lattice vectors (RLV) affect the accuracy of dispersion curves. The PWE work shows that with the larger number of RLV, the dispersion curves shift to a lower value. With the results, the FDTD method is verified and the division of the unit cell is also discussed for the trade-off of accuracy and efficiency.
Second, Fig. 6 shows the results of 2D phononic crystals with a square lattice arrangement but the inclusion is a square rod. Although the orthogonal discrete grids may limit accuracy of the inclusion shapes, the FDTD method is very convenient to study the phononic crystals with a different inclusion shape. The square inclusions are defined by designating the material constants in the corresponding grids, and furthermore the square inclusions are orientated in different directions, such as 0, 15, 30, and 45 degrees with respect to the X-axis. These cases have an equal filling fraction, but the apparent band gaps change due to the orientation of square inclusions. Fig. 6(a) shows the results of 2D phononic crystals whose filling fraction is 0.27. The filling fraction is approximately equal to, but smaller than, that in Fig. 5. The hollow rhombus “” symbols of the 0 degree orientation case, for example, represent a total band gap width of 2.58, with a normalized frequency from 3.97 to 6.55. Meanwhile, Fig. 6(b) shows the dispersion curves of the phononic crystals with filling fraction 0.36, and thus the length of the side of square rod is 30 grids, equal to the diameter of the circle in Fig. 5. A total band gap of the 0 degree orientation case appears to have a normalized frequency from 3.73 to 7.67 and the width of the band gap is 3.94. The other symbols, hollow circle, triangle, and square, represent the cases with orientation 15, 30, and 45 degrees, respectively, and the band gaps are included in Fig. 6(c).
When the filling fractions are equal, we find the total band gap widths of the phononic crystal with the square inclusion in the square lattice arrangement approximately equal to that with the circular inclusion, but the interval changes to a higher frequency range with the increase of the orientation angle. In addition, if a general radius is defined as the vertical distance from the center of the ge-ometric shape to the boundary such as the half length of a square side and the radius of a circle, then the cases of equal general radius show that the crystal with the square rod inclusion exhibits a wider total band gap.
B. Three-Dimensional Phononic Crystals
By replacing the material constant of each grid, the FDTD method is appropriate to study the 3D phononic crystals. We analyze three types of 3D phononic crystals in this subsection, including SC lattice, BCC lattice, and FCC lattice arrangements. Due to the considerable
com-(a)
(b)
(c)
Fig. 6. Dispersion relation of a square inclusion phononic crystal with different orientations: (a) The case of filling fraction 0.27; (b) The case of filling fraction 0.36; (c) The distribution of total band gaps.
(a)
(b)
Fig. 7. A 3D simple cubic lattice phononic crystal with a sphere in-clusion: (a) Lattice arrangement and the first Brillouin zone; (b) Dis-persion relation.
putation of 3D cases, we introduce the PC cluster system for the numerical work. In addition, the unit cell of the 3D phononic crystal is divided into 30× 30 × 30 grids. It is believed that the accuracy is good enough from the study of the 2D cases in the previous subsection. The time step is chosen as 0.002a/Ctto satisfy the convergence condition
of the FDTD method in wave propagation, and the calcu-lation is with the time evolution of 100,000 time steps.
Fig. 7(a) shows the unit cell of the 3D phononic crystal with a sphere inclusion in the SC lattice arrangement. The first Brillouin zone of an SC lattice is also shown in the figure. Due to the symmetric property of geometry and material, the irreducible part is a tetrahedron with corners labeled Γ, X, M, and R. The radius of the sphere is 9 grids, and thus the filling fraction of this SC phononic crystal is 0.113. The dispersion relation for a BAW in the first Brillouin zone is represented in Fig. 7(b) and a total band gap is observed from a normalized frequency of 4.68 to 4.90 with a width of 0.22.
(a)
(b)
Fig. 8. A 3D simple cubic lattice phononic crystal with a cube inclu-sion: (a) Lattice arrangement; (b) Dispersion relation.
Next, we replace the sphere inclusion in the SC lattice with a cube inclusion, as shown in Fig. 8(a). The length of the cube side equals 16 grids and the filling fraction is 0.152. The dispersion relation is analyzed along the irre-ducible volume again and a band gap from 4.48 to 5.38 is obtained. The width of the band gap is 0.9, which is wider than that of the sphere inclusion. The band gap be-comes larger when the side length of cube inclusion equal 18 grids. Therefore a cube inclusion can create a larger band gap than a sphere one for an equal general radius.
Further, the 3D phononic crystal with a BCC lattice arrangement and its first Brillouin zone are shown in Fig. 9(a). The tips of the irreducible part are labeled Γ, H, N, and P. The radius of the sphere is 9 grids, the same as in the SC case. Since each cubic unit cell in the BCC lat-tice arrangement contains two spheres, the filling fraction is 0.226. The dispersion curves for a BAW are calculated along the selected direction and shown in Fig. 9(b). A to-tal band gap appears from a normalized frequency of 5.37 to 8.07 and the width is 2.7.
(a)
(b)
Fig. 9. A 3D body-centered cubic lattice phononic crystal with sphere inclusions: (a) Lattice arrangement and the first Brillouin zone; (b) Dispersion relation.
The final 3D phononic crystal case is the FCC lattice ar-rangement in Fig. 10(a). The tips of the irreducible part of first Brillouin zone are labeled Γ, X, W, K, L, and U. With the same dimension of the lattice constant and sphere, each cubic unit cell contains four spheres and the filling fraction is 0.452. There is a noticeable total band gap from 6.79 to 14.60 with a width of 7.81.
Thus we find that the 3D phononic crystal with the cube inclusion has a larger band gap than that with the spherical inclusion; this is also true in the 2D case. It is apparent that in the 3D phononic crystal the band gap phenomenon is affected by the geometric shape of the in-clusion. Furthermore, the 3D phononic crystals, including SC, BCC, and FCC lattice arrangements, show that the band gap width in the FCC case is larger than those of the BCC and SC cases and that the band gap width of the BCC case is larger than that of the SC case. We can also inspect the phononic band gaps in terms of relative band width ∆ω/ω0, defined as ratio of the band gap to the
central frequency of gap. The relative band width values of the SC, BCC, and FCC lattice cases of spherical inclu-sions are 0.046, 0.40, and 0.73, respectively. So the band gap enlarges as the filling fraction increases, and a high-density package arrangement has a larger band gap in the
(a)
(b)
Fig. 10. A 3D face-centered cubic lattice phononic crystal with sphere inclusions: (a) Lattice arrangement and the first Brillouin zone; (b) Dispersion relation.
3D phononic crystal. The FCC lattice has a considerable total band gap of its BAW and thus can potentially be applied to the design of 3D waveguides and filters.
V. Conclusions
This paper proposes a FDTD theory of the BAW in 3D phononic crystals consisting of materials with orthotropic symmetry. The FDTD method is very suitable for deal-ing with various periodic phononic structures and easily satisfies the convergence condition of materials with large elastic mismatch. With the periodic boundary condition derived from Bloch’s theorem, the FDTD method can be used to analyze the dispersion relations of 2D and 3D phononic crystals efficiently. We also develop and execute parallel FDTD programs on a PC cluster system to mini-mize the calculating time for the 3D structures. The paral-lel programs are realized by distributing jobs with specific wave vectors, and they have increased our work efficiency many times.
In the discussion of phononic crystals, we generalize the relation between the band gap phenomenon and the geo-metric shape by extensive study of phononic crystals of different geometric designs. Some band gap characteris-tics of both 2D and 3D phononic crystals drawn from the aforementioned study are as follows:
1. In the 2D cases of square lattice arrangement, phononic crystals with the square inclusion arrange-ment have a larger band gap. The range of total band gaps varies with the orientation angle of the square inclusion while the widths do not show obvious differ-ences.
2. The geometric shapes of inclusions in 3D phononic crystals also change the appearance of total band gaps. A cube inclusion creates a larger band gap than a sphere inclusion of equal general radius.
3. The band gaps of 3D phononic crystals with a constant sphere size are calculated for SC, BCC, and FCC lat-tice arrangements. The largest band gap is obtained in the FCC case and the next largest for the BCC case. This can be understood by the increase of the filling fraction. These analyses are the fundamental steps toward development of further applications. For example, the BCC lattice phononic crystal with real dimensions of the 10-mm lattice constant and the 3-mm spherical radius supports a total band gap from 99.2 kHz to 149.1 kHz. This property can be of use for an acoustic wave filter and further to accordingly design a 3D acoustic waveguide.
References
[1] M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. Djafari-Rouhani, “Acoustic band structure of periodic elastic compos-ites,” Phys. Rev. Lett., vol. 71, no. 13, pp. 2022–2025, 1993. [2] M. S. Kushwaha, P. Halevi, G. Martinez, L. Dobrzynski, and B.
Djafari-Rouhani, “Theory of acoustic band structure of periodic elastic composites,” Phys. Rev. B, vol. 49, no. 4, pp. 2313–2322, 1994.
[3] J. O. Vasseur, B. Djafari-Rouhani, L. Dobrzynski, M. S. Kush-waha, and P. Halevi, “Complete acoustic band gaps in periodic fibre reinforced composite materials: The carbon/epoxy compos-ite and some metallic systems,” J. Phys.: Condens. Matter, vol. 6, pp. 8759–8770, 1994.
[4] C. Goffaux and J. P. Vigneron, “Theoretical study of a tunable phononic band gap system,” Phys. Rev. B, vol. 64, no. 075118, 2001.
[5] F. Wu, Z. Liu, and Y. Liu, “Acoustic band gaps in 2D liquid phononic crystals of rectangular structure,” J. Phys. D, vol. 35, pp. 162–165, 2002.
[6] F. Wu, Z. Liu, and Y. Liu, “Acoustic band gaps created by rotating square rods in a two-dimensional lattice,” Phys. Rev.
E, vol. 66, no. 046628, 2002.
[7] X. Li, F. Wu, H. Hu, S. Zhong, and Y. Liu, “Large acoustic band gaps created by rotating square rods in two-dimensional periodic composites,” J. Phys. D: Appl. Phys., vol. 36, pp. L15– L17, 2003.
[8] T.-T. Wu and Z.-G. Huang, “Level repulsions of bulk acoustic waves in composite materials,” Phys. Rev. B, vol. 70, no. 214304, 2004.
[9] Y. Tanaka and S. Tamura, “Surface acoustic waves in two-dimensional periodic elastic structures,” Phys. Rev. B, vol. 58, no. 12, pp. 7958–7965, 1998.
[10] Y. Tanaka and S. Tamura, “Acoustic stop bands of surface and bulk modes in two-dimensional phononic lattices consisting of aluminum and a polymer,” Phys. Rev. B, vol. 60, no. 19, pp. 13294–13297, 1999.
[11] T.-T. Wu, Z.-G. Huang, and S. Lin, “Surface and bulk acoustic waves in two-dimensional phononic crystals consisting of mate-rials with general anisotropy,” Phys. Rev. B, vol. 69, no. 094301, 2004.
[12] Z.-G. Huang and T.-T. Wu, “Temperature effect on the bandgaps of surface and bulk acoustic waves in two-dimensional
phononic crystals,” IEEE Trans. Ultrason., Ferroelect., Freq.
Contr., vol. 52, no. 3, pp. 365–370, 2005.
[13] T.-T. Wu, Z.-C. Hsu, and Z.-G. Huang, “Band gaps and the electromechanical coupling coefficient of a surface acoustic wave in a two-dimensional piezoelectric phononic crystal,” Phys. Rev.
B, vol. 71, no. 064303, 2005.
[14] M. Kafesaki and E. N. Economou, “Multiple-scattering theory for three-dimensional periodic acoustic composites,” Phys. Rev.
B, vol. 60, no. 17, pp. 11993–12001, 1999.
[15] I. E. Psarobas, N. Stefanou, and A. Modinos, “Scattering of elastic waves by periodic arrays of spherical bodies,” Phys. Rev.
B, vol. 62, no. 1, pp. 278–291, 2000.
[16] Z. Liu, C. T. Chan, P. Sheng, A. L. Goertzen, and J. H. Page, “Elastic wave scattering by periodic structures of spherical ob-jects,” Phys. Rev. B, vol. 62, no. 4, pp. 2446–2457, 2000. [17] J. Mei, Z. Liu, J. Shi, and D. Tian, “Theory for elastic wave
scattering by a two-dimensional periodical array of cylinders: An ideal approach for band-structure calculations,” Phys. Rev.
B, vol. 67, no. 245107, 2003.
[18] D. Garica-Pablos, M. Sigalas, F. R. Montero de Espinosa, M. Torres, M. Kafesaki, and N. Garcia, “Theory and experiments on elastic band gaps,” Phys. Rev. Lett., vol. 84, no. 19, pp. 4349– 4352, 2000.
[19] Y. Tanaka, Y. Tomoyasu, and S. Tamura, “Band structure of acoustic waves in phononic lattices: Two-dimensional composites with large acoustic mismatch,” Phys. Rev. B, vol. 62, no. 11, pp. 7387–7392, 2000.
[20] A. Khelif, B. Djafari-Rouhani, V. Laude, and M. Solal, “Cou-pling characteristics of localized phonons in photonic crystal fibers,” J. Appl. Phys., vol. 94, no. 12, pp. 7944–7946, 2003. [21] J.-H. Sun and T.-T. Wu, “Analyses of mode coupling in joined
parallel phononic crystal waveguides,” Phys. Rev. B, vol. 71, no. 174303, 2005.
[22] M. S. Kushwaha and B. Djafari-Rouhani, “Complete acoustic stop bands for cubic arrays of spherical liquid balloons,” J. Appl.
Phys., vol. 80, no. 6, pp. 3191–3195, 1996.
[23] I. E. Psarobas and M. M. Sigalas, “Elastic band gaps in a fcc lattice of mercury spheres in aluminum,” Phys. Rev. B, vol. 66, no. 052302, 2002.
[24] R. Sainidou, N. Stefanou, and A. Modinos, “Formation of ab-solute frequency gaps in three-dimensional solid phononic crys-tals,” Phys. Rev. B, vol. 66, no. 212301, 2002.
[25] X. Zhang, Z. Liu, Y. Liu, and F. Wu, “Elastic wave band gaps for three-dimensional phononic crystals with two structural units,” Phys. Lett. A, vol. 313, pp. 455–460, 2003.
[26] J. D. Achenbach, Wave Propagation in Elastic Solids. New York: North-Holland Publishing Company, 1976.
[27] W. Gropp, E. Lusk, N. Doss, and A. Skjellum, “A high-performance, portable implementation of the MPI message pass-ing interface standard,” Parallel Comput., vol. 22, pp. 789–828, 1996.
[28] J. O. Vasseur, P. A. Deymier, B. Chenni, B. Djafari-Rouhani, L. Dobrzynski, and D. Prevost, “Experimental and theoretical evidence for the existence of absolute acoustic band gaps in two-dimensional solid phononic crystals,” Phys. Rev. Lett., vol. 86, no. 14, pp. 3012–3015, 2001.
Po-Feng Hsieh received his B.S. degree from
the Department of Mechanical Engineering, National Chiao Tung University, and his M.S. degree from the Institute of Applied Mechan-ics, National Taiwan University, in 2002 and 2004, respectively. Currently, he is a research engineer in the Behavior Design Corpora-tion, Science-Based Industrial Park, Hsinchu, Taiwan, R.O.C. His research work is mainly on nature language processing and machine translation.
Tsung-Tsong Wu received his doctorate in
theoretical and applied mechanics from Cor-nell University, Ithaca, NY, in 1987. Then he joined the National Taiwan University faculty and is now a professor in the Institute of Ap-plied Mechanics. He was awarded the distin-guished research prizes of the National Sci-ence Council (NSC) three times for six years from 1995 to 2001 and is currently a distin-guished research fellow of the NSC. He is the executive board director of the Taiwanese So-ciety of Nondestructive Testing and the board director of the Quartz Industry Association of Taiwan. He is a mem-ber of the American Society of Mechanical Engineers.
Jia-Hong Sun received his B.S. degree from
the Department of Aeronautics and Astronau-tics, National Cheng Kung University, in 1995 and his M.S. degree from the Institute of Ap-plied Mechanics, National Taiwan University, in 1997. Currently, he is a Ph.D. candidate in the Institute of Applied Mechanics, Na-tional Taiwan University. His research work is mainly on phononic crystals and the calcu-lation of elastic wave propagation using the FDTD method.
