伴随着人们对水下仿生推进技术的研究，扑翼推进技术逐渐成为国内外的研究热点，目前国内外对扑翼的水动力特性已经展开了很多研究工作并取得了许多定性和定量的成果[1－3]。Karbasian和Kim[4]对振荡翼失速时的流场结构及涡特性展开数值研究；Hover等[5]基于实验方法对扑翼的不同攻角变化规律对扑翼推进特性的影响规律进行对比，发现扑翼攻角变化规律为余弦时产生的推力和推进效率都是最高；Abbaspour和Ebrahimi[6]基于计算流体力学(CFD)方法对三种常见的扑翼运动方式(纯拍动、纯俯仰、拍动结合俯仰，后文称为扑动)的水动力特性展开数值计算并进行对比分析，发现扑翼进行扑动运动时水动力特性最佳；Wang等[7]对不同系列翼型的扑翼能量提取性能进行数值研究；国内也开展一系列关于扑翼推进特性的研究。杨璞[8]基于有限体积法(FVM)及动网格技术对扑翼的非定常水动力特性展开数值研究；丁浩等[9]关于仿生扑翼推进特性及运动性能展开研究，研究了影响单翼推进性能的各种参数；曹永辉等[10]对混合翼型的推进性能展开研究，发现混合翼型与标准NACA翼型相比具有更优的推进性能；朱建阳[11]对不同扑动轨迹的扑翼气动特性进行数值研究；Jia和Gong等[12]设计了一套基于LabVIEW平台的串列扑翼的力学测量实验装置，同时利用数字粒子图像测试技术(DPIV)技术分析其流场特性。
通过上述研究可发现扑动方式是最优的扑翼运动方式，但是根据相位角及初始姿态的不同扑翼扑动方式可分为以下四种：方式1，初始位置无攻角，扑翼拍动与俯仰运动方向相同；方式2，初始位置无攻角，扑翼拍动与俯仰运动方向相反；方式3，初始位置处有最大俯仰角，扑翼拍动与俯仰运动的方向相同；方式4，初始位置处有最大俯仰角，扑翼拍动与俯仰运动的方向相反。本文主要开展四种扑动方式对扑翼推进性能影响规律的数值研究，从流场角度分析其影响机理。采用计算流体力学(CFD)方法，结合动网格技术，通过用户自定义函数(UDF)模拟扑翼运动。对计算域进行了非结构化三角形网格划分，并划分边界层更好地捕获流场边界层信息。基于FVM，利用计算机程序分别对四种扑动方式下的扑翼的非定常水动力性能进行了数值计算，对流场结构进行分析，得到四种扑动方式对扑翼推进特性的影响规律，为特定情况下扑动方式的选择提供参考。
1 运动模型
1.1 运动描述
扑翼的扑动是一种结合上下拍动及俯仰运动的二自由度运动方式，主要根据上下拍动与俯仰运动的方向分为四种运动方式，图1所示为一个周期内四种扑动方式的二维示意图。
图1中，Hmax为扑翼的拍动幅值，θmax为扑翼的俯仰幅值，运动周期为T，可定义扑翼扑动的基本运动方程如下式[13]：
图1 四种扑动方式示意图
Fig.1 Schematic diagram of four flapping modes
式中，y(t)和θ(t)分别为扑翼的拍动位移和俯仰位移，
为扑动频率，定义拍动频率跟俯仰频率一致。φ为两种运动的相位角，通过改变系数的正负与相位角φ即可满足四种扑动方式。对式(1)求导，即可得到任意时刻扑翼的速度：
式中，
分别为扑翼的拍动速度与俯仰速度。
1.2 运动参数与力学模型
本文中，扑翼假设为刚性翼，翼型选择NACA0012简化翼型，弦长
拍动幅值
俯仰幅值
扑翼的选择中心距前缘距离l=0.2c，对于扑翼运动，定义斯特劳哈尔数
式中，U为无穷远处来流速度，本文中
对于二维扑翼运动，瞬时推力及升力系数分别为：
式中：Ft(t)为水平方向瞬时推力；Fy(t)为竖直方向瞬时升力。平均推力系数及升力系数定义为：
式中，n为计算周期数，本文中n=10，取一个计算周期计算平均推力。
2 数值方法
2.1 控制方程与湍流模型
本文采用 CFD方法展开数值计算，设定流场为二维不可压缩流动的湍流流场，其运动控制方程为：
式中：
为流体的运动速度；
为空间坐标；p为流体压强；t为时间；γ为运动黏性系数；
为湍流黏性系数，k为湍流动能，ε为湍动能耗散率，cμ为常数。
扑翼在扑动过程中，前后缘均会产生涡，同时还伴随着涡的移动与脱落。考虑到Realizablek-ε湍流模型有很好的捕获复杂流场信息的能力，特别适合于这类有涡生成与脱落的复杂流场，因此选择Realizablek-ε作为本文的湍流模型来封闭URANS方程。相应的方程如下：
湍流动能k方程：
耗散率ε输运方程：
式中，右端项分别为生成项、耗散项以及壁面项。式中各常数定义为
近壁衰减函数
壁面项
2.2 计算域与网格生成
计算域尺寸及远场边界条件对计算结果的影响不可忽略，如表1所示为不同计算域尺寸对计算结果的影响，故本文选择计算域尺寸为40c×20c。由于扑翼需要按规律进行扑动，故需要结合动网格技术进行数值计算。为了提高计算效率，特将计算域分为内外域两部分，内域为半径为5c的圆形区域。整个计算域均进行非结构化三角形网格划分，为了提高计算精度，对内域进行更小尺度的网格划分。在翼型周围划分边界层以便更好地捕获边界层流场，第一层网格厚度为0.0001c。计算域及初始时刻有俯仰角的网格划分如图2所示。
表1 不同计算域尺寸对计算结果的影响
Table 1 The influence of different computational domain sizes on numerical results
2.3 数值计算
图2 计算域及网格划分
Fig.2 Computational domain and mesh division
本文利用 ANSYS FLUENT软件进行数值计算，基于Realizablek-ε湍流模型，结合动网格技术展开二维扑翼非定常水动力性能数值计算，入口边界设为速度入口，出口边界为压力出口，壁面为滑移壁面，翼型表面为无滑移壁面，内外域交界处设为Interior。在建立动网格模型时，综合采用弹性光顺法以及局部重构法，对所标记的不满足要求的网格进行自动删除及重新生成。对于弹性光顺，为不影响较远区域的网格，设定弹性系数为0.8。为不影响边界处网格节点分布，设定边界点松弛因子为0.3。对于局部重构，开启 Size Function和Must Improve Skewness，为了保证重新生成的网格尺度较小，需要将最大长度尺度(Maximum Length Scale)调小，最大网格扭曲率取0.65。为保证边界层网格随着扑翼一起运动，因此需要在FLUENT中对网格进行分区，首先将四边形网格标记出来，然后对所标记内域进行分割，最后对所有运动区域进行指定，通过网格预览功能对动网格参数进行合理设置，运动方式1不同时刻动网格更新如图3所示。利用SIMPLE算法来耦合压力场与速度场，利用二阶迎风格式来离散时间[14]，步长为最小网格尺度与流场速度的比值。
图3 不同时刻动网格更新
Fig.3 Mesh update at different moments
3 计算结果
3.1 模型验证
在进行 CFD非定常数值计算过程中，流场网格划分疏密程度以及时间步长的选择都会对计算结果造成很大的影响，因此需要对网格无关性及时间步长无关性进行验证，使计算结果更加准确与可靠。表2为不同网格数目对应的平均推力系数，可以看出当网格数目达到27万的时候，计算结果已经不受影响，因此后文计算中采取27万网格对应的设置参数进行网格划分。表3为不同时间步长所对应的平均推力系数，可以看出当时间步长为0.0005 s时，计算结果收敛，故本文所选时间步长为0.0005 s。基于上述数值无关性验证，选择与麻省理工学院(MIT)[15]相同的工况进行数值计算，对比结果如表4所示，发现随着St的增大相对误差逐渐减小，最大相对误差为12.212%，最小相对误差为5.232%，由于St=0.1时推力系数较小，故出现较大误差在许可范围内，因此所建计算模型的计算结果是可靠的。
表2 网格无关性测试
Table 2 The grid-independence tests
表3 时间步长无关性测试
Table 3 The time step-independence tests
表4 计算结果与实验结果对比
Table 4 The comparison between numerical and experimental results
3.2 扑动方式对推力及升力的影响
随着St的增大，扑翼所产生的推力系数是逐渐增大的。本文选择斯特劳哈尔数St=1对应的工况展开数值计算，即拍动频率f= 2 Hz 。
图4为两个周期内4种不同扑动方式对应的推力系数随时间变化曲线，可以看出：对于方式 1，扑翼从平衡位置向上开始扑动，推力逐渐减小，随后产生阻力。在扑翼运动到最大拍动幅值之前，阻力达到最大值时，随后扑翼开始产生推力，在扑翼回到平衡位置之前，推力达到最大值。随后扑翼开始向反方向运动，推力又开始减小，随后产生阻力。在扑翼达到平衡位置之前，推力再次产生峰值，但是比前半个扑动周期产生的峰值略小。该运动方式扑翼总体产生推力，平均推力系数为5.9；对于方式 2，扑翼总体不产生推力，随着扑翼向最大拍动幅度扑动，扑翼开始产生推力并逐渐增大，随后开始下降，最后产生阻力。在扑翼达到一次平衡位置时，扑翼阻力达到最大值；对于方式 3，可以看出推进特性最好，在整个运动周期都不产生阻力。在扑翼下拍过程中推力第一次达到峰值，随后下降到推力最小值。在后半个扑动周期，推力再次达到峰值，但是比前一个峰值小。通过对比峰值出现时刻发现，方式3的推力变化相对于方式1有明显的滞后；对于方式 4，总体产生阻力，且推力变化较方式3有明显滞后。
图4 不同扑动方式推力系数变化曲线
Fig.4 Instantaneous thrust coefficient curves with different flapping modes
为了更深入探索推力系数变化规律，对推力系数进行快速傅氏变换(FFT)转换，为了对比明显，取分析频率段为0 Hz～15 Hz，分析其频域特性，如图5所示。可以看出四种扑动方式对应的推力系数变化频率几乎一致，主频率均为4 Hz，是扑动频率的2倍，这也与已有文献的结果是一致的。方式1与方式4对应的主幅值最大，且近乎相等，方式2幅值最小，这也与图4的结果是对应的。由此可以看出，对于推力特性，四种扑动方式的推进特性对比明显，推力变化均近似呈余弦型变化规律，且变化频率近乎一致，主频率为扑动频率的2倍。方式3推力特性最好，方式1次之，方式2与方式4总体不产生推力。
图5 不同扑动方式推力系数频域特性分析
Fig.5 Frequency domain analysis of thrust coefficients with different flapping modes
图6所示为2个周期不同扑动方式对应升力系数随时间变化曲线。可以看出：四种升力系数在一个周期只有两个峰值(一正一负)，而且远大于推力系数，除了方式 3，其他扑动方式都产生随时间余弦型变化的升力，四种扑动方式产生的平均升力都几乎为0。对于方式 1，升力系数峰值最大，在初始平衡位置时，扑翼产生负值最大的升力系数，前半个扑动周期结束后，扑翼产生正值最大的升力系数；对于方式2，峰值约为方式1的一半，且峰值出现时间与方式1略有提前；对于方式3，升力系数的峰值最小，约为方式 1峰值的1/4，且峰值出现时间最为靠后；对于方式 4，最大峰值出现时间最靠后，最大拍动位置附近出现峰值，且峰值也较大。
图6 不同扑动方式升力系数变化曲线
Fig.6 Instantaneous lift coefficient curves of different flapping modes
图7所示为对升力系数进行频域特性分析，分析频率段为0 Hz～10 Hz。可以看出四种扑动方式对应的升力系数变化主频率为2 Hz，这与扑翼的扑动频率是相同的，也与已知文献是一致的。与推力系数相比，升力系数的幅值较大。对于扑翼运动，我们通常更关心扑翼向前产生的推力大小，但是对于有些航行器，扑翼所产生的竖直方向力也有实际应用价值，因此合理选择扑翼运动方式可能会产生意想不到的效果。
3.3 扑动方式对流场结构的影响
由于四种扑动方式在一个扑动周期产生的推力与升力对比非常明显，方式 1、3总体会产生推力，方式 2、4整体不产生推力。由于通常情况下扑翼所产生的推力更能引起我们的兴趣，故本小节仅对方式1和方式3的流场进行深入对比分析，揭示其对推力影响机理。
图7 不同扑动方式升力系数频域特性分析
Fig.7 Frequency domain analysis of lift coefficient of different flapping modes
如图8所示为两种扑动方式1个周期不同时刻的流场流线图，图9为不同时刻涡量等值线图。在前 1/4个周期，对于方式 1，随着扑翼从初始位置向上拍动及俯仰运动，拍动及俯仰速度减小，前缘逐渐形成前缘逆时针涡，后缘也逐渐形成后缘逆时针涡，因为起初新的涡尚未形成，而前面产生的涡已经向下游运动且强度越来越弱，逆卡门涡的推进效果变弱，同时扑翼的迎流面积逐渐增大，扑翼对流体表现更多地为阻挡效果，故前段时间扑翼产生的推力不足以抵消阻力，而随着涡的形成，强度逐渐变大，扑翼所产生的推力已经足以克服阻力，故表现为推力，因此推力先减小后增大；对于方式3，扑翼从初始位置向下拍动，同时向下做俯仰运动，此时拍动速度减小，但俯仰速度增大。尾缘的逆时针涡随着俯仰运动增强并向后分离，前缘逐渐形成顺时针涡，此时推力有所增加。随着扑翼拍动速度的减小，涡的强度也逐渐减小，故推力也逐渐减小，但由于在拍动速度最小时翼是水平的，迎流面积最小，此时逆卡门涡街的推力效果虽然最弱但足以抵消扑翼的阻力，故推力始终为正值。
图8 方式1、3不同时刻流线图(上为方式1，下为方式3)
Fig.8 The streamline of mode one and three at different moments (upper: mode one; bottom: mode three)
图9 方式1、3不同时刻涡量等值线图(上为方式1，下为方式3)
Fig.9 The vorticity contor of mode one and three at different moments (upper: mode one; bottom: mode three)
在第二个 1/4周期，对于方式 1，扑翼向下拍动并做俯仰运动，拍动速度及俯仰速度逐渐增大，此时前缘逆时针涡形成并从下表面至尾缘脱落，尾部逆压梯度增强，同时扑翼迎流面积减小引起阻力减小，因此推力持续增加。最后前缘开始形成逆时针涡，而尾涡已经脱落，因此推力又有所减小；对于方式4，扑翼开始向上拍动同时做向下俯仰运动，此时拍动速度增大，前缘形成的逆时针涡逐渐增强并从上表面向后运动，扑翼尾流场形成的逆卡门涡街的推力效果明显，同时尾缘形成顺时针涡，逆压梯度增加，故推力持续增加。
在后半个扑动周期，两种扑动方式所产生的推力均跟前半个扑动方式有相同的变化规律，故不再赘述。但是通过观察流线图及涡量等值线图可发现，扑翼在前后两个扑动周期中所形成的前缘涡并不一样，从尾缘的脱落位置也不一样，前半个扑动周期前缘形成逆时针涡对逆卡门涡街有增强作用，而后半个扑动周期前缘形成顺时针涡对逆卡门涡街有抵消作用，且导致涡量大小不同，故前后两个半周期推力峰值不同，扑动方式3较为明显。对比两种方式同一时刻，由于扑翼的姿态不同，故涡的形状及强度也不一样。总体来讲，扑动方式3在扑动过程中翼型的姿态及所形成的涡对扑翼的推进效果更为明显，因此其产生更大的平均推力。
4 结论
基于所建立二自由度水下刚性扑翼运动模型，结合动网格技术，基于FVM求解URANS方程，对水下扑翼的四种扑动方式的推进性能进行对比研究，对其流场影响机理进行对比分析。结果表明：
(1) 四种扑动方式所产生的推力在一个扑动周期均近似为余弦规律变化，主变化频率均为扑动频率的2倍。扑动方式3在整个扑动周期均产生推力，能够产生最大的平均推力，方式1次之，且方式1跟方式3的推力峰值几乎相等；扑动方式2与方式4不产生平均推力，且推力峰值几乎相等；扑动方式2跟扑动方式3产生推力峰值的时刻几乎相等。
(2) 四种扑动方式所产生的升力在一个扑动周期均为一个变化周期，四种方式均近似为余弦变化规律，主变化频率等与扑动频率；四种扑动方式所产生的平均升力均几乎为0，方式1产生的升力峰值最大，方式3产生的升力峰值最小。
参考文献：
[1]Olivier M, Dumas G.A parametric investigation of the propulsion of 2D chordwise-flexible flapping wings at low Reynolds number using numerical simulations [J].Journal of Fluids and Structures, 2016, 63: 210―237.
[2]Moriche M, Flores O, García V M.Three-dimensional instabilities in the wake of a flapping wing at low Reynolds number [J].International Journal of Heat and Fluid Flow, 2016, 62: 44―55.
[3]刘健, 蒋永, 吴金华.基于 iDDES 的双三角翼大迎角非定常涡破裂特征分析[J].工程力学, 2016, 33(4):241―249.Liu Jian, Jiang Yong, Wu Jinhua.Analysis on characteristics of unsteady vortex breakdown flows around a double-delta wings at high incidence based on iDDES method [J].Engineering Mechanics, 2016, 33(4): 241―249.(in Chinese)
[4]Karbasian H R, Kim K C.Numerical investigations on flow structure and behavior of vortices in the dynamic stall of an oscillating pitching hydrofoil [J].Ocean Engineering, 2016, 127: 200―211.
[5]Hover F S, Haugsdal ?, Triantafyllou M S.Effect of angle of attack profiles in flapping foil propulsion [J].Journal of Fluids and Structures, 2004, 19(1): 37―47.
[6]Abbaspour M, Ebrahimi M.Comparative numerical analysis of the flow pattern and performance of a foil in flapping and undulating oscillations [J].Journal of Marine Science and Technology, 2015, 20(2): 257―277.
[7]Wang Y, Sun X, Huang D, et al.Numerical investigation on energy extraction of flapping hydrofoils with different series foil shapes [J].Energy, 2016, 112: 1153―1168.
[8]杨璞.扑翼的非定常水动力特性数值研究[J].计算机仿真, 2014, 31(8): 372―376.Yang Pu.Numerical research on unsteady hydrodynamic property of flapping-foil [J].Computer Simulation, 2014,31(8): 372―376.(in Chinese)
[9]丁浩, 宋保维, 田文龙.水下仿生扑翼推进性能分析[J].西北工业大学学报, 2013, 31(1): 151―155.Ding Hao, Song Baowei, Tian Wenlong.Exploring propulsion performance analysis of bionic fapping hydrofoil [J].Journal of Northwestern Polytechnical University，2013, 31(1): 151―155.(in Chinese)
[10]曹永辉, 朝黎明, 丁浩.混合翼型对扑翼推进性能影响分析[J].上海交通大学学报, 2016, 50(8): 1228-1232.Cao Yonghui, Chao Liming, Ding Hao.Propulsion performance of complex flapping foils [J].Journal of Shanghai Jiao Tong University, 2016, 50(8): 1228―1232.(in Chinese)
[11]朱建阳.扑动轨迹对扑翼气动特性影响的数值研究[J].工程力学, 2016, 22(1): 246―251.Zhu Jianyang.Numerical study on the effect of kinematics of a flapping wing on its aerodynamic performance [J].Engineering Mechanics, 2016, 22(1):246―251.(in Chinese)
[12]Jia B B, Gong W Q.Force and flow field measurement system for tandem flapping wings [J].Experimental Techniques, 2016, 40(6): 1495―1509.
[13]Lighthill M J.Aquatic animal propulsion of high hydromechanical efficiency [J].Journal of Fluid Mechanics, 1970, 44(2): 265―301.
[14]Lu K, Xie Y H, Zhang D.Numerical study of large amplitude, nonsinusoidal motion and camber effects on pitching airfoil propulsion [J].Journal of Fluids and Structures, 2013, 36: 184―194.
[15]Triantafyllou M S, Techet A H, Hover F S, et al.Review of experimental works in biomimetic foils [J].IEEE Journal of Oceanic Engineering, 2004, 43(3): 558―594.