0 引　言

在自然界中，蜻蜓、蜜蜂等昆虫的扑翼结构是最常见的串列扑翼系统，学者们从中获得灵感，采用多翼串列的方式来提升扑翼系统的推进性能。于宪钊等^[1]分析了前翼拍动后翼静止的串列翼系统水动力；Jones等^[2]分析了串列水翼的推进性能并进行了相关试验研究；Broering等^[3]采用数值模拟方法研究了串列扑翼的气动力性能；石成明等^[4]采用气动力试验方法研究了串列翼型所产生的推力与前后翼距离之间的变化规律；李永泽等^[5]研究了3种串列双翼的气动布局，分析了最优的前后翼布局方案。许国冬等^[6]针对典型推进运动模式与能量采集模式，分析了串列双翼之间涡的干扰规律。

尽管国内外针对串列翼的气动力和水动力性能已开展许多研究，但过往研究中的串列翼运动形式大多为正弦运动。自然界生物的扑翼运动轨迹往往是更为复杂的非正弦运动，目前尚未看到针对非正弦运动串列水翼推进表现的研究。在串列扑翼式水下航行器的研发驱动下，本文采用计算流体力学方法研究非正弦升沉运动下串列水翼的推进性能，为水下航行器的扑翼运动模式选择提供指导。

1 计算方法 1.1 扑翼运动的描述

本文中非正弦运动通过椭圆形轨迹实现，其中，扑翼的升沉运动位移方程定义如下^[7]：

$ h(t) = {h_0}\frac{{S\cos (\omega t)}}{{\sqrt {{S^2}\cos {{(\omega t)}^2} + \sin {{(\omega t)}^2}} }}。$

(1)

其中：S为扑翼升沉运动的非正弦参数；h₀为升沉运动幅值； $\omega = 2{\text{π }}f$ ，f为拍动频率；t为时间。 $S = 1$ 时，扑翼运动即为正弦运动，不同S值对应的扑翼运动轨迹如图1所示。

当两扑翼前后串列排布时，前翼和后翼的升沉运动可分别表示为：

$ {h_1}(t) = {h_0}\frac{{S\cos (\omega t)}}{{\sqrt {{S^2}\cos {{(\omega t)}^2} + \sin {{(\omega t)}^2}} }}，$

(2)

$ {h_2}(t) = {h_0}\frac{{S\cos (\omega t - \varepsilon )}}{{\sqrt {{S^2}\cos {{(\omega t - \varepsilon )}^2} + \sin {{(\omega t - \varepsilon )}^2}} }}。$

(3)

其中：h₁(t)为串列排布两翼中前翼的运动方程，h₂(t)为后翼的运动方程；S为扑翼升沉运动非正弦参数；ε为前后两翼升沉运动的相位差。

1.2 物理量的定义

扑翼在运动时会受到不同方向上的水动力和水动力矩。本文对串列翼系统的前后水翼力学特性分别讨论，作用在单独进行升沉运动水翼上的瞬时推力系数C_T和瞬时升力系数C_Y可表达为：

$ {C_T} = \frac{{{F_X}(t)}}{{\dfrac{1}{2}\rho c{U^2}}}，$

(4)

$ {C_Y} = \frac{{{F_Y}(t)}}{{\dfrac{1}{2}\rho c{U^2}}}。$

(5)

其中：F_X (t)为作用在单独水翼上沿x轴的瞬时推力；F_Y (t)为作用在单独水翼上沿y轴的瞬时升力；ρ为流体的密度；U为流场来流速度。

则平均推力系数C_T,mean和平均功率系数C_P,mean可表达为：

$ {C_{T,mean}} = \dfrac{{\dfrac{1}{T} \displaystyle \int_0^T {{F_X}(t){\rm{d}}t} }}{{0.5\rho c{U^2}}}，$

(6)

$ {C_{P,mean}} = \frac{{{P_{input}}}}{{0.5\rho c{U^3}}}。$

(7)

其中：T为扑翼升沉运动周期；P_input为扑翼运动时所需功率。

对于升沉运动水翼，一个运动周期所需的功率可以通过下式计算：

$ {P_{input}} = \frac{1}{T}\displaystyle \int_0^T {{F_Y}(t) \cdot v(t)} {\rm{d}}t。$

(8)

其中：v(t)为水翼沿y轴轴向速度。

单独水翼的推进效率 $ \eta $ 可表达为：

$ \eta = \frac{{{C_{T,mean}}}}{{{C_{P,mean}}}}。$

(9)

1.3 数值方法与计算模型

本文采用基于有限体积法的商用CFD求解器STAR-CCM+进行串列水翼的水动力分析。其中，水翼的运动利用重叠网格技术来实现。数值模型的计算域及边界条件如图3所示。其中，分析的水翼为NACA0012标准翼型，翼弦长 $c = 0.1\;{\text{m}}$ 。方形区域为背景域，尺寸设置为70c×70c；中间圆形区域为重叠域，半径为1c。根据Chao等^[8]的研究，当扑翼和计算域边界的距离大于20c时，流域尺寸对水翼水动力的影响可以忽略不计，因此，可以判断本文中的计算域尺寸已足够大。

1.4 数值收敛性分析

针对数值模型中的时间步长和网格尺寸进行收敛性分析，具体计算参数取值为 $ {h_0} = 0.7c $ ， $Re = 1 \times {10^4}$ ， $ St = 0.35 $ ， $ S = 0.25 $ ，相位差 $ \varepsilon = 0 $ 。其中，Re为雷诺数，其定义为：

$ Re = \frac{{Uc}}{v}。$

(10)

式中：U为来流速度，v为水的动力粘度。St为斯特劳哈尔数，反映流体非定常运动的相似性，其定义为：

$ St = \frac{{2{h_0}f}}{{{U_0}}}。$

(11)

式中，f为运动频率。

图4为不同时间步长下前后水翼推力系数时历曲线的比较。计算中选取的时间步长分别为dt=T/1000，T/2500，T/4000，对应计算模型的网格数目为24.5万。可以看出，不同时间步长下水翼的推力曲线吻合较好，故时间步长dT=t/2500用于后续计算。图5为不同网格密度下前后水翼推力系数时历曲线的比较。对应3个不同的网格密度，扑翼表面分别具有200，400，600个节点，对应计算模型的网格数目分别为6.8万、24.5万和54.1万，计算中采用的时间步长dT=t/2500。可以看出，不同时间步长下水翼的推力曲线总体吻合较好，表明当前数值模型具有较好的网格收敛性，在后续计算中选取网格数目为24.5万的计算模型。

1.5 数值算法的验证

图6为本文数值结果与Heathcote等^[9]试验结果以及Y.L.Wang等^[10]和Young等^[11]数值结果的对比。其中，计算翼型为NACA0012水翼，弦长 $c = 0.1\;{\text{m}}$ ，雷诺数 $Re = 2 \times {10^4}$ ，运动幅值 ${h_0} = 0.175c$ ，水翼运动形式为单翼的正弦纯升沉运动。横坐标为无量纲所减频率k，定义为 $ k = 2{\text{π}} f c/{U_0} $ 。可以看出，本文数值结果与文献[9]试验结果吻合较好。图7为本文数值结果与Kinsey等^[12]结果的对比。其中，文献中的计算结果也被许国冬等^[6]采用，用于验证数值模型。其中，计算翼型为NACA0015水翼，弦长 $c = 0.1\;{\text{m}}$ ，雷诺数 $Re = 2 \times {10^4}$ ，运动幅值 ${h_0} = 1.0c$ ，水翼运动形式为串列水翼的俯仰升沉复合运动。可以看出，本文数值结果与文献数值结果吻合较好。总体而言，当前数值方法在本文中的应用是合理可行的。

2 结果与分析

主要讨论非正弦运动对串列升沉水翼推进性能的影响，其中无特殊说明，计算参数 ${h_0} = 0.7c$ ， $Re = 1 \times {10^4}$ ， $St = 0.35$ ， ${L_x} = 7.5c$ ，保持不变。

2.1 非正弦升沉系数S对瞬时水动力的影响

采用前后翼相位差 $\varepsilon = 90^\circ $ 的结果进行分析。如图8（a）所示，对于前翼的瞬时推力系数曲线：从峰值来看，非正弦系数从S=1增大或减小，都会使前翼瞬时推力系数曲线峰值增加；从周期来看，增大非正弦系数时周期基本不变，而减小非正弦系数时原有周期发生变化，峰值出现位置由水翼运动的平衡位置（0.25T，0.75T）变化为水翼运动的极端位置（0T，0.5T，1T）。如图8（b）所示，对于后翼的瞬时推力系数曲线：从峰值来看，随着非正弦系数S的减小，后翼推力系数曲线峰值增大，随着非正弦系数S的增大，峰值先减后增；从周期来看，由于相位差的存在，后翼周期应比前翼推迟0.25T，而当S=1时，峰值却提前出现，这是前翼的尾流作用所致。

如图9（a）所示，对于前翼的升力系数曲线：从峰值和谷值来看，从S=1增大或减少非正弦系数均可增大或减小升力系数曲线的峰值或谷值，且从S=1减小时相比于增大时峰值增大或谷值减小的程度更大；从周期来看，当S>1时周期基本不变，而当S<1时，峰值或谷值出现位置由水翼运动的平衡位置（0.25T，0.75T）变化为水翼运动的极端位置（0T，0.5T，1T）。如图9（b）所示，对于后翼的瞬时升力系数曲线：从峰值和谷值来看，从S=1增大或减小非正弦系数同样会使升力系数曲线峰值更大，谷值更小；从周期来看，由于相位差的存在，后翼的周期性变化比前翼慢0.25T。

2.2 非正弦升沉系数S对涡结构的影响

图10和图11分别为前后两翼在 $t = 0.5{{T}}$ 时刻不同非正弦升沉运动时周围流场的涡量云图。此时前翼处于最大负向位移位置，后翼处于向下运动的平衡位置。将水翼的上表面称为吸力面，下表面称为压力面，水翼形成升力的原因是吸力面与压力面的压强差。

前翼在非正弦系数S=1时，一组尾涡由2个正向旋转、1个负向旋转的涡组成。在脱落过程中，正向旋转的涡迅速变弱至消失。水翼首部有一对旋向相反的涡，且带着一个未与水翼表面接触的负向旋涡。当这组涡向水翼尾部运动时，涡强度逐渐减弱。由于上表面的涡弱，上下表面压力差不大，升力无明显变化。当S<1时，尾涡组成与S=1时相似，当S减小时，尾涡中一个正向旋转的涡强度增强，推力明显增大。当S>1时，一组尾涡由2对旋向相反的涡组成，且脱落后负向旋涡迅速减弱，当S增大时，上表面脱落的负向旋涡强度增大。

对于后翼在非正弦系数S=1时，有一个负向旋涡在首端脱落，有一个正向旋涡在尾端脱落。负向旋涡比正向旋涡强度更强，有一组旋向相反的涡在尾流中，观察到有前翼的尾涡运动到后翼附近。当S<1时，首端负向旋涡分为2个，尾端正向旋涡和尾涡中的负向旋涡随S的减小而减弱。当S>1时，随着S的增大，尾端的正向旋涡强度增强，首端负向旋涡强度减弱。

2.3 非正弦升沉系数S对推进表现的影响

从图12（a）可知，对于前翼，从非正弦系数S=1增大或减少，均能增加前翼的平均推力系数。在S=0.5时前翼平均推力系数在相位差为90°处变化明显，在此非正弦系数下，后翼对前翼的推力影响较大。从图12（b）可知：在大多数相位差下，从S=1增大或减小非正弦系数可以增强后翼平均推力系数，但在某些相位差下，非正弦运动的平均推力系数小于正弦运动的平均推力系数。如相位差为0时，S=0.5的平均推力系数小于S=1的平均推力系数。总体来看，随着相位差的增加，平均推力系数先增后减。

从图12（c）可知，对于前翼的推进效率，当S>1时，推进效率随S的增大而减小，当S<1，推进效率随S的减小而增大。可以观察到在S=0.5时，前翼推进效率曲线在不同的相位差下，数值发生较明显变化，推测是此非正弦系数和相位差时后翼对前翼的影响较大所致。从图12（d）可知，对于后翼的推进效率同样在S>1时，推进效率随S的增大而减小，随相位差的增大先减后增，在S<1时，随S的减小而增大，随相位差的增大先增后减。

3 结　语

本文采用计算流体力学方法和重叠网格技术研究了非正弦运动对串列升沉水翼推进性能的影响，得到如下结论：

1）通过不同时间步长和网格尺寸的数值结果比较，验证了数值模型的收敛性；通过与公开发表试验和数值结果的比较，验证了数值模型的计算精度。

2）非正弦升沉运动相比于正弦升沉运动会增大扑翼的瞬时推力峰值和升力峰值，而且，非正弦系数S的量值与S=1的偏差越大，瞬时峰值增大程度越明显。

3）非正弦系数可以改变水翼表面涡的脱落模式，当S<1时，一组尾涡由2个旋向相同、一个旋向相反的涡组成；当S>1时，一组尾涡由2对旋向相反的涡组成。

4）除个别少数相位差对应的运动状态，与正弦运动相比，当非正弦系数从S=1增大或减小时，非正弦运动下的前翼和后翼平均推力增大。

5）除个别少数相位差对应的运动状态，与正弦运动相比，当非正弦系数从S=1增大时，前翼和后翼的推进效率减小；当非正弦系数从S=1减小时，前翼和后翼的推进效率增大。除此之外，在此基础上改变双翼相位差可以进一步使串列水翼达到最佳推进效率，这将为串列双翼水下航行器推进系统的设计提供指导。