Analytical approach to piezoelectric model synthesis with the use of Cauer’s method for system design

Background Piezoceramic materials have unique property which enables direct and bilateral conversion between mechanical and electrical energy. This ability facilitates significant miniaturisation of technology and opens many opportunities in design of new actuators and energy harvesters. Mathematical modelling of piezoelectric modules is notoriously hard due to complex constitutive equations defining mechanical and electrical energy conversion. Methods The article presents research on a new synthesis method based on the Cauer’s method. Mechanical damping is introduced with the use of Rayleigh’s approximation. A discrete electromechanical model is formed based on the Mason’s piezoelectric model. The proposed approach allows modelling of piezoelectric systems based on a set of characteristic frequencies. The method allows a more general approach to the problem of modelling new systems, as opposed to application-oriented methods seen in literature. A non-standard model analysis method using edge graphs and structural numbers is also verified as a potential alternative for matrix-based method. The authors compare their precision and computation requirements. Results The structural method of mechanical model analysis gave identical results as the reference matrix method. However, the non-classical algorithm took much longer to calculate and was using more memory. The electromechanical model analysis has shown an error of 5% in comparison to resonance frequencies taken from a reference plate specification. The calculated magnitude of displacement was well above the capability of a 3.5mm thick piezoelectric plate. Conclusions The synthesis method presented in this paper allows synthesizing piezoelectric cascade models based on limited information in form of characteristic frequencies. Currently this method allows a coarse approximation of the real piezoelectric parameters with limited number of inputs. The additional method of analysis based on structural numbers offers a promising alternative to matrix calculations but requires a more thorough investigation of the computational power required to determine whether it can compete with existing algorithms.


Introduction
One of the purposes of mathematical modelling in case of physics is to give a clear and systemic picture of relations and energy flow between elements of a system.Physicists use mathematical modelling to describe phenomena and problem solve the mechanics of a wide variety of systems. 1 Mathematical synthesis supported by previous empirical analysis of underlying phenomena helps in creating models which describe real behaviour of materials and components with mathematical correlations.Deep understanding gained by thorough analysis can be later used in development of highly sophisticated systems that are used for precise control of various system variables. 2][9][10] Laplacian and Lagrangian algebra rely on polynomial and exponential equations to systemically represent components of force and energy having influence on successive derivatives of displacement or charge.Polynomial equations are among the main set of tools in control theory 2,[7][8][9][11][12][13][14][15][16] and were used during the development of signal filters composed of RLC (Rresistor, Lcoil, Ccapacitor) components. 17,18Wilhelm Cauer was among the pioneers of electrical model synthesis, creating algorithms based on transfer function polynomials, defining large range of electrical systems by using mathematical decomposition.Cauer synthesis 17,18 relies on decomposing transfer function polynomials into continued fractions.The set of coefficients of a continued fraction represents values defining each system component properties (like impedance or capacitance).This procedure allows any transfer function of a complex system to be broken down into a series of components represented by a lumped model. Moreo, because of existing analogies between electrical and mechanical elements, this method can be applied to both electrical and mechanical models.[7][8][9]11 Piezoelectric materials exhibit properties that can be described with a combination of interactions between both electrical and mechanical components.19 A set of constitutive equations for piezoelectricity 12 has been established to describe electromechanical coupling for simple and inverse piezoelectric effect.Because of the complex set of dimensional relations, described by tensors, mathematical modelling of piezoelectric components is being focused mainly on the analytical aspect of this phenomenon.Researchers largely focus on modelling simple piezoelectric components with lumped models 13,14 and matrix tensor equations 20 or try to approach the problem using FEM (finite element method) analysis.21 There are also analytical methods that can be labelled as "non-classical" which rely on AI (artificial intelligence) algorithms or graphical representations of the system like "edge graphs". Many f the developed models are used on a case-by-case basis where system variables defining piezoelectric behaviour are fully or partially known or they are being empirically determined.15,[21][22][23][24][25][26][27] There is a research gap when it comes to the subject of a full piezoelectric system synthesis where majority of the initial parameters are unknown.This article is focused on the development of a new method for piezoelectric system synthesis where the only initial information comes in the form of established resonant and anti-resonant frequencies.Presented method allows a more general approach to the problem of modelling new systems, as opposed to application-oriented methods focused on modelling systems for very specific use REVISED Amendments from Version 1 • The abstract has been extended to highlight the main merit of the manuscript as a response to both the reviewers' suggestions.• I've also put stronger emphasis on the highlight of the novelty presented in this article in the introduction section.
Corrections have also been done further down to better clarify the purpose of introducing new algorithms and models in the research work.• More English publications have been added to make it easier to understand the non-classical method used in the article as per reviewer request.• A fix was made to the content of the "Algorithm testing" paragraph regarding the source of the resonant frequencies chosen as a reference point.• More additions were made to clarify the background of the comparison.
• Another paragraph was added, explaining a test conducted to further analyse the computational load of both algorithms for systems with a larger number of degrees of freedom.Also, a memory criterion was introduced to further analyse the computational load of each algorithm by not relying solely on just the computation time criterion.• Results of the second, more thorough study were added to further explore the efficiency of a non-classical method versus the matrix method.• Conclusions were modified to account for the results of another efficiency study and to clarify the conclusions made on the basis of the electromechanical model study.• Another algorithm repository was added with all the algorithms used in the second algorithm efficiency study.
• Figures 6 and 7 were added to the manuscript, along with Tables 4 and 5. Previous Figure 6 has been changed to Figure 8.
• New data repository with upgraded algorithms has been added with a CC0 license.
Any further responses from the reviewers can be found at the end of the article cases, proposed by cited researchers.The algorithm based on Cauer synthesis method is used to construct an initial discrete mechanical model which serves as the base for piezoelectric model synthesis.Mechanical damping is introduced to the model with the help of Rayleigh's method. 28The transition between purely mechanical model and a discrete electromechanical model of the piezoelectric plate is done with the combination of Mason's model and analogies between mechanical and electrical components.Because the initial information is restricted only to a desired set of characteristic frequencies, there are many unknown piezoelectric parameters which have to be determined in this process.To solve this problem, authors have chosen an approach where most of the parameters are being approximated by using the obtained mechanical properties of an already synthesized system and a set of equations derived from the constitutive equations of piezoelectricity.
Part of this article has been devoted to a study on the potential alternative for matrix-based analysis of mechanical models with damping, by employing edge graphs and structural number algebra. 7,8,11,29,30,31The topic of analysis of mechanical systems is indirectly related to the topic of synthesis of piezoelectric systems.The studied analysis methods are used to verify synthesized piezoelectric systems.The article presents a comparison of results obtained with both methods on a variety of mathematical models.A study was also carried out to test the performance of computer algorithms developed on the basis of the non-classical method against standard matrix-based computer algorithms.Other methods based on edge graphs bare close naming resemblance like the ones used in 15,22 which may cause confusion.As a matter of clarification, this article is focusing on the method using edge graphs in combination with structural number algebra.Structural numbers visually resemble matrix calculation as each structural number construct consist of rows and columns of numbers.The number of rows and columns, however, is tightly connected to the shape of an edge graph which a given structural number relates to.Each node in an edge graph and all its connections are represented by a single structural number vector (e.g., |1 2 3 4| for a node connecting four edges numbered from 1 to 4).A structural number is then created by multiplying such node vectors by themselves according to structural number algebra.Arithmetic and logic operations based on structural number algebra were described in depth by professor Stanisław Bellert in his book 29 and in articles. 30,31More details can also be found in previous article published by authors. 32The article also directly relates to adaptations of this method used by other researchers for analysis and synthesis of various systems. 2,11,16

Mechanical model synthesis
The synthesis process begins with a selection of predetermined resonant and anti-resonant frequencies.The set of frequencies may be freely defined according to future applications or may result from a vibration spectral analysis of the object of interest.In this method, the identified set of frequencies is the only input information.It is used to generate a transfer function polynomial, which can later be transformed into a chain fraction equation.Cauer's algorithm mentioned in Refs.8, 11, 33, 34 was used to create four different types of mechanical models depending on the allocation of coefficients in the chosen transfer function polynomial.Depending on the polynomial, the cascade model resulting from the chain fractioning process can be constrained on one end or not constrained at all.Additionally, the model can be excited by either a kinetic or dynamic force.A combination of those factors creates four different types of models with their own transfer functions (detailed examples of each type and associated transfer function structures are given in Refs.8, 11).The model chosen for the purposes of this article is a cascade model with one of its ends constrained and which is being excited by the dynamic force.The choice is motivated by the close resemblance of such a model to a fixed piezoelectric stack, where individual plates can be described as successive stages of the model.According to the models proposed in Refs.8 and 11, a hypothetical transfer function for an infinite degree cascade model constrained on one end and excited with a dynamic force takes the factorial form: where: T(s)response of a system, Hresponse amplification factor, ωconsecutive frequencies (even are resonant and odd are anti-resonant), s -Laplacian variable.
The factorial form is transformed into a standard polynomial form to define the coefficients for successive powers of the equation: where: After the initial calculation of the polynomial coefficients, the transfer function is converted into a continued fraction by using the Cauer's method in its first form. 17The method involves breaking down a polynomial into a chain fraction using simple mathematical operations.The first form of Cauer synthesis is implemented by dividing the coefficient of the greatest power of a numerator by the coefficient of the greatest power of a denominator.The resulting fraction is being multiplied by the denominator and subtracted from the polynomial function, leaving the rest, which is then inverted.This process is repeated until all powers of the transfer function are reduced: I i pfraction extracted by division of parameters next to the highest power "s", R i prest of the original polynomial after subtraction.
A continued fraction resulting from the process can be used to determine the values of discrete component parameters in the synthesized system: Obtained values can be used in either mechanical or electrical cascade systems which gives this method a high degree of flexibility.An example of an infinite stage mechanical system resulting from the calculations have been shown in Figure 1.
The newly obtained mechanical model consists of springs and inertial components which can be used to determine the stiffness and mass/density of the piezoelectric material.Damping elements can be added to the model to account for energy dissipation which is an inherent property of piezoelectric materials.The mechanical quality factor is often used to indicate piezoelectric mechanical damping ability.The factor is inversely proportional to the amount of damping generated by the material and can be derived from Rayleigh's method. 28,35However, damping introduces complex derivatives into the equations of motion for any mechanical or electrical system.Every system can be considered under the case of subcritical, critical, and over-critical damping.Depending on the strength of damping elements, a system can be in a Figure 1.A mechanical cascade system composed of an infinite number of stages composed of mass "m" and spring "c" pairs, synthesized from the transfer function polynomial with the Cauer's method.
state of damped oscillations (subcritical damping) or the oscillations may stop immediately after the harmonic excitation disappears (critical damping), they do not occur (over-critical damping).Because of that, there are three separate integral solutions that describe each case respectively.A resonance (additive interference of harmonic oscillations) occurs only under the case of subcritical damping.In order to limit the range of possible solutions, it was assumed that all models would only be considered under subcritical damping, thus eliminating cases where vibrations cannot resonate.The decision was made because the project was focused on piezoelectric applications where vibration generation was their main purpose.
Another assumption limited the range of solutions to those with a Rayleigh damping coefficient related to stiffness.According to literature, 24 most cases of piezoelectric damping are restricted to damping with respect to element stiffness only.The two assumptions simplify the term for the damping coefficient.Based on the obtained equations, it is then possible to determine the damping parameters for each stage in relation to the spring stiffness: where: To verify the model in its current state, an intermediate electromotive force was also added at each stage.Its role was to simulate the actuation of piezoelectric plates by an electrical stimulus, similar to. 15,25The electromotive force acts as a simplified force generated by the electromechanical coupling used to describe the conversion of energy conversion between the mechanical and electrical part of a piezoelectric plate.The mechanical model with added damping and electromotive forces is shown in Figure 2.

Mechanical model verification
The synthesized model had to be verified after the addition of new damping elements and additional forces.An analysis of the model response was carried out in terms of its magnitude along the frequency axis.Two analysis methods were used to verify the existing mechanical model.The classical matrix method served as a reference for an additional method using structural numbers. 2,11,16,29,30,31The use of a second method is considered as part of the ongoing research 32 aimed at evaluating the accuracy and computational load of structural numbers in comparison to the matrix calculation.Current state of the algorithm based on structural numbers should still not be considered as final as there are more improvements to be made in terms of its compatibility with computation software.Starting with the classical method, a set of equations of motion was created to describe relations between parts of the model.To avoid overcomplicating the example, only two degrees of freedom were considered.The equations of motion for a two-stage model derived from the Lagrangian equations were as follows: Figure 2. A discrete cascade mechanical model with added damping elements.The model masses are marked with "m" symbols, springs with "c" symbols, dampers with "b" symbols and electromotive forces with "G" symbols.The displacement of each mass was marked with "x".
Under the assumption into account, that the system is only considered in the case of subcritical damping, an integral of a second-order differential equation 36 for the variable displacement x can be described by the equation: where: A i1real part of the i th element amplitude response, A i2imaginary part of the i th element amplitude response.
Combining (9) with (10) and separating the components with sine and cosine coefficients, then reordering the equations by combining all the coefficients standing next to consecutive amplitudes, yields a coefficient matrix: Based on the obtained coefficient matrix (11) and a vector of actuating forces, a vector of amplitudes can be calculated.To obtain absolute displacement, the sine and cosine coefficient must be combined using Pythagoreans' theorem: To calculate phase-shift of each stage along the frequency spectrum, a tangent between the real and imaginary component of displacement can be calculated: The set of derived equations (12,13) allows the amplitude response and phase shifts to be determined for any cascade system built using standard discrete elements representing model inertia, stiffness, and damping.Calculations are relatively quick and easy when considering systems with a small number of degrees of freedom.However, calculations become much more complex and computationally intensive for systems with a large number of degrees of freedom.This is an effect of the number of elements being multiplied in the determinant of a matrix.The numbers get big really fast because of the accumulation of the powers when all the matrix components are being combined.This problem prompted the search for alternative solutions for rapid analysis of more complex systems.
One of the alternatives to matrix calculations is being developed by the researchers at the Silesian University of Technology. 2,7,8,11,16The non-classical method of system analysis involves the use graphs and structural numbers to illustrate relations between forces in a system.The structural number method offers an alternative solution for deriving equations of motion by using a graphical representation of relations in the system.Each edge of the graph corresponds to a single discrete relation between two elements of a system, represented by nodes.All the edges and nodes in the graph are assigned numbers called "structural numbers".These numbers form structures similar to matrices based on edge connections with nearby nodes.Using the algebra of structural numbers, it is then possible to derive a system of equations similar to equations of motion.A detailed description of this method in the context of solving mechanical cascade systems without damping was presented in Ref. 32.The structural number algebra is an entirely separate method of performing mathematical operations and therefore requires extensive explanation.It is strongly recommended to consult the book of the original author of this method 29 or related articles 30,31 for the details on its usage.
In the context of this article, several necessary changes had to be introduced to account for damping elements in the model.The basic method does not allow a simple decomposition of the sine and cosine coefficients when damping is present in the system.A standard approach includes resolving the equations with motion derivatives replaced with Laplacian transform.The issue arises when the final equation obtained with structural number algebra is turned back into its derivative form.When motion derivatives are replaced with the integral solution for sub critically damped motion as in (10), the equation becomes extremely tangled.To overcome this problem, an intermediate method has been developed, which realises the mentioned decomposition of sine and cosine components (10) prior to the application of the structural number algebra, Effectively, that means that the entire graph of relations was duplicated, separating the edges corresponding to sine and cosine coefficients.Each case is then being calculated separately until the last step.Both resulting graphs are shown in Figure 3.
The structural numbers of both graphs are synonymous, the only difference being the opposite sign in the first derivatives of the sine and cosine function, which relate to the damping coefficients.As with the matrix method, the sine graph corresponds to the real part of the calculation and the cosine part to the imaginary part.This time, in order to obtain absolute values of the displacements, the merging (similar to ( 12)) has to be performed for the whole structural number and all the simultaneity functions. 11,29,32The amplitude response of each stage can be calculated by solving the equations: where:

number for the cosine graph,
A ijamplitude calculated with simultaneity function.
The process of obtaining the necessary structural numbers to derive equations (14-16) and the meaning of a simultaneity function has been explained in more detail in. 29,32The computer algorithm which was later used for this purpose is described in detail under the section dedicated to the algorithm testing.The nodes of each graph represent the degrees of freedom of the model they are describing.Each edge corresponds to forces acting between each degree and is being described by the formula of elements creating each force."ω" is the frequency, "m" is the mass of each element, "c" corresponds to spring loads and "b" to damping.

Electromechanical model synthesis
The synthesized mechanical model serves as a basis for the construction of a combined electromechanical model of a piezoelectric system.Constitutive equations and analogy between mechanical and electrical components were used to derive parameters of the piezoelectric model. 14,19,24,25The electromechanical coupling coefficient for piezoelectric coupling perpendicular to the plate surface was determined based on resonant, anti-resonant frequency pairs: where: k 33electromechanical coupling factor in the 3-3 direction, f aanti-resonant frequency of a given stage, f rresonant frequency of a given stage.
A broader term defining electromechanical coupling for a general case, called an effective coupling coefficient was used to approximate the capacity of piezoelectric plates 37 : Another property of a piezoelectric material frequently mentioned in specifications is the mechanical quality factor.This property describes the amount of energy that is being lost during dynamic operation of a piezoelectric plate.It corresponds to mechanical losses in a system.An approximate value of a mechanical quality factor can be calculated using the equation derived from the established Rayleigh's method: Electrical permittivity can be also approximated if piezoelectric capacitance is known by using an equation: where: ε T 33electrical permittivity of a piezoelectric material in 3-3 direction under constant stress, hplate thickness, Aplate cross section field.
Lastly, based on the lumped stiffness of the mechanical model, a piezoelectric material stiffness can be obtained using Hooke's law: where: c E 33piezoelectric stiffness in 3-3 direction under constant electrical field.

Piezoelectric discrete model verification
15]38 The basic form of constitutive equations for piezoelectric effect in a 33 mode was used: where: The equation (24) was modified using Hooke's law and by converting equations for underlying variables: where: U pvoltage across the piezoelectric plate thickness, where: Q pcharge across the piezoelectric plate surface.
By combining the constitutive equations with (25,26) and plugging the equations of motion for the initial mechanical model ( 9) into the mechanical input of the electromechanical system, a lumped piezoelectric model was created.Relations inside a single piezoelectric plate are described by a system of equations: where: Fexternal force applied to a piezoelectric plate, The discrete piezoelectric model has been shown in Figure 4. Inputs and outputs were extracted from ( 27) by grouping parameters and converting the equation into a Laplace form.The resulting matrix was used to create an algorithm for piezoelectric model analysis: The top left coefficient of the coefficient matrix (28) refers to purely mechanical transformations within the piezoelectric plate.Top right and bottom left coefficients refer to the electromechanical coupling.The bottom right coefficient refers to purely electrical transformations.The matrix can be easily expanded by adding more stages to the model as well as it can be calculated using the matrix method mentioned in (10).

Algorithm testing
A series of computer programs were prepared to test the accuracy and computational load of used algorithms.All the programs were written using mathematical software called MATLAB in version R2019b (RRID:SCR_001622).
The code was prepared using basic functions offered by the standard MATLAB suite.
The polynomial equation was created using "for" loops which added consecutive elements to its product form, depending on the number of frequencies given as an input.The MATLAB software operates interchangeably between the product form and exponential form of polynomials, so the conversion was done automatically.The Cauer synthesis was realised with basic program loops which continuously calculated the remainders of the polynomial function, as was explained with equations (3,4,5).Based on the order of elements in the resulting chain fraction, consecutive parameters in form of spring loads and masses were extracted from the equation in each iteration of the loop.
The analysis method using matrix calculations was realised by first creating the equations of motion (9) depending on the number of masses given as an input.The displacement was described using a substitute parameter "X" which then was replaced in bulk with derivatives of the sine cosine integral (10) by using MATLAB function "subs".Next, the equations were transformed into a matrix form and solved with the inbuilt MATLAB solver.The input vector of external forces was used to manipulate the output amplitudes of the matrix formula.
The structural number method was constructed from ground up, using vector transformation operations inside MATLAB.The structural numbers were represented by matrices consisting of simple numbers corresponding to numbered graph edges.Each number was bound to a set of equations representing forces acting inside the system.The method of assigning numbers to equations was explained in. 32The structural number was formed by adding vectors (each vector consists of numbers corresponding to edges connected to a single graph vertex).The process was done using inbuilt MATLAB function called "combvec".After that, several sorting operations were done in a loop to eliminate any repeats in the matrix representing the graph structural number.Further steps done to determine all simultaneity functions were also done by appropriately sorting and pruning the structural number matrix to create each structural number operation.To obtain final equations from the graph algebra, an intermediate step was introduced, which replaced each number with the equation that was bound to it.The step was also done using the "subs" function.
The eletromechanical system analysis was an extended matrix algorithm that was based on the equation (27).No new functions were introduced to make the algorithms and every step was made in analogy to the matrix method used to analyse purely mechanical systems.
The code has been saved in a proprietary MATLAB format ".mlx" and placed in repositories of Mendeley Data site. 42,43he code placed in the repository 42 contains algorithms used for the synthesis of mechanical systems and their comparison through the matrix and structural number methods.The repository 43 contains the algorithm used to synthesise and analyse a full eletromechanical model, simulating the real behaviour of a piezoelectric module.MATLAB is using an open standard for their file formats, so it is possible to open files using standard Windows tools.For reproduction purposes a user has to convert the extension of the file format into ".zip" and open its contents through the ".xml" file stored inside the package.It must be noted however, that the precision and calculation times may vary, depending on the type of software used to reproduce the results and even on the hardware used to process the computation tasks.It is also necessary to note that some of the code may have to be rewritten because of the use of MATLAB specific syntax that may be different in other software.For the sake of consistency and repeatability, all calculations were done on a single computer with a Ryzen 7 2700X processor and a 4 Â 8 GB kit of DDR4 CL16 RAM clocked at 3000 Mhz.
Algorithms used to compare the precision and computation time of both matrix and structural number methods were published in an on-line repository. 42To compare the precision of both algorithms the chosen input data consisted of a collection of resonant (41 kHz and 55 kHz) and anti-resonant (47 kHz) frequencies.The frequencies were chosen arbitrarily based around piezoelectric product specifications regarding different piezoelements in the shape of cubes and cylinders, found on the market.Both algorithms were analysing a model synthesized on the basis of chosen frequencies which do not relate to any specific real piezoelectric element.The necessary system component parameters were calculated using the synthesis algorithm.The damping ratio was left at 100% and the scaling factor H was set to a value of 0.01.The dynamic forces were set to 1N of harmonic force with a cosine wave applied to both degrees of freedom.The program named "comparison" 42 was used to compare the precision of both algorithms and produce the resulting graphs and calculate deviations for the extrema.The results obtained with a matrix method were chosen as the reference point for the results obtained with the non-classical method, using structural numbers.
A simple test was also carried out to check how computationally demanding are both methods in comparison to each other.For that purpose, two programs called "Matrix_method" and "Structural_number_method" were written in MATLAB R2019b.Both algorithms were written with a fixed set of model parameters and a fixed set of equations.
To avoid long computation times, parameters were chosen for a two-degree system resonating at frequencies of 10 rad/s and 22 rad/s with an antiresonant frequency of 17 rad/s.The computation time was measured for both methods using the inbuilt MATLAB timer function "tic/toc".The timer measured seconds that the program took to form the equation and calculate the amplitudes in both cases.Each program was run 10 times with the exact same settings and all calculation time measurements were noted down with an average time derived additionally.Before each run, all the variables generated in the previous runs were cleared to avoid any potential changes in computation speed because of reusing the stored results.
To further approximate the computational load of each algorithm, both were expanded to cover a bigger variety of available systems.Necessary automations were introduced to enable program scalability for systems with larger number of degrees of freedom.Another test was conducted with an aim to determine the calculation time for each algorithm depending on the increasing number of degrees of freedom of the analysed model.The time criterion is insufficient to conclusively determine the computational load of algorithms, however.An additional test was made using the inbuilt Matlab profiler module to verify the program RAM usage during calculations and check how much memory is used by each algorithm during execution.The test was done to roughly estimate the "space requirement" of each algorithm along with the "time requirement", which are used to estimate algorithm complexity. 39Because of increased algorithm size and because of tests being made on models with higher number of degrees of freedom, a decision was made to take five time and memory measurements for each test run.Both algorithms were rewritten and placed in a new online repository. 44r the analysis of the electromechanical system, model with a single degree of freedom was chosen.The resonant frequency was set to 590kHz, based on a thickness vibration mode of a piezoelectric plate listed on the market (https:// www.steminc.com/PZT/EN/piezoelectric-plate-45x45x35mm-55-khz).Much higher resonant frequency was chosen this time, to test model behaviour with higher frequencies.The scaling factor H was set to 0.0025.The damping factor β was calculated from an existing mechanical quality factor Q m of 83 using equation ( 7) which was inversed.The remaining input parameters such as plate dimensions, plate stiffness and capacitance were either taken from an existing piezoelectric material specification or derived using the known physical properties (17-26) (http://www.steminc.com/piezo/PZ_property.asp).The aim of this study was to verify the precision in which a model is capable to simulate the real resonant frequency based on provided material and geometry data.The amplitude response was not a subject of this test, because of insufficient information provided in the piezoelement specification and the inability to conduct real tests at that time.All parameters used in the case of this article are also contained in the program called "analysis_piezoelectric_system" in the on-line repository. 43

Results
A comparison of the efficiency and accuracy of the two methods was done with the comparing program.A two-stage mechanical system with two resonant frequencies of 41kHz and 55kHz (257611rad/s and 345575rad/s respectively) was used for the comparison.The graphs in Figure 5 show the displacement amplitudes for both degrees of freedom of the considered case, obtained by matrix calculation and structural numbers.
Differences were spotted between the results obtained by the two methods.To gain a better understanding, derivatives were taken from the model response equations to calculate the local extremes obtained by both functions.Table 1 gives a direct and relative comparison between the frequencies obtained with both functions.Table 2 gives the calculated values of the displacement amplitude for each characteristic frequency.
The nodal frequencies calculated with the structural method match perfectly with the reference matrix method.There were significant differences in the obtained amplitudes, especially around the second resonant frequency of the system in both stages.The alternative method using structural numbers proved sufficient to check the nodal frequencies.However, when it comes to accuracy in terms of amplitude response, this method was found to be subject to considerable error.Calculations were performed for models with more than two degrees of freedom and the discrepancy between amplitudes of each stage determined with both methods was found to be random.It was expected that the error may result from a mistake performed during the formulation of the structural number algorithm.Another problem with the structural number method was that the authors were unable to establish the equations necessary to calculate phase shifts in the system's response.
Table 3 shows the computation time and an average time from all measurements done on both algorithms based on matrix and structural number methods. 42The test showed that matrix method took an average of 13.972 s to calculate the amplitude response of a two-stage system.The alternative method took only 8.637 s on average, which is 38% faster than the classical method.
The first study was done on algorithms in which the full parameterisation of mathematical relationships has not yet been applied.Both algorithms were rewritten and extended to increase versatility of the program.A second study was carried out, in which the computation time and required memory were measured for algorithms used to calculate systems with different number of degrees of freedom.Each test was run five times, and an average result was taken.Systems with up to six degrees of freedom were tested with each method.Results of the study were placed in Table 4.
The results of the second study indicated a significant increase in the time and required memory for mathematical operations performed by the structural number method relative to the classical matrix method.The difference in calculation time and memory usage of the modified algorithms increased significantly with the number of degrees of freedom of the analyzed systems.The data on memory usage was collected using an inbuilt profiler module inside Matlab which has been shown in Figure 6.An error made in the first iteration of the program was also largely removed during modification of the algorithms and the results obtained with both methods in the second study were matching to the 0.001 μm.An example of results obtained for a system with three degrees of freedom was shown in Figure 7 and placed in Table 5.The comparisons between two algorithms for other studied models are also available in the data repository. 44e discrete electromechanical model 43 was verified on a system simulating the behaviour of a single piezoelectric plate with specifications very close to the real ones taken from (http://www.steminc.com/piezo/PZ_property.asp).The plate had to vibrate at a resonance frequency of 590kHz (3707079rad/s).Mechanical parameters of the plate were determined with the Cauer synthesis method.Piezoelectric plate dimensions were taken from (https://www.steminc.com/PZT/EN/piezoelectric-plate-45x45x35mm-55-khz).Resulting graph of an amplitude response for piezoelectric plate displacement and voltage were shown in Figure 8.
The model had a resonance near the frequency of 560 kHz (3516624 rad/s) which gives an error of 5% in relation to the real plate specification.The calculated magnitude of displacement at the resonant frequency is 2.6 mm which is well above the capability of a 3.5 mm thick piezoelectric plate.The voltage shows a significant increase in piezoelectric voltage in the resonant region and a complete reduction near the frequency of 590 kHz which was supposed to be the original resonant frequency.

Conclusions
The aim of this work was to demonstrate a method that allows a piezoelectric model to be synthesized using only on a set of resonant/anti-resonant frequencies as the input data.The paper has shown that combining Cauer's method of synthesis for mechanical or electrical networks with Rayleigh's damping approximation and a variant of the Van Dyke-Butterworth model can yield a discrete electromechanical cascade model which parameters can be approximated to simulate the behaviour of a real piezoelectric component.This method differs from most other methods found in the literature [13][14][15]19,24 in that it focuses on developing a new set of possible solutions based on a very limited information, rather than trying to break down a specific case under consideration. The utimate goal of this method, which is under constant development, is to provide the ability to synthesize piezoelectric systems capable of resonating in a defined range of frequencies required for a target solution.An example application of such a method could be an attempt to scale up the active vibration cancellation technology to dampen vibrations generated by motors and engines.This technology is currently used in many audio solutions and new, larger piezoelectric resonators could be used for more industrial oriented projects.This method offers the possibility of constructing a piezoelectric stack of plates which have different material properties or geometrical sizes to create multiple resonant frequencies within a single module.
The method presented in this article relies heavily on mechanical model synthesis in its initial stage.In order to proceed with confidence, it was necessary to verify the mechanical model before using it in later stages of the work.However, for more complex systems the classical matrix calculation creates a significant computational load, which has the effect of significantly extending the computation time.An additional set goal was to develop an alternative method for mechanical system analysis using graphs and structural numbers.The presented alternative solution gave promising initial results in terms of computation time.The accuracy of calculated characteristic frequencies between the standard matrix method and the proposed alternative using structural numbers was on par in terms of calculated peaks near resonant frequencies.
Based on the results of the second study conducted on expanded and parametrized algorithms for both methods, the nonclassical method using structural numbers was found much less efficient in terms of both computation time and memory required.It is possible that the outcome of the second study has been an effect of excessive use of "for loops" inside the code.The second study has also proved that results obtained with both methods are comparable in terms of precision.The benefit of classical methods like the matrix method is that they are much better optimized inside calculation software because of their popularity.The structural number method definitely needs further fine-tuning and refactoring to improve its performance but it was shown that it can be an alternative to matrix calculations.The synthesized electromechanical model has shown a small deviation from the target resonant frequency.An initial hypothesis is that this error could be corrected by applying a correction factor to the input frequency.A bigger set of tests on the algorithm has to be conducted to verify it.Another highlighted problem was the calculated amplitude response for the resonant frequency which was several orders of magnitude greater than the expected value of displacement for a plate of given thickness.It is expected that the error originates from the differences between piezoelectric material properties at its resonance and material specification which is measured for a frequency of 1kHz.Authors in Refs.40, 41 mention significant differences in piezoelectric material properties measured at their resonance in comparison to their tion which usually contains properties measured at much lower frequencies.It is clear that presented mathematical model needs more tests before it can be considered acceptably accurate.The reference plate specification (more information here) also didn't contain enough data to reliably compare the results obtained from the model analysis.Only a rough comparison of material properties could be conducted.The data of specific plates such as their exact resonant or antiresonant frequencies was sometimes missing as well as information about their intended mode of vibration (longitudinal, thickness wise) and related coupling factors.A further research will contain a series of empirical tests to obtain reference data more reliably for model validation.

Guangtao Lu
Wuhan University of Science and Technology, Wuhan, Hubei, China Thank the authors for their response and detailed revisions.The current version has been greatly improved, and however there are still some concerns that need be addressed.Moreover, impressed by the innovative effort.
If the name of the proposed method is included in the title, the novelty of this manuscript will be more easily got by the readers. 1.
Figure 3: what do these symbols, 1 X 0s , 1 X 1s and so on, stand for? 2.
The material type and some key parameters of PZT should be clearly addressed instead of just putting a website link.

3.
The resolution of the figures is low and their quality should be further improved.4.
Table 2 and 5: how to confirm the accuracy of the data? 5.
Based on the results, the proposed method can be used to analyze the characteristics of piezoelectric systems accurately.However, the computation cost is much larger than the traditional method.Therefore, what is the merit or advantage of this method?

6.
Competing Interests: No competing interests were disclosed.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Reviewer Report 09 August 2024 https://doi.org/10.5256/f1000research.164274.r307986© 2024 M. Ferreira P.This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Pedro M. Ferreira
Universidade Nova de Lisboa, Lisbon, Lisbon, Portugal The manuscript presents a novel approach to synthesizing piezoelectric systems, starting with only known resonant and anti-resonant frequencies.This method offers a broader perspective on system modelling utilizing an algorithm based on the Cauer synthesis method, an initial discrete mechanical model is constructed, forming the foundation for piezoelectric model synthesis.The manuscript is well-composed and systematically organized.However, there are a few concerns and suggestions that should be addressed before it can be accepted for indexing.In my view, these issues require more detailed clarification before the manuscript can be considered for indexing.
The introduction to the present literature review is poor, particularly in the following excerpt: "Researchers largely focus on modelling simple piezoelectric components with lumped models and matrix tensor equations or try to approach the problem using FEM (finite element method) analysis."There are several models that could be referenced in the introduction, such as Eshelby-Type Models, Micromechanics-based models, or Asymptotic Homogenization-Based Models, which are widely used to model the behaviour of piezoelectric composites.Additionally, this excerpt needs to be supported with more current references.Keywords that could help in finding up-to-date references include: "Piezoelectric composites, Electro-mechanical properties, Analytical models, Numerical models." 1.
In the last paragraph of the introduction, the excerpt "The number of rows and columns, however, is tightly connected to the shape of an edge graph which a given structural number relates to," is very unclear.It should be rewritten for better clarity and explained in more detail.

2.
The author compares the computation time and memory usage of parameterized algorithms, highlighting a significant increase in both when using the structural number method compared to the classical matrix method (as shown in Table 4).The authors should clarify why the calculation time and memory usage of the modified algorithms increase so dramatically as the number of degrees of freedom in the analysed systems rises.

3.
The conclusion is overly long and needs to be more concise, with an emphasis on the limitations of the authors' analytical approach.For better clarity, consider presenting the main conclusions and limitations as bullet points. 4.

Jacek Harazin
Thank You for your insightful review.In line with Your remarks, I have rewritten most of the introductory section of the article, adding more up-to-date literature references and including links to the numerical methods mentioned in the review.I also rewrote the research gap statement to make it more readable.The conclusion section has also been revised.Below I have included corresponding point by point responses to the various comments made in the review.
1.After rereading the introduction section, I must agree that it was too vague and contained too much information that could be considered irrelevant in the scope of the article.To amend that, I have completely rewritten most of this section, while also adding mentions and new references to current works (2020-2023) related with Eshelby-Type/micromechanical models.I must admit that I was so focused on the mathematical synthesis aspect in the context of piezoelectric models that I didn't encounter works on models mentioned in your review.Thank You for bringing this issue up because that also helped me in finding new potential paths to expand my research in the future.3. I agree that this issue was left with a vague explanation.It was already pointed out that the conclusions were overly long and realising that I was trying to make the explanation concise.
Unfortunately, I made it too vague instead.I have tried to amend this by adding more detailed explanation in the conclusions.To summarize the explanation, the uncovered issue may be due to the structure of the written code for the method based on structural number algebra.It was made using closed loops which were combined with quite complicated substitution operations and complex algebraic representations of each graph edge.Combined, these three elements could contribute to the perceived great increase in both computation time and RAM requirement.The matrix method was largely an unmodified module readily available as an integral part of the Matlab software which may have been the reason for its much better performance.The new method had to be written from ground up which probably left a lot of the performance on the table due to my lack of programming skill and also the inability to better access program's memory and CPU handling operations.
4. Most of the conclusions section was reworked according to raised issues.I have tried to form several bullet points to streamline the thought process and make the text more accessible to readers.I wasn't able to put the entire section in a set of bullet points, but I made an attempt to better separate the conclusions from both experiments described in the article to avoid confusion.I hope that will satisfy the requirements.
The manuscript proposes a new approach for parameter estimation of piezoelectric structures based on the Cauer's method, and compares its computational performance with the traditional matrix method.The authors claim that their method can synthesize a piezoelectric model using only a set of resonant and anti-resonant frequencies as the input data, and approximate the piezoelectric parameters using the synthesized mechanical model and some physical equations.
The authors also use a non-classical method based on structural numbers and edge graphs to analyze the mechanical model and reduce the computational load.The manuscript is well-written and organized, and the mathematical derivations are clear and rigorous.However, I have some concerns and suggestions that need to be addressed before the manuscript can be accepted for publication.
The main merit and novelty of the manuscript are not clearly stated in the abstract, I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 25 Mar 2024

Jacek Harazin
Thank you for the critical review.I would also like to apologise for a long break before my reply.I have prepared several corrections to the article regarding your review and they will be published with the new version.I've also tried to address your concerns directly: 1. Let me repeat my answer to the second reviewer here, as it addresses the same concerns.
The article showcases a new method to synthesize cascading electromechanical models using the Cauer method coupled with Van Dyke-Butterworth models to represent piezoelectric stacks with clumped mechanical and electrical components.Each of those methods have been already used by researchers but i have found no literature on the particular approach presented in this article.The novelty revolves around the general approach to synthesis with very limited input information.Most researchers focus on specific applications where piezoelectric parameters and the analysed structure is largely known or predetermined experimentally.Here, only the frequencies are known and there is an infinite number of possible models that could realise those predetermined resonant frequencies.Coupled with further steps presented in the article, this method offers a possibility to design piezoelectric stacks with different material properties or geometric dimensions, capable of resonating in predetermined frequencies.The limitation of this method is that it can only be used to synthesize models composed of components connected in series.So far, it can only be used to synthesise specific piezoelectric plate configuraions, where each component in a stack is responsible exclusively for a single resonant frequency set as the input.I have modified the article to better reflect those statements.
2. Following up from my response to the second reviever, I agree and i can already attest that the method is currently being verified empirically.The synthesis method has been verified already with three real piezoelectric stacks with different geometrical dimensions which has shown some discrepancy between model results and real plate behaviour.The source of those differences was traced back to how some of the piezoelectric parameters are determined in the factory (under 1kHz wave stimulation) versus plate properties in the resonance proximity.It has to be noted as well, that this method relies on a lumped representation of an otherwise continous system.The main issue of this approach is that the lumped parameters are approximations of the real parameters and an electromechanical model currently can't perfectly match the behavior of a piezoelectric plate, where the interactions and their magnitude changes between the layers of a single piezoelement.There are also boundaries between piezoelectric elements, made with adhesives and isolating layers, which also introduce many hidden variables in the real model behaviour.The presented method is still in its infancy and requires more testing and refinement.I would like to provide a consice report on the empirical results in a follow up article though, as including all the information about the procedure and test stand would make this article very bloated.
3. I was made aware that a comparison of calculation time used exclusively is a bad metric for algorithm complexity.I have tried looking for more literature on the matter and added more references to the article.Given my abilities and available resources, i have prepared a second study which focused more on the memory criterion.I've used the available Matlab evaluation modules to estimate the average RAM usage between both algorithms on a variety of models, along with a comparison between the average time it took both algorithms to complete.New suite of tests has uncovered more information about the structural number algorithm optimisation and shed some light on its shortcomings.I have modified the manusript to showcase the second study.I am aware that this method may also not be the best way to evaluate algorithm complexity and presents a rough estimation and i will try to consult a specialist on the topic with an intent of providing some collaboration with more indepth information about algorithm complexity of the proposed method.
4. I have provided additional literature references to publications of the authors involved in the development of the structural number method where they descibe mentioned models in english.I hope that will provide the necessary background which was lacking previously.
To compare the computational performance of the new model, the model response and computation time are compared with the matrix method.What are the errors of these two models?The two frequencies of the new model are compared with the ones of matrix method.Are the results of the matrix method accurate?What are the actual frequencies?
It is the same during the amplitude comparison.Moreover, to compare the computational performance of the algorithm, there are many metrics, such as time cost, the number of floating points (FLOPs) multiple (MAdds), and so on.In addition, what is the meaning of only decreasing the computation time?

3.
Experimental results: the results should be further analyzed.4.

Are sufficient details of methods and analysis provided to allow replication by others? Partly
If applicable, is the statistical analysis and its interpretation appropriate?Not applicable Are all the source data underlying the results available to ensure full reproducibility?Yes

Are the conclusions drawn adequately supported by the results? Yes
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Structural health monitoring by uisng lamb waves I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
clarity.I would also like to apologise for a long break before my reply.I have prepared several corrections to the article regarding your review and they will be published with the new version.I've also tried to address your querries to the best of my ability: 1.The article showcases a new method to synthesize cascading electromechanical models using the Cauer method coupled with Van Dyke-Butterworth models to represent piezoelectric stacks with clumped mechanical and electrical components.Each of those methods have been already used by researchers but i have found no literature on the particular approach presented in this article.The novelty revolves around the general approach to synthesis with very limited input information.Most researchers focus on specific applications where piezoelectric parameters and the analysed structure is largely known or predetermined experimentally.Here, only the frequencies are known and there is an infinite number of possible models that could realise those predetermined resonant frequencies.Coupled with further steps presented in the article, this method offers a possibility to design piezoelectric stacks with different material properties or geometric dimensions, capable of resonating in predetermined frequencies.
2. The predetermined frequencies and amplitudes were chosen arbitrarily and do not reflect any real piezoelectric system.A comparison between both analysis methods was done on the same model synthesized with the Cauer method.The matrix method was considered as a reference method because its a known and established way to calculate model response.The structural number method was compared against it and it was a pure comparison between two analytical methods.I was also made aware that a comparison of calculation time used exclusively is a bad metric for algorithm complexity.I have tried looking for more literature on the matter and added more references to the article.Given my abilities and available resources, i have prepared a second study which focused more on the memory criterion.I've used the available Matlab evaluation modules to estimate the average RAM usage between both algorithms on a variety of models, along with a comparison between the average time it took both algorithms to complete.I am aware that this method may also not be the best way to evaluate algorithm complexity and only presents a rough estimation.I will need to consult a specialist on that matter and conduct more tests but that would probably go beyond the scope of this article and should be contained in a separate paper focused only on the comparison between both methods.
3. I have provided additional literature references to publications of the authors involved in the development of the structural number method where they descibe mentioned models in english.I hope that will provide the necessary background which was lacking previously.4. I agree and i can already attest that the method is currently being verified empirically.The synthesis method has been verified already with three real piezoelectric stacks with different geometrical dimensions which has shown some discrepancy between model results and real plate behaviour.The source of those differences was traced back to how some of the piezoelectric parameters are determined in the factory (under 1kHz wave stimulation) versus plate properties in the resonance proximity.However, i would like to envelope the empirical analysis, along with all the necessary test procedures and the test stand description in a second article, serving as the continuation for this study.

Figure 3 .
Figure 3. Graph of the cascade model where: a) shows the sine coefficients b) shows the cosine coefficients.The nodes of each graph represent the degrees of freedom of the model they are describing.Each edge corresponds to forces acting between each degree and is being described by the formula of elements creating each force."ω" is the frequency, "m" is the mass of each element, "c" corresponds to spring loads and "b" to damping.
: k effeffective electromechanical coupling factor, C pplate capacitance, C mequivalent capacitance from mechanical stiffness, c mmechanical stiffness of a lumped piezoelectric model.

d 33 -
piezoelectric charge constant, s E 33mechanical compliance constant (inverse of mechanical stiffness constant).

Figure 4 .
Figure 4. Discrete piezoelectric model composed of mechanical lumped model combined by a transformer with an electrical network."F"represents an external force, "c" is the mechanical stiffness, "b" corresponds to damping, "N" is the electromechanical coupling, "C" is the electrical capacitance of the piezoelement and "U" represents the voltage exciting the circuit."a,b,h" correspond to geometric dimensions of the piezoelectric element, while "x" represents the displacement.

Figure 5 .
Figure 5. Model response calculated for a) first vibrating mass, and b) second vibrating mass.Blue line represents a response calculated with matrix method; red line represents the response calculated with structural number algebra.

Figure 6 .
Figure 6.An inbuilt Matlab profiler module used to collect the data on the algorithm memory usage.

Figure 7 .
Figure 7.A graph showing model response obtained with analysis done using both methods.Green line represents model response calculated with the matrix algorithm.Red dots show datapoints calculated with structural numbers.The obtained results show a perfect match.

Figure 8 .
Figure 8. Model response calculated for a) displacement curve, and b) voltage curve along the frequency axis.

2 .
I have decided to extend the explanation with a graphical clue in form of a new figure 1, which shows an example of a graph, from which a new structural number is being created.The example of a structural number has also been placed under new equation 1.I think that a simple figure might provide much more insight than my extensive explanations, while also conserving more space.In hindsight, I should have provided a descriptive figure in the first place.

2 .
Page 4: "A combination of those factors creates four different types of models…" It refers to Refs. 8 and 11.However, Refs.8 and 11 are published in non-English.Therefore, the details of the four models are suggested to be introduced in the manuscript.

Table 1 .
Comparison of frequencies between model responses obtained with both methods.

Table 2 .
Comparison of amplitude responses obtained with both methods on the nodes.

Table 3 .
Computation time for both methods used to analyse a simple mechanical system.

Table 4 .
Computation time and memory usage for parametrized algorithms.

Table 5 .
Results obtained for both methods in the resonance frequency data points.

the work clearly and accurately presented and does it cite the current literature? Partly Is the study design appropriate and is the work technically sound? Yes Are sufficient details of methods and analysis provided to allow replication by others? Yes If applicable, is the statistical analysis and its interpretation appropriate? Not applicable Are all the source data underlying the results available to ensure full reproducibility? Yes Are the conclusions drawn adequately supported by the results? Partly Competing Interests:
No competing interests were disclosed.