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所示。

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

$ \alpha ={\alpha }_{{\rm{max}}}{{\rm{sin}}}\left(2{\text{π}} f+\varphi \right) $

(1)

其中， $ \alpha $ 为扑动角， ${\alpha }_{{\rm{max}}}$ 为扑动角振幅， $ f $ 为扑动频率， $ \varphi $ 为扑动初始相位角。双翼同步扑动，扑动角振幅 ${\alpha }_{{\rm{max}}}=30^{ \circ }$ ，扑动频率 $ f=32\;{\rm{Hz}} $ ，扑动初始相位 $\varphi =0 ^{ \circ }$ 。扑动角速度用 $ \omega $ 表示。

雷诺数 $Re={U}_{{\rm{ref}}}D/\upsilon \;=100$ ，无量纲 $ D=18 $ ，其中 ${U}_{{\rm{ref}}}=4{\alpha }_{{\rm{max}}}fr$ ， ${U}_{{\rm{ref}}}$ 为扑翼运动参考速度， $ r $ 为扑翼的旋转半径。本文分别计算了 $ \theta =0^{ \circ }\sim 60^{ \circ } $ 时不同迎角的流场情况。流体的计算区域为 $ {L}_{x}\times {L}_{y}\times {L}_{z} $ ，每个方向取 $ 16D $ ，即 $ 288\times 288\times 288 $ ，总计2389万网格。扑翼的升力系数及阻力系数定义为：

$ {C}_{L}=\frac{{F}_{z}}{0.5\rho {U}_{{\rm{ref}}}^{2}A},\;{C}_{D}=\frac{{F}_{x}}{0.5\rho {U}_{{\rm{ref}}}^{2}A} $

(2)

其中， $ {F}_{z} $ 为扑翼受到的沿z方向升力， $ {F}_{x} $ 为扑翼受到的沿x方向阻力， $ A $ 为扑翼的面积。

1.2 数值方法 1.2.1 格子Boltzmann方法

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

$ {f}_{i}\left({\boldsymbol{r}}+{{\boldsymbol{e}}}_{i}{\rm{\delta }}_{t},t+{\rm{\delta }}_{t}\right)-{f}_{i}\left({\boldsymbol{r}},t\right)=\frac{1}{\tau }\left[{f}_{i}^{{\rm{eq}}}\left({\boldsymbol{r}},t\right)-{f}_{i}\left({\boldsymbol{r}},t\right)\right] $

(3)

式中： $ {f}_{i} $ 与 ${{f}_{i}}^{{\rm{eq}}}$ 分别为密度分布函数与其对应的平衡态分布函数， $ {{\boldsymbol{e}}}_{i} $ 为粒子离散速度， $ {\rm{\delta }}_{t} $ 为时间步长， $ \tau $ 为无量纲松弛时间。

$ {f}_{i}^{{\rm{eq}}}=\rho {w}_{i}\left[1+3\left({{\boldsymbol{e}}}_{i} \cdot {\boldsymbol{u}}\right)+\frac{9}{2}{\left({{\boldsymbol{e}}}_{i} \cdot {\boldsymbol{u}}\right)}^{2}-\frac{3}{2}{{\boldsymbol{u}}}^{2}\right] $

(4)

其中，每个方向上的分布函数的权系数 $ {w}_{i} $ 为：

${w_i} = \left\{ {\begin{array}{*{20}{l}} {1/3,}&{i = 0}\\ {1/18,}&{i = 1,2, \cdots ,6}\\ {1/36,}&{i = 7,8, \cdots ,18} \end{array}} \right. $

(5)

离散速度模型 $ {{\boldsymbol{e}}}_{i} $ 为：

$ {{\boldsymbol{e}}}_{i}=\left\{ {\begin{array}{*{20}{l}} {\left( {0,0,0} \right)}，\quad i = 0\\ {\left( { \pm 1,0,0} \right),\left( {0, \pm 1,0} \right),\left( {0,0, \pm 1} \right)}，\quad i = 1,2, \cdots ,6\\ {\left( { \pm 1, \pm 1,0} \right),\left( {0, \pm 1, \pm 1} \right),\left( { \pm 1,0, \pm 1} \right)}，\quad i = 7,8, \cdots ,18 \end{array}} \right. $

(6)

宏观密度 $ \;\rho $ 和速度u分别为：

$ \rho ={\sum }_{i}{f}_{i} $

(7)

$ {\boldsymbol{u}}={\sum }_{i}{{\boldsymbol{e}}}_{i}{f}_{i}/\rho $

(8)

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

$ {f}_{i}\left({{\boldsymbol{r}}}_{s},t\right)={f}_{i}^{{\rm{eq}}}\left({{\boldsymbol{r}}}_{s},t\right)+\left[{f}_{i}\left({{\boldsymbol{r}}}_{f},t\right)-{f}_{i}^{{\rm{eq}}}\left({{\boldsymbol{r}}}_{f},t\right)\right] $

(9)

其中， ${f}_{i}^{{\rm{eq}}}({{\boldsymbol{r}}}_{s},t)$ 为边界节点 $ {{\boldsymbol{r}}}_{s} $ 的平衡态分布函数， ${f}_{i}({{\boldsymbol{r}}}_{f},t)、{f}_{i}^{{\rm{eq}}}({{\boldsymbol{r}}}_{f},t)$ 分别为边界节点临点 $ {{\boldsymbol{r}}}_{f} $ 的分布函数与平衡态分布函数。

1.2.2 LBM移动固壁边界条件 1.2.2.1 曲面边界非平衡态外推格式

对于移动的固体壁面，本文采用非平衡态外推格式^[23]确定固体边界。如图3所示，非平衡态外推格式将固体节点 $ {{\boldsymbol{r}}}_{w} $ 的分布函数 $ {f}_{i}\left({{\boldsymbol{r}}}_{w},t\right) $ 分解为平衡态与非平衡态两部分。

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

$ q=\frac{\left|{{\boldsymbol{r}}}_{f}-{{\boldsymbol{r}}}_{b}\right|}{\left|{{\boldsymbol{r}}}_{f}-{{\boldsymbol{r}}}_{w}\right|} $

(10)

其中， $ {{\boldsymbol{r}}}_{b} $ 为固体壁面节点， $ {{\boldsymbol{r}}}_{f} $ 、 $ {{\boldsymbol{r}}}_{ff} $ 为流体节点， $ {{\boldsymbol{r}}}_{w} $ 为固体节点。固体节点 $ {{\boldsymbol{r}}}_{w} $ 的宏观密度直接由流体节点 $ {{\boldsymbol{r}}}_{f} $ 处的密度 $ \;{\rho ({\boldsymbol{r}}}_{f}) $ 代替，速度 ${{{{\bar {\boldsymbol{u}}}}}_{w}}$ 为：

$ {{{{\bar {\boldsymbol{u}}}}_w}} = \left\{ {\begin{array}{*{20}{l}} {q{{\boldsymbol{u}}_{w1}} + \left( {1 - q} \right){{\boldsymbol{u}}_{w2}},}&{q < 0.75}\\ {{{\boldsymbol{u}}_{w1}},}&{q \geqslant 0.75} \end{array}} \right.$

(11)

其中， ${{\boldsymbol{u}}}_{w1}\equiv \dfrac{{{\boldsymbol{u}}}_{b}+\left(q-1\right){{\boldsymbol{u}}}_{f}}{q}$ ， ${{\boldsymbol{u}}}_{w2}\equiv \dfrac{2{{\boldsymbol{u}}}_{b}+\left(q-1\right){{\boldsymbol{u}}}_{ff}}{1+q}$ ， $ {{\boldsymbol{u}}}_{b} $ 为扑翼给定的在b点的壁面速度，由此可计算出 $ {{\boldsymbol{r}}}_{w} $ 的平衡态分布函数 ${f}_{i}^{{\rm{eq}}}\left({{\boldsymbol{r}}}_{w},t\right)$ 。非平衡态部分也可类似地给出：

$f_i^{{\rm{neq}}}( {{{\boldsymbol{r}}_w},t} ) = \left\{ {\begin{array}{*{20}{l}} {qf_i^{{\rm{neq}}}( {{{\boldsymbol{r}}_f},t} ) + ( {1 - q} )f_i^{{\rm{neq}}}( {{{\boldsymbol{r}}_{ff}},t} ),}&{q < 0.75}\\ {f_i^{{\rm{neq}}}( {{{\boldsymbol{r}}_f},t} ),}&{q \geqslant 0.75} \end{array}} \right.$

(12)

综上， ${{\boldsymbol{r}}_w}$ 沿 ${{\boldsymbol{c}}_i}$ 方向碰撞后的分布函数为：

$ {f}_{i}\left({{\boldsymbol{r}}}_{w},t\right)={f}_{i}^{{\rm{eq}}}\left({{\boldsymbol{r}}}_{w},t\right)+\left(1-\frac{1}{\tau }\right){f}_{i}^{{\rm{neq}}}\left({{\boldsymbol{r}}}_{w},t\right) $

(13)

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

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

公式如下所示:

$ \begin{split} {f}_{p}\left({{\boldsymbol{r}}}_{p},t+{\rm{\delta }}_{t}\right)=&{f}_{p}^{{\rm{eq}}}\left({{\boldsymbol{r}}}_{p},t\right)+{f}_{p}^{{\rm{neq}}}\left({{\boldsymbol{r}}}_{p},t\right)\\ \;\;\;\;\;\;\;\;\;\; =&{f}_{p}^{{\rm{eq}}}(\overline {\rho },{{\boldsymbol{u}}}_{b})+\left[{f}_{g}({{\boldsymbol{r}}}_{g},t)-{f}_{g}^{{\rm{eq}}}({{\boldsymbol{r}}}_{g},t)\right] \end{split}$

(14)

其中：在 ${f}_{p}^{{\rm{eq}}}$ 的计算中，速度近似取壁面b点速度 $ {{\boldsymbol{u}}}_{b} $ ，密度取平均密度 $\; \bar {\rho }$ ； ${f}_{p}^{{\rm{neq}}}$ 近似取临点相应的值，由点g的分布函数 ${f}_{g}^{{\rm{n}}}$ 与非平衡态分布函数 ${f}_{g}^{{\rm{neq}}}$ 的差获得；点g的位置通过外推得到，外推方向 $ {c}_{n} $ 由 ${\boldsymbol{n}}_{{\rm{wall}}} \cdot \dfrac{{\boldsymbol{c}}_{\boldsymbol{i}}}{{|{\boldsymbol{c}}}_{\boldsymbol{i}}|}$ 的最大值确定。

1.2.3 Level-Set动边界识别方法

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

$ d\left({\boldsymbol{x}}\right)={{\rm{min}}}\left(\left|{\boldsymbol{x}}-{{\boldsymbol{x}}}_{I}\right|\right){}{{\boldsymbol{x}}}_{I}\in \varOmega $

(15)

其中， $ {\boldsymbol{x}} $ 为空间节点坐标， $ {{\boldsymbol{x}}}_{\boldsymbol{I}} $ 为边界点坐标， $ \varOmega $ 为边界。

符号距离函数表示的 $ \phi \left({\boldsymbol{x}}\right) $ 定义为：

$ \phi \left({\boldsymbol{x}}\right)=\left\{ {\begin{array}{*{20}{l}}d\left({\boldsymbol{x}}\right)&,{\boldsymbol{x}}\in {\varOmega }^+\\ 0&, {\boldsymbol{x}}\in \varOmega \\ -d\left({\boldsymbol{x}}\right)&, {\boldsymbol{x}}\in {\varOmega }^-\end{array}} \right. $

(16)

其中， $ \varOmega $ 为边界， $ {\varOmega }^{+} $ 为边界外部，即流体部分， $ {\varOmega }^{-} $ 为边界内部，即固体部分。

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

1.2.4 GPU程序移植

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

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

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

1.3 程序验证 1.3.1 正确性验证

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

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

1.3.2 网格无关性验证

本文使用不同三种网格分辨率模拟了图1所示的扑翼运动。计算参数 $\theta =0^{ \circ }$ 、 $ Re=100 $ 、 $S t={2{A}_{z}f}/{{U}_{\infty }}= 0.15$ ，其中 $ {A}_{z} $ 为拍动幅值。三个不同的网格分辨率分别为200×200×200 （Case A）、245×245×245 （Case B）、288×288×288 （Case C），三种网格分辨率对应的扑翼所占的网格点数分别为4513、5528和6498。

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

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

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

2 结果与讨论 2.1 迎角对升阻力系数的影响

如图11所示，以 $ \theta =0^{ \circ } $ 为例，在一个周期内，扑翼在 $ t=0 $ 时向下拍动，此时扑动角速度 $ \omega $ 最大； $t=T/4$ 时转为向上拍动，此时 $ \omega =0 $ ； $t=T/2$ 时向上拍动，此时 $ \omega $ 最大； $t=3T/4$ 时转为向下拍动，此时 $ \omega =0 $ 。

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

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

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

2.2 气动特性机理分析 2.2.1 压力场分析

图14为不同迎角下， $t=0$ 时上下翼面的压力云图，图中p表示LBM格子单位下的压力。从各个压力场情况可以看出，下翼面压力大于上翼面，此时 $ {C}_{L} $ 达到最大值。在 $t=0$ 时，随着迎角的增大，上下翼面的负压、正压区域面积逐渐增大，并从翼尖逐步扩展到翼中，从翼前缘扩展到后缘，导致此时上下翼面压差增大，尤其当 $ \theta =0^{ \circ } $ 与 $ 20^{ \circ } $ 时。而当θ = 40°与60°时，上下翼面正负压区面积增大不明显，而 $ \theta $ 增大使得z方向（图1（c））受力面减小，反而使 $ {C}_{L} $ 变小。

图15为不同迎角下， $t=T/2$ 时上下翼面的压力云图。从各个压力场情况可以看出，下翼面压力小于上翼面，此时 $ {C}_{L} $ 达到最小值。在 $t=T/2$ 时，随着迎角的增大，上下翼面的正压、负压区域逐渐变小，上下翼面的压力差逐渐趋于0，导致了 $ \theta =40^{ \circ } $ 与 $ 60^{ \circ } $ 的情况下， $ {C}_{L} $ 接近于0的现象。

图16为不同迎角下、 $t=0$ 、 $t=T/2$ 时上、下翼面沿展长压力系数 $ {C}_{p} $ 分布和压力差 $ {\Delta }p $ 分布。由图16（a、c）可以看出，扑翼沿翼根到翼尖方向， $ {\Delta }p $ 绝对值逐渐增大，到翼尖处时 $ {\Delta }p $ 绝对值达到最大值，同时随着迎角的增加， $ {\Delta }p $ 绝对值逐渐增大，当 $ \theta =40^{ \circ } $ 与 $ 60^{ \circ } $ 时， $ {\Delta }p $ 绝对值增大值较小，而迎角增大使得z方向（图1（c））受力面减小，反而使 $ {C}_{L} $ 变小。

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

2.2.2 流场及涡结构分析

为了更直观地展现三维流场中由扑翼产生的涡系结构，图17给出了在一个周期内4个时刻的Q准则云图。Q准则定义为 $Q=\dfrac{1}{2}({\varOmega }_{ij}{\varOmega }_{ij}-{S}_{ij}{S}_{ij})$ ，其中 ${\varOmega }_{ij}=\dfrac{1}{2}\Bigg(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}-\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\Bigg)$ 为旋转张量， ${S}_{ij}=\dfrac{1}{2}\Bigg(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}+\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\Bigg)$ 为应变率张量。

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

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

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

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

3 结　论

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

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

2） $\theta =0^{ \circ }\sim 40^{ \circ }$ 范围内，随着迎角的增加，扑翼上下翼面正负压区面积增大，上下翼面压差增大，扑翼获得的升力增加；翼尖处压差最大，给扑翼提供了主要的升力。

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

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