Bloch-wave expansion technique for predicting wave reflection and transmission in two-dimensional phononic crystals Jason A. Kulpe, Karim G. Sabra,a) and Michael J. Leamy School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332

(Received 4 October 2013; revised 10 January 2014; accepted 21 January 2014) In this paper acoustic wave reflection and transmission are studied at the interface between a phononic crystal (PC) and a homogeneous medium using a Bloch wave expansion technique. A finite element analysis of the PC yields the requisite dispersion relationships and a complete set of Bloch waves, which in turn are employed to expand the transmitted pressure field. A solution for the reflected and transmitted wave fields is then obtained using continuity conditions at the half-space interface. The method introduces a group velocity criterion for Bloch wave selection, which when not enforced, is shown to yield non-physical results. Following development, the approach is applied to example PCs and results are compared to detailed numerical solutions, yielding very good agreement. The approach is also employed to uncover bands of incidence angles whereby perfect acoustic reflection from the PC occurs, even for frequencies outside of stop bands. C 2014 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4864457] V PACS number(s): 43.35.Gk, 43.20.El, 43.20.Fn [ANN]

I. INTRODUCTION

Periodic arrays of scattering objects, embedded in a homogeneous host medium, constitute a phononic crystal (PC). Multiple scattering effects combined with the resonances of the scattering objects lead to a complex dispersion relationship and the possibility of frequency bandgaps and negative refraction.1,2 A variety of PCs have been studied such as periodic arrays of fluid cylinders3,4 or spheres5,6 in a water background. Waveguides with periodic profiles,7 locally resonant periodic structures,8 and PCs designed to filter acoustic signals9 are examples of structures engineered to exploit the frequency bandgap behavior. Regardless of the system, knowledge of the dispersion relationship and associated acoustic wave fields are crucial for characterizing material and/or device response. A multitude of methods have been developed for dispersion computations, often focusing on wave dynamics within a single unit cell. Plane wave expansion methods, such as those detailed in Ref. 4, use Fourier expansions for the wave field and periodic material properties to compute the dispersion behavior, yet can experience convergence issues for highly contrasting material properties.9 Chen and Ye3 and Krynkin et al.10 considered spherical harmonics via separation of variables to compute the dispersion relationship for given fluid, and respectively, fluid-solid crystals. Alternatively, transfer matrices7 are well-suited for one-dimensional systems whereas finite difference techniques can also be used to numerically compute the dispersion of complex or finite periodic structures.9,11 Homogenization techniques12,13 are limited to low frequency systems and thus will not be considered here. Finite element (FE) based techniques14–17 allow calculation of the dispersion and acoustic wave field for arbitrary material and geometric configurations. The Bloch theorem15,18 a)

Author to whom correspondence should be addressed. Electronic mail: [email protected]

1808

J. Acoust. Soc. Am. 135 (4), April 2014

Pages: 1808–1819

allows one to consider a single unit cell of the periodic medium thereby enabling direct calculation of the dispersion behavior and a set of Bloch waves which describe the wave field within the entire medium. In the past, authors have discussed the interaction of an external incident wave with a PC, to understand the effects of multiple scattering and PC dispersion on the wave field. The layer Korringa-Kohn-Rostoker method1,19 is a common approximation method in the photonics literature where slabs of scattering objects, of infinite height, are correlated using specialized scattering matrix methods. Similar layer scattering methods have also been considered in elastic PCs.20,21 Leroy et al.22,23 considered a layer method for modeling the reflection coefficients from successive thin layers three-dimensional (3D) of bubbles under low frequency limits; their analysis is not applicable here as we only consider two-dimensional (2D) PCs at arbitrary frequencies. Moiseyenko et al.24,25 conducted numerical FE studies that focused on diffraction of the incident wave but without a quantitative independent verification study of the pressure field results. Sanchis et al.26 compared pressure fields using doubly multiple scattering theory versus experiment for an array of 35 aluminum cylinders periodically spaced in air. The work in Ref. 26, which did not exploit the periodicity of the structure via the Bloch theorem, is the only work comparing analytical pressure fields to an independent verification study known to the authors. Laude et al.27,28 have surmounted the reflection problem by considering an expansion in terms of the Bloch eigenmodes of the PC system, at a given wave incidence angle and frequency. Laude et al. also considered wave “deafness”27 by the symmetrical or asymmetrical form of the Bloch wave on the boundary, a topic discussed in the excitability of particular Bloch waves. Laude et al.27 did not, however, validate their expansion approach of the wave fields with an independent verification study nor explicitly include a crucial group velocity selection criterion.

0001-4966/2014/135(4)/1808/12/$30.00

C 2014 Acoustical Society of America V

The motivation of this paper is to derive the equations governing the dispersion and reflection/transmission, starting with a Bloch mode expansion for a lossless PC given any frequency and incidence wave angle. This research contributes the following to the predictions of reflection and transmission of waves incident upon a PC: (1) A Bloch wave selection criterion based on group velocity considerations; (2) formulation of a dispersion technique, implemented using commercial FE software, designed for direct applicability to a Bloch wave expansion; and (3) complete solution and verification of the acoustic wave fields for several frequencies and incidence angles. II. MODEL AND DISPERSION RELATION

In a 2D periodic PC the Bloch theorem18 states that everywhere in the medium the pressure field p(x, t) takes the form of a function p~(x), hereby termed a wave mode, that is periodic on the unit cell and modulated by a harmonic plane wave of frequency x pðx; tÞ ¼ p~ðxÞeiðkxxtÞ :

(1)

In fact, every field quantity in the PC obeys a similar relationship. The Bloch wave vector k is related to x by a

dispersion relationship. The 2D PC under study is displayed in Fig. 1(a) and the unit cell in Fig. 1(b) where the lengths aj (j ¼ 1,2) denote the dimensions of the unit cell, R is the radius of the cell-centered inclusion, q is the density, and c is the sound speed of the background (subscript 1) and inclusion material (subscript 2). For the sake of simplicity, we only consider 2D square unit cells (a1 ¼ a2 ¼ a) with cell centered inclusions in this work, which could be extended in a straight forward manner to non-square unit cells through a lattice transformation.12 Periodicity of the PC imposes periodicity of the real wavenumber space (reciprocal k-space), shown in Fig. 1(c), where bj ¼ 2p/aj is the length of the cell in k-space. For notational clarity let subscript j ¼ 1,2 (or x, y) correspond to a length or vector component. Unique Bloch waves are found with real part components in the range p/aj  kj  p/aj, an extent termed the first Brillouin zone, where kj is the jth component of the wavevector. However, the imaginary components of the wavevector are unbounded.28,29 For clarity we define propagating Bloch waves as having a purely real wavevector, evanescent Bloch waves as having a purely imaginary wavevector, and decaying Bloch waves as having complex wavevectors. Implicit in this field representation is that there is an infinite number of Bloch waves p(x, t) each with corresponding wavevector k that satisfy the Bloch theorem of Eq. (1) and the acoustic

FIG. 1. (Color online) (a) A 2D PC. (b) A single unit cell is highlighted with periodicity a1 and a2 in the e^x, e^y, directions, respectively, and a cell-centered circular inclusion. (c) The first Brillouin zone, in wavenumber k-space, for a 2D unit cell. (d) A representative FE mesh of a single unit cell with reduced nodal degrees of freedom depicted as dots. The nodal degrees of freedom represented by squares are related by periodic boundary conditions via matrix T. pi is one node contained within the set pi.

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

Kulpe et al.: Bloch-wave expansion for phononic crystals

1809

wave equation. Hence, at a fixed frequency and direction in k-space there are a finite number of propagating waves, due to the limitations imposed by the Brillouin zone, yet an infinite set of decaying and evanescent waves. The complex waves of a PC constitute a complete basis and justify a wave expansion; see the Appendix for more details. Decaying and evanescent Bloch waves are necessary near the boundary of the incident medium and account for the near field wave refraction behavior.28 Traditional methods for solving the dispersion relation for a particular PC involve specifying real k and finding x as an eigenvalue. The exponentially decaying and evanescent Bloch waves and the associated band structure are not revealed by these x(k) methods,30 as k is restricted to be real. A different k(x) method will be introduced, for complex k, with specific application to wave transmission into the PC half-space. A. Derivation of dispersion technique

We next present a method which uses the finite element method (FEM) to form equations leaving a single parameter k as a complex eigenvalue and the discretized mode vector p as an eigenvector for a particular frequency x. The method presented here is formulated specifically to easily obtain assembled FEM matrices from commercial software, such as COMSOL MultiphysicsV. Through the FEM, this dispersion technique remains quite general and can adequately accommodate an arbitrary fluid inclusion. We start with the time dependent wave equation, with periodically varying density q(x) and sound speed c(x)   @ 2 pðx; tÞ 1 1 rp x; t  r  (2) ð Þ ¼ 0: qðxÞ qðxÞcðxÞ2 @t2 R

In this paper, the functions q(x) and c(x) are equal to constant values q1 and c1 for the water host (q1 ¼ 1000 kg/m3 and c1 ¼ 1500 m/s) and q2, c2 and for the circular inclusion with radius R, as indicated in Fig. 1(b). Insertion of the Bloch theorem Eq. (1) and simplification results in    1 1 r  ðr~ p þ ik~ p Þ  ik  r~ p  ðk  kÞp~ q q 

x2 p~ ¼ 0: qc2

(3)

We seek to formulate Eq. (3) with a FE model, using triangular quadratic Lagrange elements. Multiplication by a test function w(x) and integration over the unit cell results in   ð  1 rw  ðr~ p þ ik~ pÞ q cell  x2 w p  ðk  kÞp~  2 w~ p dV ¼ 0: (4)  ik  r~ qc q In anticipation of the follow-on development we zeroed, without loss of generality, the right-hand side boundary integral terms in Eq. (4) resulting from application of the divergence theorem [see discussion following Eq. (8)]. Tessellating the domain with triangular elements and invoking a Galerkin 1810

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

discretization of the pressure field p~ and test function w(x) ! wm(x) using quadratic Lagrange interpolation functions results in five elemental matrices with components ð 1 element ¼ w w dV (5a) Mmn 2 m n V qc ð 1 element rwm  rwn dV ¼ (5b) Kmn q V ð 1 Belement wm wn dV ¼ (5c) mn q V  ð  1 @wm @wn element  wm dV (5d) wn Amn;x ¼ @x @x Vq  ð  1 @wm @wn  w dV: (5e) w Aelement ¼ n m mn;y @y @y Vq Here m, n ¼ 1 to 6 for quadratic triangular elements, corresponding to the individually numbered nodes of an element, and V is the volume of an element. The six interpolation functions wm can be found in a typical FE textbook.31 Assembly of the elemental matrices [Eq. (5)] over all elements in the domain results in five assembled global matrices M, K, B, Ax, Ay and the equations of motion read

2  (6) x MþKþikx Ax þ iky Ay þ ðk  kÞB P ¼ 0: Dispersion techniques leaving the periodic mode vector p as an eigenvector require periodic boundary conditions.27 Unit cell periodicity implies that the nodal degrees of freedom on the bottom “B” are equal to the nodal degrees of freedom on the top “T” (pB ¼ pT), the left “L” to the right “R,” etc. [see Fig. 1(d)]. Implementation of periodic boundary conditions is easily achieved using a transformation matrix T17,32,33 by preferential ordering of the nodal degrees of freedom in the global nodal vector p, with pI referring to all nodes on the mesh interior. The reduced set of unique nodal ~ (note: the choice of the values is contained in the vector p left and bottom degrees of freedom in the reduced set is arbitrary) which is related to the vector p by p ¼ T~ p 2 3 pL 6p 7 7 ~¼6 p 6 B 7 4 pLB 5 pI 2

I 6I 6 6 60 6 60 6 T¼6 60 60 6 6 60 6 40 0

0 0 I I 0 0 0 0 0

(7a)

(7b)

0 0 0 0 I I I I 0

3 0 07 7 7 07 7 07 7 07 7; 07 7 7 07 7 05 I

(7c)

Kulpe et al.: Bloch-wave expansion for phononic crystals

where I and 0 are identity and zero block matrices. We impose the periodic boundary conditions on the equations of motion in Eq. (6) by premultiplication of TT and substitution of Eq. (7) yielding

~ Kþik ~ ~ ~ ~ ~ ¼ 0; x2 Mþ (8) x A x þ iky A y þ ðk  kÞB p ~ ¼ TTKT, etc., and the superscript T ~ ¼ TTMT, K where M denotes a matrix transpose. The FEM right-hand side boundary terms for the cell, omitted in Eq. (4), equate to zero after premultiplication by the matrix TT.33 Transformation of Eq. (8) to a scalar eigenvalue problem requires a wavevector parameterization leaving a single parameter k as an eigenvalue. To this end, we write k as a k ¼ k0 þ k^

(9)

and leave the vectors k0 and a^ as problem specific parameters; see Fig. 1(c) and Ref. 30. This general scheme allows the complex dispersion relationship to be calculated along any direction a^ in k-space. Inserting Eq. (9) and manipulation of the equations of motion Eq. (8) results in a quadratic eigenvalue problem for the dispersion relationship with eigenvalue ~ k and eigenvector (corresponding to the Bloch mode) p

0  ~ ¼ 0; (10) D ðxÞ þ kA0 þk2 B0 p where ~ K ~ ~ ðxÞ ¼ x2 Mþ D

~ ðxÞ þ i k0x A ~ x þ k0y A ~ y þ k2 B ~ D 0 ðx Þ ¼ D 0

(11b)



~ x þ ay A ~y A0 ¼ i ax A

(11c)

~ B0 ¼ a2 B:

(11d)

(11a)

~ K, ~ and B ~ are symmetric and A ~ j is Note that the matrices M, not, therefore the eigenvalues of Eq. (10) are not guaranteed to be real and positive. Numerical solution efficiency is obtained by converting the quadratic eigenvalue problem, of size NR  NR [NR being the size of the square matrices appearing in Eq. (10)] to a generalized eigenvalue problem ~ ¼ k~ p and writing of size 2NR  2NR by introducing u   ~ p ðS  kGÞ ¼0 (12) ~ u   0 I S¼ (13a) D0 ðxÞ A0   I 0 G¼ : (13b) 0 B0 The FE matrices are typically sparse and efficiency is obtained by using any one of the sparse eigenvalue solution routines. Solution of the discretized generalized eigenvalue problem Eq. (12) leads to a finite set of periodic Bloch modes and, possibly, complex wavenumbers k. The Bloch J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

wavevector k is reconstituted from the eigenvalue k by Eq. (9). An iterative solution for many values of x will lead to a k(x) dispersion relation for the arbitrary 2D PC, in any direction a^. The eigenvalue problem previously outlined is perfectly general, given the vectors k0 and a^, and can be tailored to identify Bloch waves for a transmitted field expansion. The geometry of the reflection-transmission problem is shown in Fig. 2 where a semi-infinite 2D PC half-space is adjacent to a homogeneous water medium with density qi and sound speed ci. For reflection-transmission problems the parallel wavevector component must be conserved along the PC interface.18,27 We require the set of Bloch waves used in our expansion to have the parallel component equal to a specified constant. Wavenumber conservation along the interface, given by unit vector ^e y, will be maintained by searching for Bloch waves along a direction a^ ¼ ^e x.30 Note, regardless of the interface’s orientation with respect to the coordinate axes, the vector a^ is always parallel to the PC interface normal. The vector k0 is composed of the parallel component of the incident wavevector, e.g., k0 ¼ [0 x/ci sin(h)]T where h is the incidence angle with respect to the interface normal ^e x. Note k  ^e y ¼ const 8 k and therefore the relationship allows for wavenumber conservation along the PC interface. Geometrically this constraint implies searching for Bloch waves along a line, given by Eq. (9) and at fixed x, through the 3D complex dispersion surface. See Ref. 18 for a further geometrical discussion of the wavevector relationship. Calculation of the group velocity for a particular Bloch wave can be directly obtained from the eigenvalue problem of Eq. (12). First we write Eq. (12) as W(x)q ¼ kq where ~ ]T is a right eigenvector of p u W(x) ¼ G1S(x) and q ¼ [~ W. The derivative of the nth eigenvalue kn with respect to x can be obtained34 as @kn @WðxÞ ¼ hH qn ; n @x @x

(14)

where now the vector hH n is the nth left eigenvector of W and q is scaled such that hH n n ¼ 1. Here the superscript H denotes the matrix Hermitian operation. Finally, the group velocity

FIG. 2. (Color online) A half-space consisting of a homogeneous water medium and a PC. Three reflected wave orders and the incident plane wave are included for an example. Kulpe et al.: Bloch-wave expansion for phononic crystals

1811

for the nth Bloch wave into the PC medium cgn  e^x is obtained by "  # 1 @kn cgn ^e x ¼Re @x 20 2 3 11 3 0 0 6B 6 7 C 7 6 H6 7q n C 7 ; (15a) ¼Re6B h 5 A 7 4@ n 4 01 @ 0 5 B D ðxÞ 0 @x 2 @ 0 ~ y þ 2x sin h B ~ 0: ~ þ i sin h A D ðxÞ ¼ 2xM 2 @x ci ci

(15b)

With the two preceding equations the group velocity can be efficiently determined after solution of the eigenvalue problem in Eq. (12). B. Verification of dispersion technique

Up to this point we have formulated a technique specifically for easy implementation in commercial FE programs. To verify the preceding dispersion technique the system parame3 a2ffiffi¼ ters are set to {a1 ¼ p ffi 1 m,T R ¼ 0.4 m, q2 ¼ 7800 kg/m , c2 ¼ 4591 m/s, a^ ¼ 1/ 2[1 1] , k0 ¼ 0}. This parameter choice allows direct comparison for the complex dispersion relationship from Laude et al. (see Fig. 6 in their work).27 Inspection of the results in Fig. 3 show that the dispersion curves calculated from the present technique agree very well with that reported in the literature. In the dispersion diagram of Fig. 3 we display the wavenumbers (as small dots) according to the real and imaginary parts. The results here indicate most pffiffiffi decaying wavenumbers have a real component equal to 2p=a (denoted as the M point in the first Brillouin zone18). At frequencies over 16 000 rad/s there is a branch with wavenumbers of complex values; these are wavenumbers not easily identified by a conventional x(k) search method. III. BLOCH WAVE EXPANSION

In this section we detail the use of a Bloch wave expansion to model the transmitted wave field in a PC half-space. The dispersion technique discussed in Sec. II will be used to identify the Bloch waves needed to compose the expansion for a particular choice of frequency and incidence wave angle. The expansion equations and solution procedure are taken from the work of Ref. 27 and are repeated here for completeness. However, we add an additional constraint pertaining to wave group velocity for mode selection. We also add discussion of the mode selection and associated energy calculations quantifying the amount of acoustic energy reflected from the PC half-space. The PC considered consists of air inclusions with R ¼ 2 cm, q2 ¼ 1.22 kg/m3, and c2 ¼ 343 m/s in a water host medium with a1 ¼ a2 ¼ 40 cm. The physical parameters chosen correspond to the parameters discussed in Ref. 35. Consider an incident harmonic plane wave pi(x, t) propagating at frequency x in a homogeneous incident medium with density and speed (qi ¼ 1500 m/s and ci ¼ 1000 kg/m3) and with wavevector ki and angle h relative to the PC normal, with h ¼ 0 corresponding to normal incidence. Standard 1812

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

FIG. 3. (Color online) Verification of the complex dispersion relation for a circular steel rod in water. The large dots represent the dispersion relationship from Ref. 27 (see Fig. 6 in their work). The real wavenumbers, imaginary wavenumbers, and complex wavenumbers (small dots) were found by the dispersion technique demonstrated in this paper.

boundary conditions for reflection problems apply: conservation of frequency and parallel wavevector, and the matching of pressure and normal particle velocity at the interface.36 Recall the discussion in Sec. II A pertaining to discussion of the choice of the parallel wavevector. To begin identifying the Bloch waves that compose the transmitted wave field we first solve the dispersion relationship, outlined previously, for a particular unit cell giving a set of Bloch wavevectors and modes. Solving Eq. (12) produces many eigenvalues and eigenvectors to choose from. A few notes are in order regarding the numerical solution of Eq. (12). First, a k(x) method does not automatically obey periodicity of the reciprocal k-space.30 We can identify unique Bloch waves by selecting waves with the real part within the first Brillouin zone: p/a  Re k  p/a;30 no constraint is placed on the magnitude of the imaginary part of the wavevector. Second, eigenvalues k which would produce an exponentially growing Bloch wave are discarded. Third, all propagating waves are chosen such that the wave’s group velocity vector cg satisfies cg  ^e x > 0:

(16)

The group velocity vector always lies in the same direction as the acoustic intensity vector averaged over the unit cell.37 Therefore, the constraint in Eq. (16) implies that at steadystate acoustic power flows into the PC medium.18 Identified propagating Bloch waves that do not satisfy this third criteria are discarded. This group velocity condition places a wavenumber selection criterion on 6k for some real wavenumber k. Symmetry of the dispersion relation about the vertical frequency axis implies a choice, based on group velocity, of real wavenumbers from a pair of wavenumbers. Through the discretized eigenvalue problem, after considering the three solution criteria, we are left with N Bloch waves and wavevectors; often in practice N  NR where NR is the number of degrees of freedom from the FEM Kulpe et al.: Bloch-wave expansion for phononic crystals

discretization in Eq. (8). The following outline is helpful in summarizing the method of Bloch wave selection: (1) Compute the dispersion relation, as defined in Sec. II, for a chosen x and h. (2) Numerically evaluate Eq. (12) for Bloch wavevectors and waves and Eq. (15) for determining group velocity into the PC half-space. (3) Remove very small erroneous real and imaginary wavenumber parts, left as numerical artifacts. (4) Select any and all propagating Bloch waves with cg  ^e x > 0 and p/a  Re k  p/a. (5) Select any and all decaying and evanescent Bloch waves, subject to a predefined acceptable maximum decay rate, such that the wave exponentially decays in a direction away from the interface. This set of waves constitute the Bloch expansion set. As an example of the Bloch wave selection, consider identifying suitable Bloch waves for the expansion at normal incidence at three frequencies. As described in Sec. II A Bloch waves with a conserved wavenumber can be found along a horizontal line at frequency x, parameterized by k. Fig. 4 depicts the complex dispersion diagram at normal incidence with lines drawn at x ¼ 5000, 10 000, and 18 000 rad/s and the selected Bloch waves indicated by squares. At 5000 rad/s there are no propagating Bloch waves and only decaying and evanescent Bloch waves constitute the expansion basis. At 10 000 rad/s there are two propagating Bloch waves (with wavenumbers at 6k), one being selected (square) and the other being discarded () as this Bloch wave’s group velocity is directed out of the PC half-space and violates the group velocity selection criterion. Similarly, at 18 000 rad/s the propagating Bloch wave on the Re k < 0 side of the diagram is chosen as this wave satisfies the group velocity criterion. In practice, additional decaying and

evanescent waves would be desired for an accurate expansion; the imaginary wavenumber diagram here is truncated to 20 rad/m for clarity. The above solution procedure applies equally to a variety of incidence angles, each with a corresponding complex dispersion diagram, as further described in Sec. IV. We next write the pressure field in the PC half-space, pt(x, t), as a summation of N Bloch waves pt ðx; tÞ ¼

N X

tn pn ðx; tÞ;

(17)

n¼1

where the modal transmission coefficients tn are obtained based on the interface boundary conditions. Similarly, the reflected pressure field, pr(x, t), is written as an expansion of N plane waves with unknown modal reflection coefficients rm as N=2 X

pr ðx; tÞ ¼

rm pmr ðx; tÞ;

(18)

m¼N=2

where pmr(x, t) ¼ exp [i(km  x  xt)] represents a reflected plane wave of order m. Periodicity of the unit cell creates a diffraction grating on the interface27 and therefore the two components of the reflected wavevector km satisfy qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 with b ¼ 2p/a kmy ¼ kiy þ b2m and kmx ¼  ðw=ci Þ2  kmy 2 2 and integer m. The sign of the square root in kmx is chosen to have reflected propagating or evanescent waves travel away from the PC medium. In Ref. 27 the summation limit N is not explicitly indicated; we choose N reflected waves such that we will have square matrices for wave amplitudes and have a unique solution in the entire field of interest. To determine the modal amplitudes of the N Bloch waves and N reflected wave orders we write the pressure continuity equations, valid at all times on the interface, as N=2 X

pi ð0; y; tÞ þ

rm pmr ð0; y; tÞ ¼

m¼N=2

N X

tn pi ð0; y; tÞ:

n¼1

(19) An equation enforcing continuity of normal particle velocity in the ^e x direction is written as tix ð0; y; tÞ þ

N=2 X m¼N=2

rm vmrx ð0; y; tÞ ¼

N X

tn tnx ð0; y; tÞ;

n¼1

(20) where from the linearized Euler’s relationship the velocities of the incident and reflected plane waves are tix ¼

FIG. 4. (Color online) The dispersion diagram for the PC at normal incidence. Horizontal lines at 5, 10 and 18 krad/s. Squares indicate selected Bloch waves. Here an “” indicates propagating Bloch waves discarded as they radiate energy out of the PC half-space (cg  e^x < 0). The first two bandgaps are denoted by the frequency ranges outlined in shaded boxes. J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

1 p Z0 i

tmrx ¼

1 pmr : Z0

(21) (22)

The impedance of the mth plane wave order is given by Zm ¼ qix/kmx. The particle velocity of the nth Bloch wave is Kulpe et al.: Bloch-wave expansion for phononic crystals

1813

tnx ¼ ~t nx exp½iðkn  x  xtÞ   1 @~ pn ~t nx ¼ þ iknx p~n : ixq1 @x

(23) (24)

Note that the Bloch wave velocity in Eq. (23) also satisfies the Bloch theorem. Solution of the expansion equations for the amplitude coefficients results from multiplication of both equations Eqs. (19) and (20) by exp(ikmyy) and integration along the boundary of the unit cell to exploit orthogonality of complex exponentials. Following the integration, Eq. (19) reads (in matrix form) i þ r ¼ Mp t

(25)

(26)

where i is a vector of the incident plane wave order amplitudes (im ¼ dm0), t ¼ [t1 t2 … tN]T, and r ¼ [rN/2 … r0 … rN/2]T, and dmn is the Kronecker Delta. For compactness, we can then define27 three matrices with components defined by ð 1 a2 p~ ð0; yÞexpð imb2 yÞdy (27) Mp;mn ¼ a2 0 n ð 1 a2 ~v nx ð0; yÞexpð imb2 yÞdy Mv;mn ¼ (28) a2 0 Zmn ¼ Zm dmn :

(29)

From a mathematical perspective, the m, n element in Mp (or Mv) represents the Fourier coefficient of the mth harmonic of the nth Bloch pressure (velocity) wave. Equations (25) and (26) can then be solved for the wave amplitude vectors via t ¼ 2ðMp þ ZMv Þ1 i

(30)

r ¼ ðMp  ZMv ÞðMp þ ZMv Þ1 i:

(31)

In summary, the solution of Eqs. (30) and (31) yields vectors t and r which, respectively, contain the Bloch wave and reflected order amplitudes and thereby allow the pressure and velocity fields to be obtained by the modal expansion [e.g., see Eqs. (17) and (18)]. Using the solution for the wave amplitudes we can consider generalized power reflection RW and transmission TW coefficients, describing the power reflected by or transmitted into the PC half-space. Note that these coefficients also depend on the frequency x and h. Following the solution for the expansion amplitudes, the time averaged intensity vector fields are defined (* denotes complex conjugation) as

1814

 1  It ðxÞ ¼ Re pt ðxÞvt ðxÞ 2

(32a)

 1  Ir ðxÞ ¼ Re pr ðxÞvr ðxÞ : 2

(32b)

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

RW ¼

Sr Sr þ St

(33a)

TW ¼

St Sr þ St

(33b)

with Sr ¼

ð a2

^  Ir ð0; yÞdy n

(34a)

^  It ð0; yÞdy: n

(34b)

0

St ¼

ð a2 0

and velocity continuity from Eq. (20) Z1 ði  rÞ ¼ Mt t;

The reflected and transmitted velocity fields follow from Euler’s equation: v ¼ (1/ixq)rp. The RW and TW coefficients are computed by

Note that by construction, RW þ TW ¼ 1, which is mandated by energy conservation in a lossless medium. If, for a particular frequency and incidence angle h, there are no propagating Bloch waves then energy must be perfectly reflected back to the incident medium and RW ¼ 1 and TW ¼ 0. IV. VERIFICATION AND RESULTS

In this section we aim to verify our Bloch wave expansion in several test cases using an independent FE model implemented in COMSOL. To the author’s knowledge there is no previous literature presenting a verification of reflected and transmitted pressure fields of a PC by an independent numerical model. Three verification studies will be performed in Sec. IV A with three different forcing frequencies and incidence angles. Modeling a PC half-space in a FE model requires a 2D geometry with a large array of air-filled circular inclusions. The half-space model requires non-reflecting boundary conditions on the edges of the workable simulation domain, akin to perfectly matched layers for wave propagation modeling (but which were developed for homogeneous media only).38 It is essential, to the numerical verification accuracy, to have the damping occur in a PC, and not a homogeneous medium as in the case for perfectly matched layers, so we employ here additional PC domains but incorporating damping n(x, y) that monotonically increases with distance from the workable domain. The following damping expression was found suitable for the simulated cases to approximate nonreflecting boundary conditions:  s H2 a2 2 s1 : (35) nðx; yÞ ¼ n0 ð 6x  H1 a1 Þ 6y7 2 Here Hj is the discrete number of unit cells in the jth direction encompassing the workable PC domain, sj is an integer value herein set to 3, and n0 is a constant value between 1000 and 10 000 depending on the selected frequency. The damping parameters are arbitrarily tuned for each frequency to minimize reflection. Depending on the location of the individual damping domains, the first or Kulpe et al.: Bloch-wave expansion for phononic crystals

second polynomial terms in Eq. (35) will be removed and the appropriate þ/ sign will be selected to ensure the damping function increases away from the domain. This damping function is of arbitrary nature and other functions, such as an exponential, could be considered as well; the important damping function feature is a slow increase away from the real PC domain (see Fig. 5). The FEM verification model, depicted in Fig. 5, can be used for either normal or oblique incidence. All numerical simulations are executed using COMSOL with a monochromatic plane wave, of unit amplitude, incident from the far left, specified using built-in radiation conditions. The number of unit cells to include in all numerical comparison models is tuned to minimize computational effort yet yield a sufficiently “infinite” domain; typically the number of cells ranges from 100 for normal incidence (H1 H2) and up to 4000 for oblique incidence (H1 H2). A. Verification studies

The three verification models will be executed at selected angular frequencies: x ¼ 10 000, 5000, and 39 000 rad/s; the first two models at normal incidence (i.e., h ¼ 0), and the last at h ¼ 30 . In all three cases the geometric and material parameters are carried over from Sec. III.

It was previously asserted that a necessary Bloch wave criterion for an accurate result is cg  ^e x > 0. As an example, we intentionally swap the propagating Bloch wave with k ¼ 5.037 rad/m, a wave with group velocity of opposite sign, and use the same evanescent waves. Here the computed coefficients are RW ¼ 5.934 and TW ¼ 4.934 which is physically incorrect as we must have RW, TW  1. The incorrect Bloch wave selection is further examined, by comparing the real part of the pressure field versus the correct real part of the expansion and numerical model, in Fig. 6(b). The observed mismatch reinforces our assertion that the expansion must contain propagating Bloch waves with group velocity pointing into the medium. 2. Verification study 2

A plane wave encounters the PC half-space at normal incidence at x ¼ 5000 rad/s, which is a frequency within the first bandgap (see Fig. 4). There are no propagating waves at this frequency and the expansion only consists of two decaying waves and eight evanescent waves (see Fig. 4). Figure 7 displays a comparison between the real and imaginary parts

1. Verification study 1

The dispersion and expansion details discussed previously in the paper are used to compute the reflected and transmitted pressure fields Eqs. (17) and (18) and power coefficients Eq. (33) in each of the three cases; only pressure fields will be compared against our numerical model. At the frequency x ¼ 10 000 rad/s and h ¼ 0 the numerical procedure described in Sec. III, led to the selection of one propagating Bloch (k ¼ 5.037 rad/m) wave and eight evanescent waves. Both real and imaginary parts of the pressure field are seen in Fig. 6(a) and here RW ¼ 0.168 and TW ¼ 0.832. Excellent agreement can be seen between the numerical model using an array of 180 unit cells. In this study, and in the next two studies, we report agreement with modal deafness as discussed in Ref. 27, namely, that Bloch modes that are antisymmetric with respect to the incident plane wave are not excited and thus ti ¼ 0 for these waves.

FIG. 5. (Color online) The geometry of a FEM model designed for normal and oblique incidence with a PC domain, a homogeneous incident domain, and five damping PC domains to replicate a semi-infinite half-space. For each damping domain the arrows indicate the direction of damping, given by Eq. (35). J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

FIG. 6. (Color online) A comparison of the (a) real and imaginary components of the pressure field from the Bloch expansion (line) and the numerical model (dots) field, along a horizontal line centered at y ¼ 0.2 m. Here x ¼ 10 000 rad/s and h ¼ 0 (normal incidence) and the interface is demarcated by a vertical dashed line. Plot (b) compares the real pressure field from (a) with the result of an incorrect Bloch wave expansion intentionally disregarding the group velocity selection criterion (dashed curve). Kulpe et al.: Bloch-wave expansion for phononic crystals

1815

FIG. 7. (Color online) A comparison of the real (top) and imaginary (bottom) components of the pressure field from the Bloch expansion (line) and the numerical model (dots) field, along a horizontal line centered at y ¼ 0.2 m. Here x ¼ 5000 rad/s and h ¼ 0 (normal incidence) and the interface is demarcated by a vertical dashed line.

of the pressure field versus FEM simulations. The numerical verification model presented in Fig. 5 consists of 15 unit cells within the workable domain. Again, excellent agreement is seen between the results from the Bloch expansion and the numerical model results. Here we compute a power reflection coefficient equal to unity, as expected. 3. Verification study 3

This study for oblique incidence raises the frequency to x ¼ 39 000 rad/s and h ¼ 30 . Calculations proceed as before and the real and imaginary parts of the pressure field are compared between the expansion and numerical FEM results in Fig. 8(a). Again excellent agreement is seen in the results utilizing two decaying Bloch waves and seven evanescent waves. At this frequency 3760 unit cells for the infinite PC and over 4.5  106 nodal degrees of freedom were employed in the comparison. The comparison simulation time exceeded 2 h versus approximately 10 min elapsed for the Bloch wave expansion. The numerical results indicate perfect acoustic energy reflection for this study. To further explain the acoustic reflection, the plot in Fig. 8(b) shows the intensity field Ix(x, 8.2). Clearly, after 10 unit cells (x ¼ 4 m) the decaying intensity decreases to a very small value. The decaying intensity field, along with Eq. (31), result in perfect acoustic reflection (RW ¼ 1) despite the frequency not residing in a bandgap, a feature discussed in Sec. IV B. B. PC critical angle

In the last verification study it was shown that perfect reflection was recorded for a frequency not in a bandgap. Analogous to the zero acoustic wave energy transmission into a homogeneous fluid we refer to this incidence angle as a critical angle. At a critical angle hc and a fixed frequency, there will be no propagating Bloch waves. The Bloch 1816

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

FIG. 8. (Color online) A comparison of the (a) real and imaginary components of the pressure field from the Bloch expansion (line) and the numerical model (dots) field, along a horizontal line centered at y ¼ 8.2 m. Here x ¼ 39 000 rad/s and h ¼ 30 and the interface is demarcated by a vertical dashed line. Plot (b) is the time averaged intensity (note the difference in x axis limits) in the e^x direction in the PC, plotted along a horizontal line through the PC.

expansion consists solely of decaying and evanescent Bloch waves and the pressure field in the PC will always be attenuated. Accordingly, all acoustic energy will be perfectly reflected from the PC. Romero-Garcia et al.39 provided theoretical and experimental evidence of an attenuated field within the PC. They refer to this dispersion behavior as angular bandgaps. This reflection effect is an intrinsic property of the PC, and can be studied by simple geometrical analysis of the dispersion behavior in k-space. To illustrate the result we display an isofrequency contour at x(kx, ky) ¼ 39 000 rad/s in Fig. 9(a). This dispersion data is calculated by varying h over the interval 6arcsin(cip/xa2). Horizontal regions in k-space, illustrated in Fig. 9(a), correspond to incidence angles that yield complex wavenumbers. Let the set of critical wavenumbers, kc, contain the wavenumbers ky in the first Brillouin zone that correspond to critical angles [shaded regions in Fig. 9(a)]. Periodicity of the dispersion relationship in k-space amounts to several critical wavenumber regions, all shifted by integer multiples of b2. Critical angles can be calculated in the range (90 , 90 ] from Kulpe et al.: Bloch-wave expansion for phononic crystals

FIG. 10. (Color online) The power reflection coefficient computed versus frequency and incidence wave angle. Specific regions of high reflection are annotated on the plot, as bandgaps or critical angle effects. Recall 0 is normal incidence.

C. Parametric study FIG. 9. (Color online) (a) An isofrequency contour in k-space at 39 000 rad/s with real wavenumbers as dots and shaded regions, where no real wavenumbers exist, indicating wavenumbers corresponding to critical angles. (b) A graphical display of the critical angles where shaded regions indicate ranges of critical angles. (c) Validation of the critical angle observation versus direct numerical calculation (dots). Solid lines are from wedge regions from (b).



ci hc ¼ arcsin ðkc þ mb2 Þ x

 (36)

for any integer m such that the arcsin argument is 1. The critical angle data is displayed graphically in Fig. 9(b) such that any acoustic plane wave with incident wavevector ki within a shaded region will be perfectly reflected by the PC half-space. The assertion that there exist multiple incidence angles where perfect acoustic reflection occurs, as defined by Eq. (33), can be verified by direct numerical evaluation. In Fig. 9(c) a sweep of angles from 0 to 60 is performed and the power reflection coefficient is computed from the Bloch wave expansion. Symmetry of the PC implies angular symmetry of this data about h ¼ 0 . At critical angles, derived solely from the geometrical analysis of Fig. 9(a), the power reflection agrees favorably from the separate direct calculations. Note, the presented geometrical analysis only indicates angles of perfect reflection. The reflection at intermediate angles can only be determined by direct expansion computations and is not applicable to the analysis in this section. Subsequently, the reflection values in excess of 55 are not of present interest. The results in this figure corroborate our assertion of identification of PC critical angles via inspection of the isofrequency contour. We note that the identified critical angles differ from the angles defined by Bragg’s diffraction law of a2 cos hc ¼ mpci/x for some integer m. J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

In this section we study the relation between the power reflection and parameters such as frequency and incidence angle. This parametric study will employ the presented, general, expansion technique and serve to detail the acoustic reflection for a variety of incident angles and frequencies. The harmonic frequency is varied from 500 to 20 000 rad/s and the incidence angle is swept from 0 to 75 . The computed reflection, seen in Fig. 10, displays large regions, particularly at low frequencies, of perfect reflection. The regions of perfect reflection, attributed to the first or second bandgap or critical angles, are annotated on the plot in Fig. 10. As the frequency increases, the number of wavenumber dispersion surfaces in k-space increases, providing diminishing bandgaps and leading to the low reflection at frequencies exceeding 16 000 rad/s. At these higher frequencies the regions of high reflection are attributed to the critical angle effect, discussed earlier, and not complete bandgaps as in the low frequencies below 7000 rad/s.

V. CONCLUSION

In this paper the viability of a Bloch wave expansion method was demonstrated, with proper selection of wave modes based on a group velocity criterion, for calculating the reflected and transmitted wave fields of a semi-infinite PC half-space. The method is formulated for seamless integration with FEM software, resulting in efficient and general calculation of Bloch waves and their group velocities. The results herein show agreement of the expanded pressure fields with direct numerical simulations obtained with commercial FEM software for a variety of frequencies and incidence angles; the expansion procedure was significantly faster than the direct simulation in the oblique case. We have also provided a geometrical analysis of the critical angles of Kulpe et al.: Bloch-wave expansion for phononic crystals

1817

a PC, yielding perfect reflection, which also agrees well with direct numerical calculations. This simple geometrical analysis should be useful for future studies needing to estimate acoustic reflection from any 2D PC. Finally, we presented a detailed parametric study and observed reflection as a strong function of incident angle and frequency. In conclusion, the presented method and its results should have application in modeling wave reflection from general 2D lossless PCs subjected to plane waves with arbitrary incidence angles and frequencies.

ACKNOWLEDGMENTS

This work is supported by ONR Grant No. N000141110259. The authors are also very grateful for fruitful discussions with Dr. Vincent Laude regarding clarification of the Bloch wave expansion. FIG. 11. (Color online) A test for Bloch wave completeness by considering an expansion of Bloch waves to replicate an arbitrary function. The difference between the expansion and g(y) is plotted versus y. The inset shows the arbitrary function g with the FEM nodes shown as dots.

APPENDIX

It was stated in Refs. 29 and 30 that the Bloch waves, at a particular frequency and plane wave incidence angle, form a complete function basis. The work later27 implements a wave expansion, as is done in this paper, without thoroughly examining the completeness. We propose a simple numerical experiment to illustrate that the set of Bloch waves is a suitable complete basis for an expansion. To establish completeness we concoct an arbitrary function g(y) which is valued on the FEM nodes on the PC interface at x ¼ 0 and take an expansion of Bloch waves to replicate this function g. We take x ¼ 10 000 and h ¼ 0 and after processing propagating waves for group velocity we are left with 16 Bloch waves (following the five step procedure outlined in Sec. III). Note several more evanescent waves are included, versus the expansion for the same parameters in Sec. IV A 1, as we are not concerned with the rate of decay for evanescent waves for the sake of this section. The expansion is similar to Eq. (17), 16 X

an pn ð0; y; tÞ ¼ gð yÞeiwt ;

(A1)

n¼1

where now the unknown coefficients an are to be determined. Taking an inner product of both sides of Eq. (A1) with the mth Bloch wave gives a set of matrix equations for the wave amplitude with components given by Fmn an ¼ Cm ð a2 pm ð0; y; tÞpn ð0; y; tÞdy Fmn ¼

(A2a) (A2b)

0

Cm ¼

ð a2 0

gðyÞeiwt pm ð0; y; tÞdy:

(A2c)

The results of the expansion are shown in Fig. 11. One will note the arbitrary function is exactly represented, between the nodes, by the Bloch wave expansion. This demonstration is not meant to be a rigorous 1818

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

completeness proof for all model parameters and frequencies, but rather a numerical experiment illustrating, in fact, that the set of Bloch waves used for the expansion described in this paper is complete. 1

B. Gralak, S. Enoch, and G. Tayeb, “Anomalous refractive properties of photonic crystals,” J. Opt. Soc. Am. A 17, 1012–1020 (2000). 2 A. Sukhovich, L. Jing, and J. H. Page, “Negative refraction and focusing of ultrasound in two-dimensional phononic crystals,” Phys. Rev. B 77, 014301 (2008). 3 Y. Y. Chen and Z. Ye, “Theoretical analysis of acoustic stop bands in two-dimensional periodic scattering arrays,” Phys. Rev. E 64, 036616 (2001). 4 M. S. Kushwaha and P. Halevi, “Giant acoustic stop bands in twodimensional periodic arrays of liquid cylinders,” Appl. Phys. Lett. 69, 31–33 (1996). 5 M. Devaud, T. Hocquet, and V. Leroy, “Sound propagation in a monodisperse bubble cloud: From the crystal to the glass,” Eur. Phys. J. E: Soft Matter Biol. Phys. 32, 13–23 (2010). 6 M. S. Kushwaha, B. Djafari-Rouhani, and L. Dobrzynski, “Sound isolation from cubic arrays of air bubbles in water,” Phys. Lett. A 248, 252–256 (1998). 7 C. E. Bradley, “Time harmonic acoustic Bloch wave propagation in periodic waveguides. Part I. Theory,” J. Acoust. Soc. Am. 96, 1844–1853 (1994). 8 Z. Liu, X. Zhang, Y. Mao, Y. Zhu, Z. Yang, C. T. Chan, and P. Sheng, “Locally resonant sonic materials,” Science 289, 1734–1736 (2000). 9 M. Sigalas, M. S. Kushwaha, E. N. Economou, M. Kafesaki, I. E. Psarobas, and W. Steurer, “Classical vibrational modes in phononic lattices: Theory and experiment,” Z. Kristallogr. 220, 765–809 (2005). 10 A. Krynkin, O. Umnova, S. Taherzadeh, and K. Attenborough, “Analytical approximations for low frequency band gaps in periodic arrays of elastic shells,” J. Acoust. Soc. Am. 133, 781–791 (2013). 11 J.-H. Sun and T.-T. Wu, “Propagation of acoustic waves in phononiccrystal plates and waveguides using a finite-difference time-domain method,” Phys. Rev. B 76, 104304 (2007). 12 S. Gonella and M. Ruzzene, “Homogenization and equivalent in-plane properties of two-dimensional periodic lattices,” Int. J. Solids Struct. 45, 2897–2915 (2008). 13 A. A. Krokhin, P. Halevi, and J. Arriaga, “Long-wavelength limit (homogenization) for two-dimensional photonic crystals,” Phys. Rev. B 65, 115208 (2002). 14 M. I. Hussein, “Reduced Bloch mode expansion for periodic media band structure calculations,” Proc. R. Soc. London, Ser. A 465, 2825–2848 (2009). 15 F. Farzbod, “Analysis of Bloch formalism in undamped and damped periodic structures,” Thesis, Georgia Institute of Technology (2010). Kulpe et al.: Bloch-wave expansion for phononic crystals

16

D. Duhamel, B. R. Mace, and M. J. Brennan, “Finite element analysis of the vibrations of waveguides and periodic structures,” J. Sound Vib. 294, 205–220 (2006). 17 P. Langlet, A. Hladky-Hennion, and J. Decarpigny, “Analysis of the propagation of plane acoustic waves in passive periodic materials using the finite element method,” J. Acoust. Soc. Am. 98, 2792–2800 (1995). 18 J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Modeling the Flow of Light, 2nd ed. (Princeton University Press, Princeton, NJ, 2008), pp. 27–43, 221–227. 19 Z. Ruan, M. Qiu, S. Xiao, S. He, and L. Thylen, “Coupling between plane waves and Bloch waves in photonic crystals with negative refraction,” Phys. Rev. B 71, 045111 (2005). 20 C. Qiu, Z. Liu, J. Mei, and M. Ke, “The layer multiple-scattering method for calculating transmission coefficients of 2d phononic crystals,” Solid State Commun. 134, 765–770 (2005). 21 Z. Hou, X. Fu, and Y. Liu, “Calculational method to study the transmission properties of phononic crystals,” Phys. Rev. B 70, 014304 (2004). 22 V. Leroy, A. Bretagne, M. Fink, H. Willaime, P. Tabeling, and A. Tourin, “Design and characterization of bubble phononic crystals,” Appl. Phys. Lett. 95, 171904 (2009). 23 A. Bretagne, A. Tourin, and V. Leroy, “Enhanced and reduced transmission of acoustic waves with bubble meta-screens,” Appl. Phys. Lett. 99, 221906 (2011). 24 R. P. Moiseyenko, N. F. Declercq, and V. Laude, “Numerical investigation of diffraction of acoustic waves by phononic crystals,” AIP Conf. Proc. 1433, 319–322 (2012). 25 R. P. Moiseyenko, S. Herbison, N. F. Declercq, and V. Laude, “Phononic crystal diffraction gratings,” J. Appl. Phys. 111, 034907 (2012). 26 L. Sanchis, F. Cervera, J. Sanchez-Dehesa, J. V. Sanchez-Perez, C. Rubio, and R. Martinez-Sala, “Reflectance properties of two-dimensional sonic band-gap crystals,” J. Acoust. Soc. Am. 109, 2598–2605 (2001). 27 V. Laude, R. P. Moiseyenko, S. Benchabane, and N. F. Declercq, “Bloch wave deafness and modal conversion at a phononic crystal boundary,” AIP Adv. 1, 041402 (2011).

J. Acoust. Soc. Am., Vol. 135, No. 4, April 2014

28

V. Laude, B. Aoubiza, Y. Achaoui, S. Benchabane, and A. Khelif, “Evanescent Bloch waves in phononic crystals,” in SPIE OPTO: Integrated Optoelectronic Devices, San Jose, CA (International Society for Optics and Photonics, 2009). 29 Y.-C. Hsue, A. J. Freeman, and B.-Y. Gu, “Extended plane-wave expansion method in three-dimensional anisotropic photonic crystals,” Phys. Rev. B 72, 195118 (2005). 30 V. Laude, Y. Achaoui, S. Benchabane, and A. Khelif, “Evanescent Bloch waves and the complex band structure of phononic crystals,” Phys. Rev. B 80, 092301 (2009). 31 J. N. Reddy, An Introduction to the Finite Element Method, 3rd ed. (McGraw Hill, New York, 2006), pp. 410–523. 32 A. S. Phani, J. Woodhouse, and N. A. Fleck, “Wave propagation in two-dimensional periodic lattices,” J. Acoust. Soc. Am. 119, 1995–2005 (2006). 33 F. Farzbod and M. J. Leamy, “The treatment of forces in Bloch analysis,” J. Sound Vib. 325, 545–551 (2009). 34 N. van der Aa, “Perturbation theory for eigenvalue problems” (2005), http://www.win.tue.nl/casa/meetings/seminar/previous/_abstract051019_ files/Presentation.pdf. (Last accessed January 10, 2014). 35 J. A. Kulpe, M. J. Leamy, and K. G. Sabra, “Modeling the acoustic scattering from large fish schools using the Bloch-Floquet theorem,” Proc. Meet. Acoust. 19, 005026 (2013). 36 A. D. Pierce, Acoustics: An Introduction to its Physical Principles and Applications, 3rd ed. (Acoustical Society of America, Melville, NY, 1994), pp. 100–125. 37 P. Yeh, “Electromagnetic propagation in birefringent layered media,” J. Opt. Soc. Am. 69, 742–756 (1979). 38 O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 6th ed. (Elsevier, Burlington, MA, 2005), 637 pp. 39 V. Romero-Garcia, R. Pico, A. Cebrecos, K. Staliunas, and V. J. SanchezMorcillo, “Angular band gaps in sonic crystals: Evanescent waves and spatial complex dispersion relation,” J. Vib. Acoust. 135, 041012 (2013).

Kulpe et al.: Bloch-wave expansion for phononic crystals

1819

Copyright of Journal of the Acoustical Society of America is the property of American Institute of Physics and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.