1. 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 two-degree-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: lumped-parameter 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,
Section 7 and
Section 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.
2. 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: : the relative transverse deflection of the beam with respect to its base; : the rotation angle of the beam cross-section; and : 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:
is the flexural rigidity of a beam,
is the mass density,
S is the cross-sectional area of beam, and
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),
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, , in Equation (2) is replaced by an additional term in the right-hand side of the boundary condition in Equation (10). Proof (Proof of Proposition 1). To show that the two problems are equivalent, we check that their weak (variational) formulations coincide. Consider the following class of functions:
Multiplying Equation (
1) by
and integrating, we have
Taking into account that
, and the first condition of Equation (
10), we obtain that
Now, we multiply Equation (
2) by
and integrate to have
Taking into account that
and
, we reduce Equation (
13) to
Thus, Equations (
3), (
12) and (
14) constitute the weak formulation of the problem given by Equations (
1)–(
3). To obtain the weak formulation of the problem given by Equations (
7)–(
8), we multiply Equation (
7) by
and integrate, which immediately yields Equation (
13). Finally, we multiply Equation (
10) by
and integrate to have
If we take into account that
and the second boundary condition from Equation (
10), we immediately reduce Equation (
15) to the form of Equation (
14). The fact that weak formulation in Equations (
14), (
15), (
4), and (
5) is the same for both problems in Equation (
1)–(
3), (
4), (
5), and (
7)–(
8) with Equations (
9) and (
10) implies that these two problems are equivalent. ☐
4. 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
. Let
Rewriting Equation (
46) component-wise and eliminating
and
, 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
and
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
. It is technically convenient to deal with new functions
and
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
and then differentiating, we get a new equation
Using the notation in Equation (
52), we rewrite Equation (
55) in terms of
and
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
. Now, we have to rewrite the boundary conditions in terms of
.
The boundary conditions for the function . 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
in terms of
(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
. From Equation (
66), we have
We use Equation (
63) for
and Equation (
64) for
to have
Modifying this equation, we get
Finally, we rewrite the condition in Equation (
51) in terms of
. From Equation (
46), it follows that
Since
, 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
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
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.
Theorem 2. (1) Under the assumptions in Equations (24), (25) and (74), Equation (72) has four linearly independent solutions with the following asymptotic approximations as :where All estimates in (75)–(78) are uniform with respect to . The functions , are analytic functions of λ. (2) Each function has four continuous derivatives with respect to and these derivatives can be approximated as follows when :where All estimates in Equation (80) are uniform with respect to . 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
and
. 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
. For the eigenvalues, we have the following equation:
Due to Equation (
74), Equation (
86) has four different roots given by the formulas
The eigenvector
, corresponding to the eigenvalue
, has the form
Let
be the following Vandermonde matrix:
The matrix
reduces
to the diagonal form, i.e.,
One can show by a straightforward calculation that
Since due to Equation (
74),
, there exists
. Let us look for a solution of Equation (
84) in the form
Rewriting Equation (
84) with respect to
Z, we have
Applying matrix
from the left, we obtain from Equation (
92) the following equation for
:
with
being introduced in Equation (
90). It is convenient to look for
Z in the following form [
36]:
where
matrix
is chosen below. Equation (
93) can be rewritten for
W as
If we assume that
and its first derivative exist and are bounded as functions of
, then for large |
|, the matrix of the system in Equation (
95) has the following asymptotical form:
where
is the commutator of
and
. Now, we are in a position to choose the matrix
. We do this in such a way that the matrix
defined by
is a diagonal matrix. Let us fix the choice of
by the following conditions:
where
are given in Equations (
87). Note, since
, the matrix
from Equation (
98) is well-defined. It can be verified by a straightforward calculation that
, defined by Equation (
97), is diagonal if
is defined by Equation (
98).
We use the notation
for the matrix
with the aforementioned choice of
. With this choice of
, Equation (
95) for
W can be modified to the following equation:
where the estimate
is uniform with respect to
and
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
of the system in Equation (
99) are asymptotically close to the solutions
of the following system:
(Explicit formulas for
are given in Equation (
104) below.) Namely, the following relations are valid:
and the estimates
are uniform with respect to
. To complete the proof, we have to calculate the diagonal elements of the matrix
. 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
:
Using Equations (
102), (
104) and (
105), we immediately obtain the asymptotic approximations for the solutions of Equation (
99). Since
and
is given by Equation (
94), we obtain
Taking into account the explicit Equation (
89) for
, we obtain that the first components of the vectors
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
. ☐
7. The Harvester Model with Constant Physical Parameters. The Left Reflection Matrix
As shown in
Section 1,
Section 2,
Section 3,
Section 4,
Section 5 and
Section 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
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
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
. 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
and
,
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
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
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
can be represented as a linear combination of the coefficients
and
:
Rewriting Equations (
185) in the form
we obtain
as a linear combination of the coefficients
and
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 .
Asymptotic approximations for the entries of as .
In the sequel, we need the asymptotic approximation for
as
. Using Equations (
172), (
175) and (
177), we obtain
Using Equations (
172), (
176) and (
178), we obtain
Equations (
190) and (
191) yield the following result:
Now, we evaluate
Q defined by
First, we use Equation (
178) and evaluate
and have
Next, we use Equation (
177) and evaluate
and have
Using Equations (
195) and (
196), we obtain the following asymptotic approximation for
Q:
Finally, using Equations (
184), (
193) and (
197), we obtain the approximation for
:
Let
then
and the entries of the left reflection matrix
can be represented as follows when
:
Substituting Equations (
201) (i)–(iii) into the left reflection matrix of Equation (
189), we obtain