Complex piezoelectric material parameters of hard PZT determined from a single disc transducer

A poled piezoelectric ceramic such as lead zirconate titanate (PZT) belongs to a C ∞ crystal class, with ten independent material constants. If we consider losses, all the material parameters can be represented in complex form as all the ten material parameters also have ten loss components associated with them. Therefore, a total of 20 independent material constants are required to obtain the entire model. In this work, we present a procedure for determining the complete reduced matrix of a piezoelectric ultrasonic disc transducer using one-dimensional resonance analysis in thickness and radial modes and parametric sweep using finite element analysis to numerically calculate the impedance while comparing it with the measured impedance. The objective is to obtain the material parameters that minimize the error between the measured impedance data and the numerical impedance curve without using any optimization techniques. We can follow the standard procedure of resonance analysis for radial and thickness modes and then conduct a sensitivity analysis using parametric sweeps in COMSOL to evaluate the unknown material values. In addition, we present the transform equations in reduced matrix form to evaluate the critically sensitive parameters concerning other consistent evaluated material constant sets. To validate the procedure, we used the impedance responses of two piezoelectric discs with diameter-to-thickness ratios of 25 and 22.5. The comparison between the numerical and experimental results shows an excellent agreement for electrical impedance curves.


Introduction
Piezoelectric ceramics are found in a variety of electronic and electromechanical devices, and there has been an increasing interest in new applications over the last few decades. In the development process of a piezoelectric device, computer simulations have become an integral part; as such a system's behaviour can be investigated with a high degree of accuracy using finite element methods (FEM). However, accurate modelling requires an accurate prior knowledge of the piezoelectric material parameters. Piezoelectric ceramics with their typical perovskite structure fall in the 6 mm symmetry class; consequently, the symmetric reduced matrix needs ten independent parameters in the linear regime: five elastic, three piezoelectric, and two dielectric parameters. Furthermore, poled piezoelectric materials exhibit nonnegligible losses, i.e. dielectric, mechanical, and piezoelectric losses. Henceforth, all the ten material parameters also have ten independent loss components associated with them. However, many manufacturers do not provide all the loss components, or they are not accurate enough; hence, the unavailability of precise material parameters poses a challenge in implementing trustworthy simulations [1].
Further, lead zirconate titanate (PZT) materials display a variety of nonlinearities and interdependencies that obscure the simple direct comparisons needed to support the user in implementing these materials into design. This can pose a challenge, particularly for high-power applications as the transducers themselves are prone to self-heating, which may alter the material parameters [2]. Even with a slight change in temperature, the behaviour can diverge significantly away from the design objectives [3]; therefore, knowing the initial linear piezoelectric material parameters and their loss components is paramount for a successful design.
Different methods have been presented in literature over the years to characterize piezoelectric material parameters accurately, and consequent improvements have been made. In many ways, the easiest method to characterize piezoceramics is by conducting a resonance analysis. This can be done by measuring electrical impedance/admittance magnitudes at critical frequencies over the relevant range showing resonances and anti-resonances using a low-amplitude frequency sweep as suggested by the IEEE standard on piezoelectricity [4]. This standard uses different analytical expressions based on geometry and poling direction of the piezo resonators for material characterization. All the material parameters can be evaluated with this method, using several specimens showing one dominant resonance due to its geometry and poling. In [5], Sherrit et al evaluated the ten independent constants, including the losses, by analyzing five standard resonance modes using five samples, each with a different geometry. However, this method still provides some drawbacks: First, it might not be possible to get all the needed geometries of the desired material from the manufacturer. Second, inaccuracies in evaluated data may arise if the required specimens come from different fabrication batches and material parameters vary from one specimen to another.
Another method to characterize piezoelectric materials is to conduct circuit analysis. In the case of high losses in the piezoelectric material, the IEEE standard recommends the circuit analysis techniques as the material constants in the impedance equations are considered to be lossless. One-dimensional (1D) circuit models were used by fitting electrical impedance to the measurement curves. Thereby, the Butterworth Van Dyke model is applicable for radial mode [6], and the KLM (Krimholtz-Leedom-Mattaei) model for thickness mode [7].
The continuous improvement of computer processing ability has recently allowed us to use numerical and iterative methods combined with FEM to find the piezoelectric constants. In [8], Pérez et al developed an FEM based optimization method to determine all the material parameters with the losses for a soft PZT material. It was further improved to reduce the complexity and computation time by Kiyono et al [9] again for soft PZT material. In [10], Rupitsch and Lerch developed an inverse method based on FEM to estimate the real material parameters also for a soft piezoelectric disc transducer, and the dissipative losses are considered according to the Rayleigh damping model. However, we could not find similar techniques involving hard PZT transducers with complex material parameters. Further, these FEM based simulations have drawbacks. They may be computationally expensive and may result in parameters not converging due to the lack of iterations [11], or the optimization algorithms converge to parameters that do not reproduce the experimental impedance if the initial guess is not very close to the final solution [12]. Finally, it poses the challenge that the iteration does not react to some less sensitive parameters.
The most commonly used transducers for ultrasonic applications found in the literature are pure hard PZT discs, and a vast majority of them rely on thickness extension mode (TEM) [13], or the fundamental radial mode (RM) [14]. However, the hard PZT disc transducer that we used, as in figure 1, shows a complicated spectrum of thickness, radial, edge and coupled vibrations. Only the RM can be clearly identified and characterized using standard methods. In the TEM the disc surface is not uniformly displaced due to the presence of spurious modes, and a one-dimensional analysis is no longer applicable [15]. We can achieve the goal of a clear thickness dominant resonance, i.e. the disc is operated as a simple piston, using large diameter to thickness ratios [16]. However, in most real-world applications, the operating frequency, power transfer capability, availability of space, and desired beam size dictate the dimensions of the transducer.
In figure 1, we see the measured impedance spectra of two specimens and the preliminary FEM model fit by using the available data provided by the manufacturer and initial guessing of other parameters. We see a difference between the impedance curves, as the manufacturers have tolerance limits in their fabrication process; even two specimens of the same material, the same geometry, and the same batch of fabrication exhibit a slight dispersion in the material properties.
To avoid this limitation, in this paper, we evaluated all the material properties using a single specimen. We used the KLM model for TEM [16] and non-iterative techniques [17] for the radial mode to accurately determine the available complex piezo parameters. In the next step, we conducted a sensitivity analysis of the unknown parameters and evaluated the most sensitive parameters for different resonance peaks using a parametric sweep in FEM. Subsequently, we used matrix inversion and equations as given in IEEE standard on piezoelectricity [4] to evaluate the critically sensitive parameters. In total, this study aims to gain a deeper understanding of the sensitivity of the different parameters and assess the complete reduced matrix with their loss components, a total of 20 material constants accurately using a single hard PZT disc specimen.

Evaluation of thickness resonance
We analyzed the TEM of a circular piezoelectric disc, on a free boundary condition having thickness t and radius a, with a high aspect ratio that satisfies diameter Φ > 20t, to have a clear radial resonance and also a dominant thickness resonance [5]. The piezoelectric disc is polarized along the thickness direction, and the two opposite planar faces are covered with silver electrodes. We applied a harmonic AC signal to the electrodes of the piezoceramic disc resonator to actuate mechanical oscillations at different resonance modes. We then measured the impedance of the piezo resonator as a function of the frequency using an impedance analyzer. The characterization of the piezoelectric, dielectric and elastic constants of the piezoelectric ceramic can be performed by resonance analysis either as real or complex quantities, depending on whether losses are ignored or taken into consideration. The impedance equation of a piezoelectric disc transducer operating in TEM is derived from the linear piezoelectric equations and the wave equations: further, The material constants in the above equations are the clamped permittivity ε S 33 , the elastic stiffness at constant electric displacement field c D 33 , the coupling factor k t , the piezoelectric pressure constant h 33 , the resonance frequency f r , anti-resonance frequency f a , Z is the impedance spectra which is also a function of frequency ω, density of the piezoelectric material ρ and geometric parameters of the disc, i.e. thickness t and disc area A.
In the IEEE standard on Piezoelectricity [4], the three material constants c D 33 , ε S 33 and h 33 are assumed to be real; however, they should be complex numbers with the imaginary part representing the loss or the out of phase component [18]. Many authors have used the complex representation to fit equation (1) and obtained good results [19,20]. Nevertheless, to investigate the influence of losses in detail, we used the KLM model. We hereby updated the method and used the resonance analysis according to the IEEE standard only to determine the real part of the material constants in equation (1), and the losses associated with them were evaluated by using KLM model.

Incorporation of transducer losses
In order to achieve the miniaturization of piezoelectric devices in some applications, such as biomedical implants, higher power density is required. The heat generation in the form of losses can limit the power density; henceforth, it is necessary to include transducer material losses for modeling of piezoelectric devices. So we included the piezoelectric material losses by introducing complex material parameters. The wellknown piezoelectric equations in contracted notation, under adiabatic and linear response conditions are used for a wide range of applications [21]: where the mechanical parameters are strain S and stress T, the electrostatic parameters are electric field E and electric displacement field D. In equation (2), the piezoelectric material parameters s E , d, and ε T can be expressed as complex numbers to represent the elastic, piezoelectric, and dielectric losses, respectively. Consequently, we can define their corresponding loss tangent as: We used the KLM model to calculate the loss factors associated with them. Figure 2(a) displays the KLM model description of a piezoelectric disc transducer operating in TEM representing it as a three-port network with frequency-dependent elements. The model includes one electric port and two mechanical ports, and standard electrical components represent all the model parameters. The circuit parameters of the KLM model are defined as follows: where Z o is the specific acoustic impedance of the piezoelectric material relating the force F and particle speed u, and c is acoustic speed. The load impedance on the front and back acoustic port is given by: , When the back acoustic port is closed by a medium, i.e. air for our measurement, the KLM model can be reduced to a twoport reciprocal network, as shown in figure 2(b). The KLM equivalent circuit can be reduced to a single matrix by concatenating the matrices of all the components. Using the KLM model, we can obtain any relation between the four remaining unknowns V, I, F, and u (voltage, current, force, and particle velocity, see figure 2(b)) using circuit theory, relating the parameters from the electric port to the acoustic port and vice versa. A detailed explanation of how to mount the matrices can be found in [22].
We first, introduced the dielectric loss straightforwardly to the KLM model by introducing the dielectric loss factor tan δ tn in the capacitance C 0 shown in equation (6) as follows [16], Further, mechanical losses are defined with the addition of an imaginary part in the form of loss tangent tan δ m to the open circuit elastic stiffness c D 33 [23], and for our piezo resonator we expressed the loss tangent tan δ m as the inverse of the quality factor, defined as in [24]: There has been disagreement in the literature about the importance of additional piezoelectric losses. According to [22], piezoelectric losses play an essential role in total losses and must also be taken into account. If both elastic and dielectric coefficients are complex, it is reasonable to consider that the piezoelectric coefficients may also be complex. In fact, Holland [20] argues the necessity of using an imaginary part for the piezoelectric constants due to the limits imposed by thermodynamical considerations. So we consider the piezoelectric pressure constant h 33 as a complex quantity and introduce the third loss tangent tan δ h [23]: This loss tangent produces a phase shift between the input electrical signal and the generated stress in the piezoelectric material. Equation (4) shows further that a complex h 33 leads to a complex turn's ratio Φ in the KLM model, which is only found in very recent publications [16]. We established the theoretical validation of all the parameters of the model, including the losses for 1D resonance analysis. In figure 3, we see the 1D analytical model of a piezoelectric disc transducer showing a proper TEM resonance. However, the hard PZT disc transducer we used did not show a clear thickness resonance, as shown in figure 4. The disc is not uniformly displaced due to the presence of spurious modes; consequently, a 1D analysis is actually no longer applicable [15]. 1D modeling of such a transducer disc with coupled resonances has been a long-standing problem, and for a disc to be operated as a simple piston in the TEM, we needed a very large diameter to thickness ratio. Nevertheless,  we investigated this limitation and conducted a 1D analysis of the chosen transducer disc having an aspect ratio of 25, despite the fact that the transducer failed to show an apparent thickness resonance.
The loss factors are closely related to the figure of merit of the transducer in terms of quality factor at resonance and anti-resonance. Figure 3 shows the definition of the quality factors at resonance Q a and anti-resonance Q b , as follows: The imaginary part of the coupling factor k t is proportional to the ratio of ∆f r and ∆f a , which is a measure of change in the half-width of the resonance and anti-resonance peak, respectively. If Q a and Q b are equal, then by definition, the imaginary component of k t is zero. However, in most cases of PZTs Q b is higher than Q a , in other words the anti-resonance operational mode has a higher efficiency. We assumed that the half-width of the resonance peaks is not affected by the presence of spurious modes and evaluated the two quality factors, as shown in figure 4.
To evaluate the imaginary parts of the material constants, we attempt to numerically fit the measured impedance of  the transducers with the 1D KLM model, while keeping the already evaluated real parts unchanged. To obtain a proper fit, it was also essential to investigate how the three different loss tangents influence the two resonance peaks and, subsequently, the quality factors.
The effect on the admittance spectra (G-spectra) by changing the mechanical loss factor tan δ m is shown in figure 5. Both the resonance and anti-resonance peaks change simultaneously with a change in tan δ m and get wider with an increase in mechanical losses; thus, leading to a lower quality factor.
The influence of the dielectric loss tan δ tn , is shown in figure 6. We varied tan δ tn , while keeping the other two loss tangents constant and checked for the quality factor at resonance Q a and anti-resonance Q b . For a disc operating in TEM, the mechanical anti-resonance corresponds to the maximum impedance, i.e. electrical resonance. Thus, by increasing tan δ tn we observe an increase in the width of the resonance peak, leading to a lower quality factor at resonance Q a , while the mechanical anti-resonance peak remained unchanged. So tan δ tn has no impact on Q b .
In a similar way, the influence of piezoelectric losses tan δ h , while keeping mechanical and dielectric losses constant is shown in figure 7. In contrast to the mechanical loss factor tan δ m , tan δ h also only affects the mechanical resonance. Moreover, increasing the losses in this case produces an opposite effect to dielectric loss. This effect is more critical in materials with high losses, such as composite transducers. Uchino  et al [25] reported this effect, which could also be masked due to the opposite effect of the dielectric losses. With our chosen transducer, we extracted the material parameters while following the above procedure, which can be seen in table 1. Figure 8 shows the magnitude of the measured impedance Z and the 1D KLM model fit. As seen in figure 8, the 1D model can only predict the TEM of the piezoelectric disc, and the fit at anti-resonance frequency is not optimal. However, the experimentally determined quality factors (figure 4) should be valid and more accurate, as the quality factor is less sensitive to the frequency dependent interpolation fit.

Evaluation of radial resonance
In the second step of our characterization process, we investigated the radial resonance of the transducer disc using noniterative methods [17] to determine further parameters. In [26], Amarande simplified the technique and developed a new technique for radial mode resonance analysis, that involves solving a system of two equations to calculate complex material parameters. For our hard PZT disc transducer, both the methods produced similar results. The electrical admittance Y, of a disc resonator oscillating in the radial mode under a sinusoidal signal, as a function of the oscillating frequency ω: where the constant A is the radius, t is the thickness, ρ is the density of the piezoelectric material. The Poisson's number σ and elastic constant c p 11 can be defined by elastic compliances s E 11 and s E 12 at constant electric field as follows: where e p 31 and further k p is the electromechanical coupling factor given by: where ε T 33 and ε p 33 are clamped and stress free permittivity respectively, and d 31 is the piezoelectric charge coefficient; and finally the function J(x) for variable x is given by: where J 0 and J 1 are Bessel functions of first kind, zeroth and first order respectively. These are infinite series, since at resonance, the value of the argument shown in equation (10) is ≈2, which makes the sum to converge very fast and with the first 10 elements of the series it is possible to calculate J 0 and J 1 accurately up to nine significant figures [21]. The characterization in the radial mode, starts with the evaluation of Poisson's number σ and a dimensionless constant η. These two parameters were calculated by a polynomial fit of the data given in the IEEE standard on Piezoelectricity [4] and were initially developed by Meitzel et al [27], We calculated σ and η with the values of a i and b i , as given in the IEEE standard and r which is the ratio of second resonance frequency f r2 and first resonance frequency f r1 . In the next step the mechanical constant c p 11 and subsequently the electromechanical coupling factor k p factor is determined as follows: Now the elastic compliances s E 11 , s E 12 and the piezoelectric constant d 31 can be can be calculated from the values of k p , η, c p 11 and ε p 33 using the non-iterative methods from equations (10) and (11) [17].

Incorporation of transducer losses
In a similar fashion to the thickness mode, the mechanical, dielectric, and piezoelectric losses were introduced by including imaginary components. However, in case of a radial resonator there are two elastic compliances, s E 11 and s E 12 , involved, and in order to separate their contribution we looked into the band width at three critical frequencies [17]: the first resonance f r1 , the second resonance f r2 , and the first anti-resonance f a1 . With our definition of bandwidth of the resonance peaks shown in figure 4, we defined them in complex form as: where ∆f r1 , ∆f r2 , and ∆f a1 are the half band width frequencies. Now we can proceed in the same way as for a lossless parameter evaluation explained in the previous section, except that now all the parameters are complex to evaluate the losses associated with all the material constants in radial mode. Finally, the value of s E 66 can be determined by using the equation presented in the IEEE standard on piezoelectricity [4] with all parameters in complex form: .
We systematically investigated the radial resonance of a hard PZT disc transducer. The radial mode resonance analysis requires to analyze three critical frequencies f r1 , f r2 , and f a1 . In [17], Sherrit et al defined the resonance and antiresonance frequencies corresponding to maxima of the real parts of Y( f ) f and fZ(f ) respectively, and not to the real part of Y(f ) and Z(f ) as defined in IEEE standard on piezoelectricity [4]. Further, they also defined the three half band frequencies ∆f r1 , ∆f r2 , and ∆f a1 we see in equation (16), using Y( f ) f instead of Y(f ). Our hard PZT disc produced almost identical results for both the definitions as the difference between these two definitions is only significant for materials with low quality factor.
We simulated the resonance spectra by giving complex material constants ε p 33 ,c p 11 and k p into equation (10). This method was then applied to generate the impedance of the PZT disc transducer around the fundamental radial resonance.  In figure 9, we see a good agreement between the measured and calculated impedance, and in table 2, we see the evaluated material constants and planar coupling factor k p with complex notation from the radial resonance analysis.

FEM sensitivity analysis
After a thorough investigation of the TEM and the radial mode of the transducer disc, we do not yet have all the material parameters of the reduced matrix. So, to conduct a sensitivity analysis of the unknown material parameters, we built a FEM model using COMSOL Multiphysics. To reduce the complexity and computational time, instead of a 3D model, we used a 2D-axisymmetric space for our PZT disc transducer.
In figure 10, we see the symmetry in the geometry of the PZT disc transducer, the boundary conditions, poling, and the mesh grid.
In COMSOL Multiphysics, we modeled the piezoelectric effect using the piezoelectric devices module, which couples the solid mechanics and electrostatics physics to combine the mechanical behaviour (equation (18)) and electric behaviour (equation (19)) of the PZT transducer.
where ρ is the density, x is the displacement, s is the stress, and F v is the volume force.
where E is the electric field, and V is the electric potential. For a small electric field and mechanical deformation, the piezoelectric effect can be modeled using strain charge form are shown in equation (20). The combined behavior of the PZT disc is modeled by coupling equations (18) and (19) into equation (20).
where S p and T q denote the mechanical strain tensor and the mechanical stress tensor in Voigt notations, respectively. E stands for the electric field intensity and D relates to the dielectric displacement. For the elastic coefficient tensor s E pq , the dielectric constants ε T mn and the piezoelectric coefficients d pm for a C∞ crystal we refer to the appendix. All the coefficients were defined in complex form, and we conducted a sensitivity analysis of the real parts and imaginary parts of the six unknown parameters s E 33 , s E 13 , s E 44 , d 33 , d 15 , and ε T 11 by taking the datasheet values provided by the manufacturer or an initial guess and performed a parametric sweep in COMSOL. Most manufacturers do not provide a complete material data sheet with complex variables. They usually offer a quality factor Q m accounting for an isotropic mechanical damping and a loss tangent tan δ accounting for an isotropic dielectric factor. Nevertheless, the initial guesses can be made in an appropriate way by examining already known parameters from radial and thickness mode resonance analysis and referring to the piezoelectric constants of similar material from literature.
Only one parameter is varied at a time while keeping the other parameters of the model fixed. Even though the numerical electrical impedance does not exactly match the experimental impedance curve with the datasheet values, it helps to identify the vibrational modes in the experimental measurement, and it can further be used to perform a sensitivity analysis. The goal is to estimate how each of the unknown parameters affects the resonance frequency and damping of the different vibration modes and the impedance values at critical points.
We numerically verified that the real and imaginary parts influence the resonances differently. The imaginary parts take into account the material losses (mechanical, dielectric, or piezoelectric losses) and influence the quality factor of the resonance peaks without any substantial influence on the resonance frequency. For example, we can see in figure 11 the influence of imag(s E 13 ) on the impedance curve, and it influences the quality factor of the CM and TEM resonance peaks. We can also observe that the frequency of the resonance peaks does not change when the imaginary part of the selected parameters is varied. The maximum of the impedance curves remains  on a vertical straight line, indicating that the resonance frequency remains unchanged in the range of variation of the parameters. So, the influence of the imaginary parts of the six unknown parameters can be neglected for this part of the analysis. Figure 12 shows the simulated spectrum across the whole frequency range, showing different resonance modes. We continued with the different vibrational modes denoted as RM1, RM2, RM3, RM4, EM, TEM, and CM in figure 12 and conducted sensitivity analysis. Further, we consider the influence on the base impedance at a frequency beyond twice that of TEMa (thickness extension mode anti-resonance) for our sensitivity analysis. Each unknown parameter was changed from −25% to +25%; always only one at a time, this means the other five unknown parameters were fixed. We conducted the sensitivity analysis of six unknown parameters (s E 33 , s E 13 , s E 44 , d 33 , d 15 , and ε T 11 ), as well as the three known that we had obtained from radial mode resonance analysis (s E 11 , s E 12 , and d 13 ). Figure 13 shows that none of the six unknown parameters has any influence on RM1 resonance frequency. 1D radial mode resonance analysis is, therefore undoubtedly valid.    In contrast, we see in figure 14 that higher radial modes are heavily influenced by s E 13 ; further, that on RM3 only s E 13 has any influence among the parameters under investigation. Figure 15 shows the influence of different material parameters on RM3, since none of the material parameters have any influence except s E 13 , so its value can be adjusted to get a perfect fit of RM3 with the measurement data.

Evaluation of real parts of the material parameters
The differentiation between d 33 and s E 33 seems to be critical. We took the analogy of the coupling factor to differentiate the influence of these two parameters. In the low-frequency limit the coupling factor is a measure of the conversion efficiency between the energy density of the applied mechanical excitation and the available energy density of the electrical signal or vice versa. The electromechanical coupling constants of the different piezoelectric resonators have been determined for the case of no loss materials [28], and in the case of the radial mode, the coupling factor is defined as: In fact, the other coupling factors k t , k 33 , k 15 and k 13 for the thickness mode, length extension mode, thickness-shear mode, and length thickness extension mode, respectively can be defined in a very similar way with analogous material parameters for different modes [29]. Further, the coupling factor increases with an increasing piezoelectric constant. In another way, it can also be calculated in terms of resonance and antiresonance frequency. In thickness mode, the real part of the coupling factor k t is defined as: So we considered RM4 and looked into the difference between the resonance frequency f RM4r and anti-resonance frequency f RM4a . In figure 16 we see how RM4 is affected by d 33 . The influence on the f RM4r is small; however f RM4a is shifted to the right with increasing value of d 33 leading to a higher difference between f RM4r and f RM4a . Further, figure 17 shows how s E 33 affects RM4. Increasing the value of s E 33 shifts both f RM4r and f RM4a to lower frequencies; however, does not affect the difference between them. We conducted similar investigation on all the six material parameters s E 33 , s E 13 , s E 44 , d 33 , d 15 , and ε T 11 respectively. Figure 18 shows how the span of the RM4 (f RM4r − f RM4a ) is affected by each parameter under investigation. We observed that only s E 13 , and d 33 have any influence on it. Since we already evaluated s E 13 by Figure 17. Influence of the real part of s E 33 on RM4. All the other material constants are kept constant and only s E 33 is varied from −25% to +25% with 5% increment. fitting RM3, d 33 can be determined from the measured span of the RM4 (f RM4r − f RM4a ). Now, as s E 13 , s E 33 and d 33 have been determined from the higher RM resonances, we considered another aspect for investigation purpose, namely the base impedance at a point far away from thickness resonance 2f TEMa (twice the antiresonance frequency of the TEM). Figure 19 shows the influence of all the six material parameters on the base impedance at 5 MHz. We see that only s E 33 and d 33 have a significant influence on it. However, the designers must be careful, as a little error in the measurement of the impedance can lead to a very large error in the calculation of the material parameters. Now, we continue and evaluate the values of s E 44 , d 15 , and ε T 11 . These three parameters did not influence the RM resonances, but the EM, CM, and TEM resonances instead. It is more precise to look into the lower EM resonances. We found that only ε T 11 has an influence, as can be seen from figure 20 showing the influence of s E 44 , d 15 , and ε T 11 on the chosen lower EM resonance.
We can further look into the higher EM resonances to identify the values of s E 44 and d 15 . It can be computationally expensive to evaluate these two parameters. Figure 21 shows   how these final two parameters affect the higher EM resonance just before the TEM.
To summarize this part of the paper, the six unknown material parameters can be divided into two categories. One category involves s E 13 , s E 33 , and d 33 which affect the resonance frequencies in a major way and need to be determined by fitting the higher RM oscillations. The second category involves s E 44 , d 15 , and ε T 11 , and these three parameters affect only the EM and CM resonances.

Evaluation of imaginary components
In the first part of the sensitivity analysis, we adjusted the real parts of the unknown parameters with an initial guess of the  imaginary parts. In the second part, we evaluate the imaginary parts of the unknown parameters by fitting the damping at different resonance modes.
In the following, we varied the imaginary parts of the remaining six unknown parameters s E 33 , s E 13 , s E 44 , d 33 , d 15 , and ε T 11 in a range of −25% to +25% from its initial guess, and we evaluated the change in the quality factor Q m at the chosen resonance peaks. Figure 22 displays the influence of imag(s E 13 ) on the resonance peak at RM3, and we observe with increasing imag(s E 13 ), Q m decreases leading to higher damping. Further, in figure 23 we see the percentage change in the quality factor Q m with the six unknown parameters. Observing that only imag(s E 13 ) influences Q m at RM3 and as expected the Q m decreases with an increasing value of imag(s E 13 ). So we adjusted it by comparing the measurement data to minimize the error in the simulation at RM3.
It can be computationally expensive to isolate the influence of imag(s E 33 ), imag(s E 44 ), imag(d 33 ), and imag(ε T 11 ) as high resolution is required to find the accurate quality factors at the EM, TEM and CM resonances; in addition, we observed imag(d 15 ) is the least sensitive parameter. So to solve these problems, we went back to the basic piezoelectric relations for TEM and the equations described by Smits [30]:  Using equation (23), the only unknown parameter d 33 is evaluated in complex form, and we only considered imag(d 33 ), as the real part of d 33 is already calculated. For a pure PZT resonator operating in TEM, the loss tangent tan δ m of the resonator corresponds to the imaginary part of the open circuit elastic stiffness c D 33 as described in equation (7), which we already determined using the KLM model from the TEM resonance analysis. It can also be expressed as: Hence, imag(s E 33 ) can be approximately calculated using equation (24). Since we already evaluated imag(s E 13 ) and imag(d 33 ), the full 9×9 [s E , ε T , d] matrix can be inverted and that will give us the 9×9 [c D , β S , h] matrix. Thereby, the value of imag(s E 33 ) may be adjusted until the value of imag(c D 33 ) determined from the inverted matrix is equal to the value of imag(c D 33 ) determined from the 1D TEM analysis. The two unknown imag(s E 44 ) and imag(ε T 11 ) do not influence imag(c D 33 ) in the inversion process.
The impact of imag(s E 44 ) and imag(ε T 11 ) can be isolated from EM oscillations, see figures 24 and 25. In a similar fashion to the real parts calculation; we can first fix the imag(ε T 11 ) by fitting Q m at a lower EM resonance with the impedance measurement and after that we can calculate imag(s E 44 ) by fitting Q m at a higher EM oscillations just before TEM.  15 ) until an acceptable value of imag(ε S 11 ) determined from the inverted matrix.

Results and discussion
In this manuscript, we proposed a methodology to evaluate the complete reduced matrix of elastic, dielectric, and piezoelectric material constants, including the losses of hard PZT 181 material, using a single disc transducer. The diameter and thickness of the sample were d ≈ 25 mm and t ≈ 5 mm . We introduced complex material parameters to account for the losses associated with each parameter. The real and imaginary parts of the material parameters influenced the resonance curves differently and were calculated in two steps. The imaginary parts cause damping without an appreciable change in the resonance frequency. Table 3 shows the evaluated material parameters for strain-charge constitutive equations, including the resonance modes and the methods used for calculation. We can convert one set of material parameters into other standard forms, depending on the boundary conditions using the transform equations mentioned in the IEEE standard on piezoelectricity [4]. We added a detailed explanation to the appendix.
We verified the accuracy of the proposed method by comparing the measured and simulated electrical impedance. Figure 26 shows the measured and simulated electrical impedance using the parameters from table 3. We observed the two resonance curves show good agreement for the entire bandwidth.
Further, figures 27 and 28 show a magnified view of figure 26, around the radial modes and fundamental TEM. We also compared the initial model fit in the magnified figures. We minimized the deviations in the final simulated resonance  curve compared to the initial model fit. However, it can be seen that still some regions around the TEM are not perfectly adjusted. Minor imperfections in the sample can lead to such discrepancies in the model fit as FEM assumes an ideal sample. Nevertheless, we can further increase the accuracy of the model fit by implementing an optimization algorithms in FEM. We followed a similar procedure for another hard PZT transducer with wraparound electrodes with 45 mm diameter and 2 mm thickness, and the results are shown in figure 29.
In our analysis, we presented the results that were obtained from a single piezoelectric disc. This accounts to the fact that,  the piezoelectric, dielectric, and mechanical constants may exhibit dispersion in their material parameters from one fabrication batch to the next and even from one sample to another. This dispersion is crucial to consider if we apply the same material constants to other specimens and particularly to other geometries to conduct FEM analysis of piezoelectric devices, and it might lead to a less precise simulation. Hence, since we are using only one specimen for the process, it provides a valuable tool for quality control in fabrication processes. Further, to characterize a new type of PZT material, we must repeat the whole process, including the resonance analysis and FEM sensitivity analysis.

Conclusion
This paper described an accurate characterization procedure of a piezoelectric transducer disc that uses 1D radial and TEM resonance analysis and FEM. We determined the complete reduced matrix of the piezoelectric material parameters, including the losses, a total of 20 independent material constants from the impedance measurement of a single sample. A comparison of the calculated impedance with the corresponding measurement showed excellent agreement. The evaluated parameters were successfully tested with transducers of two different aspect ratios. In the existing literature, an isotropic mechanical loss factor Q m and an isotropic dielectric loss factor tan δ tn are used to model piezoelectric materials, and some ignore the piezoelectric loss factor tan δ h . These assumptions can give misleading results; as we saw in our sensitivity analysis, piezoelectric material exhibits ten independent loss factors and influences different resonance modes differently.
We used the KLM model and resonance analysis to characterize the transducer in TEM by adequately accounting for the three loss factors; this method is very reliable even for high-loss materials like composite transducers. For 1D radial resonance analysis, we found the involved material constants using the non-iterative techniques found in [17]. However, the value of the evaluated losses associated with d 31 , s E 11 , and s E 12 can vary based on the technique used, especially for high-loss materials. Nevertheless, for the hard PZT181 disc transducer we used for our work, the proposed method was appropriate, and the 1D model fits the measurement optimally.
We provided an alternative method to calculate d 33 using Smits equations found in [30] and the matrix inversion technique for the least sensitive parameter d 15 . In closing, we would like to emphasize that, the true accuracy of the evaluated parameters highly depends on a well-performed measurement of the electrical impedance and geometry of the transducer disc. We did not use any optimization techniques in our method, as the analysis threshold for uncertainty in the initial guess of material constants is low. If the initial guess is not close to the final solution, it can lead the algorithm to an incorrect set of material constants. Nevertheless, we can still implement a refinement algorithm to improve the model fit further. In this work, we did not look into the computation cost of the FEM; additional investigation is required regarding the nonlinear shape functions and accuracy of the FEM to reduce the computation time and cost. Further, piezoelectric materials show various nonlinear behaviour, and they are a function of the applied electric field, applied stress, operating temperature, time, and frequency. For establishing nonlinearities, it is essential to know initial small-signal parameters correctly for model prediction, as they can significantly vary from the initial values depending on the conditions, they are subjected to during operation.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this field of study. The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgment
The University of Freiburg funded the article processing charge in the funding program Open Access Funding. The relationship described by each of the linear piezoelectric equation can be represented in matrix form as shown below: We can transform these matrices to other standard sets of constants using the transform notations from IEEE standard on piezoelectricity [4]: