带压缩性修正的离散涡方法在结冰数值模拟中的应用

2021-06-28周泽堃周峰王刚

航空工程进展 2021年3期

周泽堃，周峰，王刚

（1.西北工业大学航空学院，西安710072）

（2.中国商用飞机有限责任公司上海飞机设计研究院，上海201210）

0 引言

飞行器结冰是危害飞行安全的严重问题之一，轻则影响飞行器的飞行性能，重则造成机毁人亡的惨重后果［1-2］。而对飞机进行结冰数值模拟则成为一种研究飞行器结冰问题的有效途径。国际上商用的结冰数值模拟软件包括美国NASA的Lewice［3］和加拿大的FENSAP-ICE［4］等。

在结冰数值模拟研究中，国内起步较晚，目前也做了很多深入研究，易贤等［5］采用欧拉法计算了三维复杂外形表面的水滴收集率；周峰等［6］对部分典型二维翼型进行了结冰数值模拟；Xie L等［7］提出一种自适应插值方法将拉格朗日法计算得到的水滴收集率插值到三维复杂外形上；桑为民等［8］在进行大粒径水滴状态下的结冰模拟时考虑了水滴变形破碎的影响；杜雁霞等［9］对结冰过程涉及的热力学行为的研究现状进行了总结和展望；王海涛等［10］对机翼进行结冰数值模拟并设计了一套防除冰系统。

飞行器结冰数值模拟主要包括网格生成、计算流场、水滴撞击特性计算、结冰模型计算和结冰冰形生长计算五个模块。在流场计算时，通常采用经典的有限体积法求解N-S方程。相比于经典CFD方法，离散涡方法［11］可以不考虑网格，数值耗散低，适用于结冰这类易产生分离流动问题的分析。国内外研究者对离散涡方法进行了应用和发展。国外，K.Ramesh等［12］运用离散涡方法模拟翼型前缘产生的间歇涡运动；E.G.A.Antonini等［13］提出一种改进的离散涡方法来模拟二维翼型俯仰中的动态失速现象。国内，马明宪等［14］使用离散涡方法计算了湍流混合层的问题；董婧等［15］运用离散涡方法模拟了圆柱绕流问题；吴文权等［16］采用拉格朗日粒子法对胀量项计算，从而实现离散涡方法在可压缩流动下的计算；刘佳等［17］将离散涡方法应用于计算结冰数值模拟，与经典CFD下的结冰数值模拟过程相比，在保证精度的同时，计算效率得到提升。但是离散涡方法主要针对不可压流动求解，对于可压缩流动下的冰形模拟的应用仍有所欠缺。

本文在基于离散涡方法结冰模拟的基础上，引入普朗特—格劳尔特压缩性修正，将该方法的计算结果拓展到可压缩区间，结合改进的Messinger热力学模型［18］实现可压缩状态下翼型的结冰模拟。通过部分二维翼型的结冰数值模拟，将模拟结果与基于有限体积法的经典CFD方法下模拟结果及实验值进行对比，并分析引入压缩性修正前后对流场、结冰过程和结冰冰形的影响。

1 结冰数值模拟方法

1.1 离散涡方法基本原理及改进

二维不可压N-S方程的控制方程［19］可表示为

式中：ν为动力黏度；u为速度矢量；ω为涡量。

离散涡方法是将流场涡量离散化，在拉格朗日框架下模拟流场涡量。通过采用算子分裂法，式（2）可以分成两部分，即对流项式（4）与扩散项式（5）。

在二维不可压无黏条件下，得到：

流体微团在运动过程中涡量始终不发生改变，并且不随流体质点迁移运动。

式中：r为涡元距离原点的位移矢量；x s为涡元所处位置的坐标。

为了避免点涡在距离其无穷近处出现诱导速度无穷大的奇性问题，通过A.J.Chorin［20］提出的涡团法，涡量分布方程可变为

式中：N为涡团数；Γj为涡团的环量；δσ(x-x j)为涡量的分布函数。

在得到涡量分布后，依据Biot-Savart定律［21］得到速度场：

计算得到流场的速度分布后，就可以获得涡元的运动轨迹，再计算出下一时刻的涡量场。通过随机走步法［22］来模拟涡元的扩散运动，即：

式中：μ为满足高斯分布的随机数，其均值为零，方差为

在大多数工程问题中，必须满足特定的边界条件，流场才能求解。根据亥姆霍兹（Helmholtz）定理，得出在均匀流中不会产生涡量，只有在物体边界或者流场受非保守力作用时才会有涡量产生。本文通过壁面无穿透边界条件确定流场中生成的涡量。整个流场在满足Kelvin环量定理的前提下，通过速度场与涡量场之间交替迭代计算，就能模拟出Lagrangian框架下的绕流流动。

上述是离散涡方法的基本原理，在基本原理的基础上本文对离散涡方法进行一些改进，包括涡元合并、控制区高度修订及压力计算［23］等，具体参见文献［17］。

1.2 二维结冰热力学模型

基于改进的Messinger热力学模型［18］，将结冰表面分为若干控制体，每个控制体的质量流动情况如图1所示。在一些假设条件下，该模型通过求解质量和能量守恒方程［24］最终求解出结冰质量。

图1 控制体内质量流量示意图Fig.1 Sketch of mass variation in control volume

建立质量守恒方程（11）和能量守恒方程（12），定义每个控制体内的冻结系数f如式（13）所示。

式中：为控制体内的水滴撞击质量；为液态水溢流进入控制体的质量为液态水溢流出控制体的质量；为升华和蒸发质量；为结冰质量；为控制体内水滴撞击的动能；为控制体内水滴撞击的焓变；为液态水溢流进控制体的焓变为气动摩擦热；为结冰所释放的潜热；为液态水溢流出控制体的焓变；为液态水蒸发潜热；为对流换热热能。

联立式（11）、式（13），通过提取表面温度作为试探参数，对方程组求解，从而确定壁面各控制体内的壁温、冻结系数和各个质量项。在得到结冰表面各控制体内的每个质量项后，就可以依据结冰增长模型，计算出单次时间步长内冰层的厚度。

1.3 普朗特—格劳尔特压缩性修正

离散涡方法适用于不可压流动，在当下的研究进展中，已经实现可压缩流动下的离散涡方法。但是引入压缩性修正，将离散涡方法的结果扩展到可压缩区间仍是一种简易的方法。

对二维翼型壁面附近的绕流速度采用小扰动假设，如图2所示。

图2 二维翼型绕流［25］Fig.2 Flow around two-dimensional airfoil［25］

假设流动的总焓不发生改变，构建式（14）：

式中：H0为流动的总焓；h为当地流体的焓值。

在式（14）中，给定远场的静温，求解建立温度场，结合理想气体假设下的状态方程对压力场修正，并计算出当地密度：

壁面参数采用式（17）进行修正：

通过式（15）～式（17）对采用离散涡方法计算得到的压力、密度和压力系数进行一次修正。

对二维不可压N-S方程的动量守恒式进行简化得到：

式（18）等价于：

在壁面处通过涡量对壁面压力进行计算，而在壁面处u=0，式（2）简化为

涡量流量为-ν∇ω，通过边界的流量为，即壁面涡量的生成比率。

壁面的法向和切向用s和n表示，在二维表面得到：

由式（19）～式（21）得到：

即：

因为在每个时间步均有新涡生成，所以上述环量生成率是已知的。根据梯形法则，由式（23）得到：

式中：pj+1和p j为沿壁面单元节点上的压强；和为在一个时间步长Δt下的新生涡环量。

通过式（15）对壁面压强进行修正，利用式（24）中修正的壁面压强对涡团环量进行修正，将修正后的涡团环量分别代入式（8）和式（9），求解壁面涡量和速度场，最后通过式（14）～式（17）对压强、密度和压力系数再次修正。

2 算例验证与分析

2.1 离散涡结冰模拟结果及验证分析

基于离散涡方法进行不可压流动下的结冰数值模拟，以NLF-0414翼型和Business Jet机翼翼型作为算例并与实验值［26］进行对比，以验证本文中未采用压缩性修正下算法和程序的正确性。主程序是将翼型表面分成220块主板块，每块主板块分成5个次板块。由于前缘易发生流动分离，壁面离散单元在前缘加密。计算流程如图3所示，单次结冰模拟时长为60 s，若未满足最终结冰模拟时长则再循环计算。

图3 结冰数值模拟计算流程Fig.3 Workflow of icing numerical simulation

NLF-0414翼型计算条件如下：来流速度66.9 m/s，攻角0.3°，弦长0.9 m，静压93 182 Pa，水滴直径20μm，液态水含量0.44 g/m3，来流静温-15.0℃（霜冰状态），结冰总时长240 s。该工况下结冰模拟结果与实验值［26］的对比如图4所示。

图4 NLF-0414翼型结冰冰形比较Fig.4 Comparison of ice shape for NLF-0414 airfoil

Business Jet机翼翼型计算条件如下：来流速度90 m/s，攻角6.1°，弦长0.9 m，静压94 562 Pa，水滴直径20μm，液态水含量0.54 g/m3，来流静温-5.0℃（明冰状态），结冰总时长360 s。该工况下结冰模拟结果及与实验值［26］的对比如图5所示。

图5 Business Jet机翼翼型结冰冰形比较Fig.5 Comparison of ice shape for Business Jet wing airfoil

在未考虑压缩性修正的情况下，基于离散涡方法的结冰模拟结果与实验结果轮廓贴近，说明采取的方法在不可压流动下可以较好地模拟翼型结冰过程，从而证明了本文求解流程和计算程序的正确性。然而，当前方法未考虑流动的压缩性，可以在此基础上进行压缩性修正，以提高程序在可压缩流动下的预测能力。

2.2 引入修正后离散涡结冰模拟结果及验证分析

选取NACA 0012翼型和Business Jet机翼翼型进行计算验证。在选取相同结冰模型的情况下，分别采用未考虑压缩性修正的离散涡方法（DVM without correction）、引入压缩性修正的离散涡方法（DVM with correction）以及经典CFD方法（Classical CFD）计算流场并最终得到结冰冰形。

算例1：NACA 0012翼型，来流速度129 m/s，攻角4°，弦长0.3 m，静压90 750 Pa，水滴直径20μm，液态水含量0.5 g/m3，来流静温-12.6℃，结冰总时长120 s。

计算结果与文献［27］中的Lewice软件模拟结果和实验值进行比对，如图6所示，可以看出：基于离散涡方法引入修正前后结冰模拟结果上冰角接近，下冰角处有差异，考虑压缩性修正的离散涡方法最终冰形的下冰角处与经典CFD方法结果、文献结果均符合的较好。

图6 算例1结冰冰形比较Fig.6 Comparison of ice shape for Case 1

通过单次结冰模拟结果分析离散涡方法引入压缩性修正对流场和结冰的影响，该算例第60 s结冰模拟结果如图7所示。

图7 算例1第60 s结冰冰形比较Fig.7 Comparison of ice shape for Case 1 in 60 s

从图7可以看出：第60 s未考虑压缩性修正的离散涡方法与经典CFD方法的最终冰形轮廓不一致。相比未考虑压缩性修正的离散涡方法，引入修正后方法与经典CFD方法的最终冰形更符合。

在计算第60 s冰形时需要初始干净翼型的流场数据。计算算例1的干净翼型时，流场压力系数对比如图8所示，流场马赫数对比如图9所示，可以看出：修正后离散涡方法与经典CFD方法计算得到的压力系数和马赫数更加接近，修正后离散涡方法可以较好地模拟流场。

图8 算例1流场压力系数对比Fig.8 Comparison of pressure coefficient in flow field for Case 1

图9 算例1流场马赫数对比Fig.9 Comparison of Mach number in flow field for Case 1

分析离散涡方法引入修正前后对结冰模型计算过程的影响，算例1第60 s通过结冰模型得到的表面水含量如图10所示，横坐标为沿驻点两侧翼型弧长s与弦长c比值的无量纲距离，0为驻点位置；纵坐标表面水含量表示单位表面积上撞击和溢流入当前单元的质量流量与流出蒸发的质量流量的差值。

图10 算例1第60 s表面水含量对比Fig.10 Comparison of water quality on surface for Case 1 in 60 s

算例1第60 s通过结冰模型计算得到的冻结系数如图11所示。

图11 算例1第60 s表面冻结系数对比Fig.11 Comparison of freezing factor on surface for Case 1 in 60 s

从图10～图11可以看出：引入压缩性修正后经结冰模型得到的表面水含量和冻结系数与经典CFD方法经结冰模型得到的结果更为接近，说明引入修正后的离散涡方法相比未考虑修正的方法模拟结冰过程更准确。

算例2：Business Jet机翼翼型，来流速度129 m/s，攻角1.5°，弦长0.9 m，静压89 040 Pa，水滴直径20μm，液态水含量0.31 g/m3，来流静温-5℃，结冰总时长360 s。

计算结果与文献［26］中实验结果对比如图12所示，可以看出：未考虑修正的离散涡方法最终模拟冰形与实验值相比，最终冰形轮廓相差过大，而引入修正后最终模拟冰形与实验值贴近。

图12 算例2结冰冰形比较Fig.12 Comparison of ice shape for Case 2

算例2第60 s结冰模拟结果如图13所示，可以看出：引入修正后最终模拟冰形与经典CFD方法的结果更符合。

图13 算例2第60 s结冰冰形比较Fig.13 Comparison of ice shape for Case 2 in 60 s

计算算例2采用干净翼型时，流场压力系数对比如图14所示。

图14 算例2流场压力系数对比Fig.14 Comparison of pressure coefficient in flow field for Case 2

从图14可以看出：相较修正前离散涡方法计算结果，修正后计算得到的压力系数与经典CFD方法结果更加接近。

将算例2第60 s通过结冰模型得到的表面水含量和冻结系数进行对比，如图15～图16所示，可以看出：引入修正后经结冰模型得到的表面水含量和冻结系数与经典CFD方法经结冰模型得到结果更为接近；未考虑修正时经结冰模型计算出的表面水量在驻点附近基本没有，此时驻点附近无流入水量，而撞击到壁面的水量全部流入附近的单元内，最终冻结系数在驻点附近接近为0。引起这一问题的主要原因是未考虑压缩性的离散涡方法在结冰模型中壁面单元附近的壁温与引入修正后结果有明显差异，未修正方法计算驻点附近的壁温大于0℃，液态水无法冻结而流向其他单元。而引入修正后方法与经典CFD方法在结冰模型中的驻点附近壁温小于或等于0，使液态水在驻点附近结冰。