翼面结冰对飞机飞行安全会产生严重影响,较强结冰条件下形成的角状或楔状冰将破坏翼型的前缘形状,使其绕流特性显著改变,发生非定常分离和再附现象,形成时均流动分离泡,导致升力系数下降和失速迎角提前[1, 2]。对翼型结冰后流场和气动特性变化情况的精确预测依赖于准确把握与分离现象相关的流动特征。特别对于结冰翼型失速点附近的分离流动问题而言，其流场结构较为复杂，存在丰富的旋涡相互作用，非定常效应占主导地位[3, 4]，这对湍流模拟的精度和可靠性提出了较高的要求。但由于工程领域常用的雷诺平均(RANS)方法对湍流脉动信息采用时均化处理方式，难以对此类复杂分离流动的细节特征进行合理描述，需要应用更符合物理事实、精度更高的湍流数值模拟方法。
近年来发展的DES[5](Detached Eddy Simulation)类混合方法能够在一定程度上兼顾分离流动的计算精度和求解效率。其基本思想是利用RANS湍流模型中包含的长度尺度项构造某种混合长度尺度，以使湍流模型在壁面附近区域体现RANS方法的性质，在使用薄层网格的前提下避免近壁面雷诺应力损失，同时降低对计算资源的要求；在以大涡输运为主要特征的远离壁面流动分离区域则体现大涡模拟(LES)方法的特点，使当地雷诺应力快速降低，对大尺度湍流结构进行解析求解，保证空间湍流运动模拟的精度。
目前，利用DES类方法对翼型结冰后流动进行分析的相关研究工作已经取得了一定进展。Mogili[6]等运用DES方法研究了失速点附近结冰翼型的气动特性，证明对于此类问题DES方法较RANS方法能够取得更好的气动力预测结果；Thompson[7]等对不同迎角下的翼型结冰后流动问题做了较为系统的DES分析研究，阐述了网格尺度对数值模拟结果的影响，认为在较大迎角下DES方法对于流动分离特性的模拟精确程度尚有所欠缺；Lorenzo[8]等分别采取DES和DDES[9](Delayed Detached Eddy Simulation)方法对平尾翼型结冰失速前后的流场特征进行了分析，同样指出以上数值方法对于结冰后分离流动的精确预测尚存在一定不足。
针对DES类方法存在的相关问题，Travin[10]等将DDES方法与LES壁面模型(Wall-Modelling in LES，WMLES)相结合，提出了IDDES方法。该方法有利于分离区域内中小尺度湍流结构的充分解析，同时在附着和分离流动并存的过渡区域能够取得更为满意的结果。已经证明该方法对于干净翼型失速分离问题的分析是适用的[10]，但与结冰失速分离流动相关的应用研究还并不充分。
本文基于IDDES方法就翼型前缘角状冰导致的失速点附近流动分离问题进行数值模拟研究，以期检验数值方法对此类分离流动的模拟能力，探讨翼型结冰影响流场特征的相关物理机制，为DES类方法在结冰后翼型/机翼气动力分析中的实际应用提供参考。
计算模型及网格针对一种通用航空公务机使用的翼型GLC305开展计算分析研究，翼型前缘具有结冰时间为22.5min的典型双角状积冰[11]。积冰上下角状结构与翼型表面近似垂直，高度分别约为翼型弦长c的3%和2%。图 1描述了该角状积冰的几何形状。
Broeren[12]等在NASA Langley低湍流度风洞(LTPT)针对该翼型前缘结冰后的气动力变化问题进行了多组试验。本文选取的计算状态来流马赫数Ma=0.12，基于翼型弦长的雷诺数Re=1.8×106，来流迎角α=6°，对应风洞试验失速点附近状态。
|
| 图 1 GLC305翼型前缘双角状积冰 Fig. 1 Horn ice on the leading edge of GLC305 airfoil |
对该带冰翼型三维几何模型生成多块点对接结构化网格，采用图 2所示的C-H型网格拓扑结构，展向长度取0.5c。根据Spalart[13]的观点，应用DES类方法进行计算时，计算域可分为欧拉区域(ER)、大涡模拟区域(LR)及雷诺平均区域(RR)三类，需要根据当地流动特征设置不同的网格尺度，各区域的详细定义见文献[13]。图 2给出了计算域分区及相应的网格布置情况。
|
| 图 2 计算域分区及相应的网格布置 Fig. 2 Computational region and corresponding grid (a) 计算域分区情况 (b) 网格布置情况 |
为达到节省计算资源的目的，对ER区域和LR区域使用不同周向密度的网格，在区域交界面处采用面搭接技术进行流场信息传递，图 2中的红色实线表示网格搭接面位置。搭接面附近的网格布置情况如图 3所示。
|
| 图 3 搭接面附近网格空间截面 Fig. 3 Computational grid near the patch surface |
在翼型上表面LR区域内布置三向同性网格单元，单元尺度取0.5%c。在物面附近RR区域内布置薄层网格单元，首层网格到壁面的法向距离为1×10-5c，增长率取1.1，以保证壁面附近y+值不大于1。计算域内网格结点总数约为1.6×107。流动关注区域附近计算网格如图 4所示。
|
| 图 4 流动关注区域附近计算网格 Fig. 4 Computational grid of the focus region |
计算域远场给定无反射边界条件，物面采用绝热、无滑移和法向零压力梯度条件，展向设置周期性边界条件。
2 数值模拟方法 2.1 控制方程及其离散在有限体积法基础上，对三维可压缩非定常N-S方程进行求解：
(1) 其中Q为守恒变量，F、G、H为三向无粘通量，FV、GV、HV为三向粘性通量。对无粘通量项采用Roe-MUSCL三阶迎风通量差分分裂格式，粘性通量项采用二阶中心差分格式进行空间离散；时间推进采用二阶隐式近似因子分解方法。
2.2 湍流模拟方法根据文献[10]对 IDDES方法进行构造，实现湍流数值模拟。该方法基于k－ω SST两方程湍流模型[14]，相对DDES方法的主要改进内容包括通过引入壁面距离和网格间距参数对亚格子尺度进行重新定义；以及将DDES方法用作一种WMLES构造了RANS/LES混合长度尺度。
2.2.1 亚格子尺度重新定义LES方法中常将亚格子尺度取为网格单元体积的立方根，而在DES方法中将其取为网格单元三向长度的最大值。但在以上两种亚格子尺度定义方式中，对各向同性和自由剪切等不同性质的湍流流动，最优亚格子应力模型常数需要取不同的值，对一般情况下的湍流数值模拟缺乏普遍适用性。针对该问题，文献[10]中重新定义了一种亚格子尺度：
(2) 式(2)中hwn为沿壁面法向当地网格单元尺度，Cw为由LES解得到的经验常数，hmax为网格单元三向最大尺度，dw为网格单元与壁面距离。
2.2.2 RANS-LES 混合长度尺度构造基于上述重新定义的亚格子尺度，构造一种新的湍流混合长度尺度
对于本文使用的SST湍流模型，该混合长度尺度具备如下形式：
(3) 式中lRANS为RANS长度尺度，Δ为式(2)所定义的亚格子尺度，CDES为亚格子应力模型常数。函数frestore和fhyb的具体形式见文献[10]。
3 计算结果与分析非定常计算设置无量纲时间步长Δt*=0.0025，等价于物理时间步长Δt=0.067ms，每时间步内取15步牛顿子迭代。首先基于非定常RANS方法获得充分发展的初始流场，在初场基础上进行后续IDDES计算至非定常流场基本稳定，以有效滤除分离区域的RANS效应。
以下分别从时均角度和瞬态角度对IDDES计算结果分别进行分析，并与LTPT风洞试验结果[12]及相同网格布置下的非定常RANS方法计算结果进行比较。以对IDDES方法关于此类分离流动的模拟能力进行较为全面的评估和考核。
3.1 时均结果图 5给出了IDDES方法计算基本稳定阶段所得气动力系数的时间序列，图中实线表示各气动力系数瞬时值，虚线表示气动力系数时间平均值。序列的不规则锯齿状分布反映流场中存在强烈的非定常脉动现象,表明数值方法具备对流场内部中小尺度的湍流细节变化进行捕捉的能力。
|
| 图 5 气动力系数时间序列 Fig. 5 Time history of aerodynamic force coefficients |
表 1将计算所得的时均宏观气动力与试验值进行了对比。由表 1可知IDDES方法较RANS方法在升力量值预测上有一定提升,表明 IDDES方法能够更为准确地反映翼型表面的压力分布特征。阻力量值预测结果精度则与RANS方法相当。
| CL | Error | CD | Error | |
| EXP | 0.660 | / | 0.105 | / |
| RANS | 0.701 | 6.21% | 0.118 | 12.38% |
| IDDES | 0.688 | 4.24% | 0.093 | -11.43% |
图 6给出了计算所得时均流场压力分布与试验结果的对比情况。由图 6可知RANS方法无法反映因大尺度时均分离泡存在而形成的压力平台。IDDES方法计算所得压力平台长度能够与试验值良好吻合；但压力恢复过程相对试验结果则略为陡峭，该特点在Alam[16]等的DDES方法计算结果中也有所体现。在翼型下表面的局部分离区域内，两种方法同试验结果的吻合程度均较为良好，显示了IDDES 方法对近壁面弱分离流动计算的可靠性。不过IDDES方法在翼型后缘的压力分布计算结果较试验值有所偏移，Thompson[7]等利用DES方法也得到了相似的压力分布情况，这与翼型干净无冰情况下的试验结果比较接近。
|
| 图 6 时均流场压力分布对比 Fig. 6 Comparison of time-averaged pressure distribution |
图 7展示了时均流场速度分布对比情况，其中分离泡平均外廓和流动平均再附位置由0速度等值线描述。PIV试验测量得到的再附线位置约为53%c附近。由图可知，RANS计算所得分离泡形状和体积与试验结果存在较大差异，流动直到翼型后缘附近才发生再附。与之相较，IDDES计算所得流动再附位置约为50%c处，与试验结果大体相同；分离泡形状与试验结果较为相似，体积与之相当。不过在10%c附近，计算所得回流区域的强度与试验结果仍存在一定差距。
|
| 图 7 时均流场U向速度分布对比 Fig. 7 Comparison of time-averaged U velocity contours |
图 8以流线形式直观给出了不同数值方法计算所得的时均分离泡，并指示了分离泡的时均再附位置。可见IDDES方法有效描述了分离泡的几何形状和再附形态。
|
| 图 8 计算所得时均分离泡形状对比 Fig. 8 Comparison of time-averaged separation bubble |
图 9对比了翼型上表面附近不同位置的时均速度型计算结果。图 9(a)给出的速度型位于分离泡起始产生区域，描述了壁面附近的强回流趋势；IDDES方法略为保守地估计了流动速度沿法向的变化，这与Alam[16]等利用动态RANS/LES混合方法和DDES方法给出的计算结果较为相似。图 9(b)反映了翼型中部再附区域(0.4<x/c<0.6)的速度分布；IDDES计算结果与试验值吻合良好，显示了精确预测流动再附位置的能力。图 9(c)和图 9(d)体现了再附区域下游充分发展的附着流动速度分布，IDDES方法预测结果趋势与试验值相一致；但在壁面附近过于乐观地估计了流动速度的恢复情况。总体而言IDDES方法能够在壁面附近区域较RANS方法提供更为良好的速度分布预测结果，这与图 7提供的宏观速度分布结果相吻合。
|
| 图 9 时均流场近壁面不同位置速度型对比 Fig. 9 Comparison of time-averaged velocity profiles at different stations near the wall |
3.2 瞬态结果
图 10给出了计算所得瞬态流场展向涡量分布情况。RANS计算结果中角状冰后方的脱落剪切层一直延伸到翼型中部位置，并未显示析出旋涡的趋势，这与真实流动现象是不相符的，表明分离区域湍流涡粘系数过高。而IDDES计算结果中剪切层失稳破碎的起始位置靠近翼型前缘，能够体现外部流动与回流区域切向速度差导致的强烈剪切效应。表明在该算例当中IDDES方法能够有效和快速地降低壁面附近分离区域的涡粘系数水平，减少其对流场细节求解的影响。注意剪切层失稳大约起始于25%弦长左右，与试验结果中的压力恢复起始位置相当，表明压力恢复的直接原因是角状冰后方脱落剪切层的破碎和旋涡析出过程。
|
| 图 10 瞬态流场展向涡量分布对比 Fig. 10 Comparison of instantaneous spanwise vortices distribution |
此外，IDDES方法在分离区域能够获得比较丰富的旋涡结构。表明计算网格的单元尺度基本能够满足LR区域精细解析湍流结构的要求。由图 10可知，旋涡向壁面的输运促进了外部流动与回流区域流动的掺混，使得流动在翼型表面发生再附。由于旋涡的析出和输运过程具备相当程度的随机性，使得再附位置并不稳定于翼型表面，而是在一定范围内振荡，形成再附区域，令分离泡的体积和形状都随时间发生变化。这种再附位置的低频变化是导致文献[3]中提到的结冰翼型气动力特别是升力波动的主要原因。
图 11给出了不同时刻计算所得瞬态Q等值面分布，Q准则由Hunt[17]等提出，用于表示流场中旋涡旋转量和拉伸量间的关系。由图 11可知，IDDES方法在关注区域内不仅能够获得流场中主要的三维大尺度湍流结构，也能够捕捉到更加细微的中小尺度结构，精确描述了剪切层沿弦向逐步失稳破碎，内部旋涡脱落滚转进入下游的过程；并且处于网格加密区域内的流动尾迹仍大体能够被较为完整地加以解析。由图中同样能够观察到分离泡的形状及再附位置随时间的变化情况。
|
| 图 11 不同时刻IDDES计算所得瞬态Q等值面分布 (Q=0.1) Fig. 11 Instantaneous Q criterion of IDDES in different times (Q=0.1) |
4 结 论
在IDDES方法的基础上针对GLC305翼型失速点附近前缘角状冰诱导产生的分离流动进行了数值模拟研究，并与NASA LTPT风洞试验结果进行了对比。数值模拟结果表明：
1) IDDES方法适用于前缘角状冰导致的翼型失速分离问题的精细数值模拟，较之非定常RANS方法可更好地描述分离流动的物理特征。能够在壁面附近取得良好的速度预测结果，有效解析分离区域内的中小尺度旋涡结构，较为准确地反映大尺度时均分离泡的再附位置和形态特点。
2) 对于失速点附近翼型结角状冰引起的流动分离问题而言，由于脱落剪切层内部的不稳定旋涡析出和输运过程促进了外部流动与回流区域流动的掺混，将会导致流动发生非定常再附现象。
| [1] | Gurbacki H M, Bragg M B. Unsteady flowfield about an iced airfoil[R]. AIAA 2004-0562. |
| [2] | Broeren A P, Bragg M B, Addy H E, et al. Effect of high-fidelity ice accretion simulations on the performance of a full-scale airfoil model[R]. AIAA 2008-434. |
| [3] | Ansell P J, Bragg M B. Measurement of unsteady flow reattachment on an airfoil with a leading-edge horn-ice shape[R]. AIAA 2012-2797. |
| [4] | Ansell P J, Bragg M B. Characterization of ice-induced low-frequency flowfield oscillations and their effect on airfoil performance[R]. AIAA 2013-2673. |
| [5] | Spalart P R, Jou W H, Strelets M, et al. Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach[M]. Los Angles: Greyden Press, 1997. |
| [6] | Mogili P, Thompson D S, Choo Y, et al. RANS and DES computations for a wing with ice accretion[R]. AIAA 2005-1372. |
| [7] | Thompson D S, Mogili P. Detached-eddy simulations of separated flow around wings with ice accretions: year one report[R]. NASA CR-2004-213379. |
| [8] | Lorenzo A, Valero E, De-Pablo V. DES/DDES post-stall study with iced airfoil[R]. AIAA 2011-1103. |
| [9] | Spalart P R, Deck S, Shur M, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J].Theoretical and Computational Fluid Dynamics, 2006, 20(3):181–195. |
| [10] | Travin A K, Shur M L, Spalart P R, et al. Improvement of delayed detached-eddy simulation for LES with wall modeling[C]//European Conference on Computational Fluid Dynamics. The Netherlands: TU Delft, 2006. |
| [11] | Addy H E. Ice accretions and icing effects for modern airfoils[R]. NASA TP-2000-210031. |
| [12] | Broeren A P, Bragg M B, Addy H E. Flowfield measurements about an airfoil with leading-edge ice shapes[J].Journal of Aircraft, 2006, 43(4):1226–1234. |
| [13] | Spalart P R. Young-person's guide to detached-eddy simulation grids[R]. NASA CR-2001-211032. |
| [14] | Menter F R. Two-equation eddy viscosity turbulence models for engineering applications[J].AIAA Journal, 1994, 32(8):1598–1605. |
| [15] | Nikitin N V, Nicoud F, Wasistho B, et al. An approach to wall modeling in large-eddy simulations[J].Physic of Fluids, 2005, 12(7):1629–1632. |
| [16] | Alam M F, Thompson D S, Walters D K. Hybrid Reynolds-averaged Navier-Stokes/large-eddy simulation models for flow around an iced wing[J].Journal of Aircraft, 2015, 51(1):244–256. |
| [17] | Hunt J, Wray A, Moin P. Proceedings of the 1988 summer program[R]. Stanford: Center for Turbulence Research, 1988. |