0 引　言

近年来，微型飞行器（Micro Aerial Vehicle, MAV）因其潜在的军用和民用前景，成为航空领域研究的热门^[1-2]。MAV的雷诺数很小，通常在100～10 000之间。低雷诺数下流体的黏性效应十分显著，导致MAV的阻力大幅增加。然而，在相近的雷诺数下，自然界中的鸟类、昆虫和鱼类依赖其特有的扑翼飞行或游动方式实现自主推进^[3]。为优化MAV设计，人们围绕扑翼的推力特性及其背后的流动机理开展了大量研究。

目前已开展的大部分工作基于大展弦比假设对绕扑翼流动进行了二维近似。Koochesfahani^[4]利用激光多普勒测速（Laser Doppler Velocimetry, LDV）技术对NACA0012翼型作小幅度俯仰运动的尾迹流场进行了测量，发现缩减频率较小时，尾迹形态为卡门涡街，其中顺时针旋转的涡处于逆时针旋转的涡上方；随着缩减频率的增大，顺时针涡与逆时针涡的纵向位置交换，前者位于后者下方，该尾迹形态称为反卡门涡街。Jones等^[5]及Von Ellenrieder和Pothos^[6]对NACA0012翼型作浮沉运动的实验研究表明，当斯特劳哈尔数大于某一临界值时，出现尾迹不对称现象，即反卡门涡偏离水平线向上方或向下方倾斜。Godoy-Diana等^[7]实验研究了俯仰翼型尾迹不对称性机制，建立了反卡门涡对称性破坏准则。Andersen团队^[8-9]对俯仰翼型和浮沉翼型的实验及数值研究发现，在改变振荡频率和幅度时，除观察到卡门涡街（2S BvK）和反卡门涡街（2S RBvK）外，还出现2P型、2P+2S型、4P型、4P+2S型等更复杂的尾迹形态，其中“S”表示孤立涡，“P”表示由两个旋转方向相反的涡组成的涡对，“2”和“4”为每个振荡周期内从翼型后缘脱落到尾迹中的孤立涡或涡对的数目。

尾迹形态随扑翼运动学、动力学的演变是扑翼流场的主要特征，其与推力产生的关系是研究的重点。Von Karman和Burgers^[10]在对平板横向振荡的研究中发现，反卡门涡街的“喷流”现象能够产生推力。Koochesfahani^[4]的研究表明，当俯仰翼型尾迹形态由卡门涡街转变为反卡门涡街时，对应的尾迹流场由速度损失型变为速度增加型或“喷流”。他进一步基于积分动量定理证实，反卡门涡街是一种推力产生型尾迹。Jones等^[5]在对浮沉翼型的研究中也发现，随着无量纲浮沉速度的增大，发生阻力型尾迹(卡门涡街)到推力型尾迹(反卡门涡街)的转变。Lai和Platzer^[11]在此基础上研究了浮沉翼型尾迹的“喷流”特性，指出最大喷流速度是无量纲浮沉速度的线性函数。

尽管反卡门涡街的“喷流”理论揭示了扑翼产生推力的流动机制，但随着研究的深入，人们发现反卡门涡街的出现与推力产生并不是完全对应的。Godoy-Diana等^[12]采用粒子图像测速（Particle Image Velocimetry, PIV）技术观察了俯仰翼型的尾迹形态（ $Re = 1\;173$ ），首次发现随着斯特劳哈尔数增大，卡门涡街-反卡门涡街的转变先于阻力-推力转变，原因是翼型扑动产生的推力不足以克服其受到的阻力，Deng等^[13]基于计算流体力学（Computational Fluid Dynamics, CFD）方法对此作了类似的解释（ $Re = 1\;700$ ）；他们均未对阻力的来源进行深入分析。Das等^[14]数值研究了雷诺数对俯仰翼型推力特性的影响，在10～2000范围内均观察到了这一滞后现象，但未就其原因进行讨论。Bohl和Koochesfahani^[15]对原有积分动量定理作了改进，并用于估计NACA0012俯仰翼型的推力大小（ $Re = 12\;600$ ），发现在推力为零时尾迹速度型分布已经呈“喷流”状，认为该额外的动量通量用以克服俯仰翼型下游的压力损失。Ashraf等^[16]进一步对NACA0015俯仰翼型的流场进行了PIV测量（ $Re = $ $ 2\;900$ ），应用改进的积分动量定理估计了推力的大小，结果显示反卡门涡街的形成不能保证推力的产生，原因是其导致的尾迹流场动量增加被翼型运动引起的压力损失抵消了。积分动量定理仅涉及位于翼型后缘下游某一特定位置的控制面上的速度分布，因此文献[15-16]并未考虑控制体内部的非定常尾迹流场对推力产生的影响。

目前对扑翼推力特性的研究取得了一些重要结论。Sarkar和Venkatraman^[17]研究发现，俯仰翼型的时均推力随缩减频率的增大而增大，随平均迎角的增大而减小。Mackowski和Williamson^[18]的研究表明俯仰翼型阻力-推力转变对应的临界缩减频率随雷诺数的增大而减小。Das等^[14]进一步指出，俯仰翼型阻力-推力转变对应的临界斯特劳哈尔数 $S\!{t_{\rm{tr}}}$ 与雷诺数之间存在简单的比例关系St_tr～ $R{e^{ - 0.37}}$ 。Deng等^{[13, 19-20]}通过对二维基准流动施加展向周期性小扰动，结合Floquet理论分析了俯仰翼型尾迹流场的稳定性，发现随着振荡频率或幅度增大，主模态依次经历了由稳定模态、临界稳定模态到不稳定模态的演变过程，其中不稳定模态反映了流场的固有三维线性不稳定性，相应的尾迹称为三维尾迹，主模态为临界稳定模态的尾迹为二维尾迹-三维尾迹的转变边界；文献[13]进一步开展了三维直接数值模拟，验证了其与二维Floquet分析的一致性；文献[20]讨论了尾迹形态与推进效率的关系，指出俯仰翼型的推进效率峰值与二维尾迹-三维尾迹转变边界吻合较好。田伟等^[21]的实验研究发现，俯仰轴位置对尾迹涡演化及推力特性有较大的影响。戴龙珍和张星^[22]对比研究了扑翼在连续式和间歇式俯仰运动下的能耗问题，指出间歇式俯仰运动在推进速度较高时能耗更低，连续式俯仰运动则在推进速度较低时更节省能耗。

尽管上述工作对扑翼的流场特征和推力特性的认识已相对深入，但目前还没有相关文献针对在低雷诺数下，随着斯特劳哈尔数增大，俯仰振荡翼型的推力产生滞后于反卡门涡街形成的现象进行系统性研究，尤其对其背后的流动机制的认识存在不足。为此，本文采用CFD方法对NACA0012翼型在1 000雷诺数下作简谐俯仰运动的非定常流场进行了数值模拟，考察了阻力-推力转变临界点附近的尾迹形态特征和推力特性。采用翼型表面积分方法和基于有限控制体的气动力估计方法分别研究翼面分布力特性和尾迹流场特性变化对推力产生的影响，以探究阻力-推力转变滞后于卡门涡街-反卡门涡街转变的物理机制。

1 数值计算和验证 1.1 控制参数

考虑NACA0012翼型以正弦规律绕1/4弦长位置作纯俯仰运动：

${\alpha} ={\alpha _0}\sin{ (2 {\text{π}} ft)}$

(1)

其中， ${\alpha _0}$ 为俯仰幅度， $f$ 为俯仰频率。利用雷诺数，无量纲幅度和斯特劳哈尔数表征绕扑翼流动，其中无量纲幅度和斯特劳哈尔数分别定义为：

${A_D}{\rm{ = }}A{\rm{/}}D{\rm{ = 1}}{\rm{.5}}c\sin {\alpha _0}{\rm{/}}D$)

(2)

$S\!t{\rm{ = }}fD{\rm{/}}{U_\infty }\;\;{\rm{ or }}\;\;S\!{t_A}{\rm{ = }}fA{\rm{/}}{U_\infty }{\rm{ }}$

(3)

其中，D为翼型厚度，A为在一个俯仰周期内翼型后缘通过的垂直距离。

本文选取俯仰幅度（ ${\alpha _0}$ ）为2°～10°，对应的无量纲幅度（ ${A_D}$ ）为0.44～2.17；对每个 ${\alpha _0}$ ，通过改变俯仰频率使斯特劳哈尔数（ $S\!t$ ）在0.08～0.49范围内（对α₀ = 2°，St = 0.08～0.9）。由此形成包含81个算例的数据集，内容涵盖了无量纲幅度和斯特劳哈尔数在不同取值下的推力系数及其分量，以及相应的流场信息。

1.2 求解设置

计算域为足够大的圆形区域，直径为100倍翼型弦长，圆心位于1/4弦长位置，如图1所示。采用481（周向）×241（径向）的O型结构网格对计算域进行离散，第一层网格高度为0.001c，满足y⁺ < 1。对直径为10倍弦长的内部圆形区域加密（径向），使其具有较高的流场分辨率；该圆形区域之外为较粗化的网格。两块网格的周向节点数一致，且在交界面上的节点一一对应。

采用Fluent 17.2压力基求解器数值求解二维不可压N-S方程。远场设置为速度入口边界条件，来流速度为0.146 m/s，相应的雷诺数为1 000。选择该雷诺数的原因是，一方面 $Re = 1\;000$ 处于MAV的雷诺数范围内^[3]，另一方面已开展的大部分研究中雷诺数的选择都集中在 ${10^3}$ 量级^{[8-9, 12-14, 16, 20]}。翼型为无滑移壁面边界，两块网格交界面为内部网格面（Interior）。该雷诺数下流体黏性较强，选用层流模型进行模拟。翼型运动由用户自定义功能（User Defined Function, UDF）文件指定，整个流体域固连在翼型上，随翼型作相同的刚性运动。

压力-速度耦合求解采用SIMPLEC算法，空间离散中对流项和黏性项分别采用二阶迎风格式和二阶中心格式，时间项采用二阶隐式双步法，对每个俯仰周期分N个时间步进行时间推进，每个时间步内进行25次内迭代。为获得周期性稳定解，每个算例均至少计算60个俯仰周期。

1.3 相关验证

采用四套疏密程度不同的网格验证网格无关性，网格信息见表1。验证算例中选取α₀ = 10°， $S\!{t_A}{\rm{ = 0}}{\rm{.5}}$ 。如图2（a）所示，除粗化网格Grid 1外，两套中等网格Grid 2、Grid 3与细化网格Grid 4的阻力系数和力矩系数的时间历程曲线十分接近，Grid 3和Grid 4的周期平均阻力系数和力矩系数的相对误差低于1%。为验证时间步长无关性，选取4个依次加倍的N值，即时间步长依次减半。图2（b）显示，N为720和1 440时，计算结果差异较小。为保证较高的计算精度同时适当减小计算量，后续所有计算中都采用细化网格，时间步长取为T/ 720。

为进一步验证求解器的可靠性，将计算结果与文献中的结果进行对比。Deng等^[20]的CFD计算中 $Re = $ $ 1\;700$ ，选取其中 ${A_D}=2、St={\rm{0.5}}$ 的算例进行对比验证，为保持一致性，将模型替换为NACA0015翼型，俯仰轴选在前缘位置。图3显示，阻力系数和力矩系数的时间历程曲线与Deng等的结果十分吻合。

2 基于有限控制体的气动力分解方法

为分析俯仰翼型诱导形成的非定常尾迹对推力产生的影响，采用基于有限控制体的气动力估计方法建立局部尾迹与推力产生的联系。对不可压黏性流动，由动量定理，作用在翼型上的气动力为：

${{F}} = - \mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{\partial B} {( - p{{n}} + {{\tau }}){\rm{d}}S}$

(4)

${{F}} = - \rho \iiint_{{V_{\rm{f}}}} {{{a}}{\rm{d}}V}{\rm{ + }}\mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{\varSigma} {( - p{{n}} + {{\tau}} ){\rm{d}}S} $

(5)

其中， $ p{\text{、}}\rho {\text{、}}{{\tau }}$ 分别为压力、流体密度和黏性应力矢量， ${{a}} = {\rm{d}}{{u}}{\rm{/}}{\rm{d}}t$ 为流体加速度， ${{u}}$ 为流体速度。 ${V_{\rm{f}}}$ 为包围翼型的有限控制体，其内边界为翼面 $\partial B$ ，外边界为 $\varSigma $ ， ${{n}}$ 为指向控制体外部的单位法向量。方程（4）是在翼面直接积分求解气动力的常用方法，被大部分CFD求解器所采用，根据该式可以分析压力和黏性应力分布对气动力产生的贡献。方程（5）通过包围翼面的控制体间接求解翼型所受的气动力。Wang等^[23]在方程（5）基础上利用积分变换推导得到：

$ \begin{split} { F}\!=\! &\underbrace{\rho \iiint_{{{V}_{\rm{f}}}}{{ u}\times { \omega }{\rm{d}}V}}_{(1){\text{涡力}}({\rm{V}})}\underbrace{\!-\! \rho \iiint_{{{V}_{\rm{f}}}}{\frac{\partial { u}}{\partial t}{\rm{d}}V}}_{(2){\text{局部流体加速}}({\rm{A}})}\underbrace{-\mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{\partial B } {\frac{1}{2}\rho |{{u}}|^2}{ n}{\rm{d}}S }_{(3){\text{固体占据的虚拟流体}}({\rm{B}})}\!\!\!+\! \\ & \underbrace{-\mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{{\varSigma}} {\frac{1}{2}\rho |{{u}}|^2}{ n}{\rm{d}}S }_{(4){\text{动量诱导}}({\rm{M}})}\underbrace{- \mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{{\varSigma}} p{ n}{\rm{d}}S }_{(5){\text{压力诱导}}({\rm{P}})}\underbrace{+ \mathop{{\int\!\!\!\!\!\int}\mkern-18mu \bigcirc}\nolimits_{{\varSigma}} { \tau }{\rm{d}}S }_{(6){\text{粘性剪切}}({\rm{S}})} \\[-18pt] \end{split} $

(6)

其中， ${{\omega}} $ 为涡量。上式中每一项都有明确的物理意义^[23-24]：第（1）项为lamb矢量或涡力项；第（2）项和第（3）项分别为由翼型运动的非定常惯性效应引起的局部流体加速和翼型占据的虚拟流体对气动力的贡献，在无黏无旋流动中二者之和称为附加质量力；第（4）至第（6）项依次为外控制面上的动量诱导力、压力诱导力和黏性力。

通过实验方法得到压力场的空间分布并不方便，而对速度场的测量相对容易且测量技术较为成熟^[25]。Wang等结合不可压N-S方程进一步将压力场转化为速度场，同时选取特殊的矩形控制体以消除方程（6）中的压力项。在CFD计算中容易得到压力场，因此本文直接应用方程（6）进行气动力分解且不对控制体的选取作任何限制。为避免后处理中的数值插值和积分误差，将方程（4）和方程（6）写入UDF文件，基于原始网格在CFD求解计算的过程中同时完成对各气动力分量的计算。为减小由涡量耗散引入的积分误差和不确定性，控制体与图1中计算网格的加密区一致，该范围内涡量耗散程度较小。

3 结果和讨论 3.1 推力产生滞后于反卡门涡街出现

图4给出了零推力点附近不同俯仰幅度下周期平均推力系数随斯特劳哈尔数的变化，图中曲线为对各离散点的二次函数拟合，零推力点由拟合函数近似给出，并经计算验证满足推力系数在 ${10^{{\rm{ - 4}}}}$ 量级。图4（a）显示，在各俯仰幅度下，周期平均推力系数随 $S\!t$ 的增大而增大，在 $S\!t$ 减小至零时趋于某一定值，对应静止翼型所受到的阻力（图中未标出）；发生阻力-推力转变的临界斯特劳哈尔数 $S\!{t_{T{\rm{ = 0}}}}$ 随俯仰幅度的增大而减小。图4（b）显示不同俯仰幅度下推力系数相对于 $S\!{t_A}$ 的变化曲线较为接近，均在 $S\!t_A$ = 0.3附近出现零推力点。考虑到 $S\!t_A={A_D}S\!t$ ，因此阻力-推力转变临界斯特劳哈尔数 $S\!{t_{T{\rm{ = 0}}}}$ 与无量纲幅度 ${A_D}$ 近似为反比例关系。

图5对比了俯仰幅度为10°时不同斯特劳哈尔数下的尾迹形态。 $S\!t$ = 0.066时为典型的卡门涡街，其中顺时针旋转的涡位于逆时针旋转的涡上方。随 $S\!t$ 增大，顺时针涡和逆时针涡沿纵向靠拢，在 $S\!t$ = 0.082时基本处于一条直线上。此后随 $S\!t$ 增大顺时针涡与逆时针涡的纵向位置交换，前者位于后者下方，呈反卡门涡街排列。随 $S\!t$ 进一步增大，逐渐出现不对称尾迹和正反涡配对现象（ $S\!t$ = 0.185～0.247）。注意到在 $S\!t$ = 0.103时反卡门涡街就已经出现，但并未产生推力，直到 $S\!t$ = 0.143 时才达到零推力点。因此推力产生明显滞后于反卡门涡街的出现。

为较准确判断卡门涡街-反卡门涡街转变临界点，根据涡量极值点确定了后缘下游0.25～3.25倍弦长范围内所有尾迹涡的涡核位置 $({x_{{\rm{c,}}i}},{y_{{\rm{c,}}i}})$ ，定义涡核纵向平均偏离量为该范围内所有涡核纵向位置绝对值的平均值：

$ {\left|\: {{\overline y_{\rm{c}}}}\: \right|} {\rm{ = }}\frac{1}{M}\sum\limits_{i = 1}^{i = M} {\left| {{y_{{\rm{c,}}i}}} \right|} $

(7)

${\left|\: {{\overline y_{\rm{c}}}}\: \right|}$ 最小时的尾迹定义为临界尾迹。图6显示了不同俯仰幅度下 ${\left|\: {{\overline y_{\rm{c}}}}\: \right|} /D$ 随 ${S\!t_A}$ 的变化，除 ${\alpha _0}$ = 2°外，其余俯仰幅度下 ${\left|\: {{\overline y_{\rm{c}}}}\: \right|}$ 在 ${S\!t_A} = 0.18$ 附近达到最小值，由于正反旋转的涡在向下游运动时存在相互扰动，导致所有涡核的纵向位置并不全部为零，从而该最小值并不为零。 ${\alpha _0}$ = 2°时临界尾迹对应的 ${S\!t_A}$ 值为0.25左右，明显与其它俯仰幅度不一致，原因是更小幅度俯仰运动产生的尾迹涡强度较弱，在向下游运动时很快就耗散了，因此对涡核的提取存在较大误差，如图7所示。图7给出了不同俯仰幅度下由该方法确定的临界尾迹， ${\alpha _0}$ = 10°的临界尾迹已在图5中显示。

基于以上分析，图8绘制了尾迹形态随无量纲幅度 ${A_D}$ 和斯特劳哈尔数 $S\!t$ 的演变过程。由于 $S\!{t_A} = $ $ {A_D}S\!t$ ，因此当 $S\!{t_A}$ 取值一定时， ${A_D}$ 与 $S\!t$ 的变化关系对应 ${A_D}$ - St参数空间中的一条反比例函数曲线。可以发现，除 ${A_D}=\:0.44$ ( ${\alpha _0}$ = 2°）外，其他 ${A_D}$ 取值下卡门涡街（BvK）-反卡门涡街(RBvK)转变临界点即临界尾迹点（Neutral Wake）大致分布在曲线 $S\!{t_A}=\:0.18$ 上；类似地，曲线 $S\!{t_A}=\:0.37$ 基本对应反卡门涡街-不对称尾迹（Asymmetic Wake）转变临界。曲线 $S\!{t_A}=\:0.18$ 和 $S\!{t_A}=\; $ $ 0.37$ 可以将1 000雷诺数下NACA0012翼型俯仰运动诱导产生的尾迹形态划分为三个区域： $S\!{t_A}=\:0.18$ 左下方为卡门涡街区域， $S\!{t_A}=\:0.37$ 右上方为不对称尾迹区域，介于二者之间的带状区为反卡门涡街区域。图中同时绘制了阻力-推力转变临界点（ $\overline {{C_T}} = 0$ ），大致分布在曲线 $S\!{t_A}=\:0.3$ 上。曲线 $S\!{t_A}=\:0.3$ 处于反卡门涡街区域中间靠后位置，其左下方为阻力区，右上方为推力区。以上临界线位置与Godoy-Diana等 $({Re =} $ $ {1\;173})$ ^[12]及Andersen团队 $\left( {Re = 1\;320,{\rm{ 2\;640}}} \right)$ ^[8-9]的实验及数值结果接近，尽管选择的翼型形状不同。

图8显示，在每个俯仰幅度取值下（α₀ = 2°～10°或 $ {A_D}\:{\rm{ = }}\:0.44 $ ～2.17），卡门涡街-反卡门涡街转变临界斯特劳哈尔数 ${S\!t}^{\rm{*}}$ 和阻力-推力转变临界斯特劳哈尔数 $S\!{t_{T{\rm{ = 0}}}}$ 均不相等，且 $S\!{t_{T{\rm{ = 0}}}}$ 大于 ${S\!t}^*$ ，从而推力产生滞后于反卡门涡街出现。定义无量纲数 $\Delta {S\!t}^{\rm{*}}{\rm{ = }}S\!{t_{T{\rm{ = 0}}}} - {S\!t}^{\rm{*}}$ 来表征二者之间的相对滞后程度， $\Delta {S\!{t}}^{\rm{*}}$ 越大表明推力产生越滞后于反卡门涡街形成。图9给出了 $\Delta {St}^{\rm{*}}$ 随俯仰幅度的变化，可以看出， $\Delta {S\!t}^{\rm{*}}$ 随 ${\alpha _0}$ 增大而减小，从而推力产生与反卡门涡街出现的相对滞后程度随俯仰幅度的增大而减小。 $\Delta{S\!t}^{\rm{*}}$ 在0.06～0.14之间，对应俯仰振荡频率范围为0.7～1.7 Hz。

3.2 黏性效应对阻力-推力转变临界点的影响

不可压黏性流动中气动力的主要来源是翼型表面的压力和黏性应力分布，当压力分布不均产生的推力与剪切层的黏滞阻力相平衡时出现零推力点。为分析造成俯仰翼型推力产生滞后于反卡门涡街出现的主要因素，首先利用方程（4）考察翼型表面非定常压力和黏性应力分布对推力产生的贡献。图10给出了俯仰幅度为6°和10°时周期平均推力系数及其分量随斯特劳哈尔数的变化。发现基于方程（4）计算的总推力系数随 $S\!t$ 的变化曲线与CFD软件直接给出的结果几乎完全一致，表明该气动力分解中引入的积分误差很小。

图10显示，在两个俯仰幅度下，翼面压力分布产生的推力分量仅在 $S\!t$ 很小时为负值，为压差阻力，当 $S\!t$ 大于某一临界值时均为正值，为压力诱导推力，定义该临界值为 $S\!{t_{TP{\rm{ = 0}}}}$ 。在整个 $S\!t$ 取值范围内，翼面黏性应力分布产生的推力分量均在零推力线下方，为黏滞阻力，并随 $S\!t$ 增大呈缓慢线性增大的趋势。推力系数随 $S\!t$ 的变化曲线与其压力分量曲线接近，尤其在更大的俯仰幅度下（α₀ = 10°），表明该雷诺数下俯仰翼型的推力主要是由翼面附近的非定常压力场产生的。

为更直观地显示，图11将图10中零推力点附近进行了放大，对黏滞阻力加了负号。图中阴影区为反卡门涡街对应的斯特劳哈尔数范围，其左侧为卡门涡街区域，右侧为不对称尾迹区域。俯仰幅度为6°时，推力的压力分量在卡门涡街-反卡门涡街转变临界点之前就变为正值，但无法克服黏性剪切层产生的阻力；之后在反卡门涡街区域中间靠左部分，压力诱导推力随 $S\!t$ 增大迅速增大，但仍不足以抵消剪切层的黏滞阻力，从而总的推力依旧为负值，表现为阻力；直到 $S\!t$ 增大到临界值 $S\!{t_{T{\rm{ = 0}}}}$ 时，由翼面非定常压力分布贡献的推力恰好与剪切层的黏滞阻力相等，总的推力为零。类似地，α₀ = 10°时总推力及其分量的变化特征与α₀ = 6°一致，此时推力的压力分量恰好在卡门涡街-反卡门涡街转变临界点为零，即 $S\!{t_{TP{\rm{ = 0}}}} \approx {S\!t}^{\rm{*}}$ 。

图8中绘制了推力的压力分量为零时对应的临界斯特劳哈尔数 $S\!{t_{TP{\rm{ = 0}}}}$ 随无量纲幅度的变化，即图中曲线 $\overline {C_T^{{\rm{(PRE)}}}} = \:0$ 。曲线 $\overline {C_T^{{\rm{(PRE)}}}} = \:0$ 左下方为压差阻力区，右上方为压力诱导推力区。可以看出，各俯仰幅度下，压力分量临界点均在卡门涡街-反卡门涡街转变临界点之前或附近出现，表明推力产生滞后于反卡门涡街形成不是由翼面的压差阻力造成的。在曲线 $S\!{t_A}=0.18$ 和 $S\!{t_A}=0.3$ 之间的反卡门涡街区域，压力诱导推力无法克服剪切层的黏滞阻力，是造成零推力点后移的主要原因。从图10和图11中发现，对给定的 $S\!t$ 值，压力诱导推力大小严重依赖于俯仰幅度，且在 ${\alpha _0}$ 增大时大幅增大，而黏滞阻力的大小对 ${\alpha _0}$ 的变化并不敏感，因此俯仰幅度更大时更容易产生足够抵消黏滞阻力的压力诱导推力，从而零推力点更靠近卡门涡街-反卡门涡街转变临界点。这就解释了为什么推力产生相对于反卡门涡街出现的滞后程度随俯仰幅度增大而减小。

对给定的俯仰幅度，在推力为零和其压力分量为零时可以分别唯一确定一个临界斯特劳哈尔数St_T=0和St_TP=0，对应一个数对（St_T=0，St_TP=0）；通过改变俯仰幅度就可以确定一系列数对。图12反映了由这些数对确定的点在St_T=0- St_TP=0参数空间中的分布。发现这些点全部落在了St_T=0- St_TP=0空间中的一条直线上（实线），拟合的相关系数几乎为1。在雷诺数为无穷大或忽略黏性效应影响时，推力全部由翼面压力分布产生，因此St_T=0=St_TP=0，对应图12中的虚线。这表明低雷诺数的黏性效应将无黏时的零推力临界斯特劳哈尔数St_TP=0按照线性比例规律放大为St_T=0，且当俯仰幅度在2°～10°范围内时放大因子（系数）不依赖于俯仰幅度大小。对低雷诺数流动，翼型表面的压力分布可以通过实验测量，而对黏性应力分布的测量比较困难。因此根据St_T=0和St_TP=0的线性依赖性可以直接由压力分量临界点确定总推力的临界点，从而忽略对黏滞阻力的测量。St_T=0和St_TP=0之间的线性相关规律仅是针对NACA0012翼型在1 000雷诺数下作较小幅度（α₀ = 2°～10°）俯仰运动得到的，对不同雷诺数、不同翼型形状及更大的俯仰幅度是否成立还需要进一步探究。

3.3 非定常尾迹对阻力-推力转变临界点的影响

尾迹形态随斯特劳哈尔数的演变是扑翼流场的主要特征，为分析其对推力产生的影响，基于方程（6）定量估计非定常尾迹对推力的贡献。关于该方程中各分量的物理意义及控制体的选取已在第2节说明。图13对比了俯仰幅度为10°时由CFD直接计算的各瞬时气动力系数和基于有限控制体的尾迹流场积分结果。发现对不同的 $St$ 值，两种方法得到的升力系数和力矩系数的时间历程曲线几乎完全一致，阻力系数除在峰值位置有微小差异外也吻合得很好。图14对比了由CFD得到的周期平均推力系数和尾迹流场积分结果。可以发现，俯仰幅度为6°时，两种方法得到的结果基本一致，尾迹积分结果略微超估了推力； ${\alpha _0}$ = 10°时两种方法的计算结果吻合得很好。基于控制体的积分结果对更大的俯仰幅度来说精确度更高；在俯仰幅度较小时，尾迹涡自身强度较弱，在向下游运动时涡量耗散得较快，从而导致了一定的积分误差。

图15对比了俯仰幅度为10°时不同斯特劳哈尔数下各推力系数分量的时间历程曲线，图中推力系数上标V表示涡推力，A和B分别为翼型运动的非定常惯性效应引起的局部流体加速和翼型占据的虚拟流体对推力的贡献，M、P、S分别代表动量诱导推力、压力诱导推力和黏性剪切力。

图15显示，外控制面上的黏性应力分布对推力的贡献（ $C_T^{{\rm{(S)}}}$ ）在 ${10^{ - 4}}$ 量级，可以忽略不计。翼型占据的虚拟流体产生的推力（ $C_T^{(B)}$ ）均位于零推力线下方，呈现为阻力， $C_T^{{\rm{(B)}}}$ 的幅度较小。除 $S\!t\:=\:0.185$ 外，其它 $S\!t$ 值下局部流体加速贡献的推力（ $C_T^{{\rm{(A)}}}$ ）的幅度随 $St$ 增大而增大，但在一个俯仰周期内总是以正弦规律分布，从而其周期平均值为零，主要原因是翼型的上下俯仰运动是对称的； $S\!t=0.185$ 时， $C_T^{{\rm{(A)}}}$ 在一个周期内的波峰或波谷高度明显不一致，其它推力系数分量也出现类似的情况，这可能与该 $St$ 值下的尾迹不对称现象有关（图5）。因此，在零推力点 $S\!{t_{T{\rm{ = 0}}}}=0.143$ （ ${\alpha _0}$ = 10°）附近， $C_T^{{\rm{(A)}}}$ 和 $C_T^{{\rm{(B)}}}$ 对推力的贡献也可以忽略，后边的讨论将进一步对此进行说明。

基于以上分析，非定常尾迹对推力的贡献主要包括三部分，即涡推力（ $C_T^{{\rm{(V)}}}$ ）、动量诱导推力（ $C_T^{{\rm{(M)}}}$ ）和压力诱导推力（ $C_T^{{\rm{(P)}}}$ ）。图16分别给出了相应的时均流场，其中速度场为沿来流方向的速度分量；图中分别对速度场和压力场以来流速度和局部最大压力（主流的压力）为参考进行了无量纲处理。图16（a）中蓝色和红色区域分别是顺时针涡和逆时针涡的时间平均结果，与图5中的瞬时涡量场对应。图16（b）中无量纲速度小于1的区域为速度损失区（蓝色），大于1的区域为速度增加或“喷流”区（红色）。结合图5，发现卡门涡街导致尾迹流场速度损失( $S\!t=\:0.041,\; $ $ 0.066$ )，反卡门涡街使尾迹流场速度增加或产生“喷流”现象(St = 0.103～0.164)，介于二者之间的临界状态( $S\!t= 0.082$ )对尾迹速度场的影响很小。图16（c）中无量纲压力小于1的区域为压力损失区（蓝色），等于1的区域为主流区（红色）。结合图5，发现卡门涡街和反卡门涡街尾迹流场都出现压力损失现象。

图15显示，St = 0.041～0.164范围内， $C_T^{{\rm{(V)}}}$ 在一个俯仰周期的大部分时间内小于零，为涡致阻力；结合图5和图16（a），发现在零推力点 $S\!{t_{T{\rm{ = 0}}}}=0.143$ （ ${\alpha _0}\:{\rm{ = }}\:$ 10°）附近，卡门涡和反卡门涡自身都产生涡致阻力。St = 0.082～0.164范围内， $C_T^{{\rm{(M)}}}$ 均在零推力线上方，为动量诱导推力，其幅度随 $S\!t$ 增大逐渐增大并更加远离零推力线。因此动量诱导推力随 $S\!t$ 增大而增大，这可以通过图16（b）中反卡门涡街的“喷流”强度随 $S\!t$ 增大而增强解释。 $S\!t\;{\rm{ = }}\;0.041$ 时， $C_T^{{\rm{(M)}}}$ 均为负值，是由卡门涡街尾迹速度损失产生的动量诱导阻力； $S\!t\;{\rm{ = }}\;0.066$ 时， $C_T^{{\rm{(M)}}}$ 在一个周期内分布在零推力线上下两侧，之后的周期平均结果显示其对应小量的动量诱导推力。总体上， $C_T^{{\rm{(M)}}}$ 随斯特劳哈尔数的变化趋势与图5中尾迹形态变化和图16（b）中尾迹速度损失或增加的变化特征一致。 $C_T^{{\rm{(P)}}}$ 在一个俯仰周期的绝大部分时间内小于零，为尾迹流场压力损失引入的附加阻力； $C_T^{{\rm{(P)}}}$ 的幅度随 $S\!t$ 的增大而增大，原因是 $S\!t$ 更大时，尾迹压力损失程度更深，因此其导致的附加阻力更大，如图16（c）。

综上分析，在零推力点 $S\!{t_{T{\rm{ = 0}}}}\:{\rm{ = }}\:0.143$ （ ${\alpha _0}$ = 10°）附近，卡门涡和反卡门涡自身的涡致阻力及尾迹流场压力损失引入的附加阻力是俯仰翼型阻力的主要来源，此外卡门涡街导致的尾迹速度损失也贡献小部分阻力；反卡门涡街的“喷流效应”产生的动量诱导推力是推力的主要来源。该结论是针对俯仰幅度为10°得到的，对其它俯仰幅度也成立。为达到零推力点，涡致阻力、压力损失的附加阻力和动量诱导推力需达到平衡。

图17显示了 ${\alpha _0}$ = 6°～10°范围内各周期平均推力系数分量随斯特劳哈尔数的变化，其中 $\overline {C_T^{{\rm{(Eq6)}}}} $ 为基于尾迹积分得到的总推力系数；图中的阴影区域为反卡门涡街对应的 $S\!t$ 范围；箭头标注的数字为基于尾迹分析预测的零推力点（左侧）与CFD给出的结果（右侧）对比，二者存在小量差异，表明尾迹积分对推力临界点的预测存在一定误差。对不同俯仰幅度，在零推力点附近，外控制面上的黏性应力分布贡献的周期平均推力（ $\overline {C_T^{{\rm{(S)}}}} $ ）几乎为零，翼型运动惯性效应和翼型占据的虚拟流体对周期平均推力的贡献也很小可以忽略不计（ $\overline {C_T^{{\rm{(A)}}}} $ 和 $\overline {C_T^{{\rm{(B)}}}} $ ），这与前面的分析一致。

图17（c）显示，对α = 10°，St < 0.164时， $\overline {C_T^{{\rm{(V)}}}} $ 均小于零并随St的变化程度很小，为较稳定的涡致阻力； $\overline {C_T^{{\rm{(M)}}}} $ 在St很小时（0.041）为负值，为动量诱导阻力，之后在St增大时变为正值并大幅增大，为较大的动量诱导推力。St > 0.164时， $\overline {C_T^{{\rm{(V)}}}} $ 急剧增大并最终变为涡推力；动量诱导推力大小基本保持不变。在图中显示的整个St范围内， $\overline {C_T^{{\rm{(P)}}}} $ 均为负值，且随St增大而减小，表明由尾迹压力损失引入的附加阻力随St的增大而增大。在卡门涡街-反卡门涡街转变临界点（St^* = 0.082），即图中阴影区左侧边界位置，小量的动量诱导推力无法克服涡致阻力和尾迹压力损失引入的附加阻力，总推力为负值，表现为阻力；之后进入反卡门涡街区域，动量诱导推力随St增大逐渐增大，但仍不足以抵消涡致阻力和尾迹压力损失的附加阻力，直到 $\overline {C_T^{{\rm{(M)}}}} $ 与 $\overline {C_T^{{\rm{(V)}}}} $ 和 $\overline {C_T^{{\rm{(P)}}}} $ 相抵消时达到零推力点。因此，反卡门涡街的“喷流效应”产生的动量诱导推力无法克服反卡门涡自身的涡致阻力和尾迹流场压力损失引入的附加阻力，是俯仰翼型推力产生滞后于反卡门涡街出现的主要原因。该结论对其它俯仰幅度也成立，如图17（a）和图17（b）。α₀ = 6°时尾迹流场压力损失引入的附加阻力随St增大先略有减小后急剧增大；α₀ = 8°时各周期平均推力系数分量随St的变化规律与α₀ = 10°时基本一致。

4 结　论

本文针对在低雷诺数下( $Re = 1\;000$ )，随着斯特劳哈尔数增大，俯仰振荡翼型推力产生滞后于反卡门涡街出现的现象及其背后的物理机制进行了分析，主要结论如下：

1） 俯仰幅度在2°～10°范围内时，随 St增大，卡门涡街-反卡门涡街转变临界点与阻力-推力转变临界点均不重合，且前者先于后者出现；二者之间的相对滞后程度随俯仰幅度的增大而减小。

2） 俯仰翼型的推力主要由翼型表面压力分布产生，黏性应力分布对推力大小的影响有限。不同俯仰幅度下，翼面压力分布积分得到的推力分量在卡门涡街-反卡门涡街转变临界点之前或附近就由负值变为正值，因此零推力点后移不是由翼面的压差阻力造成的。

3） 当翼型振荡参数进入对应反卡门涡街尾迹形态区域时，翼面压力分布产生的推力分量无法克服剪切层的黏滞阻力，是导致推力产生滞后于反卡门涡街出现的主要原因。

4） 基于控制体的气动力分析表明，尽管反卡门涡街会产生“喷流效应”，但在 $S\!t$ 较小时，其产生的动量诱导推力无法克服反卡门涡自身的涡致阻力和尾迹流场压力损失引入的附加阻力，因此即使存在反卡门涡街也不能产生净推力。

本文系统地讨论了俯仰翼型推力产生与反卡门涡街形成的相对滞后特征，较深入地探究了该滞后现象的发生机制，但也存在其局限性，如暂未考虑来流速度及雷诺数的影响。此外，翼型形状、俯仰轴位置的影响还需进一步探究。尾迹形态演变是扑翼流场的主要特征，深入认知反卡门涡街与推力产生的关系有利于揭示扑翼推力产生机理，在扑翼式MAV设计方面具有潜在的应用价值。