*Article* **Modeling, Validation, and Performance of Two Tandem Cylinder Piezoelectric Energy Harvesters in Water Flow**

**Rujun Song 1,2 , Chengwei Hou 1,3 , Chongqiu Yang <sup>1</sup> , Xianhai Yang <sup>1</sup> , Qianjian Guo 1,\* and Xiaobiao Shan 3,\***

<sup>1</sup> School of Mechanical Engineering, Shandong University of Technology, Zibo 255049, China;

	- <sup>2</sup> Shenzhen Research Institute of City University of Hong Kong, Shenzhen 518057, China
	- <sup>3</sup> School of Mechatronics Engineering, Harbin Institute of Technology, Harbin 150001, China

**\*** Correspondence: guoqianjian@163.com (Q.G.); shanxiaobiao@hit.edu.cn (X.S.)

**Abstract:** This paper studies a novel enhanced energy-harvesting method to harvest water flowinduced vibration with a tandem arrangement of two piezoelectric energy harvesters (PEHs) in the direction of flowing water, through simulation modeling and experimental validation. A mathematical model is established by two individual-equivalent single-degree-of-freedom models, coupled with the hydrodynamic force obtained by computational fluid dynamics. Through the simulation analysis, the variation rules of vibration frequency, vibration amplitude, power generation and the distribution of flow field are obtained. And experimental tests are performed to verify the numerical calculation. The experimental and simulation results show that the upstream piezoelectric energy harvester (UPEH) is excited by the vortex-induced vibration, and the maximum value of performance is achieved when the UPEH and the vibration are resonant. As the vortex falls off from the UPEH, the downstream piezoelectric energy harvester (DPEH) generates a responsive beat frequency vibration. Energy-harvesting performance of the DPEH is better than that of the UPEH, especially at high speed flows. The maximum output power of the DPEH (371.7 µW) is 2.56 times of that of the UPEH (145.4 µW), at a specific spacing between the UPEN and the DPEH. Thereupon, the total output power of the two tandem piezoelectric energy harvester systems is significantly greater than that of the common single PEH, which provides a good foreground for further exploration of multiple piezoelectric energy harvesters system.

**Keywords:** piezoelectric energy harvester; tandem; energy harvesting; vortex-induced vibration; flowing water

#### **1. Introduction**

Global environmental and energy policies stress the need to increase the share of renewable resources and to enhance the efficiency of energy conversion plants and commit to developing advanced solutions for power production [1–4]. In addition to the conversion of biogas to electricity, which has attracted many scholars, the use of piezoelectric materials to convert kinetic energy in the environment into electrical energy has also become a hot spot in electrical energy conversion. Fluid kinetic energy is one of the most important clean energy resources and is widely distributed. It is of great value to convert fluid kinetic energy to electricity for human use. There are several efficient methods for energy transformations, such as turbines for hydro-plants, windmills for wind farms, and so on. Both turbines and windmills are suitable for large-scale power generation, which is inconvenient for portable, wearable, low-energy device. The piezoelectric effect can be applied for energy harvesting [5–7] and ultrasonic transduction [8–11]. Especially, piezoelectric energy harvesters (PEHs) have been widely used for energy harvesting from ambient vibration [12–16]. Hence, the piezoelectric conversion mechanism can be used to capture energy from fluid kinetic energy, which has attracted much attention among scholars [17–22].

**Citation:** Song, R.; Hou, C.; Yang, C.; Yang, X.; Guo, Q.; Shan, X. Modeling, Validation, and Performance of Two Tandem Cylinder Piezoelectric Energy Harvesters in Water Flow. *Micromachines* **2021**, *12*, 872. https://doi.org/10.3390/mi12080872

Academic Editors: Qiongfeng Shi and Huicong Liu

Received: 21 June 2021 Accepted: 23 July 2021 Published: 25 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Flow-induced vibration (FIV) and biogas are renewable and alternative sources of energy [23–25]. FIV is a common natural phenomenon, such as vortex-induced vibration of smokestacks and power transmission lines, fluttering flags in airflow, and the flow-induced vibration of seaweed and submarine cables in ocean currents. Biogas and biomass gas (syngas) have been used to feed solid oxide fuel cells [26,27]. Bochentyn et al. [28] proposed using praseodymium and samarium co-doped with ceria as an anode catalyst for DIR-SOFC fueled by biogas, and investigated the catalytic performance. De Lorenzo et al. [29,30] investigated the electrical and thermal analysis of an intermediate-temperature solid oxide fuel cell system fed by biogas. Energy harvesting from the above FIV phenomenon will be an effective method for converting fluid kinetic energy to electric energy, which can be regarded as wake-induced vibration (WIV) [31–34], vortex-induced vibration [35–40], fluttering [41–44] and galloping [45–48], etc. Particularly, PEHs based on the WIV have been extensively investigated in the past decade years. Taylor et al. [49] and Allen et al. [50] first proposed an eel-shaped PEH, based on the WIV to convert water kinetic energy into electric energy. Akaydin et al. [51] studied an airflow PEH by WIV and found that the maximum power output could be obtained when the tip of the piezoelectric beam was located downstream about twice the cylinders' diameter. Weinstein et al. [52] proposed that PEH caused by the awakening of the upstream cylinder, achieving output power of 3 mW at air velocity of 5 m/s. Li et al. [53] conducted a numerical study on the feasibility of the blunt body obtaining energy from steady state flow when *Re* = 100 and found that the performance was closely related to the local vortex dissipation and pressure gradient of large vortex in the wake near the PEH. Yu et al. [54] proposed four structures of cylindrical and PEH membranes. The beat dynamics of PEH membranes were extensively studied, which is helpful to comprehensively understand the energy collection process of elasticity, kinetic energy, and continuous energy transfer. It can be inferred that an enhanced performance of PEH can be achieved by using the WIV phenomenon.

Inspired by the WIV, some researchers have proposed a combined piezoelectric energy harvesting system to improve space utilization and output power performance [55]. Li et al. [56] investigated the wake effect of the two tandem cylinders and found that the aerodynamic forces of the two cylinders were largely determined by the spacing ratio between them. Abdelkefi et al. [57] studied two tandem PEHs in airflow. It was found that when the distance between two PEHs exceeds a critical value, wake galloping of the upstream PEH would occur, and this wake galloping can improve the average output power of the total energy harvesting system. Hu et al. [58] evaluated the performance of twin, adjacent PEHs, based on mutual interference, under wind excitation, and found that the total output power of the two wind energy harvesters could achieve up to 2.2 times that of two single harvesters at the optimal relative position. All above works indicate that the PEHs with tandem arrangement can improve the energy-harvesting capacity of the overall system. Therefore, Shan et al. [59] proposed an energy-harvesting system with two identical piezoelectric energy harvesters in a tandem configuration and investigated the energy-harvesting performance by experimentation. However, the vibration response mechanism of fluid–structure interaction and the relationship between the output power of the harvester and flow field distribution were not studied. Hence, this work proposes a simulation method to study the performance of two tandem cylinder piezoelectric energy harvesters thoroughly and systematically. The validity of the proposed simulation method is verified through the experimental results. Through the simulation, the vibration response and energy-harvesting performance under wake excitation, including the variation rules of vibration frequency, vibration amplitude, power generation, and the distribution of flow field, can be easily obtained.

#### **2. Physical Model and Simulation Method**

Figure 1 shows the physical model of a two tandem PEHs system. Each PEH consisted of a piezoelectric beam and a cylinder. The piezoelectric beam was fixed on the upper end and consisted of a piezoelectric layer (PZT-5H) and a substrate layer. The cylinder was attached to the free end and vertically immersed into incoming water. *U* is the water velocity. The water flowed over the cylinders and stimulated the PEHs to vibrate periodically. The energy harvesting system can be equivalent into two tandem independent single degree of freedom PEHs with separation spacing *W*. As shown in Figure 1, *u* is the vibration displacement of the PEH. *Meq*, *K<sup>b</sup>* are the equivalent mass and equivalent stiffness, respectively. *C* is the damping coefficient of the single PEH. *C<sup>p</sup>* is the equivalent capacitance of the PEH. Θ is the electromechanical coupling coefficient of mechanical vibration and electrical output. *V* is the generated voltage across the electrical load resistance *R*. Θ

**Figure 1.** Schematic diagram of two tandem PEHs system: (**a**) three-dimensional system; (**b**) equivalent model of two tandem PEHs.

According to the equivalent model in Figure 1b, each of the two PEHs can be described as an individual-equivalent single-degree-of-freedom model, written as

$$\begin{cases} M\_{eq}\ddot{u} + \mathbf{C}\dot{u} + \mathbf{K}\_b u + \Theta V = F\_f \\ \frac{V}{R} + \mathbf{C}\_p \dot{V} - \Theta \dot{u} = 0 \end{cases} \tag{1}$$

*ω ξ ξ* where, *F<sup>f</sup>* is the flow-induced force, which should be obtained by the computational fluid dynamics (CFD) method (Fluent® software) due to the complex interactions between the two tandem cylinders in water. *C* is the damping equaling to 2*ωnξ*, and the *ξ* is the damping ratio and can be obtained by free decay oscillation method in experimental text. The second equation is based on the Kirchhoff's law, which express that the sum of all the currents entering a node is equal to the sum of all the currents leaving the node. A key point of the measurement in the experiment is that the cylinder must be immersed into water when the free decay oscillation occurs. As a result, the added fluid damping of cylinder must be taken into account.

Figure 2a illuminates the schematic of the individual PEH. The length of piezoelectric beam and cylinder are *L<sup>b</sup>* , *Lc*. The thickness of piezoelectric and substrate layer are *hp*, *h<sup>s</sup>* , respectively. The width of piezoelectric beam is *b*. The depth of immersion into water flow of the PEH is *L<sup>f</sup>* . Point *O* is the center of flow-induced force action on cylinder. Figure 2b shows the finite element analytical model of the PEH. Elements SOLID5, SOLID45, and CIRCU94, in ANSYS® 10, were used for finite element simulation of the piezoelectric layer, substrate lay, and load resistance circuit, respectively. The sign of the piezoelectric constant was applied for indicating the polarization direction of the PZT-5H. The node on the upper piezoelectric layer surface was named common node "1", and the node on the lower surface was named common node "2". Then the load resistance *R* can be connected between the common nodes "1" and "2".

**Figure 2.** Physical model of the individual PEH: (**a**) Dimension diagram, (**b**) Finite element analytical model.

According to dimension of PEH in Figure 2a, the deflection *w<sup>b</sup>* (*x*, *t*) of the cantilever beam can be expressed as

$$w\_b(\mathbf{x}, t) = \frac{F\_f \mathbf{x}^2}{6EI} (\Re L\_b - \mathbf{x}) + \frac{M\_{Force} \mathbf{x}^2}{2EI} \qquad \mathbf{x} \le L\_b \tag{2}$$

where *EI* is the stiffness of piezoelectric cantilever beam, expressed as

$$EI = b \left[ \frac{E\_p \left( h\_c^{-3} - h\_b^{-3} \right) + E\_s \left( h\_b^{-3} - h\_a^{-3} \right)}{3} \right] \tag{3}$$

where, *ha*, *h<sup>b</sup>* , and *h<sup>c</sup>* are expressed as

$$\begin{aligned} h\_{\mathfrak{a}} &= -\frac{h\_{\mathfrak{k}}}{2} - \frac{\left(h\_{\mathfrak{p}} + h\_{\mathfrak{s}}\right)E\_{\mathfrak{p}}h\_{\mathfrak{p}}}{2\left(E\_{\mathfrak{p}}h\_{\mathfrak{p}} + E\_{\mathfrak{s}}h\_{\mathfrak{s}}\right)}\\ h\_{\mathfrak{b}} &= \frac{h\_{\mathfrak{s}}}{2} - \frac{\left(h\_{\mathfrak{p}} + h\_{\mathfrak{s}}\right)E\_{\mathfrak{p}}h\_{\mathfrak{p}}}{2\left(E\_{\mathfrak{p}}h\_{\mathfrak{p}} + E\_{\mathfrak{s}}h\_{\mathfrak{s}}\right)}\\ h\_{\mathfrak{c}} &= \frac{h\_{\mathfrak{s}}}{2} + h\_{\mathfrak{p}} - \frac{\left(h\_{\mathfrak{p}} + h\_{\mathfrak{s}}\right)E\_{\mathfrak{p}}h\_{\mathfrak{p}}}{2\left(E\_{\mathfrak{p}}h\_{\mathfrak{p}} + E\_{\mathfrak{s}}h\_{\mathfrak{s}}\right)} \end{aligned} \tag{4}$$

2 2 where, *Ep*, *E<sup>s</sup>* are the Young modulus of the piezoelectric and substrate materials.

The moment *MForce* acts on the tail of cantilever beam is expressed as

$$M\_{\text{Force}} = F\_f \left( L\_c - \frac{L\_f}{2} \right) \tag{5}$$
 
$$\text{all cylinder } w(x, t) \text{ is}$$

The displacement of terminal cylinder *w*(*x*, *t*) is

$$w(\mathbf{x},t) = w\_b(L\_{b\prime}t) + w\_b'(L\_{b\prime}t)(\mathbf{x} - L\_b) \qquad L\_b \le \mathbf{x} \le L\_b + L\_c \tag{6}$$

2

where, *w* ′ *b* is the first derivative of *x*. Then equivalent stiffness *K<sup>b</sup>* can be obtained by the following equation, using Hooke's law,

$$\mathbf{K}\_b = \frac{\mathbf{F}\_f}{w(L\_{bo}, t)} \tag{7}$$

− where, *Lbo* = *L<sup>b</sup>* + *L<sup>c</sup>* − *Lf*/2, and the *K<sup>b</sup>* is obtained at last.

The total kinetic energy of the PEH (TPEH) is the sum of kinetic energy of the piezoelectric layer and substrate layer (*Tbeam*), kinetic energy of cylinder (*Tcylinder*), and kinetic energy of its added fluid (*Tfluid*), expressed as follow

$$T\_{PEH} = T\_{beam} + T\_{cylinder} + T\_{fluid} \tag{8}$$

where,

$$\begin{cases} T\_{beam} = \frac{1}{2} \int\_{V\_p} \rho\_p \left[ \frac{\partial w\_b(\mathbf{x}, t)}{\partial t} \right]^2 dV\_p + \frac{1}{2} \int\_{V\_s} \rho\_s \left[ \frac{\partial w\_b(\mathbf{x}, t)}{\partial t} \right]^2 dV\_s\\ T\_{cylinder} = \frac{1}{2} \int\_{L\_b}^{L\_b + L\_c} \frac{M\_c}{L\_c} \left[ \frac{\partial w(\mathbf{x}, t)}{\partial t} \right]^2 dx\\ T\_{fluid} = \frac{1}{2} \int\_{L\_b + L\_c - L\_f}^{L\_b + L\_c} \frac{M\_{cf}}{L\_f} \left[ \frac{\partial w(\mathbf{x}, t)}{\partial t} \right]^2 dx \end{cases} \tag{9}$$

where, *ρ<sup>s</sup>* is the density of substrate material, *ρ<sup>p</sup>* is the density of piezoelectric material, *M<sup>c</sup>* is the mass of cylinder, and *Mcf* is the fluid added mass, written as

$$M\_{cf} = C\_M \frac{D^2}{4} \pi L\_f \rho\_f \tag{10}$$

where *C<sup>M</sup>* is the additional mass coefficient, *C<sup>M</sup>* = 1.

Meanwhile, the kinetic energy TPEH can be re-expressed as

$$T\_{PEH} = \frac{1}{2} \mathcal{M}\_{eq} \left[ \dot{w} (L\_{bo}, t) \right]^2 \tag{11}$$

and the *Meq* can be obtained based on Equations (8) and (11). Thus the equivalent mass *Meq* can be obtained by the equation:

$$\mathcal{M}\_{\mathsf{dq}} = \begin{bmatrix} \int\_{V\_{p}} \rho\_{p} \left[ \frac{\partial w\_{b}(\mathbf{x},t)}{\partial t} \right]^{2} dV\_{p} + \int\_{V\_{b}} \rho\_{s} \left[ \frac{\partial w\_{b}(\mathbf{x},t)}{\partial t} \right]^{2} dV\_{\mathbf{s}} \\ + \int\_{L\_{b}}^{L\_{b} + L\_{c}} \frac{M\_{c}}{L\_{c}} \left[ \frac{\partial w(\mathbf{x},t)}{\partial t} \right]^{2} d\mathbf{x} + \int\_{L\_{b} + L\_{c} - L\_{f}}^{L\_{b} + L\_{c}} \frac{M\_{cf}}{L\_{f}} \left[ \frac{\partial w(\mathbf{x},t)}{\partial t} \right]^{2} d\mathbf{x} \end{bmatrix} / \left[ \dot{w}(L\_{bw}, t) \right]^{2}$$

The Θ of electromechanical coupling coefficient of the single PEH in energy harvesting system can be expressed as

$$
\Theta = \sqrt{(w\_{\rm noc}^2 - w\_{\rm nsc}^2)M\_{eq}\mathcal{C}\_p} \tag{12}
$$

where, *ωnoc* is open circuit resonance frequency, and *ωnsc* is short circuit resonance frequency, which will be obtained by *R*<sup>0</sup> = 0 and *R*0→∞ in ANSYS® software, respectively. Equivalent capacitance *C<sup>p</sup>* can be

$$\mathcal{C}\_p = \frac{\varepsilon\_{33}^S b L\_b}{h\_p} \tag{13}$$

where *ε S* <sup>33</sup> is the dielectric constant of piezoelectric material.

Furthermore, the output power of the PEH can be expressed as

$$P = \frac{1}{T} \int\_0^T \frac{V^2(t)}{R} dt\tag{14}$$

where *T* is the alternating period of output voltage.

The simulation results can be obtained by Equations (1) and (14) using the User Defined Function (UDF) in Fluent® software. The simulation process can be summarized as:


Figure 3 illuminates the schematic of the flow field, including the whole grid arrangement and partial grid diagram of two tandem PEHs. The dynamic grid was applied in this CFD process. The domain of flow field grid was 20*D* × (40*D* + *W*), where *D* is the cylinder diameter and *W* is the spacing of the two tandem PEHs. The spacing of inlet to upstream piezoelectric energy harvester (UPEH) was 10 times that of *D*. The spacing of outlet to downstream piezoelectric energy harvester (DPEH) was 30 times that of *D*. In order to ensure the grid quality near the cylinders during the simulation process, the computational domain was divided into three parts, as shown in Figure 3a. The first part was the static grid. The second part was the dynamic grid, where the third part could move with the two tandem PEHs during the CFD process to reduce grid distortion, as shown in Figure 3b.

**Figure 3.** Schematic of the flow field grid of the tandem arrangement of energy harvesting system: (**a**) grid diagram of flow field; (**b**) partial grid diagram of third part.

#### **3. Experimental Setup and Simulation Verification**

Figure 4 shows the experimental platform and the prototype setup. The experiment platform was made up of an open channel, fixtures, and measurement system. In order to obtain a steady flow, some sections in the open channel, such as the contraction section, damping screens, and a cellular device were set up to reduce stream turbulence. The speed of the water was controlled by a pump driven by a frequency converter. The fixtures were used to fix the two tandem PEHs, the spacing *W* was changed by relative position of the fixtures. The measurement system was composed of an adjustable resistor, a Data Acquisition system (DAQ) from National Instruments (NI), and a computer. The output voltage of the resistor was monitored by DAQ and displayed on the computer.

**Figure 4.** Energy harvesting experimental platform and prototypes of two tandem PEHs.

Figure 5 shows the output voltage evolution curve of the upstream and downstream PEHs with time when the spacing ratio *W*/*D* is 3.3 at the flow velocity of 0.23 m/s, 0.31 m/s and 0.37 m/s, respectively. The dimension and parameters of the two tandem PEHs are listed in Table 1. The comparison shows that when the flow velocity is 0.23 m/s, the output voltage of the upstream PEH was greater than that of the downstream PEH. At a flow velocity of 0.37 m/s, the output voltage of the upstream PEH was smaller than that of the downstream PEH. By comparing Figure 5a,c,e, it can be found that the output voltage of the upstream PEH at the velocity of 0.31 m/s is higher than that at the flow velocity of 0.23 m/s and 0.37 m/s, respectively. By comparing Figure 5b,d,f, it can be seen that the output voltage of the downstream PEH increases with the increase in velocity. It can be seen that under the tandem arrangement of the two PEHs, due to the strong coupling effect, the performance of the two PEHs will change greatly, especially the downstream one.

**Figure 5.** Output voltage of two tandem PEHs versus time for *W*/*D* = 3.3. (**a**) Output voltage of UPEH at *U* = 0.23 m/s; (**b**) Output voltage of DPEH at *U* = 0.23 m/s; (**c**) Output voltage of UPEH at *U* = 0.31 m/s; (**d**) Output voltage of DPEH at *U* = 0.31 m/s; (**e**) Output voltage of UPEH at *U* = 0.37 m/s; (**f**) Output voltage of DPEH at *U* = 0.37 m/s.


**Table 1.** Dimension and material properties of the PEHs.

In order to research the energy-harvesting performance and verify the effectiveness of simulation method, Figure 6 illustrates the output power comparison between simulation and experimental results of the two tandem PEHs when the spacing ratio *W*/*D* is 3.3, 4.58, and 6.25, respectively. As can be seen from Figure 6a–c, when the spacing was smaller, the coupling effect was stronger, and the output power of the UPEH was slightly improved. With the increase in spacing ratio, the coupling effect was weakened and the energy harvesting performance of the UPEH decreased and gradually tended to resemble the performance of a single PEH. The experimental results showed that peak output power (145.4 µW of UPEH and 371.7 µW of DPEH) was generated at the velocities of 0.31 m/s and 0.41 m/s when *W/D* = 3.3, respectively. The maximum output power of DPEH was 2.56 times of that of the UPEH. When *W/D* = 6.25, the maximum output power of the UPEH was 117.9 µW and the maximum output power of the DPEH was 239.9 µW. It can be seen that the output performance of DPEH was better than that of UPEH, mainly occurred at the high flow velocity (>0.31 m/s). Especially, the output power of DPEH (371.7 µW) was 12.8 times that of UPEH (29.2 µW) when velocity was 0.41 m/s and *W/D* = 3.3. It can be summarized that the coupling effect of each PEH decreased with the increase in the spacing ratio of *W*/*D*. Meanwhile, the wake effect of the DPEH from the UPEH was gradually decreased, and the disturbance effect on UPEH from the DPEH was also weakened. It's worth noting that the comparison between the theoretical results and the experimental data is better for high *W/D* values. This is because, in the experiment, the smaller the distance, the more serious the impact of the coupled vibration. The mutual disturbance is bound to cause disturbances in the flow field, aggravate wall disturbances, and change the depth of the cylinder's water. In the theoretical analysis, the wall of the flow field is far enough away and the immersion depth of the two cylinders is assumed to be constant. This leads to smaller distances producing greater error between the simulation analysis results and the experimental results. The experimental results also show that the output power of UPEH and DPEH increased to the maximum first, and then decreased with the increase of the flow velocity, and the optimal velocity of the maximum power output could be obtained. In addition, the energy harvesting performance of a single piezoelectric energy harvester (SPEH) is given in Figure 6d. It can be noted that the simulation flow velocity corresponding to the maximum output power is lower than the experimental results. This is mainly because when the vortex-induced resonance occurs, the column immersion depth is slightly reduced, and the additional mass of water to the cylinder is slightly reduced too, but the immersion depth is assumed to be unchanged during the simulation, so that the natural frequency of the energy vibrator during simulation is slightly lower, which causes the flow velocity corresponding to the vortex-induced resonance to be slightly smaller than the experimental result. By comparing Figure 5a–c,

it can be found that the power generation capacity of the UPEH and DPEH were greater than that of the SPEH. Therefore, it can be seen that the two tandem PEHs can improve the energy-harvesting performance of the system. What's more, the overall results of the experiment and simulation are mild, and the error is within the allowable range. The main sources of error can be summarized as: in the experiment, the left and right swing caused a certain change in the depth of the cylinder's water inlet. In the simulation, it is treated as quantitative, and some assumptions in the simplification process of the two-dimensional model will also have certain contrast errors.

**Figure 6.** Output power of simulation results (SRs) and experimental results (ERs) of two tandem PEHs (**a**–**c**) and single PEH (**d**).

Therefore, the simulation results were consistent with the experimental results, as shown in Figure 6; the simulation method was verified.

#### **4. Simulation Results and Discussion**

In order to further investigate the performance of the two tandem PEHs, Figure 7a,b shows changes of vibration frequencies versus velocity. It can be concluded that frequency increases with velocity as a whole. However, as for the DPEH, the vibration performance is complicated. As shown in Figure 7a, the frequency increases with the increase of water velocity when water velocity (*U*) is less than 0.211 m/s, and then decreases with the increase of water velocity when 0.211 m/s < *U* < 0.23 m/s, and lastly increases with the increase of water velocity when *U* > 0.23 m/s. When the *W/D* = 12.5, the vibration frequency of DPEH was smaller than that of UPEH at low-speed flow, as shown in Figure 7b. Figure 7c–f are plotted for when *U* = 0.272 m/s. Due to the coupling effect of the vortex-induced vibration (VIV) itself and the wake-induced vibration (WIV) from the UPEH, there was inevitable unbalanced vibration in the DPEH, with small vibration frequency differences, such as the beat vibration signal, containing two frequency components of dominant frequency and non-dominant frequency when *W/D* = 6.25, which can be found that Figure 7c,d. When *W/D* = 12.5, the coupling effect becomes light, so that the frequency is unique as shown in Figure 7e,f.

**Figure 7.** Vibration frequencies of the PEHs versus velocity when *W*/*D* = 6.25 and *W*/*D* = 12.5. (**a**) Vibration frequencies when *W*/*D* = 6.25. (**b**) Vibration frequencies when *W*/*D* = 12.25. (**c**) Vibration response when *W*/*D* = 6.25 and *U* = 0.272 m/s. (**d**) Fast Fourier transform (FFT) when *W*/*D* = 6.25 and *U* = 0.272 m/s. (**e**) Vibration response when *W*/*D* = 12.25 and *U* = 0.272 m/s. (**f**) FFT when *W*/*D* = 12.25 and *U* = 0.272 m/s.

In order to further analyze the vibration performance of the two tandem PEHs system, the vibration amplitudes of the upstream and downstream PEH with spacing ratio *W/D* and water speed were numerically discussed, specifically. As seen in Figure 8, Part I is the area of water velocity and Part II is the area of high-speed water velocity. As for the UPEH, the vibration amplitude was not mainly affected by the spacing ratio *W*/*D*, but rather by the water velocity. The vibration amplitude of the UPEH first increased and then decreased with the increase in velocity, and reached its maximum value when the velocity was about 0.27 m/s. As for the DPEH, the vibration amplitude first decreased to the minimum value when the velocity was small, as shown in Part I in Figure 8b. This phenomenon was attributed to the flow field transition of the shear layer around the DPEH; the vibration response was weakened during the transformation of the water flow from laminar to turbulent. With further increase in the water velocity, the amplitude of vibration on the DPEH increased to its maximum point, and then gradually decreased to a stable value, as shown in Figure 8b. Due to the excitation effect of upstream vortex shedding

from the UPEH, the vibration response of DPEH was obviously greater than that of UPEH, especially at the high water speed, as shown in Part II in Figure 8a,b respectively.

**Figure 8.** Comparison of vibration amplitudes of the upstream and downstream PEH: (**a**) UPEH, (**b**) DPEH.

In order to study the energy-harvesting performance of the system, the relationship among the output power of tandem PEHs, spacing ratio *W*/*D,* and water velocity are illustrated in Figure 9. It can be found that the output power of UPEH increased first and then decreased with the increase in water speed. The maximum output power was obtained at the speed of 0.27 m/s, at which point vortex-induced resonance occurred, combining with the vibration responses in Figure 8a. As for the DPEH in Figure 9b, the relationship between *W/D* and the energy-harvesting performance of the DPEH was the same as that of the UPEH. Comparing with the UPEH, the output power of the DPEH decreased slowly with the increase in velocity, due to wake stimulation from upstream, which was different from the performance of the UPEH. In addition, it can be concluded that in the range of low water velocity, the output power of the UPEH was higher than that of the DPEH, but in the range of high water velocity, the energy-harvesting performance of the DPEH was obviously better than that of the UPE,H due to wake stimulation by vortex shedding, which is in accordance with the vibration response results in Figure 8.

**Figure 9.** Energy harvesting performance of the two tandem PEHs versus spacing ratio *W*/*D* and water velocity: (**a**) UPEH, (**b**) DPEH.

Further analyzing the coupling effect of the two tandem piezoelectric energy harvester, Figure 10 shows the velocity flow field contour in *x* direction with different spacing ratios and water velocities. It can be seen from Figure 10a,c that at a low flow speed (*U* = 0.19 m/s), a long low-speed flow region could be formed in the downstream after the water flowed through the upstream PEH (cylinder). At this time, the downstream PEH would be in the low-speed flow region. Therefore, the vibration response and power generation performance of the downstream PEH were weak at the low-speed water flow,

as shown in Figure 9b. At the same time, the vibration response of the upstream PEH was similar to that of a single PEH at low speed flow, because the downstream PEH had a weak vibration response and its influence on the upstream PEH was negligible. According to Figure 10b,d, when the flow velocity was high (*U* = 0.41 m/s), the downstream PEH was located on the vortex path where the upstream PEH fallen off, and the surrounding vortices were relatively strong. Under stimulation of the vortex, the vibration response and power output performance of the downstream PEH would be greatly improved. In addition, as the vortex moved forward with the flow field, the vortex diffused and its turbulence strength decreased, so the vibration response and energy harvesting performance of the downstream PEH gradually decreased with the increase in spacing ratio.

**Figure 10.** Flow field distribution of two tandem PEHs for various *W*/*D* and water velocities.(**a**) Flow field distribution at *W*/*D* = 3.33, *U* = 0.19 m/s; (**b**) Flow field distribution at *W*/*D* = 3.33, *U* = 0.41 m/s; (**c**) Flow field distribution at *W*/*D* = 3.33, *U* = 0.19 m/s; (**d**) Flow field distribution at *W*/*D* = 8.33, *U* = 0.41 m/s.

#### **5. Conclusions**

In this paper, a tandem arrangement of piezoelectric energy harvesters was studied to scavenge the water flow vibration energy through simulation modeling and experimental validation. Through the simulation analysis, the variation rules of vibration frequency, vibration amplitude, power generation, and the distribution of flow field were obtained. The experiment results verify the accuracy of the simulation results. The effects of spacing ratio and water flow velocity on the vibration response and the output power of two tandem PEHs were studied numerically and experimentally. It could be concluded that the vibration response and energy-harvesting performance of the two tandem PEHs were enhanced due to the coupling effect induced by upstream vortex vibration. The vibration response of UPEH was vortex-induced vibration, and maximum power could be achieved at the resonant point of the PEHs and the vortex flow. When the spacing was small, the performance was enhanced due to enhancement of the DPEH. The vibration response of the DPEH was raised under the stimulation of the vortices shedding from the UPEH, especially under high-speed water flow; the energy-harvesting performance of the DPEH as better than that of the UPEH accordingly. However, when the flow speed was low, it was noticeable that the superiority of the DPEH energy harvesting was not remarkable and the power output of the UPEH was better than that of the DPEH. The results show

that the energy-harvesting performance of the tandem PEHs in flowing water could be significantly improved, especially at the high speed flow region, which provides good support for further exploration of energy harvesting systems with multi-piezoelectric harvesters, parallel or serially arranged.

**Author Contributions:** R.S. contributed to the whole work including design, fabrication, testing of the device, and writing—original draft preparation; C.H. contributed to the simulation analysis and writing—review and editing; C.Y. contributed to fabrication and simulation analysis; X.Y. contributed to the testing guidance; Q.G. contributed to the writing—review and editing and X.S. contributed to the research idea and gave guidance to the work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant numbers 51705296, 51875116 and 52075306), Key Research and Development Project of Shandong Province under Grant 2019GGX104033, Shenzhen Science and Technology Innovation Committee under Grant No. JCYJ20200109143206663, Key Research and Development Project of Zibo City under Grant No. 2020SNPT0088 and Young Innovative Talents Introduction and Training Program Project of Shandong Provincial Department of Education.

**Acknowledgments:** The authors are also grateful to the experimental support of Harbin Institute of Technology.

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

#### **References**


**Yang Huang 1,2 , Zhiran Yi 1,3, Guosheng Hu 1,2 and Bin Yang 1,2,\***


**Abstract:** A data-driven optimization strategy based on a generalized pattern search (GPS) algorithm is proposed to automatically optimize piezoelectric energy harvesters (PEHs). As a direct search method, GPS can iteratively solve the derivative-free optimization problem. Taking the finite element method (FEM) as the solver and the GPS algorithm as the optimizer, the automatic interaction between the solver and optimizer ensures optimization with minimum human efforts, saving designers' time and performing a more precise exploration in the parameter space to obtain better results. When employing it for the optimization of PEHs, the optimal length and thickness of PZT were 6.0 mm and 4.6 µm, respectively. Compared with reported high-output PEHs, this optimal structure showed an increase of 371% in output power, an improvement by 1000% in normalized power density, and a reduction of 254% in resonant frequency. Furthermore, Spearman's rank correlation coefficient was calculated for evaluating the correlation among geometric parameters and output performance such as resonant frequency and output power, which provides a data-based perspective on the design and optimization of PEHs.

**Keywords:** piezoelectric; energy harvester; optimization; pattern search; FEM; PZT

#### **1. Introduction**

With the urgent demand of sustainable power supplies for low-power electronic applications such as wireless sensor network systems [1–4] in Internet of Things, implantable medical devices [5,6] and other devices in some extreme environments [7], energy harvesting from the ambient environments has attracted broad attention and provided potential solutions to the periodical replacement of batteries during the last few decades. Vibration energy is ubiquitous and robust mechanical energy that exists widely, including bridges [8], roads [9], human body [10–12], cars [7], etc. Therefore, a variety of vibrationbased piezoelectric energy harvesters (PEHs) with different structures have been proposed and studied [13–15]. Among these, the cantilever structure is widely used due to structure simplicity and the high average strain obtained by a given input force [2,16–19].

For cantilever PEHs, a multi-parameter coupling problem exists for obtaining highefficiency energy conversion. To improve the output performance of cantilever PEHs, researchers have studied the effects of different geometric parameters on output performance [20,21]. He et al. [22] and Jia et al. [23] proposed that the optimal mass-beam length ratio is 0.6~0.7 within linear response. Hu et al. [18] investigated the optimal length of the piezoelectric layers based on theoretical analysis, FEM simulation and experimental verification, from which they discovered that the optimal length ratio of piezoelectric layers and the beam is approximately 0.2. Furthermore, Hu et al. [18] also reported that

**Citation:** Huang, Y.; Yi, Z.; Hu, G.; Yang, B. Data-Driven Optimization of Piezoelectric Energy Harvesters via Pattern Search Algorithm. *Micromachines* **2021**, *12*, 561. https://doi.org/10.3390/mi12050561

Academic Editor: Ju-Hyuck Lee

Received: 14 April 2021 Accepted: 13 May 2021 Published: 15 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the optimal PZT layout length decreases as cantilever width increases. However, investigations on various geometric parameters in the studies above constitute single-variable optimization. Moreover, the geometric parameters were manually set and tested at a fixed interval based on the researchers' experience and intuition, which required the researcher to spend more time in trying each possible combination of geometric parameters to obtain the optimal output performance. To overcome these shortcomings, several data-driven optimization strategies combining different algorithms to maximize the output performance of PEHs have been proposed recently. For example, data-driven optimization can solve the optimization problems that are difficult to formulate or solve, in that they are non-convex, multi-objective or multi-modal; it can solve optimization problems based on derivative-free data. Ghoddus et al. [24] presented an optimization approach based on the Particle swarm optimization (PSO) algorithm to maximize the output power of PEHs with four different structures, which simultaneously optimized multiple geometric parameters. In addition, Nabavi and Zhang [25,26] proposed an analytical model for PEHs with proof mass and used a genetic algorithm to simultaneously optimize multiple objectives, including resonant frequency, output power and device volume. However, finding the optimal solution to complex high-dimensional problems such as the multi-parameter coupling problem for PEHs may require expensive objective function evaluations. Therefore, as population-based stochastic optimization techniques, the PSO and GA algorithm implemented in the previous works [24–26] may encounter low efficiency when dealing with the multi-parameter coupling problems. Moreover, parameter tuning is needed to obtain a better convergence for the PSO algorithm and GA [27], such as particle number, accelerate constant, inertia weight and population size. Although the global search ability of GA and PSO may be better, the generalized pattern search (GPS) algorithm is more suitable for the time-consuming multi-parameter coupling problem due to its effortless parameter tuning and relatively fewer objective function evaluations [28].

In this paper, we proposed a data-driven optimization strategy based on the GPS algorithm. The GPS algorithm is one of the direct search methods, which can be effectively used in derivative-free optimization. By implementing the proposed scheme, multiple geometric parameters (*lp*, *w*, *tp*, and *t<sup>b</sup>* ) were simultaneously optimized for maximum output performance based on the data without building an analytical model. Using the proposed optimization method can not only efficiently optimize the geometric parameters of PEHs to achieve high output performance but also save time for analytical modeling and complex parameter tuning. Furthermore, Spearman's rank correlation coefficient was calculated for evaluating the correlation among geometric parameters and output performance, providing a data-based perspective on the design of PEHs.

#### **2. Methods**

The unimorph PEH is composed of a piezoelectric layer, a structural layer, and a proof mass. The piezoelectric layer partially covers the flexible structural layer along the constraint end, as shown in Figure 1a. A tungsten proof mass is fixed at the free end of the structural layer, increasing average strain under given excitation and lowering the resonant frequency to meet the low-frequency ambient vibration. Figure 1b illustrates the annotation of this unimorph PEH. The data-driven optimization scheme mainly includes a FEM solver and a GPS optimizer. Figure 1c depicts the overall function of the presented strategy. In this paper, FEM simulation was carried out using COMSOL Multiphysics (Version 5.4) and the material properties used in FEM simulation are listed in Table 1. The GPS algorithm was implemented in MATLAB (Version 2020a). The communications between solver and optimizer were implemented using COMSOL LiveLinkTM for the MATLAB module, which facilitated the data-driven optimization strategy. As direct search methods, different types of pattern search algorithms have been utilized for solving engineering optimization problems [28–30]. The convergence analysis of GPS algorithm has been performed by Audet and Dennis [31].

**Figure 1.** Schematic diagram of the piezoelectric energy harvester (PEH) and the proposed strategy. (**a**) Three-dimensional view. (**b**) Annotation for geometric parameters of the cantilever PEH. (**c**) Overview of the presented optimization working process.


 () . ≜ { ∈ ℝ ∣ = ′ + Δ ⋅ νk, ∈ {1,2, ⋯ , } } At each iteration, the GPS searches a set of points around the current point called mesh, finding a better point whose value of the objective function is lower than the value before. Then, the better point is set as the current point at the next iteration. Based on this procedure (Polling), the GPS finds a sequence of points that approaches the optimum, and it does not stop until the convergence is met. Specifically, let *M<sup>k</sup>* denote the mesh at *k*-th iteration, and *x* (*i*) *k* denote the *i*-th point in the *M<sup>k</sup>* . The mesh is defined as Equation (1):

$$M\_k \triangleq \left\{ \mathbf{x} \in \mathbb{R}^n \mid \mathbf{x} = \mathbf{x}'\_k + \Delta\_k \cdot \mathbf{v}\_{\mathbf{k}'} \ k \in \{1, 2, \dots, n\} \right\} \tag{1}$$

() ( () ) () ( () ) < ( ′ ) +1 ′ = () Δk+1 = 2 ⋅ Δ ℎ +1 ′ = ′ Δk+1 = 1/2 ⋅ Δ where *x* ′ *k* denotes the current point at the *k*-th iteration while ∆*<sup>k</sup>* denotes mesh size, and ν<sup>k</sup> denotes pattern. Pattern is a set of vectors used to determine the generation of mesh. Suppose *F*(*x*) denotes the objective function in the optimization problem. At each iteration, GPS computes *F x* (*i*) *k* and looks for a better point *x* (*j*) *k* so that *F x* (*j*) *k* < *F*(*x<sup>k</sup>* ′ ); then, *x* ′ *<sup>k</sup>*+<sup>1</sup> = *x* (*j*) *k* and ∆*k*+<sup>1</sup> = 2 · ∆*<sup>k</sup>* are set. Otherwise, if the polling fails to find a better point

νk

at *k th* iteration, then *x* ′ *<sup>k</sup>*+<sup>1</sup> = *x* ′ *k* and ∆*k*+<sup>1</sup> = 1/2 · ∆*<sup>k</sup>* . The computation in the mesh at each iteration is called polling. In addition, the search method runs before polling can select a different current point, which may accelerate the optimization if tuned well. Various search methods can be set, including the genetic algorithm, Latin hypercube search, etc.

The proposed optimization strategy uses FEM simulation as the solver and GPS algorithm as the optimizer, whose workflow is depicted in Figure 2. The initialization step is mainly used for setting parameters in the GPS algorithm, including initial point, the searching and polling method, stopping criteria and other parameters used to accelerate the optimization. Moreover, a new set of geometric parameters generated by GPS is used for FEM modelling, and then, the FEM simulation is performed in order to extract the solutions for feedback to the GPS algorithm. In order to efficiently harvest energy in daily life such as cars, bridges and the human body, the development of PEHs tends to have high output power density and low resonant frequency. Therefore, the optimization objective function here is defined as normalized power density or output power. The definition of normalized power density is shown as Equation (2), which is a function of various geometric parameters and external excitation:

$$P\_n = \frac{P}{a^2 \cdot f\_1 \cdot V\_{eff}} \tag{2}$$

where *P<sup>n</sup>* denotes the normalized power density, *P* denotes the output power, *f*<sup>1</sup> denotes the first-modal resonant frequency, *Ve f f* denotes the effective volume of the specific structure, and *a* represent the excitation acceleration. The maximizing of normalized power density may indicate the trend of maximizing output power and minimizing the first-modal resonant frequency and effective volume. Therefore, the optimization problem is defined as below: <sup>1</sup> 

**Figure 2.** Workflow of the proposed data-driven optimization strategy. (**a**) Overall working mechanism. (**b**) Detailed workflow of generalized pattern search.

#### **3. Results and Discussion**

Utilizing the well-described micromachining processes and experimental setup in our previous work [17], we first fabricated and tested four devices with different length PZT layers to validate the effectiveness of the proposed data-driven optimization strategy, and then, individually optimized the PZT length and the proof mass length intending to maximize output power and compared the optimized result with the previous works [18,23]. The experimental setup is shown in Figure 3a. Controlled by the vibration controller through the power amplifier (YE2706A, Sinocera Piezotronics, Inc., Yangzhou, China), the force applied on different PEHs is generated by a shaker (JK-2, Sinocera Piezotronics, Inc.) and is monitored by a force sensor (208C02, PCB Piezotronics, New York, NY, USA). Firstly, the bulk PZT (300 μm) and beryllium bronze (50 μm

**Figure 3.** (**a**) Photograph of the experimental setup. (**b**) Photograph of the fabricated harvesters with the varied PZT layer length.

Next, the bulk PZT was thinned to around 50 μm through chemical mechanical thin-Firstly, the bulk PZT (300 µm) and beryllium bronze (50 µm) were polished so as to increase the bonding strength. Then, the bulk PZT was bonded on the beryllium bronze using conductive silver epoxy in a vacuum oven at 175 ◦C for 3.5 h.

varied from 1 kΩ to 1000 kΩ to determine the optimum load resistance and output power Next, the bulk PZT was thinned to around 50 µm through chemical mechanical thinning and polishing. After that, 20/200 nm Cr/Au was sputtered on the polished surface of the bulk PZT thick film as an electrode. Finally, the cantilevered PEHs were patterned using the ultraviolet laser method, and the proof mass was assembled at the free end. The fabricated harvesters with different PZT layer length are illustrated in Figure 3b. Under 1.0 g acceleration, the PEHs were connected in serial with the external resistance, which varied from 1 kΩ to 1000 kΩ to determine the optimum load resistance and output power. For performing FEM simulation, a geometrical configuration was considered where the free end was mounted with a proof mass, and the upper surface of the PZT was grounded. In addition, chamfer was added between the proof mass and structural layer to avoid stress singularities at the reentrant corners. Meshes were created according to the shape of the geometry. Then, the detailed meshes were determined by carrying out a mesh convergence study. Figure 4a presents the meshes of the PEHs constructed in COMSOL software. Here, we used skewness, the default quality measure in COMSOL software, to evaluate the mesh quality, which is a suitable measure for most types of meshes. To extract the optimum output power of each PEH, the FEM simulation firstly ran an eigenfrequency study to obtain the first-modal eigenfrequency of the PEHs, and then ran a frequency domain study under the first-modal eigenfrequency with an auxiliary sweep of different external resistance. The excitation was set as body load with an acceleration of 1.0 g (9.8 m/s 2 ). Using the proposed scheme, the procedures mentioned above will not stop until the GPS algorithm meets convergence, and the structure of PEHs will keep updated to seek the optimum geometric parameters.

> Figure 4b shows a comparison of normalized output power for the presented PEHs with the varied PZT layer length under different external resistance between FEM simulation and experimental results. Compared with the experimental data, the FEM simulation

showed agreement in the trend of maximum output power with varied PZT layer length, which proved the validation of the FEM simulation in determining the optimal geometric parameters. Especially when the length of the PZT layer is 3.0 mm, the FEM simulation fits well with the experimental data. The difference in output power corresponding to external resistance between FEM simulation and experimental results may come from manufacturing and testing errors. As the length of the PZT layer increases, the manufacturing errors accumulate. The difference between FEM and the experiment decreases as the length of PZT decreases. Moreover, the mechanical damping (set as 0.015) was assumed to be a constant in the FEM simulation, which may also lead to errors in power prediction for different structure [32]. Then, two single-variable optimizations (PZT length and mass length) for maximum output performance were performed and compared with the previous works [18,23]. For the cantilever harvester with proof mass, Hu et al. [18] and Jia et al. [23] have, respectively, optimized the PZT length and the mass length to maximize output power, from which they have concluded that the optimal ratio of PZT length and total length is approximately 0.2 (3 mm/15 mm), while the optimal mass-beam length ratio is 0.6~0.7.

**Figure 4.** (**a**) The mesh of the PEHs in FEM simulation. Validation of the proposed scheme by (**b**) comparison of normalized output power for the presented PEHs with the varied PZT layer length under different external resistance between FEM simulation and experimental results. (**c**) Convergence curve of output power for each single-variable optimization given by the proposed scheme. (**d**) All results searched by the proposed scheme for each single-variable optimization. The applied acceleration amplitude is 1.0 g.

Figure 4c shows the convergence curve of mass and PZT length optimization. It can be observed that the GPS algorithm can efficiently optimize the proposed optimization problems within two and four iterations, respectively, for PZT length and mass length optimization. Figure 4d presents all results evaluated by the GPS algorithm, which demonstrates the trend of output power versus varied length of the mass and PZT. The optimal PZT length and mass length are 3.437 mm and 10.090 mm corresponding to the ratio of the

total length of 0.229 and 0.672, respectively. The agreement of the optimized results and that of the previous works [18,23] proves the efficiency and effectiveness of the proposed data-driven optimization method and the FEM model. Furthermore, unlike the optimization approaches in the previous works [18,23], by which researchers manually varied the geometric parameters at fixed intervals to explore the parameter space for maximum output performance, the proposed method automatically performs the optimization task once set, saving researchers time compared to trial-and-error approaches. In addition, a more precise exploration in the parameter space can be carried out by the proposed scheme compared with the manual trial-and-error approaches. Respectively labeled as 1OPT-1 and 1OPT-2, the geometric parameters and output performance of the above two single-variable optimizations are summarized in Table 2.

**Table 2.** Geometric parameters and results of our calculations.


*l<sup>b</sup>* = 15.00 mm; *t<sup>m</sup>* = 0.20 mm; Acceleration = 1 g.

After verifying the effectiveness of the proposed scheme, further optimization for the unimorph PEHs was implemented to obtain maximum output power and minimum resonant frequency. The optimization problem is depicted as Equation (3). The geometric parameters to be optimized were set as *lp*, *w*, *t<sup>p</sup>* and *t<sup>b</sup>* . As the thickness of the whole beam becomes thinner, the higher the probability of device failure. Based on the experimental data in the previous work [18], the upper and lower bound of *t<sup>b</sup>* was set as 30–50 mm from a conservative consideration. The optimal results and the initial point in the optimization, labeled as 4OPT-1 and Ref, are summarized in Table 2.

Figure 5a presents the trajectory of the geometric parameters and output power during optimization, which directly shows how the GPS algorithm optimizes the defined optimizable parameters to obtain better values of objectives function. Each point, respectively, represents the best point at every iteration. As shown in Figure 5a and Table 2, the optimal length and thickness of PZT are 6.0 mm and 4.6 µm, respectively. In addition, the reduction in thickness of the structural layer also contributes to the improvement of output power. Figure 5b demonstrates the comparison of the output performance between the reference structure (Ref) and optimized structure (4OPT-1). It can be observed that the resonant frequency of 4OPT-1 is 2.54 times lower than that of Ref, while the output power and normalized power density were, respectively, about 3.71 and 10.10 times larger, showing a substantial improvement in output performance. Numerically, the reduction in resonant frequency greatly contributed to the improvement of normalized power density according to the definition given by Equation (2). By calculating the second derivative of the tip displacement at the first-modal vibration, the strain of the piezoelectric layer along the arc length (S1(x, y, t) = −y *∂* <sup>2</sup>z(x,t) *∂*x 2 ) is presented in Figure 5c, which demonstrates a higher average strain distribution of 4OPT-1 than that of Ref, leading to an improvement in the output performance of 4OPT-1. However, it is noted that the effects of a reduction in material properties such as piezoelectricity with the decrease in the thickness of PZT are neglected theoretically. Although the output performance of 4OPT-1 may be augmented without the consideration of the effects of material properties, the optimal result (4OPT-1) given by the proposed data-driven optimization strategy provides the trend of improving output performance by tuning geometric parameters.

Furthermore, a comparison between the GPS algorithm and GA for the four-variable optimization problem has been implemented and is depicted in Figure 6, which shows the convergence curve versus the number of function evaluations using the GPS algorithm and GA with different population size. The parameters used in GA are the default settings S1

(x, y, t) = −y

∂ <sup>2</sup>z(x,t) ∂x 2

in MATLAB; that is, the mutation method and crossover method are constraint dependent, and the selection method is stochastically uniform. Each configuration was limited to run approximately 250 times for comparison.

, ,

–

mm and 4.6 μm, respectively. In addition, the

**Figure 5.** Evaluation of the proposed data-driven optimization strategy from various perspectives. (**a**) Trajectory of the parameters during optimization. (**b**) Output power versus external resistance for previous reported high-output PEH and optimized structure. (**c**) Strain distributions along the arc length of previous reported high-output PEH and optimized structure.

**Figure 6.** Comparison among the efficiency of generalized pattern search algorithm and genetic algorithm with different population size.

To some extent, the number of function evaluations can reflect the running time due to having the same solver (FEM simulation) for each configuration. Although the genetic algorithm can obtain a higher normalized output power at the beginning, the GPS algorithm succeeds in achieving better performance than the GA algorithm with a population size of 10 or 20 after evaluation 250 times. Furthermore, the GPS algorithm requires fewer function evaluation than the GA. The reasons why the pattern search algorithm is more efficient than the GA algorithm with a population size of 10 or 20 are numerous, including the characteristics of the optimization problem, the complex parameter tuning of GA algorithm, and the size of parameter space. Theoretically, the performance of GA may be improved by increasing the population size, whereas the running time may increase for the proposed optimization problem in this work, possibly leading to a low efficiency of optimization.

Figure 7a shows a scatterplot among the geometric parameters *lp*, *w*, *tp*, *t<sup>b</sup>* , resonant frequency and output power. It is observed that the data are not completely randomly distributed in the parameter space due to the minimum searching characteristic of the GPS algorithm and GA. The peak in each histogram in the first four rows indicates the preferable geometric parameters given by the GPS algorithm and GA that may generate better output performance. The contour of the scatterplot of the power and thickness of PZT, the power and thickness of the structural layers, and the power and resonant frequency may present non-linearity among the aforementioned variables. Based on the data in the scatterplot, Spearman's rank correlation coefficient was calculated using Equation (4) to evaluate the strength and direction of the monotonic relationship among the geometric parameters, resonant frequency and output power.

$$r\_s = 1 - \frac{6\sum d\_i^2}{n(n^2 - 1)}\tag{4}$$

where *d<sup>i</sup>* = *rg*(*Xi*) − *rg*(*Yi*), representing the difference in rank between *X<sup>i</sup>* and *Y<sup>i</sup>* ; *n* is the number of observations; *X<sup>i</sup>* and *Y<sup>i</sup>* are variables to be evaluated.

) Spearman's rank correlation matrix of the geometric parameters and output performance **Figure 7.** (**a**) Scatterplot and (**b**) Spearman's rank correlation matrix of the geometric parameters and output performance (resonant frequency and output power).

The correlation matrix with normalized values is presented in Figure 7b. The greater the absolute value of the correlation coefficient, the stronger the monotonic relationship between the evaluated variables. After filtering the identical data, the amount of data available for calculating the Spearman's correlation coefficient is 4506. The maximum value in the correlation matrix is 0.8744 and the corresponding *p*-value under t-distribution is 0, suggesting a strong positive correlation between the thickness of the structural layers and resonant frequency, which means that the designers may need to reduce the thickness of the structural layers to lower the resonant frequency. The *p*-value under t-distribution is calculated using Equation (5):

$$t = r\_s \sqrt{\frac{n-2}{1-r\_s^2}}\tag{5}$$

' In addition, the coefficient between power and width is 0.7488, suggesting that the designers may increase the width to improve the output power. However, it is noted that the coefficient of resonant frequency and width is −0.3573, indicating a weak negative correlation between them. Expanding the width may increase the resonant frequency. Designers may need to balance the impact of increase on width. Moreover, although the coefficient of power and length of PZT is close to 0, its *p*-value (0.4475) is larger than the

— —

general threshold of 0.05, which indicates that this correlation may fail to reject the null hypothesis. As shown in Figure 7a, although the scatterplot of the power and length of PZT presented no monotonic relationship, the upper contour in the scatterplot showed a non-linearity relationship between them, in which the optimal length of PZT is 6.0 mm when the thickness of PZT is 4.6 µm. In addition, in the scatterplot of the power and the length of PZT, although the length of PZT is set as the optimum (6.0 mm), the output power may encounter poor performance without proper setup of other geometric parameters. Therefore, it is necessary to simultaneously optimize multiple parameters for improving the performance of PEHs using the proposed data-driven optimization strategy.

#### **4. Conclusions**

In summary, a data-driven optimization strategy based on FEM simulation (solver) and a GPS algorithm (optimizer) was proposed and implemented, which can not only optimize the geometric parameters of PEHs to achieve high performance but also save time for analytical modeling and complex parameter tuning. In addition, the proposed strategy can search the variable space more precisely than manually varied geometric parameters at fixed intervals. The effectiveness and efficiency of the proposed scheme were validated by comparing the optimization results with previous works [18,23] and experimental results. Then, a four-variable (*lp*, *w*, *t<sup>p</sup>* and *t<sup>b</sup>* ) optimization was further investigated with the aim of maximum output performance. The optimal ratio of the length of PZT and structural layer increases as the thickness of PZT decreases compared with previous works. The optimization results showed that the optimal length ratio is 0.40 (6.0 mm/15.0 mm) when the thickness of PZT is 4.6 µm, while the optimal length ratio is 0.23 (3.4 mm/15.0 mm) when the thickness of PZT is 50 µm. Furthermore, Spearman's rank correlation coefficient was calculated for evaluating the correlation among geometric parameters and output performance such as resonant frequency and output power, providing a data-based perspective on the design and optimization of piezoelectric energy harvesters.

Our evaluation showed that the GPS algorithm exhibited better performance than GA in terms of efficiency, requiring fewer function evaluations than GA. Solving the multiparameter coupling problem such as that of PEHs may require expensive computation cost. Therefore, utilizing a GPS algorithm can shorten the total optimization time for complex high-dimensional optimization problems in real life. In addition, owing to the extensive application of the FEM simulation, the proposed strategy is not limited to the optimization of PEHs, but can also be used for the optimization of other structures that are challenging to formulate, facilitating the employment of the proposed strategy in other fields.

**Author Contributions:** Conceptualization, Y.H.; methodology, Y.H.; software, Y.H.; validation, Y.H., Z.Y.; formal analysis, Y.H., Z.Y., G.H.; investigation, Y.H., G.H.; resources, B.Y.; data curation, Y.H.; writing—original draft preparation, Y.H.; writing—review and editing, Y.H., Z.Y., G.H., B.Y.; visualization, Y.H., Z.Y.; supervision, B.Y.; project administration, B.Y.; funding acquisition, B.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is funded by the Pre-research Foundation of Equipment (No. 6142601180202) and the Joint Foundation of Pre-research of Equipment and the Ministry of Education (No. 6141A02022637).

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors are also grateful to the Center for Advanced Electronic Materials and Devices (AEMD), Shanghai Jiao Tong University.

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

#### **References**

