吴世其, 洪梅, 陈希, 毛科峰, 刘科峰
国防科技大学气象海洋学院, 江苏 南京 211101
波浪滑翔机是一种新型的无人海洋航行器, 最早由美国Liquid Robotics 公司研发。近年来我国也研制了“黑珍珠”、“海哨兵”、“海鳐”和“蓝鲸”等产品样机(Qi et al, 2013; 杨燕 等, 2014; 张森 等, 2014; 胡滕艳 等, 2018; 李灿, 2018; 桑宏强 等, 2018; 孙秀军 等, 2019)。波浪滑翔机的关键创新之处在于能够从海浪中获取大量能量, 为前进提供动力(Manley et al, 2010b)。同时, 位于母船上的太阳能电池板能够不断补充电量, 供给波浪滑翔机的控制系统、通信系统以及搭载的各类传感器(Mullison et al, 2011; Goebel et al, 2014; Carlon, 2015; O’Reilly et al, 2015; Penna et al, 2018)。
波浪滑翔机由水面母船通过绳缆与水下滑翔体连接而组成, 如图1 所示。推进系统是纯机械式的, 既不产生电力, 也不消耗电力(Hine et al, 2009; Manley et al, 2009)。如图2 所示, 波面上升时, 波浪滑翔机母船随波面上升带动水下滑翔体上升, 翼板偏转并受到向前及向下两个方向的分力; 同理, 波浪滑翔机随波面下沉时, 翼板偏转并受到向前及向上两个方向的分力。在波面的上下起伏中, 翼板带动波浪滑翔机整体向前运动(Manley et al, 2010a)。
由此可见, 水下翼板的各项参数在很大程度上影响着波浪滑翔机的航行效能。关于波浪滑翔机翼板参数设计, 国内外进行了较多研究。Kraus 等基于船舶六自由度模型, 建立了第一个适用于波浪滑翔机的运动模型, 通过模拟得到翼板最大偏转角度为20°最优(Kraus et al, 2011; Kraus, 2012)。Salari 等(2015)模拟了不同海况等级(根据海况等级表, 海况分为0~9 级, 海况越高所对应波高越高), 翼板最大偏转角度为10°或20°时, 波浪滑翔机翼板的运动状 况, 发现海况越高翼板所产生的推力越大。贾立娟(2014)对稳态情况下波浪滑翔机翼板的水动力特征进行了仿真, 并比较了翼板在不同间隔分布以及最大偏转角度下的升阻比, 进而确定翼板最优的参数设计。胡合文(2015)研究了翼板的数量、间距以及最大偏转角度对串列翼板平均推力系数的影响。严日华(2017)通过仿真试验分析了翼板数量、间距、最大偏转角度、转轴位置、运动幅值、运动周期以及来流速度对推进性能的影响。吕元博等(2018)分析了波浪滑翔机翼板的最大偏转角度、波高、波周期对翼板推力系数的影响。
图1 波浪滑翔机(Manley et al, 2010b) Fig. 1 The view of Wave Glider from below
图2 波浪滑翔机工作原理示意图(Manley et al, 2010a) a. 上升状态; b. 下降状态 Fig. 2 The operational principles of Wave Glider. (a) on the rise, (b) on the fall
这些研究都说明了翼板参数以及波浪要素对波浪滑翔机的推进效能影响较大, 并求得了在单一波浪要素条件下翼板的各项参数的最优设计。但是这些研究并没有指出针对不同的波浪特征是否存在不同的翼板参数设计。
针对这一问题, 本文首先基于再分析数据集, 对南海不同季节以及台风影响下的波高波周期特征进行统计。再利用FLUENT 软件, 研究不同波浪要素条件下, 翼板最大偏转角度、转轴位置以及翼板间距对波浪滑翔机水下翼板平均推力的影响。
本文选用的数据为欧洲中尺度天气预报中心大气和海洋全球再分析数据集 ERA5 (https://cds.climate. copernicus.eu/cdsapp#!/search?type=dataset) 中有效波高和平均波周期数据, 数据的空间分辨率为0.5°× 0.5°。选取1981—2010 年逐月时间分辨率数据, 对南海(4°—24°N, 104°—124°E)不同季节的波高和波周期进行统计。选取逐小时时间分辨率的数据, 查看2017 年第23 号台风“DAMREY”最强盛时南海波高波周期分布。该台风在2017 年11 月3 日6 时(UTC)最为强盛, 最大风速达到36m·s–1, 中心气压970hPa。
我们将每年3 至5 月的月平均数据按时间平均, 作为春季的波高和波周期分布情况。同理, 6 至8 月为夏季, 9 至11 月为秋季, 12 月至次年2 月作为冬季。所有季节时间上的平均为全年平均的数值, 结果如图3a~3e 和图4a~4e 所示。图3f 和图4f 为台风最强盛时刻的波高和波周期分布。
我们将南海海域空间内各季节及全年的波高、波周期的平均值作为对应时段内的代表值, 而台风状态下则是选取台风眼附近的波高、波周期数据作为代表值。具体数值如表1 所示。
翼板的运动分为两个部分, 一是由于波浪作用带动的升沉运动, 二是升沉运动导致的被动摆动运动。
图3 南海有效波高分布 a. 春季; b. 夏季; c. 秋季; d. 冬季; e. 全年; f. 台风最强盛时刻 Fig. 3 Significant wave height distribution in the South China Sea. (a) spring, (b) summer, (c) autumn, (d) winter, (e) annual, and (f) during typhoon
图4 南海波周期分布 a. 春季; b. 夏季; c. 秋季; d. 冬季; e. 全年; f. 台风最强盛时刻 Fig. 4 Wave period distribution in the South China Sea. (a) spring, (b) summer, (c) autumn, (d) winter, (e) annual, and (f) during typhoon
表1 南海海域波高和波周期数值 Tab. 1 Wave height and period in the South China Sea
水翼被动摆动的控制方程为:
式中: ω 为角速度(单位: rad·s–1); θ 为翼板当前偏转角度(单位: rad);maxθ 为翼板最大偏转角度; M 为水翼绕转轴力矩(单位: N·m), 通过对翼板表面压强积分得到; I 为翼板绕转轴的转动惯量(单位: kg·m2)。
本文采用二维瞬态不可压缩流进行求解。对于二维不可压缩流, 连续性方程和RANS 方程(雷诺平均方程)可以表示为:
本文采用SST k-ω 湍流模型(剪切应力输运k-ω湍流模型)计算涡动粘度。SST k-ω 湍流模型相比于标准的k-ω 湍流模型, 较好地避免了模型对自由流动的敏感性, 并能精确计算光滑表面的流动分离(Menter, 1994), 其方程式可以表示为:
选取NACA0012 型翼板, 弦长C 为0.2m。采用重叠网格, 背景网格区域为6m×5m 的长方形(台风浪研究则选取背景网格区域为6m×8m 的长方形)。选用四边形对该区域进行划分, 左边界设置为速度入口, 上下边界和右边界设置为压力出口, 翼板设置为无滑移壁面。单个翼板的重叠网格区域为直径0.6m 的圆, 串列翼板的重叠网格区域为6 个直径0.22m 的圆, 选用三角形对该区域进行划分。利用GAMBIT 软件进行网格划分, 翼板头部尾部适当加密。网格划分情况如图5 所示, 其中图5a 为单个翼板网格划分情况, 图5b 为串列翼板网格划分情况。
图5 单个翼板网格划分(a)和串列翼板网格划分(b) Fig. 5 Meshing single wing (a), and tandem wings (b)
为确保仿真计算的准确性和高效性, 计算了在4种网格尺寸下, 单个翼板在南海年平均的波浪要素条件下、最大偏转角为20°和转轴位置处于翼板前端1/3 处时, 翼板摆动产生的平均推力。表2 记录了4种网格尺寸下平均推力的变化, 本文最终采用A3 网格。串列翼板的网格可以按相同的密度进行划分。
表2 网格尺寸验证 Tab. 2 Mesh size validation
首先, 通过对单个翼板的仿真, 分别确定6 种波浪环境下最佳的偏转角度以及转轴位置。其次, 将计算得到的结果, 设置为串列翼板中每个翼板的参数, 改变翼板间隔, 比较翼板的平均推力。最终得到适用于南海的波浪滑翔机翼板的各项参数。
取转轴位置为翼板前端1/3 处, 改变翼板被动摆动的最大偏转角度(maxθ ), 分别对6 种不同波浪环境下翼板运动一个周期的平均推力进行比较。首先计算maxθ =15°、20°、25°、30°、35°、40°、45°时的平均推力, 再分别对极大值附近进行加密仿真计算。台风浪情形下,maxθ 在[15°, 45°]区间内, 翼板受到推力单调递增, 所以添加了maxθ =50°、55°、60°、65°、70°时的情况, 再对极大值附近进行加密仿真计算。
图6 显示在不同波浪要素条件下, 翼板获得的平均推力随最大偏转角度的变化情况。仿真结果表明, 对于同样的最大偏转角, 波高越大, 波浪滑翔机翼板系统能够获得的推力就越大。并且波高越大, 最大偏转角度也必须相应增大才能保证翼板获得更大的推力。针对南海全年平均的波浪要素, 最大偏转角度取27°为最佳, 而春季应为20°, 夏季应为22°, 秋季应为32°, 冬季应为38°, 台风影响下应为58°。
对NACA 型翼板的摆动研究大多选取翼板前端1/3 处(Anderson et al, 1998; Read et al, 2003), 本文再添加3/10、2/10 (1/5)、1/10 3 处转轴位置进行仿真。4 种转轴位置①、②、③、④如图7 所示, 编号①, 转轴位于翼板前端1/3 处; 编号②, 转轴位于翼 板前端3/10 处; 编号③, 转轴位于翼板前端1/5 处; 编号④, 转轴位于翼板前端1/10 处。
图7 翼板转轴位置示意图 ①为1/3 转轴位置; ②为3/10 转轴位置; ③为1/5 转轴位置; ④为1/10 转轴位置 Fig. 7 Schematic diagram of rotating axis
翼板最大偏转角度分别设置为27° (全年)、20° (春季)、22° (夏季)、32° (秋季)、38° (冬季)和58° (台风), 翼板转轴位置改变为上述4 处翼板转轴位置。仿真计算单个翼板运动一个周期的平均推力, 结果如图8 所示。
图8 显示了不同波浪要素条件下, 翼板获得的平均推力随转轴位置的变化情况。仿真结果显示, 在6 种波浪环境下, 翼板的转轴位置随推力变化的趋势基本一致, 并且获得推力的极大值点都发生在翼板前端1/5 处。
翼板最大偏转角度分别设置为27° (全年)、20° (春季)、22° (夏季)、32° (秋季)、38° (冬季)和58° (台风), 转轴位置取在翼板前端1/5 处, 翼板对数为6对时, 分别仿真翼板间距为40mm、80mm、120mm、160mm 4 种不同的情况, 结果如图9 所示。
图9 显示了在不同的波浪要素条件下, 串列翼板的平均推力随翼板间间隔的变化情况。图中可以看出, 翼板间距越大, 翼板获得的平均推力越大。但 由于翼板间间隔越大, 意味着水下滑翔体的质量越大。所以推力的增加, 并不能直接意味着推进效能的提高。
本文基于ERA5 再分析数据集, 对南海不同季节以及台风经过时的波高波周期特征进行统计, 得到代表值。再利用FLUENT 软件, 研究在这些波浪要素条件下, 翼板最大偏转角度、转轴位置以及翼 板间距对波浪滑翔机水下翼板的平均推力的影响。针对波浪滑翔机翼板的数值仿真试验得到主要仿真趋势如下:
图9 翼板平均推力随翼板间隔的变化情况 a. 四季和全年平均; b. 台风 Fig. 9 Average thrust obtained by the wing under different wave condition as a function of wings spacing. (a) four seasons and annual, and (b) during typhoon
1) 波高越大, 波浪滑翔机翼板系统获取的推力越大;
2) 波高越大, 最大偏转角度也必须相应增大才能保证翼板获得更大的推力。针对南海的波浪要素, 最大偏转角度在6 种波浪情况下分别取27° (全年)、20° (春季)、22° (夏季)、32° (秋季)、38° (冬季)、58° (台风影响下)为最佳;
3) 在南海波浪要素条件下, 翼板转轴位置选取在翼板前端1/5 处较佳;
4) 适当增加翼板间距能够提高翼板获得的推力。
本文中翼板转轴位置的改变过于理想化, 实际中还应考虑翼板的弹簧结构、最大承受压力等诸多因素。同时, 文中翼板运动的控制方程较为理想化, 实际中翼板的运动应结合波浪滑翔机整体动力学特征分析, 这些都是我们下一步工作的方向。