Nonlinear Dynamic Analysis of a Piezoelectric Energy Harvester with Mechanical Plucking Mechanism

In this study, we propose an analytical approach based on the modified differential transform method to investigate the dynamic behavior of a plucking energy harvester. The harvester consists of a piezoelectric cantilever oscillator and a rotating plectrum. The analytical approach provides a closed-form solution that helps determine the starting and ending points of the contact phase between the piezoelectric cantilever and the plectrum. This analytical approach is valuable for simulating complex dynamic interferences in multiple or periodic plucking processes. To evaluate the effects of plucking speed and overlap length of the plectrum on single and periodic plucking, a series of simulations were carried out. The output voltage of the piezoelectric energy harvester increases as the overlap length of the plectrum increases. On the other hand, increasing the plucking speed tends to amplify the magnitude of the contact force while reducing the duration of the contact phase. Therefore, it is crucial to optimize the plucking speed to achieve the maximum linear impulse. For periodic plucking, successful synchronization between the motions of the piezoelectric energy harvester and the rotating plectrum must occur within a limited contact zone. Otherwise, dynamic interferences often cause the plectrum to fail to pluck the energy harvester exactly within the contact zone. Additionally, reducing the plucking speed of the plectrum and increasing the overlap length would be more advantageous for successful periodic-plucking energy harvesting.


Introduction
Energy harvesting involves a series of processes aimed at extracting ambient energy, converting it into electrical energy, and utilizing or storing it to power micro-sensors or devices. Various ambient energy sources are available for energy harvesting, including sunlight, wind, heat, and vibration. Among these sources, mechanical vibration is commonly present in our daily lives. The energy conversion mechanisms for mechanical vibration include piezoelectric, electrostatic, and electromagnetic transductions [1][2][3], which are relatively straightforward compared to other energy sources. As a result, vibration energy harvesters  have been developed for diverse applications, including small-scale portable devices, MEMS sensors for the Internet of Things, and autonomous wireless sensor networks for smart environments.
In the early stages of development, vibration energy harvesters were primarily designed using linear electromechanical oscillators, which had a limited resonant frequency range. These linear resonant-type harvesters were only capable of generating high energy within their narrow resonant frequency band. Consequently, they proved to be unsuitable for complex ambient vibration sources that involved multiple frequency components and/or time-varying frequency components [4,5].
Extensive research efforts have been dedicated to expanding the frequency bandwidth of vibration energy harvesters. To achieve this, researchers have proposed and developed various energy harvesting structures, employing different technologies such as multi-modal oscillators or oscillator arrays [6][7][8], frequency-tunable oscillators [9,10], nonlinear monostable oscillators [11][12][13], and multi-stable oscillators [14][15][16][17][18][19][20][21][22][23][24]. However, the utilization of multi-modal oscillators or oscillator arrays often resulted in larger and more complex structures for energy harvesting. Nonlinear monostable or multi-stable oscillators showed promise for broadband energy harvesting, but the presence of hysteresis due to multiple solutions often resulted in decreased energy harvesting efficiency. Conversely, an alternative approach known as frequency-up conversion has emerged for harvesting energy from low-frequency ambient vibrations [25][26][27][28][29][30]. This technique involves utilizing mechanical/magnetic impact or a plucking force to stimulate the high-frequency natural vibration of the energy-harvesting oscillator, which can then be converted into electrical energy.
Mechanical plucking energy harvesters have been developed to extract energy from the periodic motions of rotor blades or rotating plectra in various applications, such as wind turbines, human knee joints, and automobile wheels [18][19][20][21][22][23]. The fundamental design of mechanical plucking energy harvesters consists of an energy-harvesting oscillator and an external plectrum. When the plectrum impacts the energy-harvesting oscillator that has a relatively high natural frequency, it imparts a linear impulse to the oscillator. Subsequently, the oscillator undergoes free vibrations, which are then converted into electrical energy. The main challenges in designing such mechanical plucking energy harvesters arise from the contact mechanism between the energy-harvesting oscillator and the external plectrum, involving complex dynamic interactions [31][32][33][34][35][36][37]. Mathematical models based on Hertzian contact theory [36][37][38][39] and methodologies utilizing finite element analysis or other numerical methods have been employed to predict the complex dynamic behaviors of the mechanical plucking energy harvesters, including strong dynamic interferences and multiple or periodic plucking effects.
In this study, we propose an analytical approach based on the modified differential transform method to investigate the dynamic behaviors of a mechanical plucking energy harvester, which comprises a piezoelectric cantilever oscillator and a rotating plectrum. The analytical approach provides a closed-form solution that facilitates the determination of the starting and ending points of the contact phase between the piezoelectric cantilever and the plectrum. This analytical approach is valuable for simulating complex dynamic interferences in multiple or periodic plucking processes. We conduct a series of simulations to evaluate the effects of plucking speed and overlap length of the plectrum on single or periodic plucking, and we discuss the results. Figure 1 shows the schematic diagram of a piezoelectric energy harvester with a mechanical plucking mechanism. The harvester consists of a piezoelectric bimorph cantilever beam and a rotating plectrum. The bimorph cantilever beam comprises a stainless-steel substrate partially covered by two identical piezoelectric films which are polarized oppositely in the vertical direction. The natural frequency of the bimorph cantilever beam oscillator can be adjusted by altering the thickness of the stainless-steel substrate while keeping the dimensions of the piezoelectric layers unchanged. Both major surfaces of each piezoelectric layer are fully coated with conductive electrodes. These two piezoelectric lay-ers are connected in series to an external load resistance R. Table 1 provides the geometric parameters and material properties of a typical piezoelectric bimorph cantilever.  When the rotating plectrum plucks the tip of the bimorph cantilever, the cantilever structure is deflected, and the resulting mechanical strain in the piezoelectric materials is converted into the electrical current that flows through the external load resistance. Through this energy conversion mechanism, the piezoelectric cantilever energy harvester can extract useful electrical energy from rotational motions, such as individuals passing through turnstile gates, rotating vehicle wheels, and wind turbines, to power autonomous micro-sensors or devices.

Electromechanical Duffing Oscillator Model of the Mechanical Plucking Energy Harvester
The piezoelectric bimorph cantilever beam is modeled based on the Euler-Bernoulli beam theory and linear piezoelectricity. The Kelvin-Voigt damping model is employed to account for the material damping effect of the cantilever beam. Firstly, the governing field equations and boundary conditions are derived using D'Alembert's principle and the moment balance method [13]. Following a series of derivation processes, including modal analysis, discretization process, and single-mode approximation [13][14][15], the differential equations of motion for the mechanical plucking energy harvester and the initial conditions can be obtained as follows: ..
where the dot over the symbols indicates the time derivative; w is the tip displacement of the cantilever beam; ζ and ω n are the natural angular frequency and damping ratio, respectively; β is the coefficient of the cubic nonlinear term; θ is the electromechanical coupling constant; V is the output voltage; R and C p are the external load resistance and equivalent capacitance; and α f p is the generalized plucking force. In Equation (1), the cubic nonlinear term is added to describe some nonlinear effects such as nonlinear magnetic interactions [39] and structural nonlinearity [13]. The derivation process for the equivalent circuit equation, given by Equation (2), is provided in Appendix A. Additionally, to validate the oscillator model, a comparison is made between the numerical results obtained using Equations (1) and (2) and the experimental results reported in Ref. [37], which are provided in Appendix B. For generality, introducing the following dimensionless variables and parameters to Equations (1)-(3), the dimensionless electromechanical oscillator model for the plucking energy harvester and the initial conditions are obtained in the form: .. .
Note that the dimensionless mathematical model is used for all the dynamic simulations conducted in this study. Henceforth, unless specified otherwise, the variables and parameters will be treated as dimensionless.
In general, the plucking process consists of two distinct phases: the contact phase, during which the tips of the bimorph cantilever beam and the plectrum remain in contact, and the release phase, in which the beam undergoes free vibration after the separation from the plectrum. Although the contact phase has a much shorter duration compared to the release phase, it holds greater significance as it provides the kinetic energy required to initiate the subsequent free vibration of the energy harvester. The contact phase involves complex dynamic interactions between the beam and plectrum tips, encompassing sliding, shock, and vibration [38]. Furthermore, it is well-known that when the plucking speed is small (quasi-static plucking process), excessive surface friction significantly impacts the dynamic interaction, giving rise to complicated effects such as the stick-slip phenomenon. In this study, we focus on relatively high plucking speeds (dynamic plucking process), assuming the negligible influence of surface friction during the contact phase.
In the dynamic plucking process, the contact period is very short, and no plastic deformation occurs. This implies that the contact force can be considered as a spring force, dependent on local elastic deformations. Under these assumptions, the plucking force, denoted as f in Equation (6), is modeled based on the Hertzian contact theory [37,39]. Additionally, assuming a linear relationship between the contact force and the indentation, the plucking force can be expressed as follows: where κ is the linear contact stiffness; δ is the indentation; u is the vertical displacement of the plectrum; R is the radius of curvature of the circular motion of the plectrum; and ω and φ 0 are the angular frequency and phase angle, respectively. In Equation (9), δ ≥ 0 for the contact phase, whereas δ < 0 for the release phase after the loss of contact. The mechanical linear impulse, which is the integral of the plucking force over the time interval of the contact phase, determines the change in the linear momentum of the cantilever beam and initiates its free vibration. This linear impulse depends on both the duration of the contact and the contact area. The contact duration, denoted as ∆t, is measured from the initial contact time, t 0 , to the separation time, t s , and by setting t 0 = 0, it becomes equivalent to the separation time itself, i.e., ∆t = t s . If the overlap length is very small, the geometric separation distance can be approximated as equal to the overlap length between the cantilever beam and the plectrum at a trivial position, as shown in Figure 1. Therefore, the contact duration or the separation time can be evaluated as follows: where w d is the overlap length and v p is the plucking speed.
In previous studies, Equation (10) was used to determine the end of the contact phase in the plucking process. However, it is important to note that the separation time obtained from the geometric configuration may or may not coincide with the actual end of the contact phase. This is because contact loss can occur within the separation distance due to excessive relative motion between the cantilever beam and the plectrum. In this study, for all simulations conducted, it is determined whether contact loss occurs within the separation distance. To accomplish this, the analytical solution of the initial value problem described by Equations (6)-(8) is obtained and used to determine the condition of zero indentation, which represents the actual end of the contact phase prior to the separation distance.

Differential Transform Method
In this study, we utilize an analytical technique known as the differential transformation method (DTM) to derive a closed-form approximate solution of the nonlinear electromechanical model presented in the previous section for the mechanically plucked piezoelectric energy harvester.
The differential transformation of an analytical function The differential inverse transformation F(k) → f (t) can be performed as follows: As presented in Equations (11) and (12), the transformed function F(k) is a sort of the coefficient of power series of the original function f (t). The fundamental operations of the DTM are provided in Table 2. Refer to References [40,41] for the detailed derivation of the fundamental DTM operations. Table 2. Fundamental operations of the differential transform method.

Original Function
Transformed Function To solve the initial value problem of the Duffing oscillator forced by mechanical plucking, we apply the differential transformation method (DTM), at t = 0, to the initial conditions and the system of ordinary differential equations given in Equations (6)- (8), which yields the transformed initial conditions and the transformed system of equations where As shown in Equations (14) and (15), the DTM transforms the original differential equations into a system of algebraic equations involving the unknowns X(k) and Y(k), which are the transformed functions of displacement x(t) and voltage y(t), respectively. Specifically, Equations (14) and (15) form a recursive system of algebraic equations that allow us to determine the solution by substituting the values of X(k), X(k + 1), and Y(k) into the equations. This process begins with the transformed initial conditions, for k = 0, given in Equation (13).
Finally, using the differential inverse transformation, we obtain the approximate solution functions x(t) and y(t) for the initial-value problem in the following forms of truncated Maclaurin series up to the Nth order: Note that the solution functions are locally analytic at t = 0, allowing them to be approximated by truncated power series centered around the initial point. The approximate solutions obtained from the DTM demonstrate improved convergence to the exact solution with reduced error, but this convergence is guaranteed only in the vicinity of the initial point. To expand the region of convergence beyond this limited range, we will apply the Laplace-Padé resummation method to the truncated series solutions in the next section.

Laplace-Padé Resummation Method
By applying the Laplace transform with respect to t to the truncated power series solutions given in Equation (17), we obtain the transformed solution functions in terms of s: where L denote the Laplace transform. The Padé approximants for the solution functions given in Equation (18) are then calculated. These approximants are known to provide the best approximation using a rational function of a specific order [42,43]. The Padé approximation is applied to obtain the solutions for displacement and voltage, which are given as follows: n=0 a n s −n ≈ 0, (19) Sensors 2023, 23, 5978 (20) where L and M are the upper bounds of summations, b 0 is set to be zero for normalization, a n and c n are the L + 1 numerator coefficients of the displacement and voltage solutions, respectively, and b n and d n are the M denominator coefficients. In Equations (19) and (20), by equating the coefficients of like powers of s, we obtain a linear system of L + M + 1 algebraic equations for the unknowns a n and b n , which can be solved. Finally, we apply the inverse Laplace transformation with respect to s to obtain the approximate solutions in the following forms: Note that the numbers of the coefficients, L and M, can be arbitrarily chosen within the range of L + M + 1 ≤ N + 1. It is known that the error of the Padé approximation is minimized when L = M + 1 or L = M for a fixed value of L + M + 1. The Padé approximants mentioned above help expand or improve the region of convergence of the series solutions and provide better approximations compared to the truncated Taylor series solutions (Equation (17)) which have limited convergence and are only accurate in the vicinity of the initial point.

Numerical Simulation
In the previous section, the DTM was enhanced by incorporating the Laplace-Padé resummation method to improve the accuracy of the approximate solution functions. These modified DTM (MDTM) solutions are then compared with numerical results for validation purposes. For this purpose, a series of numerical simulations are conducted.
By introducing the state vector x = [x 1 , x 2 , x, x, y into Equations (6)-(8), the governing differential equations are reformulated into a system of first-order differential equations in the following vector form: The transient response of the piezoelectric energy harvester is obtained by numerically integrating Equation (22)   x(0) = 0, and y(0) = 0, ensuring that the oscillator is at rest before being plucked by a rotating plectrum. The initial angular position of the rotating plectrum is specified as φ 0 = 2π/3. As the plectrum reaches a trivial position (i.e., u = x = 0), it makes initial contact with the tip of the piezoelectric cantilever.  The plucking process of the piezoelectric energy harvester involves two distinct time scales of dynamic motion: (i) the contact phase between the cantilever and plectrum, which occurs for a very short duration, and (ii) the subsequent free vibration process of the piezoelectric cantilever, which generates electrical energy. Therefore, the MDTM solution procedures need to be applied separately to these two different time-scale motions.
The comparison between the MDTM solutions and the RK results in Figure 2 demonstrates strong agreement and confirms the validity of the analytic technique employed in this study. The analytical MDTM solution offers a simpler means of determining the duration of the contact phase compared to numerical methods. In Figure 2b, the starting and ending points of the contact phase correspond to the plucking and separation points, respectively. Since the piezoelectric cantilever is initially at rest, the plucking time t p can be determined by satisfying the condition of zero indentation at a trivial position, expressed as follows: After the initial plucking event, the separation point can be determined by satisfying the conditions of zero indentation within the geometric contact zone (0 < u < w p ). This can be expressed as follows: where τ is the time defined as τ = t − t p that is introduced for illustration purposes, so that the starting and ending points of the contact phase become at τ = 0 and τ = τ s , respectively. In Figure 3a, the contact force between the piezoelectric cantilever and plectrum during the contact phase is depicted for different plucking speeds: 15, 24, 32, and 40. The ending points of the contact phase are indicated by the open square points S1-S4 corresponding to the respective plucking speeds. As observed in Figure 3a, the duration of the contact phase, denoted as ∆τ(= ∆t), or the separation time, τ s , exhibits significant dependence on the plucking speed. When the plucking speed is relatively high, the contact phase duration or separation time is short, as the geometrical separation point is fixed at u = w d . Consequently, as the plucking speed decreases, the separation time tends to increase, as evidenced by the open square points S1-S3 in Figure 3a,b.
However, this trend is interrupted at a critical point, indicated by the solid red circle in Figure 3b. This occurs because, at low plucking speeds, the relative motion of the cantilever beam can sufficiently induce the loss of contact, even before the plectrum reaches its geometrically designated separation distance. Moreover, it is interesting to note that the separation time reaches a saturation point as the plucking speed reaches a certain critical value. The saturation value of the separation time may be influenced by the rotational frequency of the plectrum rather than its plucking speed.
As discussed in Figure 3a, increasing the plucking speed leads to a shorter duration of the contact phase, but it also results in a higher maximum contact force. Consequently, the mechanical linear impulse does not exhibit a monotonic relationship with the plucking speed. Figure 3c presents the mechanical linear impulse of the contact force, integrated over the time interval of the contact phase, with respect to the plucking speed. This figure includes an optimal plucking speed that maximizes the linear impulse, thereby yielding the largest free vibration and highest output voltage for the piezoelectric energy harvester. Figure 4 shows the time responses of the free vibration of the piezoelectric energy harvester and the corresponding output voltage at points S1, S3, and S4 (as shown in Figure 3). It is evident from the figure that the displacement obtained at point S3 is greater than those at points S1 and S4. The resulting RMS output voltages are 0.0484, 0.0572, and 0.0493 at points S1, S3, and S4, respectively.  Figure 5a displays the variations in contact force between the piezoelectric cantilever and plectrum during the contact phase. These variations are evaluated for different overlap lengths: 0.5, 1.3, 2.0, and 2.8. The overlap length plays a crucial role in determining the geometrical separation point in the contact phase, but it does not impact the dynamic interaction between the piezoelectric cantilever and plectrum. Therefore, as depicted in Figure 5a, the profiles of the contact force for the different overlap lengths are identical, except for the separation time.
As the overlap length increases, the geometrical separation point also increases proportionally. Consequently, when the overlap length is relatively small (or large), the time duration of the contact phase (or separation time) is also small (or large). The separation time tends to increase as the overlap length reaches the critical point S4, as observed in Figure 5a,b. The red circle point S4 represents the critical condition of the maximum overlap length, where the separation time and linear impulse reach their saturation points and are maximized, as shown in Figure 5b,c. As discussed in Figure 5a, as the overlap length increases, both the contact force and the time duration of the contact phase also increase. Consequently, the mechanical linear impulse tends to increase monotonically with the overlap length, as depicted in Figure 5c. Therefore, the critical point S4 indicates the optimal overlap length for maximizing the linear impulse and output voltage of the piezoelectric energy harvester. The effects of overlap length and plucking speed on the contact force and separation time in the plucking processes of the piezoelectric energy harvester have been recognized as important factors [38]. The analytical results presented in Figures 3 and 5 are in good agreement with the numerical findings reported in the previous work [38], which further supports the validity of the analytical approach proposed in this study.  Figure 6 illustrates the time responses of the piezoelectric energy harvester's free vibration and the corresponding output voltage evaluated at points S1, S2, and S3 (indicated in Figure 5). It is evident that the displacement of the free vibration during the release phase (Figure 6a-c) increases with the overlap length, and this trend is accompanied by a corresponding increase in the associated output voltages (Figure 6d-f) and output powers (Figure 6g-i). The calculated RMS output voltages are 0.0018, 0.0122, and 0.0287 at points S1, S2, and S3, respectively.

Periodic Plucking
In this section, we investigate the dynamic characteristics of the piezoelectric energy harvester when subjected to periodic plucking. The first plucking cycle corresponds to the free vibration of the energy harvester excited by a single plucking, as discussed in the previous section. Subsequent plucking cycles occur with a period determined by the rotational frequency of the plectrum. In the case where the piezoelectric cantilever is initially at rest, the plucking point for the first cycle is located at a trivial position. However, for subsequent cycles, the plucking point may vary if the free vibration of the piezoelectric cantilever is not fully damped out before the next plucking. Hence, the plucking time needs to be determined by solving the equation of zero indentation within the geometrical contact zone (|u| ≤ w p ), as shown below: If a solution to Equation (25) cannot be found within the geometrical contact zone, it indicates that the plectrum does not make contact with the piezoelectric cantilever. If a single solution exists, it represents the plucking time that triggers the subsequent plucking cycle. In this case, the separation time is determined by the geometric separation point and can be evaluated as τ s = w d /Rω presented in Equation (24). Lastly, if two distinct solutions are obtained, the first solution corresponds to the plucking time, while the second solution represents the separation time. This implies that the loss of contact between the piezoelectric cantilever and the plectrum occurs before the plectrum reaches the geometrically designated separation point.  In Figure 7a, the plucking period given as T = 628.3 is relatively long. As a result, the free vibration of the piezoelectric beam is almost damped out before the beam is plucked at rest by the plectrum. Therefore, it appears that the free vibration excited by a single plucking cycle (as discussed in the previous section) is repeated in each plucking cycle. However, from the perspective of energy harvesting efficiency, it would be more desirable to minimize the portion of small-amplitude free vibration within each plucking cycle. This can be achieved by reducing the plucking period.
In Figure 7b, the plucking period is given as T = 125.7. When comparing Figure 7a,b, it can be observed that as the rotational frequency of the plectrum increases, the plucking period becomes shorter, although the dynamic interaction during the contact phase becomes more complex. The first free vibration appears to be relatively larger than the subsequent ones because the relative motion of the cantilever beam in the subsequent plucking cycles potentially reduces the contact force magnitude due to decreased indentation of the contact surface. However, despite this observation, the RMS output voltage in Figure 7b is still higher than that in Figure 7a. Additionally, it is noteworthy in Figure 7b that the plectrum fails to pluck the piezoelectric cantilever in the second cycle. While a shorter plucking period may improve energy harvesting efficiency if the free vibration is successfully triggered in every plucking cycle, it often leads to the failure of the plectrum in periodic plucking.  Figure 8 shows an example of continuous plectrum failure in periodic plucking, obtained with a plucking frequency of ω = 0.17. In this simulation, the plectrum successfully plucks the piezoelectric cantilever at rest only during the first plucking cycle, while subsequently failing to do so for four consecutive cycles. The intentional choice of the rotational frequency in this case causes a lack of synchronization between the motions of the piezoelectric cantilever and the plectrum within their geometrically designated contact zone. In contrast, the plucking frequencies used in Figure 7a,b are specifically chosen to ensure exact synchronization: ω n = mω or T = mτ n , m = 1, 2, 3, · · · (26) where ω n and τ n are the natural angular frequency and natural period of the piezoelectric cantilever oscillator, respectively, which are given by ω n = 1 and τ n = 2π in this study, and m is the positive integer. Equation (26) indicates that the reciprocal of the frequency ratio, ω n /ω, must be a positive integer to achieve an exact synchronization between two periodic motions. This condition is essential as a necessary requirement and a fundamental design strategy, even though nonlinear effects also need to be considered. In Figure 7a,b, the inverse frequency ratios used are 100 and 20, respectively, both of which are positive integers. However, in Figure 8, the inverse value is not an integer but rather 5.88. It should be noted that reducing the plucking speed of the plectrum and increasing the overlap length can offer further advantages in the context of periodic-plucking energy harvesting.

Conclusions
In this study, we proposed an analytical approach based on the modified differential transform method to investigate the dynamic behavior of a mechanical plucking energy harvester. The energy harvester consisted of a piezoelectric cantilever oscillator and a rotating plectrum. Our analytical approach provided a closed-form solution for the beam displacement, allowing us to determine the starting and ending points of the contact phase between the piezoelectric cantilever and the plectrum, as evaluated by Equations (24) and (25). This approach proved valuable in simulating complex dynamic interferences during multiple or periodic plucking processes, where the contact phase varied depending on the relative motion of the beam.
To evaluate the effects of plucking speed and overlap length of the plectrum on single or periodic plucking, we conducted a series of simulations. Our observations showed that the output voltage of the piezoelectric energy harvester increased with the overlap length of the plectrum. This was because a longer contact duration resulted in a stronger linear impulse, triggering larger-amplitude free vibrations of the piezoelectric oscillator. On the other hand, increasing the plucking speed of the plectrum tended to enhance the magnitude of the contact force but shorten the contact phase duration. Therefore, optimizing the plucking speed was crucial to maximize the resulting linear impulse, leading to higher electric energy production from the piezoelectric cantilever.
For successful periodic plucking, it was essential to synchronize the motions of the piezoelectric energy harvester and the rotating plectrum within a limited contact zone. This synchronization could be achieved by matching the rotating frequency of the plectrum with an integer multiple of the natural frequency of the energy harvester. Failure to achieve this synchronization often resulted in the plectrum missing the piezoelectric cantilever in the contact zone due to dynamic interferences. Additionally, reducing the plucking speed of the plectrum and increasing the overlap length were advantageous for successful periodic-plucking energy harvesting.
It is important to note that the Duffing-type oscillator model used in this study represents a simplified approach to simulate the complex dynamics involved in multiple plucking scenarios of a piezoelectric energy harvester. Therefore, detailed design factors such as the geometric configurations of the rotational plectrum and the nonlinear curvature of the beam were not considered in the model. Further research, encompassing both experimental and theoretical studies, is necessary to gain a deeper understanding of the intricate plucking dynamics involved in piezoelectric energy harvesters.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
For a thin beam (based on the Euler-Bernoulli beam theory), the piezoelectric constitutive equations for the stress and electric displacement components are given as where T 1 and D 3 are the stress and electric displacement components, respectively, S 1 and E 3 are the strain and electric field components, respectively, c E 11 , e 31 , and ε S 33 are the elastic, piezoelectric, and permittivity constants, respectively, and the superscript E and S denote that the respective constants are evaluated at the constant electric field and constant strain. In Equation (A1), the strain component can be expressed as where h p and h s are the thicknesses of the piezoelectric layer and substrate, respectively, and κ is the curvature. The electric current output can be obtained from the integral form of Gauss's law: where D is the electric displacement vector, n is the unit outward normal to the electrode area A, and V is the output voltage across the resistance R. Substituting Equations (A1) and (A2) into Equation (A3), the equivalent circuit equation for the piezoelectric layers can be obtained as where w is the tip deflection of the piezoelectric beam, and C p and θ are the equivalent capacitance and electromechanical coupling constant, which are given as where b and L p are the width and length of the piezoelectric layer, respectively, and φ is the mode shape function. Note that the above derivation process was based on the single-mode approximation.

Appendix B
To validate the oscillator model, a comparison is made between the numerical result obtained using Equations (1) and (2) and the experimental one reported in Ref. [37]. To this end, the system parameters of the oscillator model are evaluated using the geometric parameters and material properties of a piezoelectric energy harvester described in Ref. [37], which are listed in Table A1. The transient response of the piezoelectric energy harvester is obtained through direct numerical integration of a linear system of state equations using the Runge-Kutta method implemented with the ode45 function in MATLAB.  Figure A1 shows (upper) the free vibration of the piezoelectric energy harvester excited by periodic plucking and (lower) the corresponding rotational motion of the plectrum. It is worth noting that the plectrum rotates in a counterclockwise direction, while the plucking frequency and plucking period are given as 19.5 Hz and 50.3 ms, respectively. Figure A1 presents four plucking processes with the period of T, followed by the subsequent free vibrations. The numerical result obtained in this study exhibits good agreement with the experimental result reported in Ref. [37], which supports the effectiveness of employing the oscillator model to simulate the intricate dynamics involved in multiple plucking scenarios of a piezoelectric energy harvester.  [37].