扑翼飞行器同时具有隐蔽性与机动性两种特征，符合现代战争隐蔽侦察与精确打击的概念，具有良好的军事用途前景[1]，因此以美国为首的多个国家[2]，近二十几年来掀起了一股扑翼飞行器研究的浪潮[3]。随着计算机技术的发展，研究扑翼飞行器的方法也愈来愈多。本文在前人对鸟类扑翼结构研究的基础上[4]，使用计算流体力学方法，建立三维鸽状刚性翅与柔性翅翼展模型，着重通过对比不同翼展的气动力参数，给出柔性翅扑动过程中的气动特性。
1 扑翼模型与扑动函数
1.1 扑动简介
自然界鸟类与昆虫的扑动方式各异，产生的影响也各不相同，但经过分解，扑动过程的运动形式主要有以下4种：① 翅翼绕着翅根的上下扑动; ② 翅面的俯仰或翅翼展向的扭转; ③ 翅翼的前后向运动； ④ 翅膀沿展向弯曲折叠。其运动形式如图1所示。
图1 扑翼运动的4种运动形式示意图
而鸟类的扑翼运动主要由翅翼的上下扑动与翅面的俯仰运动组成，因此本文也以这两种运动的特征作为数值计算的输入条件。
1.2 刚性翅模型
刚性翅模型的建立主要参照鸽子的翅翼形状及其特点，典型鸽子的飞行参数与翅翼特征在表1中列出。
表1 翅翼的特征参数
参数符号典型数值典型鸽子翅翼特征参数弦长c0.10m展长b0.3m翼展面积S0.06m2飞行速度v11m/s肱骨骨架直径d0.01m仿生翅翼翼型参数弦长c0.1m展长b0.3m前缘半径r0.005m最大厚度w0.016m最大厚度位置p20%前缘长度q0.13m后缘长度h0.3m展弦比λ3
根据典型鸽子翅翼特征参数与具体鸽子的翅翼形状，进行翅翼建模，模型如图2所示。
图2 鸽翼与刚性翅模型
根据鸽翅的生物学构造，鸽翅的前缘部分主要由肱骨支撑，量得肱骨直径约为0.01 m。因此仿生翅的前缘半径选择为0.005 m。翅翼的后缘部分为羽毛尖端，厚度可以忽略不计，因此构建模型时后缘厚度为0。鸽翅的厚度最大处位于肱骨与副羽的相接部分，约在弦向20%的位置，厚度为0.016 m。翼展前缘长度约0.13 m，后缘长度为0.3 m。鸽翅的下表面展开时没有明显的弧度，为平面型下表面。具体参数见表1所示。
1.3 柔性翅模型
建立柔性翅模型时考虑到后期与刚性翅气动特性的对比，保持柔性翅的外形与刚性翅完全相同，改变材料属性即可。同时，为了探究柔性翅支架对其气动特性的影响，这里建立2种不同的柔性翅模型：① 橡胶翅翼，0.001 m厚钢架支撑： ② 橡胶翅翼，0.002 m厚钢架支撑：
其中模型(1)、模型(2)的外形与刚性翅并无差异，只是改变了其材料属性，其外形如图3所示，并且多出一个支架部分对翅翼起到支撑作用。模型(1)、模型(2)的差异则主要在于钢架的厚度不同。
图3 柔性翅翼及其支架示意图
对翅翼建模完成后，构建外流域，外流域大小与流域流体的流速有关，由于是低速飞行，雷诺数较小，流场变化范围不大，因此取外流域的长宽大约是翅翼的6倍即可。构建后的外流域如图4所示。
图4 计算域
1.4 扑动函数
确定扑动函数时主要考虑的因素为扑翼飞行器的扑动幅值以及扑动频率，由于翼型是根据鸽子的翅翼确定的，为了使结果具有实际参考价值，扑动范围以及扑动频率也依照鸽子的飞行特点进行确定，鸽子的扑动频率在8 Hz左右，鸽子的主翅翼的扑扇幅度在45°左右，具体设置函数时，采用 8 Hz 的扑动频率，45°作为扑动幅值。由于扑动函数的确定过程较为简单，方法成熟[5]，这里直接给出扑动函数为：
(1)
求导得到其角速度函数：
(2)
2 计算前处理
2.1 计算方法
单纯刚性翅的计算方法较为简单，计算过程主要采用动网格技术[6]，这里不再过多赘述，其具体过程可参看文献[6]。柔性翅的计算则涉及流固耦合技术，流固耦合的核心在于计算时既要解得流场相关的物理量又要解得固体力学场的相关物理量，且不能忽略二者的相互作用[7]。简单地说，就是在每一步求解中，都需要解两类不同的方程并完成数据传递。因此流固耦合的计算量远远大于一般单项计算的情况，但其结果也是可观的，其求解出的结果往往更接近真实的物理现象[8]。本文采用双向耦合技术，即在翅翼在流场作用下发生形变，形变后翅翼的相关参数传给流体域，流体域接收流体中翅翼形变数据，根据翅翼的新状态进行下一时间步的流场参数计算。
2.2 计算模型
在本文的计算过程中，刚性翅的计算模型较为简单，其牵涉的计算模型只是柔性翅的一部分，而柔性翅计算时，既包含流体域计算，也包含固体力学计算，另外还涉及两者之间的数据交换，下面对柔性翅的计算模型作以详细说明：
1) 流体控制方程：
流体控制方程包括连续性方程与运动方程。连续性方程[9]如下。
(3)
运动方程：
(4)
式中： ρ为密度；u代表流体微团的速度矢量；F表示作用在单位质量上的质量力；P为应力张量。
2) 固体控制方程：
固体控制方程一般由牛顿第二定律导出：
ρd=▽·σ+F
(5)
式中： ρ是固体密度；σ是柯西应力张量；F是体积力矢量；d是固体域当地加速度矢量。
3) 流固耦合方程：
流固耦合的过程即为数据传递与应用的过程，该过程中则必须满足流体域与固体域的一些重要物理量的相等与守恒，即满足如下方程：
(6)
式中： τf、τs分别代表流体域与固体域力大小； nf、ns则接触表面的法向单位向量； df、ds代表流体域与固体域微团位移矢量； qf、qs分别代表流体域与固体域的热流量； Tf、Ts分别代表流体域与固体域的温度。
第一项则表示交界面上的流体域与固体域力相等，第二项表示位移相等，第三项表示热流量相等，第四项表示温度相等。
2.3 计算平台的搭建
ANSYS进行流固耦合计算时，在workbench中搭建如图5所示的平台。
图5 计算平台示意图
其中geometry模块用于翅翼模型的建立并共享给后面的两个计算模块。计算过程为瞬态过程，结构体的力学计算使用Transient structure，流体部分CFX与Fluent都具有流固耦合的能力，但Fluent进行动网格计算时更有优势，因此这里流体部分使用Fluent模块，最后采用耦合器system coupling。将结构体与流体的setup部分与system coupling相联接，计算过程中两者在耦合器中进行数据传递。
2.4 计算域网格划分
平台搭建完成后，需分别对流体与固体计算域进行网格划分。网格划分时，需要根据不同网格的特点进行考虑：六面体网格可以根据计算模型的结构进行块划分后生成，其计算速度快，准确性高，但其不适用于形状复杂的几何模型且在动网格中适应性差，不宜收敛。四面体网格生成速度快，相对内存占用较大，但适用于边界条件有巨大变化的动网格情况。因此这里采用四面体网格的划分策略。
流体域的网格划分可用Workbench中自带的mesh模块进行过操作，也可联接ICEM模块进行网格划分，之后将网格直接导入mesh模块。结构力学部分的网格划分过程相对简单，选择对应的body，输入控制网格大小参数即可，参数设置时尽量保持与流体接触面网格大小相似，这样能够适当减小插值次数，加快运算速度。网格仍然选取四面体网格，因为翅翼会发生明显变形，六面体网格不适合复杂形状结构体的运算。具体划分结果如图6所示。
图6 网格划分情况
3 计算过程
3.1 翅翼扑动升、阻力变化特点
在FLUENT中设置监测时，一般得到的是升阻力系数cl、cd。无量纲量升阻力系数可以对比分析升阻力大小。
考虑到鸽子的飞行速度大致在11 m/s，且鸟类在实际飞行时很少具有绝对的平飞状态，因此在FLUENT中设置气体来流速度11 m/s，设置攻角为8°，可以监测到刚性翅以及两种柔性翅升阻力系数的变化情况。
从图7、图8可以看出，0.001 m支架柔性翅的升力、阻力系数变化失去了刚性翅时域上的规律性，但从周期性的变化规律来看，趋势却是与刚性翅类似的，无论是升力还是阻力系数，都符合刚性翅“高-低-高”的变化特征。另外在时域上，一个周期内正升力部分所占的时域宽度要明显大于负升力。阻力系数的分布于升力部分也是类似的。
从图9、图10可以看出，0.002 m支架柔性翅翅翼的升阻力系数的峰值得到了大幅度的提升，另一方面 2 mm柔性翅翅翼的升阻力系数与 1 mm柔性翅翅翼的升阻力系数有类似的特点，即一个周期内整个变化曲线在时域上分布不均匀，依然出现了正的升阻力系数的时域范围明显大于负升阻力系数的时域。
图7 刚性翅与0.001 m支架柔性翅升力系数曲线
图8 刚性翅与0.001 m支架柔性翅阻力系数曲线
图9 刚性翅与0.002 m支架柔性翅升力系数曲线
图10 刚性翅与0.002 m支架柔性翅阻力系数曲线
3.2 波动时域特征
刚性翅与柔性翅的升阻力系数都呈周期性变化，柔性翅的峰值到来的时间略晚于刚性翅，且峰值大于刚性翅。另外，周期内正升阻力部分所占的时域宽度要明显大于负升阻力。这种现象对于平均升力系数来说，明显是有益的，下面对产生该种现象的原因作以分析，首先以两种柔性翅的升力为例，对升力曲线的时域作以划分，如图11、图12所示。
图11 0.002 m支架柔性翅升力系数曲线
图12 0.001 m支架柔性翅升力系数曲线
对比刚性翅与0.002 m支架柔性翅的时域可以发现，两者的上升与下降趋。势基本是同步的。而出现时域不对称的原因，主要是由于柔性翅部分出现了波动时域，波动时域的长度一定程度上“占用了”柔性翅的负时域，从而出现了负时域时域变窄的情况。而在0.001 m柔性翅中也同样出现了波动时域，在图12中可以看出，波动时域的长度要大于0.002 m支架柔性翅的波动时域长度，波动时域不但减小了柔性翅负时域的长度，且波动时域的延长直接导致了负升力的出现时刻发生了延迟。下面具体分析出现波动时域的成因与原理。
首先，依据柔性翅的升力系数变化图线，可知0.001 m支架柔性翅升力图线波动时域长度在0.06～0.08 s，0.002 m支架升力柔性翅图线在0.05 s与0.10 s附近，下面给出两种柔性翅翼在x=0.09 m处ZOY截面上的速度矢量图。
结合1.3节中的扑翼角度函数，分析两种柔性翅t=0.05 s时的状态。t=0.05 s时翅翼刚刚经过下极限位置，整个翅翼正在向水平位置复位，此时展项中前端由于受到惯性力与气动力的影响，翼展的大部分区域还处在大变形的状态。首先分析0.001 m支架柔性翅的流场情况(见图13)，0.001 m支架柔性翅在0.05 s时，尽管此时翅翼的翼根部已经过了极限位置并正在向水平位置复位，但结合扑动函数易知，翼根部位在复位前期，角速度极低，是整个扑动过程中，角速度最接近0的位置。而此时翅翼的前端与中部受到惯性力与气动力的影响，继续向下运动，到0.06 s时能够明显观察到翅的变形。从流场速度矢量可以明显地看出，在该截面中大部分气流在ZOY平面上的投影都指向了水平方向，只有少部分的气流在翼根处指向y轴的正方向，而此时翼根处一方面角速度较低，另一方面距离旋转中心较近，因此实际翼根处的旋转速度非常低。因此此时翅翼产生的力，大部分是水平力，如果翅翼是对称分布，那么水平方向的力是完全能够抵消掉，不会对飞行过程造成影响。且0.06～0.08 s期间，翅翼持续变形，速度矢量在ZOY平面的投影大部分都持续指向水平方向。这种持续变形的现象反映到升力系数变化图线上，即为波动时域的在周期内的时间占比较长。
类似地，分析0.002 m支架柔性翅的波动时域对应的流场速度矢量图(见图14)。产生波动时域的原因与0.001 m支架柔性翅的原因相同，即柔性翅翅翼中前端的持续运动导致了一定时间内，翅翼主要产生了水平力，而升力则基本在0点附近徘徊。但与0.001 m支架柔性翅不同的是，0.002 m支架柔性翅出现变形的时间更早，在0.05 s附近，翅翼已经明显变形，且恢复变形的速度也更快，在0.06 s时，流场速度矢量的已经指向了Y轴，或者说流场速度矢量在Y轴的分量在0.06 s时已经不容忽视，翅翼产生的主要力已经由水平力转变成了升力。0.10 s时，翅翼的变形过程、波动时域产生的原因与0.05 s时均是同理，这里便不再赘述。
图13 0.001 m支架柔性翅波动时域流场图
图14 0.002 m支架柔性翅波动时域流场图
造成两者变形过程不同的因素是显然的，二者流场相同、扑翼函数相同、翼型相同，只有支架的厚度相差了0.001 m。而支架影响翅翼恢复形状的原因则主要是由于弹性变形的回弹力矩产生了巨大差异。下面根据弹塑性力学理论，作以简要分析。
由于该支架的厚度只有0.001～0.002 m，其厚度远远小于整体支架的长度，因此可以作为薄板处理，根据弹性力学的薄板理论[10]，薄板的内力矩计算公式如下：
(7)
式中：Mx、My、Mz是作用在中面上每单位长度的弯矩与扭矩；h为薄板厚度；ω为薄板挠度；D为板件的抗弯刚度；κx、κy为变形后中面沿x、y方向的曲率；κxy为扭曲率；v为泊松比且有
(8)
在式(7)中，D为板件的抗弯刚度。且有：
(9)
其中：h为板件的厚度；E为弹性模量；v为泊松比。根据以上公式可知，若0.001 m支架柔性翅与0.002 m支架柔性翅发生同样程度的变形，那么0.002 m支架柔性翅产生的内力矩将是0.001 m柔性翅内力矩的8倍。
再有根据公式：
(10)
(11)
式中： β为角加速度； I为转动惯量，显然0.002 m支架柔性翅的转动惯量在同样的变形情况下是0.001 m支架柔性翅转动惯量的2倍。综合考虑转动惯量与内力矩，可知0.002 m支架柔性翅在同样的变形情况下的其恢复变形的角加速度是0.001 m支架柔性翅的4倍。因此也就出现了如图12所呈现的变形与恢复变形的情况，0.001 m支架柔性翅的翅翼变形恢复时间远长于0.002 m柔性翅翅翼。
综上所述，波动时域产生的原因是：柔性翅在极限位置附近发生了变形，导致翅翼的中上段延迟到达原本的既定位置，在该过程中，翅翼产生的主要力为水平力，而在双翅翼的情况下，水平力相互抵消。升力部分则主要由翼型与翼根部分对气流的作用决定，且其作用力很小。而柔性翅支架的形状与材料，决定了柔性翅的变形恢复速度，也就决定了波动时域的长短。
4 试验
本文扑翼研究始终以鸽子翼为仿生对象，展开了一系列的数值计算，得出了相关刚性翅与柔性翅的相关结论。在翅翼形变的基础上提出了波动时域理论。为了验证波动时域理论在鸟类飞行过程中是否具有实际意义，设计了高速摄影实验以拍摄鸽子的扑动过程。
由于鸽子的扑动过程并非一个固定的函数，为了方便对比，将鸽子的扑动过程以周期T来表征，结果如图15所示。
图15 鸽子的实际扑动与柔性翅仿真过程
从图13可以看出，在下扑阶段的变形过程中，鸽子翅翼的变形过程与数值计算的结果十分相似，但在翅翼上扑阶段，即0.4t～0.8t之间，鸽子翅翼的变形过程与数值计算结果不同。这是由于鸽子在翅翼的上扑过程中，中前段翅翼一直尽量保持在竖直状态，也就是说，鸽子翅翼在上扑过程中，鸽子尽量放缓了翅翼恢复成一个平面的过程，使得翅翼在上扑过程中尽量不产生负升力。如果联系0.002 m支架柔性翅升力系数曲线图，可以认为，鸽子在上扑过程中利用身体肌肉与骨骼结构延长了第一段波动时域在整个周期的时间占比，缩短了负时域的长度。
因此，实际鸟类在飞行中不但运用了柔性翅的波动时域特征，而且其自身的骨骼结构提升了波动时域的利用率，更利于扑翼飞行。
5 结论
1) 柔性翅的升力特性优于刚性翅，柔性翅产生大升力的主要原因在于波动时域的存在，波动时域延迟了负升力的到来时间，缩短了负升力在周期内的时间占比。
2) 柔性翅的支架对柔性翅的波动时域与峰值特性影响较大，主要由其结构、材料特性决定。
3) 波动时域在自然界鸟类飞行起着重要作用，并且由于鸟类的自身骨骼结构影响，鸟类对波动时域的利用率更高。
[1] 李峰，叶正寅，白存儒.微型扑翼飞行器空气动力学研究进展[J].华东交通大学学报，2007(1):63-65.
[2] SHYY W，AONO R，CHIMAKURTHI R K，et al.Recent progress in flapping wing aerodynamics and aeroelasticity[J].Progress in Aerospace ences，2010，46(7):284-327.
[3] 杨文青，宋笔锋，宋文萍，等.仿生微型扑翼飞行器中的空气动力学问题研究进展与挑战[J].实验流体力学，2015，29(03):1-10.
[4] 朱宝.扑翼飞行机理和仿生扑翼机构的研究[D].南京：南京航空航天大学，2010.
[5] 王姝歆.微型仿生飞行机器人飞行机理及其仿生翅运动系统研究[D].南京：东南大学，2004.
[6] 刘旭博，于纪言.鸽状扑翼飞行器气动特性研究[J].兵器装备工程学报，2019，40(01):130-135.
[7] TAKIZAWA K，TEZDUYAR T E，KOSTOV N.Sequentially-coupled space-time FSI analysis of bio-inspired flapping-wing aerodynamics of an MAV[J].Computational Mechanics，2014，54(2):213-233.
[8] 杨林.非线性流固耦合问题的数值模拟方法研究[D].青岛：中国海洋大学，2011.
[9] 潘光，钟如意，宋保维，等.扑翼翼型流体动力特性数值计算与分析[J].水下无人系统学报，2012，20(1):9-13.
[10] 徐秉业.应用弹塑性力学[M].北京：清华大学出版社，1995.
[11] 刘旭博.仿生扑翼飞行器的气动特性与结构研究[D].南京：南京理工大学,2019.