﻿ 基于标准粒子群优化算法的压电片位置与尺寸优化
«上一篇
 文章快速检索 高级检索

 应用科技  2017, Vol. 44 Issue (5): 62-69  DOI: 10.11991/yykj.201610005 0

### 引用本文

WEI Haixia, SHEN Jianbo, WANG Hongtao. Location and size optimization of piezoelectric patch based on standard PSO algorithm[J]. Applied Science and Technology, 2017, 44(5), 62-69. DOI: 10.11991/yykj.201610005.

### 文章历史

Location and size optimization of piezoelectric patch based on standard PSO algorithm
WEI Haixia, SHEN Jianbo, WANG Hongtao
College of Mechanical and Electrical Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
Abstract: The location and size of the piezoelectric patch on the cantilever for vibration energy harvester are the important factors that influence the recovery of the electric energy. A kind of novel method using the standard particle swarm optimization (PSO) algorithm was put forward for optimizing the location and size of the piezoelectric patch. Firstly, the energy equation of the cantilevered vibration energy harvester was analyzed and deduced. Secondly, the optimization program of standard PSO algorithm which takes the energy equation as the objective function was compiled using MATLAB software, so that the optimal location and size of the piezoelectric patch for the first three orders of modes were obtained by running the optimization program. Lastly, the natural frequency and local point voltage of the cantilevered piezoelectric vibration energy harvester for the first three orders of modes were computed by using the Abaqus software, and the open-circuit energy can be calculated using the energy equation, and then the optimal location and size of the piezoelectric patch can be obtained by analyzing the local point voltage and open-circuit energy. The results of the optimal location and size of the piezoelectric patch from implementing the standard PSO algorithm are consistent with the analysis results from Abaqus software. The detailed results are as follows that the optimal location is at the root of the cantilever and the size is about half length of the cantilever under the first mode, and the optimal location is at the middle and the size is about half the length of the beam under the second mode, and they are the two-thirds and about one-third beam length under the third mode.
Key words: cantilever    vibration energy harvesting    piezoelectric patch    location optimization    size optimization    energy equation    standard particle swarm optimization algorithm    simulation analysis

1 能量方程

 图 1 压电悬臂梁的结构示意

 $\begin{array}{*{20}{c}} {\rho A\frac{{{\partial ^2}u\left( {x,t} \right)}}{{\partial {t^2}}} + c\frac{{{\partial ^5}u\left( {x,t} \right)}}{{\partial {x^4}\partial t}} + EI\frac{{{\partial ^4}u\left( {x,t} \right)}}{{\partial {x^4}}} = }\\ {f\left( {x,t} \right)\delta \left( {x - L} \right)} \end{array}$ (1)

 $u\left( {x,t} \right) = \sum\nolimits_{k = 1}^\infty {{\eta _k}\left( t \right){U_k}\left( x \right)}$ (2)

Uk(x)具有如下表达式：

 ${U_k}\left( x \right) = \cosh \left( {{s_k}x} \right) - \cos \left( {{s_k}x} \right) + {v_k}\left( {\sinh \left( {{s_k}x} \right) - \sin \left( {{s_k}x} \right)} \right)$ (3)

 ${v_k} = - \frac{{\sinh \left( {{s_k}L} \right) - \sin \left( {{s_k}L} \right)}}{{\cosh \left( {{s_k}L} \right) - \cos \left( {{s_k}L} \right)}}$ (4)

sk可以通过悬臂梁频率方程得到：

 $\cos \left( {{s_k}L} \right)\cosh \left( {{s_k}L} \right) = - 1$ (5)

 $\begin{array}{*{20}{c}} {\rho A\sum\limits_{k = 1}^\infty {{{\ddot \eta }_k}\left( t \right){U_k}\left( x \right)} + c\sum\limits_{k = 1}^\infty {{{\dot \eta }_k}\left( t \right)U_k^{\left( 4 \right)}\left( x \right)} + }\\ {EI\sum\limits_{k = 1}^\infty {{\eta _k}\left( t \right)U_k^{\left( 4 \right)}\left( x \right)} = {f_0}\sin \left( {\omega t} \right)\delta \left( {x - L} \right)} \end{array}$ (6)

 $\begin{array}{*{20}{c}} {\rho A\sum\limits_{k = 1}^\infty {{{\ddot \eta }_k}\left( t \right){U_k}\left( x \right)} + c\sum\limits_{k = 1}^\infty {{{\dot \eta }_k}\left( t \right)s_k^4{U_k}\left( x \right)} + }\\ {EI\sum\limits_{k = 1}^\infty {{\eta _k}\left( t \right)s_k^4{U_k}\left( x \right)} = {f_0}\sin \left( {\omega t} \right)\delta \left( {x - L} \right)} \end{array}$ (7)

 $\begin{array}{l} \int_0^L {\rho A\sum\limits_{k = 1}^\infty {{{\ddot \eta }_k}\left( t \right){U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} } + \\ \int_0^L {c\sum\limits_{k = 1}^\infty {{{\dot \eta }_k}\left( t \right)s_k^4{U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} } + \\ \int_0^L {EI\sum\limits_{k = 1}^\infty {{\eta _k}\left( t \right)s_k^4{U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} } = \\ \int_0^L {{f_0}\sin \left( {\omega t} \right)\delta \left( {x - L} \right){U_n}\left( x \right){\rm{d}}x} \end{array}$ (8)

 $\begin{array}{*{20}{c}} {\int_0^L {\rho A\sum\limits_{k = 1}^\infty {{{\ddot \eta }_k}\left( t \right){U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} } = }\\ {\rho A{{\ddot \eta }_k}\left( t \right)\left[ {\int_0^L {{U_1}\left( x \right){U_n}\left( x \right){\rm{d}}x} + \int_0^L {{U_2}\left( x \right){U_n}\left( x \right){\rm{d}}x} + } \right.}\\ { \cdots + \int_0^L {{U_{k - 1}}\left( x \right){U_n}\left( x \right){\rm{d}}x} + \int_0^L {{U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} + }\\ {\left. {\int_0^L {{U_{k + 1}}\left( x \right){U_n}\left( x \right){\rm{d}}x} + \cdots } \right]} \end{array}$ (9)

 $\left\{ \begin{array}{l} \int_0^L {{U_k}\left( x \right){U_n}\left( x \right) = 0} ,n \ne k\\ \int_0^L {{U_k}\left( x \right){U_n}\left( x \right) \ne 0} ,n = k \end{array} \right.\left( {n = 1,2,3, \cdots ,k, \cdots } \right)$ (10)

 $\int_0^L {\rho A\sum\limits_{k = 1}^\infty {{{\ddot \eta }_k}\left( t \right){U_k}\left( x \right){U_n}\left( x \right){\rm{d}}x} } = \rho A{{\ddot \eta }_k}\left( t \right)\int_0^L {U_k^2\left( x \right){\rm{d}}x}$ (11)

 $\begin{array}{*{20}{c}} {\rho A{{\ddot \eta }_k}\left( t \right)\int_0^L {U_k^2\left( x \right){\rm{d}}x} + cs_k^4{{\dot \eta }_k}\left( t \right)\int_0^L {U_k^2\left( x \right){\rm{d}}x} + }\\ {EIs_k^4{\eta _k}\left( t \right)\int_0^L {U_k^2\left( x \right){\rm{d}}x} = {f_0}\sin \left( {\omega t} \right){U_k}\left( L \right)} \end{array}$ (12)

 ${{\ddot \eta }_k}\left( t \right) + 2{\xi _k}{\omega _k}{{\dot \eta }_k}\left( t \right) + \omega _k^2{\eta _k}\left( t \right) = {F_k}\sin \left( {\omega t} \right)$ (13)

 $\eta _k^ * = \frac{{{F_k}}}{{\omega _k^2\sqrt {{{\left( {1 - {{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}} \right)}^2} + 4\zeta _k^2{{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}} }}$ (14)
 ${\varphi _k} = \arctan \left( {\frac{{2{\xi _k}\left( {\frac{\omega }{{{\omega _k}}}} \right)}}{{1 - {{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}}}} \right)$ (15)

 $\begin{array}{*{20}{c}} {u\left( {x,t} \right) = \sum\limits_{k = 1}^\infty {\frac{{{F_k}}}{{\omega _k^2\sqrt {{{\left( {1 - {{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}} \right)}^2} + 4\xi _k^2{{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}} }}} \cdot }\\ {\sin \left( {\omega t + \arctan \left( {\frac{{2{\xi _k}\left( {\frac{\omega }{{{\omega _k}}}} \right)}}{{1 - {{\left( {\frac{\omega }{{{\omega _k}}}} \right)}^2}}}} \right)} \right){U_k}\left( x \right)} \end{array}$ (16)

ω=ωk时，梁体产生共振，此时压电悬臂梁吸收外界的振动能、产生较高的电能。梁体此时的变形主要为该阶固有频率所对应的主振型，其他各阶振型的影响可忽略不计。因此，压电悬臂梁的横向振动位移即为压电片的应变，由式(2)可知，压电片的应变近似表达式为

 $u\left( {x,t} \right) \approx {\eta _k}\left( t \right){U_k}\left( x \right)$ (17)

 $\varphi = - \frac{{bh}}{S}\int_x {\left( {{h_{31}}r\frac{{{\partial ^2}u\left( {x,t} \right)}}{{\partial {x^2}}}} \right){\rm{d}}x}$ (18)

ω=ωk时，将式(17)、式(3)代入式(18)，得开路电压表达式为

 $\varphi = - \frac{{bh}}{S}\int_x {\left( {{h_{31}}r\frac{{{\partial ^2}\left( {{\eta _k}\left( t \right){U_k}\left( x \right)} \right)}}{{\partial {x^2}}}} \right){\rm{d}}x} = \\ - \frac{{bh{h_{31}}rs_k^2{\eta _k}\left( t \right)}}{S}\int_x {\left[ {\cosh \left( {{s_k}x} \right) + \cos \left( {{s_k}x} \right) + {v_k}\left( {\sinh \left( {{s_k}x} \right) + \sin \left( {{s_k}x} \right)} \right)} \right]{\rm{d}}x}$ (19)

 ${\varphi ^ * } = \frac{{bh{h_{31}}rs_k^2\eta _k^ * }}{S}\int_x {\left[ {\cosh \left( {{s_k}x} \right) + \cos \left( {{s_k}x} \right) + {v_k}\left( {\sinh \left( {{s_k}x} \right) + \sin \left( {{s_k}x} \right)} \right)} \right]{\rm{d}}x}$ (20)

 $Q = \frac{1}{2}C{\varphi ^2}$ (21)

 ${Q^ * } = \frac{1}{2}{C^{\left( {\varphi * } \right)}}2 = \frac{1}{2}\frac{{{\varepsilon _{33}}bhh_{31}^2{r^2}s_k^4{{\left( {\eta _k^ * } \right)}^2}}}{{{x_2} - {x_1}}} \cdot\\ {\left[ {\int_{{x_1}}^{{x_2}} {\left( {\cosh \left( {{s_k}x} \right) + \cos \left( {{s_k}x} \right) + {v_k}\left( {\sinh \left( {{s_k}x} \right) + \sin \left( {{s_k}x} \right)} \right)} \right){\rm{d}}x} } \right]^2}$ (22)
2 位置和尺寸优化

 ${q_k}\left( {{x_1},{x_2}} \right) = \frac{1}{{{x_2} - {x_1}}}{\left[ {\int_{{x_1}}^{{x_2}} {\left( {\cosh \left( {{s_k}x} \right) + \cos \left( {{s_k}x} \right) + {v_k}\left( {\sinh \left( {{s_k}x} \right) + \sin \left( {{s_k}x} \right)} \right)} \right){\rm{d}}x} } \right]^2}$ (23)

2.1 标准粒子群优化算法

 $\begin{array}{*{20}{c}} {{V_{i,j}}\left( {t + 1} \right) = w \cdot {V_{i,j}}\left( t \right) + {c_1} \cdot {r_{1,i,j}}\left( t \right) \cdot }\\ {\left( {{P_{i,j}}\left( t \right) - {X_{i,j}}\left( t \right)} \right) + {c_2} \cdot {r_{2,i,j}}\left( t \right) \cdot \left( {{G_{i,j}}\left( t \right) - {X_{i,j}}\left( t \right)} \right)} \end{array}$ (24)
 ${X_{i,j}}\left( {t + 1} \right) = {V_{i,j}}\left( {t + 1} \right) + {X_{i,j}}\left( t \right)$ (25)

 $w\left( t \right) = \frac{{\left( {{w_{\max }} - {w_{\min }}} \right) \cdot \left( {{t_{\max }} - t} \right)}}{{{t_{\max }}}} + {w_{\min }}$

1) 参数初始化：设定粒子群规模M，粒子空间维数N；选取最大迭代次数tmax；学习因子为2；粒子飞行最大速度，粒子位置变化范围；设定t=0，在N维空间中随机选择粒子的位置与速度；选取最大惯性权重和最小惯性权重。

2) 计算种群中每个粒子当前位置Xi(t)处的适应度函数值，将适应度函数依然定义为式(23)。针对每个粒子，比较当前位置Xi(t)对应的适应度函数值和其个体历史最优位置Pi(t)对应的适应度函数值，优则替换Pi(t)的适应度函数值，并有Pi(t)=Xi(t)，反之则保持不变。

3) 针对各个粒子，分析个体历史最优位置Pi(t)和群体历史最优位置G(t)，并比较两者所对应的适应度函数值，优则替换G(t)的适应度函数值，并有G(t)=Pi(t)，反之则保持不变。

4) 根据进化方程(24)、(25)分别对所有粒子的位置与速度进行更新。

5) 若未达到最大进化代数tmax，则置t=t+1，返回步骤2)继续执行，若达到最大进化代数，则以得到的具有最大适应度的个体作为最优解输出。

2.2 算例

 ${q_1}\left( {{x_1},{x_2}} \right) = \\ \frac{1}{{\left( {{x_2} - {x_1}} \right)}}{\left[ {\int_{{x_1}}^{{x_2}} {\left( {\cosh \left( {18.751x} \right) + \cos \left( {18.751x} \right) + 0.7341\left( {\sinh \left( {18.751x} \right) + \sin \left( {18.751x} \right)} \right)} \right){\rm{d}}x} } \right]^2}$ (26)
 ${q_2}\left( {{x_1},{x_2}} \right) =\\ \frac{1}{{\left( {{x_2} - {x_1}} \right)}}{\left[ {\int_{{x_1}}^{{x_2}} {\left( {\cosh \left( {46.941x} \right) + \cos \left( {46.941x} \right) + 1.0185\left( {\sinh \left( {46.941x} \right) + \sin \left( {46.941x} \right)} \right)} \right){\rm{d}}x} } \right]^2}$ (27)
 ${q_3}\left( {{x_1},{x_2}} \right) = \\ \frac{1}{{\left( {{x_2} - {x_1}} \right)}}{\left[ {\int_{{x_1}}^{{x_2}} {\left( {\cosh \left( {78.548x} \right) + \cos \left( {78.548x} \right) + 0.9992\left( {\sinh \left( {78.548x} \right) + \sin \left( {78.548x} \right)} \right)} \right){\rm{d}}x} } \right]^2}$ (28)

 图 2 目标函数最优解随迭代次数的变化
3 仿真分析

 图 3 100 mm梁前三阶模态下局部点电压

1) 仿真分析并计算一阶模态下的能量，将长分别为5、10、20、30、40、50、60、70、80、90 mm，宽与厚为10 mm与0.1 mm的压电片从悬臂梁根部开始分别贴于10个相同规格的梁上，在梁自由端施加幅值1 mm、频率f1的正弦激励，通过仿真分析得到这10个压电片在激振力作用下各自产生的开路电压，再由式(22)计算得到各压电能量值，结果如图 4所示。

 图 4 100 mm梁一阶模态下的压电能量

2) 仿真计算二阶模态下的能量，按压电片贴于梁根部与梁中部两种情况讨论。当压电片贴于悬臂梁根部时，将长分别为5、10、20、30、40、50、60、70、80、90 mm的压电片从梁根部开始贴于10个梁上，在梁自由端施加幅值1 mm、频率f2的正弦激励，通过仿真分析与计算得到模态二下这10种长度压电片贴于梁根部时的压电能量，结果如图 5中Q21所示；当压电片贴于梁中部时，即压电片中心置于悬臂梁中部，将这10种长度的压电片分别贴于10个悬臂梁上，在梁自由端施加幅值为1 mm、频率为f2的正弦激励，通过仿真分析与计算得到模态二下这10种长度压电片贴于梁中部时的压电能量，结果如图 5中Q22所示。

 图 5 100 mm梁二阶模态下的压电能量

3) 仿真计算三阶模态下的能量，分压电片贴于悬臂梁根部、梁长的三分之一和梁长的三分之二处3种情况讨论。当压电片贴于梁根部时，将长分别为5、10、20、30、40、50、60 mm的压电片从梁根部开始贴于7个规格相同的梁上，在梁自由端施加幅值1 mm、频率f3的正弦激励，通过仿真分析与计算得到模态三下这7种长度压电片贴于梁根部时的压电能量，结果如图 6中Q31所示；当压电片贴于梁长的三分之一处时，即压电片中心置于距梁固定端35 mm处，将这7种长度的压电片贴于7个梁上，在梁自由端施加幅值1 mm、频率f3的正弦激励，通过仿真分析与计算得到模态三下这7种长度压电片贴于梁长三分之一处时的压电能量，结果如图 6中Q32所示；当压电片贴于梁长的三分之二处时，即压电片中心置于距梁固定端70 mm处，将这7种长度的压电片贴于7个梁上，在梁自由端施加幅值1 mm、频率f3的正弦激励，通过仿真分析与计算得到模态三下这7种长度压电片贴于梁长三分之二处时的压电能量，结果如图 6中Q33所示。

 图 6 100 mm梁三阶模态下的压电能量

4 结论

1) 标准PSO算法与Abaqus软件仿真的结果基本吻合：一阶模态下，压电片最优贴片位置位于悬臂梁根部，最优尺寸约为梁长一半；二阶模态下，压电片最优贴片位置位于梁中部，最优尺寸约为梁长一半；三阶模态下，压电片最优贴片位置位于梁的三分之二处，最优尺寸约为梁长的三分之一；

2) 结果验证了运用标准PSO算法实现压电振动能量回收系统中压电片位置和尺寸优化方法的有效性。

 [1] 文晟, 张铁明, 刘旭, 等. 基于压电效应的振动能量回收装置的研究进展[J]. 机械科学与技术, 2010, 29(11): 1515-1520. (0) [2] 边义祥, 杨成华. 基于压电材料的振动能量回收技术现状综述[J]. 压电与声光, 2011, 33(4): 611-622. (0) [3] 阚君武, 唐可洪, 王淑云, 等. 压电悬臂梁发电装置的建模与仿真分析[J]. 光学精密工程, 2008, 16(1): 71-75. (0) [4] 袁江波, 谢涛, 陈维山, 等. 悬臂梁压电发电装置的实验研究[J]. 振动与冲击, 2009, 28(7): 69-72. (0) [5] LIANG Zhu, XU Chundong, REN Bo, et al. Optimization of cantilevered piezoelectric energy harvester with a fixed resonance frequency[J]. Technological sciences, 2014, 57(6): 1093-1100. DOI:10.1007/s11431-014-5556-7 (0) [6] SODANO H A, INMAN D J, PARK G. Comparison of piezoelectric energy harvesting devices for recharging batteries[J]. Journal of intelligent material systems and structures, 2005, 16(10): 799-807. DOI:10.1177/1045389X05056681 (0) [7] 刘树林, 许小勇, 翟宇毅, 等. 振动模态对压电发电机陶瓷片粘贴位置的影响[J]. 光学精密工程, 2011, 19(8): 1801-1809. (0) [8] LI Hua, HU Shundi, TZOU H S. Size optimization of conical piezoelectric energy harvester[C]//Proceedings of 2011 Symposium on Piezoelectricity, Acoustic Waves and Device Applications (SPAWDA). Shenzhen, China, 2011:485-488. (0) [9] 高瑞贞, 张京军, 郑骥, 等. 基于改进遗传算法主动柔性结构压电元件位置优化[J]. 计算力学学报, 2008, 25(4): 542-546. (0) [10] DHURI K D, SESHU P. Multi-objective optimization of piezo actuator placement and sizing using genetic algorithm[J]. Journal of sound and vibration, 2009, 323(3/4/5): 495-514. (0) [11] HADAS Z, KURFURST J, ONDRUSEK C, et al. Artificial intelligence based optimization for vibration energy harvesting applications[J]. Microsystem technologies, 2012, 18(7/8): 1003-1014. (0) [12] 潘继, 蔡国平. 桁架结构作动器优化配置的粒子群算法[J]. 工程力学, 2009, 26(12): 35-39. (0) [13] 潘继, 陈龙祥, 蔡国平. 柔性板压电作动器的优化位置与主动控制实验研究[J]. 振动与冲击, 2010, 29(2): 117-120. (0) [14] 马天兵, 裘进浩, 季宏丽, 等. 基于粒子群的压电结构多目标同步优化控制[J]. 沈阳工业大学学报, 2012, 34(5): 569-575. (0) [15] 胡海岩. 机械振动基础[M]. 北京: 北京航空航天大学出版社, 2005: 38-1188. (0) [16] TZOU H S. Piezoelectric shells:distributed sensing and control of continua[M]. London: Kluwer Academic Publishers, 1993: 146. (0) [17] HU Shundi, CHUANG K C, TZOU H S. PVDF energy harvester on flexible rings[C]//Proceedings of 2010 Symposium on Piezoelectricity, Acoustic Waves and Device Applications (SPAWDA). Xiamen, China, 2010:100-105. (0)