On the effectiveness of a rod-like distributed piezoelectric controller in preventing the Hopf bifurcation of the visco-elastic Beck’s beam

The linear stability of a piezo-electro-mechanical (PEM) system subject to a follower force is here discussed. The mechanical subsystem is constituted by a linear visco-elastic cantilever beam, loaded by a follower force at the free end. It suffers from the Hopf bifurcation, whose critical load is strongly affected by damping, according to the well-known Ziegler’s paradox. On the other hand, the electrical subsystem consists of a distributed array of piezoelectric patches attached to the beam and connected to a properly designed second-order analog circuit, aiming at possibly enhancing the stability of the PEM system. The partial differential equations of motion of the PEM system are discretized by the Galerkin method. Linear stability analysis is then carried out by numerically solving the associated eigenvalue problem, for different significant values of the electrical parameters. A suitable perturbation method is also adopted to detect the role of the electrical parameters and discuss the effectiveness of the controller.

However, while several efforts have been made in the literature to enhance the response of non-autonomous (i.e., externally excited) systems, the vibration control of structures, subject to non-conservative forces (e.g., follower forces) has not been yet extensively addressed. However, it has to be remarked that problems involving follower forces have triggered in the literature a controversial debate [51][52][53]. In particular, several works investigated the existence in the real-world applications of follower actions; some relevant examples in engineering applications can be found: (a) in aerospace [54][55][56]; (b) in vehicle brakes [57,58], and in flexible pipes conveying a flowing fluid [59][60][61][62][63][64][65][66]. On the other hand, also experimental works have been performed to reproduce follower forces in order to demonstrate their existence and the theoretically derived findings, see configuration was derived in [28] and adopted in [87] for similar controlling purposes: it is referred to as the rod-like controller, since the spatial derivatives of the damping-and stiffness-like terms in the flux-linkage equation resemble that of the rod, i.e., a linear elastic beam undergoing only axial displacements. According to the definition given in [28], the piezoelectric devices, shunted to the electrical circuit, are idealized as an array of infinite in-parallel RCL elements; in particular, the circuit is characterized by a linear density of piezoelectric capacitance C, of inductance L, and of resistances R, r R , while has a piezoelectric coefficient E em , that actually provides the coupling between the electrical and mechanical subsystems.
The partial-differential equations of motion, governing the linear dynamics of the PEM system [87], are discretized by a Galerkin approach, and in non-dimensional form, they read (see "Appendix A" for details): M mqm + B mqm + K m q m + 2μHq m − γ G Tq e = 0, ν e M eqe + B eqe + κ e K e q e + γ Gq m = 0 (1) where the dot denotes differentiation with respect to the non-dimensional time t; q m , q e are the vectors collecting the modal coordinates of the mechanical and electrical subsystem, respectively; M θ , K θ (with θ = m, e) are the mass and stiffness matrices of the mechanical and electrical subsystem, respectively, H is the external action matrix, and G is the matrix associated with the gyroscopic nature of the electro-mechanical coupling. On the other hand, the damping operators are B θ = α θ K θ + β θ M θ (θ = m, e), where α θ and β θ (θ = m, e) represent the non-dimensional internal and external (mechanical and electrical) damping coefficients, respectively; ν e and κ e are here referred to as the non-dimensional electrical mass and stiffness, respectively. Finally, γ is the coupling parameter, and μ is the magnitude of the non-conservative force. All the coefficients are defined in "Appendix A."

The control strategy
The linear stability of the PEM system is addressed by solving the eigenvalue problem associated with Eq. (1): the eigenvalues are actually a function of the electro-mechanical parameters, i.e., (μ, α m , β m , γ, ν e , κ e , β e , α e ), and the aim is to investigate how the capability of the rod-like controller affects the PEM behavior. By letting q m = e λt u m and q e = e λt u e , the following eigenvalue problem associated with Eq. (1) is obtained: When γ = 0 and the load μ increases from zero, the uncontrolled system encounters a Hopf bifurcation, i.e., the first eigenvalue, namely the one with the lowest frequency, crosses from the left the imaginary axis of the complex plane, ω d , μ d being the critical frequency and the corresponding critical load, respectively (subscript d means 'damped'). When γ is larger than zero, a complex (and possibly beneficial in terms of stability) interaction between the mechanical and electrical subsystems occurs. Remarkably, and according to what discussed in [29,86], the electro-mechanical coupling γ is assumed to be small, i.e., γ << 1, to meet common practical applications requirements. However, as it was shown in [87], to significantly improve the controller performance, γ should be the largest possible, within the situation of moderately largely coupled systems. However, when γ is sufficiently small, the behavior of the PEM system can be synthesized as follows.
• The uncontrolled beam loses the stability at its critical load μ d and thus oscillates at the corresponding frequency ω d .
• The gyroscopic coupling brings the mechanical response into the controller equation as a forcing term with frequency equal to ω d : the controller starts oscillating as well. • Again the gyroscopic coupling returns the controller response back to the beam: its response is then modified.
Therefore, it is important to remark that the magnitude of the electrical response and thus its contribution to the overall stability of the PEM system depends on the coupling and on the electrical parameters. If they are correctly designed, the controller may enhance significantly the beam response. To achieve this goal, it is expected that the energy exchange between mechanical and electrical subsystems is maximized when the electrical frequency is close to ω d , thus giving rise to a so-called resonant controller, the S RC. However, as discussed in [29,86], this may not be necessarily the best control approach, in fact non-resonant controllers, S N RC, may be more effective in a wider range of the electric parameter space. In particular, it has been observed in [86] that the best performance is obtained when the controller has small ν e , κ e , and is away from resonance, and, in addition, when the electrical frequency is smaller than ω d . This aspect has also been confirmed when the controller has not only a single frequency (see [86]), but it becomes a multi (or infinite)-degree of freedom system (see [87]), as in the case at hand, where the controller is characterized by a multi-modal response whose (electrical) undamped natural frequencies are defined, thanks to the rod analogy, as: Finally, it is worth noticing that the eigenvalue problem associated with the PDEs (A.1) may be also directly attacked by BVP solvers. However, in this work, a Ritz-Galerkin discretization approach of the equations of motion (A.1) has been preferred since, even if less accurate when compared to former solvers, it allows a more straightforward comprehension of the role played by each mode of the uncoupled subsystems on the whole PEM stability.

Sensitivity analysis
To investigate the role of the electrical parameters for both the resonant and non-resonant controllers, a sensitivity analysis based on a perturbation approach is carried out. The procedure is driven from the literature [82] (to which the reader is referred to) and here is briefly recalled (see also "Appendix B" for details).
When 0 < γ < < 1, the eigenvalue governing the PEM system stability can be sought as a perturbation of the critical eigenvalue at the Hopf bifurcation, which, due to the smallness of the coupling, is that of the mechanical subsystem, namely λ 0 = iω d . However, because two different control strategies are adopted, two separate perturbation schemes need to be defined, since for the S N RC the critical eigenvalue is simple, while for the S RC it is defective [82].

SRC
The small positive scaling parameter 0 < ε ≤ 1 (to be reabsorbed at the end of the procedure) is introduced to rescale the electrical parameters as follows: γ → εγ , ν e → εν e , κ e → εκ e , β e → ε 3/2 β e , and α e → ε 3/2 α e . It is remarked that, here κ e and ν e are fixed to enforce the resonance condition and thus tune a specific controller frequency (remember Eq. (3)) such that ω e,k = ω d . Therefore, λ 0 = iω d is a double (defective) eigenvalue for the uncoupled PEM system. Accordingly, a Newton-Puiseux series is adopted: u m = u m,0 + ε 1/2 u m,1/2 + εu m,1 + O ε 3/2 , u e = u e,0 + ε 1/2 u e,1/2 + εu e,1 + O ε 3/2 (4) where λ k , u m,k , u e,k (k = 1/2, 1, 3/2, . . . ) are the eigenvalue sensitivities evaluated at μ 0 = μ d , i.e., the Hopf critical load of the mechanical subsystem. It is worth noting that not only the eigenpairs, but even the load has been expanded around the critical value of the uncontrolled system. Thus, the value of μ, which renders the controlled system marginally stable, is determined as an expansion from μ 0 . It is found that the sensitivity λ 1/2 governs the stability of the PEM system whose expression is derived from the solvability condition enforced at the ε-order problem, and it is the solution of the following second-order algebraic equation: where the coefficients c 0 , c 1 are defined in "Appendix B." The latter equation admits two distinct roots λ ± 1/2 = λ ± 1/2 μ 1/2 , ν e , β e , α e , γ that do not explicitly depend on κ e , since the ratio κ e /ν e has been fixed to enforce the resonance condition.
It is concluded that the PEM system endowed with the S RC controller is asymptotically stable when Re λ ± 1/2 < 0. The boundary of the stability domain is thus defined by Re λ ± 1/2 = 0, which, in a geometrical perspective, is a hyper-surface in the 5-D parameter space μ 1/2 , ν e , β e , α e , γ .

SNRC
The same rescaling of the electrical parameters of the S RC is taken, except for the electrical damping coefficients, that is: β e → εβ e and α e → εα e . Since λ 0 is a simple eigenvalue for the uncoupled PEM system, the following Mc Laurin series expansion is adopted: Also in this case, the load μ is expanded around the critical value of the uncontrolled system μ 0 = μ d , and the value μ 1 represents the sensitivity of the PEM critical load due to the presence of the controller. Under these assumptions, the perturbation algorithm is carried out, and the solvability conditions on the descending ε-order problem furnishes the first eigenvalue sensitivity, namely where g mm (ν e , β e , α e , κ e ) := (v 0 m,0 ) H G T S −1 e Gu 0 m,0 and S e , u m , v m , h mm , b mm , m mm are defined in "Appendix B," together with some details about the perturbation scheme.
It is concluded that the PEM system endowed with the S N RC controller is thus asymptotically stable when Re (λ 1 ) < 0. The boundary of the stability domain is identified by the equation Re (λ 1 ) = 0 which, from a geometrical point of view, represents a 6-D hyper-surface depending on the parameters (μ 1 , ν e , β e , α e , κ e , γ ).

Preliminary analysis on the controller frequency
In this Section, a preliminary sensitivity analysis on the PEM stability, with respect to the non-dimensional electrical stiffness and coupling coefficient γ (keeping fixed the other electrical parameters), is carried out with the purpose of exploring the effects of the resulting interaction between the spectra of eigenvalues of the mechanical and electrical subsystems.
The beam parameters, here and in what follows, are taken as α m = 0.01 and β m = 0.1, which correspond to a situation where the system suffers from a considerable detrimental effect due to the 'Ziegler paradox', entailing that damping reduces 30% of the critical load, i.e., μ d = 6.46 (being 10.02 that of the undamped beam) and ω d = 5.92.
The eigenvalue problem Eq. (2) is solved by adopting N m = N e = 9 (see "Appendix A"), and the stability diagrams, namely the diagrams of the region of the parameter-space which characterize the PEM system stability, are obtained. These diagrams are represented in terms of percentage deviation of the load μ from μ d , namely μ = 100(μ − μ d )/μ d . The eigenvalue problem is solved by varying two selected electrical parameters, namely κ e ∈ (0, 5] and μ ∈ [−5, 30]%, and for selected values of γ , namely γ = When γ = 0.0125, it is sufficiently small such that the PEM system can be considered as weakly-coupled. The results of the sensitivity analysis are illustrated in Fig. 2 where the stable regions of the stability diagram (i.e., all the eigenvalues have negative real part in these regions of the parameter-space) are denoted in blue, and the dashed lines identify the values of κ e that correspond to the resonance condition ω e,k = ω d of the k-th electrical mode. It is observed that, when the first electrical mode is tuned to the beam frequency (ω e,1 = ω d ), the controller has a considerable detrimental effect since the system is stable for negative μ. On the other hand, when κ e is larger, the controller effect is negligible ( μ 0), while if κ e is smaller significant stable regions can be detected, i.e., μ > 0, entailing that the controller is beneficial. In particular, the stability diagram exhibits multiple peaks when the resonance condition is attained on higher electrical modes, with the largest increase of critical load when ω e,2 = ω d . It is then concluded that for small γ the best control strategy, except for the resonance on the first electrical mode, corresponds to the S RC approach.
A second analysis is performed by increasing the coupling coefficient γ , so that the PEM system becomes moderately coupled. Analogously to what was shown before, the results are reported in Fig. 3 for the values of the coefficient: γ = 0.025 in Fig. 3a and γ = 0.05 in Fig. 3b.
Also in this case, the tuning of the first electrical mode to the the beam frequency reveals to have a strong detrimental effect, and analogously if κ e increases, the controller effect becomes negligible. On the other hand, if κ e decreases, again a large stable region at positive μ appears showing multiple peaks for specific values of κ e . However, it must be noted that, contrarily to the previous case, the largest increments of μ are attained when the electrical system is not in resonance with the beam. In particular, when γ = 0.025 (see Fig. 3a), the largest peak is reached when κ e is slightly larger than that corresponding to ω e,2 = ω d , i.e., in a nearly resonance condition. Whereas, when γ = 0.05 (see Fig. 3b), the largest peak is reached when κ e is considerably far from that corresponding to ω e,2 = ω d , as well as far from ω e,1 = ω d : the S N RC strategy is then more effective. In addition, in both cases, the increase of critical load is considerably larger than that obtained at smaller γ .
By summarizing, the preliminary analysis revealed that when γ is small S RC strategy is the best choice for improving the PEM stability, while if γ is moderately large, S N RC strategy becomes more effective.
However, it is of interest to further shed light on the complex interaction taking place between the electrical spectrum and the mechanical one, as the controller parameters vary. An attempt to discuss this aspect can be made by carefully analyzing the stability diagrams shown in Fig. 3. Preliminarily, it is remarked that only the first mode of the uncontrolled beam, i.e., that at lowest frequency, participates to the PEM stability, the higher ones being passive, for the considered load values. On the other hand, more than one electrical mode is found to be involved the PEM stability depending on the adopted "tuning," i.e., on the choice of the electrical stiffness κ e : the lower κ e , the higher the electrical mode possibly triggering instability. In addition, the stability of the PEM system is strongly affected by the choice of the coupling parameter γ (having set α e , β e ), in particular: when γ is small (see Fig. 3a), a weak interaction between the spectra of the uncoupled subsystems manifests itself, and the mode determining the Hopf bifurcation comes from the mechanical one (the background in the Figure is denoted in white); on the other hand, when γ is larger (see Fig. 3b), a stronger spectra interaction takes place, since the stability is governed by the mechanical subsystem for large κ e (white background), but as κ e decreases the eigenvalue triggering the PEM instability comes from the electrical subsystem (light gray background). It is worth to note that in some regions (see the dark gray area) an interaction between the mechanical and electrical modes occurs for large κ e , as well as for low κ e between electrical modes (not reported as gray scales region in the Figure). An in-depth discussion about this modal interaction will be given in Sect. 6.3.

Numerical results
In this Section, further sensitivity analyses are carried out to understand the role played by the electrical damping parameters. The linear stability analysis of the PEM system is here carried out via the previously discussed perturbation approach, whose results are compared with those obtained by the numerical analysis of the eigenvalue problem of the coupled PEM system.

Effects of the electrical damping
By exploiting the analytical expressions derived in Sect. 4, the sensitivity to the electrical damping is investigated for both the controllers: the S RC (see Sect. 4.1) and the S N RC (see Sect. 4.2). These are endowed with a source of electrical external, as well as internal damping represented by the parameters β e and α e , respectively. The stability diagrams are derived by fixing the electrical stiffness in correspondence of the largest peak appearing in Fig. 2 for the S RC and in Fig. 3a for the S N RC. The results are illustrated in Fig. 4, where the stability boundary is represented by a surface in the (β e , α e , μ)-space, the stable region being that below the surface. Moreover, the other electrical parameters are chosen according to the results derived in the previous Section, namely ν e = 0.1 and γ = 0.0125, κ e = 0.157 for the S RC (Fig. 4a), and γ = 0.025, κ e = 0.177 for the S N RC (Fig. 4b).
The obtained behavior is qualitatively similar for both the controllers: the increase of critical load is larger when both β e , α e are small; in particular, it is observed that the controllers are effective when β e = 0, α e = 0 or β e = 0, α e = 0. Thus, to improve the controller performance, the electrical external and internal damping should act separately.
Then, slices of the surfaces shown in Fig. 4 are taken at the same κ e , ν e and slightly different γ , to further emphasize the effect of the coupling and electrical damping parameters on the PEM stability. Results relevant to the PEM endowed with the S RC are shown in Fig. 5, by taking ν e = 0.1, κ e = 0.157, γ = 0.01, 0.0125, 0.15; the stability diagram is represented: in the ( μ, β e )-plane with α e = 0 (see Fig. 5a), and in the ( μ, α e )-plane with β e = 0 (see Fig. 5a). Similarly, the results of the PEM endowed with the S N RC are shown in Fig. 6, by taking ν e = 0.1, κ e = 0.177, γ = 0.02, 0.025, 0.03, in the ( μ, β e )-plane with α e = 0 (see Fig. 6a), and in the ( μ, α e )-plane with β e = 0 (see Fig. 6b).
In both the cases, the behavior is qualitatively similar: when the electrical damping (analogously for β e and α e ) is close to zero, the controller effect is negligible, i.e., the system is stable at μ 0, then large stable regions (at the right of the boundary) are obtained until an optimum value is found above which μ decreases. Here, the small increase of γ (close to the reference value) induces a shift of the stability boundaries toward higher μ, thus confirming that, in order to improve the control performance, higher electrical coupling coefficients are preferred, but with the caveat that if γ is strongly enlarged, the control strategy should change (remember Figs. 2, 3).

Analytical versus numerical results
Here, stability diagrams obtained via the perturbation approach are compared with those got by numerical analyses, directly carried out on the eigenvalue problem. This analysis aims to check the validity of the asymptotic approach and to confirm the qualitative findings previously discussed.  Fig. 7b. It can be observed that the behavior is qualitatively well captured by the asymptotic approach, though a discrepancy is quantitatively exhibited. This is due to the fact that the range of β e , α e actually slightly exceeds the parameters scaling grounding the perturbation method, introduced in Sect. 4. The ordering is respected in the insets of Fig. 7 that show a zoom in the region close to the origin: a better agreement is reached when the damping value is below the cusp of the domain, above which the ordering is violated (according to [82]).
An analogous comparison is performed at larger γ , thus entailing that the PEM is endowed with the S N RC. The results are illustrated with the same previous logic in Figs. 8 and 9. In Fig. 8, the stability diagrams are determined by taking γ = 0.025, ν e = 0.1, κ e = 0.177 and are shown: in the ( μ β e )-plane with α e = 0 (Fig. 8a); in the ( μ, α e )-plane with β e = 0 (Fig. 8b). Contrarily to what was shown for the S RC, when the coupling become stronger, the analytical results are affected by a considerable discrepancy that is clear in terms of β e , while is less evident in terms of α e . On the other hand, by evaluating the stability diagrams after a further increasing of gamma, namely γ = 0.05 and by taking ν e = 0.1, κ e = 0.2, a much more evident loss of accuracy of the perturbation approach is detected (see Fig. 9), but this is not only due to the fact that the electrical parameters are beyond the limit suggested by the ordering. In fact, it is expected that when the PEM is moderately coupled, the system response cannot be assumed as a small perturbation of the response of the beam alone. Thus, again confirming that the stronger the coupling is, the stronger the interaction between the spectra of electrical and mechanical subsystems is, and, accordingly, the instability of the PEM system is triggered by more than one mode (i.e., the first beam mode), and, in particular, that the electrical modes are no more passive.  It can be observed that the behavior of the mechanical subsystem qualitatively resembles the first mode of the cantilever beam (see Fig. 10a), even though quantitatively the solutions in C 1 , C 3 are fully overlapped, while C 2 is slightly different. On the other hand, the behavior of the electrical subsystem resembles the second electrical mode (according to the choice of the parameters) in all the considered conditions; however, a significant quantitative difference can be observed, suggesting a different role of the electrical mode in the PEM stability.

Discussion
To further investigate the mechanical/electrical modal interaction on the stability of the PEM system, its eigenvalues are numerically evaluated by solving the eigenvalue problem Eq. (2) over the full range of μ ∈ [0.95, 1.3]μ d and for the same parameters of Fig. 9a. The real and imaginary part of λ are reported versus μ in Fig. 11a, b, respectively. As shown in Fig. 11a, it is clear that there is more than one eigenvalue crossing the horizontal axis, confirming the fact that the PEM instability can be triggered by two modes. The nature of such modes is easily recognized by analyzing the results illustrated in Fig. 11b: the first eigenvalue attaining the instability at a lower critical load (that corresponding to C 1 ) is related to the first beam mode (ω m,1 ω d ) and is the same determining the stability in C 2 . On the contrary, the stability in C 3 is governed by the second electrical mode ω e,2 > ω d (here κ e = 0.2, remember Fig. 3b).
A better understanding of these findings merges when a projection of the modes of the coupled system onto the basis of the eigenvectors of the uncoupled one is carried out. Indeed, the solution of Eq. (2) delivers the eigenvalues and the corresponding eigenvectors, namely λ and u m , u e , of the coupled PEM system. On the other hand, the eigenvalue problem of the uncoupled mechanical and electrical subsystems (γ = 0) is governed by: where u 0 m , u 0 e, are the eigenvectors of the uncoupled mechanical and electrical subsystems, respectively. Then, the eigenvectors of the coupled PEM system are written as a linear combination of the eigenvectors of the uncoupled subsystems, namely where a j and b j are the (unknown) coefficients of the linear combination. Projection of the eigenvectors of the coupled PEM system onto the basis of those of the uncoupled subsystems calls for the left eigenvalue problems associated with Eq. (1) with γ = 0, which read: where the superscript () H denotes the conjugate transpose and v 0 m , v 0 e are the left eigenvectors of the uncoupled mechanical and electrical subsystems, respectively. Then, by premultiplying both members of Eq. (9) by the k-th corresponding left eigenvector, the projection reads: where use of the orthogonality between the left and the right eigenvectors (the scalar product in the sum is zero when k = j) has been made; moreover, by adopting the normalization criterion (v 0 θ,k ) H · u 0 θ,k = 1 (θ = m, e), the a k , b k coefficients can be obtained according to the following expressions: According to the latter definitions, a k , b k represent the participation factor of the k-th mode of the uncoupled subsystem to the response of the coupled PEM system. The evaluation of the a k , b k coefficients is performed on the points of the stability boundary denoted by C j with j = 1, 2, 3 in Fig. 9. Since they are complex-valued coefficients, their real and imaginary parts are separately illustrated, namely Re(a k ) and Re(b k ) are shown in Fig. 12a, c, respectively, whereas Im(a k ) and Im(b k ) are reported in Fig. 12b, d, respectively. As expected, the larger modal participation factors are obtained on the first beam mode, see Fig. 12a, b, meaning that the mechanical part of the PEM response, u m , is mostly influenced by the first beam mode, though nonzero components of the second mode are also present. Similarly, it is observed that the electrical part of the PEM system response, u e , is mostly driven by the second mode of the controller, see Fig. 12c, d, with smaller but not negligible components of the first controller mode. This confirms that the PEM stability cannot only be evaluated as a small perturbation of the first beam mode, since, when γ is moderately large, more than one mode may trigger the system instability.

Conclusions
The stability of a piezo-electro-mechanical (PEM) system has been investigated in this paper. The system is composed of a mechanical subsystem, i.e., a visco-elastic Beck's beam, and a piezoelectric subsystem that is a rod-like controller. The beam suffers the detrimental effect of the Ziegler paradox on the Hopf bifurcation, triggered by a follower force; thus, the PEM system behavior has been investigated to evaluate the capability of the controller to improve its stability. The discretized equations of motion have been recalled from the literature, and the eigenvalue problem has been solved numerically as well as via an analytical perturbation approach based on the smallness of the coupling and electrical parameters. The sensitivity of the PEM system to the main electrical parameters has been discussed, and the contribution of the electrical modes to the system stability has been analyzed. The major outcomes are summarized below.
1. The most efficient control strategy depends on the magnitude of γ : if it is small, the S RC is more effective, while when it increases, an S N RC strategy is preferable. However, the larger γ , the larger is the beneficial effect of the controller in increasing the system Hopf critical load. 2. The electrical external and internal damping can be both beneficial, but should act separately to avoid detrimental effects on the stability. 3. The perturbation approach gives a good agreement for weakly-coupled systems, while, as it is expected, it loses accuracy when γ grows toward moderately coupled systems. In that case only the numerical approach is able to well represent the PEM behavior. 4. In moderately coupled PEM systems, the stability is governed by the interaction of multiple modes: the first descends from the mechanical subsystem, while the other derives from the controller.
Further investigations may be performed on other types of controllers, i.e., different analog circuits, to analyze their effect on the system stability. Moreover, another interesting aspect which deserves to be investigated is the nonlinear response of the PEM system in the post-critical regime also in the presence of nonlinear electrical damping.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open access funding provided by Universitá degli Studi dell'Aquila within the CRUI-CARE Agreement. The authors received no additional funding for this research work.

Declarations
Competing interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
where the dot denotes differentiation with respect to the non-dimensional time t and the prime with respect to the non-dimensional abscissa s. Accordingly, the non-dimensional parameters governing the PEM behavior are defined as: where the subscripts m, e denote the terms relative to the mechanical and electrical subsystems, respectively. The variables v and ψ are thus expressed as a linear combination of N m and N e terms, respectively, each representing the product of the k-th shape function φ (θ ) k (s), collected in the vector θ (s), and the corresponding time-depending amplitude q (θ ) k (t), collected in the vector q θ (t), with θ = m, e. In [87], a variational approach is adopted to the derive the discretized PEM model, and to obtain Eqs. (1), the following operators are defined: In the latter, the shape functions adopted for v are chosen as the eigenfunctions of an unloaded undamped Euler-Bernoulli beam undergoing flexural motion, while for ψ the analogy with the rod is recalled; consequently, the eigenfunctions of an unloaded undamped rod undergoing axial motion are adopted.

Appendix B: Details on the perturbation sensitivity analysis
Details about the perturbation approach to the eigenvalue sensitivity analysis are given in the following. For the complete formulation of the asymptotic procedure here recalled, the reader is referred to [82]. Besides the different parameters scaling and the variable expansion, the leading-order problem possesses the same structure in both the cases of S RC and S N RC, namely rendered unique by a suitable normalization. It should be remarked that, at this order, λ 1/2 is still undetermined. The solvability condition must instead be invoked at the ε-order problem, which finally determines λ 1/2 , according to Eq. (5), where the involved coefficients are defined as: c 0 := γ 2 λ 0 g em g me + μ 1/2 h mm b ee 2ν e m ee (2λ 0 m mm + b mm ) . (B.13)

SN RC
For the S N RC, the leading-order problem takes again the expression given in Eq. (B.5), while the higher-order problem reads: Although the mechanical and electrical responses are of the same order, the electric oscillator behaves as passive, driven by the mechanical one. By substituting expressions (B.15) into the ε-order problem and enforcing its solvability the expression of λ 1 given in Eq. (7) is obtained.