- Email: [email protected]
Generalized eigenproblem of hybrid matrix for Floquet wave propagation in one-dimensional phononic crystals with solids and fluids
Ultrasonics 50 (2010) 91–98
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras
Generalized eigenproblem of hybrid matrix for Floquet wave propagation in one-dimensional phononic crystals with solids and fluids Eng Leong Tan * Nanyang Technological University, School of Electrical & Electronic Engineering, Block S2, Nanyang Ave, 639798, Singapore
a r t i c l e
i n f o
Article history: Received 8 January 2009 Received in revised form 14 September 2009 Accepted 15 September 2009 Available online 20 September 2009 PACS: 43.20.Bi 43.20.Fn 43.20.Gp
a b s t r a c t A method based on the solution to a generalized eigenproblem of hybrid matrix is presented for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method overcomes the numerical instability in the standard eigenproblem of transfer matrix, thus enabling Floquet waves to be determined reliably. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases are formulated. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids are discussed. The frequency–thickness range of phononic bandgap is determined conveniently based on the Floquet wavenumbers. Ó 2009 Elsevier B.V. All rights reserved.
Keywords: Phononic crystal Floquet wave Generalized eigenproblem Omnidirectional reflection
1. Introduction In recent years, there has been considerable interest in the study of phononic crystals structures, which are periodic composite materials with acoustic analogues of photonic crystals [1–7]. One of the celebrated techniques for analysis of such media is based on Floquet wave theory and transfer matrix method [8–13]. Although the transfer matrix method is applicable in principle, its direct implementation has been found to suffer from numerical instability when the layer is lossy, frequency is high and/or thickness is large. To circumvent the problem, many techniques have been proposed. One approach called the impedance (stiffness) matrix method has been introduced in [14–17]. Although the method maintains the numerical stability when the layer thickness grows to infinity, it is inaccurate when the thickness reduces toward zero [17]. Moreover, the Floquet waves for periodic media are really the eigenvectors of transfer matrix and not those of impedance matrix. In this paper, we present a method based on the solution to a generalized eigenproblem of hybrid matrix for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method utilizes the hybrid matrix of each layer in a recursive algorithm to deduce the stack hybrid matrix
* Tel.: +65 67906190; fax: +65 67933318. E-mail address: [email protected] 0041-624X/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2009.09.007
[18,19]. Similar to the impedance matrix method, the hybrid matrix method is able to eliminate the numerical instability of transfer matrix method. Moreover, contrary to the impedance matrix, the hybrid matrix remains to be well-conditioned and accurate even for zero or small thicknesses. Section 2 describes the standard eigenproblem of transfer matrix method. The stability and accuracy of both transfer matrix and impedance matrix methods are also discussed. Section 3 presents the generalized eigenproblem of hybrid matrix method. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases are formulated. While the standard eigenproblem of transfer matrix may have difficulties to determine Floquet waves, the generalized eigenproblem of hybrid matrix method can alleviate the problem. In Section 4, some numerical implementation examples are described. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids are discussed. The frequency–thickness range of phononic bandgap is determined conveniently based on the Floquet wavenumbers. 2. Standard eigenproblem of transfer matrix 2.1. Transfer matrix Fig. 1 shows the geometry of one-dimensional phononic crystals. There are in general n (could be 1) unit cells and the unit cell
92
E.L. Tan / Ultrasonics 50 (2010) 91–98 ðjÞ
ðjÞ
ing to Imðkz f Þ > 0 and Imðkz f Þ < 0, which remain bounded in the f and Pf upper and lower regions, respectively. Accordingly, wf , w can be decomposed into partitions as
" wf ¼
v>f
v
s>f
s
#
" f ¼ w
;
# >f w ;
" Pf ¼
P>f
0
0
P
# ;
ð6Þ
where the superscripts ‘‘>” and ‘‘<” stand for ‘‘upward-bounded” and ‘‘downward-bounded” decompositions, respectively. With the eigensolutions available for each layer of the unit cell, the field vector at the top of the layer can be related to that at the bottom by the layer transfer matrix:
"
ðZ >f Þ v
#
" ¼ Tf
sðZ >f Þ
ðZ
#
sðZ
ð7Þ
;
Tf ¼ wf Pf ðhf Þw1 f :
ð8Þ
By properly cascading the pertaining layer transfer matrices, the field vector at the top of the unit cell can be inferred from that at the bottom directly via the cell transfer matrix
Tc ¼ TN T2 T1 :
Fig. 1. Geometry of one-dimensional phononic crystals.
may comprise N homogeneous layers stratified in ^z direction. The upper and lower bounding interfaces of each layer f are denoted < by Z > f and Z f , respectively. For generality, the medium of each layer can be anisotropic solid or fluid material. Assuming a plane harmonic wave with exp(ixt) time dependence, each field vector satisfies a first-order differential system for solid as [20]
v d v ¼A : dz s s
ð1Þ
and s are the velocity vector and the normal stress vector, Here, v respectively. The system matrix A takes the (6 6) form
" A ¼ ix
sC1 33 C31
C1 33
qI s2 ðC11 C13 C1 sC13 C1 33 C31 Þ 33
# ;
ð2Þ
s is the transverse slowness along x direction, q the mass density, I the 3 3 identity matrix and C’s can be constructed from the stiffness constants using the abbreviated subscripts. For nonviscous and s reduce to their z components only and the correspondfluid, v ing system matrix A becomes 2 2 as
z d v z dz s
z v ¼A ; sz " # 0 q1 c12 s2
A ¼ ix
q
0
;
ð9Þ
2.2. Standard eigenproblem Using the cell transfer matrix, the wave propagation in a periodic stack of unit cells that constitute the one-dimensional phononic crystals can be analyzed conveniently. In particular, Floquet (F) waves are determined from the solutions to the standard eigenproblem
ðjÞ ¼ pðjÞ w ðjÞ : Tc w F F F
ð10Þ
In accordance with the boundedness association above, the eigensoðjÞ lutions are decomposed into upward-bounded ðjpF j < 1Þ and ðjÞ ðjÞ downward-bounded ðjpF j > 1Þ waves. pF can be assembled into < P> F and PF forming the matrix PF , while their associated eigenvecðjÞ > < < tors wF can also be assembled into v> F ; sF and vF ; sF forming the Floquet eigenwave matrix wF :
wF ¼
v>F s>F
"
; s
v
PF ¼
P>F
0
0
P
# :
The cell transfer matrix can be written as
ð3Þ
Tc ¼ wF PF ðhF Þw1 F ;
ð4Þ
where hF is the total thickness of unit cell. Furthermore, the Floquet wavenumber kzF can be deduced from (10) as (for integer m) ðjÞ
ð12Þ
ðjÞ
ðjÞ
c is the speed of sound in the fluid. (If the fluid is viscous, the corresponding system matrix is again 6 6 matrix similar to solid case.) Eq. (1) (or (3)) admits solutions in terms of the superposition of eigenwaves for each layer f:
kzF ¼ ½2mp þ argðpF Þ i ln jpF j=hF :
p6F trðTc Þp5F þ þ detðTc Þ ¼ 0:
v f ðzÞ sf ðzÞ
f ðzÞ: ¼ wf Pf ðzÞcf ¼ wf w
ð11Þ
ð5Þ
wf is a 6 6 (for solid, or 2 2 for fluid) eigenwave matrix comprising the eigenvectors; Pf(z) is a diagonal matrix with ðjÞ ðjÞ ðjÞ pf ðzÞ ¼ expðikz f zÞ, kz f being the jth wavenumber; cf is a 6 1 (or f ðzÞ is a 2 1) vector containing the unknown coefficients, and w 6 1 (or 2 1) vector which lumps the exponentials and coefficients together. Henceforth, the eigensolutions are decomposed into upward-bounded and downward-bounded waves correspond-
ð13Þ
Eq. (10) gives rise to a sixth-degree characteristic polynomial with its coefficients expressible in terms of the cell transfer matrix elements:
ð14Þ
When the unit cell is lossy and thick, the calculation of cell transfer matrix leads to numerical difficulties with the result being erroneous or unstable as
Tc jhF
!1
¼
0 v
s
v>F
0
s>F
0
1 :
ð15Þ
This renders the computation of the eigenvalues and eigenvectors inaccurate or wrong based on the incorrect Tc .
93
E.L. Tan / Ultrasonics 50 (2010) 91–98
There has been attempt to overcome the numerical problem of transfer matrix using impedance matrix
partitioned for several scenarios (the size of partition is indicated in the bracket subscript):
"
– layer f is solid, layer l is solid:
#
"
#
ðZ <1 Þ sðZ <1 Þ v ¼ Zc ; > ðZ >N Þ sðZ N Þ v "
Zc ¼
s
s>F s>F P>F
ð16Þ
#"
s
v>F
v
v>F P>F
v
½H
detðpF Zc12 p1 F Zc21 þ Zc11 Zc22 Þ ¼ 0;
ðZc12 Þ1 Zc11
ðZc12 Þ1
1
Zc22 ðZc12 Þ1
Zc21 Zc22 ðZc12 Þ Zc11
66 ¼
"
Zc jhF ¼1 ¼
0
0
s
s>F s>F
s
v>F v>F
v
" ½Hðf ;lÞ 44 ¼
ð20Þ
:
1
ð21Þ
;
We shall resort to an alternative method based on hybrid matrix, which is defined as [18]
ðZf Þ
¼ Hf
s
ðZf Þ v
s
"
# Hf ¼
;
Hf11 Hf21
Hf12 Hf22
# ð22Þ
:
The matrix Hf is called the layer hybrid matrix since it is a mixture of admittance (compliance), impedance (stiffness) and transfer elements. In each layer, one can determine Hf from the eigensolutions:
" Hf ¼
v>f
v
s>f P>f ðhf Þ
s
#"
s>f
s
v>f P>f ðhf Þ
v
#1 :
ð23Þ
ðZl Þ
s
#
" ¼ Hðf ;lÞ
sðZl Þ v
½H12 11
ðf ;lÞ
ðf ;lÞ
½H22 11
ðf ;lÞ
# ð25bÞ
;
ðf ;lÞ
½H12 31
ðf ;lÞ
½H22 11
½H11 33
ðf ;lÞ ðf ;lÞ
# ð25cÞ
;
ðf ;lÞ
½H12 13
ðf ;lÞ
½H22 33
½H11 11 ½H21 31
ðf ;lÞ ðf ;lÞ
ð25dÞ
:
ðf ;NÞ
ðf þ1;NÞ ðf þ1;NÞ 1 f ½I Hf22 H11 H21 ; ðf þ1;NÞ f 1 ðf þ1;NÞ f H12 ½I H11 H22 H12 ; ðf þ1;NÞ ðf þ1;NÞ 1 f H21 ½I Hf22 H11 H21 ; ðf þ1;NÞ ðf þ1;NÞ f ðf þ1;NÞ f 1 ðf þ1;NÞ H22 þ H21 H22 ½I H11 H22 H12 :
H11
¼ Hf11 þ Hf12 H11
ð26aÞ
ðf ;NÞ H12 ðf ;NÞ H21 ðf ;NÞ H22
¼
ð26bÞ
¼ ¼
ð26cÞ ð26dÞ
The algorithm proceeds until f = 1 and the cell hybrid matrix is obtained as
Hc ¼ Hð1;NÞ :
ð27Þ
For the two adjacent layers being of different phases, say layer f + 1 is solid and layer f is fluid, the interface conditions read
z ðZf Þ ¼ v v
ð28aÞ
z ðZ >f Þ ¼ z ðZ
ð28bÞ
s s
s
ð28cÞ
Let us suppose that layer N is solid and partition the ij-submatrix ðf þ1;NÞ as (i and j stand for 1 and/or 2) of stack hybrid matrix Hij
h i 33 ¼ ½Hðfij Itþ1;NÞ 32 ½Hðfij Izþ1;NÞ 31 2 3 2 ðf þ1;NÞ ðf þ1;NÞ ½H ½Hij tI 23 22 4 5 ¼ 4 ij tt ¼ ðf þ1;NÞ ðf þ1;NÞ ½Hij zI 13 ½Hij zt 12
ð24Þ
where H(f,l) is the total stack hybrid matrix from layer f to layer l. In and s are vectors of the above definitions one should recall that v 3 1 for solid and 1 1 for fluid. The stack hybrid matrix can be
#
– layers f+1 and f are of solid–solid or fluid–fluid:
ðf þ1;NÞ
# ;
ð25aÞ
;
To deduce the stack hybrid matrix for unit cell, we consider the recursion relation starting from the top layer N. Let us first assume that the two adjacent layers are of the same phase, i.e. solid–solid or fluid–fluid. Using the definitions of layer and stack hybrid matrices in (22) and (24) while noting the continuity conditions, one can obtain the recursion relation starting from the top layer N downward to the lower layer f as
½Hij
For solving multilayered problem, we also define the stack hybrid matrix as
"
ðf ;lÞ
½H11 11
½H21 13
" ½Hðf ;lÞ 44 ¼
3.1. Hybrid matrix
"
ðf ;lÞ
– layer f is fluid, layer l is solid:
3. Generalized eigenproblem of hybrid matrix
#
½H22 33
#
ð19Þ
:
<1 cf. P> ¼ I as hF = 0 in (17). Therefore, we see from above that F ¼ PF the eigenproblems for both transfer and impedance matrix methods will have numerical difficulties at large and/or small thicknesses.
"
ðf ;lÞ
½H21 11
#
One can see that both Zc21 and Zc12 are approaching zero, thus (18) is no longer a function of pF, and (19) is no longer invertible. On the other hand, when the cell thickness becomes small, the impedance matrix is ill-conditioned as
Zc jhF ¼0 ¼
ðf ;lÞ
– layer f is solid, layer l is fluid:
#
s>F ðv>F Þ1
½H12 33
½H21 33
½Hðf ;lÞ 22 ¼
However, when the cell thickness becomes large,
"
ðf ;lÞ
½H11 33
– layer f is fluid, layer l is fluid:
ð18Þ
the impedance matrix still cannot be applied directly for the periodic media considered here. This is because the Floquet waves are really the eigenvectors of transfer matrix and not those of impedance matrix. One might then consider the relation between transfer matrix and impedance matrix as
Tc ¼
ðf ;lÞ
ð17Þ
:
Although the characteristic polynomial (14) can be written in terms of the partitions as
"
"
#1
ðf þ1;NÞ
21
ðf þ1;NÞ
11
½Hij tz
½Hij zz
3 5: ð29Þ
By applying the interface conditions, we can obtain the recursion relation for the stack hybrid matrix as – layer f+1 is solid, layer f is fluid:
94
E.L. Tan / Ultrasonics 50 (2010) 91–98
ðf ;NÞ
ðf þ1;NÞ ðf þ1;NÞ ½I Hf22 H11zz 1 Hf21 ; ðf þ1;NÞ f 1 ðf þ1;NÞ f H12 ½I H11zz H22 H12zI ; ðf þ1;NÞ ðf þ1;NÞ H21Iz ½I Hf22 H11zz 1 Hf21 ; ðf þ1;NÞ ðf þ1;NÞ ðf þ1;NÞ ðf þ1;NÞ H22 þ H21Iz Hf22 ½I H11zz Hf22 1 H12zI :
H11
¼ Hf11 þ Hf12 H11zz
ð30aÞ
ðf ;NÞ H12 ðf ;NÞ H21 ðf ;NÞ H22
¼
ð30bÞ
¼ ¼
ð30cÞ ð30dÞ
Likewise, if layer f + 1 is fluid and layer f is solid, the stack hybrid matrix recursion relation can be obtained as – layer f+1 is fluid, layer f is solid: ðf ;NÞ
H11
f
f
ðf þ1;NÞ
¼ h11 þ h12 H11
ðf ;NÞ H12
¼
H21
ðf ;NÞ
¼
ðf ;NÞ H22
¼
f
ðf þ1;NÞ 1
½I h22 H11
f
h21 ;
f h12 ½I
ðf þ1;NÞ f 1 ðf þ1;NÞ H11 h22 H12 ; f ðf þ1;NÞ ðf þ1;NÞ 1 f H21 ½I h22 H11 h21 ; ðf þ1;NÞ ðf þ1;NÞ f ðf þ1;NÞ f 1 ðf þ1;NÞ H22 þ H21 h22 ½I H11 h22 H12 ;
ð31aÞ ð31bÞ ð31cÞ ð31dÞ
where f
ð32aÞ
f
ð32bÞ
h11 ¼ Hf11 Hf12It ðHf22tt Þ1 Hf21tI ; h12 ¼ Hf12Iz Hf12It ðHf22tt Þ1 Hf22tz ; f
h21 ¼ Hf21zI Hf22zt ðHf22tt Þ1 Hf21tI ;
ð32cÞ
f h22
ð32dÞ
¼
Hf22zz
Hf22zt ðHf22tt Þ1 Hf22tz :
From the above, we see that both (30) and (31) have similar form as (26), except partitioning of stack hybrid submatrices in (30) and replacement of layer hybrid submatrices by h in (31).
With the aid of stable QZ algorithms [21], the eigensolutions to (33) can be computed reliably through a(j) and b(j) (rather than ðjÞ pF ¼ aðjÞ =bðjÞ directly), where a(j) and b(j) are the diagonal elements of the triangular matrices QHAZ and QHBZ obtained by QZ factorization. Note that it is possible that b(j) ? 0 for some j, which could ðjÞ then make jpF j ! 1. In practice, such infinite eigenvalue can be avoided if one adopts the upward-bounded and downwardðjÞ ðjÞ bounded associations corresponding to P> F ¼ diagða =b Þ for < 1 ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ¼ diagðb =a Þ for jb j < ja j, respectively. ja j < jb j and PF ðjÞ are assembled into v> ; s> and v< ; s< formThen the eigenvectors w F F F F F ing the Floquet eigenwave matrix wF as before. Furthermore, the Floquet wavenumber kzF can be deduced from (33) as (for integer m) ðjÞ
kzF ¼ ½2mp þ argðaðjÞ Þ argðbðjÞ Þ iðln jaðjÞ j ln jbðjÞ jÞ=hF :
ð35Þ
In Fig. 2, we show the dispersion relation of Floquet wavenumber kzF hF =2p (real part) versus fhF as calculated by (a) standard eigenproblem of transfer matrix and (b) generalized eigenproblem of hybrid matrix. The unit cell comprises aluminium and epoxy with thickness hAl ¼ hEpoxy ¼ hF =2, and the incident angle is 80° corresponding to longitudinal wave in nylon. The materials are characterized by the parameters shown in Table 1. To better demonstrate the numerical problem, the cell materials are also arbitrarily made lossy with complex viscoelastic constants, 00 C IJ ¼ C 0IJ iC IJ , where C 00IJ ¼ 0:01C 0IJ . From Fig. 2, it is observed that the first method may give corrupted results for fhF larger than about 10 103 m=s. This is because when fhF increases, the cell
3.2. Generalized eigenproblem Using the cell hybrid matrix, Floquet waves can be determined via a generalized eigenproblem:
ðjÞ bðjÞ ¼ HB w ðjÞ aðjÞ ; HA w F F I Hc11 Hc12 ; HB ¼ HA ¼ 0 Hc21 Hc22
0 : I
ð33Þ ð34Þ
Table 1 Material parameters.
Al Epoxy Nylon
Density (kg/m3)
Longitudinal velocity (m/s)
Shear velocity (m/s)
2700 1140 1100
6380 2770 2600
3140 1300 1100
Fig. 2. Relation of kzFhF/2p (real part) and fhF as calculated by (a) standard eigenproblem of transfer matrix and (b) generalized eigenproblem of hybrid matrix. The unit cell comprises (viscoelastic) aluminium and epoxy with hAl = hEpoxy. The incident angle is 80° corresponding to longitudinal wave in nylon.
E.L. Tan / Ultrasonics 50 (2010) 91–98
transfer matrix becomes numerically unstable and inaccuracies arise in computing kzF . Although we have considered some loss here, the numerical problem may still persist even without viscos-
95
ity. On the other hand, the second method using the generalized eigenproblem of hybrid matrix can alleviate the problem, giving uncorrupted results over the range considered. (For extreme fhF , ðjÞ one may anticipate that some jpF j ! 0; 1 or aðjÞ ; bðjÞ ! 0.) 3.3. Application for unit cell of solid and fluid
Fig. 3. Structures for unit cell composed of solid and fluid.
Without loss of generality, let us consider a unit cell of infinite periodic media that is composed of one solid layer and one fluid layer (solid–fluid cell), as shown in Fig. 3a. To treat the solid, one previous approach is to model only the longitudinal wave [2]. However, this approach is correct only when the wave is incident normally as to be shown below. Meanwhile, one can notice that the cell is equivalent to the structure that is made up by a solid layer placed between two fluid layers (fluid–solid–fluid cell), as shown in Fig. 3b. Alternatively, it is also equivalent to the structure that is made up by a fluid layer placed between two solid layers (solid–fluid–solid cell), as shown in Fig. 3c. For all the structures of Fig. 3, the generalized eigenproblem of hybrid matrix method may be applied by taking into account all acoustic waves (including shear waves in solids). Note that for Fig. 3b, since HA and HB are 2 2 matrices, the generalized eigenproblem of (33) has two eigenvalues and eigenvectors. On the other hand, for Fig. 3c, HA and HB are 6 6 matrices. In this case, because HA and HB are not full rank matrices, the generalized eigenproblem of (33) still has two eigenvalues and eigenvectors after one eliminates four fictitious solutions. As for Fig. 3a, HA and HB are 4 4 matrices and (34) needs to be modified slightly (easier alternative method is described below). In Fig. 4, we plot the relation of kzFhF/2p (real part) and fhF for (a) 0° (normal) and (b) 80° incidence. The structures are those of Fig. 3a (solid–fluid hso ¼ hfl ) considering all waves in solid and only longitudinal wave in solid, Fig. 3b (fluid–solid–fluid 2hfl1 ¼ hfl ¼ hso ) and Fig. 3c (solid–fluid–solid 2hso1 ¼ hso ¼ hfl ). The solid and fluid are aluminium and water, respectively. Here
Fig. 4. Relation of kzFhF/2p (real part) and fhF for (a) 0° (normal) and (b) 80° incidence. The markers are for the structures of Fig. 3a (solid–fluid hso ¼ hfl ): ‘’ – all waves in solid, ‘o’ – only longitudinal wave in solid, Fig. 3b (fluid–solid–fluid 2hfl1 ¼ hfl ¼ hso ): ‘+’ and Fig. 3c (solid–fluid–solid 2hso1 ¼ hso ¼ hfl ): ‘x’. The solid and fluid are aluminium and water, respectively. The incident angles correspond to longitudinal waves in aluminium.
96
E.L. Tan / Ultrasonics 50 (2010) 91–98
and henceforth, no viscosity is assumed for simplicity. The incident angles correspond to longitudinal waves in aluminium. From Fig. 4, one can see that the results for all equivalent structures of Fig. 4 match very well at normal and oblique incidence. The approach that models only the longitudinal wave in solid is correct only when the wave is incident normally. Further examination shows that in the situation illustrated by Fig. 3b or Fig. 3c, the thickness hfl1 or hso1 can be set to arbitrary portion of hfl or hso. When the thickness becomes small,
the hybrid matrix has the advantage that it is still well-conditioned as
Hf jhf ¼0 ¼
0
I
I
0
:
ð36Þ
Then one can even set the thickness to be zero, i.e. hfl1 = 0 or hso1 = 0. This will reduce Fig. 3b or Fig. 3c to the initial situation in Fig. 3a, which can now be analyzed easily by the generalized eigenproblem of hybrid matrix method [without having to modify (34)].
Fig. 5. Real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. The unit cell comprises aluminium and epoxy with hAl = hEpoxy. The incident angles correspond to longitudinal waves in aluminium.
Fig. 6. Power reflection coefficients versus incident angles at fhF ¼ 2000m=s for n ¼ 1; 2; 8 cells, cf. Fig. 5.
E.L. Tan / Ultrasonics 50 (2010) 91–98
97
Fig. 7. Real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. The unit cell comprises a water layer between the two solid layers of Fig. 5: hWater = 0.9hAl, hEpoxy = 0.1hAl.
4. Omnidirectional reflection In this section, we apply the generalized eigenproblem of hybrid matrix method to analyze the omnidirectional reflection in onedimensional phononic crystals with n cells (nN layers) as shown in Fig. 1. Let the acoustic wave be incident from layer 0, and the reflection (r0,1) and transmission (t0,nN+1) coefficient matrices be defined as
<0 ðZ >0 Þ ¼ r0;1 w >0 ðZ >0 Þ; w >nNþ1 ðZ0 ðZ >0 Þ: w
ð37aÞ ð37bÞ
The stack hybrid matrix can be easily obtained using Floquet waves as
" ð1;nNÞ
H
¼
v>F
v
s>F ðP>F Þn
s
#"
s>F
s
v>F ðP>F Þn
v
#1 :
ð38Þ
< < By applying the radiation condition w nNþ1 ðZ nNþ1 Þ ¼ 0, the reflection and transmission coefficient matrices can be solved explicitly in terms of stack hybrid matrix as
r0;1 ¼ ½HS s<0 v <0 1 ½v >0 HS s>0 ; ð1;nNÞ
ð1;nNÞ > v>nNþ1 1 H21 ½s0 þ s<0 r0;1 ; ð1;nNÞ ð1;nNÞ ð1;nNÞ ð1;nNÞ HS ¼ H11 þ H12 ½s>nNþ1 ðv >nNþ1 Þ1 H22 1 H21 :
t0;nNþ1 ¼ ½s>nNþ1 H22
sponding power reflection coefficients versus incident angles at fhF = 2000 m/s for n = 1, 2, 8 cells. From the figure one can see that the reflectance tends toward unity as the number of the cells increases thus achieving omnidirectional reflection. Next we investigate the effect of fluid layer on the phononic bandgap. The unit cell now comprises a water layer between the two solid layers of Fig. 5 (keeping the same cell thickness): hWater = 0.9hAl, hEpoxy = 0.1hAl. In Fig. 7 we plot both real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. For this structure, the phononic bandgaps occur around fhF (0.6 1.5) 103 m/s and fhF (1.9–2.9) 103 m/s. It is seen here that one may obtain higher frequency bandgap by introducing a fluid layer between two solid layers in unit cell [7].
ð39aÞ ð39bÞ ð39cÞ
In Fig. 5, we plot both real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles (0°–90° step 5°). The unit cell comprises aluminium and epoxy with hAl = hEpoxy. The incident angles correspond to longitudinal waves in aluminium. We find that it is especially convenient to determine (if any) the frequency–thickness range of omnidirectional reflection using such scan plot of the imaginary parts of Floquet wavenumbers. Recall that if the imaginary parts of all kzF are not zero, the acoustic waves will decay rapidly along ^z direction. From the figure, one can clearly see that there is a phononic bandgap around fhF (1.8 2.2) 103 m/s and another narrower one at fhF 1.15 103 m/s. To verify this, we plot in Fig. 6 the corre-
5. Conclusion In this paper, a method based on the solution to a generalized eigenproblem of hybrid matrix has been presented for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method overcomes the numerical instability in the standard eigenproblem of transfer matrix. Thus, it enables Floquet waves to be determined reliably as has been demonstrated numerically. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases have been formulated. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids have been discussed. The frequency–thickness range of phononic bandgap has been determined conveniently based on the Floquet wavenumbers.
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] M.R. Shen, W.W. Cao, Acoustic bandgap formation in a periodic structure with multilayer unit cells, J. Phys. D: Appl. Phys. 33 (2000) 1150–1154.
98
E.L. Tan / Ultrasonics 50 (2010) 91–98
[3] D. Bria, B. Djafari-Rouhani, A. Bousfia, E.H. El Boudouti, A. Nougaoui, Absolute acoustic band gap in coupled multilayer structures, Europhys. Lett. 55 (6) (2001) 841–846. [4] A. Bousfia, E.H. El Boudouti, B. Djafari-Rouhani, D. Bria, A. Nougaoui, V.R. Velasco, Omnidirectional phononic reflection and selective transmission in one-dimensional acoustic layered structures, Surf. Sci. 1175 (2001) 482– 485. [5] D. Bria, B. Djafari-Rouhani, Omnidirectional elastic band gap in finite lamellar structures, Phys. Rev. E 66 (2002) 056609. [6] B. Manzanares-Martinez, J. Sanchez-Dehesa, A. Hakansson, F. Cervera, F. Ramos-Mendieta, Experimental evidence of omnidirectional elastic bandgap in finite one-dimensional phononic systems, Appl. Phys. Lett. 85 (1) (2004) 154–156. [7] Y. El Hassouani, E.H. El Boudouti, B. Djafari-Rouhani, H. Aynaou, Sagittal acoustic waves in finite solid–fluid superlattices: band-gap structure, surface and confined modes, and omnidirectional reflection and selective transmission, Phys. Rev. B 78 (2008) 174306. [8] W.T. Thomson, Transmission of elastic waves through a stratified solid medium, J. Appl. Phys. 21 (1950) 89–93. [9] N.A. Haskell, The dispersion of surface waves on multilayered media, Bull. Seismol. Soc. Am. 43 (1953) 17–34. [10] K.E. Gilbert, A propagator matrix method for periodically stratified media, J. Acoust. Soc. Am. 73 (1) (1983) 137–142. [11] M. Schoenberg, Wave propagation in alternating solid and fluid layers, Wave Motion 6 (1984) 303–320.
[12] A.H. Nayfeh, The general problem of elastic wave propagation in multilayered anisotropic media, J. Acoust. Soc. Am. 89 (1991) 1521–1531. [13] A.M.B. Braga, G. Herrmann, Floquet waves in anisotropic periodically layered composites, J. Acoust. Soc. Am. 91 (1992) 1211–1227. [14] B. Hosten, M. Castaings, Surface impedance matrices to model the propagation in layered anisotropic media, Ultrasonics 41 (2003) 501–507. [15] B. Collet, Recursive surface impedance matrix methods for ultrasonic wave propagation in piezoelectric multilayers, Ultrasonics 42 (2004) 189–197. [16] S.I. Rokhlin, L. Wang, Stable recursive algorithm for elastic wave propagation in layered anisotropic media: stiffness matrix method, J. Acoust. Soc. Am. 112 (3I) (2002) 822–834. [17] E.L. Tan, Stiffness matrix method with improved efficiency for elastic wave propagation in layered anisotropic media, J. Acoust. Soc. Am. 118 (6) (2005) 3400–3403. [18] E.L. Tan, Hybrid compliance–stiffness matrix method for stable analysis of elastic wave propagation in multilayered anisotropic media, J. Acoust. Soc. Am. 119 (1) (2006) 45–53. [19] E.L. Tan, Matrix algorithms for modeling acoustic waves in piezoelectric multilayers, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 54 (10) (2007) 2016–2023. [20] E.L. Tan, A concise and efficient scattering matrix formalism for stable analysis of elastic wave propagation in multilayered anisotropic solids, Ultrasonics 41 (2003) 229–236. [21] E.L. Tan, Generalized eigenproblem for acoustic wave propagation in periodically layered anisotropic media, J. Comput. Acoust. 16 (1) (2008) 1–10.
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras
Generalized eigenproblem of hybrid matrix for Floquet wave propagation in one-dimensional phononic crystals with solids and fluids Eng Leong Tan * Nanyang Technological University, School of Electrical & Electronic Engineering, Block S2, Nanyang Ave, 639798, Singapore
a r t i c l e
i n f o
Article history: Received 8 January 2009 Received in revised form 14 September 2009 Accepted 15 September 2009 Available online 20 September 2009 PACS: 43.20.Bi 43.20.Fn 43.20.Gp
a b s t r a c t A method based on the solution to a generalized eigenproblem of hybrid matrix is presented for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method overcomes the numerical instability in the standard eigenproblem of transfer matrix, thus enabling Floquet waves to be determined reliably. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases are formulated. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids are discussed. The frequency–thickness range of phononic bandgap is determined conveniently based on the Floquet wavenumbers. Ó 2009 Elsevier B.V. All rights reserved.
Keywords: Phononic crystal Floquet wave Generalized eigenproblem Omnidirectional reflection
1. Introduction In recent years, there has been considerable interest in the study of phononic crystals structures, which are periodic composite materials with acoustic analogues of photonic crystals [1–7]. One of the celebrated techniques for analysis of such media is based on Floquet wave theory and transfer matrix method [8–13]. Although the transfer matrix method is applicable in principle, its direct implementation has been found to suffer from numerical instability when the layer is lossy, frequency is high and/or thickness is large. To circumvent the problem, many techniques have been proposed. One approach called the impedance (stiffness) matrix method has been introduced in [14–17]. Although the method maintains the numerical stability when the layer thickness grows to infinity, it is inaccurate when the thickness reduces toward zero [17]. Moreover, the Floquet waves for periodic media are really the eigenvectors of transfer matrix and not those of impedance matrix. In this paper, we present a method based on the solution to a generalized eigenproblem of hybrid matrix for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method utilizes the hybrid matrix of each layer in a recursive algorithm to deduce the stack hybrid matrix
* Tel.: +65 67906190; fax: +65 67933318. E-mail address: [email protected] 0041-624X/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2009.09.007
[18,19]. Similar to the impedance matrix method, the hybrid matrix method is able to eliminate the numerical instability of transfer matrix method. Moreover, contrary to the impedance matrix, the hybrid matrix remains to be well-conditioned and accurate even for zero or small thicknesses. Section 2 describes the standard eigenproblem of transfer matrix method. The stability and accuracy of both transfer matrix and impedance matrix methods are also discussed. Section 3 presents the generalized eigenproblem of hybrid matrix method. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases are formulated. While the standard eigenproblem of transfer matrix may have difficulties to determine Floquet waves, the generalized eigenproblem of hybrid matrix method can alleviate the problem. In Section 4, some numerical implementation examples are described. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids are discussed. The frequency–thickness range of phononic bandgap is determined conveniently based on the Floquet wavenumbers. 2. Standard eigenproblem of transfer matrix 2.1. Transfer matrix Fig. 1 shows the geometry of one-dimensional phononic crystals. There are in general n (could be 1) unit cells and the unit cell
92
E.L. Tan / Ultrasonics 50 (2010) 91–98 ðjÞ
ðjÞ
ing to Imðkz f Þ > 0 and Imðkz f Þ < 0, which remain bounded in the f and Pf upper and lower regions, respectively. Accordingly, wf , w can be decomposed into partitions as
" wf ¼
v>f
v
s>f
s
#
" f ¼ w
;
# >f w ;
" Pf ¼
P>f
0
0
P
# ;
ð6Þ
where the superscripts ‘‘>” and ‘‘<” stand for ‘‘upward-bounded” and ‘‘downward-bounded” decompositions, respectively. With the eigensolutions available for each layer of the unit cell, the field vector at the top of the layer can be related to that at the bottom by the layer transfer matrix:
"
ðZ >f Þ v
#
" ¼ Tf
sðZ >f Þ
ðZ
#
sðZ
ð7Þ
;
Tf ¼ wf Pf ðhf Þw1 f :
ð8Þ
By properly cascading the pertaining layer transfer matrices, the field vector at the top of the unit cell can be inferred from that at the bottom directly via the cell transfer matrix
Tc ¼ TN T2 T1 :
Fig. 1. Geometry of one-dimensional phononic crystals.
may comprise N homogeneous layers stratified in ^z direction. The upper and lower bounding interfaces of each layer f are denoted < by Z > f and Z f , respectively. For generality, the medium of each layer can be anisotropic solid or fluid material. Assuming a plane harmonic wave with exp(ixt) time dependence, each field vector satisfies a first-order differential system for solid as [20]
v d v ¼A : dz s s
ð1Þ
and s are the velocity vector and the normal stress vector, Here, v respectively. The system matrix A takes the (6 6) form
" A ¼ ix
sC1 33 C31
C1 33
qI s2 ðC11 C13 C1 sC13 C1 33 C31 Þ 33
# ;
ð2Þ
s is the transverse slowness along x direction, q the mass density, I the 3 3 identity matrix and C’s can be constructed from the stiffness constants using the abbreviated subscripts. For nonviscous and s reduce to their z components only and the correspondfluid, v ing system matrix A becomes 2 2 as
z d v z dz s
z v ¼A ; sz " # 0 q1 c12 s2
A ¼ ix
q
0
;
ð9Þ
2.2. Standard eigenproblem Using the cell transfer matrix, the wave propagation in a periodic stack of unit cells that constitute the one-dimensional phononic crystals can be analyzed conveniently. In particular, Floquet (F) waves are determined from the solutions to the standard eigenproblem
ðjÞ ¼ pðjÞ w ðjÞ : Tc w F F F
ð10Þ
In accordance with the boundedness association above, the eigensoðjÞ lutions are decomposed into upward-bounded ðjpF j < 1Þ and ðjÞ ðjÞ downward-bounded ðjpF j > 1Þ waves. pF can be assembled into < P> F and PF forming the matrix PF , while their associated eigenvecðjÞ > < < tors wF can also be assembled into v> F ; sF and vF ; sF forming the Floquet eigenwave matrix wF :
wF ¼
v>F s>F
"
; s
v
PF ¼
P>F
0
0
P
# :
The cell transfer matrix can be written as
ð3Þ
Tc ¼ wF PF ðhF Þw1 F ;
ð4Þ
where hF is the total thickness of unit cell. Furthermore, the Floquet wavenumber kzF can be deduced from (10) as (for integer m) ðjÞ
ð12Þ
ðjÞ
ðjÞ
c is the speed of sound in the fluid. (If the fluid is viscous, the corresponding system matrix is again 6 6 matrix similar to solid case.) Eq. (1) (or (3)) admits solutions in terms of the superposition of eigenwaves for each layer f:
kzF ¼ ½2mp þ argðpF Þ i ln jpF j=hF :
p6F trðTc Þp5F þ þ detðTc Þ ¼ 0:
v f ðzÞ sf ðzÞ
f ðzÞ: ¼ wf Pf ðzÞcf ¼ wf w
ð11Þ
ð5Þ
wf is a 6 6 (for solid, or 2 2 for fluid) eigenwave matrix comprising the eigenvectors; Pf(z) is a diagonal matrix with ðjÞ ðjÞ ðjÞ pf ðzÞ ¼ expðikz f zÞ, kz f being the jth wavenumber; cf is a 6 1 (or f ðzÞ is a 2 1) vector containing the unknown coefficients, and w 6 1 (or 2 1) vector which lumps the exponentials and coefficients together. Henceforth, the eigensolutions are decomposed into upward-bounded and downward-bounded waves correspond-
ð13Þ
Eq. (10) gives rise to a sixth-degree characteristic polynomial with its coefficients expressible in terms of the cell transfer matrix elements:
ð14Þ
When the unit cell is lossy and thick, the calculation of cell transfer matrix leads to numerical difficulties with the result being erroneous or unstable as
Tc jhF
!1
¼
0 v
s
v>F
0
s>F
0
1 :
ð15Þ
This renders the computation of the eigenvalues and eigenvectors inaccurate or wrong based on the incorrect Tc .
93
E.L. Tan / Ultrasonics 50 (2010) 91–98
There has been attempt to overcome the numerical problem of transfer matrix using impedance matrix
partitioned for several scenarios (the size of partition is indicated in the bracket subscript):
"
– layer f is solid, layer l is solid:
#
"
#
ðZ <1 Þ sðZ <1 Þ v ¼ Zc ; > ðZ >N Þ sðZ N Þ v "
Zc ¼
s
s>F s>F P>F
ð16Þ
#"
s
v>F
v
v>F P>F
v
½H
detðpF Zc12 p1 F Zc21 þ Zc11 Zc22 Þ ¼ 0;
ðZc12 Þ1 Zc11
ðZc12 Þ1
1
Zc22 ðZc12 Þ1
Zc21 Zc22 ðZc12 Þ Zc11
66 ¼
"
Zc jhF ¼1 ¼
0
0
s
s>F s>F
s
v>F v>F
v
" ½Hðf ;lÞ 44 ¼
ð20Þ
:
1
ð21Þ
;
We shall resort to an alternative method based on hybrid matrix, which is defined as [18]
ðZ
¼ Hf
s
ðZ
s
"
# Hf ¼
;
Hf11 Hf21
Hf12 Hf22
# ð22Þ
:
The matrix Hf is called the layer hybrid matrix since it is a mixture of admittance (compliance), impedance (stiffness) and transfer elements. In each layer, one can determine Hf from the eigensolutions:
" Hf ¼
v>f
v
s>f P>f ðhf Þ
s
#"
s>f
s
v>f P>f ðhf Þ
v
#1 :
ð23Þ
ðZ
s
#
" ¼ Hðf ;lÞ
sðZ
½H12 11
ðf ;lÞ
ðf ;lÞ
½H22 11
ðf ;lÞ
# ð25bÞ
;
ðf ;lÞ
½H12 31
ðf ;lÞ
½H22 11
½H11 33
ðf ;lÞ ðf ;lÞ
# ð25cÞ
;
ðf ;lÞ
½H12 13
ðf ;lÞ
½H22 33
½H11 11 ½H21 31
ðf ;lÞ ðf ;lÞ
ð25dÞ
:
ðf ;NÞ
ðf þ1;NÞ ðf þ1;NÞ 1 f ½I Hf22 H11 H21 ; ðf þ1;NÞ f 1 ðf þ1;NÞ f H12 ½I H11 H22 H12 ; ðf þ1;NÞ ðf þ1;NÞ 1 f H21 ½I Hf22 H11 H21 ; ðf þ1;NÞ ðf þ1;NÞ f ðf þ1;NÞ f 1 ðf þ1;NÞ H22 þ H21 H22 ½I H11 H22 H12 :
H11
¼ Hf11 þ Hf12 H11
ð26aÞ
ðf ;NÞ H12 ðf ;NÞ H21 ðf ;NÞ H22
¼
ð26bÞ
¼ ¼
ð26cÞ ð26dÞ
The algorithm proceeds until f = 1 and the cell hybrid matrix is obtained as
Hc ¼ Hð1;NÞ :
ð27Þ
For the two adjacent layers being of different phases, say layer f + 1 is solid and layer f is fluid, the interface conditions read
z ðZ
ð28aÞ
z ðZ >f Þ ¼ z ðZ
ð28bÞ
s s
s
ð28cÞ
Let us suppose that layer N is solid and partition the ij-submatrix ðf þ1;NÞ as (i and j stand for 1 and/or 2) of stack hybrid matrix Hij
h i 33 ¼ ½Hðfij Itþ1;NÞ 32 ½Hðfij Izþ1;NÞ 31 2 3 2 ðf þ1;NÞ ðf þ1;NÞ ½H ½Hij tI 23 22 4 5 ¼ 4 ij tt ¼ ðf þ1;NÞ ðf þ1;NÞ ½Hij zI 13 ½Hij zt 12
ð24Þ
where H(f,l) is the total stack hybrid matrix from layer f to layer l. In and s are vectors of the above definitions one should recall that v 3 1 for solid and 1 1 for fluid. The stack hybrid matrix can be
#
– layers f+1 and f are of solid–solid or fluid–fluid:
ðf þ1;NÞ
# ;
ð25aÞ
;
To deduce the stack hybrid matrix for unit cell, we consider the recursion relation starting from the top layer N. Let us first assume that the two adjacent layers are of the same phase, i.e. solid–solid or fluid–fluid. Using the definitions of layer and stack hybrid matrices in (22) and (24) while noting the continuity conditions, one can obtain the recursion relation starting from the top layer N downward to the lower layer f as
½Hij
For solving multilayered problem, we also define the stack hybrid matrix as
"
ðf ;lÞ
½H11 11
½H21 13
" ½Hðf ;lÞ 44 ¼
3.1. Hybrid matrix
"
ðf ;lÞ
– layer f is fluid, layer l is solid:
3. Generalized eigenproblem of hybrid matrix
#
½H22 33
#
ð19Þ
:
<1 cf. P> ¼ I as hF = 0 in (17). Therefore, we see from above that F ¼ PF the eigenproblems for both transfer and impedance matrix methods will have numerical difficulties at large and/or small thicknesses.
"
ðf ;lÞ
½H21 11
#
One can see that both Zc21 and Zc12 are approaching zero, thus (18) is no longer a function of pF, and (19) is no longer invertible. On the other hand, when the cell thickness becomes small, the impedance matrix is ill-conditioned as
Zc jhF ¼0 ¼
ðf ;lÞ
– layer f is solid, layer l is fluid:
#
s>F ðv>F Þ1
½H12 33
½H21 33
½Hðf ;lÞ 22 ¼
However, when the cell thickness becomes large,
"
ðf ;lÞ
½H11 33
– layer f is fluid, layer l is fluid:
ð18Þ
the impedance matrix still cannot be applied directly for the periodic media considered here. This is because the Floquet waves are really the eigenvectors of transfer matrix and not those of impedance matrix. One might then consider the relation between transfer matrix and impedance matrix as
Tc ¼
ðf ;lÞ
ð17Þ
:
Although the characteristic polynomial (14) can be written in terms of the partitions as
"
"
#1
ðf þ1;NÞ
21
ðf þ1;NÞ
11
½Hij tz
½Hij zz
3 5: ð29Þ
By applying the interface conditions, we can obtain the recursion relation for the stack hybrid matrix as – layer f+1 is solid, layer f is fluid:
94
E.L. Tan / Ultrasonics 50 (2010) 91–98
ðf ;NÞ
ðf þ1;NÞ ðf þ1;NÞ ½I Hf22 H11zz 1 Hf21 ; ðf þ1;NÞ f 1 ðf þ1;NÞ f H12 ½I H11zz H22 H12zI ; ðf þ1;NÞ ðf þ1;NÞ H21Iz ½I Hf22 H11zz 1 Hf21 ; ðf þ1;NÞ ðf þ1;NÞ ðf þ1;NÞ ðf þ1;NÞ H22 þ H21Iz Hf22 ½I H11zz Hf22 1 H12zI :
H11
¼ Hf11 þ Hf12 H11zz
ð30aÞ
ðf ;NÞ H12 ðf ;NÞ H21 ðf ;NÞ H22
¼
ð30bÞ
¼ ¼
ð30cÞ ð30dÞ
Likewise, if layer f + 1 is fluid and layer f is solid, the stack hybrid matrix recursion relation can be obtained as – layer f+1 is fluid, layer f is solid: ðf ;NÞ
H11
f
f
ðf þ1;NÞ
¼ h11 þ h12 H11
ðf ;NÞ H12
¼
H21
ðf ;NÞ
¼
ðf ;NÞ H22
¼
f
ðf þ1;NÞ 1
½I h22 H11
f
h21 ;
f h12 ½I
ðf þ1;NÞ f 1 ðf þ1;NÞ H11 h22 H12 ; f ðf þ1;NÞ ðf þ1;NÞ 1 f H21 ½I h22 H11 h21 ; ðf þ1;NÞ ðf þ1;NÞ f ðf þ1;NÞ f 1 ðf þ1;NÞ H22 þ H21 h22 ½I H11 h22 H12 ;
ð31aÞ ð31bÞ ð31cÞ ð31dÞ
where f
ð32aÞ
f
ð32bÞ
h11 ¼ Hf11 Hf12It ðHf22tt Þ1 Hf21tI ; h12 ¼ Hf12Iz Hf12It ðHf22tt Þ1 Hf22tz ; f
h21 ¼ Hf21zI Hf22zt ðHf22tt Þ1 Hf21tI ;
ð32cÞ
f h22
ð32dÞ
¼
Hf22zz
Hf22zt ðHf22tt Þ1 Hf22tz :
From the above, we see that both (30) and (31) have similar form as (26), except partitioning of stack hybrid submatrices in (30) and replacement of layer hybrid submatrices by h in (31).
With the aid of stable QZ algorithms [21], the eigensolutions to (33) can be computed reliably through a(j) and b(j) (rather than ðjÞ pF ¼ aðjÞ =bðjÞ directly), where a(j) and b(j) are the diagonal elements of the triangular matrices QHAZ and QHBZ obtained by QZ factorization. Note that it is possible that b(j) ? 0 for some j, which could ðjÞ then make jpF j ! 1. In practice, such infinite eigenvalue can be avoided if one adopts the upward-bounded and downwardðjÞ ðjÞ bounded associations corresponding to P> F ¼ diagða =b Þ for < 1 ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ ¼ diagðb =a Þ for jb j < ja j, respectively. ja j < jb j and PF ðjÞ are assembled into v> ; s> and v< ; s< formThen the eigenvectors w F F F F F ing the Floquet eigenwave matrix wF as before. Furthermore, the Floquet wavenumber kzF can be deduced from (33) as (for integer m) ðjÞ
kzF ¼ ½2mp þ argðaðjÞ Þ argðbðjÞ Þ iðln jaðjÞ j ln jbðjÞ jÞ=hF :
ð35Þ
In Fig. 2, we show the dispersion relation of Floquet wavenumber kzF hF =2p (real part) versus fhF as calculated by (a) standard eigenproblem of transfer matrix and (b) generalized eigenproblem of hybrid matrix. The unit cell comprises aluminium and epoxy with thickness hAl ¼ hEpoxy ¼ hF =2, and the incident angle is 80° corresponding to longitudinal wave in nylon. The materials are characterized by the parameters shown in Table 1. To better demonstrate the numerical problem, the cell materials are also arbitrarily made lossy with complex viscoelastic constants, 00 C IJ ¼ C 0IJ iC IJ , where C 00IJ ¼ 0:01C 0IJ . From Fig. 2, it is observed that the first method may give corrupted results for fhF larger than about 10 103 m=s. This is because when fhF increases, the cell
3.2. Generalized eigenproblem Using the cell hybrid matrix, Floquet waves can be determined via a generalized eigenproblem:
ðjÞ bðjÞ ¼ HB w ðjÞ aðjÞ ; HA w F F I Hc11 Hc12 ; HB ¼ HA ¼ 0 Hc21 Hc22
0 : I
ð33Þ ð34Þ
Table 1 Material parameters.
Al Epoxy Nylon
Density (kg/m3)
Longitudinal velocity (m/s)
Shear velocity (m/s)
2700 1140 1100
6380 2770 2600
3140 1300 1100
Fig. 2. Relation of kzFhF/2p (real part) and fhF as calculated by (a) standard eigenproblem of transfer matrix and (b) generalized eigenproblem of hybrid matrix. The unit cell comprises (viscoelastic) aluminium and epoxy with hAl = hEpoxy. The incident angle is 80° corresponding to longitudinal wave in nylon.
E.L. Tan / Ultrasonics 50 (2010) 91–98
transfer matrix becomes numerically unstable and inaccuracies arise in computing kzF . Although we have considered some loss here, the numerical problem may still persist even without viscos-
95
ity. On the other hand, the second method using the generalized eigenproblem of hybrid matrix can alleviate the problem, giving uncorrupted results over the range considered. (For extreme fhF , ðjÞ one may anticipate that some jpF j ! 0; 1 or aðjÞ ; bðjÞ ! 0.) 3.3. Application for unit cell of solid and fluid
Fig. 3. Structures for unit cell composed of solid and fluid.
Without loss of generality, let us consider a unit cell of infinite periodic media that is composed of one solid layer and one fluid layer (solid–fluid cell), as shown in Fig. 3a. To treat the solid, one previous approach is to model only the longitudinal wave [2]. However, this approach is correct only when the wave is incident normally as to be shown below. Meanwhile, one can notice that the cell is equivalent to the structure that is made up by a solid layer placed between two fluid layers (fluid–solid–fluid cell), as shown in Fig. 3b. Alternatively, it is also equivalent to the structure that is made up by a fluid layer placed between two solid layers (solid–fluid–solid cell), as shown in Fig. 3c. For all the structures of Fig. 3, the generalized eigenproblem of hybrid matrix method may be applied by taking into account all acoustic waves (including shear waves in solids). Note that for Fig. 3b, since HA and HB are 2 2 matrices, the generalized eigenproblem of (33) has two eigenvalues and eigenvectors. On the other hand, for Fig. 3c, HA and HB are 6 6 matrices. In this case, because HA and HB are not full rank matrices, the generalized eigenproblem of (33) still has two eigenvalues and eigenvectors after one eliminates four fictitious solutions. As for Fig. 3a, HA and HB are 4 4 matrices and (34) needs to be modified slightly (easier alternative method is described below). In Fig. 4, we plot the relation of kzFhF/2p (real part) and fhF for (a) 0° (normal) and (b) 80° incidence. The structures are those of Fig. 3a (solid–fluid hso ¼ hfl ) considering all waves in solid and only longitudinal wave in solid, Fig. 3b (fluid–solid–fluid 2hfl1 ¼ hfl ¼ hso ) and Fig. 3c (solid–fluid–solid 2hso1 ¼ hso ¼ hfl ). The solid and fluid are aluminium and water, respectively. Here
Fig. 4. Relation of kzFhF/2p (real part) and fhF for (a) 0° (normal) and (b) 80° incidence. The markers are for the structures of Fig. 3a (solid–fluid hso ¼ hfl ): ‘’ – all waves in solid, ‘o’ – only longitudinal wave in solid, Fig. 3b (fluid–solid–fluid 2hfl1 ¼ hfl ¼ hso ): ‘+’ and Fig. 3c (solid–fluid–solid 2hso1 ¼ hso ¼ hfl ): ‘x’. The solid and fluid are aluminium and water, respectively. The incident angles correspond to longitudinal waves in aluminium.
96
E.L. Tan / Ultrasonics 50 (2010) 91–98
and henceforth, no viscosity is assumed for simplicity. The incident angles correspond to longitudinal waves in aluminium. From Fig. 4, one can see that the results for all equivalent structures of Fig. 4 match very well at normal and oblique incidence. The approach that models only the longitudinal wave in solid is correct only when the wave is incident normally. Further examination shows that in the situation illustrated by Fig. 3b or Fig. 3c, the thickness hfl1 or hso1 can be set to arbitrary portion of hfl or hso. When the thickness becomes small,
the hybrid matrix has the advantage that it is still well-conditioned as
Hf jhf ¼0 ¼
0
I
I
0
:
ð36Þ
Then one can even set the thickness to be zero, i.e. hfl1 = 0 or hso1 = 0. This will reduce Fig. 3b or Fig. 3c to the initial situation in Fig. 3a, which can now be analyzed easily by the generalized eigenproblem of hybrid matrix method [without having to modify (34)].
Fig. 5. Real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. The unit cell comprises aluminium and epoxy with hAl = hEpoxy. The incident angles correspond to longitudinal waves in aluminium.
Fig. 6. Power reflection coefficients versus incident angles at fhF ¼ 2000m=s for n ¼ 1; 2; 8 cells, cf. Fig. 5.
E.L. Tan / Ultrasonics 50 (2010) 91–98
97
Fig. 7. Real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. The unit cell comprises a water layer between the two solid layers of Fig. 5: hWater = 0.9hAl, hEpoxy = 0.1hAl.
4. Omnidirectional reflection In this section, we apply the generalized eigenproblem of hybrid matrix method to analyze the omnidirectional reflection in onedimensional phononic crystals with n cells (nN layers) as shown in Fig. 1. Let the acoustic wave be incident from layer 0, and the reflection (r0,1) and transmission (t0,nN+1) coefficient matrices be defined as
<0 ðZ >0 Þ ¼ r0;1 w >0 ðZ >0 Þ; w >nNþ1 ðZ
ð37aÞ ð37bÞ
The stack hybrid matrix can be easily obtained using Floquet waves as
" ð1;nNÞ
H
¼
v>F
v
s>F ðP>F Þn
s
#"
s>F
s
v>F ðP>F Þn
v
#1 :
ð38Þ
< < By applying the radiation condition w nNþ1 ðZ nNþ1 Þ ¼ 0, the reflection and transmission coefficient matrices can be solved explicitly in terms of stack hybrid matrix as
r0;1 ¼ ½HS s<0 v <0 1 ½v >0 HS s>0 ; ð1;nNÞ
ð1;nNÞ > v>nNþ1 1 H21 ½s0 þ s<0 r0;1 ; ð1;nNÞ ð1;nNÞ ð1;nNÞ ð1;nNÞ HS ¼ H11 þ H12 ½s>nNþ1 ðv >nNþ1 Þ1 H22 1 H21 :
t0;nNþ1 ¼ ½s>nNþ1 H22
sponding power reflection coefficients versus incident angles at fhF = 2000 m/s for n = 1, 2, 8 cells. From the figure one can see that the reflectance tends toward unity as the number of the cells increases thus achieving omnidirectional reflection. Next we investigate the effect of fluid layer on the phononic bandgap. The unit cell now comprises a water layer between the two solid layers of Fig. 5 (keeping the same cell thickness): hWater = 0.9hAl, hEpoxy = 0.1hAl. In Fig. 7 we plot both real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles. For this structure, the phononic bandgaps occur around fhF (0.6 1.5) 103 m/s and fhF (1.9–2.9) 103 m/s. It is seen here that one may obtain higher frequency bandgap by introducing a fluid layer between two solid layers in unit cell [7].
ð39aÞ ð39bÞ ð39cÞ
In Fig. 5, we plot both real and imaginary parts of kzFhF/2p as functions of fhF for a range of incident angles (0°–90° step 5°). The unit cell comprises aluminium and epoxy with hAl = hEpoxy. The incident angles correspond to longitudinal waves in aluminium. We find that it is especially convenient to determine (if any) the frequency–thickness range of omnidirectional reflection using such scan plot of the imaginary parts of Floquet wavenumbers. Recall that if the imaginary parts of all kzF are not zero, the acoustic waves will decay rapidly along ^z direction. From the figure, one can clearly see that there is a phononic bandgap around fhF (1.8 2.2) 103 m/s and another narrower one at fhF 1.15 103 m/s. To verify this, we plot in Fig. 6 the corre-
5. Conclusion In this paper, a method based on the solution to a generalized eigenproblem of hybrid matrix has been presented for stable analysis of Floquet wave propagation in one-dimensional phononic crystals with solids and fluids. The method overcomes the numerical instability in the standard eigenproblem of transfer matrix. Thus, it enables Floquet waves to be determined reliably as has been demonstrated numerically. The recursion relations of hybrid matrix for periodic multilayered structure of various solid and/or fluid phases have been formulated. Dispersion relation and omnidirectional reflection for one-dimensional phononic crystals with solids and fluids have been discussed. The frequency–thickness range of phononic bandgap has been determined conveniently based on the Floquet wavenumbers.
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] M.R. Shen, W.W. Cao, Acoustic bandgap formation in a periodic structure with multilayer unit cells, J. Phys. D: Appl. Phys. 33 (2000) 1150–1154.
98
E.L. Tan / Ultrasonics 50 (2010) 91–98
[3] D. Bria, B. Djafari-Rouhani, A. Bousfia, E.H. El Boudouti, A. Nougaoui, Absolute acoustic band gap in coupled multilayer structures, Europhys. Lett. 55 (6) (2001) 841–846. [4] A. Bousfia, E.H. El Boudouti, B. Djafari-Rouhani, D. Bria, A. Nougaoui, V.R. Velasco, Omnidirectional phononic reflection and selective transmission in one-dimensional acoustic layered structures, Surf. Sci. 1175 (2001) 482– 485. [5] D. Bria, B. Djafari-Rouhani, Omnidirectional elastic band gap in finite lamellar structures, Phys. Rev. E 66 (2002) 056609. [6] B. Manzanares-Martinez, J. Sanchez-Dehesa, A. Hakansson, F. Cervera, F. Ramos-Mendieta, Experimental evidence of omnidirectional elastic bandgap in finite one-dimensional phononic systems, Appl. Phys. Lett. 85 (1) (2004) 154–156. [7] Y. El Hassouani, E.H. El Boudouti, B. Djafari-Rouhani, H. Aynaou, Sagittal acoustic waves in finite solid–fluid superlattices: band-gap structure, surface and confined modes, and omnidirectional reflection and selective transmission, Phys. Rev. B 78 (2008) 174306. [8] W.T. Thomson, Transmission of elastic waves through a stratified solid medium, J. Appl. Phys. 21 (1950) 89–93. [9] N.A. Haskell, The dispersion of surface waves on multilayered media, Bull. Seismol. Soc. Am. 43 (1953) 17–34. [10] K.E. Gilbert, A propagator matrix method for periodically stratified media, J. Acoust. Soc. Am. 73 (1) (1983) 137–142. [11] M. Schoenberg, Wave propagation in alternating solid and fluid layers, Wave Motion 6 (1984) 303–320.
[12] A.H. Nayfeh, The general problem of elastic wave propagation in multilayered anisotropic media, J. Acoust. Soc. Am. 89 (1991) 1521–1531. [13] A.M.B. Braga, G. Herrmann, Floquet waves in anisotropic periodically layered composites, J. Acoust. Soc. Am. 91 (1992) 1211–1227. [14] B. Hosten, M. Castaings, Surface impedance matrices to model the propagation in layered anisotropic media, Ultrasonics 41 (2003) 501–507. [15] B. Collet, Recursive surface impedance matrix methods for ultrasonic wave propagation in piezoelectric multilayers, Ultrasonics 42 (2004) 189–197. [16] S.I. Rokhlin, L. Wang, Stable recursive algorithm for elastic wave propagation in layered anisotropic media: stiffness matrix method, J. Acoust. Soc. Am. 112 (3I) (2002) 822–834. [17] E.L. Tan, Stiffness matrix method with improved efficiency for elastic wave propagation in layered anisotropic media, J. Acoust. Soc. Am. 118 (6) (2005) 3400–3403. [18] E.L. Tan, Hybrid compliance–stiffness matrix method for stable analysis of elastic wave propagation in multilayered anisotropic media, J. Acoust. Soc. Am. 119 (1) (2006) 45–53. [19] E.L. Tan, Matrix algorithms for modeling acoustic waves in piezoelectric multilayers, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 54 (10) (2007) 2016–2023. [20] E.L. Tan, A concise and efficient scattering matrix formalism for stable analysis of elastic wave propagation in multilayered anisotropic solids, Ultrasonics 41 (2003) 229–236. [21] E.L. Tan, Generalized eigenproblem for acoustic wave propagation in periodically layered anisotropic media, J. Comput. Acoust. 16 (1) (2008) 1–10.