A Global Optimization Method to Determine the Complex Material Constants of Piezoelectric Bars in the Length Thickness Extensional Mode

Abstract

Optimization methods have been used to determine the elastic, piezoelectric, and dielectric constants of piezoelectric materials from admittance or impedance measurements. The optimal material constants minimize the difference between the modeled and measured admittance or impedance spectra. In this paper, a global optimization method is proposed to calculate the optimal material constants of piezoelectric bars in the length thickness extensional mode. The algorithm is applied to a soft PZT and a hard PZT and is shown to be robust.

1. Introduction

Piezoelectric constants $d_{31}$ and $d_{33}$ are indispensable parameters when we evaluate the performances of new piezoelectric materials. A common method to determine $d_{31}$ from capacitance and admittance measurements has been described in an IEEE standard [1]. First, the dielectric constant ${ϵ}_{33}^T$ is calculated from the capacitance measured at 1 kHz. Second, the elastic compliance constant $s_{11}^E$ is calculated from the frequency of maximum conductance $f_s$ in the length thickness extensional mode. Third, the electromechanical coupling factor $k_{31}$ is calculated from $f_s$ and the frequency of maximum resistance $f_p$. The piezoelectric constant $d_{31}$ is finally calculated from ${ϵ}_{33}^T$, $s_{11}^E$, and $k_{31}$ [1].

To circumvent the need of the capacitance measurement, researchers have developed iterative and non-iterative methods to determine the material constants ($s_{11}^E$, ${ϵ}_{33}^T$, and $d_{31}$) from only the admittance measurements in the length thickness extensional mode. Smits proposed an iterative method to calculate the complex material constants from the admittances at three frequencies [2]. According to the one-dimensional model in the IEEE standard [1], the three complex admittances suffice to determine the three complex material constants [3], whose imaginary parts represent the losses in the material [4]. Smits’s method involved the manual selection of three data points on an admittance curve, which was later automated by Alemany et al. [3]. Non-iterative methods have also been proposed to calculate the material constants from particular near-resonance frequencies and the corresponding admittances [5,6].

In contrast to the iterative and non-iterative methods that use a limited number of admittance data, fitting methods are based on the whole admittance spectra and are therefore robuster. In a fitting method, the material constants are determined by fitting the modeled admittance or impedance spectra to the measurements. The optimal material constants are those that minimize the difference between the modeled and measured admittance or impedance spectra.

When the modeled admittance or impedance spectra come from the one-dimensional model of the length thickness extensional mode, the optimal material constants can be calculated with local optimization algorithms, such as the simplex algorithm [7], the Levenberg–Marquardt algorithm [8], and the Nelder–Mead algorithm [9]. Unlike the global optimization algorithms that were recently proposed for the thickness extensional mode [10] and the radial mode [11], these local optimization algorithms may not converge to the globally optimal material constants, especially when the initial guess of the material constants is inaccurate. To propose a global optimization algorithm for the length thickness extensional mode is the first motivation of this paper.

The modeled admittance or impedance spectra in the fitting methods can also come from finite element simulations [12,13,14,15,16]. Finite element simulations are able to predict the admittance at any frequency, even when different vibration modes are coupled. On the contrary, each one-dimensional model only applies to a specific vibration mode. Furthermore, one-dimensional models are only valid for samples with particular shapes and dimensions. The constraint is unnecessary in the finite element simulations. The fitting methods based on finite element simulations, however, are too complicated to be used in everyday research for most researchers.

Although there are several fitting methods based on the one-dimensional model of the length thickness extensional mode [7,8,9], no open-source software is available. Consequently, researchers have to write their own code to realize the algorithms if they cannot afford a commercial software license [17]. The second motivation of this paper is to publish an open-source code to determine the material constants in the one-dimensional model of the length thickness extensional mode.

The global optimization method proposed in this paper has the following advantages over other methods. It does not need a capacitance measurement to determine the dielectric constant ${ϵ}_{33}^T$, in contrast to the common method in the IEEE standard [1]. Unlike the iterative and non-iterative methods [2,3,5,6], our method takes full advantage of the measured admittance spectra and is less prone to experimental errors. Compared with the local optimization methods [7,8,9], the present global optimization method is able to find the globally optimal material constants even when their initial guesses are chosen randomly. Based on the one-dimensional model, the present method is more efficient and convenient than the fitting methods based on finite element simulations [12,13,14,15,16].

This paper is arranged as follows. In Section 2, we introduce the global optimization method to fit the modeled admittance and impedance spectra to the measured ones. The algorithm is applied to calculate the optimal material constants of PZT materials and the results are presented in Section 3 and discussed in Section 4. Conclusions are given in Section 5.

2. Materials and Methods

In the one-dimensional model of the length thickness extensional mode, the electrical admittance of a rectangular bar is [1]:

$$Y_{\mathrm{mod}}(\mathit{y},f)=\frac{2\mathrm{i}\pi flw}{t}[{ϵ}_{33}^T-\frac{d_{31}^2}{s_{11}^E}+\frac{d_{31}^2}{s_{11}^E}\frac{tan\left(\pi fl\sqrt{\rho s_{11}^E}\right)}{\pi fl\sqrt{\rho s_{11}^E}}],$$

(1)

where $\mathrm{i}=\sqrt{-1}$ and f is the frequency of the driving voltage. The density ($\rho$), length (l), width (w), and thickness (t) of the bar are experimentally measured. The one-dimensional model applies when $l\gg w$ and $w>3t$ [1]. The admittance also depends on an elastic compliance constant ($s_{11}^E$), a dielectric constant (${ϵ}_{33}^T$), and a piezoelectric constant ($d_{31}$) of the material. The three material constants are assumed to be complex numbers. The imaginary parts of $s_{11}^E$, ${ϵ}_{33}^T$, and $d_{31}$ are related to the elastic, dielectric, and piezoelectric losses, respectively [4]. The six-dimensional real vector $\mathit{y}$ is defined as:

$$\mathit{y}={(\mathrm{Re}\left\{s_{11}^E\right\},\mathrm{Im}\left\{s_{11}^E\right\},\mathrm{Re}\left\{{ϵ}_{33}^T\right\},\mathrm{Im}\left\{{ϵ}_{33}^T\right\},\mathrm{Re}\left\{d_{31}\right\},\mathrm{Im}\left\{d_{31}\right\})}^{\mathrm{T}},$$

(2)

where “$\mathrm{Re}$” and “$\mathrm{Im}$” represent the real and imaginary parts, respectively. The superscript “$\mathrm{T}$” represents the transpose of the vector.

The experimental admittances are denoted as $Y_{\exp }\left(f_n\right)$, where $f_n$ ($1\le n\le N$) are N frequencies around the resonance frequency of the first length thickness extensional mode. We denote the average relative error of the admittance and impedance in the one-dimensional model as

$$E\left(\mathit{y}\right)=\frac{E_G\left(\mathit{y}\right)+E_B\left(\mathit{y}\right)+E_R\left(\mathit{y}\right)+E_X\left(\mathit{y}\right)}{4},$$

(3)

where

$$\begin{array}{c}{E}_G\left(\mathit{y}\right)=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\left[\frac{G_{\mathrm{mod}}(\mathit{y},f_n)-G_{\exp }\left(f_n\right)}{G_{\exp }\left(f_n\right)}\right]}^2},\end{array}$$

(4)

$$\begin{array}{c}{E}_B\left(\mathit{y}\right)=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\left[\frac{B_{\mathrm{mod}}(\mathit{y},f_n)-B_{\exp }\left(f_n\right)}{B_{\exp }\left(f_n\right)}\right]}^2},\end{array}$$

(5)

$$\begin{array}{c}{E}_R\left(\mathit{y}\right)=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\left[\frac{R_{\mathrm{mod}}(\mathit{y},f_n)-R_{\exp }\left(f_n\right)}{R_{\exp }\left(f_n\right)}\right]}^2},\end{array}$$

(6)

$$\begin{array}{c}{E}_X\left(\mathit{y}\right)=\sqrt{\frac{1}{N}{\sum }_{n=1}^N{\left[\frac{X_{\mathrm{mod}}(\mathit{y},f_n)-X_{\exp }\left(f_n\right)}{X_{\exp }\left(f_n\right)}\right]}^2}.\end{array}$$

(7)

Here, $G=\mathrm{Re}\left\{Y\right\}$, $B=\mathrm{Im}\left\{Y\right\}$, $Z=1/Y$, $R=\mathrm{Re}\left\{Z\right\}$, and $X=\mathrm{Im}\left\{Z\right\}$ are the conductance, susceptance, impedance, resistance, and reactance, respectively. The subscripts “$\mathrm{mod}$” and “$\exp$” represent the modeled and measured quantities, respectively. The optimal material constants correspond to the vector $\mathit{y}$ that minimizes the average relative error $E\left(\mathit{y}\right)$. The minimum average relative error is calculated with the Levenberg–Marquardt modification of Newton’s method as follows [18].

(1)

Measure the admittance, density, and dimensions of the sample. Calculate the experimental conductance, susceptance, resistance, and reactance from the measured admittance.

(2)

Randomly select a set of material constants $(s_{11}^E,{ϵ}_{33}^T,d_{31})$ as an initial guess, where:

$$\begin{array}{c}5\times {10}^{-11}{\mathrm{m}}^2/\mathrm{N}<\mathrm{Re}\left\{s_{11}^E\right\}<2\times {10}^{-12}{\mathrm{m}}^2/\mathrm{N},\end{array}$$

(8)

$$\begin{array}{c}-1\times {10}^{-12}{\mathrm{m}}^2/\mathrm{N}<\mathrm{Im}\left\{s_{11}^E\right\}<1\times {10}^{-12}{\mathrm{m}}^2/\mathrm{N},\end{array}$$

(9)

$$\begin{array}{c}1\times {10}^{-9}\mathrm{F}/\mathrm{m}<\mathrm{Re}\left\{{ϵ}_{33}^T\right\}<5\times {10}^{-8}\mathrm{F}/\mathrm{m},\end{array}$$

(10)

$$\begin{array}{c}-1\times {10}^{-9}\mathrm{F}/\mathrm{m}<\mathrm{Im}\left\{{ϵ}_{33}^T\right\}<1\times {10}^{-9}\mathrm{F}/\mathrm{m},\end{array}$$

(11)

$$\begin{array}{c}-5\times {10}^{-9}\mathrm{C}/\mathrm{N}<\mathrm{Re}\left\{d_{31}\right\}<-5\times {10}^{-9}\mathrm{C}/\mathrm{N},\end{array}$$

(12)

$$\begin{array}{c}-1\times {10}^{-9}\mathrm{C}/\mathrm{N}<\mathrm{Im}\left\{d_{31}\right\}<1\times {10}^{-9}\mathrm{C}/\mathrm{N}.\end{array}$$

(13)

(3)

Define:

$$\mathit{x}=(\frac{\mathrm{Re}\left\{s_{11}^E\right\}}{{10}^{-11}{\mathrm{m}}^2/\mathrm{N}},\frac{\mathrm{Im}\left\{s_{11}^E\right\}}{{10}^{-13}{\mathrm{m}}^2/\mathrm{N}},\frac{\mathrm{Re}\left\{{ϵ}_{33}^T\right\}}{{10}^{-8}\mathrm{F}/\mathrm{m}},\frac{\mathrm{Im}\left\{{ϵ}_{33}^T\right\}}{{10}^{-10}\mathrm{F}/\mathrm{m}},\frac{\mathrm{Re}\left\{d_{31}\right\}}{{10}^{-10}\mathrm{C}/\mathrm{N}},\frac{\mathrm{Im}\left\{d_{31}\right\}{)}^{\mathrm{T}}}{{10}^{-12}\mathrm{C}/\mathrm{N}})$$

(14)

such that the absolute value of each component of $\mathit{x}$ is close to 1.

(4)

Calculate the average relative error E from Equations (3)–(7). Note that E is also a function of $\mathit{x}$. Calculate the gradient and the Hessian matrix of E and denote them as $\mathit{g}\left(\mathit{x}\right)$ and $\mathit{F}\left(\mathit{x}\right)$, whose expressions are derived in Appendix A.

(5)

Define $\mathsf{\Delta }\mathit{x}=-{[\mathit{F}\left(\mathit{x}\right)+\mu \left(\mathit{x}\right)\mathit{I}]}^{-1}\mathit{g}\left(\mathit{x}\right)$, where $\mu \left(\mathit{x}\right)=max\{-2{\mathsf{\lambda }}_{min}\left(\mathit{x}\right),0.001\}$ and ${\mathsf{\lambda }}_{min}\left(\mathit{x}\right)$ is the minimum eigenvalue of the Hessian matrix $\mathit{F}\left(\mathit{x}\right)$. $\mathbf{I}$ is the $6\times 6$ identity matrix.

(6)

Calculate the average relative errors $E(\mathit{x}+\alpha \mathsf{\Delta }\mathit{x})$ for step sizes $\alpha =1/25$, 1/5, 1, and 5. Find the minimum average relative error $E(\mathit{x}+{\alpha }_{\mathrm{opt}}\mathsf{\Delta }\mathit{x})$, where ${\alpha }_{\mathrm{opt}}$ is the optimal step size. Take $\mathit{x}+{\alpha }_{\mathrm{opt}}\mathsf{\Delta }\mathit{x}$ to be the new $\mathit{x}$.

(7)

In the inner iteration, repeat steps 4–6 until the absolute values of all components of ${\alpha }_{\mathrm{opt}}\mathsf{\Delta }\mathit{x}$ are less than ${10}^{-8}$. Then $E\left(\mathit{x}\right)$ is the locally minimum average relative error. Calculate the locally optimal material constants $(s_{11}^E,{ϵ}_{33}^T,d_{31})$ from Equation (14).

(8)

In the outer iteration, repeat steps 2–7 for 100 times with different sets of initial material constants. The globally minimum average relative error is defined as the minimum among all locally minimum average relative errors. The corresponding material constants are the globally optimal material constants.

The flow chart of the algorithm is shown in Figure 1. In the end, we check whether the optimal material constants are reasonable. For a material to be passive, i.e., energy is dissipated instead of generated in the material, the following conditions must be satisfied: [4]

$$\begin{array}{c}\mathrm{Im}\left\{s_{11}^E\right\}\le 0,\end{array}$$

(15)

$$\begin{array}{c}\mathrm{Im}\left\{{ϵ}_{33}^T\right\}\le 0,\end{array}$$

(16)

$$\begin{array}{c}\mathrm{Im}\left\{s_{11}^E\right\}\mathrm{Im}\left\{{ϵ}_{33}^T\right\}\ge {\left(\mathrm{Im}\left\{d_{31}\right\}\right)}^2.\end{array}$$

(17)

3. Results

The algorithm proposed in Section 2 is applied to a soft PZT and a hard PZT to calculate their optimal material constants in the length thickness extensional mode. The densities of the soft PZT and the hard PZT are $7.85\times {10}^3$ kg/m³ and $7.81\times {10}^3$ kg/m³, respectively [7]. Both samples have the same dimensions. The length, width, and thickness are 12 mm, 3 mm, and 1 mm, respectively [7].

The experimental admittance and impedance of the soft PZT are obtained from Figure 1 in [7] and plotted with symbols in Figure 2. Note that the data are not evenly spaced in the frequency range between 110 kHz and 140 kHz, which corresponds to the first length thickness extensional mode. For example, the reactance data are denser when $f\approx 134.5 \mathrm{kHz}$ (Figure 2h). Consequently, we use the linearly interpolated admittance and impedance at 301 frequencies evenly spaced between 110 kHz and 140 kHz as the experimental data in our calculations. The modeled admittance and impedance that best fit the experimental data are then calculated with the present algorithm and plotted with lines in Figure 2. The globally minimum average relative error is 6.5% and is found in 40 runs among 100 runs with different initial material constants. The optimal material constants are listed in

Table 1. The optimal material constants of the soft PZT and the hard PZT. ${ϵ}_0=8.854\times {10}^{-12}F/m$ is the vacuum permittivity.

	Soft PZT	Hard PZT
$s_{11}^E$$\left({10}^{-12}\right)$ m2/N)	$13.47-0.15\mathrm{i}$	$9.719-0.012\mathrm{i}$
${ϵ}_{33}^T/{ϵ}_0$	$1307-41\mathrm{i}$	$765.1-3.4\mathrm{i}$
$d_{31}$$\left({10}^{-12}\right)$ C/N)	$-133.4+3.8\mathrm{i}$	$-63.49+0.32\mathrm{i}$

and are confirmed to satisfy Equations (15)–(17). The relative differences between our results and the previous results [7] are less than 3% and 7% for the real and imaginary parts of the material constants, respectively.

Similarly, the experimental admittance and impedance of the hard PZT in the first length thickness extensional mode at 251 frequencies evenly spaced between 148 kHz and 158 kHz are linearly interpolated from Figure 2 in [7]. The original experimental data and the best fitting curves are plotted in Figure 3 with symbols and lines, respectively. The globally minimum average relative error is 10.4% and is found in 48 runs among 100 runs. The optimal material constants are listed in Table 1 and are confirmed to satisfy Equations (15)–(17). The relative differences between our results and the previous results [7] are less than 7% for the real parts of the material constants. The imaginary parts of the present optimal ${ϵ}_{33}^T$ and $d_{31}$, however, do not agree with the previous results [7]. Compared with the previous best fitting results [7], which are plotted with red dashed lines in Figure 3, our best fitting admittance and impedance are closer to the experimental results.

4. Discussion

In Section 3, the minimum average relative errors are 6.5% and 10.4% for the soft PZT and the hard PZT, respectively. Note that these average relative errors represent the relative differences between the modeled admittance and impedance spectra and the measured ones. The errors of the optimal material constants are difficult to estimate but should be far less than 10%. Otherwise, the differences between the modeled and measured resonance frequencies would have been evident in Figure 2 and Figure 3.

The average relative error of the modeled admittance and impedance spectra cannot be reduced to zero due to the following two reasons. First, the measured admittance and impedance may not be accurate. They are sensitive to the holding position of the sample. Second, the one-dimensional model (Equation (1)) only applies when $l\gg w$ and $w>3t$ [1]. When these conditions are not satisfied or when the sample is not a perfect rectangular bar, the one-dimensional model is only an approximation.

The present global optimization method depends on the range of the initial guess of material constants. If this range is larger than that in Equations (8)–(13), the algorithm may need more runs to find the globally optimal material constants. The algorithm may even find different locally optimal material constants that have almost the same average relative error. Actually, according to Theorem A1 in Appendix B, different materials may have similar admittances around the same resonance frequency. For example, if we have several materials and the material constants of the n-th material are ${(2n-1)}^2{s}_{11}^E$, ${ϵ}_{33}^T+4n(n-1)d_{31}^2/s_{11}^E$, and ${(2n-1)}^2{d}_{31}$, where $n\ge 1$, then the n-th resonance frequency of the n-th material is close to the first resonance frequency of the first material. Furthermore, the admittance spectrum of the n-th material near its n-th resonance frequency is similar to that of the first material near the first resonance frequency. This phenomenon is illustrated in Figure 4, where the first material is chosen as the soft PZT.

Although the present algorithm may converge to different optimal material constants, we can easily check whether the vibration mode in the given frequency range is the first mode. According to Figure 4, those conductance spectra with local maxima at lower frequencies outside the given frequency range should be neglected. Only the optimal material constants that correspond to a modeled admittance spectrum with the first mode in the given frequency range are the ones we need.

5. Conclusions

In this paper, we proposed a global optimization method to determine the complex material constants of piezoelectric bars in the length thickness extensional mode. The calculated optimal material constants and the best fitting admittance and impedance are consistent with the previous results [7] for a soft PZT and a hard PZT. Unlike the previous algorithms [7,8,9], the present algorithm involved 100 runs from random initial material constants and was able to find the globally minimum average relative error of the modeled admittance and impedance. Even when the initial material constants were randomly chosen from a wide range, at least 40% of runs found the globally minimum average relative error, which shows the robustness of the current algorithm. On the contrary, the local optimization methods [7,8,9] depend on well-chosen initial material constants, otherwise they may find a local minimum instead of the global minimum of the average relative error.

The present algorithm can also be modified to determine the optimal material constants in other vibration modes with one-dimensional models.