A new algorithm to calculate complex material parameters in piezoelectric stacks

Abstract In this paper, the hybrid differential evolution and symbiotic organism search (HDS), is the first-time developed for general solutions of a piezoelectric stack in ultrasonic transducers. The convergence and reliability of the new algorithm are verified through comparison with corresponding data from similar previous publications and differential evolution (DE) algorithm. This study also presents and discusses the calculation results using HDS for commercial piezoelectric stacks. The Matlab HDS programs for a segmented piezoelectric (PZT) model have advanced features including its applicability to any configurations, thickness and arbitrary layer numbers of PZT. Using the novel proposed technique, there is no requirement for initial data guess, no limitations for piezoelectric stacks and the convergence rate is much faster than DE. Therefore, the HDS is promising for direct evaluation of specific aging or degradation mechanisms of ultrasonic transducers. Graphical Abstract


INTRODUCTION
Piezoelectricity has been considered as one of the most important electromechanical phenomena in which materials transform a voltage into mechanical deformation or convert mechanical loads into electric charges (C. S. Feng et al. (2021)).Piezoelectric effect has been widely applied to many applications such as transducers, actuators (C. S. Feng et al. (2021), M. T. Chorsi et al. (2019)), sensors (S.Zhang and F. Yu (2011)).It is worth mentioning that piezoelectric material and piezoelectric stack parameters are extremely important in the calculation, simulation and design of ultrasonic transducers, which are popularly used in many applying, such as cleaning, welding, and ultrasonic assisted machining (M.Sako et al. (2014)).Meanwhile, these parameters are considerably inconsistent due to the diversity of material types and they made by different manufacturers around the world.Furthermore, piezoelectric materials represent a performance loss under any changes of configurations as well as working conditions.Several studies have focused on identifying a full set of such material coefficients based on experimental data by regression techniques, each of which is applied to a specified material and has certain limitations.
Generally, finite element method and the analytical method are mainly employed for modeling piezoelectricity.Although the finite element method has been devoted since the early 1970s (Y.Kagawa (2005), M. Naillon et al. (1983), Y. Amini et al. (2015), A. Ferrari and A. Mittica (2012), S. Y. Wang (2004), A. Benjeddou (2000)) as it can deal with complex models and provide high-accuracy solutions, the computational cost is still a considerable problem.To overcome this issue, the analytical models have received much interest from many researchers worldwide.Nevertheless, designers would face problematic challenges since a full set of material parameters must be identified from samples performing in different frequency ranges, in order to input into the governing equations.Because of their frequency-dependent characteristic, the analytical models seem to be never exact enough.Consequently, designers will experience problems in employing these data because of their inaccuracy or even cannot achieve expected solutions.In such cases, a comparison between the data of modeling and the experiment is really challenged.For this reason, finding a full correct set of material parameters from experimental data is required.To tackle this problem, several attempts have been proposed to estimate the piezoelectric parameters, such as IEEE's method (American National Standards Institute, 1988), Smiths' method (J.G. Smits (1976)), nonlinear regression (K.W. Kwok et al.(1997)), and gradient methods (Y.Roh and M. S. Afzal (2018), E. Heikkola et al. (2006).However, it has been found that these methods would face significant difficulties because of the highly nonlinear nature of the piezoelectricity phenomena, thus solutions are not able to be optimized or even are trapped in some local regions.To overcome this problem, non-gradient-based methods, evolutionary algorithms (EAs), could be solutions since they do not require sensitivity information of the functions because of using random searching techniques in given design spaces.Some noticeable EAs for the estimation of piezoelectricity materials are included genetic algorithms (F.Karami and R. Morsali (2019)), differential evolution (C. S. Feng et al.(2021), T. Liu et al. (2017), G. Wang et al. (2015)), and so on.
Recently, HDS has been proposed for size and shape optimization of truss structures under multiple frequency constraints (Nguyen- Van et al. (2021)).This method has also proved its efficient performance in terms of both global and local search on topology optimization (K.D. Dang et al. (2022)).However, it has not been applied in the field of piezoelectricity materials yet.Consequently, this paper aims to provide an extension of HDS in the area of piezoelectricity for computing the 3 complex parameters of piezoelectric materials in a configuration of identical layers.Moreover, from the successful application of HDS algorithm, future studies can develop for calculating 3*n complex parameters of different properties of each layer.This is especially significant in the research and evaluation of the aging or degradation mechanisms of ultrasonic transducers in practice, which for the best of the authors' knowledge, the methods mentioned have not yet been employed.

MATHEMATICAL MODEL OF PIEZOELECTRICITY
The phase-impedance models of piezoelectric materials for different modes of piezoelectric resonators are wellknown for many researchers (S.Sherrit and B. K. Mukherjee (2007)).Focusing into thickness mode, it is possible to synthesize 3 models that are mainly used in the calculation, analysis and design of transducers, including the Martin (D. A. Berlincourt et al. (1964)), the IEEE and the Kwok models (K.W. Kwok et al. (1997), D. A. Berlincourt et al. (1964), IRE Standards (IRE, 1961).These models are briefly described as below.

Martin's model
The admittance of a piezoelectric segmented system is given by: Latin American Journal of Solids and Structures, 2023, 20(3), e482 3/13 ( ) Where: and n is number of the piezoelectric wafers, N is electromechanical transformation ratio, γ is the complex propagation constant for each segment; c, K and C 0 are respectively longitudinal velocity of propagation, longitudinal wavenumber (ω/c) and longitudinally clamped capacitance for each segment, which given by [23]:

IEEE's model
Consider a low-loss piezoceramic wafer of thickness L which is relatively small in comparison with its lateral dimension, electrode area A, and polarized in the direction of its thickness.The electrical impedance of a piezoelectric is given as follows (D. A. Berlincourt et al. (1964), IRE Standards (IRE, 1961): The impedance of the piezoelectricity plate is given as follows (K.W. Kwok et al. (1997)): Where the resistance, Re, is ( ) The reactance, Im, is given as follows: Latin American Journal of Solids and Structures, 2023, 20(3), e482  4/13 ( ) and: ( ) The electrical impedance or admittance of a piezoelectric is depended on the material parameters, particularly, clamped dielectric permittivity 33 S ε , electromechanical coupling constant for the thickness mode t k , the elastic stiffness constant at constant electric displacement 33 D c .However, the piezoelectric possess dependent on mechanical, dielectric and electromechanical losses (Dahiya, R.S., Valle, M ( 2013)).These losses can be justified by considering elastic, dielectric and piezoelectric constants, which are used in these models (1) and ( 5) as complex numbers (K.W. Kwok et al. (1997)), L.F. Brown (2000)).Those parameters are rewritten as follows: ( ) Where, tan , tan  (Martin, 1964).In this paper, the HDS algorithm based on the general Martin model will allow solving problems for the general configuration of piezoelectric materials and this content will be presented in the following sections.

OPTIMIZATION PROBLEM
This paper aims to find optimal parameters of design variables corresponding to 6 piezoelectric material parameters.
D X stands for the design variable vector and is formed as follows: The optimal function of the algorithm is mean square error (MSE) of the amplitude of impedance (A) in logarithm scale and phase angular (ϕ) and can be expressed as follows: Where, nos is the number of data points, k A  and k φ  are the kth measured data of the amplitude of impedance and phase angular, respectively; k A and k φ are the kth corresponding simulated data of the amplitude of impedance and phase angular, respectively.The optimization process is established to find optimal parameters by minimizing the root mean square error (MSE) using the following objective function: Where, α =1 is the unit normalization in the unit of 1/Ω and the cost function f(X) used in the optimization problem, is minimized by using the HDS and DE which is briefly expressed in the next section.

THE HYBRID DIFFERENTIAL EVOLUTION AND SYMBIOTIC ORGANISMS SEARCH
The HDS algorithm suggested by Nguyen- Van et al. (2021) is hybrid from two classical algorithms differential evolution (DE) and symbiotic organisms search (SOS).Motivated by a symbiotic relationship in which one organism receives benefits both from mutualism and commensalism, HDS has proved its good performance in terms of good accuracy and fast convergence rate for many benchmarks and structural optimization.However, its application to fitting parameters of piezoelectricity has not been investigated yet so far.
Four main stages of HDS method, including initialization, mutation, crossover, and selection, are given as follows: • Initialization Random NP individuals are generated in a defined search space.Every i th individual (i = 1, 2, …, NP) is a vector of D design variables and is defined as follows: Latin American Journal of Solids and Structures, 2023, 20(3), e482  6/13 [ ]( ) 0 , , 0,1 , 1, 2,..., where NP stands for individual numbers, D denotes the number of design variables, L, H correspondingly symbolize for the lower and upper boundaries, rand[0; 1] is an arbitrary number in the [0,1] interval, g denotes the g th generation/population.The i th individual vector is expressed as follows.

Mutation
This stage aims to obtain the diversity of the population.Firstly, a mutant vector g i V is created by using the following mutation schemes are "DE/rand/1", "SOS", "DE/best/1" and "DE/best/2".These four operators are briefly summarized as follows: ) DE/best/1: ( ) Where R 1 ; R 2 ; R 3 ; R 4 and R 5 are random differential numbers chosen from [1, 2, 3, …, NP]; F is a random number between 0.4 and 1; MU denotes ( ) To switch from global search to local one, an adaptive parameter is used and given as follows: / 1 Where f m and f b denote the mean and best objective function values in the previous population, respectively.If

• Crossover
This stage helps to obtain better optimal solutions by using a crossover operator.Concretely, the i th trial vector

• Elitist selection
This technique aims to generate elite individuals for the next generation by choosing the best ones in a combined group of the new individual g i V and the current individuals g i X .This new population has a size of 2NP, but only NP best individuals are chosen for the next generation.
Latin American Journal of Solids and Structures, 2023, 20(3), e482 7/13 It should be noted that this study focuses mainly on the accuracy of the engineering solutions, only global search operators are used, the local search in HDS can be eliminated.Thus, the modified flowchart of HDS is illustrated in Figure 1.

Comparison of mathematical models
In this section, MATLAB codes are programmed based on mathematical models of Martin (Equations 1-4), IEEE (Equation 5), and  to calculate the phase-impedance spectra of the four samples namely PVDF, P(VDF-TrEE), 1-3/PZT/Epoxy Composite, and Lead Metaniobate.The data for material properties and constants of these materials are shown in Table 1 and Table 2, respectively.The simulation results will be use in the following sections.ε is vacuum permittivity and is equal to 8.85×10-12 (F/m).Figure 2 shows modeling data of impedance magnitude and impedance phase of the three models.It should be pointed out that these three models give the same results corresponding to the same input data from Table 2.This is reasonable because IEEE and Kwok are specified models of Martin's general solution with: (i) one piezoelectric ring and (ii) the thickness of the ring is relatively small (K*L<<π).For this reason, the following sections will employ the model of Martin (Equations 1-4) to calculate the material parameters of a piezoelectric stack with large thickness installed in ultrasonic transducers for welding, cleaning and ultrasonic-assisted machining.

Optimization
This section shows the idealness of HDS in determining material parameters by the results obtained from two MATLAB regression programs, one is based on HDS algorithm with Martin's models, the other uses DE algorithm.The first stage uses these programs to calculate the material constants from Kwok's work based on the simulation data of Kwok's model in the previous sections.The stage is conducted to evaluate the convergence rate of HDS and DE with the expectation that HDS will give good results for all materials.In the second stage, the programs are applied to calculate the complex material constants for real experimental data of a commercial PZT material (Beijing Ultrasonic, 2023).
In both cases, the number of generations, population, and other parameters in HDS are the same and can be seen in Table 3.It should be noted that HDS uses a random search technique, thus its results in each run may different, but only the best optimal results are reported.The optimal results of HDS and DE are reported in Table 5.It shows that HDS gives a good performance since it provides the same results as those in Kwok's work (K.W. Kwok et al. (1997)).However, the convergence history illustrated in Figure 3 and Figure 4 show that the convergence rate of HDS is much better than DE, then HDS's calculation number is 2-3 time less than the DE's calculation number.
The average algorithm convergence speed of HDS reached 5350 -6500, much faster than DE which is 10000-15000.The expectation optimal value (cost function f(X)) is set at 0.001.This is reasonable since the regression function is performed on the data set calculated from Kwok's, so there is no effect of noise as in the data from the experiment.
The comparison results in Table 5 and Figure 4 show that the regression method with a newly developed HDS algorithm for Martin's general model gives complete agreement with all surveyed piezoelectric materials.This allows using a new method to determine the piezoelectric material parameters with arbitrary configuration from measurement data with high accuracy and reliability.

HDS algorithm applying in experimental data
The experimental data for this stage is achieved from the authors' Ministry of Education and Training project (B2020-TNA-02).The experiment was set up to measure the impedance-phase spectra of 4 stacked PZT rings by the Keysight E5061B-3L5 ENA impedance analyzer.The PZT rings are stacked by a pair of very thin wires (0.005 mm diameter).The PZT ring dimensions and its material properties are shown in Table 6.
Latin American Journal of Solids and Structures, 2023, 20(3), e482 10/13  As the mention above, the program using HDS is applied to find variables  7.
The optimal results are provided in Table 8.The convergence history is illustrated in Figure 5.
No doubt to say that HDS gives a great performance since the data of simulation with corresponding optimal constants fitted well with those from experiments.
Latin American Journal of Solids and Structures, 2023, 20(3), e482 11/13  Similar to the above results, the convergence rate of HDS algorithm is much better than DE algorithm (the left side of Figure 5).However, the expectation optimal value is 8.4.The cause of this difference is the "non-smooth" or noise of experimentally measured data.Nonetheless, the results presented in the right side of Figure 5 show that the results from HDS regression algorithm are ideal and promising.

CONCLUSION
In this study, the HDS regression method is successfully applied for the estimation of the material parameters of piezoelectric materials with extreme accuracy for the first time.The development of a nonlinear regression algorithm based on HDS to calculate the main coefficients of piezoelectric materials has many outstanding advantages over previous methods, including (i) applicable capacity to all piezoelectric material configurations with varied thickness and unlimited layers which is particularly popular in transducers; (ii) the calculation of material parameters requires only impedance-phase-frequency spectra data in the resonant-anti-resonance frequency domain and an initial choice for data or any parameter is not necessary; (iii) the HDS can be applied to piezoelectric rings which are currently installed in the transducers, as a result, degradation of the material can be calculated as well as for the calculation and design of suitable boosters or horns, and these contents will be mentioned in the next studies.

ACKNOWLEDGMENT
This research is supported by project B2020-TNA-02.

Figure 1 .
Figure 1.The modified flowchart of HDS

Figure 2 .
Figure 2. Comparison the Martin model with the IEEE and Kwok's models

Figure 3 .Figure 4 .
Figure 3. Convergence history of HDS and DE for the data form Kwok's work PZT material (bjultrasonic.com) from the 4 stacked PZT rings.The lower and upper bounds of these parameters are shown in Table

Figure 5 .
Figure 5.The rate of convergence of the HDS and DE at the left; at the right is comparison of calculated impedance -black line and phase -red line, solid for HDS and dashed line for DE with experimental data (|Z|, black point makers and phase, red point makers. Kwok et al. (1997)Kwok model in which the factor loss tangent are tan , tan k m δ δ and tan e δ are employed to build the mathematic model, are basically the explicit representation of the Equation5in IEEE.These equations are convenient for programming nonlinear regression problems to determine the complex coefficients of piezoelectric materials in Equations 18-20, Kwok's nonlinear regression method then can overcome the limitations of methods provided by IEEE and PRAP (K.W.Kwok et al. (1997)).Equation1-4 in the Martin's model can determine complex material constants for all values of n.The IEEE model is a special case of Martin model when n = 1 and the thickness is relatively small or in other words, |K|*L<< π k m δ δ and tan e δ are the electromechanical coupling, elastic and dielectric factor loss tangent, respectively.Latin American Journal of Solids and Structures, 2023, 20(3), e482 5/13

Table 2 .
Material constants of piezoelectric materials

of the convergence and the accuracy of HDS In
this stage, the Matlab regression program HDS is firstly tried for regression of the six optimal values of of the four materials namely PVDF, P(VDF-TrEE), 1-3/PZT/Epoxy Composite and Lead Metaniobate.The target data for fitting are obtained from the simulation in the previous section.During the optimization process, the boundaries of six design variables are provided in Table4.

Table 4 .
Boundaries of design variables inKwok's model

Table 5 .
Optimal results obtained by HDS and DE

Table 6 .
Several parameters of PZT rings

Table 7 .
Bounds of design variables inKwok's model