Keywords
piezoelectricity, mathematical model, synthesis
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.
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.
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.
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.
piezoelectricity, mathematical model, synthesis
See the authors' detailed response to the review by Pedro M. Ferreira
See the authors' detailed response to the review by Saurabh K. Yadav
See the authors' detailed response to the review by Guangtao Lu
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 Systems that precisely control room temperature or humidity, systems responsible for engine control and its smooth acceleration or deceleration or systems used for active vibration cancellation are only few of the examples.3–6 Polynomial equations are great mathematical tools for approximation and interpolation of numerous phenomena which is why they are commonly used in mathematical modelling.7–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 theory2,7–9,11–16 and were used during the development of signal filters composed of RLC (R – resistor, L – coil, C – capacitor) components.17,18
Wilhelm 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 synthesis17,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. Moreover, because of existing analogies between electrical and mechanical elements, this method can be applied to both electrical and mechanical models.7–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 piezoelectricity12 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 models13,14 and matrix tensor equations20 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 of 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–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 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.28 The 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,31 The 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 in15,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 book29 and in articles.30,31 More details can also be found in previous article published by authors.32 The article also directly relates to adaptations of this method used by other researchers for analysis and synthesis of various systems.2,11,16
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:
T(s) – response of a system,
H – response 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:
ai – numerator coefficients,
bi – denominator coefficients.
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.17 The 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:
– fraction extracted by division of parameters next to the highest power “s”,
– rest 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,35 However, 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 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:
– stiffness damping coefficient,
– damping value,
– spring stiffness.
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,25 The 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.
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,31 The use of a second method is considered as part of the ongoing research32 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:
Under the assumption into account, that the system is only considered in the case of subcritical damping, an integral of a second-order differential equation36 for the variable displacement x can be described by the equation:
– real part of the ith element amplitude response,
– imaginary part of the ith 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,16 The 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 method29 or related articles30,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 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.
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,32 The amplitude response of each stage can be calculated by solving the equations:
– structural number for the sine graph,
– structural number for the cosine graph,
– amplitude 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,32 The computer algorithm which was later used for this purpose is described in detail under the section dedicated to the algorithm testing.
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,25 The electromechanical coupling coefficient for piezoelectric coupling perpendicular to the plate surface was determined based on resonant, anti-resonant frequency pairs:
– electromechanical coupling factor in the 3-3 direction,
– anti-resonant frequency of a given stage,
– resonant 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 plates37:
– effective electromechanical coupling factor,
– plate capacitance,
– equivalent capacitance from mechanical stiffness,
– mechanical stiffness of a lumped piezoelectric model.
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:
– electrical permittivity of a piezoelectric material in 3-3 direction under constant stress,
– plate thickness,
– plate cross section field.
Lastly, based on the lumped stiffness of the mechanical model, a piezoelectric material stiffness can be obtained using Hooke’s law:
– piezoelectric stiffness in 3-3 direction under constant electrical field.
To analyse piezoelectric response with calculated parameters, a discrete electromechanical model was created based on the Butterworth-Van Dyke model and other lumped piezoelectric models available in the literature.13–15,38 The basic form of constitutive equations for piezoelectric effect in a 33 mode was used:
– electric displacement field,
– mechanical strain,
– electric field,
– mechanical stress,
– piezoelectric charge constant,
– mechanical compliance constant (inverse of mechanical stiffness constant).
The equation (24) was modified using Hooke’s law and by converting equations for underlying variables:
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:
– external force applied to a piezoelectric plate,
– external voltage applied to a piezoelectric plate,
– lumped piezoelectric mass,
– lumped piezoelectric damping,
– lumped piezoelectric stiffness,
– lumped piezoelectric capacitance,
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:
“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.
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).
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.32 The 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,43 The code placed in the repository42 contains algorithms used for the synthesis of mechanical systems and their comparison through the matrix and structural number methods. The repository43 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.42 To 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.39 Because 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.44
For 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 Qm 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
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.
Blue line represents a response calculated with matrix method; red line represents the response calculated with structural number algebra.
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.
Matrix method | Structural numbers | Direct difference | Indirect difference | |
---|---|---|---|---|
1st stage | 0.0101 mm | 0.0101 mm | 0 mm | 0% |
0.0017 mm | 0.0027 mm | 0.001 mm | 58.82% | |
2nd stage | 0.0033 mm | 0.0033 mm | 0 mm | 0% |
0.0004 mm | 0.0006 mm | 0.0002 mm | 50% |
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.42 The 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.44
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.
Frequency [Hz] | Amplitude [μm] | |
---|---|---|
Matrix method | Structural method | |
45000 | 22.0713 | 22.0717 |
75000 | 6.9048 | 6.9048 |
93000 | 1.6714 | 1.6714 |
The discrete electromechanical model43 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.
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 literature13–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 ultimate 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 non-classical 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 specification 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 anti-resonant 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.
Reference data for the specific piezoelectric module has been sourced from a piezoelectric materials vendor: https://www.steminc.com/PZT/EN/piezoelectric-plate-45x45x35mm-55-khz.
A list of detailed material properties used for algorithms and the electromechanical model analysis was taken from: http://www.steminc.com/piezo/PZ_property.asp.
Mendeley Data: Model algorithms. https://doi.org/10.17632/gdt8yr9m3d.4 42
This project contains the following data:
Mendeley Data: Piezoelectric electromechanical model (first iteration), https://doi.org/10.17632/pydnfmkfdk.3 43
This project contains the following underlying data:
Data are available under the terms of the Creative Commons Zero “No rights reserved” data waiver (CC0 1.0 Public domain dedication).
Source code available from:
Archived source code available from: https://doi.org/10.17632/gdt8yr9m3d.4 42
License: CC0 1.0
Piezoelectric electromechanical model
Source code available from:
Archived source code available from: https://doi.org/10.17632/pydnfmkfdk.3 43
License: CC0 1.0
Improved analysis algorithms for structural and matrix methods
Source code available from:
Archive source code available from: https://doi.org/10.17632/xk443p7m7d.1 44
License: CC0 1.0
Views | Downloads | |
---|---|---|
F1000Research | - | - |
PubMed Central
Data from PMC are received and updated monthly.
|
- | - |
Is the work clearly and accurately presented and does it cite the current literature?
Yes
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?
Yes
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Maintenance, Piezoelectricity, Sensors
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Structural health monitoring by uisng lamb waves
Is 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.
Reviewer Expertise: Piezoelectric Composites, Piezoelectric Behaviour, Advanced Materials, Smart Materials
Is the work clearly and accurately presented and does it cite the current literature?
Yes
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?
Yes
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: Mathematical modelling
Is the work clearly and accurately presented and does it cite the current literature?
Yes
Is the study design appropriate and is the work technically sound?
Partly
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
Alongside their report, reviewers assign a status to the article:
Invited Reviewers | ||||
---|---|---|---|---|
1 | 2 | 3 | 4 | |
Version 3 (revision) 06 Sep 24 |
read | read | ||
Version 2 (revision) 02 Apr 24 |
read | read | read | |
Version 1 05 Oct 23 |
read | read |
Provide sufficient details of any financial or non-financial competing interests to enable users to assess whether your comments might lead a reasonable person to question your impartiality. Consider the following examples, but note that this is not an exhaustive list:
Sign up for content alerts and receive a weekly or monthly email with all newly published articles
Already registered? Sign in
The email address should be the one you originally registered with F1000.
You registered with F1000 via Google, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Google account password, please click here.
You registered with F1000 via Facebook, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Facebook account password, please click here.
If your email address is registered with us, we will email you instructions to reset your password.
If you think you should have received this email but it has not arrived, please check your spam filters and/or contact for further assistance.
Comments on this article Comments (0)