Coupled electro-aeroelastic energy harvester model based on piezoelectric transducers, VIV-galloping interaction and nonlinear switching circuits

Miniaturization, multi-functionality, and low power consumption are key features in the development of electronic devices today. These characteristics open up the possibility for the development of small, self-powered systems. Piezoelectric materials have emerged as a promising solution due to their ability to convert mechanical energy into electrical energy from the surrounding environment, such as wind. In this research work, we investigate an aeroelastic piezoelectric energy harvester based on the interaction between vortex-induced vibrations and galloping. The device, a bimorph cantilevered beam with a square-section bluff body attached to its tip, produces mechanical oscillations that are extracted by various electrical circuits, including a nonlinear switching interface. A numerical model, implemented in a Simulink environment, is developed to solve the coupled electro-aeroelastic problem for each electrical interface. The results are compared to experimental data collected from wind tunnel tests conducted at various wind velocities and load resistances, and a satisfactory correlation is observed. The validated model is finally used to carry out a parametric analysis aimed at optimizing the performance of the harvester in a low wind speed range. This investigation could contribute to a better understanding of an aeroelastic harvesting system and potentially assist in the optimization of design parameters and electrical interfaces for specific applications.


Introduction
Power harvesting allows the conversion of mechanical energy from the environment, which would be otherwise wasted, into usable electric energy [1]. Piezoelectric transduction is one of * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the most promising harvesting mechanisms because of its high electromechanical coupling factor compared to electrostatic, electromagnetic, and triboelectric conversion [2].
In the fields of civil and aerospace engineering, vibrations of significant amplitude may occur when a structure is subjected to fluid dynamic loads. Attempts are made to avoid these kinds of vibrations in large-scale systems such as buildings, bridges, and aircrafts to prevent possible damage or failures. In recent years, these aeroelastic phenomena have been proposed as sources of energy for small-scale systems [3]. Flutter, galloping, and vortex-induced vibrations (VIV) are the instabilities typically exploited for this purpose [4]. Since rotational modes' oscillations are not converted into electrical energy in the case of piezoelectric transduction, it is convenient to focus on galloping and VIV, which feature a single translational degree of freedom in the transverse direction to the incoming flow [5]. The main drawback of VIV is having a limited speed range in which there is a sufficient power production, the socalled synchronization region [6]. Galloping-based devices on the other hand, are not constrained to operate in narrow velocity ranges, but they can have a high instability onset speed, especially in fairly damped systems.
In the cases where the velocity associated with vortex resonance U v and the velocity associated with galloping U g are not sufficiently different from each other, an interaction between the two aeroelastic phenomena may arise [7]. This interference exhibits a velocity-unrestricted behavior like that of galloping, but it begins at the typical velocity of VIV. Since typically U v < U g , it may be advantageous to realize a harvesting device that exploits the above-mentioned interaction, allowing for energy extraction starting from very low speeds. This condition is certainly desirable for powering electronic devices in an almost continuous manner.
One of the first studies to propose utilizing the interaction between two aeroelastic phenomena for power generation was conducted by Andrianne et al [8]. This study involved experimental and numerical analyses on a square-section cylinder that used a coil-magnet as a transducer. He et al [9] demonstrated one of the earliest applications of this interaction in the context of piezoelectric harvesting. They conducted experimental investigations comparing generators that kept VIV and galloping decoupled with those in which the interference occurred. Such comparisons revealed that the energy output of the harvesting devices with the interaction was preferable because they produced an adequate energy output starting from low velocities. In a successive study by Yang et al [10] a piezo-aeroelastic mathematical model capable of predicting the behavior of both devices with and without interaction was developed. Recently, Wang et al [11] investigated various cross-sectioned bluff bodies, given as a combination of circular-and cuboid shapes to enhance the performance of a scavenger that simultaneously activated both VIV and galloping.
In the aforementioned studies, the electrical interface considered is a simple resistor. However, in order to practically implement piezoelectric harvesting for powering electronic devices, it is necessary to examine more complex circuits operating in DC. Thus, the objective of this paper is to extend previous research and investigate typical energy harvesting interfaces such as the standard DC interface and the switch synchronized harvesting on inductor (SSHI). To achieve this goal, an efficient numerical model of the coupled electric and aeroelastic domains is required. However, deriving the governing equations for complex nonlinear circuits is challenging. Although one approach to address this issue is to convert the mechanical properties of the harvester (such as mass, stiffness, damping, etc) into their electrical equivalents and use a circuit simulation software, it should be noted that this conversion is not straightforward in cases where aeroelastic phenomena are present. For this reason, we have developed a numerical model in a Simulink environment, which allows us to analyze the entire electromechanical system without any need for conversion. This investigation can improve the understanding of the harvesting system and aid in optimizing design parameters and the electrical interface for specific applications.

Mathematical model
The energy harvesting device considered in the present study is a bimorph cantilevered beam with a square prism attached to its free end. By changing the mass of the bluff body, the characteristic velocities U v and U g can be modified, enabling the generator to combine the effects of VIV and galloping. Initially, a simple resistor is used as an electrical interface, which is later replaced with more complex circuits in subsequent sections. We utilize the model of Tamura and Shimada [12] to describe the interaction between the aeroelastic phenomena, following the dissertation of [10] and generalizing it to allow its application to the bimorph prototype shown in figure 1.
The piezoelectric beam equation is based on the Euler-Bernoulli beam theory, which was chosen over Timoshenko's model due to the beam's small cross-sectional dimensions relative to its length squared: as a result, the contributions of transverse shear strain and rotary inertia can be considered negligible [13]. The governing equation is then: where m is the mass per unit length of the cantilever, c a is the viscous damping coefficient and c s is the equivalent coefficient of strain rate damping; YI is the flexural stiffness of the beam and θ is the electromechanical coupling coefficient, which are expressed as: Regarding the force and moment at the free end, the inertial contributions related to M t and I t (mass and moment of inertia of the prism) can be separated from the aerodynamic ones: (2.5) The deflection w can be expressed by means of a convergent set of eigenfunctions and make use of only the first term of the summation, since it has been shown that VIV and galloping induce oscillations around the first natural frequency of the system [14,15]: w (x, t) = ϕ 1 (x) η 1 (t). By replacing this relationship in (2.1), multiplying everything by ϕ 1 (x) and integrating over the length of the beam it is possible to derive the following equation: where ω 1 and ζ 1 are the natural pulsation and mechanical damping associated with the first eigenmode. The eigenfunction is normalized so that ϕ 1 (L b ) + 0.5ϕ ′ 1 (L b ) = 1. Consequently, the equivalent mass, piezoelectric coupling coefficient, and aerodynamic modal force are given by: The aerodynamic force is expressed by means of the model of Tamura and Shimada, which combines the quasi-stationary model of galloping with the wake-oscillator for the VIV: (2.10) The angular displacement of the wake-oscillator α instead satisfies the equation where γ and ω VIV are the natural damping and pulsation of the wake oscillator, and C L0 the rms amplitude of the oscillating lift coefficient.
To complete the set of governing equations, it is necessary to introduce an additional equation, associated with the electrical domain (h pc is the distance between the center of the piezoelectric layer and the neutral axis and ε ε 33 is the permittivity component at a constant strain), which decomposed by means of the eigenfunctions method becomes In summary, the equations of the model are: (2.14)

Numerical model
The previous section discussed a standard AC interface circuit, where a resistive load is connected in parallel to the piezoelectric generator. This simple circuit is often used in literature for  an initial assessment of the electromechanical behavior [16]. However, generally electronic devices operate on DC [17], so it is necessary to investigate other electrical interfaces. Therefore, a standard DC interface (figure 2) is considered, which includes a diode bridge rectifier, a filter capacitor, and a resistive load. The diode bridge rectifier ensures that a DC current flows through the load, while the filter capacitor smooths the voltage across the resistor.
In a standard DC interface it cannot be guaranteed that energy is continuously transferred from the mechanical to the electrical domain: during a certain amount of time in each cycle of oscillation there is a return of energy [18]. This problem is solved by using nonlinear harvesting interfaces, which modify the waveform of the piezoelectric voltage to have a one-way flow of energy. An example of nonlinear techniques is the parallel switch synchronized harvesting on inductor, or P-SSHI, which involves inserting an inductance driven by a switch in parallel to the piezoelectric element [19], as can be seen in figure 3.
The switch stays open for most of the time, and when the displacement of the piezoelectric reaches its maximum or minimum values in the cycle, the switch closes, establishing an oscillating electrical network composed of the piezoelectric capacitor and the inductance. The switch remains closed for half the period of oscillation of the inductor-capacitor (LC) circuit. During this interval, the electrical energy stored on the electrodes is transferred to the filter capacitor via the inductance, and the voltage across the piezoelectric is almost instantaneously reversed [20].
To practically implement a P-SSHI interface (figure 4), we make use of the electronic breaker proposed by Richard [21], which automatically switches when the voltage across the piezoelectric elements reaches its maximum/minimum values. In this way there is no need for a displacement sensor and a controller for the synchronization and generation of the switch signals: these devices may have a power consumption superior to the extracted power, thus the advantage of realizing a self-powered interface [22].
To evaluate the performance of the developed energy harvester, it is necessary to solve its governing equations. These are a set of nonlinear ordinary differential equations (ODEs), the solution of which can be found numerically via numerous solving algorithms, for example Matlab's ODE45 function. The use of the latter, however, implies the possibility to write explicitly the equations associated with the problem, which is certainly possible when the interface circuit is made of a simple resistor, but becomes much more difficult when circuit components of greater complexity are involved, such as diodes and transistors. For example, a P-SSHI circuit includes 1 current generator, 4 capacitors, 10 diodes, 4 transistors, 1 inductance and 6 resistors. The equations associated with this interface are numerous and non-trivial, because of the multitude of components and the non-linearity of the characteristic relations of diodes and transistors. Any change in the circuit configuration would require re-deriving the entire set of equations associated with the electrical part. To overcome these issues, we use Simulink's Simscape package, which enables the creation of models of physical components based on physical  links that integrate directly with block diagrams and other modeling paradigms. Figure 5 shows the Simulink model for the standard AC interface, which comprises a mechanical and an electrical subsystem. The mechanical subsystem is comprised of a block diagram that represents the first two equations of (2.14) in state space form. As for the circuit interface, as previously stated, the Simscape package is used, and consequently we can easily draw the circuit diagram. Within the electric circuit there is a current generator controlled by an external variable. Its output current is the numerical value presented at the input port of this circuit component, −χη 1 . In addition, there is a voltage sensor in parallel to the piezoelectric capacitor and load resistance which allows to generate the signal corresponding to the piezo-voltage to be provided as input to the equations of the mechanical domain.
To implement different electrical interfaces, it is sufficient to change the circuit configuration in the electrical domain, while the block diagram of the mechanical domain remains the same. Figure 6 shows the standard DC interface and the self-powered version of the P-SSHI. As can be seen, there are two voltage sensors in the circuits. The one in parallel with the capacitor and the current generator of the piezoelectric generates the V signal that is given as input to the mechanical subsystem, while the one in parallel to the resistive load allows to evaluate the generated power.
An advantage of this kind of modeling is that the blocks associated to the different circuit elements allow to include the non-idealities, and to calibrate the parameters associated to the characteristic equations of the components according to the datasheet of the correspondents used in the practical realization, for example h FE for transistors or the voltage drop across a diode. The adopted solving strategy allows the simulation of the whole electromechanical problem within a single environment. To the authors' knowledge the only alternative found in the literature in the case of complex circuits, consists in transposing the mechanical properties (mass, stiffness, damping, etc) to electrical quantities and use a circuit simulation software, like in [23]. However, when aerodynamics is present it is not always trivial to operate in this way (for example, it has not been found any electric correspondence for the VIV-galloping interaction in literature) and therefore it appears evident the convenience of using Simulink with Simscape.
Moreover, it is possible to notice how the methodology presented here is generalizable to harvesting problems of different nature (e.g. energy production from base vibrations, flutter, etc) provided that the mechanical equations can be written as a set of ODEs.

Model verification
To create the piezoelectric generator in practice, we manufactured a fiberglass reinforced polymer beam with two DuraAct A12 patches attached in a bimorph configuration near the fixed end. The use of a composite material allows to avoid plasticization phenomena that may arise in the proximity of the fixed end because of the high mechanical stresses. The prismatic body was manufactured in polystyrene foam, and its mass is such to allow the VIV-galloping interaction.
The first natural frequency of the system was measured from the frequency response functions obtained by exciting the piezoelectric generator with an electromagnetic shaker. Testlab's Spectral Testing environment was used to obtain the frequency response function (FRF), in conjunction with a LMS Scadas III , which handles both the input signal conditioning and the data acquisition ( figure 7). The mechanical damping coefficient was instead obtained using the logarithmic decrement method.
The constraint has been implemented by a steel beam fixed at both ends, such as to be rigid enough not to deform during the motion of the harvesting device, but at the same time sufficiently thin to affect the aerodynamic field in a limited way. To ensure the first condition, we carried out a finite element analysis on the whole experimental system to verify that the steel beam does not participate significantly to the first mode shape. The finite element method (FEM) model was created with MSC Marc and it is shown in figure 8 on the left. As can be seen, the model is three-dimensional and makes use of 3D Solid 8-nodes type elements. The inertial properties of the prismatic body were enclosed in a single node located at its center of gravity, which was connected via RBE2s (rigid body The electrical circuits were realized on a breadboard: a coaxial cable carries the piezoelectric terminals to the breadboard, where they are combined in pairs to achieve a parallel configuration and then connected to the different circuit elements of the various interfaces.
A Keysight 34470A multimeter with data-logging functionalities is used for voltage measurements, with a sampling time of 299 µs. To avoid any kind of aerodynamic interference we make use of a Fastcam Mini AX100 high framerate camera for displacement measurements, which allows the acquisition of a large number of images per cycle: the camera  was set to record with a frame rate of 750 fps, a shutter speed of 1/5000 s and a resolution of 1024 × 1024. The whole experimental setup in a wind tunnel of Sapienza University of Rome is shown in figure 9.
Since the camera lens was not perfectly aligned with the piezo beam (figure 10) on the vertical axis it has been necessary to correct the displacement measurements to eliminate the parallax error.
The characteristics of the piezoelectric beam are shown in table 1.
Two types of measurements were made in the wind tunnel tests: 1. Measurements at varying speed with a fixed value of resistance. 2. Measurement at varying resistance with a fixed value of the flow velocity.
In the numerical simulations where velocity is varied, the considered speed range starts at U = 1.6 m s −1 , at which an initial condition η 10 = 0.01B is imposed on the displacement. In order to simulate the continuous velocity growth observed in the experiments, we assumed that the initial displacement at each successive velocity was equal to the steady-state displacement evaluated at the previous velocity. On the other hand, in simulations with different resistance values, we trace the bifurcation curve up to the desired velocity to determine the initial displacement condition to be applied. Only then do we proceed to vary the resistance value, R.
The Tamura and Shimada model requires the choice of some parameters, which depend on the section considered. The values used in the present numerical simulations were found from sources in the literature and are listed in table 2.
In their study Tamura and Shimada used a value f m = 1.16, which is the one of a circular cylinder. Since in [10] it was found that f m = 3.48 provided more accurate results for a square cylinder, we decided to apply this value in the following analysis.
The numerical solver used in the simulations is daessc, a variable step solver designed specifically for Simscape which provides robust algorithms specifically designed to simulate differential-algebraic equations arising from modeling physical systems (in our case the electrical circuits).
Since to evaluate the mass M 1 and the electro-mechanical coupling coefficient χ it is necessary to know the eigenfunction ϕ 1 , a code capable of solving the eigenvalue problem of the piezoelectric generator was developed in Wolfram  Mathematica, because of its symbolic calculation capabilities which allowed to obtain an analytic expression of ϕ 1 over the whole domain.

Standard AC
In figure 11(a) the trend of the root mean square voltage across the load as the speed varies is shown. The curve is a spline interpolation of the values provided by numerical simulations, while the points are the experimental data. The onset speed of the instability is very well predicted, at around 2.4 m s −1 . In the 8 m s −1 speed range, the numerical and experimental values are very close to each other and almost coincident at some points. The error is maximum between 4 and 5 m s −1 , where it reaches 63% in U = 4 m s −1 , and it is possible that this is due to an incorrect estimation of the aerodynamic parameters related to VIV, in particular f m . This term has a significant impact on the behavior of the aeroelastic system as it was pointed out by Mannini et al [24]. In figure 11(b) the behavior of the maximum displacement of the free end of the beam against the flow speed is shown. As it can be seen, there is still a good agreement between experimental results (dots) and numerical (curve), but the slopes of the two are slightly different.
The generated power was chosen as the comparison parameter when resistance varies, since it is expected to exhibit a maximum and it is desired to verify that numerical simulations can adequately predict the optimal value of R, as well as the experimentally recorded levels of P.
As it can be seen in figure 11(c), there is an excellent agreement between the results, both qualitatively and quantitatively.
The predicted values are very close to the experimental ones, substantially coinciding everywhere except around the maximum, where however the error stays small and lower than 10.3%. The experimentally measured optimal resistance is 68 kΩ, while the numerically predicted one is 61 kΩ.

Effect of higher order modes.
In section 2 we reduced the continuous model to a 1 degree of freedom (DOF) lumped system by considering only the first mode shape, since both galloping and VIV induce oscillations around the first natural frequency of the system. In the present section we take into account the effect of the higher mode shapes to evaluate the error introduced by the above assumption. Thus, we compare the root mean square voltage curves obtained for the 1DOF and 3DOF systems, where the latter makes use of the first three bending modes of the beam. The governing equations for the 3DOF system are [10]: M rηr + 2M r ζ r ω rηr + M r ω 2 r η r + χ r V = f a r = 1, 2, 3 (4.1) From figure 12 it can be seen that the voltage in the case of the 3 degree of freedom system is slightly higher at all wind speeds. Overall, the trend of the curve remains unchanged and the differences in voltage always remain less than 0.9 V. The 1DOF model is advantageous for its simplicity and it is more compact to represent in a state-space form, without introducing particularly significative errors.

Standard DC
In figure 13(a) it is possible to observe the behavior of the load voltage as the flow speed varies. In this case, since the voltage is no longer sinusoidal but stays approximately constant, it no longer makes sense to consider its root mean square. However,  there are still slight oscillations due to the voltage ripple and thus we take the average of the steady state voltage values. The curve is associated with the numerical values while the points are the experimental data. Again, the onset velocity of the instability is well predicted and there is also very good qualitative and quantitative agreement between the results. The largest discrepancy still occurs between 4 and 5 m s −1 , with the maximum error of 25% at U = 5.14 m s −1 .
Also, in the case of the standard DC interface, the extracted power exhibits a maximum as a function of the resistance, as shown in figure 13(b). A qualitative comparison between the numerical (curve) and experimental (dots) results shows an acceptable agreement. From a quantitative point of view, the predicted power values are close to the experimental data, with maximum relative error of 28% for R = 1 MΩ. However, the optimal resistance shows a considerable discrepancy: numerical simulations provide an optimum load value around 130 kΩ, while experimentally this is found to be around 82 kΩ.

P-SSHI
In this case, it was necessary to add a resistor in series with the inductance within the Simscape representation of the electrical interface in order to account for the non-idealities that appeared in the real circuit. When observing the waveform of the voltage across the piezo elements, it became evident that the inversion of the voltage had limited efficiency due to losses associated with internal resistances of the components, breadboard, cables, etc.
From figure 14(a) it can be seen that, similarly to the other two circuits considered above, there is very good agreement between the numerical and experimental U − V load bifurcation curve. For this interface, the maximum error amounts to 21% in U = 4.7 m s −1 .
In figure 14(b) the behavior of the power against the resistance is shown. Generally, it is possible to observe how  the numerical model overestimates the values of P for all the considered resistances, with a discrepancy that tends to become more and more evident for high values of R. This effect may arise from the fact that the losses in the real circuit tend to increase the closer we get to the open circuit condition, and this performance degradation remains unexplained and is not captured by the proposed mathematical model.

Optimization of the harvesting device
Once the numerical model is validated, it is used to perform an optimization of the harvesting device. First, a mechanical optimization is performed by changing two geometrical features of the piezoelectric beam and obtaining the trend of efficiency and extracted power with respect to this control parameters. Next, a preferred configuration is selected, and the electrical interfaces are optimized and compared between each other.

Optimization of the mechanical system
For the sake of simplicity, the optimization of the mechanical domain is carried out by considering the standard AC interface as electrical circuit. The geometric features chosen as control parameters are the length of the beam and the side of the square prism. The prismatic body was assumed to be made of polystyrene foam: changing the side of the square section, the mass of the prism changes accordingly. The values of M 1 , ω 1 and χ are derived for each pair of geometric parameters. The length of the beam was varied between 14 and 18 cm, while the side of the square between 3.5 and 5.5 cm. It is assumed that, inside these ranges, the modal damping does not vary significantly from the condition measured on the original beam and therefore can be considered as a constant: the electromechanical system is quite sensitive to variations of ζ 1 . Each pair of parameters yielded a galloping-to-VIV velocity ratio below 2, ensuring the interaction between the two aeroelastic phenomena.
Next, we exploit the Simulink-Simscape model to obtain the maximum power and efficiency for each couple of mechanical parameters. We define the efficiency as done by Shu and Lien [25] ε = P P in (5.1)  where, P is the time-averaged power dissipated across the load resistor and P in is the time-averaged power done by the external force, that in our case is the aerodynamic one (P in = To find the maximum values of power and efficiency for each pair of mechanical parameters, we also varied the load resistance R. Thus, for each couple L b − B the maximum values of power and efficiency with the corresponding optimal resistances are extracted and this process is iterated at 4, 6 and 8 m s −1 . Since we are interested in optimizing the performance of the harvesting device over the whole speed range of 0-8 m s −1 , and since the optimal resistances do not vary significantly over this velocity range, we take the mean of the efficiency and power at these three velocities. The results are shown in figure 15. The graphs show that the highest average efficiency (ε max = 0.137) is provided by the configuration L b = 14 cm, B = 3.5 cm. Instead, the configuration that guarantees the maximum average power (P max = 22.4 mW) is L b = 14 cm, B = 5.5 cm. Although the second configuration had a slightly lower average efficiency (ε = 0.126) than the optimal value, it produced significantly higher power output, making it the selected choice.

Optimization of the electrical system
In this section, we focus on optimizing the electrical domain of the energy harvester with the new piezoelectric beam (L b = 14 cm, B = 5.5 cm). We consider the standard DC interface and the P-SSHI in their optimized configurations and compare them in order to get a sense of the maximum power that can be extracted by the generator. The standard AC interface is not discussed because, being a source of an AC voltage, it is not capable of powering typical electronic devices by itself and would still require the presence of a rectifier.
As far as the standard DC interface is concerned, there are two degrees of freedom: the capacitance of the filter capacitor and the load resistance. The filter capacitance should be kept in a range that allows for a smooth voltage across the resistor while having a short transient time. However, from our numerical simulations, we found that power output had no significant dependence on the capacity at different flow speeds. Regarding the resistive load, it was observed that the power exhibits a maximum around 160 kΩ and that such value does not vary significantly with the wind speed.
For the self-powered SSHI interface, we have as free parameters those related to the electronic breaker, the inductance, the filter capacitance, and the load resistance. From our simulations, the electronic breaker resistances and the filter capacitance indeed do not seem to affect the output power in a significant manner, and the same is true for the inductance as long as it is selected to be sufficiently large. On the other hand, the load resistance influences the output power, and its ideal value varies with wind speed, going from 300 kΩ at 4 m s −1 to 365 kΩ at 8 m s −1 . Therefore, it is not possible to choose a unique optimal value: we selected R = 330 kΩ, for which we have found an excellent performance over the entire speed range considered.
Finally, a comparison between the different optimized interfaces can be drawn, by looking at the output power as the wind speed increases ( figure 16). As can be seen, the P-SSHI interface shows significantly higher performance than the standard DC interface, with the former averaging around 2.15 times the power output of the latter and a maximum power produced of 63.42 mW.

Conclusions
In the present work a numerical model has been developed to predict the behavior of a piezoelectric generator based on the interaction between VIV and galloping even in the case of complex harvesting interfaces. A prototype of the generator has been subsequently manufactured and wind tunnel tests have been performed. The comparison between numerical and experimental results has provided good feedback, thus validating the proposed model. Finally, the developed model was used to perform a parametric analysis aimed at optimizing both the mechanical and electrical domain. From the obtained results, it was found that the P-SSHI interface guarantees on average 2.15 times the output of the standard DC interface, with a maximum power of 63.42 mW at 8 m s −1 .
It can be emphasized that the developed numerical model is relatively simple to implement and allows to obtain results in a short time. Such model can be useful in a preliminary phase of design to obtain some initial data, while a more accurate analysis may be useful in a successive stage.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).