迎角对扑翼运动气动性能影响的高性能数值研究

祝衔琦，孙祥程，王 娴，*

(1. 西安交通大学 航天航空学院，机械结构强度与振动国家重点实验室，西安 710049；2. 陕西省先进飞行器服役环境与控制重点实验室，西安 710049)

0 引 言

经过上亿年的进化，鸟类、昆虫、蝙蝠等生物均通过扑翼运动实现飞行。国内外学者受到扑翼生物的启发，研制了扑翼飞行器。扑翼飞行器主要在悬停、升力、能耗、控制等性能方面具有很大的优势。然而，目前扑翼飞行器的发展仍面临无法提供更高升力，飞行效率低，尺寸无法达到扑翼昆虫量级，无法做出更精确的飞行控制等问题[1]。因此，为了突破扑翼飞行器高升力、高效率、小尺寸等技术瓶颈，亟需开展更为细致的扑翼运动动力学研究，为设计更高效、高机动性的扑翼飞行器提供可靠的理论基础。扑翼飞行器在飞行过程中，相对来流方向与翼弦间的角度为迎角。迎角是影响扑翼飞行器飞行时升阻力系数的重要因素。当迎角超过一定范围，飞行器的飞行性能将会下降。因此，为了获得扑翼飞行器飞行时最佳的飞行姿态，开展不同迎角对扑翼运动的影响机理研究尤为重要。

目前关于迎角对扑翼运动的影响机理研究主要以实验为主：2007 年，宁琦等[2]对扑翼飞行器进行了风洞试验，研究了迎角对其气动特性的影响；2014 年，鲍锋等[3]基于色流实验与粒子图像测速（particle image velocimetry, PIV）技术，开展了迎角对单自由度扑翼运动升力的影响机理研究。此外，对扑翼运动其他运动参数的研究中，1989 年，Ennos[4]等基于高速摄影技术，对扑翼昆虫飞行时的运动学规律进行了探索；1999 年，Dickinson 等[5]通过模型实验的方法确定了昆虫飞行时的雷诺数范围为40 与400 之间。然而，实验方法普遍存在代价大、成本高等问题，且实验能够提供的数据有限，难以精细地获得能够直观、全面反映其空气动力学机制的流场细节特征。因此，高效、准确地对扑翼运动进行数值模拟研究已成为近些年的热点。国内外学者主要通过求解N-S 方程的方法对扑翼运动开展数值模拟研究。1998 年，Liu等[6]通过求解三维非定常N-S 方程的方法模拟并分析了昆虫飞行时的非定常流场结构；2000～2004 年，Wang 等利用二维计算[7-9]、三维实验[10]的方法研究了蜻蜓前飞时流场的变化与前后翼之间的影响，解释了蜻蜓飞行高升力机制；2002～2015 年，孙茂等[11-14]通过求解N-S 方程的方法，对扑翼昆虫运动时的稳定性高、能耗低等特点进行了研究，并对气动力机理进行了归纳；2008 年，Young 等[15]通过求解三维不可压N-S 方程的方法模拟了昆虫扑翼运动，研究了振幅、频率等参数对扑翼运动性能的影响。

扑翼的运动方式复杂、外形复杂，可靠的模拟结果依赖于复杂边界条件易实施的数值方法与高分辨率网格。近年来，高效准确的格子Boltzmann 方法（lattice Boltzmann method, LBM）[16]在计算流体力学（computational fluid dynamics, CFD）中得到普遍应用，由于其算法简单、边界条件易实施、质量动量严格守恒、适于并行等优点，被应用于各个领域，在扑翼运动的数值模拟中也得到快速发展。2014 年，Keisuke 等[17]采用浸没边界格子Boltzmann 方法针对类蜻蜓扑翼模型前后翼相位滞后角开展了数值模拟研究，结果表明，扑翼运动的移动方向取决于相位滞后角。2020年，彭连松等[18]采用基于LBM 的商业软件Xflow，对串列双翅扑翼和单对翅扑翼进行数值模拟，结果表明，双翅悬停效率高于单翅拍动的悬停效率。此外，对扑翼运动进行数值模拟时，对运动复杂边界的精确描述是十分重要的。动边界识别方法总体上分为拉格朗日法和欧拉法，基于拉格朗日法的动边界识别法，其重构网格计算时间长，计算效率低。而Level-Set 方法[19]基于欧拉法，即使在规则分布的欧拉网格中也能精确捕捉到边界点的位置。同时，由于网格的规则性，更适合并行计算，使得重构网格的时间大大缩短，若在较高的网格分辨率下，可精准描述复杂物体边界运动。如2005 年，李会雄等[20]基于Level-Set方法捕捉了气-液两相流运动相界面; 2006 年，陈凡红等[21]基于Level-Set 方法追踪了各种外形的物体在旋转流场中的变形还原效果，研究表明，Level-Set 方法可以准确追踪运动界面，且容易实现，具有较大通用性。另一方面，在硬件上，近些年来，图形处理器（graphics processing unit，GPU）发展速度远超中央处理器（central processing unit, CPU），GPU 开始逐渐应用于非显示的通用图形处理器（general purpose graphics processing unit, GPGPU），出现了较多基于CPU/GPU 异构系统的CFD 高性能计算在各类计算流体力学方法（如N-S、LBM 等）上的应用研究，使得CFD 计算效率更高。其中GPU 的单指令多线程模式与LBM 算法完美匹配，两者结合，较之CPU 可达百倍以上的加速比[22]。

基于以上种种，本文针对扑翼运动问题，采用LB 方法求解流场，Level-Set 方法追踪复杂运动边界，并结合GPU 并行加速技术，开发了扑翼运动三维高效数值模拟程序，对扑翼在不同迎角下，上下扑动时其气动性能开展数值研究。通过精细捕捉其流体力学特征量—涡系结构、压力场等，分析其升阻力产生的机制，获得其最佳运动迎角。同时，本文也为扑翼运动的研究提供了高效、可靠的数值方法。

1 计算模型与数值方法

1.1 计算模型

本文根据Liang[12]选用的扑翼模型进行建模，其几何模型如图1 所示。

图1 扑翼几何模型Fig. 1 Prototype of flapping-wing

翼面的展长L=4.90 cm ，平均弦长D=1.24 cm，厚度s=0.14 cm， θ为迎角。两个扑翼分别绕各自的翼根作上下扑动运动，扑动方程为：

其中，α为 扑动角，αmax为 扑动角 振 幅，f为扑动频率，φ为扑动初始相位角。双翼同步扑动，扑动角振幅αmax=30°， 扑动频率f=32 Hz ，扑动初始相位 φ=0°。扑动角速度用 ω表示。

雷诺数Re=UrefD/υ =100， 无量纲D=18，其中Uref=4αmaxfr，Uref为扑翼运动参考速度，r为扑翼的旋转半径。本文分别计算了 θ=0°～60°时不同迎角的流场情况。流体的计算区域为Lx×Ly×Lz，每个方向取1 6D，即 288×288×288，总计2 389 万网格。扑翼的升力系数及阻力系数定义为：

其中，Fz为 扑翼受到的沿z方向升力，Fx为扑翼受到的沿x方向阻力，A为扑翼的面积。

1.2 数值方法

1.2.1 格子Boltzmann 方法

本文采用LBM-D3Q19 模型[16]（图2）求解流场，其离散方程如下：

图2 D3Q19 模型Fig. 2 D3Q19 model

对于流场，宏观边界条件为：进口u=U∞,v=w=0；出口以及其余边界均采用充分发展边界条件。分布函数f边界条件为：所有边界均采用非平衡态外推格式[23]，如下：

1.2.2 LBM 移动固壁边界条件

1.2.2.1 曲面边界非平衡态外推格式

对于移动的固体壁面，本文采用非平衡态外推格式[23]确定固体边界。如图3 所示，非平衡态外推格式将固体节点rw的分布函数fi(rw,t)分解为平衡态与非平衡态两部分。

图3 曲线边界非平衡态外推格式Fig. 3 Schematic of extrapolation method for curved boundary conditions in lattice Boltzmann method

非平衡态部分中，q表示与壁面相交一个网格间距中处于流体区域的比例：

1.2.2.2 新生流体节点的重新填充

物体在流场中运动时，当前时刻是固体中的网格节点下一时刻可能会成为流体节点，此时需要重新填充新生流体节点[24]。如图4 所示，蓝色虚线为t时刻固体壁面，t+δt时刻，固体壁面运动至蓝色实线处，p点由固体节点变为流体节点。此时，新生流体节点p需要进行分布函数的计算。

图4 新生流体节点的重新填充Fig. 4 Schematic of initialization of new fluid nodes

公式如下所示:

1.2.3 Level-Set 动边界识别方法

Level-Set 算法是一种用隐函数描述边界运动的方法[19]，利用等值面函数 φ(x) ， 让 φ以一定速度运动，零等值面即为固体边界。在欧拉坐标系下， φ(x)可以用符号距离函数表示。符号距离函数表示当前节点与物体边界之间的最短距离，它在边界附近充分光滑，可以用来表示Level-Set 函数。符号距离函数d(x)的定义如下：

针对流/固动边界的实时更新问题，采用Level-Set 动边界识别方法对标准格子立方体网格的边界进行实时更新。如图5 所示，图5（a）为物体运动前状态，每个节点存储 φ(x)。 将 φ(x)跟随物体作平移、旋转运动后，如图5（b）所示，并将运动后的 φ(x)插值至欧拉背景网格中，更新每个节点的 φ(x)，物体随即发生运动，从而识别出动边界。

图5 Level-Set 动边界识别方法Fig. 5 Schematic of identification method of moving boundary of Level-Set

1.2.4 GPU 程序移植

GPU 并行模式是将数据分配于多条线程（thread），多条线程组成一个块（block），GPU 中多个运算单元同时计算每个block 中的threads，而每条thread 中的数据则串行计算。

本文自主搭建了基于LBM-GPU 与Level-Set 动边界识别方法的数值计算平台，计算流程图如图6 所示。在流场计算的每一个迭代中，经过分布函数的碰撞与迁移之后，对下一迭代步的边界位置采用Level-Set 方法进行识别，并映射至LBM 网格中，对网格节点属性，即固体内部、固体壁面及流体，进行更新，并对新生流体节点的分布函数进行填充，完成一次循环。

图6 计算流程图Fig. 6 Flow chart of numerical method

此外，本计算采用一颗NVIDIA Tesla K20M GPU加速，以2 389 万个网格计算工况为例，计算7 200 个LBM 步长，总耗时4.20 h；相同情况下使用一颗Intel Core i7-3770 CPU 计算共耗时298.2 h，加速比为71。

1.3 程序验证

1.3.1 正确性验证

本文利用圆柱沉降问题进行程序验证，并与Li[25]的模拟结果进行对比。文献采用应力积分法求得物体受力，从而获得物体运动的轨迹，实现追踪动边界的效果。如图7 所示，圆柱在一个充满流体的管道中自由下落，由于自身的重力以及流体的作用，圆柱会边旋转边平移地下沉。管道宽度L=0.4 cm，圆柱直径d=0.1 cm， 圆柱密度为 ρc=1.003 g/cm3。流体密度 ρ=1 g/cm3，运动黏度 υ =0.01 cm2/s，重力加速度G=980 cm/s2，圆柱在靠近左侧壁面x=0.076 cm位置开始释放。

图7 圆柱沉降示意图Fig. 7 Schematic of sedimentation of a circular cylinder

图8 为圆柱下沉过程中圆柱下落轨迹、水平方向速度、竖直方向速度和角速度随时间变化的曲线图。可以看出，本文结果与文献对比一致。此外，从图8（a）可以看出，圆柱运动最终达到了以恒定速度up竖直下落的状态。稳态雷诺数定义为Re=dup/υ，计算得到该工况下的稳态雷诺数为Re=1.03，文献计算值为Re=1.03，两者一致，本程序的正确性得以验证。

图8 数值计算与参考结果验证[25]Fig. 8 Validation of numerical with reference case[25]

1.3.2 网格无关性验证

本文使用不同三种网格分辨率模拟了图1 所示的扑翼运动。计算参数θ =0°、Re=100、St=2Az f/U∞=0.15， 其中Az为拍动幅值。三个不同的网格分辨率分别为200×200×200 （Case A）、245×245×245 （Case B）、288×288×288 （Case C），三种网格分辨率对应的扑翼所占的网格点数分别为4 513、5 528 和6 498。

图9 为三种网格分辨率下的翼尖处涡量云图。从图中可以看出，涡量云图基本一致。但随着网格分辨率的增加，在扑翼尾迹处的旋涡分离趋势更加明显，涡的细节也更清晰。

图9 翼尖涡量云图Fig. 9 Vorticities contour at wing-tip

图10 为不同网格分辨率下的升力系数与阻力系数随时间变化曲线图，其中t为无量纲时间：当前时间步数与一个周期时间步长的比值，扑翼上下扑动一次为一个周期。可以看出，不同的网格分辨率获得的结果基本一致。

图10 网格无关结果图Fig. 10 Grid independence results

本文计算6 个扑动周期，一个扑动周期对应1 200个LBM 步长，共7 200 个LBM 时间步长，对于A、B、C 三种工况，计算时间分别为3.68 h、3.86 h、4.20 h。三种工况的计算均采用GPU 加速，因计算时间较短，本文在之后的计算中均采取288×288×288 的网格分辨率，以捕捉到更精细的流场细节。

2 结果与讨论

2.1 迎角对升阻力系数的影响

如图11 所示，以 θ=0°为例，在一个周期内，扑翼在t=0 时 向下拍动，此时扑动角速度 ω最 大；t=T/4时转为向上拍动，此时 ω=0；t=T/2时向上拍动，此时ω最 大；t=3T/4 时转为向下拍动，此时ω =0。

图11 一个周期内扑翼运动示意图Fig. 11 Schematic of flapping-wing movement in a period

图12 为不同迎角下升力系数CL、阻力系数CD随t变化的曲线。曲线取自扑翼运动的第5、6 个周期，此时流场已经处于稳定状态，各个周期内升阻力系数变化稳定一致。如图12（a）所示，在一个周期内随着时间推进，CL在扑翼下扑过程中达到最大值，在扑翼上扑过程中达到最小值。随着 θ 的 增大，CL总体呈上升 趋 势，直 到 θ=60°时 ，CL开 始 下 降。当 θ=40°与60°时 ，CL最小值接近于0。如图12（b）所示，当θ=0°时，扑翼的迎风面积极小，CD接近于0，且CD变化较小。当迎角逐渐增大时，扑翼的迎风面积逐渐增大，CD呈上升趋势，并在扑翼下扑过程中达到最大值，在扑翼上扑过程中达到最小值。

图12 升阻力系数随时间变化曲线Fig. 12 Variation curve of lift and drag coefficient with time

图13 分别为平均升力系数C¯L、平均阻力系数C¯D以及升阻比C¯L/C¯D随迎角变化曲线图，此处，平均升阻力系数C¯L、C¯D为一个周期内的算术平均值。从图13（a）中可以看出，C¯L随迎角的增大呈先增大后减小 的 趋 势，当 θ ≥40°时 ，C¯L开 始 下 降。从 图13（b）中可以看出，C¯D随着迎角的增大逐渐增大。从图13（c）中可以看出，C¯L/C¯D随迎角增大呈先增大后减小的趋势，当 θ=10°时 ，C¯L/C¯D最大，此时扑翼的气动性能达到最佳。当 θ=60°时 ，C¯L/C¯D较 θ=0°时更小，该迎角下扑翼运动气动性能反而下降。

图13 平均升阻力系数随迎角变化曲线Fig. 13 Variation curve of average lift and drag coefficient with angle of attack

下文通过压力场和流场结构对扑翼运动的气动特性机理进行分析。

2.2 气动特性机理分析

2.2.1 压力场分析（图1（c））受力面减小，反而使CL变小。

图14 t =0上下翼面压力云图Fig. 14 Pressure contour of lower and upper surface of wing ont=0

图15 t =T/2上下翼面压力云图Fig. 15 Pressure contour of lower and upper surface of wing ont=T/2

同样地，分析图16（b、d），从翼根到翼尖方向，Δp绝 对值也逐渐增大，但是随着迎角的增大， Δp绝对值明显比下扑阶段小，定量说明了图12（a）中高迎角时CL最小值接近于0 的情况。两幅正压曲线在翼尖处下降趋势最明显，说明在翼尖处压力变化最大。

图16 压力系数和压力差沿翼型弦长方向变化曲线图（t =0 ，t =T/2）Fig. 16 Pressure coefficient and pressure difference distribution along chord length direction of airfoil on t =0 and t=T/2

2.2.2 流场及涡结构分析

图17 为不同迎角下扑翼运动涡系结构随时间的演变，图中Ω表示LBM 格子单位下的涡量。从图中可以观察到前缘涡，扑翼拍动过程中，上翼面前缘开始产生前缘涡。前缘涡具有附着于扑翼表面不脱落的特性，附着于翼面的面积越大，提供的升力越大，而一旦当前缘涡从扑翼上脱落，升力系数将会大幅下降。可以看出，当t=0与t=T/4时，在翼尖处前缘涡最大，此时，翼尖处获得的升力最大。此外，在同一迎角下，t=T/4 ～ 3T/4时，附着在翼面的前缘涡不断从翼面前缘脱落，此时CL小于该迎角下整个周期的C¯L；t=3T/4 至 下一个周期的t=T/4，附着在翼面的前缘涡开始生成，此时CL大于该迎角下整个周期的C¯L，这解释了图12（a）所示的升力变化趋势。

图17 不同迎角下扑翼运动涡系结构随时间的演变Fig. 17 Variation of vortices structure of flapping-wing for different angle of attack with time

图18 为不同迎角下，t=0时翼尖处的流线与压力图。从图18 中可以看出，在拍动的过程中，扑翼下方压强较大，上方压强较小，空气在翼尖处上下压差更大，从扑翼下表面经由翼尖绕至上表面，形成翼尖涡。随着迎角的增大，翼尖涡逐渐由翼面处攀升至翼面上方，强度与影响范围逐渐增大，从而导致涡致阻力的增加，使得CD逐渐增大。

图18 t =0不同迎角流线与压力图Fig. 18 Streamline and pressure contour for different angle of attack on t=0

如图17 所示，在t=0时，即扑翼处于下拍阶段，翼尖涡从翼尖处分离并顺着来流速度方向传输，导致t=0 ～T/2时 涡致阻力下降，此时扑翼CD下降；在t=T/2时，即扑翼处于上拍阶段，翼尖涡又开始快速生成，使得t=T/2 ～T时涡致阻力上升，此时扑翼CD上升。

如图17 所示，随着迎角的增大，前缘涡附着于翼面的面积明显增大，且与翼尖涡逐渐融合，并在扑翼尾迹处逐渐形成了封闭的涡环，不易耗散，且涡环与扑翼间的相互作用，使得涡环中的能量得以重复利用，使扑翼维持在一个高升力状态。 θ=60°时，扑翼尾迹处的旋涡开始变得不稳定并破碎形成一个个细小的涡，流场结构紊乱，旋涡的能量被分解耗散，此时，扑翼的气动性能下降。

3 结 论

本文结合LBM 与Level-set 动边界识别方法，采用GPU 加速，建立了三维含有运动边界流场的高效数值模拟方法并开发了相应程序。基于此，开展了不同迎角下扑翼运动气动特性机理分析。研究结果表明：

1） 随着迎角的增加，扑翼的平均升力系数先增大后减小，平均阻力系数逐渐增大。升阻比在10°时达到最大。此时，扑翼能够获得最佳气动性能。

2） θ =0°～40°范围内，随着迎角的增加，扑翼上下翼面正负压区面积增大，上下翼面压差增大，扑翼获得的升力增加；翼尖处压差最大，给扑翼提供了主要的升力。

3） 随着迎角的增大，前缘涡附着于翼面的面积明显增大，扑翼后方脱落的涡旋也较难耗散，提高了扑翼的升力；同时，翼尖涡的强度与影响范围逐渐增大，导致涡致阻力增加，阻力系数增加。

LBM 结合Level-Set 方法可精确捕捉具有复杂外形的运动边界，同时，采用GPU 加速能够获得极高的加速比。本文的数值方法可为扑翼运动提供高效、可靠的研究手段。