无人水下航行器在海洋探索、环境监测和水下搜救等领域都存在巨大的价值和应用前景。随着海洋探索活动的逐渐深入，人们对新型水下航行器的性能要求也随之提高。传统水下航行器多采用螺旋桨作为推进装置，螺旋桨具有推进效率高的优势，但存在噪声大、低速条件下机动能力不足等问题。因此，仿生推进水下航行器凭借高效率、高机动、低噪声等优势成为新型水下航行器研发的一大热门方向。在仿生推进机理方面，早在1926年，Breder[1]就根据推进原理的不同，将鱼类推进方式分为中央鳍条/对鳍推进（median and/or paired fin, MPF）和身体/尾鳍（body and/or caudal fin, BCF）推进。多数鱼类采用BCF推进模式，该模式让鱼类具备较强的机动能力和较高的游动速度；采用MPF推进模式的鱼类则具备更高的推进效率和稳定性。1971年，Warren[2]最早发现海龟产生推进力的主要原因在于其前肢呈“8”字形轨迹的扑动。美国麻省理工学院Licht[3]通过观察海龟的游动过程发现，柔性前肢的存在增强了海龟的机动性和稳定性，并对海龟的转向速度进行了定量测量。Cubero[4]通过研究计算得出在理想状态下，柔性扑翼推进相比于螺旋桨推进可提升近两倍效率。在仿生推进结构方面，鲼科鱼是一类典型的通过扑动胸鳍产生前进推力的鱼类。它们有一对相对身体较大的扁平柔性胸鳍，在水下游动时具备出色的机动性和稳定性，可实现如向前推进、转弯、上升和潜水等多种不同的动作，许多受鲼科鱼启发的扑翼推进水下机器人也因此诞生。基于现有的生物学资料和研究人员对海洋中鲼科鱼的实地观测结果[5-6]可以发现，此类鱼的胸鳍弦向截面可通过NACA系列标准翼型来近似描述，其胸鳍扑动频率通常为0.5~2Hz，可在扑动频率较低的情况下产生较大的推进力。鲼科类机器鱼的胸鳍通常由柔性蒙皮和驱动器组合而成，根据柔性蒙皮的变形方式可分为被动变形胸鳍[7-8]和主动变形胸鳍[9-10]；根据驱动器的种类可分为非完全柔性胸鳍[7-10]和完全柔性胸鳍[11]。Lauder等[12]研制了一种利用胸鳍推进的机器鱼，发现在扑翼过程中个别鳍条产生了较为明显的形变，且胸鳍在弦向截面上存在一定程度的面积变化。Gao等[7]研制了一种模仿牛鼻鲼的机器鱼，对比了柔性胸鳍在无肋、肋展向布置、肋弦向布置时的推进能力，发现无肋的完全柔性胸鳍表现出最佳的推进性能。西北工业大学丁浩[13]系统分析了不同翼型厚度、翼型弯度、旋转轴位置下单自由度柔性扑翼的水动力性能，给出了具备最佳推进效果的柔性翼参数。在仿生推进水动力数值模拟方面，Blake[14]采用经典叶片单元理论将柔性翼分解成多个刚性小单元，进而分析柔性翼对应的流体动力特性。Frank等[15]利用浸入边界法求解了蝠鲼机器人的胸鳍流体动力学，发现大部分推力是由胸鳍在展向方向的远端部分产生的，且推进力随着时间产生周期性变化。Schultz等[16]利用水动力数值仿真方法模拟了鱼类的前游过程，得到了鱼前游过程对应的流场分布。Bozkurttas等[17]通过水动力分析模拟了不同展弦比胸鳍对运动性能的影响。结果表明，高展弦比胸鳍具有相对更大的推进力和更高的推进效率。值得注意的是，目前对于鲼科类机器鱼的流体动力学分析大多局限于只对胸鳍本身进行定性分析，很少有将胸鳍和水下航行器结合为一个整体进行仿真，有关胸鳍流体动力对水下航行器整体运动参数影响的研究较少。本文以融合软体机翼的水下航行器为研究对象，提出一种高可靠性、高稳定性的软体可扑动机翼设计方案，使水下航行器具备扑翼推进模式，从而增强其机动性。根据实际工作情况分别在二维和三维层面搭建软体机翼翼面波动数学模型。基于数学模型通过流体动力学仿真分析不同工况和扑动参数下的扑翼运动表现，得到扑翼过程中对应的流场涡量变化云图，阐述扑翼推进性能与流场内产生的后缘涡街之间的对应关系，总结扑翼产生流场后缘涡的脱落、排列及耗散规律，从水动力仿真的角度验证软体机翼扑动推进的可行性及软体扑翼水下航行器机动性的提升。1　软体扑翼推进系统设计软体扑翼推进系统由软体机翼单元、扑翼驱动单元、浮力补偿单元三部分组成，整体结构如图1所示。扑翼运动是由一个驱动电机带动一根鳍条做上下往复周期性摆动，软体机翼在主动鳍条和水阻尼的共同作用下产生被动变形，形成一个传播方向与前进方向相反的推进波，推动水下航行器向前运动，其波幅由前向后逐渐递减。为保证水下航行器能够正常工作，需调整水下航行器整体的净浮力使其在水面时呈近似中性。本文采用浮力补偿单元用于平衡软体机翼单元和扑翼驱动单元带来的负浮力。10.12050/are20220211.F001图1软体扑翼整体结构Fig.1Soft flapping-wing layout与传统水下滑翔机平板机翼相比，本文软体扑翼结构的优势在于具备更强的机动能力，同时兼具水下滑翔机低能耗滑翔运动模式，更能适应多样化的水下任务需求。与常见的仿蝠鲼柔性翼和多鳍条波动翼结构相比，本文软体扑翼结构只采用一个驱动电机作为动力源，扑翼运动时软体机翼在水阻尼和鳍条的共同作用下自然产生柔性变形。相比于模仿柔性翼多关节协同运动而使用多个驱动器或线驱动的机构，本文软体扑翼具有结构更简单可靠、工作更稳定和能耗更低等优势。在传统水下滑翔机基础上，用软体扑翼替换刚性平板机翼，即得到本文提出的软体扑翼水下航行器。软体扑翼共包括通信天线、前导流罩、三段密封舱、软体机翼、尾部导流罩和尾舵几个主要部分，其整体布局如图1所示。2　翼面波动数学模型搭建本文软体机翼尺寸及主要扑动参数见表1。软体扑翼结构在扑翼运动时，驱动电机带动鳍条做上下往复周期性摆动，可将鳍条运动近似看成相对中心截面对称的正弦运动，鳍条的运动学方程如式(1)所示10.12050/are20220211.T001表1软体机翼尺寸及主要扑动参数Table 1Dimensions and main flapping parameters of soft wing参数名称取值展弦比1.6鳍条长度(l0)/mm450机翼展长(D)/mm500最薄厚度/mm1.5横截面形状NACA0003扑动频率(f)/Hz0.5~1.5扑动振幅(θ)/(°)10~30波动参数o1.8转换系数ε2.0雷诺数1.91×105θ(t)=Asin(2πft)+φ (1)式中：θ(t)为t时刻驱动鳍条偏离中心位置的角度；A为扑翼运动的振幅；f为扑翼运动的频率；φ为初始时刻驱动鳍条偏离水下航行器中心位置的角度。当扑翼运动的对称中心与水下航行器中心位置重合时，其大小等于0。软体机翼在鳍条的带动下产生扑动，并与周围流体接触产生相互作用力，从而推动水下航行器向前移动。根据鳍条运动学方程搭建扑翼运动数学模型。为了便于给出扑翼运动数学模型，先将模型简化到二维层面。取软体机翼弦向截面带入运动学方程观察扑翼运动时软体机翼的形变效果。软体机翼的弦向截面可通过NACA0003标准翼形近似表示。基于鳍条运动学方程和实际扑翼工作原理。二维层面下，软体机翼弦向截面翼面波动数学模型可以近似表示为y2d(t)= (1- 1/(1+10t))Asin(2πox- 2πft) (2)式中：t为当前时间；y2(t)为t时刻弦向截面上各点的Y轴坐标；o为波动参数；x表示弦向截面上各点的X轴坐标，f为扑动频率；A为扑动振幅；(1- 1/(1+10t))项是为了让t=0时刻的初始形变由NACA0003标准翼型开始产生。将摆动幅值转化为角度，引入扑翼特征变量k等参数得到近似模型的简化表达式为y2d(t)=Δtl0tanθsin(kx- ωt) (3)式中：k=2πo；Δt = (1- 1/(1+10t))；ω=2π f；l0为鳍条的长度；θ为驱动鳍条偏离中心位置的角度。将f =1Hz、θ=20°带入式(3)，在扑翼运动达到稳定即Δt约等于1时，得到扑翼弦向截面在一个周期内的运动情况如图2所示。图中T表示扑翼运动周期。10.12050/are20220211.F002图2扑翼弦向截面一个周期内运动情况Fig.2Movement of the cross section of flapping-wing in chordwise in one period在二维扑动方程的基础上，引入Z轴坐标将其转化为三维扑动方程，机翼扑动是由于鳍条绕着X轴上下摆动而产生的。因此，机翼表面上各点在Z轴方向的运动幅度沿Z轴正向逐渐增加。在三维层面下，软体机翼扑动对应的翼面波动方程可以表示为y3d(t)=Δtl0tanθsin(kx-  ωt)(εz) (4)式中：ε为转换系数，其大小等于1/D，D为软体机翼展长；z为翼面上各点的Z轴坐标。将f =1Hz、θ=20°带入式(4)，在扑翼运动达到稳定即Δt约等于1时，得到水下航行器在一个周期内的扑翼运动情况如图3所示。10.12050/are20220211.F003图3扑翼水下航行器一个周期内运动情况Fig.3Movement of underwater vehicle with flapping-wings in one period3　水下扑翼运动流场显示与分析3.1　二维扑翼运动瞬态流场仿真本文水下扑翼运动瞬态流场仿真选用Realizable k- ε湍流模型求解黏性不可压缩雷诺方程，采用动网格方法控制扑翼边界运动。取扑翼的弦向截面作为二维仿真主体，求解水下扑翼运动对应的流场数值变化。扑翼弦向截面网格划分如图4所示，其网格规模为8.3万，网格类型为三角形网格，共包含速度入口、压力出口、长方形水域、前景域、背景域、扑翼壁面、重叠交界面等部分。在仿真过程中，通过Fluent软件的UDF（User Defined Function）接口定义扑翼运动模式，使扑翼壁面在运动方程的控制下产生柔性形变。机翼在形变过程中与周围流体接触产生相互作用力，整个前景域由于扑翼推进力的作用在长方形背景域内逐渐向前移动。最后，通过UDF运动方程迭代输出仿真过程中扑翼壁面在每个迭代步下的运动参数和推进力变化趋势。10.12050/are20220211.F004图4扑翼弦向截面网格划分Fig.4Grid layout of the cross-section of flapping-wing in chordwise本文采用网格无关性方法对仿真准确性进行验证。在工况和扑翼参数相同时，扑翼运动在三种网格规模下的平均推进力见表2。10.12050/are20220211.T002表2网格无关性验证Table 2Verification of grid independence网格规模/万平均推力/N算例a8.312.4算例b1012.1算例c1511.8由表2可得，当网格规模达到一定数量后，继续增加网格密度对本文扑翼平均推力影响较小。令扑翼参数f=1Hz、θ=20°，得到二维层面下软体扑翼运动在5s内产生的流场涡量变化云图如图5所示。10.12050/are20220211.F005图5扑翼运动在5s内产生的二维流场涡量云图Fig.5Vorticity contours of 2D flow field generated byflapping-wing movement in 5s由图5可得，软体机翼在发生变形时与周围流体相互接触，导致流体在机翼后缘脱落产生涡旋。在扑翼由静到动的初始阶段，流场内产生涡旋相较于扑翼运动达到稳定后更为明显。将初始时刻产生的涡旋称为启动涡，由于启动涡的存在，导致扑翼运动在初始时刻产生相较于稳定时更大的推进力。当扑翼运动达到稳定后，由图5观察到流场内的涡旋呈上下两列有序分布并缓慢沿着前进方向的反方向移动。观察涡旋的旋转方向可以发现，上排涡旋沿逆时针方向旋转；下排涡旋沿顺时针方向旋转，以此规律分布的涡旋正与卡门涡街[18]分布规律相反，即为诱导额外推力产生的逆卡门涡街。在涡旋脱落方面，在0~5s时间范围内由图5观察到扑翼后缘每隔1s脱落两个完整涡旋。扑翼运动的周期T为1s，即软体机翼每完成一次上扑或下扑动作（T/2），机翼尾部则脱落一个完整的后缘涡。这表明在每个T/2时间内，扑翼运动对应的流体涡旋先在机翼尾部聚集并逐渐增大；在T/2时刻，涡旋大小达到峰值并完成脱落，表明扑翼运动产生的推进力呈周期性变化趋势，且在每个T/2时刻达到一个推进力峰值。在涡旋耗散方面，以图5中1s时刻和4s时刻的涡旋分布情况为例，通过对比观察可以发现，图5中由右至左方向的第二个后缘涡在1s时刻整个涡旋的能量聚集在涡旋中心，在涡量云图中大部分呈现红色分布，而在4s时刻该涡旋的边缘在涡量云图中部分呈现黄色和绿色分布，表明随着时间的推移该涡旋的能量正逐渐向外耗散。进一步观察1s和5s时刻的涡旋分布情况可以发现，当扑翼运动在初始阶段时，机翼尾部刚刚脱落而下的后缘涡在涡量云图中大部分为红色。而当扑翼运动达到稳定即在5s时刻时，机翼尾部刚刚脱落而下的后缘涡在涡量云图中大部分为绿色，这再次表明当扑翼运动处于由静到动的初始阶段时，其产生的推进力相较于扑翼运动达到稳定时更大。令扑翼参数f=1Hz、θ=20°，分别取来流速度vi=0、0.3m/s、0.6m/s，得到扑翼运动过程中其推力系数和运动主体的运动学参数随时间变化趋势，如图6所示。由图6可得，扑翼运动产生的推力系数呈周期性变化趋势。当计算结果收敛后，扑翼在一个周期内的上扑和下扑阶段推进力变化趋势基本相同。在来流速度位于0~0.6m/s范围内，当来流速度较小时，扑翼运动在初始时刻产生的推进力更大，但来流速度对于稳定后的扑翼推进力影响较小，扑翼整体的位移和速度平均值在来流速度较小时取得较大值。10.12050/are20220211.F006图6不同来流速度对扑翼运动的影响Fig.6Influence of different inflow velocities on the motion of flapping-wing令来流速度vi=0.3m/s，扑动频率f =1Hz，分别取扑动振幅θ=10°、20°、30°，得到扑翼运动过程中其推力系数和运动主体的加速度、位移、速度变化趋势，如图7所示。由图7可得，在扑翼振幅位于10°~30°范围内，当扑翼振幅较大时，其产生的推力系数也较大，在5s的仿真时间内，扑翼整体的位移、速度平均值及加速度平均值随扑动振幅的增大而增大。10.12050/are20220211.F007图7不同扑动振幅对扑翼运动的影响Fig.7Influence of different amplitude on the motion of flapping-wing令来流速度vi=0.3m/s，扑动振幅 θ=20°，分别取扑动频率f =0.5Hz、1Hz、1.5Hz，得到扑翼运动过程中其推力系数和运动主体的加速度、位移、速度变化趋势如图8所示。由图8可得，在扑翼频率位于0.5 ~1.5Hz范围内，扑翼运动推力系数随频率的增加而增大，运动主体的位移、速度平均值及加速度平均值随频率的增加而增大。10.12050/are20220211.F008图8不同扑动频率对扑翼运动的影响Fig.8Influence of different frequency on the motion of flapping-wing综上所述，总体仿真结果表明，软体机翼在本文扑翼模式下可产生稳定的前进方向推进力，具备良好的运动性能。3.2　水下航行器扑翼运动瞬态流场仿真在二维流场仿真基础上进一步完成三维水下航行器扑翼运动瞬态流场仿真，其网格类型为Tetra/Mixed非结构网格，网格规模为153.3万，且通过网格无关性验证表明网格数量达到此规模后仿真已具备较高程度的准确性。令扑翼参数f =1Hz、θ=30°，仿真得到三维层面下水下航行器软体扑翼运动在水平面内一个运动周期产生的流场涡量变化情况，如图9所示。由图9进一步验证，在扑翼水下航行器由静到动的初始阶段，其表面的流场涡量最大，当扑翼运动逐渐达到稳定后，水下航行器表面的涡量减小，同时扑翼后方的流场产生移动方向与水下航行器前进方向相反的推进涡旋。具体地说，涡旋首先沿着扑翼边缘聚集，并随着扑翼运动逐渐向后均匀流去，在扑翼运动达到稳定后，推进涡旋呈半圆形排列包裹在扑翼后缘周围并随着软体机翼的扑动逐渐均匀向后移动。该仿真结果验证了本文扑翼结构可产生带动水下航行器向前移动的推进力，水下航行器在扑翼模式下具备较好的机动性。10.12050/are20220211.F009图9扑翼运动在水平面对应的流场涡量变化云图Fig.9Vorticity contours of flow field generated by underwater vehicle with flapping-wings in horizontal plane仿真得到扑翼水下航行器在纵垂面内一个运动周期产生的流场涡量变化情况如图10所示。由图10可得，扑翼运动在纵垂面内产生方向周期性变化的涡旋，当扑翼下扑时涡旋流向向上；当扑翼上扑时涡旋流向向下，两个方向的涡旋大小近似相同。该仿真结果表明，扑翼运动过程中水下航行器受到方向周期性变化的升力，其变化频率与扑翼扑动频率一致且在两个相反方向上的升力峰值基本相等，即在整数个扑翼运动周期内水下航行器受到的升力在0值上下呈对称性波动，波动平均值近似为0。该仿真结果验证了水下航行器扑翼模式在纵垂面内具备一定的稳定性。10.12050/are20220211.F010图10扑翼运动在纵垂面内的流场涡量变化云图Fig.10Vorticity contours of flow field generated by underwatervehicle with flapping-wings in vertical plane4　结论通过研究，可以得出以下结论：（1）本文提出一种软体可扑动机翼设计方案，主要包括软体机翼单元、扑翼驱动单元和浮力补偿单元。可将该扑翼结构集成于传统水下滑翔机得到本文软体扑翼水下航行器。近似搭建了软体扑翼运动数学模型。给出了软体扑翼水下航行器整体布局，其与传统水下滑翔机相比具有更高的机动性，同时保留了低能耗滑翔模式；与常见的仿生多关节扑翼和波动翼相比，其在结构上更加简单可靠、能耗更低。（2）本文分别完成了二维和三维层面下的扑翼运动瞬态流场仿真，得到了扑翼运动过程中对应的流场内涡量变化情况，阐述了流场后缘涡与扑翼推进特性之间的对应关系，总结了不同扑动参数和工况对扑翼推进性能的影响规律。仿真结果表明，扑翼产生的推力系数随扑动振幅、频率的增加而增大，而来流速度对稳定后的扑翼推力系数影响较小。仿真结果从水动力分析的角度验证了水下航行器软体机翼扑动推进的可行性，以及其机动性的提升。