大幅度非定常运动流场分析的新方法及其在柔性扑翼上的应用

图２果蝇翼半个扑动周期（下扑过程）的几个瞬时位置
扑动翼在半个扑动周期中（上扑或下扑），在中间阶段一般保持迎角ｎ和扑动速度不 变．转动角速度为０；在扑动开始阶段经历快加速和快翻转；在扑动结束阶段则经历快减 速和快翻转。参考文献［８］有详细运动规律的数学表达式，在此不再赘述。
探索创新
交流（第三集）——第三届中国航空学会青年科技论坛文集
１９１
根据Ｄｉｃｋｉｎｓｏｎ的试验¨纠和孙茂的数值计算一Ｊ，翼的扑动和转动运动之间的相位差是 影响其扑动翼气动力特性的一个关键因素。据此将悬停飞果蝇翼的扑动模式分为三类：超 前模式、对称模式和滞后模式。定义如下：当转动运动和扑动加减速运动相位同步时称为 对称模式，当转动运动相位提前于扑动加减速运动时称为超前模式，反之称为滞后模式。 本文计算采用的参数为：ｐ＝１．２２５ｋｇ／ｍ３，扑动翼长Ｒ＝２．５ｍｍ，平均弦长ｃ＝ ０．８７９ｍｍ，扑动频率２７２Ｈｚ，参考速度Ｕ＝２．２ｍ／ｓ。其雷诺数１３３，和试验值相近。扑动 幅度、转动角幅度以及上下扑动过程中间阶段的迎角分别为咖＝１６０。，Ａａ＝１００。，Ｏｔ＝ ４０。。翻转的时刻和翻转的持续时间分别为△Ｌ＝０．ＩＴ，Ａｒ，＝０．３２Ｔ。超前模式的转动运 动相位比挥动运动相位超前８％，对称模式两者一致，滞后模式则落后８％。图３为３种
对有运动边界的问题，计算网格需随着边界做适当的调整。使用动态网格时，为避免 网格运动引入额外的计箅误差，控制体表面面积、体积、网格运动速度等几何参数的计算
须满足几何守恒定律
戋￡ｄ以一屯％。出＝ｏ．
４
（６）
非结构动态变形嵌套网格方法
本文将双重Ｄｅｌａｕｎａｙ图映射的网格变形技术¨１’’２ ３和非结构嵌套网格方法结合起来， 用于解决既有整体大幅运动，又有局部小变形的非定常流场问题。飞鸟、昆虫等动物扑翼 飞行时，翅膀上下大幅扑动的同时又有弦向和展向的柔性变形，翅膀的动态柔性变形可能
法［Ｉ旬】很好地解决了这一问题。它通过对可压缩流控制方程的时问导数项的预处理，使方
程的系数矩阵的特征值保持在同一量级而不至于相差太大，使传统的可压缩流方法具备了 求解低速流场的能力。但日前，该方法主要用于定常流场的计算。 飞鸟、昆虫等扑动飞行时，翅膀既有大幅度的扑动又有局部的柔性变形１４．５ Ｊ。对含有
探索创新交流（第三集）——第三届中国航空学会青年科技论坛文粟
『｛ｍ≤旧
ａ
ｍ
Baidu Nhomakorabea鬻鬻鋈纛 貉骶鬯
探索创新交流（第三集）——第三届中国航空学会青年科技论坛文集
１９３
常有效。图６（ｃ）是一个周期内柔性翼和剐性翼产生的升力阻力系数的比较，可以看出， 在大部分时间内，柔性翼产生的升力要稍大于、而阻力又稍小于刚性翼．说明翼的柔性有
待。本文综合参考文献［１３］总结的针对定态网格的当地最大预处理和当地并全局截断
预处理的两种方法，提出对ＡＬＥ格式的预处理参数的选取办法。 （４）
其中：Ｍａ。。为当地最大的相对马赫数，定义为
探索创新交流（第三集）——第三届中国航空学会青年科技论坛文集
肘ｎ眦ｌ：ｍａｘ
ｆ卫啦，卫上止１
、 Ｃ Ｃ ，
１８９
５计算结果与分析
采用ｉ维果蝇翼和二维柔性蜻蜒翼的扑动流场作为本文方法的应用实例，考察算法的
效率和精度。 ５．１三维果蝇翼的扑动流场
该算例采用的果蝇翼平面形状和Ｄｉｃｋｓｏｎ等人””的动态比例模型形状相同（如图１ （ａ）所示）。等厚度平板，边缘光滑过渡，平板厚度为弦长的５％。模型为刚性翼，但流 场的计算过程仍采用嵌套阿格结合网格变形的方法来实现。 计算区域为整个流场的一半，嵌套网格分为３层，如图Ｉ（ｂ）和图Ｉ（ｃ）所示。里层 为贴近扑动冀的非结构混合阿格，冀表面划分为３２１８４个三角形．黏性层内为拉伸比大的
Ｐｐ’０
ｕｐｐ’Ｐ
０ ０
Ｐ
０ ０ ０
Ｐ ｐｗ
ＰＴ ｕｐｒ ｖｐｒ ｗｐｒ
Ｆ＝
ｖｐＰ’０ ｗｐｐ’０
ｐｐ
（３）
０
ｐｖ
式中：对理想气体，以’定义为
Ｐｐ．，＝萨１一掣＝毒∥＝南，．ｃ为声速。
卢＝ｍｉｎ（ｍａｘ（讹２一，ｐ。沁），１．０）
７Ｈｏ一（１一ｐｈＰ）ｐｕ
ＰｒＨｏ＋ｐ＾ｒ
此处，对ＡＬＥ格式，预处理参数ｐ的选取不同于定态网格下的方程求解，需小心对
利于气动性能的改善。
圈
。，嗽．…。……，，儡
ｆ１１黍性■
。，魄…｛…“……濡
ｍ剐性一
ｆ０时间壹他的升阻力量教
田６二维螭蜒翼模型柔性翼和刚性翼的计算结果比较
６结论
本文发展了一种用可压缩流控制方程模拟扑翼低雷诺数非定常流场的数值方法。通过
律，并和Ｄｉｃｋｉｎｓｏｎ等¨５１的动态比例模型试验结果的比较。由图４可见，计算出的升力系
数随时间变化的趋势和试验结果吻合得比较好，在部分区域，两者在量值上有较小的差 别；实际上，Ｄｉｃｋｉｎｓｏｎ等人的试验中扑动和转动角度有±４．５０的系统误差，因此，在比 较计算和试验结果时，不必也不可能做到两者在量值上完全相同。 通过本算例验证表明，本文发展的数值方法对三维低雷诺数扑翼流场计算有效，计算 结果合理可靠。
２方程离散与求解
式（１）在引入预处理的伪时间导数项后变成
ｆ旦＿０
７＂Ｌｏ Ｑｄｙ＋生Ｏ Ｌ。）Ｗｄｙ＋‰。）（Ｆ（形）一％形）ｄ５ ２虹。）ＦｖｄＳ
ｔ
（２）
式中：ｒ、ｔ分别为伪时间和物理时间；Ｆ为伪时间导数项的预处理矩阵；Ｑ为压力、速度
和温度的原参变量，Ｑ＝（Ｐ，ｕ，ｔ，，ｔ‘７，ｒ）１。。注意到，当ｒ一∞时，式（２）左边的第１
不同运动模态，扑动角沙和转动角Ｏｔ在一个扑动周期中随时间的变化，以及无量纲扑动
速度ｔ，Ｉ＋和转动角速度ａｌ＋随时间的变化。
ｔ／ｒ
ｔ／Ｔ
（ａ）扑叻和转动的角度
（”扑动和转动的角度
图３果蝇翼３种不同运动模态的扑动规律
计算中，每个扑动周期分成２００个物理时间步，每时间步用１５步内迭代推进，各种 模态各计算了４个周期。图４是果蝇翼在３种不同扑动模态下的升力系数随时间的变化规
是影响高升力产生的重要冈素Ｈｊ Ｊ。计算这类扑翼的非定常流场，网格变形技术和嵌套网 格方法结合使用，尤为必要也非常有效。双重Ｄｅｌａｕｎａｙ图映射方法和非结构嵌套网格方 法结合使用的步骤如下： （１）将流场划分为多个重叠的计算域并各自生成独立的网格； （２）为每套含有运动壁面边界的网格生成Ｄｅｌａｕｎａｙ背景罔，并建立网格点和背景图 单元一一对应的映射关系； （３）根据壁面边界的整体运动轨迹，将整个背景图刚性移到新的位置； （４）根据壁面边界的主动变形规律或者流同耦合计算出的壁面变形量，将背景图的
自然界中，飞鸟、昆虫等飞行生物无一例外采用扑翼方式飞行，其特点是将举升、悬
停和推进功能集于一个扑翼系统，可以用很小的能量进行长距离飞行，同时具有较强的机
动性。关于扑翼飞行原理的深入研究将有助于提高现有飞行器的性能，推动微型飞行器的 开发和研制。对扑翼的低雷诺数非定常流场的数值研究是其巾的重要内容。 传统的可压缩流计算方法在中等亚声速到超声速的流场模拟中应用广泛，但众所周 知，在求解低速流场时会面临收敛速度慢、精度低的所谓“刚性”问题；预处理方
理方法引入到Ａｒｂｉｔｒａｒｙ—Ｌａｇｒａｎｇｉａｎ．Ｅｕｌｅｒｉａｎ（ＡＬＥ）格式的非定常可压缩Ｎ—Ｓ方程中，在 非结构的动态混和网格基础上求解有运动边界的低速流场。控制方程在物理时间方向上用
预处理的伪时间步推进，空间上用二阶精度迎风格式有限体积法离散。离散后的方程用隐 式、ＬＵ—ＳＧＳ迭代求解。双重Ｄｅｌａｕｎａｙ图映射的网格变形技术ｕｈｌ２１和非结构嵌套网格方 法结合使用，解决扑翼既有整体大幅运动，又有局部柔性变形的复杂运动边界问题。最后 作为应用，选取两个典型的扑翼算例验证本文发展的算法。
项消失，方程恢复到式（１）的形式。双时间步推进的方法在每个物理时间步的外循环中
包括伪时间步的内循环，亦即伪时间步将两个物理时间之间的流场当做定常问题求解。因
此，物理时间步不受方程系统矩阵的刚性影响，而内迭代的无需时间精确的伪时间步可以 用预处理、当地时间步长等加速收敛办法来求解。 本文采用参考文献［３］的预处理矩阵，定义ｆ为
后向隐式欧拉方法离散，两物理时间步之间用预处理的伪时间步推进求解，并采用当地时 间步长等加速收敛的方法。由于预处理改变了系统矩阵的特征系统，所有的黎曼不变量的 边界条件不再适用，因此，所有基于特征量的边界条件都必须做相应的改变以适应预处理
后的特征系统。详细计算方法可参考参考文献［１４］。
３几何守恒定律
１８８
探索创新
交流（第三集）——第三届中国航空学会青年科技论坛文集
１非定常流体控制方程
积分形式ＡＬＥ格式的可压缩Ｎ—Ｓ方程为
３
， ， ，
｝Ｉ ＷｄＶ＋牵 （Ｆ（形）一口印Ｗ）ｄＳ＝巾 ｏ‘ｍ“） ＪａｆＸＩ） Ｊｍ（Ｉ）
Ｆ，ｄＳ
（１）
式中：力（ｔ）为随时间变化的控制体；守恒量Ｗ＝（Ｐ，ｐｕ，ｐｖ，ｐｗ，ＰＣ．）７；Ｆ（Ｗ）、Ｆ， 分别为无黏通量和黏性通量；移。为控制体表面的法向运动速度。如上文所述，由于传统 的可压缩流方法在计算低速流场时面临“刚性”问题，因此，本文在式（１）中引入预处 理的伪时间导数项，用双时间步推进求解。
５．２二维柔性蜻蜓翼的扑动
另一应用实例是二维蜻蜒扑翼。该模型有一前一后两个扑翼，并考虑翼的弦向柔性变 形，问题较为复杂。本文仅计算该模型在悬停状态时的流场，其扑动规律为，前后两翼在 ６００倾斜的平面内以１８００相位差上下扑动和弦向俯仰运动的同时¨６ｌ，按照参考文献［５］ 测量出的规律做主动的弦向动态弯度变形，变形的最大弯度为弦长的１６％。 图５是二维蜻蜒翼扑动过程中几个时刻的网格图。不难看出，两方法的结合使用可有 效地处理这类既有整体大幅运动又有局部小变形的流场问题，省去了网格再生的繁琐工 作，保证了网格质量且效率高。
（５）
卢。ｉ。取为卢呐。＝３ 赫数。
ｍａｘ（讹：，Ｍａ］．。。），Ｍａ。为来流－ｑ赫数，Ｍａ。．。。为边界运动的平均马
经过上述预处理，方程系统矩阵的所有特征值大小保持在同一量级，使得预处理后的 方程系统在所有速度范围内都有良好的性态。 式（２）的求解，在空间上采用二阶迎风格式的有限体积法离散，物理时间采用二阶
壁面边界点移至新位置；
（５）根据第（２）步建立的映射关系，将网格点一一移动到新的位置； （６）建立各套网格的嵌套关系进行流场计算；
探索创新交流（第三集）——第三届中国航空学会青年科技论坛文集
（７）下一物理时间步，重复执行（３）一（５）步。 两种方法的结合，省去了同格再生的繁琐工作，保证了网格质量，而且效率较高。
运动边界的非定常流场计算，每时间步都需要调整网格来适应计算域的变化。网格动态变 形技术多用于较小变形或较小运动幅度问题，对较大相对幅度运动的问题则需要网格再生
来保证网格的质量帕１；嵌套网格技术处理大幅边界运动的问题非常有效，广泛用于战机 投弹等非定常流场的计算。７棚Ｊ，但目前仅限于固壁的刚性运动，如果壁面有变形，也需要 重新生成网格来适应新的壁面边界¨…。 鉴于以上问题，本文发展一种适于扑翼低雷诺数非定常流场的数值计算方法。将预处
探索创新
交流（第三集）——第三届中国航空学会青年科技论坛文集
１８７
大幅度非定常运动流场分析的新方法 及其在柔性扑翼上的应用
肖天航 昂海松周新春
（南京航空航天大学微型飞行器研究中心，江苏南京２１００１６）
摘要：将预处理引入到ＡＬＥ格式的非定常可压缩Ｎ—Ｓ方程中．发展了一一种能有效模拟大幅度边界运动 的非定常流场的数值方法。控制方程在非结构的动态网格上用双时间步推进并用隐式、ＬＵ—ＳＧＳ迭代求 解；双重Ｄｅｌａｕｎａｙ圈映射的网格变形技术和非结构嵌套网格方法结合使用。解决既有整体大幅运动，又有 局部柔性变形的复杂运动边界问题。扑翼应用实例的结果表明，该方法将可压缩流方法的计算能力拓展到 了低码赫数范围．用于计算大幅度非定常流场能得到满意的计算精度和效率。 关键词：非定常流；预处理；双时间步；嵌套网格；网格变形；扑翼；低雷诺数
婚霄蘸 懈
定常流场。往外为逐渐变稀的第３层网格。
（ａ）果蝇置模型 Ｏ）嵌套弼格 （ｃ）嵌套网格局部放大圈 图１
＝棱柱，往外为四面体，共ｌｏｌ万网格单元。第２层为较密的笛卡尔网格以准确地捕获非
三维果蝇翼模型和计算网格
图２是半个扑动周期（下扑过程）几个等时间间距的果蝇翼的瞬时位置。坐标轴固
定在翼根处且：轴垂直扑动平面。冀绕ｚ轴的方位角转动称为“扑动”，扑动角记为廿； 翼绕Ｙ轴的翻转称为“转动”，转动角记为“，卫称为迎角。