Numerical results

In this section, we firstly conduct numerical results to validate the conver-gence for our finite element model. Secondly, we report the numerical compar-isons with our structure-preserving method and the traditional Arnoldi method for solving the GEP (11) to demonstrate the accuracy and efficiency of the pro-posed eigenvalue solver. All computations are carried out in MATLAB 2010b on a HP workstation with an Intel Quad-Core Xeon X5570 2.93GHz and 60 GB main memory, using IEEE double-precision floating-point arithmetic.

To make our numerical computation reliable, first, we scale the displacement field q and potential field ϕ by 10⁻⁵ and 10⁵, respectively, and scale the mass density ρ accordingly. The entries of the stiffness matrices K^qq, K^qϕ, K^ϕq and K^ϕϕ, and the mass matrix M^qq are about the same order (close to O(1)) after the scaling. Secondly, we reorder the interior nodes so that the matrix M1 in

(a) Sparse pattern of matrix M1 without permutation from locally refined meshes with mesh length p/80.

(b) Sparse pattern of the permuted matrix M1 from locally refined meshes with mesh length p/80.

Figure 3: Sparsity of matrix M1.

(13) has narrower band structure. As a result, the LU factorization of M1 can be computed easily and solutions of the linear systems in (25) and (26) can be obtained efficiently from (29) and (30). The sparsity patterns of the matrices M₁ obtained from FEM descritization on a locally refined mesh are shown in Figure 3.(a) without nodal reordering, and shown in Figure 3.(b) with nodal reordering. A mesh that are locally refined twice near the interface over an uniform mesh is shown in Figure 5.

The configuration of our computational domain shown in Figure 2 is as follows. The domain width AB and height CD are set to be 10⁻⁶and 3× 10⁻⁶, respectively. The ratio of the electrode width EF versus the domain width is set to be ¹₂ and the ratio of the electrode thickness DE versus the domain height is

1

15. The material constants of LiTaO3and LiNbO3are taken from measurements obtained by Kushibiki, Takanaga and Sannomiya [18]. Also, it has been shown that the viscous damping coefficient κ1 ≈ O(10⁻⁸) in the 10 KHZ operation range and κ1 ≈ O(10⁻¹⁰) in the MHZ operation range for a family of PZT materials [23, 26]. In general, κ1 depends on the operation frequency ω. The viscous damping coefficient is extrapolated to GHZ operation range according to the reciprocal rule κ₁∝ _ω¹ [7]. In our numerical studies, the viscous damping coefficient κ₁is set to be 10⁻¹⁴ and the mass damping coefficient κ₂is taken as 1− κ1 to account for the effect from the electrode weight.

4.1. Accuracy and convergence of finite element approximation

Firstly, we show that our finite element model gives accurate results in pre-dicting the center of stopping band of LSAW on the filters with aluminum elec-trodes on top of piezoelectric substrates 36^o YX-LiTaO₃ and 64^o YX-LiNbO₃. The dispersion diagrams of the attenuation constant α and the propagation constant β associated with the eigenvalue λ(ω), that is most close to −1 on

−0.05 0 0.05

(a) Dispersion diagram of LiTaO3 near the stopping band

(b) Dispersion diagram of LiNbO3near the stopping band

Figure 4: Dispersion diagrams

the complex plane for frequency ω around the stopping bands, are shown in Figure 4(a) and Figure 4(b) for crystals 36^o YX-LiTaO3 and 64^o YX-LiNbO3, respectively. A typical shear wave displacement associated with the eigenvector of the computed eigenvalue is shown in Figure 5.

Figure 5: A shear wave displacement.

In order to measure the convergence of the eigenvalues, tests over three successively refined meshes with initial mesh size h = ₂₀^p are performed. The dimensions of matrices M₁and M₂associated with uniform meshes and locally refined meshes are list in second and third columns of Table 1, respectively. We set λ_[i,ξ]to be the eigenvalue obtained from meshes with mesh length p/(10×i).

Here the index ξ = “u” and ξ = “ℓ” denote the mesh is uniform without and with local refinement, respectively. Using λ_[16,u] as exact value, the convergence of eigenvalues can be verified from|λ[16,u]−λ[i,u]| and |λ[16,u]−λ[i,ℓ]| for i = 2, 4, 8.

−1.5 −1 −0.5 0 0.5 1

Real part of eigenvalue

Imaginary part of eigenvalue

λ_3,I

Figure 6: Distribution of eigenvalues for 64^oYX-LiNbO3 at frequency 2.180 GHZ.

The values of|λ[16,u]− λ[i,u]| and |λ[16,u]− λ[i,ℓ]| (i = 2, 4, 8) for 36^oYX-LiTaO3

and 64^o YX-LiNbO3 at ωs and ωe are shown in Table 2 where ωs and ωe are the frequencies for which the stopping band starts and ends, respectively. From Table 2, it can be seen that the accuracy of eigenvalue is increased as the mesh length being reduced and is improved by using locally refined meshes.

Moreover, it is known that the wave propagation velocity is about 4112m/s for 36^o YX-LiTaO3 and about 4478m/s for 64^oYX-LiNbO3. Since the domain width p = 1× 10⁻⁶, clearly, the center of the stopping band is about 2.056 GHZ and 2.239 GHZ, respectively. We compute the center of the stoping band by averaging ωsand ωeon different mesh lengths and show it in Table 3. Obviously, one can see that the central frequency is monotonically converged to a constant when the mesh length is reduced. The numerical error from our finite element simulations is less than 0.2% and 1.2% for 36^oYX-LiTaO₃and 64^oYX-LiNbO₃, respectively.

4.2. Comparison of Algorithm 1 and traditional Arnoldi method

From Tables 2 and 3, we already show that the accuracy of the computed eigenvalues and central frequency of stopping band obtained from the locally refined mesh with mesh length p/80 is almost the same as those obtained from uniform mesh with mesh length p/160. Therefore, in the following numerical computations, we only consider the coefficient matrices in the GEP (11) that are generated by the finite element discretization on the mesh that is locally refined twice over the uniform mesh with mesh length p/80.

Let the pair (λk,I, λk,O), k = 1, . . . , N , denote the reciprocal pairs of eigen-values of (11) where λk,I and λk,O lie inside and outsideU, respectively.

Fig-Uniform Local refine mesh length M₁ M₂ M₁ M₂

p/20 3554 183 4968 183

p/40 14548 363 17192 363

p/80 58856 723 63960 723

p/160 236752 1443

Table 1: Dimension information of matrices M1and M2 obtained from FEM discretization

36^o YX-LiTaO3 64^o YX-LiNbO3

ω (GHZ) ωs= 2.028 ωe= 2.075 ωs= 2.177 ωe= 2.257

|λ[16,u]− λ[2,u]| 0.0222 0.0161 0.0955 0.0797

|λ[16,u]− λ[2,ℓ]| 0.0178 0.0141 0.0857 0.0751

|λ[16,u]− λ[4,u]| 0.0076 0.0056 0.0299 0.0458

|λ[16,u]− λ[4,ℓ]| 0.0036 0.0042 0.0166 0.0329

|λ[16,u]− λ[8,u]| 0.0016 0.0015 0.0060 0.0087

|λ[16,u]− λ[8,ℓ]| 0.0001 0.0006 0.0007 0.0016

Table 2: The values of|λ[16,u]− λ[i,ξ]| for different mesh lengths at frequencies ωsand ωe.

mesh length h_u=₄₀^p h_ℓ=₄₀^p h_u= ₈₀^p h_ℓ= ₈₀^p h_u=₁₆₀^p fc

LiTaO3 2.05246 2.05222 2.05206 2.05191 2.05183 LiNbO3 2.21623 2.21464 2.21385 2.21305 2.21305

Table 3: Computed center frequency fc(GHZ) of stopping bands of LSAW on various meshes.

Here, huand h_ℓdenote the mesh length of meshes without and with local refinement, respec-tively.

ure 6 displays the eigenvalues {λ1,I, . . . , λ9,I, λ1,O, . . . , λ5,O} of the LiNbO3 at frequency ω = 2.180 GHZ in which reciprocal pairs (λk,I, λk,O) for k = 1, . . . , 5 close toU may be of interests. The notation O(λ) represents the set of all the eigenvalues that cluster at the origin of the complex plane. Suppose 2N eigen-values nearU are desired. The Arnoldi process in (20)–(21) for GTSHIRA is set to restart if the desired eigenpairs are not convergent when the dimension j of the subspace span{bY_j} grows more than 5N. The number of restarting Arnoldi process is denoted by “#Iter” in the following.

A standard iterative approach for solving the GEP (11) is to apply Arnoldi method on the equation directly. However, the reciprocal property of the eigen-values is not guaranteed to be preserved in the computation. For Algorithm 1, based on the (S + S⁻¹)-transform, if λ and µ are the eigenvalues of (11) and (16), respectively, then λ and µ satisfy the relation µ = λ + λ⁻¹. As a result, we can obtain the k-th reciprocal pair (λk,I, λk,O ≡ 1/λk,I) by solving the al-gebraic equation µk = λk,I + λ⁻¹_k,I after the k-th eigenvalue µk of Kz = µN z is computed. Hence, the reciprocity is automatically preserved. Two numer-ical comparisons on preserving reciprocal property, between Algorithm 1 and traditional Arnoldi method, are listed in the following where eigenvalues of 64^o YX-LiNbO₃at frequency ω = 2.180 GHZ are computed.

• Traditional Arnoldi method does not guarantee that half of the computed eigenvalues lie inside of the unit circle and the others are outside. For example, when we use Arnoldi method to compute four eigenvalues (i.e., 2N = 4) of (11) which are near −1, the four convergent eigenvalues are λ1,I, λ1,O, λ2,I and λ3,I. Clearly, the reciprocal property of eigenvalues is lost.

• Suppose one wants to compute the five reciprocal pairs (λk,I, λk,O) for k = 1, . . . , 5. As shown in Figure 6, no matter what shift value τ is chosen, there exists some k^⋆ ∈ {1, . . . , 5} such that the eigenvalues in O(λ) are closer to τ than the eigenvalue λk^⋆,I or λk^⋆,O. As a result, the eigenvalue λk^⋆,I or λk^⋆,Owould not be discovered by the Arnoldi method.

For example, if we take a shift value τ =−2.89, then the desired reciprocal pair (λ_5,I, λ_5,O) nearU can not be discovered by the Arnoldi method. On the contrary, in Algorithm 1, according to the relationship µ = λ + λ⁻¹, the eigenvalue µ ofKz = µN z is far away from the shift value τ +1/τ when λ is closed to the origin. Naturally, Algorithm 1 will not converge to those unwanted eigenvalues inO(λ). As a result, all the desired eigenvalues can be discovered more easily by Algorithm 1 than the traditional Arnoldi method. Our numerical results in Table 4 show that not only all the desired eigenvalues are found by Algorithm 1, even when the number of desired eigenvalues is set to 2N = 18, it also converges much faster than the traditional Arnoldi method. In fact, it only takes two restarting steps for Algorithm 1 to converge for all cases shown in Table 4. In addition,

from the rightmost column of Table 4, one can see that all the computed eigenvalues indeed preserve the reciprocal property. On the contrary, the reciprocity of the convergent eigenvalues obtained by Arnoldi method are diminished about 3 significant digits.

From the above comparison, Algorithm 1 preserves the reciprocal property of the eigenvalues of the GEP (11) effectively. For measuring the accuracy of Algorithm 1, let us define the relative residual of an eigenpair (λ, u) of (11), where u = [ψ_i, ψ_l]^⊤, as following:

where∥ ∗ ∥F is the Frobenius matrix norm. The maximal relative residuals of the ten desired eigepairs for 64^o YX-LiNbO3 with various frequency are shown in Figure 7. From these numerical results, one can see that the eigenpairs pro-duced by Algorithm 1 possess high accuracy in terms of relative residual error.

Next, let’s compare the efficiency of Algorithm 1 and traditional Arnoldi method. The CPU times in computing ten desired eigenpairs (i.e., 2N = 10) for 64^o YX-LiNbO₃ with various frequencies by using Algorithm 1 and tradi-tional Arnoldi method are shown in Figure 8. On average, Algorithm 1 only takes 476 seconds of CPU time to compute the desired eigen pairs for all fre-quency ω in the search range. Obviously, the proposed Algorithm 1 is more efficient compared to the traditional Arnoldi method which takes 527 seconds of CPU time to get all the desired eigenpairs.

Method 2N Computed eigenvalues #Iter max{|λk,Iλk,O− 1|}

Arnoldi

Table 4: Convergent reciprocal pairs and the associated errors of reciprocity for 64^o YX-LiNbO3at frequency 2.180 GHZ versus different eigensolvers with various “2N” which denotes the number of interested eigenvalues. Here, ˜λk∈ O(λ) for k = 1, . . . , 4.

2.14 2.16 2.18 2.2 2.22 2.24 2.26 2.28 2.3 2.32

Figure 7: Maximal relative residu-als of the ten desired eigenvalues pro-duced by Algorithm 1.

2.14 2.16 2.18 2.2 2.22 2.24 2.26 2.28 2.3 2.32

Figure 8: CPU times for computing ten desired eigenpairs by the tradi-tional Arnoldi method and Algorithm 1.

5. Conclusions

In this paper, we have modeled the leaky surface acoustic wave propaga-tion on a simple resonator with an interdigital transducer (IDT) where elec-trodes are arranged periodically on piezoelectric substrates (PZT) such as 64^o YX-LiNbO3and 36^oYX-LiTaO3. The energy conservation equation (4) is dis-cretized by finite element method (FEM) on a single cell domain with proper periodic boundary conditions as shown in Figure 2. Equation (4) is discretized on locally refined meshes in order to increase the accuracy of our numerical solutions. Our FEM simulation for predicting the center frequency of the stop-ping bands of the resonator is convergent and accurate within an error about 1% compared to experimental data as shown in Tables 2 and 3. For computing the dispersion diagram near the center of stopping band of the resonator, we transform the GEP (11) into the TPQEP (2) to reveal the important reciprocal relationship of the eigenvalues in which the eigenvalues appear in reciprocal pairs (λ, 1/λ). The TPQEP (2) is then solved by GTSHIRA so that the reciprocal relationship of the eigenvalues can be automatically preserved. Our numerical results show that the traditional Arnoldi method converges slowly and fails to preserve the reciprocal property of the eigenvalues near the unit circle. On the other hand, the proposed structure-preserving method in Algorithm 1 not only converges to those eigenpairs faster than that of the traditional Arnoldi method but also possesses high accuracy in terms of relative residual error. Furthermore, the reciprocal property of the eigenpairs are kept nicely under machine precision.

Searching a good crystal cut of various PZTs for high frequency filter de-sign based on LSAW is important. Our numerical studies here show that the dispersion diagram of a resonator with a prescribed crystal cut on its PZT substrate can be computed accurately and efficiently by discretizing the model equation on locally refined meshes and solving the resulted GEP by the proposed structure-preserving method in Algorithm 1. As a result, the computation time

in searching effective crystal cuts can be shorten and the computed dispersion diagrams can be more accurate and reliable.

Acknowledgments

This work is partially supported by the National Science Council, the Taida Institute of Mathematical Sciences, the Center for Advanced Study in Theo-retical Sciences and the National Center for TheoTheo-retical Sciences in Taiwan.

The corresponding author C. T. Wu thanks the support from National Science Council in Taiwan under the grant number 99-2115-M-009-001.

References

[1] H. Allik and T. Hughes. Finite element method for piezoelectric vibration.

Int. J. Numer. Methods Eng., 2:151–157, 1970.

[2] M. B. Angel, M. I. Rocha-Gaso, M. I. Carmen, and A.V. Antonio. Surface generated acoustic wave biosensors for detection of pathogens: A review.

Sensors, 9:5740–5769, 2009.

[3] S. N. Atluri and M. Nakagaki. Computational methods for plane prob-lems of fracture. Computational Method in Mechanics of Fracture, North-Holland Publishing Co, Amsterdam:170–227, 1986.

[4] B. A. Auld. Acoustic Fields and Waves in Solids. Krieger Publishing Company, Malabar, Florida, 1990.

[5] W. Bond. The mathematics of the physical properties of crystals. BSTJ, 22:1–72, 1943.

[6] M. Buchner, W. Ruile, A. Dietz, and R. Dill. FEM analysis of the reflection coefficient of SAWS in an infinite periodic array. In Proc. IEEE Ultrason.

Symp., 371-375, 1991.

[7] C. Cai, H. Zheng, M. S. Khan, and K. C. Hung. Modeling of material damping properties in ansys. Inernational ANSYS conference Proceedings, 2002.

[8] C. K. Campbell. Surface Acoustic Wave Devices for Mobile and Wireless Communications. Academic Press, INC., 1998.

[9] D. P. Chen and H. A. Haus. Analysis of metal-strip SAW gratings and transducers. IEEE Trans. on Sonics and Ultrasonics, 32, No. 3:395–408, 1985.

[10] L. C. Chin, V. V. Vardan, and V. K. Vardan. Hybrid finite element for-mulation for periodic piezoelectric arrays subjected to fluid loading. Int. J.

Numer. Meth. Engrg., 37:2987–3003, 1994.

[11] E. K.-W. Chu, T.-M. Hwang, W.-W. Lin, and C.-T. Wu. Vibration of fast trains, palindromic eigenvalue problems and structure-preserving doubling algorithms. J. Comput. Appl. Math., 219:237–252, 2008.

[12] P. Grisvard. Elliptic Problems in nonsmooth domains. Pitman Publishing INC., 1985.

[13] A Hilliges, C. Mehl, and V. Mehrmann. On the solution of palindramic eigenvalue problems. In Proceedings 4th European Congress on Com-putational Methods in Applied Sciences and Engineering (ECCOMAS), Jyv¨askyl¨a, Finland, 2004.

[14] M. Hofer, N. Finger, G. Kovacs, J. Sch¨oberl, S. Zaglmayr, U. Langer, and R. Lerch. Finite-element simulation of wave propagation in periodic piezo-electric SAW structure. IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 53 No. 6:1192–1201, 2006.

[15] T.-M. Huang, W.-W. Lin, and J. Qian. Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration on fast trains. SIAM J. Matrix Anal. Appl., 30:1566–1592, 2008.

[16] C.F. Ipsen. Accurate eigenvalues for fast trains. SIAM News, 37, 2004.

[17] M. Koshiba, S. Mitobe, and M. Suzuki. Finite-element solution of periodic waveguides for acoustic waves. IEEE Trans. Ultrason. Ferroelectr. Freq.

Control, 34, No. 4:472–477, 1987.

[18] J. Kushibiki, I. Takanaga, and T. Sannomiya. Accurate measurements of the acoustical physical constants of LiNbO3 and LiTaO3 single crystals.

IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 46:1315–1323, 1999.

[19] R. Lerch. Simulation of piezoelectric devices by two- and three-dimensional finite elements. IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 37, No.

2:1990, 1990.

[20] W.-W. Lin. A new method for computing the closed-loop eigenvalues of a discrete-time algebraic Riccatic equation. Linear Alg. Appl., 96:157–180, 1987.

[21] D.S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Structured polyno-mial eigenvalue problems: Good vibrations from good linearizations. SIAM J. Matrix Anal. Appl., 28:1029–1051, 2006.

[22] D.S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Vector spaces of linearizations for matrix polynomials. SIAM J. Matrix Anal. Appl., 28:971–

1004, 2006.

[23] J. E. San Miguel Medina, F. Buiochi, and J. C. Adamowski. Numerical modeling of a circular piezoelectric ultrasonic transduer radiating in water.

ABCM Symposium Series in Mechatronics, 2:458–464, 2006.

[24] S. Mitobe, M. Koshiba, and M. Suzuki. Theoretical determination of equiv-alent circuit parameters for reflective SAW metallic gratings. Electr. Com-mun. Jpn, 70, No. 6:37–45, 1987.

[25] M. Mohamed, EL Gowini, and W. A. Moussa. A finite element model of a mems-based surface acoustic wave hydrongen sensor. Sensors, 10:1232–

1250, 2010.

[26] G. Nadar, E. C. Nelli Silva, and J. C. Adamowski. Effective damping value of piezoelectric transduer determined by experimental techniques and numerical analysis. ABCM Symposium Series in Mechatronics, 1:271–279, 2004.

[27] R. V. Patel. On computing the eigenvalues of a symplectic pencil. Linear Alg. Appl., 188:591–611, 1993.

[28] U. R¨osler, D. Cohrs, A. Dietz, G. Fisherauer, W. Ruile, P. Russer, and R. Weigel. Determination of leaky saw propagation, reflection and coupling on LiTaO3. IEEE Ultrasonics Symposium, pages 247–250, 1995.

[29] C. Schr¨oder. A QR-like algorithm for the palindromic eigenvalue problem.

Technical report, Preprint 388, TU Berlin, Matheon, Germany, 2007.

[30] C. Schr¨oder. URV decomposition based structured methods for palindromic and even eigenvalue problems. Technical report, Preprint 375, TU Berlin, MATHEON, Germany, 2007.