6.二维高升力翼型NLR7301高精度气动数值模拟

1、中国航空学会2007年学术年会气动专题06二维高升力翼型NLR7301高精度气动数值模拟白文颜洪成水燕梁益华中国航空第一集团公司第六三一研究所，气动数值模拟航空科技重点实验室，陕西西安，710068摘要本文工作是建立大型飞机高升力设计高精度气动数值模拟能力系列工作的一部分。采用多个求解器，高精度 数值模拟了二维高升力两段翼型 NLR7301气动特性。系统研究了不同网格类型、网格密度、不同求解器、不同湍流 模型和转捩模型的计算精准度。采用的基本算法为有限体积法，流动模型为雷诺平均N-S（ RANS ）方程。与经典的试验数据确认研究表明，RANS技术可以比较精确地预测流动细节，以及阻力系数、最大升

2、力系数、失速迎角等关键气动特性。关键词高升力系统RANS方程湍流模型转捩NLR7301飞机的多段翼高升力增升系统设计及其空气动力特性研究一直为民机设计研究的前沿课题。近 十年来，雷诺平均N-S（ RANS ）方程方法已成为分析高升力系统的主要方法。针对高升力系统的CFD模拟，国际间广泛开展了高升力系统流场的计算、风洞实验以及它们之间的对比和评估工作。朱自强等1综述了高升力系统空气动力数值模拟计算已取得的成果、当前的水平、以及存在的 问题。由于高升力系统的几何复杂性及其大量的计算时耗，RANS方法目前在三维高升力系统模拟中应用较少，大多应用在二维简化模拟。在二维多段翼型绕流中也存在着各种令人困扰

3、并难以解决的 复杂流动物理现象，如分离、转捩、尾迹/边界层掺混、激波/边界层干扰等。目前二维模拟仍存在着大量问题：RANS方法的湍流模型应用仍具有最大不确定性；转捩模型对多段翼型的计算结果具有相对较大的影响，继续对转捩模型的研究是当前的一个研究热点；仍难以准确计算出CLmax值和失速迎角。Rumsey等对当前CFD计算高升力流场能力进行了评估2。2003年国内开展的“ CFD统一算例研讨”中选用了 NLR7301高升力翼型作为其研究算例之一7，当时整体计算收敛困难，计算的阻力系数相差较大。为了建立CFD在高升力翼型的工业应用能力及可信度，本文基于RANS方法，针对二维高升力两段翼型NLR730

4、1进行了 CFD高精度气动数值模拟。计算结果较为精确地预测了最大升力系数、失速 迎角，分析了计算流场的流动特性及流动现象，对比了CFD软件及其湍流模型在高升力翼型的应用情况，为CFD在民机高升力系统的应用打下了良好的基础。1算例描述本文选取典型的高升力两段翼型NLR7301作为研究内容。该模型先后在NLR Amsterdam 3m x 2m低速风洞和NLR North-East Polder 3m x 2.5m低速风洞进行了大量试验测量。为防止在襟翼上发生流 动分离，模型设计时襟翼偏转角设计为20°，是一个典型的起飞襟翼装置。风洞试验选取了两个不同缝间距的构型，两个构型的主翼和襟翼结

5、构相同，主翼襟翼搭接量为5.3%c，缝间距一个为1.3%c，另一个为2.6%c。风洞试验马赫数为0.185，雷诺数为2.51x106 （基于平均气动弦长 c=0.57m），试验测量出前缘有较小的层流分离泡，但对全局流场影响较小。压力分布的测量通过模型表面的大量测压 点获得。基于假定固壁附近的对数法则，利用“razor-blade ”技术获得了少量的表面摩擦力分布。升力和力矩系数则通过测量表面压力的积分得到，而阻力系数通过尾迹面法获得，尾迹面的位置选取 为离后缘相对一倍弦长距离。边界层剖面特性则通过边界层探针获得。模型表面的转捩位置的获得 则通过“ sublimation ”技术。所有测量数据已

6、出版34并可用于数值计算验证对比。该模型也经过大量的验证研究计算56，本文计算选择了 1.3%c构型。2数值方法2.1求解器介绍本文用于计算模拟的 CFD软件共三个： WiseManPlus、WiseMan和WoF90。这三个软件都是在有 限体积法的框架下求解可压缩RANS方程，其中前两个程序是基于格心格式求解N-S方程的多块结构网格求解器，而WoF90是基于格点格式求解N-S方程的混合网格求解器。三个求解器都包含动网格技 术，多种湍流模型，多重网格方法等，适用于飞行器低、亚、跨、超和高超声速大攻角定常、非定 常内外流动。三个求解器采用了不同的计算方法：在空间离散上，WiseManPlus和W

7、iseMan采用三阶迎风Roe格式，WoF90采用标准中心格式； 在时间推进上，WiseManPlus采用隐式近似因子分解方法，WiseMan采用隐式LU-SGS分解方法，WoF90采用三阶Runge-Kutta显式方法。同时，三个求解器都采用了多重 网格措施加速收敛。对于低速问题，可采用低马赫数预处理方法。鉴于已有的湍流模型应用经验 1，本文主要选取了三种涡粘性湍流模型和一种非线性湍流模型： Spalart-Allmaras模型，求解简单，适合于附着流动和小分离流动，为广泛应用的湍流模型之一。Wilcoxk- 3模型适合于模拟壁面剪切流，但结果过分依赖于 3的自由来流值。Men ter SS

8、T k -模型发挥了k- 3模型在近壁附近的稳定性又利用了k- &模型在边界层外部的独立性，能更精确地模拟近物面的湍流附面层，是比较出色的一个两方程模型。EASM k - 3模型是将EASM模型和k- 3两方程模型耦合到一起的非线性湍流模型。它给出了更高级的湍流剪应力表达式，包含了平均速度梯度，平均应变率和旋转应变率张量，克服了线性涡粘模型固有的缺陷。2.2边界条件本文在计算N-S方程时，三个计算软件壁面边界条件采用通常的无滑移、绝热壁条件，压力满足法向零梯度。同时WiseManPlus和WoF90在远场计算时采用了可压缩远场涡修正方法修正外部流动。 湍流模型中的参数在流入边界采用自由

9、来流值，流出边界可通过壁面上的值外推得到。3计算网格为了减少网格对计算精度的影响，参照ecarp5和文献 的网格研究，生成了本文的计算网格。计算网格按原始弦长无量纲化 （c），考虑远场的影响，网格的远场约100c左右。生成了三套网格用于网格收敛性分析，对于粗网格和细网格，在中等网格上保持网格拓扑关系和远场不变，按照每方向 约1.5倍关系修改网格点数得到。保持三套网格第一层网格点到固壁距离约为10-6数量级，满足y+ 1，为了精确模拟边界层效应，在高度为0.005c内的边界层分布网格点约 5080个。对于主翼和3中国航空学会2007年学术年会气动专题06副翼尾流区域，利用单独的网格块进行加密，保

10、证对尾流的准确捕捉。网格示意图如图1所示。图1.计算网格4中国航空学会2007年学术年会气动专题06#中国航空学会2007年学术年会气动专题064计算研究分析4.1网格收敛性分析为了排除网格对计算结果的影响，对网格进行收敛性分析。选取WiseManPlus求解器，在攻角为Alpha= 6.0 °的流动状态下对粗、中、细三套网格进行网格分析。表1网格收敛性分析计算结果网格点数(2D)Cl Cd粗网格112106 2.423060.035896中等网格231792 2.430430.025124细网格461046 2.463130.023395Richardson 夕卜推0.023015

11、试验值2.3660.0225上表给出了网格收敛性分析的计算各网格点数和计算的升阻力系数。从上表计算结果可以看出:三套网格的升阻力系数呈现出单调性，升力系数随着网格的加密单调递增，但偏离试验值；阻力系 数随着网格加密单调递减，靠近试验值。利用Richardson外推方法对三套网格计算的阻力系数进行外推，阻力系数收敛且外推得到的值与试验值相差约 力系数计算研究更加精确。+0.0005(5count)。可看出网格的收敛性分析对阻图2. Cp 曲线(Alpha=6.0 ° )图3. Cf 曲线(Alpha=6.0 ° )图2给出了三套网格计算的主/襟翼表面的Cp分布。三套网格计算的

12、Cp与试验数据吻合较好。从图可看出计算的主翼前缘的峰值都较试验值偏高，峰值随着网格点数的增加越来越大；由于不同密 度的网格对主翼的尾流捕捉情况不同，三套网格计算的Cp在襟翼表面差异较为明显，特别在襟翼的下表面表现最为明显。图3给出了三套网格计算的主/襟翼表面的Cf分布。此研究中由于采用全湍流计算，未能有效地模拟转捩现象，因此在层流区域和转捩区域计算的Cf与试验值有所差异。三套网格整体上较为准确预测了Cf的分布，细网格在附面层区域和尾流区域的加密使得在襟翼计算出现了转捩特征现象，但与试验值仍存在差异。分析看出网格加密对Cp曲线影响较小，少量的网格点即可准确地预测表面压力系数，但表面摩擦力系数Cf

13、的精确预测则需要高质量和足够密的网格。4.2求解器计算对比分析为了更加准确地评估CFD求解器在高升力翼型的应用能力，本文对三个求解器进行了对比分析 研究。三个求解器计算采用相同的中等网格，都采用SA湍流模型，对全部试验状态进行了模拟，分析了各求解器对预测G_max和失速迎角的能力。图4给出了三个求解器计算的升、阻力系数与试验数据的比较。结果表明，WiseMan计算的升力系数较试验值偏高，阻力系数较试验值偏差较大；WoF9 0计算的升力系数在小攻角情况下与试验值基本一致，阻力系数较试验值偏差较大；WiseMa nPlus计算的升、阻力系数较试验值偏大，但整体计算情况较好，计算的全部状态的升力系数

14、相差在4%以内，阻力系数相差约1020% o WiseMan和WoF90计算的失速迎角都约为 13.1 °，较试验值(约为14.1°)小1°左右。WiseMan计算的最大升力系数(CLmax= 3.045 )与试验值(试验值为 3.041 )相当，WoF90的最大升力系数(CLmax=2.95 )相对试验 值偏小；WiseMa nPlus准确地预测出失速迎角，计算的最大升力系数(CLmax= 3.153 )比试验值偏大。从曲线的整体趋势可以看出，三个求解器都较为精确的预测了升、阻力系数；在失速迎角前，三个 求解器计算的结果都与试验值吻合较好，预测的CLmax和失速

15、迎角都较为准确，失速后计算结果相差较大。综上表明结构和非结构求解器模拟高升力翼型的水平相当。图4.升阻力曲线5中国航空学会2007年学术年会气动专题06图5图8分别给出了三个不同攻角下三种求解器所得到的表面压力系数或摩擦系数与试验值的 对比曲线。从压力系数曲线图可看出，三个求解器的表面压力系数分布与试验值吻合较好，在主翼 前缘峰值都较试验值偏高，存在细微差别。WiseMa nPlus求解器计算的主翼峰值高于WoF90和WiseMan计算峰值。三个求解器的计算压力系数分布随攻角增大差异有所变大。图8表明,WiseManPlus与WoF90计算的摩擦系数结果比较接近，而WiseMan的计算结果整体

16、上与其它两个求解器偏差比较大。三个求解器都采用全湍流计算，摩擦系数模拟都未能捕捉转捩特征。综上三个求解 器都较准确地预测了翼型表面的压力系数分布，对于摩擦系数的模拟还存在细微差别。图 5. Cp 曲线(Alpha=6.0 ° )#中国航空学会2007年学术年会气动专题066中国航空学会2007年学术年会气动专题06图 6. Cp 曲线(Alpha=10.1 ° )0 04 0 03 WuMinFibiI口 JSPXJ1rJCf*二0.01.1 1 1 . .02414CISMI001-WMtan側卄1-0006:1 /广、-汀丿CM?:DQg.1 . 1 a . 1 I .

17、 Ia 1 . n 1 i I , a IU-fl 0第卿郵 11® 11 1.15 12图8. Cf 曲线(Alpha=6.0 ° )#中国航空学会2007年学术年会气动专题06#中国航空学会2007年学术年会气动专题06图 7. Cp 曲线(Alpha=13.1 ° )4.3湍流模型计算对比研究高升力系统采用RANS方法的最大不确定性之一来源于湍流模型，因此研究高升力系统适用的湍流模型成为当前重要的问题。借鉴文献1的总结，采用中等网格，针对攻角为Alpha= 6.0。的流动状态，选取WiseManPlus中含有的在高升力翼型模拟中较为广泛应用的几种湍流模型进行

18、了计算对比 研究。Sta.9Sta.12Sta.13Sta.14(x/c= 0.60Wi ng Lower) ( xf/c= 0.119Flap upper) (xf/c= 0.219Flap upper) (xf/c= 0.317 Flap upper)图11.边界层速度型分布表2湍流模型对比分析结果湍流模型Cl Cdexp 2.3660.0225SA 2.43040.0251Wilcox k-w2.43500.02651Menter SST2.407250.02542EASM k-w2.404750.02648从表2的结果比较可以看出：WiseMa nPlus求解器的各湍流模型计算的升力系

19、数相差不大，都在试验值3%以内，但阻力系数相差比较大，SA和SST模型的阻力系数更为接近试验值。图9给出了各湍流模型计算的Cp分布，从图可知各湍流模型计算的 Cp分布与试验数据吻合较好，各模型之间相差较小。图10给出了各湍流模型计算的 Cf分布，各湍流模型的计算的Cf分布存在差异。在主翼上表面各模型计算的Cf相差较小，襟翼则相差较大，表明各模型对尾流在襟翼上表面的影响不尽相同。图 11给出了不同湍流模型计算的不同站位的速度型，计算结果与试验数据吻合较好，正确地模拟出了 附面层特性。在主翼下表面都为全湍流计算，计算结果基本一致，但与试验值相差较大。各模型在 襟翼上表面的速度型较为准确地预测了尾流

20、作用，预测出了襟翼上表面在尾流影响下吸力峰值的减 少。但各湍流模型对尾流模拟能力的不同，导致相互之间计算的襟翼上表面速度型存在少量差异， SA和SST模型计算速度型更接近试验值。综合分析可知，SA、SST湍流模型计算的结果更加精确，能较好的模拟高升力翼型气动特性。4.4转捩研究表3翼型表面转捩位置主翼缝翼上表面下表面上表面下表面Alpha=6.0 °experiment0.0270.0400.6350.6751.0631.081laminarCFDx 旳.035x 旳.652x M.060laminarAlpha=10.1 °experiment0.0250.0370.68

21、00.7101.0541.072laminarCFDx «0.032x 旳.692x M.054laminarAlpha=13.1 °experiment0.0240.0350.6850.7151.0471.064laminarCFDx «0.029x 旳.705x M.048laminar转捩模型对多段翼型的计算具有相对较大的影响，使计算和实验的转捩点位置接近是计算获得成功的关键要素，对转捩模型的研究仍是当前的一个研究热点。为了验证转捩模型在高升力翼型的应用，本文采取固定转捩方法，对WiseManPlus求解器的转捩模型进行分析，湍流模型选用SA模型。三个流动状

22、态试验测得转捩位置如表3所示。图 12. Cp 曲线(Alpha=6.0 ° )图 13. Cp 曲线(Alpha=10.1 °)图 14. Cp 曲线(Alpha=13.1 ° )图 15. Cf 曲线(Alpha=6 ° )Sta.9Sta.10Sta.12Sta.13Sta.14(x/c=0.60Wing Lower) (x/c=0.84 Wing Lower) (xf/c= 0.119Flap upper) (xf/c=图5.速度型分布0.219Flap upper) (xf/c= 0.317 Flap upper)表4转捩分析计算结果ALPHA

23、=°ClCdexp2.3660.0225transition2.48440.02153no transition2.43040.0251表4给出了有无转捩计算的升阻力系数，结果表明：WiseMa nPlus求解器的转捩模型计算对结果有着显著影响。由于层流区的存在，转捩计算的升力系数都较无转捩计算偏大，而阻力系数减小更 接近试验值。转捩计算预测阻力更为精确。图12图14给出了 WiseManPlus求解器在不同攻角下有、无转捩计算的主翼/襟翼表面的Cp分布比较，计算的Cp分布与试验数据吻合较好。在主翼和襟翼上表面都准确的给出了转捩区域，转捩区 Cp变化较大。图15给出了 WiseMa

24、nPlus软件有、无转捩计算的摩擦系数分布，Cf分布与试验值吻合较好，转捩计算准确地预测出Cf在转捩区域急剧变化，模拟了真实流动的转捩特征现象。襟翼在有、无转捩计算时Cf相差较大，主要是因为层流区域对摩擦力影响较大。图5给出了有、无转捩计算的不同站位的速度型，转捩计算的速度型与试验值吻合较好。主翼的 Sta.9位置处于层流区域，Sta.10位置为转捩影响的尾流区域，转捩计算的速度型精确地模拟了转捩区 的边界层速度分布，与试验值非常接近。无转捩计算未能捕捉转捩特征，计算结果存在较大差别。 在襟翼上表面，有无转捩计算结果相差较小，但转捩计算对尾流的模拟能力更准确，计算的速度分 布更接近试验值。从整

25、体来看，转捩计算结果更能准确地预测翼型的流动特性，对边界层的模拟更 加接近试验值。通过二维高升力两段翼型 NLR7301算例的高精度气动数值模拟，得到以下几点结论：（1）对于高升力系统的CFD气动数值模拟，网格生成至关重要，网格质量对计算结果有着直接 影响。（2）从计算结果可以看出，转捩计算对结果有着显著影响，更能反映出真实流动。（3）从不同网格、不同方法和不同湍流模型的计算情况来看，计算结果与试验数据吻和较好， 验证了 CFD求解技术的准确性及可信度，为计算软件在大型飞机高升力系统设计的应用验证等后续 工作打下了良好的基础。（4）下一步工作将对网格的收敛性、三个求解器的多种典型湍流模型、多状

26、态的转捩计算及自动转捩模型的实现作进一步研究。感谢日本JAXA的雷忠博士提供相关论文和网格数据。参考文献1 朱自强，陈迎春，吴宗成，陈泽民，高升力系统外形的数值模拟计算”航空学报，Vol.26, No.3, 2005.2 Rumsey C L, Ying S X. Prediction of high lift: review of present CFD capability J . Progress in Aerospace Sci2ence,2002,38 :145 - 180. 23 Van den Berg, B.: Boundary Layer Measurements on a Two-Dimensional Wing with Flap, NLR TR 79009 U, 1979.4 B. van den Berg , J.H.M. Goode n: LOW-SPEED SURFACE PRESSURE AND BOUNDARY LAYERMEASUREMENT DATA FOR THE NLR 7301 AIRFOIL SECTION WITH TRAILING EDGE FLAP, AGARD