引　　言

尽管微型飞行器技术有了很大的发展，但昆虫、鸟类等生物的飞行能力至今仍让人造飞行器望尘莫及，因此研究扑翼气动性能具有重要意义^[1-3]。在以往的研究中，通常把扑翼翼型作为刚性翼型处理，忽略了其弹性变形的影响。刚性翼型在超过其相应的临界攻角时，会发生严重的流动分离，从而影响飞行器的飞行状态。而昆虫、鸟类等生物的翅膀在飞行中往往可以旋转较大角度且仍能保持平稳飞行。这不仅与其特殊的扑翼运动飞行方式有关，还与其柔性翼型结构有关。例如，Tang、Shyy等^[4]采用梁模型研究了低Reynolds数下柔性平板翼型的气动性能，研究表明柔性结构在空气动力载荷作用下发生变形，表现出类似扑翼的俯仰运动，产生了明显的升力和推力。Visbal等^[5]采用高保真隐式大涡模拟方法，研究了表面微柔性薄膜的SD7003翼型气动性能，研究表明柔性对膜翼型上流动结构使得升力增加和失速延迟。张兴伟、周超英等^[6]通过弱流固耦合（仅考虑结构变形对流场的影响）方法分析了柔性扑翼变形的影响，发现适当的柔性变形可以提高升力、增大升阻比。王姝歆、周建华等^[7]以昆虫的翅膀为基础，对柔性翅模型进行了合理的假设与简化，建立了柔性翅模型并进行了实验研究和分析，结果表明，柔性翅的柔性变形有助于改善飞行时的空气动力学特性。Kang等^[8]研究了前缘处局部柔性翼的绕流问题，发现局部柔性翼型对提高升力有着积极的影响。进一步，陶真新等^[9]数值分析了吸力面为三段柔性结构的NACA0012翼型的气动性能，发现翼面上表面的柔性起到提高升力、抑制流动分离的作用。

尽管以上研究已涉及柔性翼型气动性能的研究，但他们往往仅关注展向柔性变形或局部柔性变形的情况，多简化为薄膜、梁等简单结构进行建模，而缺乏对全局柔性翼型气动性能的研究。本研究以柔性椭圆扑翼翼型为例，建立有限元模型，分析不同来流风速、迎角下柔性椭圆翼型的空气动力学特性。

1.   刚性扑翼翼型空气动力学仿真分析

1.1   计算模型与网格划分

为了验证建模、仿真分析方法的有效性以及参数设置的合理性，首先对刚性椭圆翼型模型周围的流场进行仿真分析。利用ANSYS Workbench开展双向流固耦合计算，建立刚性翼型周围的流场模型，采用Meshing模块划分网格，调用计算流体力学Fluent模块，基于有限差分格式进行流场计算。

为了便于对比验证，扑翼翼型参考Pesavento和Wang所采用的果蝇翅膀尺寸^[10]，典型果蝇质量为1 mg，平均翼展长约为0.2 cm，椭圆翼型长轴长为0.068 cm，椭圆率为0.25（图1）。流场尺寸为16 mm × 10 mm，椭圆位于流场中央，翼型迎角为27.5°。流体介质为空气，速度为2.94 m/s，方向垂直于入口边界水平吹向椭圆翼型，左侧为流场的速度入口，右侧为出口，上下边界为固定壁面。

流场区域分为两部分，椭圆壁后尾流区与外围流场区，流场整体采用三角形非结构网格进行划分，如图2所示。为提高求解精度，同时也为了更好地呈现椭圆壁后流场，在尾迹区进行网格加密处理；同时为了提高计算效率，外围流场区采用较为稀疏的网格。考虑到低Reynolds数下椭圆壁边界层的黏性作用，椭圆边界层使用结构化网格进行加密处理，边界层网格膨胀率为1.1，共设置20层，网格厚度从壁面向外逐层递增，网格总数为19808。采用正交质量检查项进行网格质量检查，网格质量达到0.97。

在Fluent分析中，设置瞬态计算。由于果蝇的尺寸小，速度较低，在实际飞行中，其Reynolds数范围一般为10²~10⁴，因此在流场分析中主要分析其椭圆翼型低Reynolds数下的空气动力学特性，计算模型采用$k{\text{-}}\omega$ SST湍流模型^[11-12]。k-ω SST湍流模型湍动能k和其比耗散率ω的关系为

变量k和ω的控制方程分别为

图  1  椭圆翼型几何模型：（a）流体域几何尺寸；（b）椭圆翼型几何尺寸

Figure  1.  The geometric model for the elliptical airfoil: (a) geometric sizes of the fluid domain; (b) geometric sizes of the elliptical airfoil

下载: 全尺寸图片 幻灯片

图  2  流体网格划分：（a）整体网格划分；（b）壁面附近网格加密处理

Figure  2.  Meshing of the fluid domain: (a) meshing of the whole model; (b) mesh refinement near the wall

下载: 全尺寸图片 幻灯片

方程（2a）和（2b）中右侧前三项分别为生成项、耗散项和扩散项，方程（2b）右侧最后一项为交叉扩散项，其中

式(1)~(4)中，$\varepsilon$为湍动能耗散率；$\;\mathrm{\rho }$为密度；${u}_{i}, {u}_{j}, {u}_{k}$为速度分量；$\mu$为层流黏性；${\nu }_{t}$为湍流运动黏性；$\;{\beta }^{*}$，$\beta$，γ，${\sigma }_{k}$，${\sigma }_{\omega }$和${\sigma }_{\omega 2}$为封闭常数。

采用SIMPLE的二阶格式求解方法对N-S离散方程进行求解，提高求解精度。时间步长为0.000 1 s，共计算100个时间步，计算得到0.01 s内的结果。

1.2   模型验证

分别计算得到该模型流场的速度云图和压力云图（图3），椭圆翼型后产生了明显的周期性尾涡脱落，这一特征与Pesavento和Wang^[10]数值计算的流场特征一致。

图  3  刚性翼型流场流速分布：（a）速度云图；（b）压力云图

注　为了解释图中的颜色，读者可以参考本文的电子网页版本，后同。

Figure  3.  Characteristics of the flow field around the rigid airfoil: (a) the velocity contour; (b) the pressure contour

下载: 全尺寸图片 幻灯片

为了便于与文献结果定量对比验证，将时间无量纲化为

其中，$v$为来流速度，$v$=2.94 m/s；$t$为时间；$c$为椭圆翼型的弦长，$c$=0.68 mm。

数值仿真计算所得椭圆翼型升力随时间的变化，开始阶段由于求解过程未收敛出现较大幅度扰动，随后升力值在0.4~0.6 mg之间呈现出周期性变化，周期约为4（图4）。这是因为翼型后缘出现周期性尾涡脱落。升力平均值约为0.5 mg，这是因为单只翅膀承担了果蝇重量。计算所得升力曲线的这些特征（扰动周期、平均升力），与Pesavento和Wang^[10]的数值计算结果保持一致。表明所建立的模型、计算方法及参数设置等是有效的，可供后续研究。采用三种网格数量（1.23×10⁴，1.98×10⁴，5.23×10⁴）进行了网格无关性验证，稀疏网格对升力峰值的监测不够，中等网格和密网格结果相近。因此，后续研究采用中等网格（即网格数1.98×10⁴）即可满足计算精度。

图  4  刚性翼型升力与Pesavento和Wang^[10]的数值计算结果对比

Figure  4.  Comparison of lift forces of the rigid airfoil to the numerical results of Pesavento and Wang^[10]

下载: 全尺寸图片 幻灯片

1.3   不同风速、迎角下的流场及气动性能

选取翼型迎角为0°，6°，12°，16°，20°，24°，28°，32°，风速为2 m/s，5 m/s，10 m/s，分别计算得到椭圆翼型周围的流场（图5），压力场（图6），升、阻力及升阻比（图7）。

图  5  不同迎角下的刚性翼型流场速度云图

Figure  5.  Velocity contours of the flow field around the rigid airfoil at different attack angles

下载: 全尺寸图片 幻灯片

从速度云图（图5）可以看出，随着攻角增大，椭圆翼型后流动分离范围逐渐增大。当迎角较小时，流动为稳定的层流；当迎角达到24°，翼型后的流动开始出现扰动；当迎角达到28°以上时，翼型后出现周期性的涡脱落。从不同迎角下的压力云图（图6）可以看出，当迎角增大到一定程度，由于脱落涡的出现，出现了脱离翼型的负压中心。

从升力曲线图(图7(a))可以看出，升力随来流风速的增大而增大；当迎角较小时，升力随迎角线性增大；随着迎角进一步增大，由于椭圆翼型后流动分离加剧，升力随迎角的增大呈非线性变化，迎角越大升力增大的速率越慢。从阻力曲线图(图7(b))可以看出，当迎角较小时，阻力随着迎角缓慢增大；当迎角增大达到16°左右时，阻力随迎角迅速增大。而升阻比随迎角变化成非单调趋势(图7(c))，随迎角的增大先增大再减小，当迎角为10°~15°时达到最大值。

图  6  不同迎角下的刚性翼型流场压力云图

Figure  6.  Pressure contours of the flow field around the rigid airfoil at different attack angles

下载: 全尺寸图片 幻灯片

图  7  不同风速下刚性翼型的升力、阻力及升阻比：（a）升力；（b）阻力；（c）升阻比

Figure  7.  Lift forces, drag forces and lift-drag ratios of the rigid airfoil at different wind speeds: (a) lift forces; (b) drag forces; (c) lift-drag ratios

下载: 全尺寸图片 幻灯片

2.   柔性翼型流固耦合仿真分析

2.1   流固耦合模型

在以上刚性翼计算模型的基础上，将椭圆翼型考虑为柔性翼。果蝇翅膀上分布着沿展向翅脉，且翅脉的刚度远大于翼面。为了简化模型，在柔性扑翼椭圆翼型模型中心的空心圆柱（直径为0.02 mm）面，采用固定边界，以模拟刚性翅脉（图8）。流体作用力与柔性翼型固体结构相互耦合，一方面流体在固体结构表面产生作用力，使固体结构产生变形；结构形状的改变，引起流场结构的改变。流固耦合仿真分析在ANSYS Workbench平台中用Fluent计算流体模块进行流场分析，用Transient Structural瞬态结构模块进行结构分析，使用System Coupling耦合求解器将两个模块关联起来，双向传递流体的计算数据和固体结构的计算数据，将流场与结构变形进行耦合。在流固耦合分析时，需要建立流固界面进行数据传递，因此柔性翼型模型应为三维模型。

柔性翼型固体结构分析的动力学控制方程^[9]可写为

式中${\boldsymbol{M}}$, ${\boldsymbol{C}}$, ${\boldsymbol{K}}$分别为质量矩阵、阻尼矩阵、刚度矩阵；${\boldsymbol{F}}\left(t\right)$为流体作用在固体结构表面上的瞬态耦合力，其对结构的综合作用可以是力，也可以是力矩。柔性材料密度为1 000 kg/${{\rm{m}}}^{3}$，材料本构模型采用线弹性模型，Poisson比为0.3，弹性模量可调。

图  8  柔性椭圆翼型有限元模型：（a）网格划分；（b）流固耦合面与固定面

Figure  8.  The finite element model for the flexible elliptical airfoil: (a) the meshing; (b) the fluid-solid coupling and the fixed surface

下载: 全尺寸图片 幻灯片

2.2   网格划分

流固耦合计算在Fluent中采用动网格技术，采用非结构网格进行划分。在流固耦合计算中，当固体结构发生变形后，流场需要接收固体变形的位移数据，流场边界也相应的改变，因此需要使用动网格进行网格重构。具体为在Dynamic Mesh中设置Smoothing与Remeshing，将椭圆壁面设为System Coupling，端部壁面设为Deforming，其余设置与刚性翼型一致。固体结构也采用非结构网格进行划分。

2.3   模型验证

在验证模型中，椭圆翼型迎角设置为27.5°，来流速度为2.94 m/s。柔性翼型弹性模量为0.01 MPa，0.1 MPa，1000 MPa，分析流固耦合模型的流场和固体结构变形。

从翼型周围流场速度云图（图9）可以看出，随着椭圆翼型弹性模量的减小，翼型周围的最大流速减小，但是翼型对周围流场的扰动范围扩大。类似地，从翼型周围流场压力云图（图10）可以看出，随着椭圆翼型弹性模量的减小，翼型周围产生的最大负压值减小，但是翼型周围的负压区扩大。从柔性翼型结构位移云图（图11）可以看出，随着椭圆翼型弹性模量的减小，柔性翼型结构位移增大。

图  9  不同弹性模量的柔性翼型流场速度云图

Figure  9.  Velocity contours of the flow field around the flexile airfoil with different Young’s moduli

下载: 全尺寸图片 幻灯片

图  10  不同弹性模量的柔性翼型流场压力云图

Figure  10.  Pressure contours of the flow field around the flexile airfoil with different Young’s moduli

下载: 全尺寸图片 幻灯片

图  11  不同弹性模量的柔性翼型位移云图

Figure  11.  Deformation contours of the flexile airfoil with different Young’s moduli

下载: 全尺寸图片 幻灯片

图  12  柔性翼型的升力和阻力：（a）升力；（b）阻力

Figure  12.  Lift forces and drag forces of the flexible airfoil: (a) lift forces; (b) drag forces

下载: 全尺寸图片 幻灯片

图  13  不同迎角下的刚性翼型（左）与柔性翼型（右）的周围流场

Figure  13.  Velocity contours of the fluid field around the rigid (left) and the flexible (right) airfoil with different attack angles

下载: 全尺寸图片 幻灯片

从升力曲线图（图12（a））可以看出，柔性椭圆翼型升力开始阶段由于求解过程未收敛出现较大幅度波动，随后升力值稳定在0.4~0.6 mg之间呈现出周期性变化，总体变化趋势与刚性翼型一致。但随着弹性模量的降低，升力值稳定扰动振荡的周期增大。而随着弹性模量的降低，扰动振荡的振幅呈先减小再增大的趋势。当弹性模量降低到0.1 MPa时，升力值扰动振荡幅度几乎消失，升力值达到一稳定值（约0.5 mg）；而当弹性模量进一步降低时，升力值又呈现出周期性的震荡变化。由此可见，较刚性翼型，柔性翼型对于稳定升力扰动有着积极的作用。类似地，从阻力曲线图（图12（b））中可以看出，较刚性翼型，柔性翼型的阻力略大于刚性翼。且随着弹性模量的降低，柔性翼型的阻力逐渐出现了周期性振荡。

2.4   不同迎角下流场及气动性能

柔性翼型的迎角分别选取0°，6°，12°，16°，20°，24°，28°，32°进行仿真分析。风速为5 m/s，柔性翼型弹性模量为0.1 MPa。

由翼型周围流场（图13）可以看出，随着迎角逐渐增大，刚性翼型后的流场逐渐出现扰动，翼型后伴随有尾涡脱落。而柔性翼型后的流场与刚性翼型有很大不同，流体流过柔性翼型后尽管出现了流动分离，但流场并未出现明显的扰动。由此可见，较刚性翼型，柔性翼型延缓了尾涡脱落时间，显著抑制了尾流流场的扰动范围及幅度。

由柔性翼型的升力曲线（图14（a））可以看出，升力随着迎角增大而增大，且柔性翼型的升力略大于刚性翼型。由柔性翼型的阻力曲线（图14（b））可以看出，阻力随迎角的增大而增大，且柔性翼型的阻力略大于刚性翼型。但从升阻比（图14（c））看，柔性翼型与刚性翼型并无相显著差异。也就是说，柔性翼型并未损失飞行效率。

图  14  柔性翼型的升力、阻力及升阻比随迎角变化：（a）升力；（b）阻力；（c）升阻比

Figure  14.  Lift forces, drag forces and lift-drag ratios of the flexible airfoil vs. the attack angle: (a) lift forces; (b) drag forces; (c) lift-drag ratios

下载: 全尺寸图片 幻灯片

3.   结　　论

本文从仿生学的角度出发，通过流固耦合数值仿真分析了柔性翼型在不同风速、迎角下的空气动力学特性。研究发现，较刚性翼型，柔性翼型延缓了尾涡脱落时间，有效降低升力扰动振荡频率；柔性翼型显著抑制了尾流流场的扰动，降低升力扰动振荡幅值，合适的弹性模量翼型使得扰动振荡完全消除，这将有利于改善飞行器运动平顺性、降低气动噪声。

致谢　本文作者衷心感谢江苏省仿生功能材料重点实验基金对本文的资助。