三维可变形翼展向变形气动力特性数值模拟研究

陈 默，张 楠，郑文涛，李永成

（中国船舶科学研究中心，江苏 无锡 214082）

0 引 言

随着现代航空技术的飞速发展，传统飞行器的设计已经不能满足人类对飞行性能和飞行品质不断提升的需求。利用仿鸟类及昆虫飞行时翅膀展向变形现象而发明的可变翼飞行器一直是各国航空发展领域的研究热点[1]。可变翼飞行器甚至可利用机翼的扭转和弯曲辅助替代操纵面，实现飞行器横滚、偏航以及俯仰等不同姿态变化。

自上世纪80 年代开始，国外就开始了对变形机翼的研究，基于洛克韦尔公司提出的可充分利用飞行器在飞行时产生的扭转变形概念，美国NASA 兰利研究中心联合美国空军及洛克韦尔公司开展了主动柔性机翼（active flexible wing，简称AFW）项目，旨在研究主动柔性机翼技术[2]，且在工程设计中得到初步应用。到1996 年，在“主动柔性机翼”计划研究基础上，美国提出了利用机翼扭转使得结构载荷降低、飞行包线提高的主动气动弹性机翼（active aeroelastic wing，简称AAW）计划[3]。

2006 年，欧洲主动气动弹性飞行器（3AS，active aeroelastic aircraft structures）计划提出了以“主动内部机构”概念来实现机翼结构的主动气动弹性控制，该结构同样可以实现对机翼展向变形的控制[4]；2017 年，葡萄牙贝拉国内大学的Santos 等设计了可变展长机翼结构，采用伸缩方式，由齿轮齿条机构带动外部机翼从做成空壳状的内部机翼中伸出，机翼的最大展长可达2.5 m，且装配在无人机上进行了飞行测试[5]。

国内，燕山大学的Guo 等[6]于2016 年提出刚性滑动蒙皮的设计关键在于机翼表面蒙皮划分的概念，并相应研究了两种有关机翼表面蒙皮划分的方法；哈尔滨工业大学的王礼佳提出了一种可变翼展的伸缩翼，机翼采取多级伸缩的方式，在内部布置剪叉机构来驱动机翼进行运动，除伸缩外，机构还能够通过形状记忆合金丝实现后掠角的变形[7]；南京航空航天大学的宋哲[8]利用形状记忆合金作为驱动器，设计了一种能够改变厚度的机翼，其分析结果显示，厚度增加后的机翼阻力系数降低，升力系数和升阻比均得到了提高。

由于计算流体力学中关于建立复杂变形几何结构的流场数值模拟技术逐渐趋于成熟，基于数值模拟方法研究变形翼的气动力特性这一可行途径已得到广泛采用。Wang等[9]利用数值模拟方法研究了无人机机翼伸缩对其气动力特性的影响，结果表明机翼伸缩对升/阻力系数影响不大，但对升阻比影响显著；程春晓等[10]采用可压缩Navier-Stokes 方程和Spalart-Allmaras(S-A)湍流模型，对柔性后缘可连续变弯度二维NACA 0012 翼型升/阻力系数、表面压力分布和流场结构展开了深入的研究，基于该变形方式下的翼型，其升力系数和升阻比均得到较大改善；马洋等[11]采用RANS 方法，结合SSTk-ω两方程湍流模型对伸缩/后掠变构型机翼在亚声速和超声速阶段的气动力特性开展了研究，相对于固定机翼外形，该构型升阻比优势明显。

此外，机动性一直是水中航行体性能研究和设计优化的核心命题之一。舵作为主要操纵装置，其水动力性能（舵效）的优劣直接影响船舶和水下航行体的活动范围和机动性能，随着船舶和水下航行体尺度越来越大，高舵效舵的需求尤其迫切。对于三维可变形机翼的初步研究，可为新型翼舵及新型操纵控制技术的研发和水下航行器的优化设计提供技术支撑。

综上，国内外学者已针对三维可变形翼的气动力特性开展了研究。但目前绝大多数变形翼方案均基于改变机翼展长、后掠角和弦长来实现气动性能的改变，对于机翼展向变形的研究内容较为稀缺。因此，以展向变形作为出发点来开展三维可变形翼气动力特性的研究具有重要的研究意义。有鉴于此，本文对三维可变形翼定常状态下给定静态展向变形和非定常状态下给定动态展向变形气动力特性及变化规律进行了数值模拟研究。

1 研究对象

1.1 静态展向变形模型

静态变形计算对象为三维NACA 0012 机翼（展弦比b/c=0.5 m/0.2 m=2.5，其中b和c分别表示机翼的展长和弦长），及以其为基础的一系列给定静态展向变形的可变形翼（即展向弯曲的刚性机翼），如图1和图2所示。静态展向变形规律按轴线变形表达式（1）来确定，变形起始位置bs为0 mm，变形幅度ha分别为50 mm、100 mm和150 mm。攻角设为α=-18°～18°，间隔2°（抬头为正），局部角度加密。

1.2 动态展向变形模型

动态变形计算对象基础原型与静态展向变形一致，变形规律按轴线变形表达式（2）来确定，变形幅度ha分别为100 mm和150 mm，变形频率f分别为0.5 Hz、1 Hz、2 Hz、3 Hz和4 Hz。攻角α设为5°。

一个运动周期内动态展向变形翼的运动过程如图3所示，其中T表示运动周期。

2 数值计算方法

2.1 控制方程

不可压缩流体连续性方程与RANS方程的张量形式为

2.2 湍流模型

本文采用k-ε湍流模型。湍流动能k和耗散率ε方程为

式中，C1ε= 1.44，C2ε= 1.92，Cμ= 0.09，σk= 1.0，σε= 1.3，Gk、Gb为湍流动能生成项，σk、σε表示湍流Prandtl数，Sk、Sε为源项[12]。

2.3 网格划分与边界条件

采用多块结构化网格划分方式，计算区域网格的划分如图4所示，机翼表面周围网格如图5所示。对模型表面附近网格加密，且第一层网格间距根据y+来确定（y+平均值约为1）。动态变形过程利用径向基函数（radio base function，RBF）进行变形区域网格节点插值置换来实现，即根据边界位移定义控制点位移插值场，并通过该插值场将网格节点平移到新位置，使网格可以进行非刚性变形。RBF变形基本原理如下[13]：

每个控制节点i的已知位移可展开为

边界条件定义如下：

速度入口：机翼前缘向前3c，设定来流速度的大小与方向；

压力出口：机翼尾缘向后6c，设定相对于参考压力点的流体静压值；

壁面：机翼外表面，设定无滑移条件；

外场：距机翼表面3c～5c，速度为未扰动的主流区速度。

2.4 数值求解

控制方程采用有限体积法离散，动量方程中的对流项采用二阶迎风差分格式，扩散项采用中心差分格式，湍流模型方程采用二阶迎风差分格式，时间项采用二阶隐式格式。压力速度耦合方法选用流场计算中经典的SIMPLE算法。

3 数值计算方法验证

3.1 网格收敛性验证

由于网格数量在CFD数值计算中扮演极其重要的角色，本文首先开展网格收敛性验证工作，以确保数值计算的可靠性。以无变形三维NACA 0012 机翼为计算对象，选取五套网格划分方案，依照第22 届ITTC 推荐规程7.5-03-01-02，网格三向加细比均为 2，网格数量分别为7.53 万、20.55 万，55.58万、169.83万和442.12万，对应方案代号设为G1、G2、G3、G4和G5。不同数量网格第一层高度相同，y+均为1。五套网格方案不同攻角下，三维NACA 0012机翼升/阻力系数变化曲线如图6所示。其中，来流速度设为30 m/s，流体介质为空气，特征雷诺数Re=4.11×105。升力系数CL和阻力系数CD的定义如式（9）所示，其中，D和L分别表示模型受到的阻力和升力，ρ表示空气密度，本文取为1.225 kg/m3，V为来流速度，A为翼舵平面面积，本文A=0.1 m2。

从图中可以看出，随着网格数的增加，计算结果趋于稳定。当网格数增大到第四套网格时，不同攻角下的升/阻力系数近乎不变化。综合计算精度和计算资源考虑，网格方案G4较为适合本文的计算需求，在下文三维机翼静态和动态展向变形计算中，网格划分方案均选用G4。

3.2 计算精度验证

现存文献中关于NACA 机翼流体动力的试验数据一般针对二维（无限展长）情况，本文选用的是三维NACA 0012 机翼，缺乏公开的试验数据对计算结果进行验证。因此本文首先考虑雷诺数相同时展弦比变化对机翼数值计算结果的影响，选用的展弦比分别为1（b/c=0.2 m/0.2 m）、2.5（b/c=0.5 m/0.2 m）及无限展长（二维，c=0.2 m），特征雷诺数Re=4.11×105，网格划分方案一致。

图7 给出了不同展弦比下NACA 0012 机翼升力系数随攻角变化的对比曲线。由图中可以看出，二维翼型计算结果与试验数据吻合较好，且失速点预报结果一致，从而验证了本文网格划分方案和物理建模方法的适用性。此外，展弦比越小，失速角越往后移，对应的最大升力系数也越小，这一趋势与经典理论吻合。

同时考虑雷诺数变化对三维机翼升力系数数值计算结果的影响，选用的特征雷诺数Re分别为4.11×105、1×106和2×106，计算结果如图8 所示。从图中可以看出，在本文后续计算将采用的特征雷诺数Re=4.11×105基础上，继续增大雷诺数对升力系数计算结果几乎没有影响，且各雷诺数对应的失速角基本一致。

对于三维机翼气动力计算精度问题，本文选用某舵的舵轴力矩实验结果进行计算精度验证，该模型试验是在中国船舶科学研究中心低速风洞中完成的。其中，舵展弦比b/c=1.595，来流风速为45 m/s，舵角δb=-25°～25°，每间隔5°变一个舵角，特征雷诺数Re=3.42×105。图9 给出了舵法向力系数随舵角变化的对比曲线。从图中可以看出，数值计算结果与试验数据吻合较好，且失速区间一致，从而验证了本文网格划分方案和物理建模方法对三维机翼计算结果的可信度。

3.3 时间步长验证

考虑到时间步长在非定常计算中起到极其重要的作用，在进行三维机翼动态展向变形研究前开展时间步长选取验证工作。变形幅度ha为100 mm，变形频率f=4 Hz，即运动周期T=0.25 s，来流速度设为30 m/s,流体介质为空气。分别选取时间步长Δt= 0.01 s、0.005 s、0.0025 s 和0.001 25 s 进行计算验证。计算得到的变形翼升/阻力系数随时间变化曲线如图10 所示。从图中可以看出，随着时间步长的减小，计算结果趋于稳定。当时间步长Δt≤0.005 s时，除短暂的初始扰动时间段外，升/阻力系数随时间变化关系几乎一致。基于计算精度和计算资源的综合考虑，在下文动态展向变形计算中，时间步长Δt均选用0.005 s。

4 计算结果及分析

4.1 静态展向变形幅度的影响分析

（1）对于阻力系数，在失速发生之前，负攻角时，不同展向变形幅度下的阻力系数相对于原始无变形三维机翼均略有增大，且相对增幅均在2.8%以下；正攻角时，不同展向变形幅度下的阻力系数均小于原始无变形三维机翼，且下降幅度随着变形幅度的增大而增大，最大相对降幅达6.0%。因此从阻力产生的能耗考虑，正攻角下展向变形幅度增加对其气动力性能有利。

（2）对于升力系数，其增减趋势几乎与阻力系数一致，且增减幅度更加明显，负攻角下相对增幅最高可达6.3%，正攻角下最大相对降幅达9.8%。从图13中横剖面流线和速度梯度范围分布对比结果来看，正负攻角下升/阻力系数呈现略微不对称的现象是由变形翼的三维效果所引起。从提高升力的角度考虑，负攻角下展向变形幅度增加对其气动力性能有利。此外，不同展向变形幅度下的机翼失速点一致，均在±15o附近；结合图14 和图15 可知，展向变形的存在对于机翼迎流面压力没有明显影响，主要是改变背流面负压值，负攻角下展向变形增加时背流面压力值会更低。

（3）对于升阻比，不同展向变形幅度下最大升阻比均出现在±6°附近，与原始无变形三维机翼一致，大小略有差异。从升阻比角度考虑，正负攻角下展向变形幅度为50 mm时其气动力性能均最为优越。

（4）总体而言，静态展向变形对三维机翼气动力性能影响不大。

4.2 动态展向变形幅度的影响分析

对三维NACA 0012可变形翼模型以不同变形频率和不同变形幅度做动态展向变形运动的气动力特性进行数值模拟，表1 给出了不同变形翼组合的阻力系数和升力系数平均值的比较。图16 给出了变形幅度一定时不同变形频率下升/阻力系数时历变化过程。由以上结果可得给定动态展向变形基本规律：

表1 动态展向可变形翼阻力系数和升力系数平均值比较Tab.1 Comparison of average drag coefficient and lift coefficient of dynamically spanwise deformable airfoil

（1）展向动态变形的存在对阻力系数和升力系数平均值略有影响，阻力系数平均值均随变形频率的增大而微幅下降，而对于升力系数，展向变形幅度ha=100 mm 时其平均值随展向变形频率的增大而微幅上升，展向变形幅度ha=150 mm 时则无明显规律可循。因此从兼顾阻力能耗和提高升力角度来考虑，以小幅度高频做动态展向变形时气动力性能会略微得到改善。

（2）不同变形频率下的升力系数和阻力系数时历变化规律几乎一致，呈类余弦曲线形式，且阻力系数谷值点均正对应于升力系数峰值点，均出现在变形翼经平衡位置向下做展向变形运动的瞬间。相对应地，变形翼经平衡位置向上做展向变形运动时升力系数达到谷值点，阻力系数峰值点则出现在变形翼向上到达最大变形位置时。此外，升力系数和阻力系数变化幅度均随变形频率的增大而增大，且变形幅度越大，变化幅度越明显。因此在变形翼设计阶段，可考虑以不同变形频率和不同变形幅度来规划一个展向变形周期内的运动，以达到气动力性能最优化的结果，例如可使变形翼在由下到上半个周期内低频变形，而在由上到下半个周期内高频变形，上下半程变形幅度亦可做相应调整。

图17 和图18 分别为变形幅度ha=150 mm 时，截取的纵剖面和横剖面一个展向变形周期内的速度云图。从图中可以看出，变形翼在经平衡位置向下变形到向上变形回至平衡位置半个运动周期内，机翼迎流面附近速度逐渐增大，背流面附近速度逐渐减小，对应的升力系数随之逐渐减小；同时机翼前缘驻点附近低速区范围逐渐减小，对应的阻力系数随之逐渐增大；往后半周期内则正好相反。

5 结 论

本文基于三维NACA 0012 系列展向可变形翼，开展了定常状态下给定静态展向变形和非定常状态下给定动态展向变形气动力特性及变化规律的数值模拟研究。主要结论如下：

（1）静态展向变形对三维机翼气动力性能略有影响，影响程度与攻角有关。正攻角下升力系数略有减小，阻力系数略有增大，负攻角则正好相反。不同静态展向变形幅度下对应的失速角完全一致。

（2）展向动态变形的存在对三维机翼阻力系数和升力系数平均值略有影响。阻力系数平均值均随变形频率的增大而微幅下降，升力系数变化规律与展向变形幅度有关。从兼顾阻力能耗和提高升力角度来考虑，以小幅度高频做动态展向变形时气动力性能会略微得到改善。

（3）不同变形频率下的升力系数和阻力系数时历变化规律几乎一致，呈类余弦曲线形式。阻力系数谷（峰）值点对应升力系数峰（谷）值点，出现在变形翼经平衡位置向下（上）做展向变形运动的瞬间；升力系数和阻力系数变化幅度均随变形频率的增大而增大，且变形幅度越大，变化幅度越明显。在变形翼设计阶段，可考虑以不同变形频率和不同变形幅度来规划一个展向变形周期内的运动，以达到气动力性能最优化的结果。

上述结论可为展向变形翼的气动力性能预报以及优化设计提供技术支撑，分析思路适用于变形翼其它变形形式下的展向或弦向变形研究。