Finite element modelling (FEM) is a valuable tool for transducer construction, and may be both time- and cost-saving in practical design. This was the longterm motivation for the development of a finite element program for piezoelectric problems at UoB and CMR. The implementation and testing of a FEM code is described briefly in the present paper .
A basic function in the program is the calculation of resonance frequencies of cylindrical piezoelectric ceramic disks using FEM. The disks studied here are electroded on the major surfaces, and poled in the thickness direction, with one electrode grounded. A Ph.D. work by Guo  has been used as a basis for much of the present work, and attempts were made to reproduce some of the results reported in ref. .
Since the general finite element formulation for a piezoelectric material was first given by Allik and Hughes  in 1970, this method has been widely used for vibration analysis of many electromechanical devices at both low and high frequencies, see for example , ,  and .
In the derivation of the finite element method used here, the Hamilton principle has been applied to a general piezoelectric medium, following Tiersten . Point forces and point charges are introduced, and cylindrical coordinates (r,z) are used, due to the axisymmetric nature of the problem (where r and z are the radial and thickness coordinates, respectively). The material constant matrices of a piezoelectric ceramic are of the same form as for the crystal class 6mm. This procedure leads to a variational principle formulation. The core of the finite element method is the division of the region of analysis (in our case the half of the disk in the (r,z) plane) into nodal points and elements. 8-node isoparametric elements  are used, as in ref. . For these elements, the displacement u at any point in each local element is assumed to be a function of the displacement ûi, i=1...8 in the 8 nodes of the element, and likewise for the potential. This is achieved using 8 interpolation functions, Ni, as given in ref. .
When the disk is divided into nodes and elements, the variational formulation may be transformed to the finite element equations given on matrix form, which can be found for example in , :
where [Muu] is the mass-matrix, [Kuu] is the mechanical stiffness matrix, [Kuf] = [Kfu]^T is the piezoelectric stiffness matrix, [Kff] is the dielectric stiffness matrix, [û] is the vector of the displacement components in all of the nodes, [f] is a vector of the electric potential in all of the nodes, [F] is the force vector,and [Q] is the charge vector. The mass- and stiffness-matrices are calculated using Gaussian integration, the material constants for the piezoelectric material and the interpolation functions and their derivatives. Further explanation of the formulaes may be found for example in  or .
Mechanical and electrical boundary conditions are applied to Eq. (1) through [F] and [Q]. The finite element equations are transformed to [H]-form  , and the resonance frequencies are found by short-circuiting the electrodes (f=0) in vacuum ([F]=0). Impedance and admittance is calculated using mode superposition and structural damping  .
The program (FEMP, "Finite Element Modelling of Piezoelectric structures") is coded both in MATLAB and in Fortran 90. A significant testing of the code has been made. The first test was a comparison between calculated resonance frequencies using the MATLAB and Fortran 90 versions. The relative difference between calculated frequencies of the two programs were less than 0.2 ppm for several PZT-5A disks. Simulation results for isotropic cylinders were compared with NAFEMS-benchmarks  and with the commercial finite element program ANSYS , with good results (less than 80 ppm relative difference).
Simulated results using FEMP were also compared with the commercial finite element program ABAQUS , for PZT-5A disks with diameter to thickness ratio (D/T-ratio) of 0.1 to 20. Good agreement was found, with relative differences of less than 6 ppm between the simulated resonance frequencies, for frequency- thickness (fT) numbers in the range 0-2600. For the same disks, comparisons with simulated results reported in the Ph. D. Thesis of N. Guo  showed a considerably poorer agreement. Attempts were made to use exactly the same simulation parameters as in ref. ; however, relative differences of up to 12 % were observed. These differences have not yet been explained, although contacts have been made with Guo to identify the reasons for the poor agreeement.
Testing against other published results have also been made, where different element types than those implemented in FEMP are used. Therefore, in these comparisons a perfect agreement can not be expected, due to different convergence properties of different element types. Some of these references are: Kunkel et al.  (up to 2% relative difference), Lanceleur et al.  (up to 2%), Jensen  (up to 4%), Lerch  (up to 0.6%) and Mercer et al.  (up to 0.7 %). Comparisons with measurements of the electric impedance for a few disks have also given promising results.
A FEM program for piezoelectric ceramic disks has been implemented. A good agreement has been found with the commercial finite element program ABAQUS for resonance frequencies of several piezoceramic disks. A good agreement has also been found with the NAFEM benchmarks and the commercial finite element program ANSYS, for isotropic cylinders. A fair agreement was found with several published results where the simulation parameters were not identical. Poor agreement was obtained with results reported in ref. . The reason for this poor agreement remains to be explained, and has to be studied further.
It is planned to continue the work presented here under a Dr. Scient. fellowship by J. Kocbach, granted by the Research Council of Norway and starting autumn 1997. In this work, further testing/verification of FEMP is planned, in addition to further comparisons with measurements, implementation of better loss models, the study of more complicated transducer constructions (housing, matching layers and backing), and inclusion of radiation into a fluid medium.
 Most of the work presented here is described in: Kocbach J., "Endelig element modellering av piezoelektriske skiver", Cand. Scient. (M.Sc.) thesis, Dept. of Physisc, University of Bergen, Bergen, Norway (1996). (In Norwegian.)
 Guo N., The vibration characteristics of piezoelectric discs, Doctoral Thesis, Department of Mechanical Engineering, Imperial College of Science, Technology and Medicine, London, SW7 (1989).
 Allik H. and Hughes T.J.R., "Finite element method for piezoelectric vibration," International Journal for Numerical Methods in Engineering 2, 151-157 (1970).
 Kunkel H.A., Locke S. and Pikeroen B., "Finite-element analysis of vibrational modes in piezoelectric ceramic disks," IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control 37(4) 316 (1990).
 Lerch R., "Berechnung des Schwingungsverhaltens piezoelektrischer körper mit einem Vektorprozessor," Siemens Forschungs- und Entwicklungsberichte, Springer-Verlag 15(5), 234-238 (1986).
 Naillon M., Coursant R.H. and Besnier F., "Analysis of piezoelectric structures by a finite element method," Acta Electronica 25 4, 341-362 (1983).
 Tiersten H.F., Linear piezoelectric plate vibrations: Elements of the linear theory of piezoelectricity and the vibrations of piezoelectric plates (New York, Plenum Press, 1969).
 Zienkewitz O.C., The Finite Element Method, (McGraw-Hill, New York, 1977)
 NAFEMS, The standard NAFEMS Benchmarks, TNSB rev. 3, National Agency for Finite Element Methods and Standards, Glasgow, Scotland (October 5, 1990).
 Ansys revision 5.0, technical description of capabilities, Swanson Analysis Systems, Inc., Johnson Road, P.O. Box 65, Houston, PA 15342-0065, USA (1992).
 Abaqus/Standard Users Manual, Version 5.4, Hibbit, Karlsson & Sorensen (1994).
 Lanceleur P., Francois de Belleval J. and Mercier N., "Modeling of transient deformation of piezoelectric ceramics", IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency control 39 (2), 293-301 (1992).
 Jensen H., "Calculations for piezoelectric ultrasonic transducers", Report Risø-R-536, Risø National Laboratory, Roskilde, Denmark (1986).
 Mercer C.D., Reddy B.D. and Eve, R.A.: "Finite element method for piezoelectric media", Technical report UCT/CSIR, Applied Mechanics Research Unit No. 92, University of Cape Town, Republic of South Africa (1987).