生命诞生于海洋，经历了数十亿年的物种选择与进化，造就了数以万计的流体力学大师—海洋生物。鱼类作为常见的水生生物，遍布于江河湖海，大多数鱼类通过摆动尾鳍产生推力，在水中机动灵活，同时拥有极高的推进效率^[1-2]。人们通过观察鱼类的生物构造与运动特性，总结提取物理原理及运动规律，最终服务于仿生推进器设计，摆翼推进器就是一种模仿鱼类波状游动的新型推进装置，相比传统推进装置，摆翼推进器拥有极高的航行效率，较低的噪声污染、较强的机动性能等优势，近年来成为越来越多科研工作者关注的热点问题。

Zhu^[3]在研究中指出，柔性鱼尾是水生动物提高运动效率的重要机制，柔性水翼的流固耦合作用能进一步提高推进效率。Dylan Iverson等^[4]研究指出柔性水翼产生了比刚性水翼更高的推力值，同时也具有更高的推进效率，但柔性的水翼存在大量的控制参数，不利于变量控制。Anderson等^[5]研究发现，刚性水翼在一定条件下的推进效率仍可达到惊人的87%，远高于目前主流的螺旋桨推进器。摆翼推进器的研究可促进推进效率大幅提高，节约能源，减少碳排放，具有很高的工程、环保意义。

从研究现状来看，基于二维水翼求解，忽略水翼的三维效应，化繁为简的研究已经取得了大量的理论和实验成果。然而，真实的三维流动不同于二维的情况，由于绕流在上下翼面产生的压强差，三维水翼表面的流动存在展向速度，从而形成尾涡层，在翼尖处向下游卷成两个旋向相反的喇叭状旋涡^[6]，产生诱导阻力耗散能量。特别是小展弦比的水翼，翼尖涡与前缘涡相互作用强烈，对流动产生的影响不容忽视。因此，本文考虑了三维效应，对三维水翼进行了数值仿真，研究了多种形状的水翼水动力性能变化规律，为摆翼推进器的设计提供参考。

本文探究了水翼的展弦比、根梢比、前掠和后掠对水动力性能的影响。定义水翼展长为$ H $，翼根弦长为$ c $，翼尖弦长为$ {c}' $，展弦比$ AR=H/c $，根梢比$ TR=c/{c}' $。由于水翼左右对称，各算例均使用半模计算，翼根弦长取$ c=0.2 $m。

模型投影形状如图1所示：

对以上形状的水翼进行建模，控制不同的几何变量，生成了23种水翼模型。其中，展弦比不同的矩形翼有8种，梯形翼、前掠及后掠翼分别有5种根梢比。考虑到三维效应可能是梯形翼、前掠及后掠翼性能的重要的影响因素之一，选择比较适中且能体现三维效应的展弦比为4的水翼，作为初步探讨。每种模型将在进速$ {V}_{A} $在0~8 m/s之间的13种工况下进行数值仿真，模型详细信息如表1所示：

图  1  水翼投影形状示意图

Figure 1.  Schematic diagram of foil projection

参照李明阳^[7]的工作，水翼设置为前进（X方向平移）、横荡（Y方向振荡）、艏摇（绕Z轴旋转）三自由度耦合运动，水翼运动使用绝对坐标系$ xOy $如图2所示。水翼匀速前进，横荡、艏摇运动均为正弦规律，横荡运动方程为$ Y={Y}_{0}{sin}\left(2\pi ft\right) $，式中$ Y $为横荡位置，$ {Y}_{0} $为最大横荡摆幅，$ f $为摆动频率，艏摇运动方程为$ \theta ={\theta }_{0}{sin}\left(2\pi ft+\phi \right) $，式中$ \theta $为艏摇摆角，$ {\theta }_{0} $为最大艏摇摆角，$ \phi $为横荡运动和艏摇运动的相位差，旋转运动参考点为镜像面内翼型导缘顶点。

图  2  水翼运动、受力情况示意图

Figure 2.  Schematic diagram of foil motion and force

本文取横荡摆幅$ {Y}_{0}=0.5\;\mathrm{m} $，摆幅比$ {Y}_{0}/c=2.5 $，运动频率$ f=1\;\mathrm{H}\mathrm{z} $，最大艏摇摆角$ {\theta }_{0}=0.3\;\mathrm{r}\mathrm{a}\mathrm{d} $，横荡与艏摇运动相位差为$\varphi =\pi /2$，一周期内水翼运动、受力情况如图2所示：

Floc"H^[8]类比了螺旋桨，定义了摆翼推进器的进速系数，其选用的特征长度为等效直径$ {D}_{eq}=\sqrt{4BH/\pi } $，式中B为横荡扫掠宽度，H为水翼展长。本文使用了类似的方式，选用B作为特征长度，定义进程$ {h}_{P}=T{V}_{x} $与$ B $的比值为进速系数，以$ J $表示：

式中，$ T $为摆动周期，$ {V}_{x} $为前进速度。上式可知，在摆动频率f、扫掠宽度$ B $一定的条件下，改变$ {V}_{x} $值即可得到不同的$ J $，本文水翼运动的不同工况即指不同进速系数的工况。

许多文献中，采用St数代替进速系数来描述摆翼的不同工况，两者关系表示为：

定义绝对坐标系下，一个周期内沿前进方向力的积分均值为平均推力$ \bar {{F}_{x}} $，表示为：

参照螺旋桨性能参数的定义，以水翼扫掠面积$ BH $为参考面积，以$ fB $为参考速度，定义推力系数$ {K}_{T} $，表示为：

按相对速度计算，摆翼推进器的工作Re数可表示为：

式中，$ \nu $为运动粘度系数，本文各工况雷诺数$ Re $在$ 0.6\times {10}^{6} $~$ 1.7\times {10}^{6} $范围内。

定义一个周期上横荡运动与艏摇运动瞬时功率的积分均值为平均输入功率，表示为：

定义一个周期上平均推力功率与平均输入功率的比值，称为推进效率，表达式为：

FINE/Marine 是NUMECA 公司为船舶与海洋工程打造的专业CFD 软件包，其集成的由EMN (Equipe Modélisation Numérique)公司开发的ISIS-CFD流动求解器，能够满足水翼多自由度耦合运动的流场的模拟要求，其采用的控制方程为

连续性方程

动量守恒方程

其中，$ \varOmega $为控制体体积，$ A $为控制面，$ \overrightarrow{{V}} $为流场速度，$ {{V}}_{i} $为即$ \overrightarrow{{V}} $的i分量，方向与$ {x}_{i} $轴一致，$ \overrightarrow{{{V}}_{{\rm{d}}}} $为控制体表面运动速度，$ \overrightarrow{{n}} $为控制面法向量，$ {{\tau }}_{ij} $为黏性流体应力张量，$ p $为压强。

算例的湍流模型均为$ k $-$ \omega \left(sst\right) $模型，该湍流模型由Menter提出^[9]，是CFD中普遍应用的双方程涡流粘度模型。该模型考虑了湍流切应力传递的影响，在逆压梯度和分离流中表现优异。

1)网格划分

NUMECA系列软件中的HEXPRESS模块生成的是全六面体非结构网格，网格由体到面生成。在物面附近生成的贴体网格，能准确描述物体的外形。计算域内绝大部分网格单元都接近于长方体，网格的正交性较好，质量较高，能达到很高的计算精度，并能很好的模拟粘性边界层。

本文对计算域网格进行四级加密，在水翼表面及后方位置加大网格密度，以捕捉水翼摆动产生的尾流旋涡，计算域网格划分情况如表2：

网格分级划分区域如图3所示，由于水翼及其运动均关于$ xOy $面对称，因此算例均采用半模计算。

图  3  网格分级划分示意图

Figure 3.  Schematic diagram of grid classification

NUMECA系列软件支持Python语言编写的脚本文件，可对网格划分的操作步骤进行录制，本文所有算例采用录制脚本文件的方法，批量生成相同的网格，因为水翼表面积不同，体表网格数目不同，各种形状的水翼网格数量如表3：

为方便网格捕捉图形，同时方便后期制作，对翼型尾缘使用了直线尖化处理，水翼表面网格即翼型周围网格如图4：

图  4  不同翼型表面网格示意图

Figure 4.  Schematic diagram of grid classification

2)边界条件

算例边界条件如图5所示，水翼表面为物面边界，同时采用壁面函数。水翼前方为远场（无限远来流，速度为零），左、右侧和后方均为零压梯度（第二类Neumann边界条件，即$ \partial p/ \partial n=0 $），上方为描述压力（第一类Dirichlet边界条件，计算域重力沿Z轴的矢量方向，大小为$p=-\rho g\left({\textit{z}}\left(t\right)-{{\textit{z}}}_{0}\right)$），下方为镜像边界（速度场与镜像平面相切，即$ {\overrightarrow{V}}_{mir}\times \overrightarrow{n}=0 $）。

图  5  边界条件示意图

Figure 5.  Schematic diagram of boundary condition

1)实验对比

Schouveiler^[10] 观测了NACA0012翼型上下摆幅$ {Y}_{0}= $0.075 m，摆动频率$ f $=1.2 Hz，St =0.45，最大攻角为20°的工况下的尾涡，本文进行了同工况的数值模拟，其对比如图6所示。从对比可以看出，尾涡形状非常接近。在尾涡刚脱落的阶段，水翼尾涡面发生卷曲，出现“合并”及“缠绕”现象，局部地方的涡量发生集中，连续的涡量最终发展成一列离散的集中涡，旋向相反，交替排列在运动轨迹两侧规则，形成反卡门涡街。

图  6  尾涡对比图

Figure 6.  Diagram of trailing vortex contrast

Schouveiler同时给出了St=0.3，$ f $=0.8 Hz，最大攻角20度（对应摆角幅度约为23.3度）工况下的受力随时间变化曲线，与本文工作的对比如图7，黑线为Schouveiler的实验测量结果，红线为本文数值模拟结果。可以看出，Fx和Fy的对比中，本文结果与Schouveiler的实验结果吻合的较好，Mz由于绝对值较小，较难分辨，二者峰值大致相近。

图  7  文献与仿真数据对比

Figure 7.  Comparison of reference and simulation data

2)网格无关性验证

为研究网格数量对计算结果的影响，使用展弦比AR=2的矩形水翼进行了验证，分别计算了网格数为36万、71万的算例，对比两种网格计算结果如图8：

图  8  不同数量网格计算结果

Figure 8.  Different number of grid results

结果显示，36万网格与71万网格的效率曲线几乎一致，浮动在0.2%左右，推力系数峰值误差为1.3%。36万网格精度与71万差距较小，考虑到计算机性能及计算时长，本文算例采用36万网格划分方式进行计算。

a)效率分析

本文的水翼弦长均为0.2 m，绘制展弦比不同的矩形翼效率高于40%的部分的工况，如图9所示：

图  9  不同展弦比矩形水翼效率曲线

Figure 9.  Efficiency curves of rectangular foil with different aspect ratio

从曲线看出，展弦比对效率产生显著影响，矩形翼展弦比越大，最高效率越大。存在某一进速系数（图中为J=6）附近，该点处，部分较大展弦比的水翼效率基本相同，但效率曲线斜率突然增大，使效率有较快增长。在该点之前的工况下，小展弦比反而有一定的效率优势。该点之后，大展弦比水翼的效率实现反超，呈现出展弦比越大，效率越高的趋势。并且在某一进速区间内，水翼保持着高效运行的状态。但是，展弦比增大引起的效率的增长幅度逐渐变小（图中展弦比AR=10与12的效率曲线几乎一样）。同时，展弦比的增加，引起了最高效率点工况的变化（图中展弦比越大，最高效率点的进速系数越大，使效率曲线峰值向右上方移动），因此，水翼存在最优展弦比，过分加大展弦比对效率提高作用不大。

b) 推力系数分析

展弦比影响的推力系数曲线如图10所示，结果表明，不同展弦比水翼的推力系数曲线走势相似，均先减后增，到达峰值后再减。在效率$ \eta > $40%段(J>4)或起动阶段（J=0附近），展弦比越大，最大推力系数越大。各个曲线依旧存在一个交汇点（图中在J=3.5附近），部分大展弦比水翼在该点的推力系数基本相同，该点之后开始显现大展弦比优势。同时，小展弦比水翼推力系数曲线较为平缓，大展弦比的曲线出现了尖点，峰值位置向右上方移动。展弦比增大引起推力系数的增幅逐渐减小，过分增大展弦比对推力系数增长作用不大。

图  10  不同展弦比矩形翼推力系数曲线

Figure 10.  Thrust coefficient curves of rectangular foil with different aspect ratio

c)流场分析

参照Dong Y^[11]整理的分辨复杂流场中旋涡的判据，本文选取Lambda 2（$ {\lambda }_{2} $判据）方法来显示三维流场中的旋涡。取$ {\lambda _2} = - 0.01 < 0$，显示展弦比AR=3和6的矩形翼在不同工况的尾流场旋涡，如图11：

观察图11(a)、图11(c)发现，进速系数较小的工况，水翼周围产生了大量因摆动产生的小旋涡，与翼尖涡在尾流场中相互融合，发展成清晰的链式涡环结构，与WANG等^[12]在鱼尾鳍尾流场不同展向截面的DPIV实验，以及陈刚等^[13]刚性椭圆翼尾流场的显示结果一致。在进速系数较大的工况下（图11(b)和图11(d)），水翼后方的尾涡层明显变得平直，翼尖涡由于攻角周期性变化，也出现了周期性的消长变化。形成一个涡环的时间约为一个周期，水翼后方计算域长度为7米，因此在进速较高的工况（如图中的6m/s的工况），仅能显示1-2个涡环。

图  11  矩形水翼三维涡环示意图

Figure 11.  Schematic diagram of three-dimensional vortex ring of rectangle foil

本文中，效率曲线及推力系数曲线均出现了汇交点，而小展弦比水翼的三维效应使得效率曲线和推力系数曲线均趋于光滑，本文算例展弦比AR=2、3的水翼均出现了三维效应，与Xia Wu等^[14]在综述中提到展弦比小于4会出现三维效应完全吻合。对于本文的推力系数曲线来说，在雷诺数为$ 7\times {10}^{5} $附近出现了交汇点。此工况下，水翼在某时段攻角很大，前缘涡发生流动分离，产生了动力失速，胡健等^[15]指出，水翼的非定常运动产生的旋涡是推力产生的主要来源，尾涡的反卡门涡街式的排列，使得水翼产生了推力。水翼最大攻角随进速增大而减小，推力系数尖点与效率曲线拐点位置相近，在效率曲线拐点位置，提供推力由非线性升力部分逐步变成线性升力部分，涡系发生变化，仅尾缘涡开始周期性分离，引起推力和阻力均减小，但升力变化不大，表现为效率提高。如果此时继续增大进速，攻角进一步减小，接近最佳升阻比，在最佳升阻比附近的区间保持高效运行。进速进一步增加，升阻比下降，效率随之下降，直到某一进速后，X轴方向受力平衡，效率为零，处于巡航状态。给出本文运动状态下水翼的最大攻角随进速系数变化趋势如图12所示：

图中显示，在本文最大进速系数工况(J=8)时的最大攻角已经接近于零，如再增大进速，有效攻角将减到零甚至变为负攻角，负攻角工作状态也被称作减速状态或能量收集状态。

图  12  最大攻角随进速系数变化示意图

Figure 12.  Diagram of variation of maximum angle of attack with advance coefficients

a)效率分析

为研究根梢比对推进性能的影响，采用展弦比AR=4，而根梢比不同的水翼，绘制不同根梢比的梯形水翼效率大于40%的曲线如图13所示：

图  13  不同根梢比梯形水翼效率曲线

Figure 13.  Efficiency curves of trapezoid foil with different tip ratio

结果显示，根梢比TR=2的梯形水翼最高效率略高于矩形水翼（根梢比可认为是1），根梢比增大引起了效率曲线向右偏移，因此，梯形水翼更适合在较高进速系数的工况工作。进一步增大根梢比，效率曲线开始向右下方偏移，偏移量逐步减少，根梢比TR=5和6的水翼形状接近菱形，效率曲线基本相同。

b)推力系数影响

根梢比影响的推力系数曲线如图14所示，结果表明，梯形水翼的推力系数较矩形水翼大幅降低。曲线走势规律不变，峰值位置略微向右偏移，但增大根梢比使推力系数曲线变得平缓。增大根梢比使得水翼形状接近菱形，因此，推力系数曲线随着根梢比增大而逐渐趋近某一曲线。从推力需求的角度来说，根梢比引起了太多的推力损失，相比梯形水翼，矩形能够在相同展弦比，相同效率下产生更大的推力。当然，可以考虑加大梯形水翼尺寸来使推力达到指定值。

图  14  不同根梢比梯形水翼推力系数曲线

Figure 14.  Thrust coefficient curves of trapezoid foil with different tip ratio

c)流场分析

显示根梢比TR=2和6的梯形水翼在不同工况的旋涡如图15：

图  15  矩形水翼三维涡环示意图

Figure 15.  Schematic diagram of three-dimensional vortex ring of trapezoidal hydrofoil

Xia Wu等^[14]研究指出流体动力特性受其末端形状的影响很大，倾斜的导流边可以产生更复杂的翼尖涡，从结果来看，梯形水翼倾斜的导边使得展向速度加大，从而产生了更为复杂的尾涡层，翼尖涡卷起更加明显。水翼跨中存在的尖点使流体在其前缘发生了分离，产生的旋涡与翼尖涡与尾缘涡相互作用，涡环结构变得复杂，首缘倾斜角度越大，流场越复杂，图中根梢比TR=6的梯形水翼低速工况（图d）尾流场变得十分混乱，尾涡层卷曲十分严重，翼尖涡的强度变大，尾流场变得肥大。在高速工况下，整个尾涡层强度增大，诱导阻力、能量损耗均增大，引起效率下降。

a)效率分析

绘制展弦比AR=4，不同根梢比的前掠、后掠水翼与梯形翼对比的效率曲线如图16所示：

图  16  不同扫掠角度梯形水翼与三角翼效率对比曲线

Figure 16.  Efficiency curves of trapezoid delta hydrofoil with different sweep angles

结果显示，相比于梯形水翼，存在一定的后掠角度能使效率曲线微幅向上偏移，而前掠角度将使效率曲线大幅向下偏移。同时，翼梢的形状也对偏移量有影响，翼梢弦长较小的水翼对前掠和后掠形式较为敏感，偏移显著。

b)推力系数影响

前掠后掠影响的推力系数曲线如图17所示，结果表明，后掠水翼相比于梯形水翼，不仅可以微幅提高效率，还可提高推力系数，而前掠形式依旧不理想，推力系数损失十分明显。不过，根梢比TR=2的前掠水翼与大展弦比矩形水翼十分相似，效率曲线出现了拐点，推力系数曲线出现了尖点，说明前掠水翼仍然保持了矩形水翼的特点，只是尾部形状变化，引起了效率和推力的损失。

图  17  不同扫掠角度梯形水翼与三角翼推力系数对比曲线

Figure 17.  Thrust coefficient curves of trapezoid and delta foil with different sweep angles

c)流场分析

后掠水翼的尾流场如图18所示，相比于梯形水翼的尾涡结构，后掠水翼在低速工况下变得复杂，其倾斜程度更高的导边，使得首缘分离情况更加严重，翼尖涡形状发生改变，与首缘涡相互缠绕，图中根梢比TR=6的水翼，翼尖涡螺旋状加剧，向后发展融合形成涡环。而高速工况下，尾涡层缠绕情况减轻，尾涡强度比矩形水翼大而比梯形水翼小。

图  18  后掠水翼三维涡环示意图

Figure 18.  Diagram of three-dimensional vortex ring of back-swept hydrofoil

前掠水翼的尾流场如图19所示，前掠水翼在低速工况下的涡环结构比较清晰，相比于矩形翼，后缘倾斜以及出现跨中尖点，致使相应涡环发生了更加复杂的折叠，加剧了能量耗散，效率降低。高速工况下，水翼跨中尖点使尾涡层向翼尖卷曲更加明显，尾涡强度明显增加，致使效率降低。综合三者来看，后缘形状对推进器性能产生了更加显著的影响，设计水翼形状时需要综合考虑。

图  19  前掠水翼三维涡环示意图

Figure 19.  Diagram of three-dimensional vortex ring of forward-swept hydrofoil

本文通过数值方法对不同形状的水翼进行研究，发现形状对摆翼推进器的性能有着显著影响。对比结果，得出了不同展弦比、根梢比以及前掠和后掠水翼的效率及推力系数变化规律，同时发现，三维效应会改变水翼表面的流动状况，从而引起水动力性能的改变。研究得到以下几个结论：

1)展弦比会影响水翼的最高效率，而且存在一个摆翼推进器能够一直处于高效运行状态的最佳工作范围，与从文超等^[16]对二维水翼的研究中指出，水翼在一定St数范围内高效运行是一致的。展弦比AR<4的水翼有着明显的三维效应。可通过增加展弦比来提高水翼效率，当展弦比增加到足够长，效率将接近一个饱和状态，展弦比对效率影响减弱，再增大展弦比对提高效率的作用不大，因此存在水翼展弦比优化问题。对正常工作的摆翼推进器来讲，展弦比越大，相应的推力系数越大，且推力系数峰值出现在最高效率点之前的工况。

2)根梢比增大，水翼推力损失较大。但发现根梢比TR=2时，仅推力系数有所下降，并没有影响水翼的最高效率，而且最高效率点随根梢比增大而偏向进速系数较大的工况。如果根梢比过大，会对水翼性能产生不良影响，效率及推力系数均会下降。

3)后掠水翼的效率及推力系数均高于梯形水翼和前掠水翼，倾斜的前缘具有使效率曲线向右偏移的作用，有一定的性能优势。前掠水翼水动力性能较差，从结果对比发现，水翼后缘形状对性能有明显影响，倾斜后缘在低速工况下，涡环会发生折叠，加剧能量耗散，高速工况下后缘尖点的存在使翼尖涡强度增加，效率和推力均下降。

综合几种形状来看，矩形水翼性能最优，后掠水翼次之，梯形及前掠水翼的尾缘形状导致性能较差。本文算例设置水翼后方计算域长度为7m，高速工况下的涡环较长，不能完全展示。如果要获取更详细的尾流场结构，需要增大计算域长度，同时进一步精细化网格。此外, 对于摆翼推进器形状的研究也存在较多研究方向，需要在后续工作中更深入的研究。