INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 39,909-922 (1996) NUMERICAL ALGORITHMS FOR SOLUTIONS OF LARGE EIGENVALUE PROBLEMS IN PIEZOELECTRIC RESONATORS YOOK-KONG YONG' AND YOUNG CHO' Department of Civil and Environmental Engineering, Rutgers University, P.O. Box 909, Piscataway, N J 08855-0909, U.S.A. SUMMARY Two algorithms for eigenvalue problems in piezoelectric finite element analyses are introduced. The first algorithm involves the use of Lanczos method with a new matrix storage scheme, while the second algorithm uses a Rayleigh quotient iteration scheme. In both solution methods, schemes are implemented to reduce storage requirements and solution time. Both solution methods also seek to preserve the sparsity structure of the stiffness matrix to realize major savings in memory. In the Lanczos method with the new storage scheme, the bandwidth of the stiffness matrix is optimized by mixing the electrical degree of freedom with the mechanical degrees of freedom. The unique structural pattern of the consistent mass matrix is exploited to reduce storage requirements. These major reductions in memory requirements for both the stiffness and mass matrices also provided large savings in computational time. In the Rayleigh quotient iteration method, an algorithm for generating good initial eigenpairs is employed to improve its overall convergence rate, and its convergence stability in the regions of closely spaced eigenvalues and repeated eigenvalues. The initial eigenvectors are obtained by interpolation from a coarse mesh. In order for this multi-mesh iterative method to be effective, an eigenvector of interest in the fine mesh must resemble an eigenvector in the coarse mesh. Hence, the method is effective for finding the set of eigenpairs in the low-frequency range, while the Lanczos method with a mixed electromechanical matrix can be used for any frequency range. Results of example problems are presented to show the savings in solution time and storage requirements of the proposed algorithms when compared with the existing algorithms in the literature. KEY WORDS piezoelectric resonators; numerical algorithms;large scale eigenvdue problems 1. INTRODUCTION T h e analysis and design of piezoelectric resonator devices present a challenging set of problems to the engineer: (1) The geometry a n d boundary conditions of the devices are complex, (2) the resonating element is a composite, consisting of metal electrodes on a piezoelectric plate or substrate, a n d in thin film devices, the piezoelectric layer is grown on a dielectric layer, a n d (3) the devices usually operate at higher modes of vibration. Hence, the analysis of these devices demands good analytical models, and detailed finite element solution of such analytical models. This is especially relevant when an accurate prediction of the unwanted modes is required. * Associate professor Graduate research assistant CCC 0029-5981/96/060909-14 0 1996 by John Wiley & Sons, Ltd. Received 3 November I994 Revised 18 April 1995 910 Y.-K. YONG AND Y.CHO The numerical models would generate a set of coupled piezoelectric eigenvalue equations of the following form: M,U Kuuu K,g+ = F (1) + + K , ~ P+ K+,+ = Q (2) where the matrices Mu,, K,, K , and K,, are the mechanical mass, mechanical stiffness,piezoelectric coupling, and dielectric matrices, respectively. The vectors u, 4, F and Q are, respectively, the mechanical displacement, electric potential, mechanical force and electric charge vectors. A static condensation' of the potential 4 in equations (1) and (2) renders a single equation written as M,U + K*u = F* (3) where K * = K,, - KUgK,&,lKg, (4) F* = F - K , + K i i Q (5) and K* is the condensed electromechanical stiffness matrix and F* is the corresponding electromechanical forcing function. The components of the electric potential vectors are recovered from equation (2) by 4 = K,-,'(Q - K& (6) A generalized eigenvalue problem can be obtained from equation (3) by assuming harmonic solutions, and setting F* equal to zero to yield (K* - W ~ M , , , ) = U0 (7) where o is the natural frequency. Equation (7) is the piezoelectric eigenvalue problem based on Allik and Hughes' finite element formulation, and has been widely used in most piezoelectric eigenvalue problems for over two decades. It is found in some recent The matrix K*, while still symmetric, does not retain the sparsity structure of the original mechanical stiffness matrix K,,, and is now very densly populated. Hence, direct solution of equation (7) requires extensive storage requirements and consequently long solution time. In an effort to keep the sparsity structure of the stiffness matrix, , ~ separates the a perturbation method was proposed by Yong and Zhang' and Boucher et ~ l . which mechanical eigenvalue problem from the piezoelectric eigenvalue problem. The piezoelectric effect was then introduced into the mechanical eigenpairs as perturbation parameters. The accuracy of the perturbation solution is predicated on the assumption that the perturbation parameters were at least an order of magnitude smaller than the mechanical stiffness. Hence, the method does not provide accurate eigenvalues for resonator materials such as lithium niobate, zinc oxide, and aluminum nitride which have stronger piezoelectric coupling coefficients than, say, quartz. Guo et ul.' proposed a scheme in which the piezoelectric eigenvalue matrix equations (1) and (2) are solved directly using the Lanczos eigensolver. Their scheme retained the sparsity structures of the finite element matrices, and provided savings in both memory and solution time when compared with the direct solution of equation (7). However, it yields a very large half bandwidth in the columns corresponding to the electrical degrees of freedom because of the large differences in the degree of freedom numbers between the electrical degrees of freedom and mechanical degrees of freedom in a node. For finite element matrices with large half bandwidths, such as those in three-dimensional cubic problems, and circular plate problems, an iterative scheme which does not require complete factorization and hence stores only the non-zero terms would provide more savings in memory; and, for very large-scale problems, in computational time too. PIEZOELECTRIC RESONATORS 91 1 2. ALGORITHMS FOR THE SOLUTION OF THE PIEZOELECTRIC EIGENVALUE PROBLEM In this paper, two sets of algorithms are presented to solve the piezoelectric eigenvalue problem without static condensation of the electric potential. In the first set, the scheme proposed in Reference 7 using a Lanczos eigensolver is enhanced by a new storage scheme for the stiffness and consistent mass matrix. The stiffness matrix storage scheme uses an equation numbering system that mixes the electrical and mechanical degrees of freedom in the assembled global matrix. The equation number corresponding to the electrical degree of the freedom for any node is sequential to equation numbers of the mechanical degrees of freedom for the same node. The storage scheme for the consistent mass matrix follows that of Reference 9. This set of algorithms is denoted subsequently as the Lanczos method with a mixed electromechanical matrix (LMEM). In the second set of algorithms, an iterative scheme denoted as the multi-mesh Rayleigh quotient iteration (MRQI) is presented. The scheme does not require complete factorization and exploits the sparse matrix structure itself. Convergence of results is significantly improved by using initial eigenpairs from a coarse mesh. Results from example problems using a three-dimensional brick element are presented. 2.1. Lanczos method with a mixed electromechanical matrix ( L M E M )for the piezoelectric eigenvalue problem The piezoelectric eigenvalue matrix equations are obtained from equations (1) and (2) by setting F and Q to zero, and assuming harmonic solution, which yields the matrix equation in the form: where If the left-hand-side electromechanical matrix is stored in a skyline scheme as proposed in Reference 7, the piezoelectric coupling K,, would have a very large half bandwidth equal to (NMDOF-MDOF-I),where NMDOF is the total number of mechanical degrees of freedom in the finite element model, and MDOF is the mechanical degrees of freedom per node. Depending on the size of the finite element model, this half bandwidth will always be orders of magnitude larger than the mean half bandwidth of the mechanical stiffness matrix K,,. Since the Lanczos eigensolver requires factorization of the electromechanical matrix, all the zero values between the non-zero values of K,, and K,, have to be stored. The half bandwidth of K,, can be reduced to the same order of magnitude as the mean half bandwidth of K,, by rearranging the equation numbers for the electrical degrees of freedom: For each node, the electrical degree of freedom number is made sequential to the mechanical degrees of freedom number for the node, so that K,,,K,, and K,, in equation (8) are all mixed together consisting of one large electromechanical stiffness matrix with a mean half bandwidth of the same order of magnitude as the mean half bandwidth of K,,. The consequence of this mixing of electrical and mechanical degrees of freedom in equation (8) is great savings in memory and computational time when compared with the previous scheme of Reference 7. The rearrangement results in zero diagonal values on the right-hand mass matrix of equation (8) which correspond with the electrical degrees of freedom. 912 Y.-K. YONG AND Y.CHO These zero 'mass' terms imply infinite eigenvalues for the electrical modes which are theorectically correct since the piezoelectric eigenvalue equations assume infinite velocity for the electrical waves. The eigenmodes which are of interest in the acoustic resonator correspond only to the mechanical degrees of freedom, hence these infinite eigenvalues are neglected, and do not present a problem in the algorithm. Pivotal problems resulting from the zero diagonal mass values can be handled with a simple modification to the Cholesky's factorization of the mass matrix. We use a piezoelectric eigenvalueequation with a mixed electromechanicalstiffnessmatrix [K] instead of equation (8): CK1 { X I = n c m { 4 (10) is the vector of mechanical displacements and electric potential. The electric potential degree of freedom is numbered sequential to the mechanical displacement degrees of freedom for a node. For an eight-node bilinear brick element with three mechanical and one electrical degrees of freedom per node, the element stiffness and mass matrices are formed by the following equations: K" = sv Me = PI' CDl CBI dV (12) p[N]'[N]dV J-V where 0 0 CBI N2.X2 0 = N2.X3 0 0 0 N1 0 0 0 N1 0 0 O N 2 0 0 0 N I O O 9x32 0 Nz 0 0 0 0 0 N 2 O ... (15) The terms Ni, [C], [el and [ E ] are the shape functions, mechanical stiffness constants, piezoelectric constants, and dielectric permittivities, respectively. Since the [ N ] matrix contains zero values in every fourth column corresponding to the electrical degree of freedom, the global mass matrix [MIalso contains zero column values corresponding to every electrical degree of freedom. 913 PIEZOELECTRIC RESONATORS Hence, [MI could be factorized using Cholesky's factorization without pivotal problems by skipping the zero columns associated with the electrical degrees of freedom. An LDLT factorization is performed on the mixed electromechanicalstiffness matrix [K].The skyline of the matrices remain the same after factorization, and could be written in the form: 1 CLI CDl WIT{ X I 3, = c641 C641T{4. (17) Equation (17) could be further transformed into a standard eigenvalue matrix equation: 1 - 1 v = [A] v (18) where v = [ 9 I T {x} and [A] = [B]TIL]-TID]-l[L]-l [64] The application of the Lanczos algorithm' transforms equation (18) to the form where V = [XI V * , [ T I is a tridiagonal matrix, and [XI is the Lanczos matrix. In general, the zero values beneath the skyline of both the stiffness and mass matrices must be stored since the factorization of these matrices would replace most of these zero values with non-zero values. The consistent mass matrix however has a special sparse matrix form with . ~ us investigate the structure of the consistent mass a fill-in pattern which is known u p r i ~ r iLet matrix of an eight-node brick element by reviewing its definition. The structure of the N matrix in equation (15) causes the element mass matrix in equation (13) to be composed of subdiagonal matrices in the form of: M', = 0, if Ii - j l # lm (22) MFi=O, i f i = 4 p , p = 1 , 2 , . . . , 8 (23) where M ; and M i are elements of Me, 1 is the difference in the node numbers associated with row i and columnj, respectively, and m is equal to the number of degrees of freedom per node which is four in this case. The structure of the mass matrix for an eight-node brick element as given by equations (22) and (23) is shown in Figure 1. The assembled global mass matrix has the same structure as the element mass matrix. Hence, equations (22) and (23) also apply to the global mass matrix. Now, we investigate the fill-in pattern in the Cholesky factorization procedure: ~ 1 =Mi',' 1 ( u.. lI = M.. - (24) i- 1 1 uii)'I2 i = 2,3,. . . ,n, i + 4p, p = 1,2,. . . (25) k=I,k#4p i- 1 k=l.k#4p where uij are the elements of the factorized upper triangular matrix U,and n is the order of the matrix. 914 Y.-K. YONG AND Y.CHO * SYMMETRIC I Non-zero terms 000 00 0 0Zero terms 32x32 Figure 1. Structure of the mass matrix for an eight-node brick element Examining equations (22) and (26), the matrix U is found to be a sparse matrix and have the same structure as M ,namely, uij=O, ifli-jl#lm forjai uii=O, i f i = 4 p , p = l , 2 , . . . (27) (28) From equations (27) and (28), we observe that the fill-in's would not happen at the locations with li - j I # lm and every fourth column and row which contain zero values can be discarded without affecting other values. Hence, we could store, inside the envelope, only the terms with I i - j I = lm for the mass matrix. With this storage scheme, the percentage savings of memory over the skyline storage is about ((m- l)/m) x 100per cent when the problem size is large. The savings in memory, for example, is about 75 per cent for a problem which has 4 degrees of freedom per node. There is also a decrease in computational time due to a smaller number of floating point operations. 2.2. Multi-mesh, Rayleigh Quotient Iteration Method (MRQI) for the piezoelectric eigenvalue problem A multi-mesh, Rayleigh Quotient Iterative Method (MRQI) is proposed to gain further savings in memory and computational time for large-scale piezoelectric eigenvalue problems. This method is advantageous if only a small set of eigenpairs is desired, and the eigenvectors could be represented on a coarser mesh of the final finite element mesh. Hence, the method is useful for the extreme eigenpairs in the low-frequency range. Rayleigh quotient iteration. The Rayleigh quotient iteration" is a revision of the symmetric QR algorithm. It is useful when a few eigenvalues and eigenvectors are desired. The Rayleigh quotient is 1 = r(u) = u*K*u u Mu 915 PIEZOELECTRIC RESONATORS 2 x 2 x 2 Coarsest Mesh 4 x 4 x 4 Coarser Mesh 8 x 8 x 8 Original Mesh Figure 2. Hierarchy of the cubic mesh applied in the iterative solution which minimizes (1 (K* - 1M)u 11. The scalar r(u) is called the Rayleigh quotient of u. If u is an approximate eigenvector, then r(u) is a reasonable choice for the corresponding eigenvalue. On the other hand, if 1 is an approximate eigenvalue, then the inverse iteration theory gives the idea that the solution to (K* - 1M)u = b will also be a good approximate eigenvector. This idea constitutes the basic algorithm of the Rayleigh quotient iteration. By starting with an arbitrary non-zero vector uo, a desired accuracy level of eigenpair solutions are obtained after a moderate number of iterations. However, the number of iterations required for convergence depends on the size of the problem, which results in excessively long solution time for large size problems. We found that by starting with an initial vector which resembles the final eigenvector, rather than an arbitrary vector uo, the required number of iterations for convergence is significantly reduced. The initial eigenpairs could be obtained from a coarse mesh of the problem, and the eigenvectors could then be represented in the final mesh by interpolation. The coarse mesh problem could be solved cheaply by the Lanczos method proposed in the previous section, and would not involve much memory since the memory allocated for the final mesh could be utilized temporarily. Figure 2 shows, for example, the use of a 2 x 2 x 2 cube mesh as a coarsest mesh and 4 x 4 x 4 cube 8 mesh problem. One could also use mesh as a coarser mesh for the original 8 x 8 ~ cube a 2 x 2 x 2 mesh as a coarser mesh for 8 x 8 x 8 mesh problem by interpolating twice. However, this will result in poor initial conditions leading to more iterations and longer solution time. Table I shows the algorithm for solving a piezoelectric eigenvalue problem using the Rayleigh quotient iteration. In this study, the convergence criterion is set as where 11 11 indicates the Euclidean norm of a vector and e4 is the convergence tolerance for the required eigensolutions. Preconditioned conjugate gradient iteration. In Table I, which shows the Rayleigh quotient iteration algorithm, the step involving the solution of the linear matrix equation Bzk+ = Muk provides for the next iterate Z k + 1. The matrix B is usually not positive definite. Hence, a conjugate gradient iterative methodlo for obtaining zk+ will not normally perform well. However, and quite fortunately, the Rayleigh quotient iteration does not need accurate solutions of z k + to yield satisfactory convergence patterns in the final eigenpairs, and a conjugate gradient iterative method could yield reasonable estimates of z k + even when B is not positive definite. A preconditioned conjugate gradient iterationlo is adapted for the solution of this linear equation. 916 Y.-K. Y O N G A N D Y. CHO Table I. Algorithm for rayleigh quotient iteration Initial vector uo is obtained from the eigenvectors of a coarse mesh For k = 0,1, . . . ,n, repeat Calculate #& = r(uk) exit the loop The convergence speed of the conjugate gradient iteration method can be accelerated by applying an appropriate preconditioning matrix. Theoretically the ideal preconditioning matrix is a completely factorized matrix of B itself. A complete factorization is however not the right choice because it defeats the primary purpose of this algorithm, namely, to use minimal amounts in the of memory storage. We note that the contribution of the piezoelectric term K,,&$K& matrix B is relatively small, and is not formed explicitly. Hence, one choice for the preconditioning matrix is the incompletely factorized matrix (K,,,,- l k M ) where the piezoelectric term is neglected. The LDLT decomposition is used in the incomplete factorization. The convergence criterion for the solution of this linear equation is set as where is the convergence tolerance. Table I1 shows the algorithm for the preconditioned conjugate gradient iteration. 3. NUMERICAL EXAMPLES AND RESULTS Numerical examples and their results are presented in this section to demonstrate the efficiency of the proposed algorithms. A cubic solid and a prismatic bar are used to show the performances of each algorithm. The two models are discretized with different fineness of meshes to generate different problem sizes. The results of the proposed LMEM and MRQI schemes are compared with (a) the Lanczos method with a condensed electromechanical matrix (LCEM) from Reference 1 and equation (7), and (b) the Lanczos method with a separated electromechanical matrix (LSEM)from Reference 7 and equation (8). All the schemes except the MRQI scheme have an efficient memory storage method described in equations (22) and (23) for the consistent mass matrix. The stiffness matrix in the LCEM scheme is a dense matrix, and hence is stored in a full matrix. The stiffnessmatrix in the LSEM and LMEM schemes is stored in a skyline format. In the MRQI scheme, only the non-zero terms of both the mass and stiffness matrices are stored. The examples use a three-dimensional, eight-node bilinear element with three mechanical, and one electrical degrees of freedom per node. The material type of the examples is PZT4 and the 917 PIEZOELECTRIC RESONATORS Table 11. Algorithm for preconditioned conjugate gradient iteration Solve for Z k + in Bzk+ = Muk where B = (K, + K.,K;tK& - L,M) Set initial vector zP+ = 0 Set ro = Muk For j = 1,. . .,n,repeat Solve P y j - = rj- for y j - where P: preconditioner Calculate flj = yJ-lrj-l/yJ-2rj-2 If j = 1, fll = 0 Calculate Pi = y j - + PjPj- If j = 1, P 1 = ro Compute aj = y:- Irj- ,/[Pj'BPj] Obtain updated vector z:+~ = z { i t + a,Pj Compute rj = rj- - ajBP, Check for convergence ; coefficients of the material are obtained from Reference 6. The convergence criteria for the MRQI scheme are set as lo-' for &+ in the Rayleigh quotient iteration (Table I), and lo-' for el in the preconditioning conjugate gradient iteration (Table 11). The relatively large parameter of 10- for is justified by the fact that the marginal return for setting smaller parameters is very low, namely, the Rayleigh quotient iteration needs only a reasonable estimate of Z k + i for the good convergence rates, and the convergence rate of the preconditioned conjugate gradient iteration is decelerated near the exact solution. ' 3. I . Cubic solid example In this example, the cubic solid model is discretized using meshes of increasing fineness starting from 3 x 3 x 3 to 16 x 16 x 16 elements. Figure 3 shows the comparison of storage requirements for the four schemes versus the number of equations (number of degrees of freedom). Some savings are gained from the LSEM scheme which uses a skyline storage of the separated electromechanical matrix when compared with the LCEM scheme which requires full matrix storage. The LMEM scheme which uses a skyline storage of the mixed electromechanical matrix (equation (10)) shows big savings in the storage requirements when compared with the two previous schemes. Since the MRQI scheme needs to store only the non-zero values of the stiffness and mass matrices, it has the smallest storage requirements. For the 16 x 16 x 16 element mesh problem with approximately 20 OOO degrees of freedom, the LSEM scheme uses about 43 per cent of the storage requirement of the LCEM scheme; while the LMEM and MRQI schemes use only 15 and 1.8 per cent, respectively. Figure 4 shows the solution-time comparisons of the four schemes when two eigenpairs are calculated. The solution time for the MRQI scheme includes time for generating initial eigenpairs from a coarse mesh. For the MRQI scheme, the 6 x 6 x 6 mesh problem used initial eigenpairs obtained by the Lanczos method from a 3 x 3 x 3 coarse mesh; and the 16 x 16 x 16 mesh problem used eigenpairs obtained by the MRQI scheme from a 8 x 8 x 8 coarse mesh, which in turn used eigenpairs solved by the Lanczos method from a 4 x 4 x 4 mesh problem. The graph indicates that the LMEM scheme shows the best performance for problem sizes smaller than the 10 x 10 x 10 918 Y.-K. YONG AND Y. CHO 3 - , ,k lo3 ; --)5 Y u) c C : : lo-’ LCEM scheme [I] --A--LSEMschemeV] MRQl scheme (proposed) LMEM scheme (PrOpOSed) I I 1o3 104 1 No. of Equations Storage requirement vs. number of equations in the cube example Figure 3. Storage requirements vs. number of equations in the cube example Solution time vs. number of equations in the cube example (2 eigensolutions) Figure 4. Solution time vs. number of equations in the cube example (2 eigensolutions) 919 PIEZOELECTRIC RESONATORS Table 111. Total number of iterations required for 2 eigenpairs Type of iteration Rayleigh quot. Problem size iter. 4X4X4 6 x 6 ~ 6 8X8X8 10 x 10 x 10 12 x 12 x 12 1 6 x 1 6 x 16 6 5 4 4 4 4 Prec. conj. grad. iter. 525 135 136 178 218 316 element mesh (5323 degrees of freedom), while the MRQI scheme is clearly superior for larger problems. The reason for the poor performance of the MRQI scheme in small size problems is that the small problem size forces the scheme to use meshes which are too coarse for good initial eigenpairs, which in turn leads to large increases in the number of iterations required for convergence of the final eigenpairs. Overall, the LCEM scheme is the worst performer. For example, in the 12 x 12 x 12 element mesh problem (8787 degrees of freedom), the solution times for the MRQI, LMEM, and LSEM schemes are respectively, about 0.6,1, and 6 per cent that of the LCEM scheme. Table 111 shows the total number of iterations required in the MRQI scheme using the given convergence criteria. The table shows that the total number of iterations in the Rayleigh quotient iteration is unchanged by the problem size (except the first two small size problems). Other iteration schemes, such as, the Gauss-Seidel method, and the conjugate gradient iteration, require more iteration steps for larger problem size. This is because the finer mesh problem uses initial conditions which are more accurate, i.e. for example, the 16 x 16 x 16 element mesh problem uses initial conditions from a 8 x 8 x 8 element mesh, while the 8 x 8 x 8 element mesh problem uses initial conditions from a 4 x 4 x 4 element mesh. The preconditioned conjugate gradient iteration, however, shows a consistent increase in number of iterations with the problem size, with the exception of the 4 x 4 x 4 element mesh problem. The poor showing in the 4 x 4 x 4 mesh problem is due to poor initial conditions from the 2 x 2 x 2 element coarse mesh. 3.2. Prismatic bar example A prismatic bar example is also conducted to compare the storage requirements and solution times for the four schemes. The following meshes of cubic shaped elements are used in Figures 5 and 6: 2 x 2 x 8,4 x 4 x 16,6 x 6 x 24,8 x 8 x 32, and 10 x 10 x 40 element meshes. Figure 5 shows the storage requirements versus the number of equations (degrees of freedom) for the four schemes. The performances in terms of memory requirements of the schemes are similar to those discussed in the previous cubic solid example. The MRQI scheme uses the least amount of memory. For the 10 x 1 0 x 4 0 element mesh, the MRQI, LMEM and LSEM schemes use, respectively, 2, 6, and 24 per cent of the memory required for the LCEM scheme. Figure 6 shows the solution times of the four schemes versus the number of equations when two eigenpairs are calculated. We observe that the LMEM and LSEM schemes exhibit an improved performance over the cubic solid problem because the mesh for the prismatic bar yields a smaller 920 Y.-K. YONG AND Y.CHO / LCEM scheme [I1 LSEMscheme MRQl scheme (proposed) LMEM scheme (proposed) lo-'1 I 1o3 I I 10' No. of Equation Storage requirements vs. number of equations in the bar example Figure 5. Storage requirements vs. number of equations in the bar example 5 : 8 v) Y @105 E --+-LCEM scheme [I] - -A-- LSEM scheme [7] MRQl scheme (proposed) LMEM scheme (proposed) F - .-s :lo4 v) lo3 lo2 10' 7 loo I I 1o3 I lo4 1 No. of Equations Solution time vs. number of equations in the bar example (2 eigensolutions) Figure 6. Solution time vs. number of equations in the bar example (2 eigensolutions) PIEZOELECTRIC RESONATORS 921 half bandwidth in the stiffness matrix, which leads to smaller storage requirements and faster solution time for these two schemes. The figure shows that the LMEM scheme is the best performer for this prismatic bar example. For the 8 x 8 x 32 element mesh, the LMEM, MRQI, and LSEM scheme require approximately only 0-6,2,and 5 per cent of the solution time needed for the LCEM scheme. For practical, very large-scale finite element models, we expect that the MRQI scheme would prove to be quite competitive with the LMEM scheme. The LMEM scheme is expected to be best performer for mesh problems with narrow half bandwidths. Overall, the LCEM scheme is the worst performer. 4. CONCLUSION Two sets of algorithms for solving piezoelectric eigenvalue problems were proposed. The performances of the Lanczos method with a Mixed Electromechanical Matrix (LMEM) and the Multi-mesh, Rayleigh Quotient Iteration Scheme (MRQI) were compared with the conventional Lanczos method using a condensed electromechanical matrix’ (LCEM), and the Lanczos method using a separated electromechanical matrix’ (LSEM). The proposed MRQI scheme shows the best performance in terms of storage requirements. In the examples conducted, significant savings in storage requirements were obtained with this method in large-size problems. In terms of solution times, the MRQI scheme is advantageous for problems with a stiffness matrix of large half bandwidth, such as that in the cubic solid example. The MRQI scheme, however, has a restriction that the desired eigenmode must be reasonably well-represented on a coarser mesh. Hence, the scheme is only well-suited for finding eigenpairs in the low-frequency spectrum of very large-scale problems. The LMEM scheme, while requiring more memory than the MRQI scheme, is much less restrictive, and could be used for finding eigenpairs in any bandwidth of frequencies. The LMEM scheme shows big savings in solution times over the existing LCEM and LSEM schemes, especially for problems with a stiffness matrix possessing a compact half bandwidth, such as that in the prismatic bar example. The LCEM scheme shows the worst performance. The LSEM scheme in which the skyline storage method is maintained through the Lanczos algorithm shows an improvement in performance over the LCEM scheme. Overall, these two methods are not competitive with the proposed LMEM scheme, or the MRQI scheme. ACKNOWLEDGEMENT Support by the F A A Center for Computational Modeling of Aircraft Structures (CMAS) at Rutgers, The State University of New Jersey, is gratefully acknowledged. REFERENCES 1. H. Allik and T. J. R. Hughes, ‘Finite element method for piezoelectric vibrations’, Int. j . numer methods eng., 2, 151-157 (1970). 2. H. S. Tzou and C. I. Tseng, ‘Distributedmodel identification and vibration control of continua: piezoelectric finite element formulation and analysis’,J . Dyn. Systems, Measurement, Control, Tran. ASME, 113, 500-505 (1991). 3. J. A. Hossack and G. Hayward, ‘Finiteelement analysis of 1-3 composite transducers’, IEEE Trans. Ultrason. Ferroelec. Freq. Controls, 38, 618-629 (1991). 4. S. K. Ha, C. Keilers and F.-K. Chang, ‘Finite element analysis of composite structures containing distributed piezoceramic sensors and actuators’,A I A A J., 30, 772-780 (1992). 5. Y.-K. Yong and Z. Zhang, ‘A perturbation method for finite element modeling of piezoelectricvibrations in quartz plate resonators’,IEEE Trans. Ultrason. Ferroelec. Freq. Control, 40, 551-562 (1993). 922 Y.-K. YONG AND Y. CHO 6. D. Boucher, M. Lagier and C. Maerfeld, ‘Computation of the vibrational modes for piezoelectric array transducers using a mixed finite element-Perturbation method’, IEEE Trans. Sonics Ultrason., su-28, 318-330 (1981). 7. N. Guo, P. Cawley and D. Hitchings, ‘The finite element analysis of the vibration characteristics of piezoelectric discs’, J . Sound Vibration, 159, 115-138 (1992). 8. T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cliffs, N.J., 1987. 9. Y.-K. Yong and Z. Zhang, ‘Numerical algorithms and results for sc-cut quartz plates vibrating at the third harmonic overtone of thickness shear’, IEEE Trans. Ultrason. Ferroelec. Freq. Control, 41, 6855693 (1994). 10. G. H. Golub and C. F. Van Loan, Matrix Computations, 2nd edn., John Hopkins Univ. Press, Baltimore, MD, 1989.

1/--страниц