2021, Vol. 39
2. 中国飞行试验研究院, 西安 710089
2. Chinese Flight Test Establishment, Xi'an 710089, China
当飞行环境的温度低于零度时，云层中的过冷液滴与机体撞击会立即冻结，从而改变飞机气动外形, 对飞机的操纵性和稳定性有直接影响，甚至导致重大飞行事故[1]。根据美国联邦航空管理局(FAA)飞行安全部统计，从1990年到2000年，世界上由于结冰导致的飞行事故多达388起，约占气象飞行事故总数的12%[2]。
飞机结冰的相关研究始于20世纪30年代，考虑到飞行试验的安全性和经济性，早期的飞机结冰研究主要采用风洞实验。20世纪30年代末，随着CFD技术的迅速发展，关于积冰的数值模拟计算也开始兴起。1994年前的大部分结冰研究集中在结有较大明冰的对称翼型，而同年的美鹰航空4184号航班事故将人们的关注点转移到了不同翼型对不同冰形的敏感性分析上[3]。1999年Lee和Bragg[4]进行了不同展向冰脊模型对NACA23012和NLF0414翼型气动性能影响的研究，发现相比结冰形状，结冰位置对气动性能的影响要大得多，且两种翼型由于干净构型下不同的压力分布，在气动损失方面存在巨大差异。同年Dunn等[5]对具有简单襟翼的NACA23012翼型进行了SLD结冰条件下展向冰脊效应的数值模拟，结果表明相对较小的突起(高度为0.83%~1.39%弦长)会对升力产生显著影响。之后Chung等[6]将研究扩展到三维机翼上，发现二维和三维模拟结果趋势相似，后者比前者提前4°发生失速。2004年Pan和Loth[7]扩展了以前的工作，进一步考虑了雷诺数和马赫数的影响，发现雷诺数和马赫数对结冰翼型的影响远远小于对干净翼型的影响。
在飞行过程中，由于大气含水量、过冷液滴温度等环境因素以及不同翼型的外形影响，冰形在积累冻结的过程中除外形、结构外，在表面粗糙度上也会有很大差异。虽然目前关于冰形几何外形(如冰覆盖区域、冰最大厚度及其位置)对翼型气动力的影响已开展了大量研究[8-9]，但是，关于冰形表面粗糙度对翼型气动特性影响的研究却比较缺乏。目前对冰形表面粗糙度的研究主要包含两个方面的内容：一是对冰形粗糙表面的测量方法和表示形式进行说明，进而分析结冰条件对冰形表面粗糙度的影响[10-13]；二是分析冰形表面粗糙度对空气流传热以及表面水膜发展的影响，从而改进积冰预测模型[14-17]。
考虑到冰形表面粗糙度同样可能对翼型失速特性造成一定影响，本文参考文献[10]中对三维冰形表面粗糙度的测量及归类方法，总结了对二维冰形粗糙度影响相对较大的两个参数：冰形表面粗糙高度Ra和冰形表面粗糙分散度ρ，并在此基础上自主提出了一种粗糙度量化表示方法；之后采用该方法对光滑的冰形曲线进行了粗糙度叠加，并基于RANS方法模拟计算了不同结冰构型下结冰翼型的气动特性，最终分析总结出冰形粗糙度参数对干净翼型的失速特性影响规律，为进一步研究飞机结冰安全问题提供了更多参考依据。
1 冰形表面粗糙度 1.1 粗糙度参数冰形表面粗糙程度的增加会对气动力产生如下影响[3]：1) 增加机翼表面的摩擦力；2) 导致边界层转捩提前；3) 增加边界层厚度，提前后缘流动分离。本文参考文献[10]中对冰形表面粗糙度的测量和表示方法，结合二维冰形表面粗糙程度的变化规律，提取了如表 1所示的两个粗糙度参数。
| 表 1 冰形粗糙度参数 Table 1 Ice roughness parameters |
|
其具体物理意义如图 1所示，Ra控制冰形表面的轮廓偏距，其值为整块冰形轮廓偏距的方差，Ra越大则冰形表面越粗糙；ρ控制冰形表面轮廓单元的不平度间距，是一无量纲非负整数，其值为轮廓峰顶和谷底间隔坐标点数，ρ越小则冰形表面越粗糙。
|
图 1 冰形粗糙度参数物理意义示意图 Fig.1 Physical meaning of roughness parameters |
为真实地表示结冰状态，同时保证冰形的几何外形一致，在同一光滑的冰形表面添加粗糙度处理，通过输入冰形表面粗糙度高度Ra和冰形表面粗糙分散度ρ来控制冰形表面的粗糙程度。
通过对冰形厚度曲线上各点叠加一组随机数组B=(a1, a11, a12, …, a1ρ, a2, a21, a22, …, a2ρ, a3, …, an)(n >冰形厚度曲线数据点数)进行粗糙度修正：
| $ \bar{f}({{s}_{i}})=f({{s}_{i}})+{{b}_{i}} $ |
其中f(si)为修正后的冰形厚度，f(si)为光滑的初始冰形厚度，si对应翼型表面各点沿翼面到前缘点处的曲线长度，bi对应数组B中各元素。
当ρ=0时，B=(a1, a2, …, an)是均值为0、方差为Ra的正态分布随机数；当ρ≠0时，(aj1, aj2, …, ajρ)为(a1, a2, …, an)的三次样条插值。
输入不同Ra和ρ得到的冰形变化如图 2所示。
|
图 2 Ra和ρ变化对冰形表面粗糙程度的影响 Fig.2 Effects of Ra and ρ on the ice surface roughness |
本文的数值模拟采用ANSYS Fluent 18.0进行计算，控制方程为雷诺平均的N-S方程，空间离散格式为二阶迎风格式。根据Chi[18]和Chung[19]等对结冰翼型数值模拟时使用的各种湍流模型的研究结果，本文选用Spalart-Allmaras(S-A)湍流模型。
2.2 算例验证算例采用NACA0012结冰翼型(图 3)作为测试翼型，文中所用的结冰外形和结冰后的气动力实验数据均来自美国伊利诺伊大学厄巴纳-香槟分校的风洞实验结果[20]。
|
图 3 NACA0012结冰翼型 Fig.3 NACA0012 iced airfoil |
网格采用C型拓扑结构，网格量约85 000，壁面处y+=1，远场距壁面20个弦长。图 4给出了翼型附近的网格总览图以及前缘结冰部分的网格细节图。计算状态：Ma=0.18，Re=1.8×106。
|
图 4 NACA0012结冰翼型网格 Fig.4 NACA0012 iced airfoil mesh |
为检验网格无关性，依次将网格数沿x和y方向同时乘以1.4(2倍网格量)和0.7(0.5倍网格量)，分别计算这三种网格在4°~ 6°迎角下的气动力系数，得到结果如图 5所示。发现当网格量翻倍，在各迎角下计算得到的气动系数值变化很小，说明在现有网格基础上增加网格数量对计算结果的影响很小，可认为所采用网格已达到网格无关，因此采用该网格的参数设置进行后面的计算。
|
图 5 网格无关性检验结果 Fig.5 Results of grid independence test |
对比其在-8°~ 8°迎角范围内的升力系数和力矩系数实验数据(图 6)，发现气动力系数计算结果与实验结果基本吻合。使用基于S-A湍流模型的RANS数值模拟方法在研究结冰翼型的气动力中获得了较好的预测结果，后续将采用该方法对不同结冰情况下的翼型气动力进行计算评估。
|
图 6 NACA0012结冰翼型气动力系数比较 Fig.6 The (a) lift and (b) moment coefficients of NACA0012 iced airfoil |
采用控制变量的方法，保持其他特征参数不变，分别改变某一特定特征量，采用上一节提到的方法对不同结冰翼型进行气动影响评估。计算模型采用NACA23012结冰翼型，计算状态：Re=1.8×106，Ma=0.18。
3.1 流线冰型 3.1.1 参数设置各Case下的参数取值如表 2所示，得到的冰形如图 7所示，Ra=0与ρ=∞均为表面光滑的结冰翼型。
| 表 2 流线冰型粗糙度参数取值 Table 2 Parameters of rime ice roughness |
|
|
|
图 7 不同特征量下的流线型冰形曲线 Fig.7 Rime ice shapes with different parameters |
注意，在冰形粗糙度叠加程序中，冰形表面的粗糙度幅度变化是一组均值为Ra的随机数，为排除随机数造成的不确定性，此次分析时取粗糙度幅度变化均为一固定值Ra，即冰形曲线各点等幅振荡进而产生粗糙度变化。
3.1.2 失速特性影响分析计算不同特征参数结冰翼型在0°~ 14°迎角范围内的气动力系数，并分析其在失速迎角附近的气动特性。
比较冰形表面粗糙高度Ra的计算结果(图 8)发现，随着Ra的增大，由于冰形表面扰动幅度的增加，流动难以保持稳定状态，结冰翼型失速明显提前，在一定范围内，Ra每增加0.8，失速迎角提前3°左右；相应的最大升力系数明显减小，相比于光滑结冰翼型，Ra=0.8和Ra=1.5时的结冰翼型最大升力系数分别减小了近25%和45%。以10°迎角下的飞行状态为例(图 9)，当冰形表面光滑(Ra=0)时，流动经过前缘结冰区域仍可保持稳定流动，翼型上翼面无分离区域，前缘吸力峰明显；而当Ra增加到1.5时，由于粗糙度造成的扰动上翼面的流动近乎完全分离，吸力损失明显。
|
图 8 Case 1流线结冰翼型气动力系数曲线 Fig.8 The (a) lift and (b) moment coefficients of a rime iced airfoil in Case 1 |
|
图 9 流线型结冰翼型在10°迎角下的压力及流线图 Fig.9 Pressure contours and streamlines around a rime iced airfoil with α=10° |
而对于冰形表面粗糙分散度ρ(图 10)，发现随着ρ的减小，由于冰形表面扰动频率增加，失速性能有明显下降，翼型气动特性变化规律与Ra增大情况类似，但影响较Ra要小，ρ=1时的结冰翼型失速迎角相比ρ=2时提前了2°，最大升力系数减小约15%。
|
图 10 Case 2流线结冰翼型气动力系数曲线 Fig.10 The (a) lift and (b) moment coefficients of a rime iced airfoil in Case 2 |
采用同样的方法对角状冰形的表面粗糙度参数进行分析，各Case下的取值如表 3所示，得到的冰形如图 11所示。
| 表 3 角状冰型粗糙度参数取值 Table 3 Parameters of glaze ice roughness parameters |
|
|
|
图 11 不同特征量下的角状冰形曲线 Fig.11 Shapes of glaze ice with different parameters |
计算不同特征参数结冰翼型在0°~ 6°迎角范围内的气动力系数，并分析其在失速迎角附近的气动特性。
对于冰形表面粗糙高度Ra，通过比较角状结冰翼型在4°迎角下的压力及流线图(图 12)发现，随着Ra的增大，由于冰形表面扰动幅度的增加，流动难以保持稳定状态，导致结冰翼型上表面分离区域变大，造成失速迎角与最大升力系数均变小，气动性能损失增加。但由于角状冰型自身几何外形特征会产生较大的压差阻力，导致冰角下游流场分离，造成较大的气动力损失，因此Ra增大产生的摩擦阻力所造成的气动损失相比之下要小得多。
|
图 12 Case 1中时的角状结冰翼型在4°迎角下的压力及流线图 Fig.12 Pressure contours and streamlines around a glaze iced airfoil in Case 1 at AoA=4° |
由不同Ra在不同迎角下的气动力变化曲线(图 13)可知，相比于表面光滑(Ra=0)的结冰翼型，Ra=0.8时的结冰翼型失速迎角基本不变，最大升力系数只减小了不到7%，Ra=1.5时的结冰翼型失速迎角减小1°左右，最大升力系数减小约19%。
|
图 13 Case 1角状结冰翼型气动力系数曲线 Fig.13 The (a) lift and (b) moment coefficients of a glaze iced airfoil in Case 1 |
同样比较不同冰形表面粗糙分散度ρ在不同迎角下的气动力变化曲线(图 14)发现，ρ对气动力造成的影响要略小于冰形表面粗糙高度Ra；同样由于角状冰型自身几何外形特征造成的气动力损失相比冰形表面粗糙分散度ρ减小造成的气动损失要大得多，因此结冰翼型的失速迎角与最大升力系数呈不规则变化，ρ=∞(即光滑冰形)时结冰翼型的气动性能损失依旧最小，ρ=2和ρ=3下得到的气动系数十分接近，但是气动特性均要略劣于ρ=1，这与之前流线结冰翼型总结的规律相反。
|
图 14 Case 2角状结冰翼型气动力系数曲线 Fig.14 The (a) lift and (b) moment coefficients of a glaze iced airfoil in Case 2 |
本文初步开展了冰形表面粗糙度叠加方法的研究，并采用CFD手段评估了翼型结冰后对气动特性的影响，重点分析了冰形表面粗糙度对翼型失速特性的影响规律，得到如下结论：
1) 对于本文提出的两类冰形粗糙度参数，冰形表面粗糙高度Ra对翼型气动特性的影响要略大于冰形表面粗糙分散度ρ；
2) 对于流线结冰翼型，由于其冰形本身的几何外形对翼型的气动特性影响较小，因此冰形粗糙度参数的改变会对翼型气动力造成较明显的影响，在研究飞机结冰安全问题时需将其考虑在内；
3) 对于角状结冰翼型，由于两个凸起的羊角状外形会产生较大压差阻力，导致明显的气动力损失，因此冰形粗糙度参数的改变对翼型气动性能的影响可基本忽略，只需考虑几何外形的影响。
本文给出的冰形表面粗糙度表示方法可在一定程度上反映冰形粗糙度对翼型气动特性的影响规律，为研究飞机结冰问题提供思路。但在实际情况下，气流沿展向也会有流动变化，因此可通过增加数组B维数的方式，进一步将该方法扩展到三维冰形上，对结冰机翼的气动力进行粗糙度量化分析，更加准确地反映真实情况下冰形粗糙度对飞行器气动特性的影响。
| [1] |
易贤. 飞机积冰的数值计算与积冰试验相似准则研究[D]. 绵阳: 中国空气动力研究与发展中心, 2007, 12. YI X. Numerical computation of aircraft icing and study on icing test scaling law[D]. Mianyang: China Aerodynamics Research and Development Center, 2007, 12. (in Chinese) |
| [2] |
CEBECI T, KAFYEKE F. Aircraft icing[J]. Annual Review of Fluid Mechanics, 2003, 35(35): 11-21. |
| [3] |
BRAGG M B, BROEREN A P, BLUMENTHAL L A. Iced-airfoil aerodynamics[J]. Progress in Aerospace Sciences, 2005, 41(5): 323-362. DOI:10.1016/j.paerosci.2005.07.001 |
| [4] |
LEE S, BRAGG M. Effects of simulated-spanwise-ice shapes on airfoils-experimental investigation[C]//Aerospace Sciences Meeting & Exhibit. 2013.
|
| [5] |
DUNN T A, LOTH E, BRAGG M B. Computational investigation of simulated large-droplet ice shapes on airfoil aerodynamics[J]. Journal of Aircraft, 1999, 36(5): 836-843. DOI:10.2514/2.2517 |
| [6] |
CHUNG J, REEHORST A, CHOO Y, et al. Navier-Stokes analysis of flowfield characteristics of an ice-contaminated aircraft wing[J]. Journal of Aircraft, 2000, 37(6): 947-959. DOI:10.2514/2.2716 |
| [7] |
PAN J, LOTH E. Reynolds-averaged Navier-Stokes simulations of airfoils and wings with ice shapes[J]. Journal of Aircraft, 2004, 41(4): 879-891. DOI:10.2514/1.587 |
| [8] |
LEE S, BRAGG M B. Experimental investigation of simulated large-droplet ice shapes on airfoil aerodynamics[J]. Journal of Aircraft, 1999, 36(5): 844-850. DOI:10.2514/2.2518 |
| [9] |
朱东宇, 张付坤, 裴如男, 等. 临界冰形确定方法及其对气动特性影响研究[J]. 空气动力学学报, 2016, 34(6): 714-720. ZHU D Y, ZHANG F K, ZHOU Z H, et al. Research on critical ice shape determination and its effects on aerodynamics[J]. Acta Aerodynamic Sinica, 2016, 34(6): 714-720. (in Chinese) |
| [10] |
ANDERSON D N, HENTSCHEL D B, RUFF G A. Measurement and correlation of ice accretion roughness[C]//36th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV: AIAA, 1998.
|
| [11] |
FORTIN G, ILINCA A, LAFORTE J L, et al. New roughness computation method and geometric accretion model for airfoil icing[J]. Journal of Aircraft, 2004, 41(1): 119-127. DOI:10.2514/1.173 |
| [12] |
常士楠, 李延, 霍西恒. 飞机结冰表面粗糙度模型建立与分析[C]//中国力学学会学术大会. 郑州: 中国力学学会, 2009. CHANG S N, LI Y, HUO X H. Analysis and modeling of aircraft icing surface roughness[C]//CCTAM. ZhengZhou: CSTAM, 2009. (in Chinese) |
| [13] |
SHIN J. Characteristics of surface roughness associated with leading-edge ice accretion[J]. Journal of Aircraft, 2015, 33(2): 316-321. |
| [14] |
BRAGG M, CUMMINGS M, HENZE C. Boundary-layer and heat-transfer measurements on an airfoil with simulated ice roughness[C]//Aerospace Sciences Meeting & Exhibit. 1996.
|
| [15] |
BRAGG M B, BROEREN A P, BLUMENTHAL L A. Iced-airfoil aerodynamics[J]. Progress in Aerospace Sciences, 2005, 41(5): 323-362. DOI:10.1016/j.paerosci.2005.07.001 |
| [16] |
MATHEIS. Numerical investigation of roughness effects in aircraft icing calculations[D]. Iowa State University, 2008.
|
| [17] |
李瀚涛, 舒立春, 胡琴, 等. 考虑覆冰粗糙度影响的风力发电机叶片气动性能数值仿真[J]. 电工技术学报, 2018, 33(10): 2253-2260. LI H T, SHU LC, HU Q, et al. Numerical simulation of wind turbine blades aerodynamic performance based on ice roughness effect[J]. Transactions of China Electrotechnical Society, 2018, 33(10): 2253-2260. (in Chinese) |
| [18] |
CHI X, LI Y, CHEN H, et al. A comparative study using CFD to predict iced airfoil aerodynamics[C]//AIAA Aerospace Sciences Meeting & Exhibit. Reno. NV: AIAA, 2005.
|
| [19] |
CHUNG J J, ADDY H E. A numerical evaluation of icing effects on a natural laminar flow airfoil[C]//38th Aerospace Sciences Meeting. Reno, NV: AIAA, 2000.
|
| [20] |
GURBACKI H M. Ice-induced unsteady flow field effects on airfoil performance[D]. Urbana, IL: Univ. of Illinois at Urbana-Champaign, 2003.
|