2. 上海交通大学航空航天学院, 上海 200240
2. School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240, China
0 引言
随着人类对于海洋探索的深入，深海对人身安全造成的威胁性也逐步加剧.尽管我国已经成功研制了“蛟龙号”深海载人下潜器，但是运行成本过高、风险大，应用面较窄.因此具备复杂环境和高危环境作业的水下机器人成为深海探测领域的利器[1-5].无人水下机器人(UUV)由于系统较为简单、成本低，得到了广泛应用.作为UUV的分支，自主水下机器人(AUV)发展迅速，能智能化地完成相应任务.民用方面，AUV可完成海底管道铺设、资源勘探和化学污染物处理等任务；军事国防方面，AUV在海底排雷、军事侦察和核电站反应器检查等中也有广泛应用.
仿生机器人是AUV的一类，它模仿了自然界中海洋生物特别是鱼类的生物性状.海洋中的鱼类经过长期的自然进化，具有推进效率高、前进阻力小、机动性强的特点.因此仿生机器人具有更好的适应性，是未来智能水下机器人的重要发展方向.鱼鳍运动方式分为两种，分别为身体/尾鳍推进(BCF)和中间鳍/对鳍推进(MPF)[6-7].前者通过身体的来回波动产生推力，例如鳗鱼和鲹科鱼类；后者通过鱼鳍的拍动产生推力，例如图 1所示的蝠鲼和鲀科类[8-10].应用到仿生机器人领域，BCF形式的运动方式会使机械结构较为复杂，设计和控制难度大；而MPF形式的运动方式较为简单，在实际中更易于实现. Chen等[11]采用离子聚合物复合材料研究了仿生蝠鲼机器人. Li等[12]受蝠鲼胸鳍拍动的启发，于2017年研究了可快速运动的电子鱼.王杨威等[13]提出了一种仿生环形长鳍波动推进器及其控制方法. Wang等人[14]借鉴蝠鲼外形设计了一种仿生蝠鲼水下滑翔机，并分析了浮力等特性随重心位置的变化规律. Luo等[15]使用CST参数化方法结合遗传算法对蝠鲼自动水下机器人的翼型进行了优化，优化结果使升阻比提高了10%.
|
| 图 1 蝠鲼科鱼类 Fig.1 Manta fish |
以上这些研究虽然初步揭示了MPF类仿生机器人的推进机制和可行性，但这些经验仍然不能满足精细化设计要求，需要更为精确的参数才能保证仿生机器人发挥最大使用价值.例如拍翼的上下拍动幅度、时间比例和最大俯仰角度都会显著影响翼型的升阻力特性，因此如何设定这些参数对机器人的行为具有重要影响.虽然目前已有关于空气扑翼的相关细致研究[16-17]，但空气与水的物理特性差异较大，因此不能完全借鉴空气扑翼的相关结论. Nonthipat等[18]针对伴有扭转弹簧的半主动式水下拍翼，利用势流理论模拟了翼型周围的流场，并分析了弹簧系数对推力和效率的影响.其他关于水下拍翼的精细研究尚未多见.为顺利地开展相关研制，需要对鱼鳍拍动过程中的迎角、拍动幅度、扭转变形等环节对升阻力的影响进行深入研究.
本文着重研究蝠鲼胸鳍的拍动过程，以二维的Clark Y翼型近似为鳍的剖面.采用动态混合网格对拍动周期内的上下-俯仰运动进行了数值模拟，得到了拍动幅度、俯仰幅度和上拍过程时间比对升阻力的影响规律.为进一步的三维拍动规律研究奠定了基础.
1 数值方法及验证本文采取动态混合网格技术对计算域进行离散，网格划分工具为Pointwise[19].翼型外围采用Hyperbolic法向增长的方式生成全四边形边界层网格，第一层网格高度为10-5c，c为翼型弦长，其余流场区域采用三角形网格填充.瞬态流场求解器为Fluent，湍流模型为SA一方程模型[20].该模型适合求解外流场问题，其计算开销低、稳定性好.通过求解关于修正湍流黏度的输运方程获得流场信息：
|
(1) |
式中：Gv为湍流黏度生成项，Yv为近壁面处产生的湍流黏度损耗项，
压力—速度项采取“Coupled”格式，空间和时间离散格式均为二阶精度.时间步长取为一个周期的1/400，单个时间步内迭代次数为50，确保计算收敛.
1.1 Delaunay动态网格技术Delaunay插值法是一种基于背景网格寻址的高效节点类动网格方法，针对中等变形问题具有很高的鲁棒性，它避免了求解大型矩阵的高耗时缺点[21].基本原理为：
1) 基于原始背景网格，建立起物面网格节点与远场边界节点的Delaunay三角形；
2) 在已建立的Delaunay三角形内定位流体域内的网格节点，过程如下：
任意一个网格节点与Delaunay三角形存在唯一的对应关系，这种关系可以由面积坐标来确定.如图 2所示，△ABC内存在一点P，则点P的位置可由面积坐标(e1，e2，e3)唯一确定：
|
(2) |
|
| 图 2 Delaunay三角形面积坐标 Fig.2 Delaunay triangle area coordinates |
其中，
|
(3) |
|
(4) |
|
(5) |
|
(6) |
显然对于面积坐标有：
|
(7) |
因此若当物面网格节点发生δ的位移后，流体域内网格节点位移由下式得到：
|
(8) |
3) 当物面网格发生变形时，Delaunay三角形的拓扑结构和流体域内网格节点的定位指针不变；
4) 根据定位指针更新得到变形后的流体网格.具体过程可参见文[19]，在此不作过多叙述.
1.2 网格无关性分析为使网格因素带来的误差最小，在如图 3所示的中等网格基础上，分别在物面节点数量和加密区网格上进行调整.最后建立了5套不同等级的计算网格，网格信息如表 1.来流速度U∞=10 m/s，攻角=0°，弦长c为0.3 m，介质为20 ℃时的液态水.在保证计算条件一致的情况下，得到这4种网格下物面的升阻力系数结果(如图 4)，综合考虑计算资源和计算精度，选取中等网格(对应3号)进行后续计算.
|
| 图 3 中等密度网格示意图 Fig.3 Medium density grid schematic |
| 网格编号 | 物面节点数 | 网格节点数 | 网格单元数 |
| 1 | 61 | 17 143 | 19 292 |
| 2 | 121 | 31 542 | 39 450 |
| 3 | 141 | 37 111 | 47 708 |
| 4 | 161 | 43 094 | 56 794 |
| 5 | 181 | 49 426 | 66 578 |
|
| 图 4 网格无关性分析 Fig.4 Grid independence analysis |
这里升阻力和扭矩系数分别定义如下：
|
(9) |
|
(10) |
|
(11) |
其中Fd与Fl分别为物面受到的阻力和升力，ρ为流体密度，S为拍翼等效面积.
1.3 数值方法验证Young等人通过实验测试了在不同拍动相位角下NACA0012翼型的动态推进性能[22]. Li等人[23]也对该测试结果进行了数值验证.其关于翼型重心的运动方式为
|
(12) |
|
(13) |
其中，y0和θ0分别为上下和俯仰运动的赋值，k为拍动频率控制参数，φ为俯仰相位角.式中参数取值分别为
|
当来流雷诺数取40 000时，由本文方法得到的平均推力系数与已有文献的对比结果如图 5所示，与Young和Li等的平均误差分别为-7.5%和5.4%，均在工程误差允许范围内，可见本文关于求解拍翼推进特性的数值方法是可信的.
|
| 图 5 数值方法可靠性验证 Fig.5 Reliability verification of numerical method |
真实的蝠鲼鳍拍动过程是极其复杂的，包含上下、俯仰和沿展向的波动.正是因为蝠鲼鳍的柔性特点，仿生结构只能以多段拍动逼近真实的柔性波动，从多段运动的角度来看，最基本的上下—俯仰运动仍然具备代表性.如图 6所示，一个完整的拍动周期T由两段构成，第一段(Part 1，下区)处于基准线以下，运动周期为T1，第二段(Part 2，上区)在基准线以上，运动周期为T2.在一个拍动周期内，上拍最高点到基准线的距离为pA，这里A代表波峰到波谷的距离，取为0.5 m，参数p称为上区比例.与此同时，翼型会经历一个完整周期的俯仰运动.
|
| 图 6 蝠鲼运动的数学简化 Fig.6 Mathematical simplification of manta movement |
因此翼型的上下运动可由如下数学表达所概括：
|
(14) |
俯仰运动通过翼型绕前缘点进行转动而实现，故俯仰运动过程可由式(15)表示：
|
(15) |
这里αmax为俯仰运动的最大拍动角.从流动方向来看，当α取负值时意味着翼型为抬头状态，反之为低头状态.
为保证拍动过程的光滑性，需增加速度约束条件，即保证两段拍动的速度光滑衔接：
|
(16) |
由此可得：
|
(17) |
则拍动频率为
|
(18) |
式(17)表示，Part 1与Part 2的运动周期之比只与p有关，p的取值可以反映出两段运动的相互关系. p < 0.5说明Part 1的周期高于Part 2，p>0.5说明Part 2的周期高于Part 1，当p=0.5时则为标准的正弦运动.现取拍动频率f=1 Hz，通过调整p的取值从而确定T1与T2.
3 结果与讨论为研究参数p和αmax对拍动过程升阻力的影响，现使p取遍0.1~0.9，增量为0.1；αmax取遍-25°~30°，增量为5°.为使图表简洁明了，在此只对部分p值下的升阻力结果进行分析.但需说明的是，本文所得结论与上述简化处理不相干涉.
3.1 参数p对升阻力的影响本节控制俯仰运动的最大拍动角αmax为0，即单纯研究扑动的上下运动而不考虑其俯仰过程，针对所有p值进行瞬态数值模拟，并取部分结果进行分析.
如图 7、图 8所示，横坐标为一个拍动周期内的相位，纵坐标为升阻力系数，负的阻力系数意味着拍动产生推力.对阻力结果而言：在所有情况中，拍动过程产生的最大阻力几乎一致，但最大推力却各不相同；且一个拍动周期内均出现两次波峰和两次波谷，两个波峰对应同样的阻力系数，出现在下拍最低点和上拍最高点，这里对应拍动速度转折点；随着p的增大，Part 1所经历的时间越短，此时波峰和波谷均发生左移，但波谷并不严格对应于Part 1和Part 2交接处，而是稍微滞后，这是因为对于瞬态问题而言，前一时刻的流场结构影响着后一时刻.对升力结果而言：整个拍动周期内只出现了一次波谷和一次波峰，且波谷出现在Part 1与Part 2交接处右侧，此时也是拍动速度发生转折的时刻.
|
| 图 7 不同p值下阻力系数随相位变化情况 Fig.7 Drag coefficient changes with phase at different p |
|
| 图 8 不同p值下升力系数随相位变化情况 Fig.8 Lift coefficient with phase change under different p |
结合图 7，图 8来看，当拍动产生最大推力时，同时也会形成最大升力，说明该时刻来流速度和拍动速度形成的耦合速度使得翼型处于最佳状态.由表 2可见最佳状态发生时刻tp随p近似线性递减，但可以看到的是，该时刻对应的平均等效迎角αe约为8.6°.
| p | tp/s | αe/(°) |
| 0.1 | 0.13 | 8.3 |
| 0.3 | 0.07 | 8.5 |
| 0.5 | 0.035 | 8.7 |
| 0.7 | 0.018 | 8.8 |
| 0.9 | 0.008 8 | 8.6 |
图 9和图 10分别是不同p值下，一个周期内拍动所产生的平均阻力和升力系数.可见不论p取何值，平均阻力系数始终取负值，说明拍动过程可以持续不断地产生推力.且平均阻力系数关于p对称分布，当p=0.5时，平均阻力系数取极小值，此时拍动过程产生的平均推力最大.而平均升力系数随着p值的增加而线性递减.相较于升力系数，我们更加关注水下机器人的推进性能，因此从平均阻力系数分布来看，当p=0.5时会得到最好的推进效果，此时的上下拍动过程为标准的正弦运动.
|
| 图 9 平均阻力系数随p的变化曲线 Fig.9 Average drag coefficient as a function of p |
|
| 图 10 平均升力系数随p的变化曲线 Fig.10 Average lift coefficient as a function of p |
结合上节结论，本节取上区比例p为0.5，令αmax取遍-20°~30°，分析最大拍动角对拍动过程升阻力的影响.
如图 11所示，显示了一个拍动周期内阻力系数随相位的变化情况.当拍动角的绝对值越大(从-25°~0°)，表明拍动过程中的最大迎风面积越大，因此会受到更大的阻力.阻力对角度的敏感程度很高，其变化趋势几乎与拍动角的变化趋势保持一致，即当拍动角绝对值增大使迎风面积增大则阻力系数立即上升.随着拍动角绝对值的降低，最大阻力系数逐渐减小.当αmax=-5°时，阻力系数变化趋势开始反向，原来的波峰相位变为波谷相位，一方面因为此时迎风面积的减小使受到的阻力大大降低，另一方面是因为此时的拍鳍处于抬头状态，在来流速度和拍动速度的耦合下，使得等效迎角处于有利减阻的状态.这种有利减阻的状态在αmax=5°时更为突出，此时拍鳍处于下拍的低头状态，此刻阻力转变为推力. 图 12为一个拍动周期内升力系数的历程曲线.可见其分布与阻力系数类似，当拍动角绝对值越大时，在抬头状态下的等效迎角越大，因此其升力系数波峰越高.当αmax超过-5°时，波峰转变为波谷，原因在于此时等效迎角为负，拍鳍产生负升力.
|
| 图 11 不同拍动角下的阻力系数历程曲线 Fig.11 Drag coefficient history curve at differentαmax |
|
| 图 12 不同拍动角下的升力系数历程曲线 Fig.12 Lift coefficient history curve at differentαmax |
如图 13所示，为平均阻力系数随αmax的变化情况，可见当αmax大于10°或αmax小与-5°时，其综合阻力系数为正，产生阻力；而当αmax大于-5°小与10°时，其综合阻力系数为负，产生推力.因此在实际设计中，不可将最大拍动角的大小设置过大，以免造成阻力.从本文结果来看αmax取5°可使平均推力最大.另外，从图 14可见，αmax取5°也可使拍鳍处于良好的升力状态.
|
| 图 13 平均阻力系数随最大拍动角的变化 Fig.13 The average drag coefficient varies withαmax |
|
| 图 14 平均升力系数随最大拍动角的变化 Fig.14 The average lift coefficient varies withαmax |
由3.1小节和3.2小节可以确定最佳的拍动参数组合为p=0.5，αmax=5°.本节针对该组合的拍动过程进行详细分析.
如图 15所示，为上述参数组合下的拍翼在一个周期内的运动过程，云图显示的是拍动中的Q涡量.可见在下拍阶段(相位=0~0.5π)附体涡主要集中在拍翼头部和拍翼上部，拍翼尾部也带起一小部分的尾涡，且随着下拍进行，拍翼低头，附着在拍翼上表面的涡逐渐减少.当运动到最低点时(相位=0.5π)，拍翼呈水平状，此时拍翼上表面的涡消失.随着拍翼的上拍(相位=0.75π~1.5π)，拍翼先低头再抬头，在这个过程中涡结构开始在拍翼下表面生成，且随着上拍的进行，下表面的涡结构渐次减少直至在上拍最高点(相位=1.5π)消失.此后拍翼下降(相位=1.5π~2π)，继续抬头，涡结构又在拍翼上表面聚集直至完成一个周期的拍动.
|
| 图 15 最佳参数组合下的拍动过程 Fig.15 Flapping process under optimal parameters combination |
在此引入拍翼的拍动效率[23]：
|
(19) |
式中，
|
(20) |
|
(21) |
其中，Cm为拍动过程中的力矩系数.
故在最佳组合情况下的拍动效率为
|
(22) |
本文以Clark Y翼型近似为二维蝠鲼拍鳍，研究了上区比例p和最大拍动角αmax两个重要参数对拍动过程中升阻力的影响规律，并得到了最佳的参数组合.作为结论，以下几点需被重点提及：
1) 从平均阻力系数分布来看，取p=0.5会得到最好的推进效果，此时的上下拍动过程为标准的正弦运动.
2) 平均等效迎角约为8.6°时，拍动产生最大推力的同时也会形成最大升力.
3) 在实际设计中，不可将最大拍动角的大小设置过大，以免造成阻力激增. αmax取5°可使平均推力最大，也可使拍鳍处于良好的升力状态.
4) 文中仅对二维情况进行了分析，确定了两个主要参数的组合，对拍鳍设计提供了初步的依据.而真实的拍鳍运动会更加复杂，具有更高的自由度，例如拍鳍沿展向的柔性波动等，这是下一步更加精细化研究的重点.
| [1] | Nonthipat T, Surasak P, Varangrat J. Semi-active flapping foil for marine propulsion[J]. Oceaning Engineering, 2018, 147: 556–564. DOI:10.1016/j.oceaneng.2017.11.008 |
| [2] | Alexander P. Robot Fish:Bio-inspired fishlike underwater robots[J]. Underwater Technology, 2017, 34(3): 143–145. DOI:10.3723/ut.34.143 |
| [3] |
周浩, 姜述强, 黄海, 等.
基于视觉感知的海生物吸纳式水下机器人目标捕获控制[J]. 机器人, 2019, 41(2): 242–249.
Zhou H, Jiang S Q, Huang H, et al. Vision based target capture control for sea organism absorptive underwater vehicle[J]. Robot, 2019, 41(2): 242–249. |
| [4] |
熊俊峰, 何玉庆, 韩建达.
水面机器人主动增强LPV航向保持控制[J]. 信息与控制, 2018, 47(3): 267–275.
Xiong J F, He Y Q, Han J D. Active LPV yaw keeping control of an unmanned surface vehicle[J]. Information and Control, 2018, 47(3): 267–275. |
| [5] |
王超, 胡志强, 衣瑞文, 等.
高速水下机器人通气空化减阻技术的水洞实验研究[J]. 机器人, 2018, 40(6): 779–785.
Wang C, Hu Z Q, Yi R W, et al. Water tunnel experiment research of ventilated cavitation drag reduction technology for a high speed UUV[J]. Robot, 2018, 40(6): 779–785. |
| [6] | Launder G V. Fish locomotion:Recent advances and new directions[J]. Annual Review of Marine Science, 2015, 7: 521–545. DOI:10.1146/annurev-marine-010814-015614 |
| [7] |
王延杰, 郝牧宇, 张霖, 等.
基于智能驱动材料的水下仿生机器人发展综述[J]. 水下无人系统学报, 2019, 27(2): 123–133.
Wang Y J, Hao M Y, Zhang L, et al. Progress of biomimetic underwater robot based on intelligent actuating materials:A review[J]. Journal of Unmanned Undersea Systems, 2019, 27(2): 123–133. |
| [8] |
云忠, 温猛, 蒋毅, 等.
仿生蝠鲼胸鳍摆动推进机构设计与水动力分析[J]. 浙江大学学报(工学版), 2019, 53(5): 872–879.
Yun Z, Wen M, Jiang Y, et al. Design and hydrodynamic analysis of pectoral fin oscillation propulsion mechanism of bionic manta ray[J]. Journal of Zhejiang University(Engineering Science), 2019, 53(5): 872–879. |
| [9] | Chew C M, Lim Q Y, Yeo K S. Development of propulsion mechanism for robot manta ray[C]//IEEE International Conference on Robotics and Biomimetics. Piscataway, NJ, USA: IEEE, 2016: 1918-1923. |
| [10] | Cai Y, Bi S, Zheng L. Design optimization of a bionic fish with multi-joint fin rays[J]. Advanced Robotics, 2012, 26(1/2): 177–196. |
| [11] | Chen Z, Um T I, Bart S H. A novel fabrication of ionic polymer-metal composite membrane actuator capable of 3-dimensional kinematic motions[J]. Sensors and Actuators A:Physical, 2011, 168(1): 131–139. DOI:10.1016/j.sna.2011.02.034 |
| [12] | Li T, Li G, Liang Y, et al. Fast-moving soft electronic fish[J]. Science Advances, 2017, 3(4): e1602045. DOI:10.1126/sciadv.1602045 |
| [13] |
王扬威, 闫勇程, 刘凯, 等.
基于CPG的仿生环形长鳍波动推进器运动控制[J]. 机器人, 2016, 38(6): 746–753.
Wang Y W, Yan Y C, Liu K, et al. Motion control of a bionic circular long-fin undulating propeller based on CPG[J]. Robot, 2016, 38(6): 746–753. |
| [14] | Wang Z, Yu J, Zhang A. Hydrodynamic performance analysis of a biomimetic manta ray underwater glider[C]//IEEE International Conference on Robotics and Biomimetics. Piscataway, NJ, USA: IEEE, 2016: 1631-1636. |
| [15] | Luo Y, Pan G, Huang G Q, et al. Parametric geometric model and shape optimization of airfoils of a biomimetic manta ray underwater vehicle[J]. Journal of Shanghai Jiao Tong University:Science, 2019, 24(3): 402–408. DOI:10.1007/s12204-019-2076-4 |
| [16] |
常兴华, 马戎, 张来平, 等.
S1223翼型俯仰-沉浮运动的非定常气动特性分析[J]. 空气动力学学报, 2017, 35(1): 62–70.
Chang X H, Ma R, Zhang L P, et al. Numerical study of the phunging-pitching motion of S1223 airfoil[J]. Acta Aerodynamica Sinica, 2017, 35(1): 62–70. |
| [17] |
肖天航, 罗东明, 郑祥明, 等.
仿生飞行器非定常气动优化设计研究进展与挑战[J]. 空气动力学学报, 2018, 36(1): 81–87.
Xiao T H, Luo D M, Zheng X M, et al. Progress and challenges in unsteady aerodynamic optimization design of bionic flapping-wings[J]. Acta Aerodynamica Sinica, 2018, 36(1): 81–87. |
| [18] | Nonthipat T, Surasak P, Varangrat J. Semi-active flapping foil for marine propulsion[J]. Ocean Engineering, 2018, 147: 556–564. DOI:10.1016/j.oceaneng.2017.11.008 |
| [19] | Mesh generation software for CFD-Pointwise, Inc. |
| [20] | Beaugendre H, Morency F. Penalization of the Spalart-Allmaras turbulence model without and with a wall function:Methodology for a vortex in cell scheme[J]. Computers and Fluids, 2018, 170: 313–323. DOI:10.1016/j.compfluid.2018.05.012 |
| [21] | Xue Q L, Ning Q, Hao X. Fast dynamic grid deformation based on Delaunay graph mapping[J]. Journal of Computational Physics, 2006, 211: 405–423. DOI:10.1016/j.jcp.2005.05.025 |
| [22] | Young J. Numerical simulation of the unsteady aerodynamics of flapping airfoils[D]. Sydeny, Australian: University of New South Wales, 2005. |
| [23] | Xin T L, Yi L L, Jia Q K, et al. Reduced-order thrust modeling for an efficiently flapping airfoil using system identification method[J]. Journal of Fluids and Structures, 2017, 69: 137–153. DOI:10.1016/j.jfluidstructs.2016.12.005 |