杜小振, 柳叙强, 王国全, 朱 尚,赵 岩
(山东科技大学 机械电子工程学院,山东 青岛 266590)
振动水柱压电功能材料发电研究需完成两项关键技术:
1) 振荡水柱气室气压分析。
2) 压电功能材料形变与输出电压分析。
在振荡水柱气室仿真分析基础上对气室水动力性能进行研究,并对气室结构进行改造,加以实验验证。Luo Y等[1]结合气室结构参数和气孔阻尼系数分析结果表明,前墙厚度和气室高度对波能捕获效率影响较大。宁德志等[2]建立非线性数值模型得出,气室内捕获的波浪力与入射波成正比增加。陈帆等[3]设计了一种双圆筒沉箱结构的振荡水柱式波能发电装置,对不同参数沉箱模型的波能转换和吸收效果进行了试验研究,优化了设计参数。李宏伟[4]总结了源造波法与摇板造波法规律,研究了不同方法下振荡水柱二维数值模型的区别。Mahnamfar F等[5]通过仿真和实验建模对比数据设计气室的几何参数,提高了波能转换效率。Elhanafia A[6]研究了不同海浪状况下波浪力与振荡水柱气室的相互作用,发现波浪力和气室内气动力直接影响共振频率下气室出口的垂直气动力。压电发电研究主要根据变形与电压关系对优化压电模型结构和分析发电特性。许伟斌[7]利用气室内输气管往复气流推动铷磁铁往复运动,采用直接式电磁感应发电机作为振荡水柱波能发电系统的发电装置,并试验分析了设备可行性。王淑云等[8]提出一种基于圆弧限位压电发电装置,并制作模型,实验确定压电振子最大输出电压仅与限位圆弧半径有关,半径越小发电量越高。K.F. Wang等[9]基于哈密尔顿原理对单晶压电能量收集器的挠曲电效应分析,结果表明,较小的压电层厚度挠曲形变对输出电压和功率产生显著影响。Hajhosseini M等[10]应用广义差分正交法计算出周期性变截面的双压电晶片的电压输出,并与普通压电梁对比得出,在宽频范围内产生电压更优异。S Z Li等[11]在不同的预应力机械振动下研究压电换能器的能量收集性能,结果表明,换能器在共振频率转移到较低值时,捕获能量随着预应力的增加而增加。欧阳嘉艺等[12]提出了浮子结构形式、弹簧力与阻尼力共同作用下的结构、压电振子三方面的优化方案,指出压电材料的机电耦合优势。顾亚雄等[13]设计开发了一种可收集环境中多方向振动的压电发电装置,并进行有限元分析和实验测试,研究提高压电发电环境振动能量收集效率。当前仅对振荡水柱气室和压电发电模型分别研究。本文将压电发电特性应用于振荡水柱气室结构模型中,通过FLUENT软件内置用户自定义函数(UDF)功能建立二维数值波浪水槽模型,模拟数值水槽造波及波浪进入振荡水柱气室,并产生气室气压的波能收集过程,将气室气动力作用于气室出口处的压电片,利用仿真与理论计算模型研究振荡水柱压缩空气进而驱动压电材料形变发电过程,探索能量输出特性。
为实现模拟实验造波采用推板造波法,在造波端模拟推板往复运动,驱动水槽内水体造波。推板造波的造波效果好,且便于调节波浪波高、波长。利用FLUENT软件的动网格及UDF二次开发,编UDF程序,造波板速度u取值冲程[14]为
(1)
推板运动波面方程
(2)
式中:η为波面;T为周期;h为水深;t为时间;S为推板振幅;x为波浪运动方向坐标;H为波高;Ks(i) 是ω2=-gKstan(KsH) 的第i个根。
-ηlsin(Kx-ωt)
(3)
根据式(3)确定S和波浪振幅ηl之间的关系为
(4)
为减小水槽终端产生的反射波干扰造波在小波区采用数值沙滩法消波。设置开启FLUENT系统中的数值模拟(Numerical Beach)功能设置自由液面标高为0,河床标高为-2.2 m,并设定波浪源等。以水槽末端为消波结束点,设置为60,设置消波波长为20,波数2,该方法可减少数值分析时间,方便调节消波效果。
由于所建立压电片模型为固定边界的圆形发电结构,气室内空气压强作用于压电圆片后在电极面上形成自由电荷形成电势。机电耦合过程中自变量是应变W及电场强度E,因变量是应力O及电位移D,选用机械夹持和电学短路边界条件为基础的第二类压电方程组数学表达式为
O=cEW-eiE
(5)
D=eW+εWE
(6)
式中:cE为弹性刚度系数;e为压电应力系数;ei为e的转置;εW为加紧压电常数。
第n次正压电效应压电片产生的电压值为
(7)
加载外力Fj时,压电片厚度方向的变形量为
(8)
式中:gni为压电系数;Kp为弹性系数;A为压电片表面积;hp为压电片厚度。
造波、消波为数值波浪水槽模拟提供理论基础,计算相应波长、波高、波幅,确定压电片材料参数,分别在FLUENT和ANSYS系统的APDL中进行设置。利用ICEM CFD建立二维数值水槽模型如图1所示,水槽长60 m、宽7.1 m,水深2.2 m,气室距离左侧水槽壁25 m,气室右侧水域为消波区。选取海平面为0坐标点,在气室内部创建监测点a检测气室压力。
图1 振荡水柱发电二维数值水槽模型
在CAE前处理软件ICEM中进行网格划分,为简化网格划分步骤对气室水槽和压电片分别划分网格。气室水槽部分采用结构化网格,划分为气室、水、空气3个区域,为减少仿真时间提高仿真精度,设置网格密度气室>水>空气。压电片部分应用APDL指令进行三维六面体映射网格划分,设置整体线单元边长为0.001 m,压电片厚度方向网格划分为6段。
设置数值水槽的空气面上边界为压力入口,左壁面为速度入口,底面和右壁面为无滑移壁面边界条件,空气密度为1.225 kg/m3,海水密度为998.2 kg/m3,波浪湍流参数:湍动能为0.001 024 m2/s3,湍动耗散率为0.000 025 799 m2/s3。标记水相区域,设置压强监测点及速度检查点。气室模型数值水槽压力云图和气室速度云图分别如图2、3所示。由图可看出,出口速度在波浪压缩空气作用下变化较大,气室内气体流速变化存在梯度;从水槽底部到气室内空气部分的压力逐渐降低,模拟结果符合现场海况。
图2 数值水槽压力分布云图
图3 气室速度分布云图
压电圆片受力输出能量模拟分析需考虑其厚度和半径,取气孔有效面积为压电片的受力发电部分,在APDL中建立三维压电片模型,其半径为0.1 m,厚度为0.001 m,选用单元类型为Solid5,设置KEYOPT(1)=3激活压电自由度、位移和电压,压电材料密度为7 600 kg/m3,弹性模量为1.69×105,泊松比为0.066,压电片边缘为固支约束,压电功能材料的介电矩阵、压电应力矩阵及弹性系数矩阵分别为
(9)
(10)
(11)
求解方法为完全瞬态求解分析法,载荷设置为斜坡载荷,关闭自动时间步控制,应用缺省子步数和时间步长。
水槽气室收集的波浪能在气室内形成压强,作用在气室出口的压电片上引起压电片变形,其受力云图如图4所示。
图4 压电片变形受力分布
采用推板造波模拟振荡水柱(OWC)气室数值水槽内波浪能转化为气室空气压力的变化情况,监测点a检测到空气压力变化曲线如图5所示。气室内前期压力有2个峰值,压强p变化较大,后期趋于稳定,气室波高起伏一致且呈现拟正弦变化,波浪到达气室9 s时,在气室惯性力及气室内壁发射波作用下有较大的起伏,14 s左右达到共振,此时压强约为13 Pa,稳定后整体压强在5 Pa左右浮动。分析可知,在波浪到达气室前期(14 s前)气室内壁反射波对气室气压影响较大,可利用前期两处峰值气压提高波能转换效率,改善装置性能。
图5 监测点a气室压强与正弦拟合曲线对比
气室内波浪上、下波动引起气室气压变化,在气室出口处驱动压电片上、下往复运动变形,产生电压输出电能。模拟分析过程是将空气压强结果作为压电片下表面的动载荷,作瞬态分析。因为APDL中施加的变动载荷需采用公式化的规律载荷,因此,提取气室压强仿真结果进行公式化处理拟合为正弦曲线(见图5),拟合式为
p=-1.23+9.26cos(1.36t-3.63)
(12)
由于气室气压在14 s前变化较大,导致在负方向整体拟合气压数据偏高约-13 Pa,正方向趋于稳定值(约4.8 Pa)。
将式(12)导入到APDL指令,得到压电片变形产生的电压变化量,如图6所示。输出最大电压量约为0.47 V,压电片的最大变形量约为3 m。
图6 压电片电压变化图
由图5、6可知,受到气室内壁反射波及惯性力影响,整体拟合气室气压在负方向偏大。为进一步提高仿真准确性,优化数据取气室压强14 s后稳态压强数据拟合正弦曲线。拟合正弦压强式为
p=-0.044 6+8.24cos(-1.35t+65.59)
(13)
稳定后气室气压及其正弦拟合曲线对比如图7所示。
图7 气室压强、正弦拟合曲线对比
导入变载荷式重新计算模拟得到压电片压强变化如图8所示。电压最大值约为0.25 V,最大变形量为1.8 m,压电发电部分输出能量相对较低,但此为波浪稳定时压电片发电输出能量,反映了系统稳态输出特性。
图8 压电发电输出电压变化图
本文基于压电效应对振荡水柱波浪能发电系统进行模拟,合流体动力学理论在FLUENT分析系统中建立气室水槽模型,并实现模拟造波。分析收集波浪能在气室内转换为压强情况,进而分析气压驱动压电片形变产生电能。从波浪能到气室气压的转换过程有2个气压峰值分别由气室内壁反射波和共振作用产生,可利用此大幅度气压变化提高装置能量转换效率;模拟输出电压幅值为0.2~0.4 V。分析过程为新型微型压电发电电源的设计提供技术依据。