自然界中的鱼类有着优异的游动能力，它们经历长期的自然选择，进化出了各式各样的身体形态及游动模式[1]，例如擅长复杂环境下作高机动性游动的带状蝴蝶鱼，擅长远距离巡游的金枪鱼，擅长短时间内加速游动的梭子鱼．在身体-尾鳍推进模式下，鱼类能够进行快速地定向巡游及快速机动，这种运动模式下鱼类的身体与尾鳍之间存在涡系干扰，表明身体的脱落涡能够提升尾鳍的推进性能，但是这种机理尚未被完全掌握．针对活鱼或者仿生鱼鳍的水池试验通常能够获得运动学和动力学特征，以及流场的速度矢量切片，但是难以获得表面压力分布、水动力性能及流场精细的涡结构特征，因此难以开展机理研究[2]．例如文献[3-4]试验研究了柔性扑翼的水动力性能．随着计算流体力学(CFD)方法的发展，采用计算机编写程序能够对鱼类游动过程中的流场的速度、压力等开展详细的数值计算．文献[5]建立了仿鱼类游动模式下的摆动翼运动模型，计算结果表明仿鱼类运动模式下摆动翼的推力系数更大．文献[6]采用面元法(BEM)数值计算了仿生金枪鱼在均匀来流中的水动力性能，结果表明身体脱落的涡能够与尾鳍产生的涡相互干扰，提升尾鳍脱落涡的强度．文献[7]数值计算了仿生金枪鱼在均匀来流中游动时的水动力性能，结果表明躯体产生的涡能够与尾鳍产生的涡相融合，提升尾鳍的推进性能．考虑鱼类的游动是身体-尾鳍与流体相互作用下的自主推进游动，文献[8]研究了尾鳍运动参数对仿生鱼巡游过程中水动力性能的影响．文献[9]数值计算了仿生金枪鱼的加速-巡游过程，结果表明在巡游阶段尾鳍能够吸收躯体产生的涡中的能量，减少功率消耗．文献[10]数值计算了仿生金枪鱼的摆动-滑行游动过程．文献[11]研究了摆动翼在振动半圆柱尾流场中的水动力性能，结果表明半圆柱涡能够提升摆动翼的推力系数或者推进效率．针对整个仿生鱼的数值计算中，已有研究大都发现尾鳍能够从躯体产生的涡中吸收能量，从而提高推进性能，然而针对这种现象的机理研究较少．本研究基于CFD方法，建立了串列翼数值计算模型，包括一个大弦长的前置翼和一个后置的小弦长摆动翼，数值计算了摆动翼在均匀来流、前置翼刚性不运动及柔性波动运动的尾流场中的水动力性能，由流场涡结构的演变特征分析摆动翼水动力性能提升的机理，为提升仿生推进器的推进性能提供理论参考．1 数值计算模型1.1　数值计算模型串列翼数值计算模型的总长为1.15 m，其中上游的水翼为NACA 0018翼型，弦长Cmf为0.95 m，定义为前置翼；下游水翼为NACA 0012翼型，弦长Ccf为0.15 m，定义为摆动翼；前置翼与摆动翼的间隙为0.05 m．为了研究前置翼对摆动翼水动力性能的影响及影响机理，定义为模式1、模式2、模式3三种运动模式．模式1为摆动翼在均匀来流中的平移-摆动运动；模式2为摆动翼在刚性不运动的前置翼的尾流场中进行平移-俯仰运动；模式3为摆动翼在前置翼进行柔性波动运动的尾流场中进行平移-俯仰运动．模式2和模式3的数值计算模型相同，以模式2的数值计算模型为例，如图1所示，图中：黑色翼型为前置翼；红色翼型为摆动翼；随体坐标系o-xy位于前置翼前缘点向后0.3 m处．数值计算过程中前置翼和摆动翼的运动在随体坐标系下进行．10.13245/j.hust.230402.F001图1串列翼数值计算模型1.2　串列翼的运动模型串列翼的运动模型包含两个部分：前置翼的柔性波动运动，摆动翼的平移-摆动运动．由鱼类的身体-尾鳍(BCF)游动模式定义前置翼的柔性波动运动，设定前置翼的前0.3 m为刚性体，前置翼的其余部分为柔性体，进行侧向的柔性波动运动[6]，有y(xf,t)=a(xf)sin(kxf-2πft)，(1)式中：xf为x轴坐标；t为时间；a(xf)为振幅包络线函数；k为波数，定义为k=2π/λ，其中λ为波长，本研究取为1.25 m；f为运动频率．模式3与鱼类的身体-尾鳍游动模式相同，对于这种游动模式，a(xf)为二次多项式形式[6]，具体为a(xf)=a1xf+a2xf2，(2)式中a1和a2均为常系数．摆动翼关于前缘点的平移-摆动运动[6-7,9]为yf(t)=Apsin(2πft);θf(t)=θmsin(2πft-Δθ), (3)式中：Ap为平移幅值；θm为摆角幅值；∆θ为摆动运动与平移运动之间的相位差．针对摆动翼推进性能的研究表明：当∆θ为90°时，推进性能最优．考虑到本研究重点在于涡-翼的相互干扰，∆θ值恒为90°．在模式3中，摆动翼在前置翼的带动下进行平移运动，应对式(3)中摆动翼的平移运动进行改进．由于前置翼与摆动翼之间存在间隔，因此将摆动翼的前缘点坐标代入式(1)得到摆动翼的平移运动．此外，由于摆动翼的摆动运动与平移运动是同步的，且存在恒定的相位差，因此应对式(3)中的摆动运动进行改进，具体为yf(t)=a(x1)sin(kx1-2πft);θf(t)=θmsin(kx1-2πft-Δθ),式中：x1为摆动翼前缘点的x轴坐标；a(x1)为平移运动幅值，即Ap．由式(1)与式(2)给出的运动模型，设置参数a1=3.32×10-3，a2=-0.133，θm=25°．图2给出了模式2和模式3的一个游动周期内的前置翼和摆动翼的运动形态．10.13245/j.hust.230402.F002图2游动周期内的前置翼和摆动翼的形态变化1.3　水动力系数及推进效率在数值计算过程中，摆动翼受到流体作用的纵向力Fx、侧向力Fy和关于前缘点的力矩Mz．从图1可以看出：Fx的正方向指向摆动翼的尾缘，因此Fx表现为推力时为负值，表现为阻力时为正值，为了便于分析，定义摆动翼的推力F1=-Fx．推力F1、侧向力Fy和力矩Mz的无因次化系数为Ct=2F1/(ρCcfUin2);Cy=2Fy/(ρCcfUin2);Cm=2Mz/(ρCcf2Uin2),式中：Uin为入口速度；ρ为流体的密度．相应地，平均推力系数C1、平均侧向力系数C2和平均力矩系数C3定义为C1=∫tt+nTCtdt/(nT);C2=∫tt+nTCydt/(nT);C3=∫tt+nTCmdt/(nT),式中：n为采样周期数；T为运动周期．摆动翼的平均输入功率Pin和系数C4为：Pin=1nT∫(Fyyf'(t)+Mzθf'(t))dt;C4=2Pin/(ρCcfUin3)．摆动翼的推进效率ηcf=C1/C4．2 数值计算方法2.1　流体控制方程对于不可压缩的黏性流体，其流动满足连续性方程和雷诺平均纳维-斯克托斯(RANS)方程，即∂ui/∂xi=0；∂(ρui)/∂t+∂(ρuiuj)/∂xj=-∂p∂xj+∂∂xjμ∂ui∂xj-ρui'uj'+Si,式中：i，j=1，2，3；ui和uj为流体速度分量的时均值；p为流体压力时均值；μ为流体的动黏滞系数；ui'uj'为速度脉动分量乘积的时均值；Si为广义源项．湍流模型采用RNG k-ε模型[3]，选择基于压力的求解器，由标准壁面函数法对进壁面流体的运动进行处理，计算物体表面的y+值分布在30~80范围内．采用SIMPLE算法对压力和速度进行耦合，梯度项选择基于单元体的最小二乘法，压力项选择二阶，动量项选择二阶迎风格式，时间项采用一阶隐式格式离散．2.2　计算域及网格划分策略计算域及网格划分如图3所示．计算域为80Ccf×40Ccf的矩形，关于ox轴左右对称．计算域的入口和左右两侧设置为速度入口，出口设置为压力出口．前置翼的前缘点到计算域入口的距离为30Ccf．串列翼的四周尺寸为10Ccf×8Ccf的计算域定义为内部计算域，对其进行三角形网格划分，网格最大尺寸为0.5 mm．整个计算域中除去内部计算域的区域定义为外部计算域，对其进行四边形网格划分；摆动翼的尾流场中尺寸为10Ccf×8Ccf的区域由尺寸为0.5 mm的四边形网格进行离散；外部计算域的其他区域网格以等比数列形式快速地向四周拉伸．这种网格划分策略能够对摆动翼周围的流场特征进行有效的求解和捕捉，在提高计算精度的同时减少网格数量，降低计算成本．当对模式210.13245/j.hust.230402.F003图3计算域及网格划分策略与3进行数值计算时,整个计算域的网格数量为3.676 74×105．当对模式1进行数值计算时，须要将前置翼从计算域中去掉，因此整个计算域的网格数量(3.396 10×105)会有略微的变化．本研究采用上述网格划分策略，以模式1为验证算例，分别划分数量为1.968 92×105，3.396 10×105，5.856 68×105的三套不同网格，计算了摆动翼在Ap=0.42 Ccf，f =1 Hz，θm=25°及Uin=0.4 m/s时的C1，分别为0.835 61，0.827 12和0.825 93，呈收敛趋势，表明本研究用于数值计算的网格满足收敛性要求．时间步长设为T/1 000以保证库朗数小于1，相应的数值计算方法验证可参考作者前期的研究工作[5]．采用弹簧光顺法和局部网格重构法对串列翼运动过程中内部计算域的网格进行变形和重构，外部计算域的网格保持不变．3 数值计算结果与分析3.1　不同平移幅值时摆动翼的推进性能针对不同的Ap，三种运动模式中摆动翼的水动力性、功率消耗、推进效率进行对比分析．Ap取0.21 Ccf，0.32 Ccf，0.42 Ccf，0.52 Ccf和0.62 Ccf，每个Ap对应的a1和a2值如表1所示．其他运动参数为：f =1 Hz，θm=25°，Uin=0.4 m/s．10.13245/j.hust.230402.T001表1摆动翼不同平移幅值Ap/Ccfa1/10-3a20.210.320.420.520.625.324.323.322.321.32-0.073-0.103-0.133-0.163-0.1933.1.1　水动力性能图4给出了不同Ap下三种运动模式中的摆动翼的Ct及C1．从图4(a)可以看出：三种运动模式中，Ct的波动幅值均是随着Ap的增加而增加．除了当10.13245/j.hust.230402.F004图4不同平移幅值下摆动翼的推力系数Ap=0.21Ccf时，模式2中Ct的波动幅值最大之外，当Ap取其他值时，均是模式3中Ct的波动幅值最大，模式1的最小．图4(b)中给出了三种运动模式中摆动翼的C1，并且分别给出了模式2和模式3的C1较模式1增加的百分比ζ2-1和ζ3-1．从图4(b)中可以看出：随着Ap的增加，三种运动模式的C1均是逐渐增加的，且增长速率不断加快．模式1和模式2的C1有着相似的变化趋势和增长速率，随着Ap的增加，两者的C1之间的差值表现出小幅度的增加，当Ap=0.21Ccf时，模式2的C1比模式1的大0.11，当Ap=0.62Ccf时，模式2的C1比模式1的大0.23，表明：随着Ap的增加，刚性前置翼的尾流场对摆动翼推进性能的提升变化不大．模式3中C1的增长速率明显高于另外两种运动模式，当Ap=0.21Ccf时，模式3的C1比模式1的大0.09，当Ap=0.62Ccf时，模式3的C1比模式1的大1.36，表明：随着Ap的增加，前置翼的柔性波动产生的尾流场对提升摆动翼的推进性能有着明显的促进作用．此外，从图4(b)中可以看出：随着Ap的增加，ζ2-1呈递减趋势，ζ3-1先增加后减小，表明相比于Ap的增加所带来的提升作用，刚性前置翼的尾流场对摆动翼推进性能提升的占比逐渐降低，前置翼的柔性波动运动对摆动翼推进性能提升的占比先增加后减少．3.1.2　功率消耗及推进效率图5给出了三种运动模式中摆动翼的平均输入功率系数C4和推进效率ηcf．从图5中可以看出C4和C1有着一致的变化趋势，摆动翼的C4均是随着Ap的增加而增加，并且增加的速率也是逐渐加大的．当Ap=0.21Ccf时，三种C4没有明显的差异，随着Ap的增加，模式3的C4逐渐扩大与另外两种运动模式之间的差值；由C4的变化趋势可知前置翼的存在能够增加摆动翼的功率消耗，且前置翼进行柔性波动运动时摆动翼消耗的能量显著增加．10.13245/j.hust.230402.F005图5摆动翼的平均输入功率系数和推进效率从图5可以看出：三种运动模式中摆动翼的ηcf均是先增加后减小，转折点为Ap=0.52Ccf．对于任意的Ap，模式1中摆动翼的ηcf始终最小，但值得注意的是：随着Ap的增加，模式1与模式2的差值逐渐减小，当Ap=0.21Ccf时有最大差值7.75%，当Ap=0.62Ccf时有最小差值0.80%．当Ap=0.21Ccf时，模式2的ηcf最大，当Ap＞0.21Ccf时，模式3的ηcf始终最大，并且随着Ap的增加，模式3与模式1之间的差值先增加后减小，当Ap=0.32Ccf时有最大差值6.61%，当Ap=0.62Ccf时有最小差值1.57%．此外，值得注意的是：模式2和模式3的ηcf相比于模式1的ηcf增加的百分比均是随着Ap的增加而减小，当Ap=0.21Ccf时分别为30.5%和21.5%，当Ap=0.62Ccf时分别为1.9%和3.7%，随着Ap的增加，前置翼的存在对摆动翼推进效率的增益效果逐渐降低．以上分析表明：前置翼尾流场的存在能够提升摆动翼的推进效率，并且前置翼进行柔性波动运动的尾流场能够进一步地提升摆动翼的推进效率；同时也表明摆动翼能够从前置翼的尾流场中吸收能量，从而提升自身的推进效率．3.2　不同运动频率时摆动翼的水动力性能针对不同的f，三种运动模式中摆动翼的水动力性、功率消耗、推进效率进行对比分析．f取0.5，0.75，1.0，1.25，1.5 Hz，其他运动参数为：Ap=0.42Ccf；θm=25°；Uin=0.4 m/s．3.2.1　水动力性能图6给出了不同f下三种运动模式中的摆动翼的推力系数Ct及C1．从图6(a)可以看出：Ct的波动幅值均是随着f的增加而增加；对于相同的f，模式3的Ct的波动幅值最大，模式1的波动幅值最小．从图6(b)中可以看出：三种运动模式的C1均是随着f的增加而增加，增加的速率也是逐渐增加，并且与3.1.1节中的变化特征一致，不再赘述．图6(b)给出了模式2和模式3的C1较模式1增加的百分比ζ2-1和ζ3-1，可以看出：ζ2-1和ζ3-1的变化趋势一致，均是随着f的增加而降低，并且始终ζ2-1小于ζ3-1，表明相比于f的增加所带来的提升作用，前置翼刚性不运动或者进行柔性波动运动的尾流场对摆动翼推进性能提升的占比均是逐渐降低．10.13245/j.hust.230402.F006图6不同运动频率下摆动翼的推力系数和平均推力系数3.2.2　功率消耗及推进效率图7给出了三种运动模式中摆动翼的C4和ηcf．从图7中可以看出：C4均是随着f的增加而增加，并且增加的速率也是逐渐加大的，变化特征与3.1.2节相同，不再赘述．从图7中可以看出：ηcf均是先增加后减小，转折点为f =0.75 Hz；对于任意给定的f，模式3的ηcf始终最大，模式1的ηcf始终最小，随着f的增加，模式2和模式3的ηcf相比于模式1增加的百分比均是逐渐地减小；当f =0.5 Hz时，模式3的ηcf比模式1的高出319.4%，模式2的ηcf比模式1的高出309.4%；当f =1.5 Hz时，模式3的ηcf比模式1的高出13.5%，模式2的ηcf比模式1的高出5.1%，表明前置翼的存在能够提升摆动翼的推进效率，并且在较小的f处，提升的效果更为明显，前置翼的柔性波动运动能够进一步地提升摆动翼的推进效率．对比模式2和模式3可以看出：两者有着相似的变化规律，且相差的值不大．10.13245/j.hust.230402.F007图7摆动翼的平均输入功率系数和推进效率3.3　摆动翼水动力性能增加的机理分析从前置翼上脱落于流场中的涡与摆动翼之间相互作用的机理，取运动参数为：a1=3.32×10-3；a2=-0.133；f =1 Hz；θm=25°；Uin=0.4 m/s．3.3.1　流场压力分布云图由于在一个周期内摆动翼运动的对称性，因此以1/4周期为例，给出了流场的压力分布云图，如图8所示，图中p为压力．当t=0.0 s时，摆动翼的摆角最大，摆动翼的左右表面分别被高压和低压所覆盖，压力差产生作用于摆动翼的水动力，从模式1到模式3，摆动翼两侧低压和高压的强度均不断增加，因此作用于摆动翼上的水动力不断增加，又由于摆角相同，因此C1不断增加(见图4和图6)．随着摆角的减小，摆动翼两侧的高压和低压的强度先增加后减小，当t=0.1 s时,达到最大，但是模式3中摆动翼两侧的压力幅值始终最大．摆动翼两侧的压力变化特征决定了摆动翼产生推力的性能，因此从模式1至模式3摆动翼的推力系数逐渐增加．10.13245/j.hust.230402.F008图8流场的压力分布云图3.3.2　流场涡结构特征及轴向速度分布图9给出了三种运动模式下的尾涡结构，图中红色漩涡的方向为逆时针，蓝色漩涡的方向为顺时针．从图9中可以看出三种游动模式的尾流场均为交错排列的反卡门涡街形式．每个旋漩涡在完全发展后呈现出逗号状，存在一个强度集中的涡核，相邻两个漩涡间的诱导作用产生一个指向摆动翼侧后方的射流，增加了摆动翼周围流体的轴向速度，由动量定理可知流体受到摆动翼的加速作用之后必然产生一个反作用于摆动翼的力，这个力便是摆动翼产生推力来源．因此，纵向上的射流强度越大，摆动翼产生的推力也必然越大．摆动翼尾流场中的纵向射流是由漩涡的诱导作用产生的，因此漩涡的产生和演变特征是影响摆动翼水动力性能的根本原因．10.13245/j.hust.230402.F009图9摆动翼的尾涡结构3.3.3　涡-翼相互作用的机制通过分析三种运动模式中的漩涡产生和演变特征，由涡-翼相互作用的机制探讨摆动翼水动力性能增加的机理．图10给出了摆动翼产生的涡结构的演变过程．对于模式1，当t=0.0 s时，摆角在正向最大值，摆动翼的尾缘处有一个尾缘涡TEVf-1，10.13245/j.hust.230402.F010图10摆动翼的尾涡结构1—LVf-1；2—TEVf-1；3—MV；4—TEVm-1；5—TEVm-2；6—LVf-2；7—LVf-1a；8—LVf-1b；9—TEVf-1a；10—TEVf-1b.TEVf-1与摆动翼下表面的附着涡连接在一起，在摆动翼的上表面附着一个与下表面旋转方向相反的涡LVf-1，此时摆动翼上下表面的附着涡存在的范围较小；随着摆角的减小，LVf-1不断发展壮大，并且TEVf-1在尾缘的带动下逐渐被拉长并发展壮大，当t=0.25 s时，摆角减小至零，TEVf-1从尾缘脱落；之后摆角在负向逐渐增加的同时向上进行平移运动，LVf-1不断地向尾缘移动并逐渐脱落于尾流场中，形成一个新的尾缘涡LVf-1，当t=0.5 s时,完全脱落至尾流场中，TEVf-1在脱落于尾流场中之后，逐渐聚集成一个强度高度集中的涡核．对于模式2，在摆动翼的上游存在着静止的前置翼，从前置翼的尾缘不断地脱落两层旋转方向相反的层涡MV．当t=0.0 s时，由于摆动翼的平移运动位于平衡位置处，摆角处于正向的最大值，因此MV位于摆动翼的下方，之后随着摆动翼向下的平移运动，MV被摆动翼截断，从前置翼上脱落的涡下泄于摆动翼的上表面，MV由上下两层旋转方向相反的带状涡组成，摆动翼上表面的附着涡LVf-1与MV中旋转方向的涡源源不断地融合在一起，相比于同时刻模式1中的LVf-1，模式2中的LVf-1强度更大；尾缘涡TEVf-1的演变过程与模式1的一致，不再赘述；当t=0.25 s时，摆动翼向下平移运动至极限位置，摆角减小至零，从t=0.25 s至t=0.5 s，摆动翼开始向上进行平移运动，并向负向进行摆动运动，MV自由脱落于尾流场中，LVf-1逐渐向尾缘运动并脱落于尾流场中．MV与LVf-1的这种融合行为使得模式2中摆动翼上表面附着涡的强度更大，并进一步导致摆动翼上表面流体的纵向速度大于模式1的(见图10(a)和(b)的对比)，模式2中摆动翼上表面的压力值更低，提升了摆动翼的推进性能．对于模式3，在摆动翼的上游存在着柔性波动运动的前置翼，并且前置翼的柔性波动运动模式为鱼类的BCF游动模式．不同于模式1和模式2，模式3中前置翼的尾缘交替的产生旋转方向相反的高强度漩涡，漩涡从前置翼的尾缘脱落，从前置翼与摆动翼的间隙穿过，流经摆动翼的表面脱落于尾流场中．当t=0.0 s时，前置翼的尾缘处存在一个明显的尾缘涡TEVm-1，TEVm-1的旋转方向与摆动翼上表面的附着涡LVf-1的旋转方向相反，在TEVm-1的挤压下，模式3中LVf-1的强度明显大于其他两种运动模式的；随着前置翼和摆动翼的运动，TEVm-1挤压着LVf-1向摆动翼的尾缘移动，并且LVf-1的强度和尺寸不断增加，当t=0.175 s时，LVf-1运动至摆动翼的尾缘，并与摆动翼的尾缘涡TEVf-1相遇；当t=0.25 s时，前置翼和摆动翼的平移运动至达单侧的最大位置处，TEVm-1从前置翼的尾缘脱落，在前置翼的尾缘处产生一个新的尾缘涡TEVm-2，LVf-1和TEVf-1相互挤压的位置发生变形．之后前置翼和尾鳍的平移运动方向改变，当t=0.35 s时，TEVm-1聚集成一个强度高度集中的圆形涡，LVf-1与TEVf-1分别分裂成LVf-1a和LVf-1b与TEVf-1a和TEVf-1b，其中LVf-1b，TEVf-1a和TEVf-1b脱落于尾流场中，LVf-1a形成一个新的尾缘涡，前置翼的尾缘涡TEVm-2通过前置翼与摆动翼的间隙运动至摆动翼的下表面，并对摆动翼下表面的附着涡LVf-2产生挤压作用．之后，随着前置翼的运动，TEVm-2不断地发展壮大，对摆动翼下表面附着涡LVf-2的挤压作用不断增强，导致LVf-2快速的增强．4 结论a．在摆动翼的上游存在一个较大弦长的前置翼，前置翼的刚性不运动或者柔性波动运动均能够提升摆动翼的推进性能，包括推力系数和推进效率．b．前置翼为刚性不运动或者柔性波动运动时对推力系数的增量均是随着平移幅值和运动频率的增加而增加，其中前置翼刚性不运动时增加的幅度较小，前置翼进行柔性波动运动时增加的幅度较大．c．在三种运动模式下，摆动翼的推进效率均是随着平移幅值和运动频率的增加而先增加后减小；其中前置翼进行柔性波动运动时，摆动翼的推进效率几乎始终最大，前置翼不存在时摆动翼的推进效率始终最小；随着平移幅值和运动频率的增加，前置翼的存在对摆动翼推进效率的增益效果均是逐渐的降低．d．当前置翼刚性不运动时，从其尾缘脱落两层旋转方向相反的涡，其中与摆动翼表面附着涡旋转方向相同的涡能够与摆动翼表面的附着涡相融合，提升附着涡的强度；当前置翼进行柔性运动时，其尾缘处脱落一个强度较大的尾缘涡，对摆动翼表面的附着涡产生挤压作用，提升附着涡的强度并加速附着涡的脱落．