奚泉，唐渔滢，于通顺，于春明，曾兴井，徐琨
1.中国电建集团昆明勘测设计研究院有限公司，云南 昆明 650051
2.青岛市交通工程质量安全监督站，山东 青岛 266101
3.中国海洋大学 工程学院，山东 青岛 266400
4.中国葛洲坝集团市政工程有限公司，湖北 宜昌 443000
在“双碳”的战略目标下，海洋可再生能源迎来发展新格局[1]。波浪能因具有开发难度较低、获取方式简单等优势，被认为是品质最高的海洋能[2-3]。振荡浮子波能发电装置[4-5]因结构简单、能量转换效率高等优势是目前常见的三大波浪能利用技术之一。
振荡浮子装置大多是由波能俘获装置、质量体、系泊系统以及发电装置组成，其工作原理是利用俘获装置在波浪的作用下振荡起伏来达到汲取波浪能的目的，之后再由能量摄取（power takeoff，PTO）系统实现波浪能向电能转化。振荡浮子式波浪能发电装置的水动力性能[6]、能量转换效率[7]及性能优化[8]成为目前研究的热点。美国Ocean Power Technologies (OPT)公司研发了PowerBuoy 振荡浮子波能装置，并对其进行了设计优化[9]；
Shi[10]和Han等[11]对振荡浮子获能特性进行了理论分析，并进一步提出了获能谱的概念；
于通顺等[12-13]探究了振荡浮子前侧波浪爬升特性，为后续浮子竖向尺寸设计提供理论依据；
Agamloh等[14]利用有限体积法分别模拟单个浮子和双浮子在不同规则波作用下的运动过程，探究了浮子的能量输出功率。吴必军等[15]和余海涛等[16]均对单浮子波能装置的俘获宽度比进行了试验研究，并指出波浪频率是影响装置俘获宽度比的主要因素。Waters等[17]对驱动线性发电机的浮子进行了试验研究，提出非线性外负载是影响装置俘获宽度比的主要因素。
通过以上研究发现，前人的研究中尚未开展浮子垂荡过程中表面所受的冲击作用的研究。振荡浮子在复杂海况中通过大幅垂荡获能，垂荡过程中浮子与波浪作用过程剧烈，忽略了浮子表面受到的强烈波浪冲击作用，将会增大浮子长期运行的风险。因此本文基于FLOW 3D 建立波浪与浮子耦合作用的数值模型，开展浮子所受冲击作用过程的研究，探究浮子所受冲击荷载的作用规律。
1.1 控制方程
对于波浪和结构物相互作用的问题，波浪被认为是不可压缩黏性流体的运动。考虑湍流流动，将连续性方程和雷诺平均纳维-斯托克斯（Reynolds-averaged Navier-Stokes，RANS）方程作为流体运动的控制方程。假设右手笛卡尔坐标系O-xyz以未受扰动自由面的原点定义，则连续性方程可以写作：
动量方程为
式中：ρ为流体密度，P为流体内部压强，VF为可流动的体积分数，u、v、w分别为在x、y、z这3 个方向上的流速分量，Ax、Ay、Az分别为x、y、z这3 个方向可流动的面积分数，Gx、Gy、Gz分别为对应x、y、z这3个方向的重力加速度，fx、fy、fz分别对应3个方向的黏滞力加速度。
1.2 紊流模型
本文采用重组化RNGk-ε模型，其表达式为
式中：PT为由剪切效应引起的湍流动能的生成项，KT,Diff和 εT,Diff为扩散项，CDIS,1和CDIS,2为模型系数，k为由流速梯度引起的紊动动能的产生项。
1.3 运动与碰撞模型
Flow-3D 中运动与碰撞模型可以模拟流固耦合下结构物的运动，该运动可以是指定的运动，也可以是与流体动态耦合的运动，移动物体都可以具有6 个自由度。该模型还允许存在多个（最多100 个）线性弹簧和扭力弹簧、弹性绳索以及系泊缆绳，它们可以附接到移动的物体上，并可对其施加力或扭矩。根据刚体运动学，刚体的一般运动可分为平移和旋转运动。结构物上任何点的速度等于结构物上任意选择的基点速度加上绕该基点旋转引起的速度。选择S为刚体上的基点，R为刚体上任意点，则R点速度表达式为
式中：VS为基点的速度；
rR/S为点R相对于基点S的矢量距离；
ω是刚体的角速度矢量，与基点的选择无关。
1.4 自由表面追踪的流体体积方法
FLOW-3D 中采用了流体体积方法进行自由追踪。定义流体计算域内独个单元的流体体积函数为F（x,y,z,t)。F为流体所占单元体的比值，若单元全部充满流体，则F值为1；
若不存在流体时，F值则为0；
而当F介于0～1 时，则表明该网格内存在自由表面。F函数表示单位体积，并满足方程：
2.1 三维数值波浪水槽的建立
建立数值波浪水槽的侧面示意如图1（a）所示，水槽长度取6 倍左右波长，宽度取浮子宽度的5 倍。为便于后续进行数值模型的验证，选用Yu等[12]试验中浮子的尺寸进行模型的建立：圆柱型浮子直径为0.4 m，高度为0.4 m，吃水深度为0.12 m。为了合理描述浮子的形状，该水槽由3 个网格块组成，中间网格块为浮子运动加密区，浮子放置在距离入口的2 倍波长(L)范围内。水槽的边界条件定义为：最左端是可以实现函数造波功能的波浪，最右端为出流，水槽的底部设为壁面，上部为压力，水槽顺波向两侧对称。如图1（b）和图1（c）所示。水槽尾端部分（约1/4 水槽长度）设有阻尼消波区用来抑制波浪反射，左端和右端的阻尼系数分别设为1 和8。
图1 三维数值波浪水槽边界设置及网格块划分
2.2 网格的划分及网格独立性验证
本模型采用波高范围内3 种不同的网格数量进行了波面网格独立性试验，即设置了网格数量分别为10、20 和30 的3 个不同算例，流域内网格的大小和数量如表1 所示。
表1 无浮子时计算流域网格数量和尺寸
图2 比较了3 种网格尺寸下计算区域内相同位置（2L处）自由面高程的历时曲线，结果发现，加密网格数为20 与30 的计算结果差距较小，加密网格数为10 的计算结果与其他2 种网格数量的计算结果差距较大，但是加密网格数为30 的运行时间明显增加。综合考虑模型计算时间和精度，本文在Z方向的波高(H)范围内网格尺寸建议设置为 Δz≤H/20。
图2 不同网格尺寸的自由面高程历时曲线
为了合理描述浮子的形状，在浮子运动加密区采用3 种不同总网格数量，浮子周围网格数量如 表2所示，表中Δx、Δy和Δz分别为沿x、y和z方向的网格尺寸。所划分网格可以较好地描述浮子的形状。
表2 浮子周围的网格数量和尺寸
图3 比较了在不同网格尺寸下采用上述波浪条件捕捉浮子垂荡运动的能力。结果表明，0.01 cm和0.005 cm 的网格尺寸下浮子位移相差不大，说明以0.01 cm 划分的网格具有足够高的模拟分辨率。因此，浮子周围网格尺寸建议设置为Δx 图3 不同网格尺寸的浮子垂荡位移历时曲线 2.3.1 造波能力验证 图4 为距波浪入口2L附近（放置浮子处）规则波的数值计算结果与理论值的时程曲线（波高H=0.2 m，周期T=2 s，水深h=1 m）。在入射波传递一段时间到达稳定状态后，很明显波峰处略尖陡，波谷处略平坦，且相对误差低于4%； 图4 自由波面理论值与计算值的时程曲线 2.3.2 波浪对浮子作用的验证 为验证数值模型对浮子运动捕捉的准确性，选用Yu等[12]的试验数据与数值计算结果进行比较。表3 给出了连续5 个振动周期中浮子运动最高点与最低点的数值结果与试验数据的对比（波高H=0.2 m，周期T=2 s，水深h=1.0 m）。可以看出，数值结果与试验数据的相对误差均在10%以内，所建立数值模型能够较好地模拟波浪对浮子的作用。 表3 垂荡位移数值结果与试验数据 2.3.3 波浪冲击压力观测点设置和定义 根据对称性，在浮子迎浪方向、侧面以及后方（分别为Row A、Row B 和Row C）设置3 列压力传感器，定义s为浪高仪至浮子底端的垂直距离，d为浮子高度，浮子柱面观测点和压力传感器的布置如图5（a）和表4 所示。浮子底部设置了11 个压力传感器，θ为所在点与x轴负方向的夹角，相邻2 组测点之间的夹角为45°，浮子底面观测点和压力传感器的布置如图5（b）和表5 所示。模型中所有虚拟观测点的布置与实际位置误差小于3 mm。 表5 浮子底面压力传感器设置 图5 波浪冲击压力虚拟观测点设置 表4 浮子柱面压力传感器设置 波浪在冲击的过程中涉及到诸多因素，比如水气掺混、瞬时效应等，这些因素导致波浪冲击垂荡浮子的过程具有很强的随机性，因此冲击过程中的压力峰值具有很强的离散性。本文根据文献[18]对每个虚拟的压力传感器采集的压力数据进行了处理，即浮子上的压力测点取样时，将每个周期内的压力峰值由大到小排列，并且取前1/3 个峰值的平均值作为各测点的冲击压力峰值特征值，即统计分析的特征值。每个测点的冲击压力峰值的实测数据记为Pt，浮子静止时测点的初始压力记为P0，将实测数据与初始压力无量纲化后记为Pt/P0。浮子所受的水平和竖直方向波浪荷载峰值数据的提取则是将每个周期内的压力峰值提取出来，并且对峰值相差较小的几组数据取平均值为波浪荷载的特征值。水平和垂直波浪荷载分别用Fx和Fz来表示，对波浪荷载无量纲化分析，选取浮子的重力ρgSd为无量纲化参数，其中ρ为浮子的密度，g为重力加速度，S为浮子底面积。 2.3.4 计算工况 数值计算中选用圆柱形浮子，实际浮子高度为1.6 m、质量为3 860 kg、吃水深度为0.48 m。计算工况如表6 所示。 表6 计算工况设置 3.1.1 冲击压力峰值在浮子的立侧处分布规律 图6 为工况1 中A3 传感器记录的浮子波浪压力历时曲线，由图可以发现，浮子受到的冲击压力具有明显的周期性，周期与入射波周期保持一致，冲击压力近似为正弦脉冲形式，冲击阶段压力逐渐增大，之后随着波面下降开始逐渐回落，当测点脱离水体，压力维持在零值几乎不变。 图6 浮子立侧不同压力测点冲击压力历时曲线 图7 为工况1 条件下浮子表面不同测点上的压力分布，1#和2#测点冲击压力在不同位置处差距较小，3#测点的冲击压力差别较大。因为波面高度达不到，故4#和5#测点的冲击压力为0。距离浮子底部越大，无量纲冲击压力值越大，即相较于初始压力提升越大，3#测点的无量纲冲击压力达到最大值。距离浮子底部相同距离的测点中，Row C 侧的冲击压力最大，Row A 侧的冲击压力最小，表明浮子后侧所受到的冲击压力最大。 图7 浮子立侧处不同测点压力分布 3.1.2 冲击压力峰值在浮子底面分布规律 图8 为工况1 条件下M1#、O#和S1#测点记录的浮子的底部波浪压力历时曲线，由图可以发现，浮子底端受到的冲击压力具有明显的周期性，周期与入射波周期保持一致； 图8 浮子底面不同压力测点冲击压力历时曲线 3.2.1 波陡参数对浮子周围波浪冲击效应影响 图9 给出了不同波陡参数下，最大波浪冲击压力沿的垂直分布。 从图9 中可以发现，最大波浪冲击压力随着波陡参数的增加而增大，距离浮子底部距离越近的位置上波浪冲击压力越大。波陡较大时，浮子靠上位置处的波浪冲击压力最小值出现在Row B 侧，这是由于浮子侧面的波浪爬升高度处于最低导致，其他位置的最大冲击压力出现在浮子中间（Row C），最小冲击压力依然发生在迎浪侧（Row A）。 图9 不同波陡参数下浮子上最大压力垂直分布 3.2.2 散射参数对浮子周围波浪冲击效应影响 图10 给出了不同散射参数下，最大波浪冲击压力沿Row A、Row B 和Row C 的垂直分布。从图10 可以看出，改变浮子半径对于Row B 和Row C侧相同压力测点上波浪冲击压力峰值的改变较小。散射参数较小时，Row A 侧的相同测点压力峰值随散射参数变化较小； 图10 不同散射参数下浮子上最大压力垂直分布 3.2.3 线性PTO 阻尼参数对浮子周围波浪冲击效应影响 图11 给出了不同弹簧刚度系数下，最大波浪冲击压力沿Row A、Row B 和Row C 的垂直分布。 从图11 中可以看出，弹簧刚度系数的增加对于浮子上的波浪冲击压力影响较小，浮子上相同压力测点上的压力峰值随弹簧系数的增加几乎保持不变。通过对比发现，施加PTO 后对于浮子上的波浪冲击压力改变很小，且浮子上的波浪冲击压力垂直分布规律与浮子自由垂荡运动时保持一致。 图11 不同弹簧刚度系数下浮子上最大压力垂直分布 图12 给出了不同阻尼系数下，最大波浪冲击压力沿Row A、Row B 和Row C 的垂直分布。从图12 中可以看出，增加阻尼系数对浮子截面上波浪冲击压力垂直分布的影响较小； 图12 不同阻尼系数下浮子上最大压力垂直分布 浮子侧面的波浪冲击压力受波陡参数的影响最大，冲击压力随波陡增大而增大； 图13 和图14 给出了工况3 水平和垂直波浪荷载与Row A 压力测点时程对照。A1 和A2 处的波浪压力与水平波浪荷载的历时曲线存在较小的相位差，水平波浪荷载的最大值与最小值出现的时间较早于A1 处的波浪压力最值； 图13 水平波浪荷载与Row A 压力测点时程曲线对照 图14 垂直波浪荷载与Row A 压力测点时程曲线对照 由于水平波浪荷载受力位置与Row A 侧相似，都为迎浪面，因此历时曲线存在的相位差较小，在浮子较高位置处的相位差几乎一致。但浮子所受垂直波浪荷载的受力面与浮子Row A 侧波浪压力受力面存在较大的差异，从图14 中可以发现，垂直波浪荷载的最值出现时间和该位置处的波浪冲击压力最值出现时间存在较大差异。A1 和A2 处的波浪压力与垂直波浪荷载最值出现的时间相反，即垂直波浪荷载出现最大值时，该位置的波浪压力处于最小位置； 1）在规则波作用下，浮子所受冲击压力具有明显的周期性，浮子底部波浪冲击压力近似为正弦脉冲形式。 2）在入射波前进方向上，随着浮子底部传感器与浮子前部距离的增大，传感器所受的冲击压力峰值也增大。 3）波浪冲击压力受波陡参数的影响最大，冲击压力随波陡增大而增大。浮子迎浪侧的波浪冲击压力随散射参数的增大有增加的趋势，其他位置处的冲击压力受散射参数的影响不明显。浮子上的最大波浪冲击压力受PTO 阻尼系统的影响很小。 4）浮子迎浪侧距离浮子底部较高位置处的波浪冲击压力与波浪荷载最值出现时间几乎保持一致。 5）本研究只针对线性波浪作用下的浮子冲击压力进行了研究，随机波浪作用下浮子荷载特征的研究需进一步开展。
2.3 三维数值水槽验证
同时可以发现波浪的反射效应并不明显，这表明水槽的造波能力和消波效果都比较可靠。
3.1 波浪冲击压力特征
浮子底部中心所受的波浪压力小于波浪传播方向上中心两侧的波浪压力。
3.2 波浪参数和阻尼参数对浮子上波浪冲击压力垂直分布的影响
散射参数较大时，Row A侧的相同测点压力峰值有明显提升。总的来说，增加浮子半径对于浮子截面上相同位置的波浪冲击压力影响较小。
增加PTO 阻尼对浮子垂荡运动有影响，但浮子上的波浪冲击由于波浪参数保持一致，导致浮子上的波浪冲击压力受其他因素的影响很小。
Row A 侧的波浪冲击压力随散射参数增大而增加，其他位置处的冲击压力受散射参数的影响不明显；
浮子上的最大波浪冲击压力受PTO 系统阻尼的影响很小；
浮子侧面的压力分布在s>0.24 m 位置上的压力处于最小值，这主要是由于波浪在该位置的爬升较小，冲击强度较低所导致。3.3 最大波浪荷载与波浪冲击压力关系研究
浮子迎浪侧距离底部较高位置处（A3、A4 和A5）的波浪压力与水平波浪荷载的历时曲线几乎不存在相位差，即水平波浪荷载的最大值与最小值出现的时间和该位置处的波浪冲击压力最值同步出现。
浮子迎浪侧距离底部较高位置处（A3、A4 和A5）的波浪压力与垂直波浪荷载的最值同步出现，与水平荷载相似。