结冰数值模拟中网格收敛性验证

2020-04-08李佩哲王福新

科学技术与工程 2020年3期

陈航, 谢亮, 李佩哲, 王福新, 刘洪

(上海交通大学航空航天学院，上海 200240)

飞机穿越含有过冷水滴的云层时，撞击在机身表面的水滴会凝结成冰。积冰会改变飞机的气动外形，恶化机翼的气动特性，严重危及整架飞机的飞行安全。随着航空工业的发展，飞机性能越来越先进，但由结冰引发的飞行事故仍然屡见不鲜，必须对飞机结冰问题予以高度重视[1]。

20世纪80年代，欧美等国家开始采用数值模拟方法研究飞机结冰问题。结冰数值模拟是随着计算流体动力学的发展而兴起的，具有高效率、低成本、适用范围广等优点。近年来，随着对飞机结冰物理过程认识的深入和数值计算方法的发展，数值模拟方法在研究飞机结冰问题中的应用越来越广泛。现已开发出数个具有实用价值的结冰预测软件，其中最为著名的是美国航空航天局开发的LEWICE软件[2]。LEWICE软件以及法国航空研究中心的结冰计算软件[3]被美国联邦航空局和欧洲联合航空局所采用，作为确认飞机是否满足结冰适航要求的工具。中国的许多团队在这一领域也开展了相应研究。空气动力研究与发展中心易贤等[4]、南京航空航天大学孙志国等[5]、西北工业大学周志宏等[6]、北京航空航天大学杨胜华等[7]和常士楠等[8]均采用数值模拟的手段对飞机结冰问题开展了深入研究。

在飞机结冰数值模拟中，保证网格收敛性是冰形计算结果可信的必要条件。对此问题，中国学者们已经进行了探讨。蒋胜矩等[9]将二维可压黏流Navier-Stokes方程的解用于翼型结冰的模拟中，对霜冰和明冰冰形进行了计算，并使用三套疏密程度不同的网格对霜冰冰形进行模拟，分析了霜冰情况下的网格收敛性；朱程香等[10]也验证了霜冰情况下的网格收敛情况。前人对网格收敛情况的讨论主要集中在霜冰情况，得出目前算法可以做到网格收敛的结论，然而对于情况更为复杂的带角状冰的混合冰及明冰情况，当前算法是否可以得到网格收敛的解还未见系统的讨论。

由于角状冰对气动特性影响更为显著，因此研究目前结冰算法对于角状冰的网格收敛性能力尤为重要。现基于目前普遍使用的算法开发仿真程序，并使用温度从低到高的四种算例详细定量地验证当前结冰模拟算法的网格收敛性情况。

1 结冰数值模拟方法

结冰数值模拟流程如图1所示。在结冰数值模拟中，首先需要计算空气流场，之后根据流场状况进行水滴撞击特性计算，使用拉格朗日法计算所有水滴粒子的运动轨迹从而得出壁面的水滴收集系数，再使用考虑粗糙度影响的S-A(spalart-allmaras)湍流模型计算壁面热流，结合结冰热力学模型计算得出壁面网格的结冰量，以此进行边界重构并通过网格变形的方法更新网格，重复这个过程直到达到指定的结冰时间为止[11]。

图1 结冰数值模拟流程

1.1 空气流场计算

根据质量、动量、能量守恒关系建立二维非定常Navier-Stokes方程，控制方程为

(1)

式(1)中：Q为守恒向量；∂Ω为区域Ω的边界；Fc表示对流矢通量；Fv表示黏性矢通量；n表示边界的外法向量。

采用有限体积法离散，对流通量计算采用Roe格式，黏性通量计算采用中心格式，时间推进采用LU-SGS(lower-upper symmetric gauss-seidel)隐式时间推进方法，湍流模型采用S-A(spalart-allmaras)模型。

1.2 水滴撞击特性计算

使用拉格朗日法建立水滴运动轨迹方程[12]，控制方程如式(2)所示：

(2)

式(2)中：μ为动力黏度；ud为水滴速度；ua为流体速度；ρd为水滴密度；ρa为空气密度；g为重力加速度；D为水滴直径；CD为阻力系数；Re为相对雷诺数，表达式为

(3)

阻力系数CD采用圆球阻力模型计算，表达式为

(4)

通过求解水滴运动轨迹方程，可以得出水滴撞击到壁面的位置，从而求得壁面的水滴收集系数。

1.3 结冰量计算

1.3.1 结冰热力学模型

结冰计算基于经典的Messinger模型[13]，将结冰表面划分为若干个控制体，对每一个控制体建立质量和能量守恒方程以求解结冰量。质量守恒方程和能量守恒方程分别如式(5)、式(6)所示[14]。

根据质量守恒，控制体内结冰的水量等于进入该控制体的水量减去离开该控制体的水量，即：

mice=mim+min-mout-mev

(5)

式(5)中：mice为结冰的水量；mim为撞击到该控制体内的水量；min为从上一控制体流入的水量；mout为流出到下一控制体的水量；mev为蒸发或升华损失的水量。

根据能量守恒

eim+ein+qice=eout+eev+qc

(6)

式(6)中：eim为撞击到控制体内的水带来的能量；ein为从上一控制体流入的水带来的能量；qice为结冰释放的潜热；eout为流出到下一控制体的水带走的能量；eev为蒸发或升华的水带走的能量；qc为对流换热带走的热量。

1.3.2 对流换热系数计算

在能量守恒方程中，最重要的热流项是对流换热项，因此，对流换热系数的计算是结冰模拟中的关键问题。使用CFD (computational fluid dynamics)方法求解对流换热系数通过湍流模型来完成，其中，Spalart等[15]提出的S-A模型具有计算简单、与试验结果吻合较好等优点。

S-A湍流模型输运方程为

(7)

湍动黏度μt由式(8)求得：

(8)

式(8)中：fv1为黏性阻尼系数，定义为

(9)

式(9)中：Cv1为常数；χ为脉动速度与分子黏度的比值，即：

(10)

雷诺应力由式(11)计算：

(11)

式(11)中：ρ为流场当地密度；Sij为速度散度，ui、uj分别代表i、j方向的速度。

(12)

(13)

式中：d为到最近壁面的距离；S为涡度，定义为

(14)

式(14)中：Ωij为时均旋转率张量，定义为

(15)

输运方程中fw函数为

(16)

式(16)中：Cw3为常数；g的表达式如式(17)所示：

g=r+Cw2(r6-r)

(17)

式(17)中：r为脉动速度与到最近壁面距离关系的函数，其表达式如式(18)所示：

(18)

壁面边界条件为

(19)

上述S-A湍流模型适合于光滑表面的湍流模拟。在结冰计算中，翼面由于积冰的影响会变得粗糙。由于表面粗糙组织对气流的扰动作用，流动转捩的位置将大大前移，壁面与气流间的对流换热也将被大大强化，因此，在结冰计算中必须考虑粗糙度对壁面对流换热的影响。

Aupoix等[16]基于等效沙粒粗糙高度模型对S-A模型进行了扩展，得到了考虑粗糙度影响的S-A湍流模型的扩展模型。粗糙壁面的扩展模型在原始湍流模型的基础上作如下改动。

向到最近壁面的距离d中加入了粗糙度的影响，d变为

d=dmin+0.03ks

(20)

式(20)中：ks为等效沙粒粗糙高度。

(21)

式(21)中：n表示边界的外法向量。

调整χ的值以改变黏性阻尼函数fv1，χ改为

(22)

式(22)中：常数Cr1=0.5。

fv2函数变为

(23)

使用CFD方法求解流场后，在壁面附近，根据傅里叶定律可以求得壁面对流换热量q：

(24)

根据对流换热系数定义即可求出粗糙表面对流换热系数：

(25)

式(25)中：Twall为壁面温度；T∞为环境温度。

1.3.3 结冰表面粗糙度计算

结冰表面粗糙度的产生机理非常复杂，与水膜的流动稳定性、水滴及气流的扰动等很多因素都有关系，是一个耦合的多尺度问题，很难用解析方法模拟其特征，因此在结冰数值模拟中一般不对其做数值预测，而是利用试验得到的经验公式进行计算。采用Shin等[17]提出的等效沙粒粗糙高度模型计算粗糙度，将等效沙粒粗糙高度ks与结冰气象条件中的液态水含量(liquid water content，LWC)、环境温度T∞和水滴粒径(median volume diameter，MVD)关联起来，其中：

(26)

式(26)中：kLWC、kT∞和kMVD分别为表征含水量、环境温度、水滴粒径影响的函数，其表达式分别为

kLWC=0.571 4+0.245 7LWC+1.257 1LWC2

(27)

(28)

(29)

(30)

式(30)中：c为弦长。

1.4 网格变形

随着冰的增长，外部流场和水滴收集系数等参数会随着结冰外形变化而变化，为保证计算结果的准确性，需要对计算网格进行更新。网格变形技术实际上是将初始状态下壁面的位移量通过插值的方法均匀地分布到流场内部的节点中，达到网格更新的效果。采用网格变形技术能够最大限度地保留初始网格特征，在网格更新速度上比重新生成网格要快得多，而且，利用上一次结冰计算的流场值作为新网格流场计算的初始值，可以大幅减少网格更新后的流场计算时间[18]。

在网格变形开始之前需要判断冰的生长方向，选取的冰生长方向如图2所示[19]。