Asymptotic and Spectral Analysis of a Model of the Piezoelectric Energy Harvester with the Timoshenko Beam as a Substructure

Mathematical analysis of the well known model of a piezoelectric energy harvester is presented. The harvester is designed as a cantilever Timoshenko beam with piezoelectric layers attached to its top and bottom faces. Thin, perfectly conductive electrodes are covering the top and bottom faces of the piezoelectric layers. These electrodes are connected to a resistive load. The model is governed by a system of three partial differential equations. The first two of them are the equations of the Timoshenko beam model and the third one represents Kirchhoff’s law for the electric circuit. All equations are coupled due to the piezoelectric effect. We represent the system as a single operator evolution equation in the Hilbert state space of the system. The dynamics generator of this evolution equation is a non-selfadjoint matrix differential operator with compact resolvent. The paper has two main results. Both results are explicit asymptotic formulas for eigenvalues of this operator, i.e., the modal analysis for the electrically loaded system is performed. The first set of the asymptotic formulas has remainder terms of the order O( 1 n ), where n is the number of an eigenvalue. These formulas are derived for the model with variable physical parameters. The second set of the asymptotic formulas has remainder terms of the order O( 1 n2 ), and is derived for a less general model with constant parameters. This second set of formulas contains extra term depending on the electromechanical parameters of the model. It is shown that the spectrum asymptotically splits into two disjoint subsets, which we call the α-branch eigenvalues and the θ-branch eigenvalues. These eigenvalues being multiplied by “i” produce the set of the vibrational modes of the system. The α-branch vibrational modes are asymptotically located on certain vertical line in the left half of the complex plane and the θ-branch is asymptotically close to the imaginary axis. By having such spectral and asymptotic results, one can derive the asymptotic representation for the mode shapes and for voltage output. Asymptotics of vibrational modes and mode shapes is instrumental in the analysis of control problems for the harvester.


Introduction
Energy harvesting, an extremely popular topic in contemporary engineering literature and practice, is understood as the process and result of converting the energy available in the environment into electrical energy, which can be consumed or stored for later use.There are different types of vibration-based energy harvesters, e.g., piezoelectric, electromagnetic and electrostatic.An advantage of a piezoelectric vibration energy harvester is the ability of piezoelectric material to convert the mechanical energy directly to electrical energy without external input.This in turn allows simpler practical designs for piezoelectric energy harvesters compared to their electromagnetic and electrostatic counterparts.Piezoelectric materials create an electric charge in response to a mechanical stress (the direct piezoelectric effect) and, in the reciprocal process, these materials develop mechanical strain in response to the electrical potential (the inverse piezoelectric effect).Since the charge generated in the piezoelectric material is proportional to the applied stress, an energy harvester is designed to maximize stress under a certain mechanical load.The most commonly used geometric structure is the cantilever beam, since this structure provides the highest average strain for a given input force.Ambient vibrations are coupled to the cantilever-mass system through the base of the cantilever causing the structure to oscillate.The alternating bending strains, due to oscillations, are converted to an alternating voltage by the piezoelectric material.To maximize the voltage output, ambient vibration frequencies should be close to the natural oscillation frequency of the beam.
In the present work, we analyze a mathematical model of the harvester.This model consists of a system of three coupled partial differential equations describing the dynamics of a cantilever Timoshenko beam whose top and bottom faces are covered by very thin piezoelectric layers.The system of equations is equipped with the standard sets of the boundary and initial conditions (see Equations ( 1)-( 5) of Section 2).It is important to emphasize that in our previous works on energy harvesting models (see Shubov [1][2][3], Shubov and Shubov [4]), we considered the Euler-Bernoulli beam model as a vibrating structure.In the present work, we deal with the Timoshenko beam model, which makes mathematical treatment much more challenging.Before we turn to the mathematical statements and results, we would like to mention several research directions closely related to the topic of the present paper.
The development of micro and nano technologies, smart sensors, and health monitoring devices created the need of independent energy sources for small devices and self-powered sensors.Stoykov et al. [5] studied the electro-mechanical system of vibrational energy harvester, where the beam is excited by external and kinematic periodic forces and damped by an electrical resistor through the coupled piezoelectric transducer.Nonlinearities are introduced by stoppers limiting the transverse displacement of the beam, and the action of the stoppers is modeled as Winkler elastic foundation.
To study the composite beam, the authors used the geometrical nonlinear version of Timoshenko beam theory [5].Using their models, the authors can capture discontinuities of structural parameters such as thickness discontinuities or impact of stoppers, which is important in engineering practice.
Michelin and Doare [6] discussed the way by which self-sustained oscillations resulting from air-solid instabilities, such as flutter of a flexible flag in an axial air flow, can be used for energy harvesting.Piezoelectric patches attached to the surface of the flag transform the solid deformation into an electric current powering purely resistive output circuits.A flexible plate becomes unstable due to flutter above the critical flow velocity which depends on the plate's properties (e.g., density and rigidity) and can be chosen to be lower than the typical flow velocity, which leads to self-sustained large-amplitude flapping of the plate.The study is focused on numerical analysis of the nonlinear dynamics of the system and on determination of its harvesting efficiency and robustness of the process.
In the review paper, Abdelkefi [7] discussed different issues on the energy harvesting from aeroelastic vibrations.These harvesters can be placed in such areas as high wind spaces, ventilation outlets, ducts of buildings, lifting surfaces of aircraft structures, and biomedical implants.Qualitative and quantitative characteristics of various types of flow-induced vibrations energy harvesters are presented in the paper.Limitations and recommendations on mechanical and aerodynamic modeling, power conditioning circuits optimization, and prototypes fabrications are discussed in detail.
Anton and Inman [8] performed a concept study on piezoelectric energy harvesting in unmanned air vehicle (UAV) applications through flight testing of a remote controlled aircraft.It is stated that the choice of their particular plane is of primary importance because of its long wingspan and flexible foam wings, which provide a good environment for piezoelectric vibration harvesting using wing deflections.It is reported that a small battery can be fully charged when increasing the volume of piezoelectric material.Briant and Garcia [9] are among the first researchers who analyzed theoretically and experimentally a two-degree-of-freedom typical section model as a power harvesting device driven by aeroelastic vibrations.In [9], a time-domain switching energy extracting scheme is considered in order to increase the level of the harvested power.A wind tunnel testing demonstrates that the wing section geometry does not affect the linear flutter speed and its frequency.The authors show that the optimal wing profile depends on the operating range of wind speeds and on the required level of harvested power.
Elvin and Elvin [10] performed a linear analysis for a cantilever pipe to investigate the effects of passive piezoelectric damping with a load resistance on the flutter speed.They demonstrated that the mechanical stiffening effects of the open-circuit cause an increase in the flutter speed compared to the short-circuit.The authors showed that the larger the piezoelectric electromechanical coupling coefficient, the higher the flutter speed.For the open circuit case, it is indicated that an increase in the piezoelectric capacitance is accompanied by a decrease in the flutter speed.
Erturk et al. [11] presented a frequency domain analysis and experimental validations for a twodegree-of-freedom typical section as a wing-based piezoaeroelastic energy harvester.They focused on the problem of harvesting energy at the flutter boundary and analyze the effects of the piezoelectric coupling on the linear flutter speed.In their mathematical modeling, they introduced a piezoelectric coupling to the plunge degree of freedom and considered a load resistance in the electrical field.They presented the modified lumped-parameter aeroelastic equations.As for the aerodynamic lift and moment modeling, a two-dimensional unsteady representation is based on the Theodorsen's approach.It is indicated that a good agreement between the mathematical model and the experimental measurements is obtained.It is shown that the generated voltage increases when increasing the electrical load resistance.Furthermore, there is an optimum value of the electrical load resistance in which the harvested power is maximized.
De Marqui et al. [12] presented a time-domain piezoaeroelastic modeling and numerical simulations of a generator wing with embedded piezoceramics for continuous-and segmented-electrode configurations.Their wing-based piezoaeroelastic energy harvester model is obtained by combining an electromechanically coupled finite element model based on the classical plate theory with an unsteady vortex-lattice model [13] representing the aerodynamic loads.They reported that, when using segmented electrodes, torsional motions of the coupled modes become more important and therefore change the flutter speed which improves the broad band performance of the harvester.Doare and Michelin [14] presented theoretical analysis of harvesting energy from the flutter of a flexible plate subjected to an axial flow.They presented local and global stability analyses based on the global coupled electromechanical equations of motion.Other linear analyses are performed in [15,16].The authors presented theoretically and experimentally a guideline to design piezoaeroelastic energy harvesters that can generate energy at low freestream velocities.The effects of varying the mass and stiffness parameters on the linear flutter speed and frequency are investigated.They reported that the linear flutter speed is more sensitive to the flap center of mass location, the flap mass moment of inertia, and the flap mass.In [15,16], a semi-empirical nonlinear model is used to determine the unsteady aerodynamic loads for large flap deflections.
Dias et al. [17] considered an aeroelastic system for which they study the possibility of harvesting energy from both piezoelectric transduction and electromagnetic induction, named hybrid energy harvesting.The authors presented the linear governing equations of a hybrid aeroelastic energy harvester and investigate in detail the effects of the electrical and aeroelastic properties on the linear flutter speed and performance of the harvester at the flutter boundary.It is demonstrated that a reduction in the radius of gyration and frequency ratio and an increase in the eccentricity result in a decrease in the linear flutter speed and an increase in the harvested power at the flutter boundary.A three-degree-of-freedom hybrid (piezoelectric and electromagnetic) aeroelastic energy harvester that includes control surface in the airfoil (flap effect), was considered by Dias et al. [18].A linear coupled lumped-parameter model is presented including both the piezoelectric and electromagnetic couplings.Theoretical results show that the linear flutter speed decreases for large values of the eccentricity and plunge-to-pitch frequency ratio.
Abdelkefi et al. [19] proposed various mathematical models for beam energy harvesters: lumpedparameter and distributed-parameter models.They designed and studied a unimorph cantilever beam-based harvester that undergoes bending-torsion vibrations by creating an offset between its center of gravity and the shear center, thereby leading to a coupling between the bending and torsion vibrations.The offset is created by placing two masses asymmetrically at the tip of the beam.The beam-mass system is modeled using the Euler-Bernoulli beam theory and Hamilton's principle to derive the coupled governing equations of motion and associated boundary conditions for a base excitation.The exact mode shapes and natural frequencies of the harvester are calculated and used as basis functions in Galerkin's scheme to derive a reduced-order model.Closed-form expressions are obtained for all needed outputs.It is demonstrated that, when the lowest global frequencies are close to each other, an interesting broadband frequency harvester appears.It is also shown that, when the asymmetry between the two masses increases, the electrical harvested power as well as the produced voltage increases.
Bibo and Daqaq [20] investigated the harvester which consists of a rigid airfoil supported by nonlinear flexural and torsional springs which is placed in an incompressible air flow and subjected to a harmonic base excitation in the plunge direction.The flow field sets the elastic structure into limit-cycle oscillations that induce an alternating strain, and, hence, a charge in the piezoelectric laminates.The paper investigates the transduction of piezoaeroelastic energy harvesters under the combination of vibratory base excitations and aerodynamic loadings.An approximate analytical solution describing the harvester response is obtained using the center manifold reduction and the method of normal forms.The analytical solution is used in conjunction with numerical simulations to investigate the harvester's performance below and above the flutter speed.
Toprak and Tigli [21] presented a comprehensive review on the history and current state-of-the art piezoelectric energy harvesting.A brief theory section presents the basic principles of piezoelectric energy conversion and introduces the most commonly used mechanical architectures.The theory section is followed by a literature survey on piezoelectric energy harvesters, which are classified into three groups depending on the size of the devise.The size of a harvester affects parameters such as its weight, fabrication method, achievable power output level, and potential application.
Having in mind potential applications of piezoelectric energy harvesters, it is important to extend the thin-beam models to reasonably thick-beam configurations as well as to varying geometries of the harvesters.The dynamics of a beam-like flexible structure strongly depends on its aspect ratio and the operating frequency range.The Euler-Bernoulli model is the classical model for slender beam configurations with sufficiently high length-to-thickness aspect ratio so that the shear distortion and rotary inertia effects can be neglected.The Rayleigh model introduces the effect of rotary inertia to the Euler-Bernoulli model but it neglects the effect of transverse shear distortion.The Timoshenko beam model accounts for both the shear distortion and rotary inertia effects and is widely used for modeling the dynamics of moderate length-to-thickness ratio beams.Erturk [22] presents approximate analytical distributed-parameter electromechanical modeling of cantilevered piezoelectric energy harvesters based on the Euler-Bernoulli, Rayleigh, and Timoshenko beam theories.The technique of [22] is an electromechanical version of the assumed-modes method.After deriving the distributed-parameter energy expressions, the extended Hamilton's principle is employed to obtain the discretized electromechanical Lagrange's equations.An axial displacement variable is kept in the formulation to account for its coupling with the transverse displacement (in the Euler-Bernoulli and Rayleigh models) or cross-section rotation (in the Timoshenko model) due to possible structural asymmetry.The steady-state electromechanical response expressions are obtained for harmonic base excitation.Experimental validations are given for the thin-beam case by comparing the assumed-mode predictions with the experimental and analytical results.
The motivation of the present research is the fact that there exist many well-developed mathematical models of energy harvesters whose rigorous analytical treatment would be very desirable.Our main goal is to present rigorous mathematical analysis of a piezoelectric energy harvester with the Timoshenko beam model as a substructure.We consider two cases: (i) the beam model with variable structural parameters; and (ii) the beam model with constant structural parameters.Our goal is to derive the asymptotic approximation for the infinite set of the natural frequencies of the model.To this end, we perform asymptotic analysis of the corresponding set of three partial differential equations subject to four boundary conditions, and homogeneous initial conditions.We would like to obtain the answer to the following question: In which way do the parameters of piezoelectric layers and the electric circuit (such as the piezoelectric coupling constant θ, the capacitance of the layers C, and the resistance of the external load R) enter the asymptotic formulas describing the distribution of the natural frequencies?The answer to this question is the following: The asymptotic formula for the frequency distribution has the leading term, the second order term, and the remainder term.Based on our analysis, we obtain that the leading asymptotical term contains contribution only from the structural part of the model, the second order term contains contributions from piezoelectric and circuitry parts of the model (see Theorems 5 and 6 below).To prove such results, we have to derive refined asymptotic approximations for the vibrational modes.However, obtaining higher order asymptotics for the case of variable parameters structure is quite involved and extremely lengthy.That is why we consider two models of the harvester, one contains the structural part of the model with variable parameters, for which we derive the first order asymptotic approximation for the modes (that contains only the leading term and the remainder term), and the second model that contains constant parameter structural part of the model, for which we derive higher order asymptotic approximation for the modes (that contains the leading term, the second order term, and the remainder term).Now, we are in a position to describe the content of the present paper.In Section 2, the initial boundary-value problem is formulated in the form of a system of three coupled partial-differential equations (see Equations ( 1)-(5) below), in which Equation (2) contains the distributional terms.We represent the problem in an equivalent form so that there are no distributional terms in the equations but there is an additional term in the boundary conditions.In Section 3, the problem is reformulated and is given in the form of a first order in time evolution equation in the Hilbert state space corresponding to the model.The dynamics generator is the main object of interest for the rest of the paper.Its properties are presented in Theorem 1.The importance of the dynamics generator is that the set of its eigenvalues being multiplied by "i" (the imaginary unit) coincides with the set of the natural frequencies of the harvester model.
In Section 4, we analyze asymptotic and spectral properties of the polynomial operator pencil whose eigenvalues coincide with the dynamics generator eigenvalues.We describe the fundamental system of solutions for the pencil equation.In Section 5, we construct two special matrices, which we call the left reflection matrix and the right reflection matrix.Using these matrices, we can construct the solution of the pencil equation that satisfies four boundary conditions.The asymptotic approximation for the vibrational modes as the number of a mode tends to infinity is justified in Theorem 5 of Section 6.This theorem together with the method of its proof is the first main result of the paper.Finally, Sections 7 and 8 are devoted to the constant parameter case and refined asymptotic approximations for the modes that are proven in Theorem 6.This theorem is the second main result of the paper.The proof is given for the constant coefficient case since the derivation of the refined asymptotic formulas in the general variable coefficient case would require significant increase of the length of the paper.
In the conclusion of the Introduction, we emphasize that the goal of the present manuscript is to provide rigorous derivation of the analytical formulas for two cases, i.e., for the structure with variable parameters and for the one with constant parameters.The derivations of all the formulas are quite lengthy and complicated.The corresponding numerical results will be presented in the forthcoming work.Preliminary numerical results show very good agreement with the analytical formulas obtained in the paper.In author's opinion, analytical formulas are important since they can provide insights not available from purely numerical results.

Statement of the Energy Harvester Problem
In this section, we present a formulation of the initial-boundary value problem describing the dynamics of a specific model of an energy harvester.The electromechanical equations of motion are presented in transverse and rotational vibrations using Timoshenko beam theory.The beam is assumed to be excited by small (not necessary harmonic) transverse motion of the base.The schematic diagram of the piezoelectric, vibration-based energy harvester considered in this study is shown on Figure 1.It is a bimorph configuration consisting of a piezoelectric material layer bonded to both surfaces of a supporting, inactive core.Electrodes are assumed to cover the upper and lower surfaces of each layer and they are wired together in a parallel configuration.The voltage across each layer is assumed to be equal, and the charge displaced by each layer is additive.Because each layer experiences opposite strains (i.e., one layer is in tension while the other is in compression), they must be in the same direction to avoid charge cancellation.It is assumed that the electrodes and connecting wires have negligible resistance and that the resistivity of the piezoelectric material is significantly higher than that of the external circuitry (i.e., an insulator).In what follows, the equations of motion for the electromechanical system are derived through force, moment, and charge balances adopting the Timoshenko beam assumptions (see Inman [23], Timoshenko et al. [24], Shubov [25]).Now, we are in a position to present mathematical statement of the problem.The model is governed by a system of three coupled partial differential equations for the following unknown functions: W(x, t): the relative transverse deflection of the beam with respect to its base; Φ(x, t): the rotation angle of the beam cross-section; and v(t): the voltage across the energy harvester.
To describe the dynamics of the non-homogeneous Timoshenko beam with the piezoelectric patches, we need the following notations: EI(x) is the flexural rigidity of a beam, ρ(x) is the mass density, S is the cross-sectional area of beam, and K(x) is the shear stiffness of a cross-section (as in Benaroya [26], Rao [27], Elishakoff [28], and Elishakoff et al. [29]).For the electrical circuit part, we denote by C the net clamped (i.e., constant strain) capacitance of the piezoelectric material.The electrical load attached to the energy harvester is assumed to be presented by a simple resistor with resistance R; and θ is the electromechanical coupling coefficient.The following boundary-value problem for the electromechanical system is valid: In Equation (1), y(t) is the absolute transverse base displacement, and δ(•) is the delta function in Equation ( 2).Since we consider the clamped-free model, the boundary conditions are (see [18,30]): Note that the base excitation only directly contributes to the transverse dynamics of Equation (1); it does not appear in the rotational dynamics equation of Equation (2).Furthermore, the transverse dynamics is coupled to the rotational dynamics, and the rotational dynamics is coupled to the electrical dynamics, but the electrical dynamics is not directly coupled to the transverse dynamics.This set of pairwise coupled equations is in contrast to the Euler-Bernoulli formulation in which the transverse dynamics is directly coupled to the electrical dynamics (see Dietl [30], Erturk and Inman [31,32], Shubov [1][2][3], Shubov and Shubov [4], and Stoykov et al. [5]).In our first statement, we justify the reformulation of the boundary-value problem in Equations ( 1)-( 5) in the form that does not contain any distributional terms.
Proposition 1.The boundary-value problem given by Equations ( 1)-( 5) is equivalent to the following boundary-value problem: The boundary conditions are: Remark 1.The only difference between Equations ( 1)-( 3) and Equations ( 7)-( 8) is that the distributional term, 2) is replaced by an additional term θv(t) in the right-hand side of the boundary condition in Equation (10).

Operator Reformulation of the Problem
Since we are interested in the asymptotic approximation for the natural frequencies of the system, we place y(t) = 0 (where y(t) is an external force), and arrive at the following system of partial differential equations governing the harvester dynamics: The boundary conditions are given in Equation ( 9) and (10).
Let H be the set of five-component vector-valued functions T obtained as a closure of smooth functions satisfying the conditions ϕ 0 (0) = ϕ 2 (0) = 0 in the following energy norm: If we introduce a new vector-valued function F by the rule: T , then it can be easily checked that the problem in Equations ( 16)-( 18), ( 9) and ( 10) can be written as a first order in time linear system in the energy space, i.e., where the superscript "0" stands for the notation of the initial state.The operator L (the dynamics generator) can be given by the following differential expression: defined on the domain The operation Remark 2. We introduce the imaginary unit "i" into the definition in Equation ( 22) of the dynamics generator and into Equation (20) for convenience.As is shown below, L is a finite-rank perturbation of a selfadjoint operator.Thus, owing to this factor, we deal with a selfadjoint operator rather than with a skew-selfadjoint one.

Properties of the coefficients:
there exist two constants 0 Derivation of the formula for the adjoint operator L * .Lemma 1.The operator L * is defined by the differential expression obtained from Equation (22) in which the second column is replaced with the column 1, 0, 0, 0, θ A L T .This operator is defined on the domain Proof.We have to derive the matrix differential expression and the boundary conditions in such a way that, for any F ∈ D L and G ∈ D L * , the following relation holds: We have for any F ∈ H: Evaluating the inner product induced by norm (Equation ( 19)) yields
Theorem 1.The operator L has the following properties.

(i)
L is an unbounded closed non-selfadjoint operator in H.

(ii)
The inverse operator L −1 exists and is a compact operator in H. Therefore, L has purely discrete spectrum of normal eigenvalues (see Gohberg and Krein [33]).
Proof.The fact that L is non-selfadjoint is due to two reasons: the matrix differential expressions for L and L * are different and the domains D(L) and D(L * ) are different (see Equations ( 23) and ( 26)).
Let us prove that L −1 exists and is compact.To this end, we consider the equation and prove that this equation has a unique solution Rewriting Equation ( 38) component-wise, we obtain the following system: Integrating Equation (39) (iv) and using the condition ψ 0 (L) − ψ 2x (L) = 0, we obtain Substituting Equation (40) into Equation (39) (iii), we get Integrating this equation and taking into account that E(L)ψ 0 (L) = θψ 4 , we obtain Using Equation (39) (i) we derive from Equation (39) (v) that Substituting this result into Equation ( 41), we obtain the following representation which upon integration yields the formula for ψ 0 , i.e., Integrating ψ 2 from Equation ( 40), we derive that Collecting together Equations (39) (i), (39) (ii), and ( 42)-( 44), one can see that the following formula for the solution of Equation ( 38) is valid: As follows from the above proof, the operator L −1 is a bounded operator from H into D(L) if D(L) is equipped with the norm of H 2 .Since the embedding H 2 → H 1 is compact, we get the compactness of the operator L −1 .

Non-Selfadjoint Operator Pencil Generated by the Harvester System
We perform an asymptotic analysis of the fourth order parameter-dependent ordinary differential equation associated to the harvester system (Marcus [34]).We consider the eigenvalue-eigenvector equation for the operator L.
Rewriting Equation (46) component-wise and eliminating ψ 1 and ψ 3 , we get the system of two coupled ordinary differential equations equipped with four boundary conditions: The boundary conditions are Our goal is to eliminate the unknowns ψ 2 and ψ 4 from the system in Equations ( 47) and (48) and from the boundary conditions in Equation ( 49)-(51) and write the pencil eigenvalue problem in terms of ψ 0 .It is technically convenient to deal with new functions f 0 and f 2 defined by Before we present the spectral problem for the pencil, let us complete some preliminary steps.First, we note that Substituting Equation (53) into Equation (48), we get From this equation, we have Now, we turn to Equation (47).Dividing it by ρ(x) and then differentiating, we get a new equation Using the notation in Equation (52), we rewrite Equation (55) in terms of f 0 and f 2 as which can be represented in the form In the sequel, we will work with the following equation: It is convenient to introduce new notations: In these notations, Equations ( 54) and (56) become Substituting Equation (58) into Equation (59), we obtain the desired form of the pencil equation: where Thus, our main equation (Equation ( 60)) is given in terms of f 0 (λ, x) = K(x)ψ 0 (λ, x).Now, we have to rewrite the boundary conditions in terms of f 0 (λ, x).
The boundary conditions for the function f 0 (λ, x).Using Equations ( 49) and ( 52), we get the first left-end condition From Equation (48), we get Substituting Equation (63) into Equation (47), we obtain From Equation (64), we obtain the second condition at x = 0 in terms of f 0 (λ, x) (see Equation (53)): The right-hand side boundary conditions.From now on, we replace Equation (50) with the more general condition in Equation (66) containing a parameter α.This helps us to separate the spectral branches in the future, i.e., we consider the first right-hand side condition as The condition in Equation ( 50) is obtained from Equation (66) by setting α = 0. From Equation (66), we have We use Equation (63) for f 2 and Equation (64) for ψ 2 to have .
Modifying this equation, we get Finally, we rewrite the condition in Equation (51) in terms of f 0 .From Equation (46), it follows that Since ψ 1 = iλψ 0 , we get Now, we take into account the condition in Equation ( 51) and rewrite Equation (68) as Using Equation (53), we modify this equation to the form Hence, our goal is to find the solution of the pencil Equation (60) satisfying the left-end boundary conditions in Equations ( 62) and (65) and the right-end boundary conditions in Equations ( 67) and (69).
Fundamental system of solutions of the pencil equation (Equation (60)).
Let us rewrite Equation (60) in the form of a polynomial with respect to λ ( Marcus [34]) This equation can be written in the form Coefficients w nm (x) can be easily found from Equation (70) (see Naimark [35]).Following the terminology of Fedoruk [36] and Olver [37], the terms of Equation (70) that satisfy the condition n + m = 4 are called the leading terms of Equation (70) and the remaining terms are called the lower order terms.Let us rewrite Equation (60) keeping the leading terms at the left-hand side and moving the lower order terms to the right-hand side.We have where Now, we impose the following assumption on the structural parameters: The next statement is concerned with the fundamental system of solutions of Equation (70).To this end, we perform an asymptotic analysis when the spectral parameter λ tends to infinity.Similar analysis can be found in the works by Shubov [38,39].However, to keep the paper self-contained, we outline the proof of the next statement.
(2) Each function ψ j (λ, x) has four continuous derivatives with respect to x ∈ [0, L] and these derivatives can be approximated as follows when |λ| → ∞: where Proof.First, we consider the equation obtained from Equation (72) in which the right-hand side has been replaced with zero.For the latter equation, we prove that it has four linearly independent solutions with the approximations given precisely by Equations ( 75)-(78).Secondly, we refer to the standard methods of asymptotic analysis (see Fedoruk [36], Miller [40], Olver [37]) to prove that the lower order terms do not destroy the asymptotic behavior of the solutions given by Equations ( 75)-( 78).Thus, we consider the equation where B(x) = Sρ(x)/EI(x) and A(x) = ρ(x)/K(x).Let us reduce this equation to the first order linear system, i.e., we consider the linearization of Equation ( 82) with respect to λ (see Marcus [34] and Fedoruk [36]).To this end, we make the substitutions Equation ( 82) can be rewritten in the following form: where In the next step, we reduce the system in Equation (84) to the asymptotically diagonal form.To this end, we have to calculate the eigenvalues and the eigenvectors of the matrix A(x).For the eigenvalues, we have the following equation: Due to Equation (74), Equation (86) has four different roots given by the formulas The eigenvector F j (x), corresponding to the eigenvalue λ j (x), has the form Let T 0 (x) be the following Vandermonde matrix: The matrix T 0 (x) reduces A(x) to the diagonal form, i.e., One can show by a straightforward calculation that det T 0 (x) = −4 A(x)B(x) A(x) − B(x) 2 . (91) Since due to Equation (74), A(x) = B(x), there exists T −1 0 (x).Let us look for a solution of Equation ( 84) in the form Y(x) = T 0 (x)Z(x).
Rewriting Equation (84) with respect to Z, we have Applying matrix T −1 0 (x) from the left, we obtain from Equation (92) the following equation for Z(x): with Λ 0 (x) being introduced in Equation (90).It is convenient to look for Z in the following form [36]: where 4 × 4 matrix T 1 (x) is chosen below.Equation ( 93) can be rewritten for W as If we assume that T 1 (x) and its first derivative exist and are bounded as functions of x ∈ [−L, 0], then for large |λ|, the matrix of the system in Equation (95) has the following asymptotical form: where Λ 0 (x), T 1 (x) is the commutator of Λ 0 (x) and T 1 (x).Now, we are in a position to choose the matrix T 1 (x).We do this in such a way that the matrix B(x) defined by is a diagonal matrix.Let us fix the choice of T 1 (x) by the following conditions: where λ j (x) 4 j=1 are given in Equations (87).Note, since A(x) = B(x), the matrix T 1 (x) from Equation (98) is well-defined.It can be verified by a straightforward calculation that B(x), defined by Equation (97), is diagonal if T 1 (x) is defined by Equation (98).
We use the notation Λ 1 (x) for the matrix B(x) with the aforementioned choice of T 1 (x).With this choice of T 1 , Equation (95) for W can be modified to the following equation: where the estimate O λ −1 is uniform with respect to x ∈ [0, L] and Λ 1 (x) is a diagonal matrix whose diagonal elements are equal to: Using standard methods of asymptotic analysis (such as in Fedoruk [36], Miller [40], Olver [37]), we obtain that the solutions W j (λ, x) of the system in Equation (99) are asymptotically close to the solutions W 0 j (λ, x) of the following system: (Explicit formulas for W 0 j (λ, x), j = 1, 2, 3, 4, are given in Equation (104) below.)Namely, the following relations are valid: and the estimates O λ −1 are uniform with respect to x ∈ [0, L].To complete the proof, we have to calculate the diagonal elements of the matrix T −1 0 (x) dT 0 (x) dx .Straightforward calculations yield the following results: Substituting Equation (103) into the system in Equation ( 101) and taking into account Equation (90), we immediately obtain the following expression for the k-th component of the four-vector W 0 j (λ, x): Using Equations ( 102), ( 104) and (105), we immediately obtain the asymptotic approximations for the solutions of Equation (99).Since Y j (λ, x) = T 0 (x)Z(λ, x) and Z(λ, x) is given by Equation (94), we obtain Taking into account the explicit Equation (89) for T 0 (x), we obtain that the first components of the vectors Y j (λ, x), j = 1, 2, 3, 4 are given exactly by the formulas from the second line of (105).Statement 1 of Theorem 2 is completely shown.
The results of Statement 2 are almost straightforward consequences of Equations ( 83) and the explicit form of the matrix T 0 (x).

Asymptotics of the Spectrum of the Operator L. Reflection Matrices at the Ends of the Beam
In this section, we derive asymptotic approximations for two special matrices that we call the leftreflection matrix and the right-reflection matrix.We apply the method of reflection matrices developed in a series of our earlier papers: Shubov and Peterson [41], Shubov [25,38,39].This method is very efficient in completing mathematical calculations which are extremely lengthy when using direct approach.Using these matrices, we obtain two equations (see Equations ( 157) and (158) below) generating two branches of the eigenvalues.Since both the operator L and the pencil P (λ) have the same spectra, we will derive the aforementioned equations for the pencil.To this end, we look for a solution of Equation (70) in the form of a linear combination of the fundamental set of solutions in Equations ( 75)-( 78) where C j (λ) , and D j (λ), j = 1, 2 are arbitrary functions of λ.Note that physically, in the first sum of Equation ( 106), there are the Fourier representations for the terms corresponding to the waves moving from right to left and in the second sum-the Fourier representations for the terms corresponding to the waves moving from left to right (see Equations ( 75)-( 78)).
Proof.We look for a solution Ψ(λ, x), which satisfies two boundary conditions in Equations ( 62) and (65).Substituting Equation (106) into the boundary conditions in Equations ( 62) and (65) related to the end x = 0, we obtain two equations for four unknown coefficients C j (λ), and D j (λ), j = 1, 2. It is convenient to represent the asymptotical approximations for the derivatives of the function Ψ(λ, x) Using Equation ( 106) we obtain that the first boundary condition in Equation ( 62) at x = 0 has the form The asymptotical form for the second boundary condition in Equation ( 65) at x = 0 is Using Equations ( 108) and ( 110) we rewrite Equation (112) as follows: Simplifying this equation and collecting like terms, we have We rewrite this equation in the desired form as Taking into account the assumption in Equation ( 74) SK(0) = EI(0) , we combine Equations ( 111) and (113) in a system and have This system can be written in a matrix form as Taking the inverse of the left-hand side matrix and applying it to the right-hand side, we finally obtain the desired result (Equation ( 108)).
Using the boundary conditions at the right end of the beam, we derive the representation for the right-reflection matrix R r (λ).We need the notations with η(L) and η(L) being given in Equation (79).The following result holds.
Assume that the parameter α satisfies the following condition: For the solution Ψ(λ, x) of Equation ( 106) to satisfy the right-end boundary conditions in Equations ( 67) and ( 69), the coefficients C j (λ), D j (λ), j = 1, 2, must be related by means of the right-reflection matrix R r (λ): with R r (λ) being given by The following asymptotic approximations are valid for the entries d jk (λ), j = 1, 2, k = 1, 2, as |λ| → ∞: Proof.Let us represent the condition in Equation (67) in the asymptotical form: we modify asymptotic Equations ( 108)-( 110) taken for x = L to the form Substituting the approximations in Equations ( 125)-(127) into Equation (123), we obtain the asymptotical form of the condition in Equation ( 67) Collecting coefficients for C 1 (λ)e(λ) and using Equation (124), we get Collecting coefficients for C 2 (λ) e(λ) and using Equation (124), we get Collecting coefficients for D 1 (λ) e(λ) −1 , we get Collecting coefficients for D 2 (λ) e(λ) −1 , we get Collecting together Equations ( 129)-(132), we obtain the following asymptotic approximation for the first right-end boundary condition: Now, we turn to the condition in Equation (69) that can be given in the form Using Equations ( 75)-( 78) and (108), we represent this equation in the asymptotical form as Collecting coefficients for C 1 (λ)e(λ), we get Collecting coefficients for C 2 (λ) e(λ), we get Collecting coefficients for D 1 (λ)e −1 (λ), we get Collecting coefficients for D 2 (λ) e −1 (λ), we get Collecting together Equations ( 136)-(139) and using Equation (118), we obtain the following asymptotic approximation for the second right end boundary condition: Let Using Notations (141), we rewrite the right-end boundary conditions in Equations ( 133) and (140) as the following system: The system can be written in a matrix form as where E(λ) is defined in Equation (117).Let A(λ) be the matrix at the left-hand side of Equation ( 143) and B(λ) be the matrix at the right hand-side; then, Equation (143) becomes Now, we derive the formulas for the entries of A −1 (λ).We have Taking into account that Θ(λ) which means that the matrix A −1 (λ) exists.Direct calculations yield the following results for the entries of A −1 (λ) = A ij 2 i,j=1 : Let D(λ) be the notation for the product , then the following formulas hold: Thus, we finally obtain Taking into account Equations ( 107) and (141), we arrive at the desired result of the theorem.

The Spectral Asymptotics
In this Section, we derive the asymptotic distribution of the eigenvalues as the number of an eigenvalue tends to infinity.To this end, we derive the spectral equation whose solutions coincide with the problem eigenvalues.Let us represent Equations ( 107) and (119) in the forms where From Equation (149), it follows that This equation has nontrivial solutions if and only if where I is the identity matrix.In the sequel, we need the following result.
Lemma 2. The set of roots of Equation ( 152) is located in a strip parallel to the real axis.
Proof.Notice Equation (152) has the same set of roots as the following equation: where E(λ) is defined in Equation (117).Using the contradiction argument, assume that there exists a sequence of roots of Equation ( 153) denoted by λ m 153) for λ = λ m , we obtain (see Equation ( 119)): It can be readily checked that, when y m → ∞, we have e λ m → 0 and e λ m → 0. Hence, Equation (154) can be written in the form which means that the left-hand side approaches its limit: while the right-hand side approaches 0 as m → ∞.This contradiction means that the set of all roots that belongs to the upper half-plane can be located in a strip parallel to the real axis.Using a similar argument, we can show that the set of all roots that belongs to the lower half-plane can be located in a strip parallel to the real axis as well.
In the following, we consider Equation (153) in a strip on the complex plane parallel to the real axis.Let us multiply matrices D(λ) of ( 148) and E −1 (λ)R l (λ)E −1 (λ).As the result, we reduce the spectral equation (Equation (153)) to the form Due to Lemma 2, when λ belongs to the strip parallel to the real axis, we have e(λ) 1 and e(λ) 1; therefore, the asymptotical form of the spectral equation becomes Using Equations ( 117) for e(λ) and e(λ), and Equations (147) for d 11 (λ) and d 22 (λ), we can readily check that Equation (156) splits up into the following two equations: Each equation has a countable set of solutions.The solutions of Equation ( 157) we call the α-spectral branch and the solutions of Equation ( 158) we call the θ-spectral branch.
In the statement below (Theorem 5), we provide the formulation of the results on the asymptotical distribution of the roots of Equations ( 157) and (158).We do not present the proof of the theorem since very detailed proofs of similar results can be found in Shubov [2,3], Shubov and Shubov [4], Shubov and Peterson [41], and Shubov [25,38].

Theorem 5. (1)
The following asymptotic approximation is valid for the α-branch of the spectrum: (2) The following asymptotic approximation is valid for the θ-branch of the spectrum: Conclusion.If the harvester model is based on the cantilever beam as a substructure, then the leading asymptotical terms are not affected by the piezoelectric parameters.For the Timoshenko beam, this conclusion follows from Equations ( 159) and (160).For the Euler-Bernoulli beam, it follows from the asymptotic formulas derived in Shubov [1][2][3] and Shubov and Shubov [4].

The Harvester Model with Constant Physical Parameters. The Left Reflection Matrix
As shown in Sections 1-6, for the harvester model with the Timoshenko beam as a substructure, the asymptotic distribution of the vibrational modes has the following important property: the leading asymptotical term does not depend on the properties of the piezoceramic patches.It means that the perturbation induced by these patches is, in some sense, small.At the same time, the entire effect of the energy harvesting takes place exactly due to the physical properties of the piezoelectricity.Therefore, it is important to figure out the level of accuracy in the asymptotical formulas for the spectrum so that the contribution of the piezoelectric patches physical parameters will be revealed.To this end, we consider the same model of the harvester based on the Timoshenko beam substructure with constant (not variable) physical parameters.We derive refined asymptotic approximations for the vibrational modes and show that the second asymptotical term (next to the leading one) contains contribution from the electrical circuit of the model.Similar results can be obtained for the case of the variable structural parameters of the beam.However, preliminary calculations show that the derivation of the formulas is extremely lengthy and the result does not provide any additional insight into the physics of the problem.
We consider the following system of equations equipped with the boundary conditions: In this section, we use the standard engineering approach and look for the solutions of the boundary-value problem in the form W(x, t) = e iλt w(λ, x), Φ(x, t) = e iλt ϕ(λ, x).
Substituting these representations into Equations ( 161) and ( 162), we obtain a system for w and ϕ To eliminate w from the system and write a new equation in terms of ϕ, we have to complete some preliminary steps.First, we derive from Equation (167) that Then, we differentiate Equation (166) once and use Equation (168) to get −ρλ 2 w (λ, x) + K ϕ (λ, x) − w (λ, x) = 0, which yields the desired form of the equation for ϕ The corresponding characteristic equation is The roots of the biquadratic equation are At this moment, it is convenient to introduce new notations: and represent the roots of Equation (171) in the form For definiteness, we assume A > B. Then, the following asymptotic approximations are valid for the roots of Equation (171): Rewriting Equation (174) in terms of the physical parameters of the structure, we get Let Λ(λ) and Λ(λ) be defined by The general solution of Equation ( 169) can be given as the following linear combination of the exponentials: where C j (λ) and D j (λ), j = 1, 2, are arbitrary functions of λ.Now, we establish two constraints on the coefficients in order to satisfy the left-end boundary conditions in Equation (163).First, we rewrite these conditions ϕ(λ, 0) = w(λ, 0) = 0 in terms of ϕ.
To this end, we use Equation (166) and have and also from Equation (167) we have Differentiating this equation once and substituting the result for ϕ − w into the formula for w, we obtain Therefore, the second boundary condition at the left end becomes: Now, use the general solution of Equation ( 169) in the form of Equation ( 179) and obtain that, for the solution to satisfy the left-end boundary conditions, the coefficients must satisfy the following linear system: The second equation can be modified as Introducing a new notation we reduce the system of Equations (182) to the following form: By subtracting the first equation from the second one we immediately obtain that C 2 (λ) can be represented as a linear combination of the coefficients D 1 (λ) and D 2 (λ): Rewriting Equations (185) in the form we obtain C 1 (λ) as a linear combination of the coefficients D 1 (λ) and D 2 (λ) Rewriting Equations ( 186) and (188) as as one matrix equation, we obtain The matrix at the right of this equation we will call the left reflection matrix and denote it by R l (λ).

The Right Reflection Matrix. The Spectral Asymptotics
In this section, we derive the constraints on the coefficients C i (λ) and D i (λ), i = 1, 2, that follow from the fact that the general solution in Equation (179) has to satisfy two right-hand side boundary conditions.As we already know, the boundary conditions at the right end of the beam are To rewrite the first condition in terms of ϕ, we use Equations (180) for w and (166) for w to get the equation Using Equations ( 177) and (178), we derive the general solution derivatives: Substituting Equations (206) into the condition in Equation (205), we obtain It is convenient to introduce new notations.Namely, let Rewriting Equation (207) in terms of Equation (208), we have Now, we turn to the second boundary condition in Equation (204).Let Then, Equation (204) becomes The second boundary condition generates the second relation between Let us introduce two two-vectors: Equations ( 209) and (212) can be written in a matrix form as where and E(λ) is a diagonal matrix: E(λ) = diag e(λ), e(λ) .The matrix that relates vectors C(λ) and D(λ) of Equation (213), given by R r (λ) = E −1 (λ)A −1 (λ)B(λ)E −1 (λ), is called the right reflection matrix.The next statement is the main result of this Section.

Theorem 6.
(1) The initial-boundary value problem in Equations ( 161)-(164) has a countable set of complex vibrational modes, which belongs to a strip parallel to the real axis and can accumulate only at infinity.There may be only a finite number of multiple modes of a finite algebraic multiplicity each.
(2) The entire set of modes asymptotically splits into two disjoint subsets: we call them the α-branch and the θ-branch and denote them by λ α n n∈Z and λ θ m m∈Z , Z = Z 0 .If α ≥ 0, then the α-branch is asymptotically close to a horizontal line in the upper half-plane and the second branch is asymptotically close to the real axis.

Conculsions
In the paper, a mathematical model of a piezoelectric energy harvester with a Timoshenko beam as a substructure is considered.The model is represented in the form of three coupled partial differential equations equipped with the boundary conditions at the left and right ends of the flexible structure.Two cases are considered: (i) the first case is concerned with the beam with spatially variable structural parameters; and (ii) the second case is concerned with the beam with constant structural parameters.In both cases the asymptotic distributions of the natural frequencies have been derived for electrically loaded harvesters.In Case (i). the leading asymptotical term and the remainder term have been obtained.The derivation of the higher order asymptotical terms for a variable parameter model becomes technically very complicated and lengthy.In Case (ii), the two-term asymptotical approximation for the natural frequencies of the model has been obtained.It is shown that the leading asymptotical terms do not contain the piezoelectric parameters, while the second order asympotical term contains the piezoelectric coupling coefficient.