基于Myers模型的三维结冰数值仿真

2018-09-29雷梦龙常士楠杨波

航空学报 2018年9期

雷梦龙，常士楠，杨波

北京航空航天大学航空科学与工程学院，北京 100083

近年来，大量的飞行事故导致飞机的结冰研究成为一个热点。根据最近的飞行数据，9%的飞行事故是由于飞机在严酷天气下结冰所致，所有气象因素带来的飞行事故中，12%是飞机结冰导致的，且92%是在飞行中发生结冰[1]。

飞机上结的冰，就其外形可以分为明冰、霜冰、混合冰3种。霜冰与机翼外形较为接近，对机翼的气动性能影响较小。明冰是在气温较高、液态水含量较大时，撞击到气动面上的水滴没有立即全部结冰，而是产生水膜沿机翼表面流动，形成溢流冰。严重情况下可能形成冰角，破坏气动性能[2]。目前对机翼结冰问题的研究方式主要分为理论、实验和数值仿真，数值仿真作为一种经济快速的模拟方法，一直是结冰领域最常用的方式之一。

机翼表面结冰问题的数值仿真过程中，由于复杂的结冰机理及界面处气-液-固耦合传热传质特性，准确预测明冰表面的水膜流动一直是一个难点[3]。早期对飞机结冰问题的研究中，通常使用Messinger结冰模型[4]求解溢流水和结冰量，如现在仍在使用的LEWICE[5]和ONERA[6]等。Messinger模型引入了冻结系数，在一个控制体内，撞击到机翼表面的水一部分会冻结成冰，未冻结的水沿气流方向全部流入位于下游的控制体。Messinger模型相对于实际结冰过程过于简陋：① 没有考虑到空气压力梯度和剪切力对水膜流动的影响，它假设某个控制体中没有冻结的水全部流入下游的控制体中，与实际情况不符，尤其是在有回流流场的冰角位置；② 在水膜结冰过程中，没有考虑空气、水膜、冰层之间的导热；③ Messinger模型在提出之初是用于计算二维结冰问题的，虽然被扩展到三维结冰情况，但在三维情况下，Messinger模拟更为困难，每个控制体的上游和下游都不止一个控制体，水流流入和流出的方向和流量难以准确确定。国内应用Messinger模型计算三维结冰的文献已经有很多，如文献[7]计算了MS-317后掠翼的结冰，文献[8]计算了发动机进气道的结冰，文献[9]计算了转子的结冰。

为了克服Messinger模型的缺点，水膜的流动必须考虑在内，现有的商业软件中，FENSAP-ICE就考虑了水膜流动[10]。此外，Myers提出了适用于任意三维表面的水膜流动模型[11]，并基于此流动模型，提出了Myers结冰模型[12]。Myers模型在计算明冰的结冰增长时，考虑了水膜、冰层之间的导热，在计算水膜流动时，考虑到了空气剪切力、空气压力梯度和曲面特性[13]。目前对Myers模型的研究还较少，其中，比较深入的有文献[14-15]。

对机翼结冰的数值仿真通常是单步法，在一次计算空气和水滴流场的基础上，认为机翼结冰后的冰形对空气和水滴流场没有影响，不再根据新的冰形进行流场计算。单步法需要的计算量更少，且算法更简单，因此大部分结冰计算都使用单步法计算结冰[16]。Myers模型计算结冰[17]时，需要考虑此处的结冰厚度、水膜厚度等信息，但在冰形改变后重新计算结冰时，机翼上的结冰厚度和水膜厚度会初始为0，因此Myers模型在使用多步法计算时需要复杂的算法。

为了充分考虑水膜流动对结冰的影响，提高结冰计算的准确性，本文采用了Myers结冰模型，并在其基础上加以改善，采用单步法和多步法计算结冰冰形。为了验证结冰模型的正确性，本文将NACA0012平直翼和GLC-305后掠翼的计算结果与商业软件以及实验结果对比，并提供了三维发动机进气道的计算结果。本文分析了仿真结果的合理性，并讨论了Myers结冰计算模型的优势和不足之处。

1 Myers水膜流动模型

水膜在曲面上的流动受到空气压力梯度、重力分量、空气剪切力等影响，为了简化水膜流动模型，需要进行以下假设：

1) 水的黏度、密度在不发生相变时，随温度变化很小，在计算水膜流动时将其视为常数。

2) 忽略机翼表面粗糙度对水膜流动的影响。

图1为任意三维曲面上的贴体坐标系示意图，贴体坐标系用(s1，s2，η)表示，其中s1、s2与曲面相切，η是曲面在该点的法线，不同控制体的s1、s2和η的方向不同。

曲面可以用第1类和第2类基本量描述，对于不同类型的曲面，曲面基本量的取值不同，实际计算中，可以假设曲面为平板来简化计算[18]，如ICECREMO软件就是以平板假设简化水膜流动模型。Myers等以润滑理论简化水膜在平板上流动的Navier-Stokes控制方程，得到[19]

(1)

(2)

(3)

式中：u为s1方向的速度；v为s2方向的速度；μw为水的黏度；p为空气压力；ρw为水的密度；ε为水膜高度和长度的比值；Re为雷诺数；gs1、gs2、gη分别为重力在s1、s2和η方向的分量。根据润滑理论，O(ε,ε2Re)的数值大小可以忽略。

水膜在三维表面上s1和s2方向单位长度的体积流量为水的流动速度在高度上的积分，即

(4)

(5)

式中：h为水膜厚度；b为冰层厚度；As1、As2为空气剪切力在s1、s2方向的分量。水膜流动的质量守恒方程为

(6)

(7)

式中：LWC为液态水含量；V∞为空气远场的速度；β为水滴收集系数；me为水膜的蒸发质量;ρi为冰的密度。

2 Myers结冰增长模型

飞机表面的结冰厚度和水膜厚度相比于机翼的整体尺寸通常很小，因此，Myers模型在分析机翼表面、冰层、水膜之间的导热时，将结冰过程视为准静态过程，忽略s1、s2方向的导热，只考虑与机翼表面垂直的η方向导热，并且忽略机翼与其表面冰层或水膜之间的接触热阻，认为机翼表面温度与冰层或水膜底部温度相同。

图 2为机翼表面的冰层结构，从下到上依次为冰层、水膜、空气。冰层的底部温度为机翼表面温度Ts，顶部温度为水的凝固点Tm，水膜底部温度为Tm，顶部温度为Th，空气温度为Ta。通常一个控制体内的热流分为向空气的对流换热热流Qc、蒸发热流Qe、升华热流Qs、水滴动能Qk、结冰释放的潜热Ql、水滴显热Qd、气动加热Qa、水膜流动带入的热流Qin、水膜流出带走的热流Qout、冰层内的导热热流Qcond。

Qc=hcv(Tm-Ta)

(8)

Qe=mele

(9)

Qs=msls

(10)

(11)

Ql=micelf

(12)

Qd=mimpcw(Td-Tm)

(13)

(14)

Qin=mincw(Tin-Tm)

(15)

Qout=moutcw(Tout-Tm)

(16)

(17)

式中：hcv为对流换热系数；me为蒸发质量；ms为升华质量；mimp为撞击的水滴质量；mice为结冰质量；vd为水滴速度；le、ls、lf分别为蒸发、升华、结冰潜热；cw为水的比热；Td为水滴温度；r为恢复系数；ue为边界层速度；ca为空气的比热；Tin为流入温度；Tout为流出温度；ki为冰的导热系数;min和mout为流入和流出质量。对于霜冰和明冰，热流项的个数以及其中一些参数的取值略有不同。

2.1 霜冰结冰模型

结冰开始阶段，过冷水滴撞击到机翼表面时，由于机翼表面温度低于水的凝固点，结冰类型为霜冰。结霜冰时，机翼表面的热流平衡方程为

Qk+Qa+Ql-Qc-Qd-Qs-Qcond=0

(18)

在每一个时间步长对式(18)迭代求解可得出霜冰结冰情况下机翼表面的温度Ts。

霜冰的结冰量为撞击到机翼表面的水量，此结冰增长微分方程为

(19)

将结冰过程视为准静态过程，忽略s1、s2的导热和冰层本身的储热能力，根据控制体内的热流平衡，可得霜冰表面的导热微分方程为

(20)

式(20)的右侧是空气向机翼表面方向传热的热流密度，将其写为E0d-E1d·Ts的形式，即

E0d=LWC·V∞·β·[cw(Td-Tm)+ciTm]+

(21)

E1d=LWC·V∞·β·ci+hcv

(22)

解式(20)，可得出冰层内的温度分布为

(23)

霜冰顶部温度为Tm时对应的冰厚为临界厚度bf，超过bf，冰的表面会出现水膜，霜冰转化为明冰。令式(23)中T=Tm、η=b，可得临界冰厚bf为

(24)

若bf>0且结冰厚度b≤bf，冰层的表面温度小于水的凝固点Tm，此时结冰类型为霜冰；当结冰厚度升至b>bf时，冰层的表面温度等于Tm，结冰类型为明冰；若bf<0，则无论结冰厚度为多少，冰层表面温度都小于Tm，结冰类型为霜冰。

2.2 明冰结冰模型

明冰结冰时，水膜除向空气导热外，还向冰层导热，机翼表面每个控制体内的热流与霜冰情况有所不同，其热流平衡方程为

Qk+Qa+Ql+Qin-Qc-Qd-Qe-

Qcond-Qout=0

(25)

在每一个时间步长内解此方程，可得出明冰结冰情况下的机翼表面温度Ts。

Myers模型认为，只要冰厚b≤bf，结冰类型就是霜冰。但是实际情况下，若水的流量过大，结冰释放的潜热不足以使水膜全部结冰时，即使冰厚b≤bf，也可能产生水膜，继续向后流动。因此，本文在计算结冰时，将明冰分为两种类型。

当临界冰厚bf>0且冰厚b

根据控制体内的热流平衡，水膜与空气接触面的导热微分方程为

(26)

式(26)右侧是空气向水膜方向传热的热流密度，将式(26)的右侧写为E0f-E1f·Tm的形式，可以得出水膜内η方向的温度分布，即

mincwTin-moutcwTout

(27)

E1f=LWC·V∞·β·cw+hcv+mincw-moutcw

(28)

导热微分方程式(26)的解为水膜内的温度分布，即

(29)

冰层的温度分布为

(30)

忽略s1、s2方向的导热，只考虑η方向的导热，水膜的结冰速率微分方程为

讯问笔录通常是由监察机关调查人员制作，内容通常是一问一答的记录双方的语言，笔录制作本身可能带有倾向性，讯问笔录不能完全真实地反映讯问过程中发生的情况；此外，被调查人的伤痕也并不意味着他一定受到了刑讯逼供，造成伤痕的可能性很多，比如打架斗殴，甚至不排除被调查人以自伤的方式来诬陷调查人员，况且非法取证手段多样，并不一定会在身体上留下痕迹，若是采取精神摧残或者恶劣讯问环境的方式，很难证明当时的真实情况。录音录像不仅可以反映讯问时双方的问答，还可以记录当时所发生的一切情况，可以起到证明是否有非法取证行为的作用。

(31)

本文对Myers模型的改进主要是对结冰类型判断标准的修改。Myers模型认为，水滴撞击到机翼表面发生结冰，若结冰厚度b≤bf，结冰类型为霜冰，结冰厚度增长至b>bf后，结冰类型转变为明冰。此推论是假设结冰过程为准静态得出的，但是实际结冰情况下，即使结冰厚度b≤bf，若上游流入和水滴撞击的流量Q足够大，也有可能不会全部冻结，从而使结冰类型由霜冰转化为明冰。因此，本文对结冰类型判断依据进行修改，结霜冰时，判断流量Q是否超过最大结冰质量，若超过最大结冰质量，则水不会全部冻结，冰的表面会产生水膜，结冰类型转变为第2类明冰。令式(31)中的水膜厚度h=0 m，此时的结冰量即是当前位置的最大结冰量。

为了比较本文对Myers模型的改进对结冰冰形的影响，本文提供了改进的Myers模型与原Myers模型在两种不同条件下的计算结果对比，两个计算结果都使用单步法。

表1给出了NACA0012翼型的环境条件，图3是Case 1的结冰冰形计算结果，图 4是Case 2的结冰冰形计算结果。从图3和图4可以明显看出，本文改进的Myers模型与原Myers模型计算结果的区别主要体现在结冰的上下极限附近，这两个位置的结冰量主要来自上游流入的水膜。由于其结冰厚度没有超过临界结冰厚度，原Myers模型将这两个位置处的冰作为霜冰计算，积聚在这里的水膜会全部冻结，形成冰角，造成此处的结冰量过大；而本文改进的Myers模型会通过计算这两个位置处的水膜流量大小，将其分类为完全冻结的霜冰和部分冻结的明冰区别处理，因此结冰冰形更符合实际情况。由于采用单步法，结冰冰形会与多步法计算结果产生较大差距，多步法的计算结果将在4.1节给出。

表1 NACA0012翼型结冰条件Table 1 Icing condition of NACA0012 airfoil

图 5和图 6分别是Case 1和Case 2的水膜厚度计算结果对比。以机翼驻点为分界线，机翼上部和下部各有3个水膜厚度峰值。在驻点处，水滴大量撞击，水膜厚度逐渐增加，形成第1个厚度峰值，在空气压力和剪切力作用下，水膜流速加快，厚度减小；水膜流量继续累积，厚度再次增加，形成第2个厚度峰值，水滴撞击量开始小于冻结量，水膜厚度再次减小；最后，空气压力梯度和剪切力减小，使水膜流速减缓，厚度增加，形成第3个厚度峰值。Case 1相比于Case 2，由于气温较高，水膜覆盖的区域更大。原Myers模型对明冰的判断标准过于简单，把上下结冰极限处的冰归类为霜冰，使水膜全部冻结为冰，因此水膜覆盖范围比改进的Myers模型小，结冰厚度更大。

3 数值计算方法

数值计算的第1步是确定s1、s2和η的方向，本文中，η的方向为网格的法向，指向机翼外部；s1的方向为网格表面空气流动方向在网格面上的投影，在此方向上压力梯度和空气剪切力很大，因此，计算水膜流动时，s1的方向为水膜流动的主要方向；s2的方向为s1和η方向的垂向。

(32)

式中：t为时间；A为网格的面积；n为网格边的法向，指向网格外，方程右边中括号内的各项依次为流入流出项、冻结项、撞击项。

在结冰初始阶段，结冰厚度为0 m，此时的结冰类型为霜冰，霜冰结冰的数值计算式为

(33)

无论是第1类明冰还是第2类明冰，冰层顶部温度都等于水的凝固点Tm。明冰结冰量的数值计算式为

(34)

在单步法的计算中，需要在从计算开始到计算结束的时间内，以一定的时间步长Δt为间隔，计算结冰量和水膜厚度，最终得出结冰的冰形；在多步法的计算中，需要每隔一定的时间间隔更新冰形及其流场，并重复上述计算步骤。

4 计算条件及结果

本文使用的流场数据均由流体计算软件Fluent得到。流场数据需要包含机翼表面的坐标信息、空气压力、空气剪切力、水滴撞击系数、对流换热系数等数据。然后，用本文改进的Myers模型结冰计算程序读取数据并求解。

4.1 NACA0012平直翼

本文使用的是弦长为0.533 4 m的NACA0012平直翼，水滴直径为20 μm，速度分别为60.1、102.8 m/s，空气液态水含量分别为1.3、0.55 g/m3，空气温度分别为266.45、262.04 K，结冰时间分别为480、420 s，如表 1所示。NACA0012翼型对比使用的商业软件计算结果与实验结果来自文献[20]。

两组计算条件均为明冰结冰类型，不同的是Case 1的气温高，液态水含量高，结冰时间长，相比于Case 2会产生更多水膜在气流作用下流动，形成明显的溢流冰。在结冰计算中，结冰造成的机翼形状变化会对流场产生不可忽视的影响，因此在整个过程中，本文使用多步法根据冰形变化重新计算流场。

图 7为Case 1的仿真结果与实验结果及ONERA软件仿真结果的对比。在结冰的冰角处，本文的结冰极限与实验结果及ONERA软件仿真结果非常接近。上结冰极限的下游处没有结冰是因为此处的空气回流使得水膜在压力梯度的作用下不会向下游位置流动。因此，上冰角的位置可能是由于流场计算软件对空气流场计算结果的差异所致。在下结冰极限位置处，本文的计算结果与实验结果十分吻合。在结冰的上半部分，本文的结冰量略小于实验结果，原因可能是Myers模型在计算水膜流量时，把机翼表面作为光滑表面。整体上本文计算的结冰量与实验结果相近。

图 8为Case 2的仿真结果与实验结果及LEWICE软件计算结果的对比。相比于实验结果和LEWICE软件计算结果，本文计算的结冰冰形较为圆润，原因可能与Case 1相同，是由于没有考虑粗糙度对水膜流动的影响所致。本文计算的上冰角的位置以及下结冰极限与实验结果吻合很好。

本文提供了单步法和多步法(5步)的计算结果。多步法的上下结冰极限处的结冰厚度较低，原因可能是因为随着结冰的进行，冰形对流场的影响逐渐突显出来，在上下结冰极限处，受到结冰冰形的影响，水滴收集量相比于没有结冰的机翼要小，因此结冰厚度也减小。在水滴的主要撞击位置，由于空气-水滴流场受到结冰冰形的影响，多步法的结冰厚度稍高于单步法的结冰厚度。单步法相对于多步法的优势不明显，原因可能是计算网格在每一步的网格重构时会产生误差，随着结冰过程的进行，这些误差累积起来，最终对冰形产生较大影响。

4.2 GLC-305后掠翼

本文的目的是对三维机翼的结冰进行仿真。为此，选取三维GLC-305后掠翼。表 2为GLC-305后掠翼的计算条件，大气压力为101 300 Pa，水滴直径为20 μm。图 9为GLC-305后掠翼的几何外形，RLE为可变机翼前缘，MAC是平均气动弦长，截面 A、截面 B与机翼前缘垂直，截面 C为机翼翼根,图中1 in=25.4 mm。图 10为GLC-305后掠翼的表面网格。本文对比使用的实验结果和LEWICE计算结果来源于文献[21]。

表2 GLC-305后掠翼的结冰条件Table 2 Icing condition of GLC-305 swept wing

GLC-305后掠翼的结冰计算条件为明冰条件，从计算结果中可明显观察到明冰的特征。图 11为翼根处的结冰厚度云图，图 12为翼梢处的结冰厚度云图。翼梢处的结冰厚度比翼根处更大，原因是翼梢处的水滴收集系数和对流换热系数比翼根处大，因此结冰情况更严重。

本文将计算结果与实验结果以及LEWICE软件计算结果对比。图 13(a)为截面A处的结冰冰形对比，图 13(b)为截面B处的结冰冰形对比，图 13(c)为截面C处的结冰冰形对比。在截面A处，本文计算的冰角不明显，但结冰量与实验结果大致相同。在截面B处，本文计算的上冰角稍小于实验结果，但是对于上下结冰极限以及结冰量总体趋势，本文与实验结果符合得很好。在截面C处，本文的计算结果与实验结果十分接近。在截面A和截面C中，LEWICE过度预测了冰角的大小，并错误预测了冰角的位置。

4.3 发动机进气道

发动机进气道结冰也是威胁飞机飞行安全的重要问题。发动机进气道在进口部分通常做成翼型形式，其内表面结冰范围和结冰强度通常比外表面大得多，原因是内表面气动特性恶化，速度场分布不均匀以及气流局部分离。

本文计算了发动机进气道在明冰结冰情况下的结冰冰形，结冰条件如表 3所示，空气压力为101 325 Pa，水滴直径为20 μm。发动机进气道的外形和表面网格如图 14所示。

本文选取的结冰条件较为严重，速度为150 m/s，温度为263 K，液态水含量为2 g/m3，属于明冰结冰条件。图 15为发动机进气道表面的水滴收集系数(β)云图，图 16为发动机进气道的结冰厚度云图，图 17为发动机进气道顶部截面的结冰外形，图 18为底部截面的结冰外形。可以看出，本文中发动机进气道处的冰形有明显的溢流冰特征。结合图17和图18可以发现，发动机进气道的结冰主要集中在进气道的内侧，而在外侧结冰量就少了很多。原因是水滴主要撞击在内侧，在外侧撞击量少很多，造成进气道内部的水滴收集量相对外部较大。