基于联合仿真的波浪能发电装置PTO系统参数优化研究

何广华，杨凯博，栾政晓，张志刚，刘朝纲，荆芃霖

（1.哈尔滨工业大学（威海）海洋工程学院，山东威海264209；2.山东船舶技术研究院，山东威海264209；3.哈尔滨工业大学机械工程学院，哈尔滨150001）

0 引 言

由于能源消耗和环境污染等问题，海洋可再生能源的开发和利用逐渐受到各国的关注[1]。其中波浪能以储量大、无污染、可重复开发利用等优点，成为国内外海洋能开发利用研究的热点。波浪能发电装置根据其工作原理可分为振荡体式、振荡水柱式、聚波越浪式等[2]。其中，振荡体式又可分为振荡浮子式（点吸收式）、筏式和摆式。近年来，振荡浮子式WEC 因其体积小、制造安装方便、不受波浪方向影响等优势受到了广泛的关注[3]。常见的振荡浮子式WEC 的PTO 系统可分为机械式和液压式，其中液压式PTO 利用液压缸和液压泵等元件将液压能转化成电能，机械式PTO 利用齿轮齿条、滚珠丝杠[4-5]等机械装置将机械能转化成电能。

波浪能发电装置的PTO 系统是波浪能发电的关键技术之一，不同研究者对PTO 系统进行了不同程度的简化研究。Rico 等[6-8]利用数学公式计算了摇臂圆柱形浮子的运动学方程，分析了线性PTO 阻尼对浮子俘获功率的影响；张亚群等[9]利用物模试验分析了振荡浮子的水动力性能，实验显示在最优PTO 负载下浮子可获得最大俘获宽度比；宋文杰等[10]利用实验分析了双行程工作形式PTO 的输出特性，但未对单行程工作形式进行对比分析；冯辉等[11]基于Simulink和AMESim 搭建了联合仿真模型，该模型包括了各液压元件与发电机，较好地模拟出了能量转换系统的发电特性曲线，但其未根据仿真结果对PTO 元件进行优化选型；Li等[12]利用CFD 工具研究了波物相互作用与机械元件对浮子的影响，发现机械元件产生的力和力矩对单个浮子的运动响应影响巨大；刘延伟等[13]对超越离合器进行了理想建模，并分析了在切换瞬间超越离合器的动态特性。

然而，对于采用机械式PTO 的振荡浮子波浪能发电装置，PTO 系统的传动方式、发电机选型及电机的电子负载等参数都会对波浪能发电装置的获能产生影响，导致将PTO 负载简化为线性阻尼的方法不够精确。针对该问题，本文采用水动力模型与PTO系统联合仿真的方法，使所添加的负载更接近真实情况，其中，在STAR-CCM+中建立波浪能发电装置数值模型，在AMESim 中建立PTO 系统仿真模型，通过二者的联合仿真进行水动力与PTO系统的耦合研究；研究PTO负载的不同行程工作模式与不同传动比对浮子运动响应和输出功率的影响，并针对目标海况对PTO 参数进行优化。本文使用的联合仿真方法可应用于机械式PTO 的振荡浮子式波能装备设计，为不同海况与各类振荡浮子式波浪能发电装置的机械参数优化与选型工作提供参考。

1 浮子水动力模型的建立

1.1 控制方程

将水视为粘性不可压缩流体，其控制方程为连续性方程和纳维-斯托克斯（N-S）方程，可表示为

式中，vx、vy和vz分别表示x、y和z方向的速度分量，fx、fy和fz分别表示x、y和z方向受到的质量力，∇2为拉普拉斯算子，t表示时间，ρ表示流体密度，p表示流体压力，ν表示流体的运动粘度。

由于流域包括气、液两相，因此本文以流体体积分数法（VOF）对自由液面进行追踪和捕捉。其中，相的分布和液面的位置由相体积分数描述，相体积分数αi定义为

式中，Vi为网格单元中第i相的体积，V为网格单元总体积。网格单元中所有相的体积分数满足下式：

式中，N为网格单元中相的总个数，文中N=2。

1.2 三维数值波浪水池的建立与验证

采用的数值水池如图1 所示，流域左侧为速度入口边界，右侧为压力出口边界，两侧和底部均为无滑移壁面边界，顶部为速度入口边界。其中上游和下游均设置3倍波长，两侧设置2.5倍波长，四周均设置消波区，沿计算域边界设置了力消波区域[14]。

图1 数值水池模型示意图Fig.1 Schematic diagram of numerical pool model

采用切割体网格进行网格划分。由于涉及到浮子的运动，需采用重叠网格。为了更好地捕捉自由液面，对两相交界处进行加密处理。网格的划分如图2所示。

图2 摇臂浮子网格划分示意图Fig.2 Diagram of float-arm buoy mesh division

为验证数值水池的可行性，对上述建立的三维波浪水池进行网格收敛性验证和时间步收敛性验证。首先选取了三种网格尺寸，具体网格数量与尺寸如表1所示。

表1 网格尺寸与数量Tab.1 Mesh size and quantity

由图3 可以看出，与网格A 的垂荡振幅对比，网格B 和网格C的相对误差分别为1.43%和6.14%，h代表浮子垂荡。综合考虑数值精度和计算成本，最终采用网格B。

图3 网格收敛性验证Fig.3 Grid convergence verification

接下来根据网格B 尺寸的数值模型进行时间步验证。由图4可以看出，时间步Δt=0.02 s时，浮子垂荡曲线明显与Δt=0.01 s 和Δt=0.005 s 有偏移，而Δt=0.01 s 和Δt=0.005 s 的垂荡振幅相对误差为1.28%。故在后续计算中采用Δt=0.01 s的时间步长。

图4 时间步收敛性验证Fig.4 Time step convergence verification

2 PTO系统仿真模型

2.1 波浪能发电装置运行模式简介

本文设计的波浪能发电装置可分为三部分，分别是摇臂浮子、机械传动装置和发电机，如图5所示。

图5 波能发电装置组成示意图（1-摇臂浮子，2-机械传动装置，3-发电机；俯视图）Fig.5 Schematic diagram of WEC compo⁃sition(1-Buoy,2-Mechanical trans⁃mission,3-Motor;top view)

摇臂浮子作为发电装置的波能捕获单元，在波浪激励力的作用下绕铰接轴做纵摇运动，浮子的纵摇通过机械传动装置将运动传递给发电机。其具体的工作流程为：浮子的纵摇通过同步带将运动传递给两根带轮轴，两根带轮轴上分别安装有工作方向相反的超越离合器，单个入射波周期内，两根带轮轴将通过齿轮组交替向发电机的齿轮轴传递动力，经齿轮箱增速后提高发电机的转速。由于浮子的上升与下降行程分别对应两个超越离合器的单程做功，当研究单个行程的做功特性时，也可以在本机械装置上实现。

2.2 PTO系统的仿真模型建立

AMESim 是系统工程高级建模和仿真平台，用以建立复杂的多学科多领域的系统模型，包括流体、机械、电磁和控制等，并在此基础上进行仿真计算和分析，研究元件或系统的稳态与动态性能。本文集成机械、电气与信号元件搭建出PTO 系统的完整模型，并参照仿真结果对系统进行优化。根据上述机械传动装置的工作原理，通过AMESim 搭建了波能发电装置PTO 系统的仿真模型，如图6所示。

图6 波能发电装置PTO系统仿真模型Fig.6 PTO system simulation model of WEC

上述AMESim 仿真模型中，机械部分包括超越离合器、齿轮系、带轮系与齿轮箱。当超越离合器处于接合状态时，将超越离合器视作一个无转动惯量的刚体，根据工作方向传递转矩；当超越离合器处于超越状态时，动力传递中断。以其主动端与从动端的转角差作为模型中工作状态的判断依据，则超越离合器传递转矩M(t)的表达式如下：

式中，Kc为超越离合器扭转刚度，μ为等效粘性摩擦系数，θm(t)为主动端扭转角度，θs(t)为从动端扭转角度。

忽略超越离合器主从动件之间的间隙与变形因素的影响，则机械传动装置的运动方程可表示为

式中，ωb为浮子的角速度，ωM为电机转子的角速度，γ为等效传动比，ig为增速齿轮箱传动比，rb为同步带轮的外径比，rg为齿轮组的分度圆半径比。

改变AMESim 内信号部分的触发方式，分别建立了下降行程做功与上升行程做功两种单行程的动力传递方式，波能发电装置的仿真模型如图7～8所示。

图7 浮子下降行程做功的PTO系统仿真模型Fig.7 PTO system simulation model of downward working buoy

图8 浮子上升行程做功的PTO系统仿真模型Fig.8 PTO system simulation model of upward working buoy

AMESim 电气部分以直流电机作为PTO 系统的发电机，电阻作为PTO 系统的电子负载，用电阻功率来衡量PTO 系统的发电功率。上述三种仿真模型的电气部分与机械部分的参数设置完全相同，只根据做功方式改变了信号部分的触发方式，用以执行联合仿真。

3 水动力-PTO联合仿真设置

本文开展了水动力模型与PTO 系统联合仿真研究，计算出了在PTO 系统作用下，浮子角速度、电机功率与电机转速等相关数据。其仿真过程如图9所示。

图9 STAR-CCM+与AMESim联合仿真示意图Fig.9 Co-simulation diagram of STAR-CCM+and AMESim

AMESim 将PTO 的负载转矩传递到STAR-CCM+中的水动力计算模块，同时考虑了电机转子转动惯量的惯性转矩与电机阻尼的电磁转矩，两者共同组成了电机的负载转矩。电机的负载转矩通过超越离合器和增速齿轮箱等机械传动系统转变为作用在STAR-CCM+模型上浮子的负载力矩。因此，本联合仿真方法相较于直接添加阻尼系数的计算方法更加接近真实PTO系统的工作状态。

本文中的计算模型为带摆臂的圆柱形浮子，所用海况来源于威海国家浅海综合试验场，Star-CCM+内的浮子参数、海况参数设置如表2所示。

表2 Star-CCM+计算模型参数Tab.2 Calculation parameters of the model in Star-CCM+

浮子后端的PTO 系统仿真模型通过AMESim 搭建，发电机选用Z4-100-1 型号直流电机，AMESim内的电机参数、机械元件参数设置如表3所示。

表3 AMESim计算模型参数Tab.3 Calculation parameters of the model in AMESim

下文将往返双行程做功的浮子称为1 号浮子，上升行程做功的浮子称为2 号浮子，下降行程做功的浮子称为3号浮子。图10为双行程做功浮子在一个周期内的运动响应。为研究做功模式以及等效传动比对浮子的角速度与电机功率的影响，对工况进行设置，如表4所示，并根据仿真结果，在特定区间完成PTO系统的负载匹配与优化选型。

表4 工况设置Tab.4 Setting of working conditions

图10 单个波周期内浮子的运动响应Fig.10 Motion response of the buoy in a single wave period

4 结果分析

1 号、2 号和3 号浮子在不同工况下的浮子角速度和电机瞬时发电功率PI如图11～15 所示。由于采用直流电机，故电机瞬时功率与电机转子的角速度成正比。图中，当角速度小于0 时，浮子逆时针转动，对应上升行程做功；角速度大于0时，浮子顺时针转动，对应下降行程做功。

图11 1:20传动比的浮子角速度和电机瞬时功率Fig.11 Buoy angular velocity and instantaneous power under 1:20 transmission ratio

图12 1:40传动比的浮子角速度和电机瞬时功率Fig.12 Buoy angular velocity and instantaneous power under 1:40 transmission ratio

可以看出，在低传动比的情况下（图11～12），三种浮子的电机功率与浮子角速度曲线吻合较好，说明低传动比情况下PTO 负载力矩小，动力传递状态下PTO 系统对浮子运动的影响几乎可以忽略。随着传动比的增大（图13～15），2号浮子对应的1号浮子上升行程段与3号浮子对应的1号浮子下降行程段，角速度曲线前半部分吻合较差，后半部分吻合良好。其原因为：相较于双行程的做功，浮子在单行程做功状态下只有单行程内存在负载力矩，故从自由运动切换至负载状态时的负载力矩发生突变，导致前半段吻合较差，待后半段运动稳定后，功率曲线与角速度曲线吻合良好。在高传动比情况下，随着PTO 系统的负载力矩增大，会显著影响浮子的运动状态，工作行程内浮子角速度的幅值将显著降低。

图13 1:60传动比的浮子角速度和电机瞬时功率Fig.13 Buoy angular velocity and instantaneous power under 1:60 transmission ratio

图14 1:80电机瞬时功率传动比的浮子角速度和电机瞬时功率Fig.14 Buoy angular velocity and instantaneous power under 1:80 transmission ratio

图15 1:100传动比的浮子角速度电机瞬时功率Fig.15 Buoy angular velocity and instantaneous power under 1:100 transmission ratio

图16 为1 号、2 号和3 号浮子在不同传动比工况下一个入射波周期内电机的平均功率图。平均功率PT的计算公式如下：

图16 单个波周期内三种浮子的平均功率Fig.16 Average power of three buoys in a single wave period

从图中可以看出，1 号浮子的输出功率始终大于2 号与3 号浮子，这是因为1 号浮子在一个波周期内对PTO 系统全程做功。在传动比小于40 时，1 号浮子的平均功率是2号和3号浮子的近两倍。但随着传动比的增大，1号浮子的平均功率衰减严重，2号与3号浮子的功率逐渐超过1号浮子的单程功率，在传动比为80到100之间时，1号浮子的输出功率仅为3号浮子的1.3倍。这是因为在高传动比情况下，浮子受到的负载力矩过大，1号浮子在双行程均受到较大的负载，导致其运动受阻，无法完全吸收波浪的能量。而2号和3号浮子，只有一段行程受到了较大的负载力，故单段输出功率较1号浮子更高。

由于计算模型中考虑了摆臂的重力影响，且3号浮子在上升过程不受PTO 的负载力矩，在一个波周期内吸收了更多波浪的能量，使得其输出功率会略大于2号浮子。

此外，随着传动比的增大，三种浮子的平均功率均先上升后下降，说明PTO 系统存在最佳传动比，在60 到80 优化区间内等差选取10 组增速比工况进行计算，其平均功率如表5 所示。据表可知，电机的最大平均功率出现在70 至75 增速比之间，故在波能装置设计中其最佳传动比应取70 至75 之间，并在此区间内配合选择齿轮箱的参数与齿轮带轮的外径。

表5 优化区间1号浮子平均功率Tab.5 No.1 buoy average power in optimized interval

5 结 论

（1）本文基于发电机和机械装置的工作原理，利用AMESim 搭建了PTO 负载模型，相比设置线性阻尼系数更加贴近真实情况。

（2）低传动比状态下，负载力矩对浮子的运动影响很小。随着传动比的增大，浮子的角速度将逐渐衰减。联合仿真模型可通过改变传动比来寻找PTO 系统输出功率的最佳参数区间，进而确定元件的最佳选型，在本文海况中，1 m 的圆柱形浮子机械装置的最佳传动比在72.5 附近，最大平均功率为35.1 W。

（3）三种做功方式中，浮子双行程做功的平均功率要大于单行程的平均功率。随着传动比的增大，其平均功率先增大再减小，传动比增大到80之后，PTO系统对浮子的负载力矩过大，浮子的纵摇角速度衰减严重，浮子无法充分吸收波浪的能量，导致双行程的平均功率降低为单行程的1.3倍。

（4）联合仿真的计算结果能指导装置最大平均发电功率的参数优化工作，本文的联合仿真模型可用于不同海况与不同浮子尺寸的振荡浮子式波浪能发电装置的设计、机械传动系统的参数优化和发电机的选型。

（5）联合仿真模型涉及水动力学与机械等学科领域，可用于指导实际工程中振荡浮子式波浪能发电装置的功率优化工作，确定最优功率下的发电机型号、齿轮箱增速比、浮子尺寸等数据，降低实际开发成本。在未来多浮子波浪能发电装置的开发中，联合仿真模型可以计算更复杂的PTO负载，从而简化多浮子波能装置中PTO 的机械设计与浮子的数量和布置等复杂的优化问题，缩短实际工程中波能装置的开发周期，为不同海况下的振荡浮子式波浪能发电装置的设计提供借鉴。