- Email: [email protected]
Presentation and investigation of a new two dimensional heterostructure phononic crystal to obtain extended band gap
Physica B 489 (2016) 28–32
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Presentation and investigation of a new two dimensional heterostructure phononic crystal to obtain extended band gap Mohammad Bagheri Nouri n, Mehran Moradi Department of Mechanical Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran
art ic l e i nf o
a b s t r a c t
Article history: Received 1 January 2016 Received in revised form 15 February 2016 Accepted 17 February 2016 Available online 18 February 2016
In this paper, a new heterostructure phononic crystal is introduced. The new heterostructure is composed of square and rhombus phononic crystals. Using finite difference method, a displacement-based algorithm is presented to study elastic wave propagation in the phononic crystal. In contrast with conventional finite difference time domain method, at first by using constitutive equations and straindisplacement relations, elastic wave equations are derived based on displacement. Then, these forms are discretized using finite difference method. By this technique, components of stress tensor can be removed from the updating equations. Since the proposed method needs less elementary arithmetical operations, its computational cost is less than that of the conventional FDTD method. Using the presented displacement-based finite difference time domain algorithm, square phononic crystal, rhombus phononic crystal and the new heterostructure phononic crystal were analyzed. Comparison of transmission spectra of the new heterostructure phononic crystal with those creating lattices, showed that band gap can be extended by using the new structure. Also it was observed that by changing the angular constant of rhombus lattice, a new extended band gap can be achieved. & 2016 Elsevier B.V. All rights reserved.
Keywords: Heterostructure phononic crystal Band gap extension Finite difference time domain Displacement-based formulation
1. Introduction Much attention has been attracted to the phononic crystals in recent years. These non homogenous structures are periodic arrangement of inclusions in an elastically different host material. The frequency ranges that propagation of elastic waves is prohibited are called band gap. By this feature, these crystals can control propagation of elastic waves [1,2]. Since the capability of the crystals arise from phononic band gaps, it is essential to achieve the desired phononic band gaps. Acoustic filter and waveguide are some applications of phononic crystals. Phononic crystals are studied theoretically and experimentally in several researches [3–6]. To achieve the desired band gap, various works have been done on these crystals such as introduction of nesting phononic crystals [7], construction of a hybrid triangular lattice phononic crystal [8], presentation of twinned-square periodic structures [9] and introduction of hierarchical phononic crystals [10]. Heterostructure phononic crystals have been also investigated in some studies recently such as: presenting of a sandwich type of heterostructure in the 2D square lattice by replacing two columns of circular scatter cylinders with two columns of square scatter n
Corresponding author. E-mail address: [email protected] (M. Bagheri Nouri).
http://dx.doi.org/10.1016/j.physb.2016.02.023 0921-4526/& 2016 Elsevier B.V. All rights reserved.
cylinders [11], investigation of the interface-guided mode of Lamb waves in a phononic crystal heterostructure plate [12], analyzing of heterostructure phononic crystal waveguides [13], construction of an one-dimensional multilayer phononic heterostructure by combining finite periodic phononic crystals and Fibonacci quasiperiodic phononic crystals [14] and design of a phononic heterostructure by merging two one-dimensional phononic crystals [15]. By combining two-dimensional square and rhombus lattices, a new heterostructure phononic crystal is presented in this paper. Then, square and rhombus phononic crystals and the heterostructure phononic crystal were analyzed using displacementbased finite difference time domain method. Band gap of the new heterostructure phononic crystal was computed and compared with band gaps of its square and rhombus lattices. Effect of angular constant of the rhombus lattice on band gap of the heterostructure phononic crystal was studied too.
2. Model and calculation method 2.1. The model The proposed heterostructure phononic crystal has been shown in Fig. 1. The proposed heterostructure is composed of square and rhombus phononic crystals. This heterostructure is finite in horizontal direction and extends infinitely in vertical
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
29
u2l + 1/2, m + 1/2; n + 1 = u2l + 1/2, m + 1/2; n + v2l + 1/2, m + 1/2; n + 1/2 Δt
(8)
l + 1/2, m; n l + 1/2, m σ11 = C11 (u1l + 1, m; n − u1l, m; n )/Δx l + 1/2, m + C12 (u2l + 1/2, m + 1/2; n − u2l + 1/2, m − 1/2; n )/Δy
l, m + 1/2; n l, m + 1/2 ⎡ l, m + 1; n σ12 = C44 − u1l, m; n ⎣ u1
(
+ Fig. 1. The proposed hetero-structure phononic crystal that is composed of square arrangement (left) and rhombus arrangement (right).
a′ = a/(2 sin (γ /2))
(1)
where a′ and γ are constants of the rhombus lattice and a is constant of square lattice. 2.2. Calculation method Two-dimensional (2D) phononic crystals are considered here. The formulation can be easily derived for 3D phononic crystals. The inclusions are arranged along the z direction and are repeated in the xy plane. The equations of elastic waves can be written as:
ρvi̇ = σij, j
(2)
vi = u̇i
(3)
σij = Cijmn um, n
(4)
In above equations, ρ = ρ (x, y ), Cijkl (x, y ), vi and ui stand for the density, elastic stiffness tensor, i component of velocity and displacement of the structure, respectively. The summation agreement over dummy indices is considered, too. Since propagation of elastic waves in the xy plane is assumed, the displacement, velocity and stress tensor of the structure are independent of z, i.e., ui = ui (x, y, t ), vi = vi (x, y, t ) and σij = σij (x, y, t ). Consider that the inclusion and the host materials are isotropic. Then the discretized form of Eqs. (2–4) for study of the mixed mode is given as below: l + 1/2, m; n l − 1/2, m; n v1l, m; n + 1/2 = v1l, m; n − 1/2 + Δt /ρl, m ⎡⎣ σ11 − σ11
(
+
(
l, m + 1/2; n σ12
v2l + 1/2, m + 1/2; n + 1/2
=
−
l, m − 1/2; n σ12
v2l + 1/2, m + 1/2; n − 1/2
)
Δy⎤⎦
+
Δt /ρl + 1/2, m + 1/2
(5)
l + 1/2, m + 1; n l + 1/2, m; n × ⎡⎣ σ22 − σ22
(
(
l + 1, m + 1/2; n l, m + 1/2; n + σ12 − σ12
u1l, m; n + 1 = u1l, m; n + v1l, m; n + 1/2 Δt
) Δx
) Δy
) Δx⎤⎦
(6)
(7)
−
) Δy
u2l − 1/2, m + 1/2; n
(
) Δx⎤⎦
l + 1/2, m; n l + 1/2, m σ22 = C11 u2l + 1/2, m + 1/2; n − u2l + 1/2, m − 1/2; n
+ direction. The relation between constants of square and rhombus lattices can be written as:
(
u2l + 1/2, m + 1/2; n
l + 1/2, m C12
(
u1l + 1, m; n
−
(9)
u1l, m; n
) Δx
(10)
) Δy (11)
where (l,m) and n stands for location of a typical node and n-th time step, respectively. As Eq. (5) reveals, updating of x component of velocity at node (l,m) requires calculation of components of stress tensor at four l + 1/2, m; n l − 1/2, m; n l, m + 1/2; n l, m − 1/2; n coordinates, i.e., σ11 , σ11 , σ12 , σ12 . As Eqs. (9– 11) implies, computation of each of these stress components needs sizeable elementary arithmetical operations. To calculate transmission spectra of a phononic crystal, Fourier transformation of displacement data are taken. To compute transmission spectra of a phononic crystal, the calculated stress components have no direct contributions and are only used for updating components of displacement. By obtaining displacement-based forms of elastic wave equations and discretization of resultant equations, components of stress tensor can be eliminated from the updating equations. This course of action gives efficient updating equations that need less elementary arithmetical operations than conventional FDTD method. So the needed computational cost of the efficient updating equations is less than that of conventional FDTD method. The procedure of obtaining these updating equations is described below. For points far apart the interface of the inclusion and the host material (elastic constants of the structure are not discontinuous in these points), by substituting divergence of Eq. (4) in Eq. (2) the displacement-based form of elastic wave equation can be derived. By discretization of the resultant equations, displacement-based finite difference time domain (DBFDTD) formulation can be obtained. Consider that the inclusion and the host materials are isotropic. The displacement-based form of elastic wave equations for study of the mixed mode in this case can be given as below:
ρv1̇ = C11u1,11 + C44 u1,22 + ( C12 + C44 ) u2,12
(12)
ρv2̇ = C11u2,22 + ( C12 + C44 ) u1,12 + C44 u2,11
(13)
In Eqs. (12) and (13) C11, C12 and C44 are elastic constants of an isotropic material. To discretize Eqs. (12) and (13), all derivatives are replaced with central difference estimates. The consequent equations can be applied to update the velocity and displacement of the non interfacial grid points. For the interfacial grid points, the discretized form of the Eqs. (2–4) are used. For a suitable presentation of the discretized form of Eqs. (12) and (13), the subsequent coefficients are defined:
30
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
l, m λ1l, m = C11 /Δx2 l, m λ2l, m = C44 /Δy2
(
)
l, m l, m λ3l, m = C12 + C44 /(ΔxΔy)
Fig. 2. The computational model used for finite square phononic crystal analysis.
l + 1/2, m + 1/2 η1l + 1/2, m + 1/2 = C11 /Δy2 l + 1/2, m + 1/2 η2l + 1/2, m + 1/2 = C44 /Δx2
(
l + 1/2, m + 1/2 l + 1/2, m + 1/2 η3l + 1/2, m + 1/2 = C12 + C44
) (ΔxΔy)
(14)
The coefficients introduced by Eq. (14) may be calculated before the time evolution loop of the DBFDTD method and saved in appropriate matrixes in order to avoid computing them in each time step. Discretizing Eq. (12) using the finite difference technique leads to:
v1l, m; n + 1/2 = v1l, m; n − 1/2 + Δt/ρl, m Σ
(15)
where:
(
)
(
Σ = − 2 λ1l, m + λ2l, m u1l, m; n + λ1l, m u1l + 1, m; n + u1l − 1, m; n +
λ2l, m
+
λ3l, m
( (u
u1l, m + 1; n
+
u1l, m − 1; n
l + 1/2, m + 1/2; n 2
−
Fig. 3. Transmission spectra through the square phononic crystal obtained by DBFDTD calculation along the ГX direction: for longitudinal (solid line) and transverse polarization (dashed line) of incident waves.
)
)
u2l + 1/2, m − 1/2; n
− u2l − 1/2, m + 1/2; n + u2l − 1/2, m − 1/2; n
)
(16)
Discretized form of Eq. (13) can be given as:
v2l + 1/2, m + 1/2; n + 1/2
= v2l + 1/2, m + 1/2; n − 1/2 + Δt/ρl + 1/2, m + 1/2 χ
(17)
where:
(
)
χ = − 2 η1l + 1/2, m + 1/2 + η2l + 1/2, m + 1/2 u2l + 1/2, m + 1/2; n
( (u (u
+ η1l + 1/2, m + 1/2 u2l + 1/2, m + 3/2; n + u2l + 1/2, m − 1/2; n +
η2l + 1/2, m + 1/2
+
η3l + 1/2, m + 1/2
l + 3/2, m + 1/2; n 2
l + 1, m + 1; n 1
+
u2l − 1/2, m + 1/2; n
) )
− u1l + 1, m; n − u1l, m + 1; n + u1l, m; n
)
(18)
Subsequently, Eqs. (7) and (8) can be used to update components of displacement of the structure. The efficiency of DBFDTD formulation can be easily understood by comparison of Eqs. (15) and (16) with Eqs. (5), (9) and (10). This comparison shows that Eqs. (15) and (16) have about 50 percent less elementary arithmetical operations.
3. Results and discussion In following examples, the inclusion and the matrix are steel and epoxy, respectively. The density and elastic constants C11 and C44 of steel are assumed to be 7780 kg/m3, 264 GPa and 81 GPa respectively, and those for epoxy are 1142 kg/m3, 7.54 GPa and 1.48 GPa. To validate result of the DBFDTD method, transmission spectra of a square phononic crystal was computed and compared with result of the FDTD method. Diameter of steel cylinders and the lattice constant are 5 and 7 mm, respectively. To calculate transmission spectra of the crystal, a sample with five unit cells in x direction and one unit cell in y direction was oriented along ГX direction and surrounded by two homogenous regions (Fig. 2). The perfectly matched layer (PML) [16] was adopted as the absorbing boundary condition for the left and right boundaries of the computational domain. For the top and bottom boundaries the periodic boundary condition was used.
A Gaussian wave packet was launched along the x direction at the left homogenous region. The x component (y component) of displacement was recorded at middle of right homogenous region and averaged over a lattice constant. By taking fast Fourier transform of this average, transmission spectra through the crystal for longitudinal (transverse) polarization of incident waves was obtained. Fig. 3 exhibits DBFDTD results of the transmission spectra through the crystal. To obtain the transmission spectra, a grid of 60*60 points in the unit cell was considered and the equations of motion were solved over 216 time steps with each time step lasting 7.08 ns. Fig. 3 shows a normalized frequency range extending from 3.8 to 8.7 that the crystal forbids propagation of elastic waves along ГX direction for any polarization of incident waves. Normalized frequency is defined as ωa/vt where ω , a and vt are respectively angular frequency, the lattice constant and transverse velocity of elastic waves in epoxy. This means that there is a band gap along ГX direction extending from 3.8 to 8.7. The computed band gap is in good agreement with that calculated by the FDTD method (3.8–8.5) [17]. 3.1. The heterostructure composed of square and triangular lattices First, the heterostructure composed of square and triangular phononic crystals was studied. Diameter of steel rods is 4 mm. Lattice constant of both square and triangular Lattices are 6 mm. Fig. 4 shows the sample was used to calculate transmission spectra of the heterostructure. The sample was surrounded by two homogenous regions (not shown here). The PML was used as the absorbing boundary condition for the boundaries along the x direction. For the boundaries along the y direction the periodic boundary condition was applied. A Gaussian wave packet was launched along the x direction at the left homogenous region. The x component of displacement was recorded at middle of right homogenous region and averaged over a lattice constant. By taking fast Fourier transform of this
Fig. 4. The computational model used for heterostructure phononic crystal analysis.
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
31
Fig. 5. Transmission spectra through the square lattice (dot line), the triangular lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction computed for longitudinal polarization of incident wave. The longitudinal component Ux of the displacement is given in arbitrary units.
Fig. 7. Transmission spectra of square lattice (dot line), rhombus lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction of propagation computed for longitudinal incident wave. The longitudinal component Ux of the output displacement field is given in arbitrary units.
average, transmission spectra through the heterostructure was obtained. Fig. 5 exhibits DBFDTD results of the transmission spectra through the heterostructure for longitudinal polarization of incident waves. To obtain the transmission spectra, a grid of 60*60 points in the unit cell of square lattice was considered and the equations of motion were solved over 216 time steps with each time step lasting 6.09 ns. The CPU times to calculate transmission spectra are 2 min and 16 s for the DBFDTD method and 3 min and 34 s for the conventional FDTD method on twenty 3 GHz Intel Xeon CPUs. So the computations time of the DBFDTD method is 36 percents less than that of the conventional FDTD method. The numerical examples were calculated by a developed FORTRAN code based on parallel processing. Fig. 5 shows that there is a band gap extending from 110 to 238 kHz that the square phononic crystal forbids propagation of elastic waves along the x direction. The band gaps of the triangular phononic crystal locate at 120–288 and 313–367 kHz. In the frequency range of Fig. 5, two band gaps were found for the heterostructure, i.e., 113–288 and 313-365 kHz. As Fig. 5 indicates the heterostructure only allows propagation of the waves that both the square and triangular lattices allow this propagation. So the band gap can be extended by combination of square and triangular lattices and construction of the heterostructure. Due to proximity of two constituent lattices’ band gap ranges, heterostructure transmission spectra are approximately same as triangular lattice spectra.
and 400 kHz for the rhombus phononic crystal and the band gap of the heterostructure locates at 110–400 kHz. Fig. 7 implies that by combination of the square and the rhombus lattices and creation of the heterostructure, a new extended band gap can be achieved. Comparison of these transmission spectra shows that heterostructure band gap has been extended extensively. As the last example, the heterostructure composed of square and rhombus lattices was considered. The angular constant of the rhombus phononic crystal is 30° here. The diameter of steel rods and the lattice constant of the square phononic crystal are 4 and 6 mm, respectively. According to Eq. (1), the lattice constant of rhombus phononic crystal is 11.59 mm. Fig. 8 shows the sample was used to calculate transmission spectra of the heterostructure. The transmission spectra through the heterostructure along x direction for longitudinal polarization of incident waves were shown in Fig. 9. To calculate the spectra, the same procedure of the previous examples was complied. In the frequency range of Fig. 9, five band gaps appear for the rhombus phononic crystal, i.e., 67–148, 173–246, 274–282, 341– 362 and 382–387 kHz. Also there are five band gaps for the heterostructure in the frequency domain of Fig. 9, i.e., 68–246, 274– 286, 339–362 and 382–387 kHz. As Fig. 9 indicates, a new different extended band gap can be attained using the heterostructure composed of the square and the rhombus lattices. According to the obtained results, the band gap of the heterostructure phononic crystal would be the wider than those of other studied heterostructures if the angular constant of the rhombus lattice is 80°.
3.2. The heterostructure composed of square and rhombus lattices 4. Conclusion As second example, the heterostructure composed of square and rhombus lattices was considered. The angular constant of the rhombus phononic crystal is 80°. The diameter of steel rods and the lattice constant of the square phononic crystal are 4 and 6 mm, respectively. According to Eq. (1), the lattice constant of rhombus phononic crystal is 4.67 mm. Fig. 6 shows the sample was used to calculate transmission spectra of the heterostructure. The transmission spectra through the heterostructure along x direction for longitudinal polarization of incident waves were shown in Fig. 7. To obtain the spectra, the same procedure of the previous example was obeyed. In the frequency range of Fig. 7, a band gap occurs between 183
Fig. 6. same as Fig. 4 but here, the triangular lattice was replaced with a rhombus lattice whose angular constant is 80°.
In this paper a new heterostructure phononic crystal was presented. The heterostructure was composed of the square and the rhombus phononic crystals. Also displacement-based finite difference time domain method was introduced and applied for simulation of elastic waves propagation through phononic crystals. Comparison of the transmission spectra computed by the DBFDTD method with the results of FDTD method, proved the efficiency of the DBFDTD method. Also Comparison of the computational time showed that the computational time of the DBFDTD method is 36 percents less than that of the FDTD method. Propagation of elastic waves through the proposed heterostructure phononic crystal was studied using the DBFDTD method. Investigation of three calculated examples demonstrated that the
Fig. 8. same as Fig. 6 but the angular constant of rhombus lattice is 30° here.
32
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
Fig. 9. Transmission spectra of square lattice (dot line), rhombus lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction of propagation computed for longitudinal incident wave. The longitudinal component Ux of the output displacement field is given in arbitrary units.
band gap can be extended using the heterostructure phononic crystal. Effect of the angular constant of the rhombus phononic crystal on the band gap of the heterostructure was studied too. This investigation showed that the band gap of the heterostructure phononic crystal would be the wider than those of other studied heterostructures if the angular constant of the rhombus lattice is 80°.
References [1] M.S. Kushwaha, P. Halevi, L. Dobrzynski, B. Djafari-Rouhani, Acoustic band structure of periodic elastic composites, Phys. Rev. Lett. 71 (1993) 2022–2025. [2] R. Martinez-Sala, J. Sancho, J.V. Sanchez, V. Gomez, J. Llinares, et al., Sound attenuation by sculpture, Nature 378 (1995) 241-241.
[3] F.R. Montero de Espinosa, E. Jime´nez, M. Torres, Ultrasonic band gap in a periodic two-dimensional composite, Phys. Rev. Lett. 80 (1998) 1208–1211. [4] A. Khelif, A. Choujaa, R. laihem, M. Wilm, S. Ballandras, V. Laude, Experimental study of band gaps and defect modes in a two-dimensional ultrasonic crystal, IEEE Ultrason. Symp. 377 (2003) 377–380. [5] Y. Pennec, B. Djafari-Rouhani, J.O. Vasseur, A. Khelif, P.A. Deymier, Tunable filtering and demultiplexing in phononic crystals with hollow cylinders, Phys. Rev. E 69 (2004) 046608. [6] B. Djafari-Rouhani, J.O. Vasseur, A.C. Hladky-Hennion, P. Deymier, F. Duval, B. Dubus, Y. Pennec, Absolute band gaps and waveguiding in free standing and supported phononic crystal slabs, Photonic Nanosrtuct. 6 (32) (2008) 32–37. [7] J.G. Hu, W. Xu, Band structures of phononic crystal composed of lattices with different periodic constants, Physica B 441 (2014) 89–93. [8] R. Wei, B. Wu, C. He, H. Zhao, Phononic band structure in a two-dimensional hybrid triangular graphite lattice, Physica B 404 (2009) 3795–3798. [9] M. Torres, F.R. Montero de Espinosa, Ultrasonic band gaps and negative refraction, Ultrasonics 42 (2004) 787–790. [10] Y.L. Xu, X.G. Tian, C.Q. Chen, Band structures of two dimensional solid/air hierarchical phononic crystals, Physica B 407 (2012) 1995–2001. [11] W. Zhao, Y.C. Zhao, L.B. Yuan, Characteristics of heterostructures in sandwich type of 2D phononic crystal, Physica Scr. T129 (2007) 112–115. [12] P.P. Huang, Y.W. Yao, F.G. Wu, et al., Interface-guided mode of Lamb waves in a two-dimensional phononic crystal plate, Chin. Phys. B 24 (2015) 054301. [13] Y.W. Yao, F.G. Wu, Z.L. Hou, Y.Y. Liu, Propagation properties of elastic waves in semi-infinite phononic crystals and related waveguides, Eur. Phys. J. B 58 (2007) 353–360. [14] Z.J. CHEN, Design and optimization of omnidirectional band gap for one-dimensional periodic and quasiperiodic phononic heterostructures, Chin. Phys. Lett. 32 (2015) 014301. [15] X.Y. Xiao, R.P. Chen, Acoustic band gap extension in one-dimensional solid/ fluid phononic crystal heterostructure, J. Comput. Acoust. 22 (2014) 1450010. [16] W.C. Chew, Q.H. Liu, Perfectly matched layer for elastodynamics: a new absorbing boundary condition, J. Comput. Acoust. 4 (1996) 341–359. [17] Y. Tanaka, T. Yano, S. Tamura, Surface guided waves in two-dimensional phononic crystals, Wave Motion 44 (2007) 501–512.
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Presentation and investigation of a new two dimensional heterostructure phononic crystal to obtain extended band gap Mohammad Bagheri Nouri n, Mehran Moradi Department of Mechanical Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran
art ic l e i nf o
a b s t r a c t
Article history: Received 1 January 2016 Received in revised form 15 February 2016 Accepted 17 February 2016 Available online 18 February 2016
In this paper, a new heterostructure phononic crystal is introduced. The new heterostructure is composed of square and rhombus phononic crystals. Using finite difference method, a displacement-based algorithm is presented to study elastic wave propagation in the phononic crystal. In contrast with conventional finite difference time domain method, at first by using constitutive equations and straindisplacement relations, elastic wave equations are derived based on displacement. Then, these forms are discretized using finite difference method. By this technique, components of stress tensor can be removed from the updating equations. Since the proposed method needs less elementary arithmetical operations, its computational cost is less than that of the conventional FDTD method. Using the presented displacement-based finite difference time domain algorithm, square phononic crystal, rhombus phononic crystal and the new heterostructure phononic crystal were analyzed. Comparison of transmission spectra of the new heterostructure phononic crystal with those creating lattices, showed that band gap can be extended by using the new structure. Also it was observed that by changing the angular constant of rhombus lattice, a new extended band gap can be achieved. & 2016 Elsevier B.V. All rights reserved.
Keywords: Heterostructure phononic crystal Band gap extension Finite difference time domain Displacement-based formulation
1. Introduction Much attention has been attracted to the phononic crystals in recent years. These non homogenous structures are periodic arrangement of inclusions in an elastically different host material. The frequency ranges that propagation of elastic waves is prohibited are called band gap. By this feature, these crystals can control propagation of elastic waves [1,2]. Since the capability of the crystals arise from phononic band gaps, it is essential to achieve the desired phononic band gaps. Acoustic filter and waveguide are some applications of phononic crystals. Phononic crystals are studied theoretically and experimentally in several researches [3–6]. To achieve the desired band gap, various works have been done on these crystals such as introduction of nesting phononic crystals [7], construction of a hybrid triangular lattice phononic crystal [8], presentation of twinned-square periodic structures [9] and introduction of hierarchical phononic crystals [10]. Heterostructure phononic crystals have been also investigated in some studies recently such as: presenting of a sandwich type of heterostructure in the 2D square lattice by replacing two columns of circular scatter cylinders with two columns of square scatter n
Corresponding author. E-mail address: [email protected] (M. Bagheri Nouri).
http://dx.doi.org/10.1016/j.physb.2016.02.023 0921-4526/& 2016 Elsevier B.V. All rights reserved.
cylinders [11], investigation of the interface-guided mode of Lamb waves in a phononic crystal heterostructure plate [12], analyzing of heterostructure phononic crystal waveguides [13], construction of an one-dimensional multilayer phononic heterostructure by combining finite periodic phononic crystals and Fibonacci quasiperiodic phononic crystals [14] and design of a phononic heterostructure by merging two one-dimensional phononic crystals [15]. By combining two-dimensional square and rhombus lattices, a new heterostructure phononic crystal is presented in this paper. Then, square and rhombus phononic crystals and the heterostructure phononic crystal were analyzed using displacementbased finite difference time domain method. Band gap of the new heterostructure phononic crystal was computed and compared with band gaps of its square and rhombus lattices. Effect of angular constant of the rhombus lattice on band gap of the heterostructure phononic crystal was studied too.
2. Model and calculation method 2.1. The model The proposed heterostructure phononic crystal has been shown in Fig. 1. The proposed heterostructure is composed of square and rhombus phononic crystals. This heterostructure is finite in horizontal direction and extends infinitely in vertical
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
29
u2l + 1/2, m + 1/2; n + 1 = u2l + 1/2, m + 1/2; n + v2l + 1/2, m + 1/2; n + 1/2 Δt
(8)
l + 1/2, m; n l + 1/2, m σ11 = C11 (u1l + 1, m; n − u1l, m; n )/Δx l + 1/2, m + C12 (u2l + 1/2, m + 1/2; n − u2l + 1/2, m − 1/2; n )/Δy
l, m + 1/2; n l, m + 1/2 ⎡ l, m + 1; n σ12 = C44 − u1l, m; n ⎣ u1
(
+ Fig. 1. The proposed hetero-structure phononic crystal that is composed of square arrangement (left) and rhombus arrangement (right).
a′ = a/(2 sin (γ /2))
(1)
where a′ and γ are constants of the rhombus lattice and a is constant of square lattice. 2.2. Calculation method Two-dimensional (2D) phononic crystals are considered here. The formulation can be easily derived for 3D phononic crystals. The inclusions are arranged along the z direction and are repeated in the xy plane. The equations of elastic waves can be written as:
ρvi̇ = σij, j
(2)
vi = u̇i
(3)
σij = Cijmn um, n
(4)
In above equations, ρ = ρ (x, y ), Cijkl (x, y ), vi and ui stand for the density, elastic stiffness tensor, i component of velocity and displacement of the structure, respectively. The summation agreement over dummy indices is considered, too. Since propagation of elastic waves in the xy plane is assumed, the displacement, velocity and stress tensor of the structure are independent of z, i.e., ui = ui (x, y, t ), vi = vi (x, y, t ) and σij = σij (x, y, t ). Consider that the inclusion and the host materials are isotropic. Then the discretized form of Eqs. (2–4) for study of the mixed mode is given as below: l + 1/2, m; n l − 1/2, m; n v1l, m; n + 1/2 = v1l, m; n − 1/2 + Δt /ρl, m ⎡⎣ σ11 − σ11
(
+
(
l, m + 1/2; n σ12
v2l + 1/2, m + 1/2; n + 1/2
=
−
l, m − 1/2; n σ12
v2l + 1/2, m + 1/2; n − 1/2
)
Δy⎤⎦
+
Δt /ρl + 1/2, m + 1/2
(5)
l + 1/2, m + 1; n l + 1/2, m; n × ⎡⎣ σ22 − σ22
(
(
l + 1, m + 1/2; n l, m + 1/2; n + σ12 − σ12
u1l, m; n + 1 = u1l, m; n + v1l, m; n + 1/2 Δt
) Δx
) Δy
) Δx⎤⎦
(6)
(7)
−
) Δy
u2l − 1/2, m + 1/2; n
(
) Δx⎤⎦
l + 1/2, m; n l + 1/2, m σ22 = C11 u2l + 1/2, m + 1/2; n − u2l + 1/2, m − 1/2; n
+ direction. The relation between constants of square and rhombus lattices can be written as:
(
u2l + 1/2, m + 1/2; n
l + 1/2, m C12
(
u1l + 1, m; n
−
(9)
u1l, m; n
) Δx
(10)
) Δy (11)
where (l,m) and n stands for location of a typical node and n-th time step, respectively. As Eq. (5) reveals, updating of x component of velocity at node (l,m) requires calculation of components of stress tensor at four l + 1/2, m; n l − 1/2, m; n l, m + 1/2; n l, m − 1/2; n coordinates, i.e., σ11 , σ11 , σ12 , σ12 . As Eqs. (9– 11) implies, computation of each of these stress components needs sizeable elementary arithmetical operations. To calculate transmission spectra of a phononic crystal, Fourier transformation of displacement data are taken. To compute transmission spectra of a phononic crystal, the calculated stress components have no direct contributions and are only used for updating components of displacement. By obtaining displacement-based forms of elastic wave equations and discretization of resultant equations, components of stress tensor can be eliminated from the updating equations. This course of action gives efficient updating equations that need less elementary arithmetical operations than conventional FDTD method. So the needed computational cost of the efficient updating equations is less than that of conventional FDTD method. The procedure of obtaining these updating equations is described below. For points far apart the interface of the inclusion and the host material (elastic constants of the structure are not discontinuous in these points), by substituting divergence of Eq. (4) in Eq. (2) the displacement-based form of elastic wave equation can be derived. By discretization of the resultant equations, displacement-based finite difference time domain (DBFDTD) formulation can be obtained. Consider that the inclusion and the host materials are isotropic. The displacement-based form of elastic wave equations for study of the mixed mode in this case can be given as below:
ρv1̇ = C11u1,11 + C44 u1,22 + ( C12 + C44 ) u2,12
(12)
ρv2̇ = C11u2,22 + ( C12 + C44 ) u1,12 + C44 u2,11
(13)
In Eqs. (12) and (13) C11, C12 and C44 are elastic constants of an isotropic material. To discretize Eqs. (12) and (13), all derivatives are replaced with central difference estimates. The consequent equations can be applied to update the velocity and displacement of the non interfacial grid points. For the interfacial grid points, the discretized form of the Eqs. (2–4) are used. For a suitable presentation of the discretized form of Eqs. (12) and (13), the subsequent coefficients are defined:
30
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
l, m λ1l, m = C11 /Δx2 l, m λ2l, m = C44 /Δy2
(
)
l, m l, m λ3l, m = C12 + C44 /(ΔxΔy)
Fig. 2. The computational model used for finite square phononic crystal analysis.
l + 1/2, m + 1/2 η1l + 1/2, m + 1/2 = C11 /Δy2 l + 1/2, m + 1/2 η2l + 1/2, m + 1/2 = C44 /Δx2
(
l + 1/2, m + 1/2 l + 1/2, m + 1/2 η3l + 1/2, m + 1/2 = C12 + C44
) (ΔxΔy)
(14)
The coefficients introduced by Eq. (14) may be calculated before the time evolution loop of the DBFDTD method and saved in appropriate matrixes in order to avoid computing them in each time step. Discretizing Eq. (12) using the finite difference technique leads to:
v1l, m; n + 1/2 = v1l, m; n − 1/2 + Δt/ρl, m Σ
(15)
where:
(
)
(
Σ = − 2 λ1l, m + λ2l, m u1l, m; n + λ1l, m u1l + 1, m; n + u1l − 1, m; n +
λ2l, m
+
λ3l, m
( (u
u1l, m + 1; n
+
u1l, m − 1; n
l + 1/2, m + 1/2; n 2
−
Fig. 3. Transmission spectra through the square phononic crystal obtained by DBFDTD calculation along the ГX direction: for longitudinal (solid line) and transverse polarization (dashed line) of incident waves.
)
)
u2l + 1/2, m − 1/2; n
− u2l − 1/2, m + 1/2; n + u2l − 1/2, m − 1/2; n
)
(16)
Discretized form of Eq. (13) can be given as:
v2l + 1/2, m + 1/2; n + 1/2
= v2l + 1/2, m + 1/2; n − 1/2 + Δt/ρl + 1/2, m + 1/2 χ
(17)
where:
(
)
χ = − 2 η1l + 1/2, m + 1/2 + η2l + 1/2, m + 1/2 u2l + 1/2, m + 1/2; n
( (u (u
+ η1l + 1/2, m + 1/2 u2l + 1/2, m + 3/2; n + u2l + 1/2, m − 1/2; n +
η2l + 1/2, m + 1/2
+
η3l + 1/2, m + 1/2
l + 3/2, m + 1/2; n 2
l + 1, m + 1; n 1
+
u2l − 1/2, m + 1/2; n
) )
− u1l + 1, m; n − u1l, m + 1; n + u1l, m; n
)
(18)
Subsequently, Eqs. (7) and (8) can be used to update components of displacement of the structure. The efficiency of DBFDTD formulation can be easily understood by comparison of Eqs. (15) and (16) with Eqs. (5), (9) and (10). This comparison shows that Eqs. (15) and (16) have about 50 percent less elementary arithmetical operations.
3. Results and discussion In following examples, the inclusion and the matrix are steel and epoxy, respectively. The density and elastic constants C11 and C44 of steel are assumed to be 7780 kg/m3, 264 GPa and 81 GPa respectively, and those for epoxy are 1142 kg/m3, 7.54 GPa and 1.48 GPa. To validate result of the DBFDTD method, transmission spectra of a square phononic crystal was computed and compared with result of the FDTD method. Diameter of steel cylinders and the lattice constant are 5 and 7 mm, respectively. To calculate transmission spectra of the crystal, a sample with five unit cells in x direction and one unit cell in y direction was oriented along ГX direction and surrounded by two homogenous regions (Fig. 2). The perfectly matched layer (PML) [16] was adopted as the absorbing boundary condition for the left and right boundaries of the computational domain. For the top and bottom boundaries the periodic boundary condition was used.
A Gaussian wave packet was launched along the x direction at the left homogenous region. The x component (y component) of displacement was recorded at middle of right homogenous region and averaged over a lattice constant. By taking fast Fourier transform of this average, transmission spectra through the crystal for longitudinal (transverse) polarization of incident waves was obtained. Fig. 3 exhibits DBFDTD results of the transmission spectra through the crystal. To obtain the transmission spectra, a grid of 60*60 points in the unit cell was considered and the equations of motion were solved over 216 time steps with each time step lasting 7.08 ns. Fig. 3 shows a normalized frequency range extending from 3.8 to 8.7 that the crystal forbids propagation of elastic waves along ГX direction for any polarization of incident waves. Normalized frequency is defined as ωa/vt where ω , a and vt are respectively angular frequency, the lattice constant and transverse velocity of elastic waves in epoxy. This means that there is a band gap along ГX direction extending from 3.8 to 8.7. The computed band gap is in good agreement with that calculated by the FDTD method (3.8–8.5) [17]. 3.1. The heterostructure composed of square and triangular lattices First, the heterostructure composed of square and triangular phononic crystals was studied. Diameter of steel rods is 4 mm. Lattice constant of both square and triangular Lattices are 6 mm. Fig. 4 shows the sample was used to calculate transmission spectra of the heterostructure. The sample was surrounded by two homogenous regions (not shown here). The PML was used as the absorbing boundary condition for the boundaries along the x direction. For the boundaries along the y direction the periodic boundary condition was applied. A Gaussian wave packet was launched along the x direction at the left homogenous region. The x component of displacement was recorded at middle of right homogenous region and averaged over a lattice constant. By taking fast Fourier transform of this
Fig. 4. The computational model used for heterostructure phononic crystal analysis.
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
31
Fig. 5. Transmission spectra through the square lattice (dot line), the triangular lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction computed for longitudinal polarization of incident wave. The longitudinal component Ux of the displacement is given in arbitrary units.
Fig. 7. Transmission spectra of square lattice (dot line), rhombus lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction of propagation computed for longitudinal incident wave. The longitudinal component Ux of the output displacement field is given in arbitrary units.
average, transmission spectra through the heterostructure was obtained. Fig. 5 exhibits DBFDTD results of the transmission spectra through the heterostructure for longitudinal polarization of incident waves. To obtain the transmission spectra, a grid of 60*60 points in the unit cell of square lattice was considered and the equations of motion were solved over 216 time steps with each time step lasting 6.09 ns. The CPU times to calculate transmission spectra are 2 min and 16 s for the DBFDTD method and 3 min and 34 s for the conventional FDTD method on twenty 3 GHz Intel Xeon CPUs. So the computations time of the DBFDTD method is 36 percents less than that of the conventional FDTD method. The numerical examples were calculated by a developed FORTRAN code based on parallel processing. Fig. 5 shows that there is a band gap extending from 110 to 238 kHz that the square phononic crystal forbids propagation of elastic waves along the x direction. The band gaps of the triangular phononic crystal locate at 120–288 and 313–367 kHz. In the frequency range of Fig. 5, two band gaps were found for the heterostructure, i.e., 113–288 and 313-365 kHz. As Fig. 5 indicates the heterostructure only allows propagation of the waves that both the square and triangular lattices allow this propagation. So the band gap can be extended by combination of square and triangular lattices and construction of the heterostructure. Due to proximity of two constituent lattices’ band gap ranges, heterostructure transmission spectra are approximately same as triangular lattice spectra.
and 400 kHz for the rhombus phononic crystal and the band gap of the heterostructure locates at 110–400 kHz. Fig. 7 implies that by combination of the square and the rhombus lattices and creation of the heterostructure, a new extended band gap can be achieved. Comparison of these transmission spectra shows that heterostructure band gap has been extended extensively. As the last example, the heterostructure composed of square and rhombus lattices was considered. The angular constant of the rhombus phononic crystal is 30° here. The diameter of steel rods and the lattice constant of the square phononic crystal are 4 and 6 mm, respectively. According to Eq. (1), the lattice constant of rhombus phononic crystal is 11.59 mm. Fig. 8 shows the sample was used to calculate transmission spectra of the heterostructure. The transmission spectra through the heterostructure along x direction for longitudinal polarization of incident waves were shown in Fig. 9. To calculate the spectra, the same procedure of the previous examples was complied. In the frequency range of Fig. 9, five band gaps appear for the rhombus phononic crystal, i.e., 67–148, 173–246, 274–282, 341– 362 and 382–387 kHz. Also there are five band gaps for the heterostructure in the frequency domain of Fig. 9, i.e., 68–246, 274– 286, 339–362 and 382–387 kHz. As Fig. 9 indicates, a new different extended band gap can be attained using the heterostructure composed of the square and the rhombus lattices. According to the obtained results, the band gap of the heterostructure phononic crystal would be the wider than those of other studied heterostructures if the angular constant of the rhombus lattice is 80°.
3.2. The heterostructure composed of square and rhombus lattices 4. Conclusion As second example, the heterostructure composed of square and rhombus lattices was considered. The angular constant of the rhombus phononic crystal is 80°. The diameter of steel rods and the lattice constant of the square phononic crystal are 4 and 6 mm, respectively. According to Eq. (1), the lattice constant of rhombus phononic crystal is 4.67 mm. Fig. 6 shows the sample was used to calculate transmission spectra of the heterostructure. The transmission spectra through the heterostructure along x direction for longitudinal polarization of incident waves were shown in Fig. 7. To obtain the spectra, the same procedure of the previous example was obeyed. In the frequency range of Fig. 7, a band gap occurs between 183
Fig. 6. same as Fig. 4 but here, the triangular lattice was replaced with a rhombus lattice whose angular constant is 80°.
In this paper a new heterostructure phononic crystal was presented. The heterostructure was composed of the square and the rhombus phononic crystals. Also displacement-based finite difference time domain method was introduced and applied for simulation of elastic waves propagation through phononic crystals. Comparison of the transmission spectra computed by the DBFDTD method with the results of FDTD method, proved the efficiency of the DBFDTD method. Also Comparison of the computational time showed that the computational time of the DBFDTD method is 36 percents less than that of the FDTD method. Propagation of elastic waves through the proposed heterostructure phononic crystal was studied using the DBFDTD method. Investigation of three calculated examples demonstrated that the
Fig. 8. same as Fig. 6 but the angular constant of rhombus lattice is 30° here.
32
M. Bagheri Nouri, M. Moradi / Physica B 489 (2016) 28–32
Fig. 9. Transmission spectra of square lattice (dot line), rhombus lattice (dash line) and the heterostructure phononic crystal (solid line) along the x direction of propagation computed for longitudinal incident wave. The longitudinal component Ux of the output displacement field is given in arbitrary units.
band gap can be extended using the heterostructure phononic crystal. Effect of the angular constant of the rhombus phononic crystal on the band gap of the heterostructure was studied too. This investigation showed that the band gap of the heterostructure phononic crystal would be the wider than those of other studied heterostructures if the angular constant of the rhombus lattice is 80°.
References [1] M.S. Kushwaha, P. Halevi, L. Dobrzynski, B. Djafari-Rouhani, Acoustic band structure of periodic elastic composites, Phys. Rev. Lett. 71 (1993) 2022–2025. [2] R. Martinez-Sala, J. Sancho, J.V. Sanchez, V. Gomez, J. Llinares, et al., Sound attenuation by sculpture, Nature 378 (1995) 241-241.
[3] F.R. Montero de Espinosa, E. Jime´nez, M. Torres, Ultrasonic band gap in a periodic two-dimensional composite, Phys. Rev. Lett. 80 (1998) 1208–1211. [4] A. Khelif, A. Choujaa, R. laihem, M. Wilm, S. Ballandras, V. Laude, Experimental study of band gaps and defect modes in a two-dimensional ultrasonic crystal, IEEE Ultrason. Symp. 377 (2003) 377–380. [5] Y. Pennec, B. Djafari-Rouhani, J.O. Vasseur, A. Khelif, P.A. Deymier, Tunable filtering and demultiplexing in phononic crystals with hollow cylinders, Phys. Rev. E 69 (2004) 046608. [6] B. Djafari-Rouhani, J.O. Vasseur, A.C. Hladky-Hennion, P. Deymier, F. Duval, B. Dubus, Y. Pennec, Absolute band gaps and waveguiding in free standing and supported phononic crystal slabs, Photonic Nanosrtuct. 6 (32) (2008) 32–37. [7] J.G. Hu, W. Xu, Band structures of phononic crystal composed of lattices with different periodic constants, Physica B 441 (2014) 89–93. [8] R. Wei, B. Wu, C. He, H. Zhao, Phononic band structure in a two-dimensional hybrid triangular graphite lattice, Physica B 404 (2009) 3795–3798. [9] M. Torres, F.R. Montero de Espinosa, Ultrasonic band gaps and negative refraction, Ultrasonics 42 (2004) 787–790. [10] Y.L. Xu, X.G. Tian, C.Q. Chen, Band structures of two dimensional solid/air hierarchical phononic crystals, Physica B 407 (2012) 1995–2001. [11] W. Zhao, Y.C. Zhao, L.B. Yuan, Characteristics of heterostructures in sandwich type of 2D phononic crystal, Physica Scr. T129 (2007) 112–115. [12] P.P. Huang, Y.W. Yao, F.G. Wu, et al., Interface-guided mode of Lamb waves in a two-dimensional phononic crystal plate, Chin. Phys. B 24 (2015) 054301. [13] Y.W. Yao, F.G. Wu, Z.L. Hou, Y.Y. Liu, Propagation properties of elastic waves in semi-infinite phononic crystals and related waveguides, Eur. Phys. J. B 58 (2007) 353–360. [14] Z.J. CHEN, Design and optimization of omnidirectional band gap for one-dimensional periodic and quasiperiodic phononic heterostructures, Chin. Phys. Lett. 32 (2015) 014301. [15] X.Y. Xiao, R.P. Chen, Acoustic band gap extension in one-dimensional solid/ fluid phononic crystal heterostructure, J. Comput. Acoust. 22 (2014) 1450010. [16] W.C. Chew, Q.H. Liu, Perfectly matched layer for elastodynamics: a new absorbing boundary condition, J. Comput. Acoust. 4 (1996) 341–359. [17] Y. Tanaka, T. Yano, S. Tamura, Surface guided waves in two-dimensional phononic crystals, Wave Motion 44 (2007) 501–512.