CN102682145A - 飞行结冰的数值模拟方法 - Google Patents

Abstract

本发明是一种用于飞行器飞行结冰的数值模拟方法，来模拟飞行器在空中飞行期间遭遇结冰时的状态。主要特征是包括用于模拟空气-超冷水滴的运动的两相流流动的单流体模型中壁面水膜表面的速度分解和水膜厚度的算法、在计算冰层形状和内部的温度分布的水膜结冰状态模型中采用网格加密方法追踪结冰界面的算法、基于固定计算网格利用上述模型和算法进行飞行结冰数值模拟计算的流程。

Description

飞行结冰的数值模拟方法

技术领域

本发明涉及航空工程领域，是计算流体力学在航空工程领域的应用，具体地讲，是一种用于模拟飞行器的飞行结冰的数值方法。这种数值模拟技术可以用计算机高级程序语言实现，并通过计算机的运行来模拟飞行器在空中飞行期间遭遇结冰时的状态。

背景技术

飞行器在一定的飞行高度范围内穿过云层时，大气中的超冷液态水滴(温度在冰点以下，但以液态水滴形式存在)会碰撞在飞行器多个部件表面形成水膜，如机翼、机身、驾驶舱、尾翼、发动机进气口等部件的表面都容易形成大量积聚的液态水膜。通常用超冷液态水含量(单位体积内的超冷液态水滴的总重量，量纲kg/m³)来表示结冰条件。如果超冷液态水含量较高，飞行器一些部件表面上积聚的水膜会在结成冰层。这种现象被称为飞行结冰。大量的飞行结冰会增加飞行器的重力、改变重心、改变飞行器外形和表面粗糙度，引起阻力增加、升力减小和失速角减小。同时，结冰会阻碍飞行器表面一些运动部件的功能，如襟翼、平衡器的运动，危害飞行器的稳定性和操纵性。

为保证飞行器在遭遇结冰条件下飞行的安全性，飞行器制造者必须表明飞行器能在各种结冰飞行条件下满足飞行包线以取得适航证。适航取证的过程是通过空中飞行测试、风洞测试、计算机数值模拟三种手段共同实现的。空中飞行是最直接的测试手段，完全在自然条件下进行。但是这种方法不仅成本昂贵，而且自然条件无法完全达到飞行包线上的所有条件，不可能进行逐个条件下的验证。在陆地上唯一替代空中飞行测试方法就是倚重结冰风洞。在风洞中产生高空飞行时遇到的超冷水含量和水滴尺度等结冰条件。结冰风洞制造成本昂贵，而且风洞中无法产生的大气中超冷水滴尺度的分布。此外，由于风洞测试中采用缩小的飞行器几何模型，所以流动雷诺数不准确，以致难以准确预报真实的飞行器飞行结冰状态。

自从上世纪80年代，国际上逐步开始了飞行器飞行结冰过程的数值模拟的研究，不仅在理论上对这一复杂的现象有了深刻的认识，而且已经将数值模拟的结果应用到工程上。飞行结冰的计算机数值模拟已经成为一种新产品开发设计和适航取证的支持工具。先进的数值模拟技术还可以弥补结冰风洞实验中的误差和缺陷。飞行结冰过程的数值模拟技术经过20多年的发展后，出现了具有代表性的相关产品。例如，来自美国的LEWICE软件、来自法国的ONERA软件、来自加拿大的Fensap软件。

飞行器飞行结冰的数值模拟技术中所采用的流程都是基于三个主模块的调用。三个主模块及其各自主要功能是分别为：

空气流动模块：用于求解飞行器外部流场(包括飞行器表面)信息，即求解流体流动(此处为空气运动)控制方程；

超冷水滴运动模块：用于求解大气中的超冷水滴与飞行器表面的碰撞过程，以获得飞行器表面的液态水的状态(用液态水收集率表示)，即求解水滴运动方程；

结冰状态模块：用于求解液态水的结冰过程，获得结冰后的几何形状。

上述各种软件中所采用的数值模拟技术也存在一些差别，主要体现在外部流场求解的方法和水滴运动的描述方法。例如，LEWICE和ONERA软件中利用面元法求(Panel Method)解外部流动，再进行可压缩修正。这意味着将飞行器外部流体流动视为不可压缩势流流动，是对真实情况的近似。尽管FENSAP软件中按照两相流湍流流动求解，但通过大量简化性假设，分别求解空气和水滴的运动方程，仅考虑了空气对水滴的阻力。

其次，在LEWICE软件采用拉格朗日方法为框架描述水滴运动问题，追踪水滴的运动轨迹，在处理复杂几何表面的结冰问题时，受到一定限制。而ONERA和FENSAP软件采用欧拉方法为框架描述水滴运动问题，并将空气和水滴的运动视为两相流体的流动。至于结冰状态模型，各种软件均以著名的Messinger结冰模型为基础。该模型是一个零维模型，认为冰层内部的特性是均等的，从结冰过程的能量守恒形式出发，并结合质量守恒关系，建立了常微分方程。LEWICE和ONERA软件中将结冰过程视为一个准稳定过程，即在一个结冰计算时间间隔内，外界空气和水滴的运动是不变化的。这种准稳定假设在冰层外部的流动有分离的情况下显然是不正确的。Fensap软件中建立了冰层增长的时间相关项，解决了这个问题，但是需要求解额外两个偏微分方程。上述各个软件中均需要在每一个时间段因为冰层形状的改变而进行网格重构。该过程中需要对网格进行局部坐标插值、光顺、正交处理，以提高网格质量，同时，要对网格上的变量进行插值运算。这个过程不仅耗时，而且插值运算会降低整体计算精度。

下面以FENSAP软件模拟飞行结冰的求解过程为例(见图1)，说明该软件中各个模块之间的关系和求解过程。从图1中可见，在某一个时间段首先进行外部流场的空气运动、水滴运动的单独求解，或结合在一起求解。然后将获得的参数，如壁面集水率，壁面剪切力、壁面传热量，带入结冰模型，计算壁面各个点上的冰层的厚度，获得了下一时刻飞行器表面的形状。然后围绕结冰后的飞行器外形将计算网格重构，并进行坐标插值、光顺、正交处理，再进入下一个时间段进行计算。该软件生产商介绍一个二维NACA0012机翼的结冰模拟计算，共47万个网格点，8个CPU，需要3.5个小时完成。其中，网格重构过程占据总体计算时间的15％。此外，这种计算方法在求解空气运动流动的控制方程时，没有考虑水滴对空气流动的影响，很明显，会产生一定误差。

发明内容

本发明是一种用于飞行器飞行结冰的数值模拟方法，这种数值模拟技术可以用计算机高级程序语言实现，并通过计算机运行来模拟飞行器在空中飞行期间遭遇结冰时的状态。该数值模拟方法目的是为更接近真实的飞行条件，兼顾计算精度、效率和功能。本发明提出的方法的主要特征是包括用于模拟空气-超冷水滴的运动的两相流流动的单流体模型中壁面水膜表面的速度分解和水膜厚度的算法、在计算冰层形状和内部的温度分布的水膜结冰状态模型中采用网格加密方法追踪结冰界面的算法、基于固定计算网格利用上述模型和算法进行飞行结冰数值模拟计算的流程。其中模拟空气-超冷水滴的运动的单流体两相流流动的模型是一组描述飞行器外部流场流体运动的偏微分方程组；计算冰层内部的温度分布的模型是一组描述水膜流动和冰层内温度分布和相变的偏微分方程，该方程组的解也可以用来识别结冰界面。本发明中使用单流体两相流流动的模拟方法，即将空气-超冷水滴的两相流动归一为单一物质的流动，只建立一组流体控制方程求解飞行器外部流场，只是在边界结冰的位置进行两相流速度的分解。

大气中的超冷水滴统计平均尺度较小，一般在50μm以下，比较均匀地分布在大气的对流运动中和空气中一起运动，成为一种空气-超冷水滴的混合流体。飞行器外部流场中空气-超冷水滴组成的二元混合物，混合均匀，热力学性质接近。在飞行结冰的数值模拟中，考虑到大气对流运动和飞行器的速度相比是小量，飞行器外部流场中的两相速度滑移是因为飞行器的亚音速飞行产生的向上游扰动，是一个小量。所以，就相当于某种等效的单流体的流动，在一个很小的流体微团内可视为连续介质。等效混合物的物理性质，如密度、比热、粘性系数、导热系数等都可从两种组分中的对应参数按照质量或体积分数进行加权平均获得。例如二元混合物中，体积分数为α、粘性系数μ、密度ρ和速度

分别用下标1表示空气、下标2表示超冷水滴、下标m表示混合物，显然有

α₁+α₂＝1。                (1)

混合物的密度ρ_m和粘性系数μ_m按照体积分数进行加权平均

ρ_m＝α₁ρ₁+α₂ρ₂，        (2)

μ_m＝α₁μ₁+α₂μ₂。        (3)

速度的体积分数进行加权平均分别是

${\stackrel{\→}{U}}_m=\frac{1}{{\mathrm{\ρ}}_m}({\mathrm{\α}}_1{\mathrm{\ρ}}_1{\stackrel{\→}{U}}_1+{\mathrm{\α}}_2{\mathrm{\ρ}}_2{\stackrel{\→}{U}}_2).---\left(4\right)$

所以，空气-超冷水滴的单流体两相流流动控制方程的形式完全与已知的单一组分的可压缩流体的控制方程一样，使得两相流问题简化。两相间的影响(如水滴受空气的阻力和浮力)通过使方程组封闭的体积分数扩散方程体现。同时混合流体的湍流脉动也可以按照单组份流体的湍流模型求得。空气-超冷水滴的单流体两相流流动控制方程包括：两相体积分数的扩散方程、连续方程、动量方程、能量方程。如果来流是湍流，还需要增加湍流模型。上述方程组除了扩散方程，其他方程的形式均与公知的单一组分的可压缩流体的控制方程的形式一致。同时控制方程组的空间、时间离散方法也与单组份流体的一致。空气-超冷水滴的单流体两相流流动的模拟结果给出飞行器外部流场的流动信息，即各个计算网格点上(或是网格单元内部)的空气和超冷水滴的密度、压力、速度，以及由上述独立变量推导出的温度、动力粘度等信息。边界网格内的液态水滴在壁面形成一定厚度的水膜，结冰将首先在水膜和飞行器干净的壁面之间发生，之后将在水膜和已经结成的冰层之间发生。

本发明开始于壁面网格内混合物流动的速度分解，从两相混合物流动中分解出超冷水滴的速度。其原理是：先假设网格内的超冷水滴全部形成假想水膜(实际上只有部分形成水膜)，其厚度由当地的超冷水滴的体积分数α₂决定。而假想水膜的速度就是超冷水滴的速度，所以，求得假想水膜的速度后再按照积分时间计算真实的水膜厚度和速度。图2给出假想水膜表面速度的分解原理图。图2(a)是壁面网格内的空气-超冷水滴的分离和假想水膜形成原理图。以二维流动为例，按照壁面网格中的超冷水滴的体积分数折算成壁面上的假想水膜高度h_f。网格中其余部分是空气，几何中心距离壁面高度是h₁；水膜的几何中心距离壁面高度是h₂。这里的壁面可以指飞行器未结过冰的干净的壁面，也可是已经结成的冰层的表面。形成的假想水膜仍然是流动的，因为流体粘性的作用，在壁面形成边界层。图2(b)不可压缩流边界层示意图。水膜是不可压缩流，按照公知的不可压缩流的边界层理论，流体在壁面上的速度为零，边界层中的速度在X-方向的分布为u(x，y)，在不考虑内部压力梯度时可由公知的Blasius公式给出，即

$\frac{u(x,y)}{U_{\∞}}=\frac{2y}{\mathrm{\δ}\left(x\right)}-{\left(\frac{y}{\mathrm{\δ}\left(x\right)}\right)}^2.---\left(5\right)$

上式中和图2(b)中，x是欲求速度的位置，y是距壁面的高度，U_∞是来流速度。边界层的厚度δ(x)，由下式给出，

$\mathrm{\δ}\left(x\right)=5\sqrt{\frac{\mathrm{xv}}{U_{\∞}}},---\left(6\right)$

其中，v是流体的运动粘度系数。所以，来流速度U_∞和运动粘度v确定后，公式(5)和(6)可以确定边界层中某点(x，y)的速度u的大小。

一种情况是，边界网格中的假想水膜被包含在边界层内，假想水膜和空气共同形成的边界层的速度分布。在假想水膜和空气接触面上，由于两相物性差异，存在滑移速度。为了求得假想水膜表面速度。该值需要从两相流单流体模型中的混合物速度中分解出来。图2(c)给出按照两相流单流体模型的边界层速度分布。其中H是网格高度，h_m是几何中心，u_m是几何中心处的两相混合速度。图2(d)是假想水膜-空气边界层的速度分布。其中，u₁是空气的速度、u₂是假想水膜的速度、u_2f是假想水膜表面的速度，其他符号和图2(a)一致。此外，需要建立假想水膜-空气边界层的速度分布积分与两相流单流体的边界层速度分布积分相等的假设，即图2(c)中的ABC的面积等于图2(d)中的ABCD的面积。目的是为了将假想水膜表面的速度u_2f从混合流动中分解出来。图2(e)给出了水膜表面速度的分解算法流程图，其具体步骤是：

(1)设定一个壁面上超冷液态水的运动粘度v为初始值。

(2)按照不可压缩边界层速度分布公式(5)求得假想水膜表面的速度u_2f。公式(5)中假想水膜厚度h_f代替y，来流速度是无穷远处混合流体的速度，即飞行器的飞行速度。

(3)按照质量平均求液态水边界层的速度u₂。该速度是指h₂处的速度，其表达式为，

$u_2=\frac{1}{2}{u}_f.---\left(7\right)$

(4)求空气的速度u₁。因为混合流体的速度由公式(4)表示，所以单相空气速度为，

$u_1=\frac{{\mathrm{\ρ}}_m{u}_m}{{\mathrm{\α}}_1{\mathrm{\ρ}}_1}-u_2.---\left(8\right)$

(5)空气-假想水膜边界层的速度分布积分S与两相流单流体模型的边界层速度分布积分S_m进行比较，求二者的差值。S和S_m的表达式分别为，

$S=\frac{1}{2}{u}_2{h}_f+u_1(H-h_f),---\left(9\right)$

$S_m=\frac{1}{2}{u}_{m}H.---\left(10\right)$

(6)比较步骤(5)的结果S和S_m，如果在误差范围内则结束，否则调整运动粘度，回到步骤(2)，直到求得精确的假想水膜表面速度u_2f。

在获得假想水膜表面速度u_2f后，同样按照不可压缩流边界层理论求假想水膜沿着壁面的法向速度，而该速度也是超冷水滴在壁面法向速度分量v₂，

$v_2=\frac{0.8604\·v}{U_{\∞}\·x}{u}_{2f}.---\left(11\right)$

至此可以将假想水膜厚度调整至真实的水膜厚度，即

h_f＝v₂·Δt，(12)

式中，Δt代表积分时间长度。真实的水膜表面流动速度

也可由边界层内部流动速度的线性分布的假设获得，即

与u₂之比等于h_f与h₂之比。

下面的计算流程将进入结冰状态模型。本发明提出的水膜的结冰状态模型是一个计算冰层形状和内部的温度分布的模型。其中，需要使用网格加密方法追踪结冰界面。

因为飞行结冰是由于超冷水滴积聚在飞行器表面上形成的水膜的内部的相变产生的，所以，将整体计算区域按照水膜的边界位置重新划分为外部流场区和内部结冰计算区。结冰过程的计算区域应该包括水膜和与水膜接触的飞行器壁面或是以前已经生成的冰层。水膜和冰层的接触面，也称作结冰界面。本发明提出的模型需要在水膜中生成至少三层计算网格，在其下方的冰层中至少生成三层网格。这六层网格形成相变区，是在原计算网格中加密的计算网格，水凝结成冰的相变过程将发生在其中，并形成新的结冰界面。图3给出了结冰模型使用的网格划分方法和加密方法的原理图。其中，图3(a)中表示在t时刻，冰层的上部是超冷液态水滴形成了水膜，在水膜中生成三层计算网格，即在原计算网格内加密，可见重新划分网格后，内部区覆盖水膜和已经结成的冰层。图3(b)中表示在t+Δt时刻，水膜中的下面两层已经结成了冰，同时，在冰层上方有新的水膜形成，同样在其中生成三层加密网格，与原来的网格共同形成了内部结冰状态计算的计算域。如果水膜是在未结冰的边界上，则直接使用传统的Messinger结冰模型，不对水膜进行网格加密，直接求取水膜的结冰量，即飞行器壁面形成的冰层。Messinger结冰模型是零维模型，不计算水膜和冰层内部的温度分布，该模型的形式和求解方法是公知的。

水膜内部划分网格后开始进行结冰过程的计算。水膜的结冰过程可以视为带有相变的不可压缩液-固两相流流动问题。其控制方程包括连续方程、动量方程和以温度变化表示的能量方程，求解后给出水膜和冰层内部的温度分布和结冰量。外部速度、温度已由流场的解获得，成为已知的边界条件，其求解数值方法为公知的。图4给出了该模型的输入、输出变量关系。输入参数是水膜表面速度、温度，输出变量是冰层的高度既是结冰界面的位置和冰层内部温度分布。

如果该模型得出结冰界面等于水膜高度，说明水膜完全结成冰，是一种霜状冰；如果结冰界面的位置小于水膜的高度，说明水膜并未完全凝结，形成瘤状冰；此时，需要将剩余的水膜转变成网格内的水滴的体积分数，作为下一个时间间隔的外部流场计算的壁面边界条件。如果剩余水膜厚度为h′_m、上标’表示剩余水膜折算后的体积分数，则从图2(a)可知，

$\frac{{2h}_1-h_m^{\′}}{h_m^{\′}}=\frac{{\mathrm{\α}}_1^{\′}}{{\mathrm{\α}}_2^{\′}},---\left(13\right)$

显然仍有α′₁+α′₂＝1，同样满足公式(1)，。

完成一个时间间隔内的飞行结冰计算后，按照有新的结冰界面重新划分为外部流场区和内部结冰区，进入下一个时间间隔的计算。图5表示了本发明提出的飞行器飞行结冰的数值模拟方法计算流程的详细步骤是，

(1)未结冰的飞行器周围计算网格的生成；

(2)规定计算的开始时刻；

(3)进行飞行器外部空气-超冷水滴的运动的单流体两相流流动模拟；

(4)求飞行器壁面网格内假想水膜的厚度；

(5)进行空气-水膜的两相流速度分解，获得假想水膜表面速度；

(6)求假想水膜表面运动的法向速度；

(7)求正是水膜的厚度和表面速度；

(8)检测壁面是否已经有结冰。如果没有结冰，用Messinger结冰模型计算结冰量，否则进入下个步骤；

(9)将水膜和其下方的冰层进行网格加密，构成相变区，进入结冰状态模型；

(10)在冰层内部计算温度分布，同时求得飞行器壁面结冰量并构成新的结冰界面；

(11)将剩余的水膜折算到外部计算域的边界中的超冷水滴的体积分数中；

(12)重新划分外部流场计算域和内部结冰计算域；

(13)回到步骤(2)，进行下一个时刻的流体两相流单流体流动模拟计算。

本发明提出的数值方法的流程的特点是：

1.计算开始生成的计算网格将作为整个计算的背景网格，在以后时刻的计算中，该网格不做移动，只是在局部进行加密处理，是非结构化网格；

2.飞行器外部流场按照单流体两相流处理，其输出结果是混合流体的速度、压力、密度、温度等参数，在壁面网格内进行空气-水膜速度的两相流分解，在水膜-冰层内部计算结冰的相变过程。

3.在每个计算时刻内的整体计算域分为外部计算域和内部计算域，根据上个时间段内构成的结冰边界线判定。外计算域的网格用于空气-超冷水滴的运动的单流体两相流流动方程组的求解；内计算域的网格用于水膜-冰层内温度分布和相变的偏微分方程组的求解。求解后获得新的结冰界面，即结冰的位置和形状。

本发明的计算方法不需要在每个时间间隔内因为结冰造成的飞行器外部形状的改变而重新构造计算网格、无需进行由于网格变动而进行的任何插值运算，保证了计算精度，节约了计算时间。此外，结冰量和冰层内部的温度分布由结冰模型计算结果同时给出，即求得壁面结冰量和冰层几何形状的同时，获得了冰层内的温度分布，是对飞行结冰过程的认识、分析、以及除冰系统设计的重要信息。

附图说明

图1Fensap软件的计算流程图

图2水膜表面速度分解的原理图

(a)壁面网格内的空气-超冷水滴的分离和假想水膜形成原理图

图中，1：飞行器壁面、2：超冷液态水膜的几何中心h₂、3：空气的几何中心h₁；

(b)假想水膜表面速度的分解使用的不可压缩流边界层示意图

(c)边界网格内的两相流的单流体混合边界层的速度分布

图中，1：水膜速度u₂、2：水膜表面h_f处速度u_2f(点D)、3：水膜表面空气速度点(点C)、4：在h₁处空气速度u₁、5：壁面网格边界的空气速度(点B)、6：点A、7：水膜速度为零处(点O)；

(d)边界网格内的水膜-空气边界层的速度分布

图中，1：混合流体在h_m处的速度u_m、2：混合流体在H处的速度(点B)、3：点A、4：混合流体速度为零处(点O)；

(e)水膜表面速度的分解算法的流程图

图3结冰状态模型使用的网格划分方法和加密方法的原理图

(a)t时刻水膜-冰层和计算网格

图中，1：飞行器壁面、2：时刻t以前生成的冰层、3：t时刻三层水膜中的第一层、4：t时刻三层水膜中的第二层、5：t时刻三层水膜中的第三层；

(b)t+Δt时刻的水膜-冰层和计算网格

图中，1：飞行器壁面、2：时刻t以前生成的冰层、3：时刻t的第三层水膜生成的冰层、4：时刻t的第三层水膜生成的冰层、5：t+Δt时刻三层水膜中的第一层、4：t+Δt时刻三层水膜中的第二层、5：t+Δt时刻三层水膜中的第三层：

图4结冰状态模型的输入、输出变量关系图

图5本发明提出的飞行器飞行结冰的数值模拟方法计算流程图

图6围绕二维NACA0012机翼的初始计算网格

图7计算结束时的结冰状态

具体实施方式

以下以一个具体实施方式进一步说明本发明，一种用于飞行器飞行结冰的数值模拟方法，这种数值模拟技术可以用Fortran90计算机高级程序语言实现，并通过计算机运行来模拟二维NACA0012机翼在空中飞行期间遭遇结冰时的状态。

数值模拟的条件是：

来流速度：57m/s

来流温度：243K

大气压力：100kPa

攻角：0°

液态水含量：2.58g/m³

超冷水滴直径：50μm

冰的密度：917kg/m³

根据图5给出的本发明提出的飞行器飞行结冰的数值模拟方法计算流程的详细步骤是，

1.未结冰的飞行器周围计算网格的生成。首先在NACA0012周围生成计算网格。这一个C-型二维结构化网格，如图6所示，围绕机翼256个网格，沿壁面垂直方向96个网格，并采用正交处理。

2.规定计算的开始时刻t＝0，此时机翼表面尚未结冰。同时，规定时间步长Δt。

3.进行飞行器外部空气-超冷水滴的运动的单流体两相流流动模拟。首先给出模拟空气-超冷水滴运动的单流体两相流模型。以下变量的下标1表示空气，下标2表示超冷水滴，下标m表示混合物，下标t表示湍流。密度、速度、压力温度、时间，分别用ρ，

p，T，t表示。速度矢量在笛卡尔坐标x，y方向上的分量分别为u，v。物性参数粘性系数、定压比热分别用μ，Cp表示。二元混合物中，空气的体积分数为α₁、超冷水滴的体积分数为α₂。在单流体两相流模型中，上述混合物参数以及密度和速度均可以由从两种组分中的对应参数按照体积分数进行加权平均获得。空气-超冷水滴运动的二维单流体两相流模型是：组分扩散方程

$\frac{{\∂\mathrm{\α}}_2}{\∂t}+\▿\·\left({\mathrm{\α}}_2{\stackrel{\→}{V}}_m\right)=0;---\left(14\right)$

连续方程

$\frac{{\∂\mathrm{\ρ}}_m}{\∂t}+\▿\·\left({\mathrm{\ρ}}_m{\stackrel{\→}{V}}_m\right)=0;---\left(15\right)$

动量方程

$\frac{\∂\left({\mathrm{\ρ}}_m{\stackrel{\→}{V}}_m\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_m{\stackrel{\→}{V}}_m{\stackrel{\→}{V}}_m\right)=-\▿p_m+{\mathrm{\λ}}_m\▿(\▿\·\stackrel{\→}{V})+({\mathrm{\μ}}_m+{\mathrm{\μ}}_T)({\▿}^2{\stackrel{\→}{V}}_m+\▿(\▿\·{\stackrel{\→}{V}}_m))+{\mathrm{\ρ}}_m\stackrel{\→}{g};---\left(16\right)$

能量方程

$\frac{\∂\left({\mathrm{\ρ}}_m{E}_m\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_m{\stackrel{\→}{V}}_m{E}_m\right)=---\left(17\right)$

$-\▿\·\left(p_m{\stackrel{\→}{V}}_m\right)+[{\mathrm{\λ}}_{m}I(\▿\·\stackrel{\→}{V})+({\mathrm{\μ}}_m+{\mathrm{\μ}}_T)({\▿}^2{\stackrel{\→}{V}}_m+{\stackrel{\→}{V}}_m\▿)]\·\▿{\stackrel{\→}{V}}_m-\▿\·{\stackrel{\→}{Q}}_m+{\mathrm{\ρ}}_m\stackrel{\→}{g}\·{\stackrel{\→}{V}}_m$

其中，E_m是混合流体内能、Q_m是热流量、

是重力、λ_m由Stokes定律给出。

上述方程组的求解采用有限体积法，变量存储在计算网格中心，方程组的空间离散采用二阶Roe方法，时间离散采用LU-SGS方法，均是公知的这里不再叙述。

4.在当地坐标下进行空气-水膜的两相流速度分解，获得水膜表面速度和厚度。根据图2关于假想水膜表面速度的分解算法流程图。其步骤是：

(1)设定一个壁面上超冷液态水的运动粘度v为初始值。

(2)按照边界层速度分布公式(5)求得假想水膜表面的速度u_2f。公式(5)中假想水膜厚度h_f代替y来流速度是无穷远处混合流体的速度，即飞行器的飞行速度。

(3)按照质量平均求液态水边界层的速度u₂。

(4)求空气的速度u₁。

(5)空气-水膜边界层的速度分布积分S与两相流单流体模型的边界层速度分布积分S_m进行比较，求二者的差值。

(6)比较步骤(5)的结果S和S_m，如果在误差范围内则结束，否则调整运动粘度，回

到步骤(2)，直到求得精确的假想水膜表面速度u_2f。

5.求假想水膜表面运动的法向速度。按照公式(11)求假想水膜沿着壁面的法向速度，而该速度也是超冷水滴在壁面法向速度分量v₂。

6.求真实水膜的厚度和表面速度。按照公式(12)将假想水膜厚度调整至真实的水膜厚度。

7.用线性插值方法获得真实的水膜表面流动速度

8.因为在t＝0时刻没有已经形成的冰层。可以用Messinger结冰模型计算结冰量，其过程是公知的，不在叙述。在获得结冰形状后，将冰层占据的网格标示出来，以备下个时刻使用。在t+Δt，如果水膜覆盖在被标示的网格上，则进入下个步骤。

9.将水膜和其下方的冰层进行网格加密，构成相变区。

10.在冰层内部计算温度分布，同时求得飞行器壁面结冰量并构成新的结冰界面。

结冰状态模型的方程组是在当地坐标系下建立。当地坐标系由ξ，ζ两正交方向组成，其中方向ζ-是壁面外法线方向。结冰状态模型将水膜和冰层视为一个整体计算域，在其中求解不可压缩流的流动，其中用液相率表示发生在冰层和水膜的交界面处的结冰，用特大阻力表示冰层中水膜的运动特征。水膜在ξ，ζ两正交方向流动速度分别用u，v表示、压力用p表示，此外，密度ρ、定压比热C_p、导热系数k、结冰潜热L表示。下表w代表水膜、i代表冰层。具体地，

连续方程

$\frac{\∂u}{\∂\mathrm{\ξ}}+\frac{\∂v}{\∂\mathrm{\ξ}}=0;---\left(18\right)$

能量方程

$\frac{\∂\left({\mathrm{\ρ}}_{*}{C}_{p^{*}}T\right)}{\∂t}+\frac{\∂\left({\mathrm{uC}}_{p^{*}}T\right)}{\∂x}+\frac{\∂\left({\mathrm{vC}}_{p^{*}}T\right)}{\∂y}=\frac{\∂}{\∂x}\left(k_{*}\frac{\∂T}{\∂x}\right)+\frac{\∂}{\∂x}\left(k_{*}\frac{\∂T}{\∂y}\right)-L[\frac{\∂f}{\∂t}+\frac{\∂\left(\mathrm{uf}\right)}{\∂x}+\frac{\∂\left(\mathrm{vf}\right)}{\∂y}];---\left(19\right)$

其中，液相率f定义为

$f=\frac{{\mathrm{\ρ}}_{*}{C}_{p^{*}}T-{\mathrm{\ρ}}_i{C}_{p^i}T}{{\mathrm{\ρ}}_w{C}_{\mathrm{pw}}T-{\mathrm{\ρ}}_i{C}_{\mathrm{pi}}T};---\left(20\right)$

式中，

ρ_*＝fρ_w+(1-f)ρ_i； (21)

k_*＝fk_w+(1-f)k_i；    (22)

$C_{p^{*}}=f{C}_{\mathrm{pw}}+(1-f)C_{\mathrm{pi}}+L\frac{\∂f}{\∂T};---\left(23\right)$

动量方程

$\frac{\∂u}{\∂t}+\frac{\∂\left(\mathrm{uu}\right)}{\∂\mathrm{\ξ}}+\frac{\∂\left(\mathrm{uv}\right)}{\∂\mathrm{\ξ}}=-\frac{1}{{\mathrm{\ρ}}_w}\frac{\∂p}{\∂\mathrm{\ξ}}+[\frac{\∂}{\∂\mathrm{\ξ}}\left(v_{*}\frac{\∂u}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(v_{*}\frac{\∂u}{\∂\mathrm{\ζ}}\right)]-g_{\mathrm{\ξ}}-D;---\left(24\right)$

$\frac{\∂v}{\∂t}+\frac{\∂\left(\mathrm{vu}\right)}{\∂\mathrm{\ξ}}+\frac{\∂\left(\mathrm{vv}\right)}{\∂\mathrm{\ζ}}=-\frac{1}{{\mathrm{\ρ}}_w}\frac{\∂p}{\∂\mathrm{\ζ}}+[\frac{\∂}{\∂\mathrm{\ξ}}\left(v_{*}\frac{\∂v}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(v_{*}\frac{\∂v}{\∂\mathrm{\ζ}}\right)]-g_{\mathrm{\ζ}}-D;---\left(25\right)$

其中，特大阻力

$D=\left\{\begin{array}{cc}0,& \mathrm{\&epsiv;}<f\&le;1\\ A,& f\&le;\mathrm{\&epsiv;}\end{array}\right.,---\left(26\right)$

式中，A是一个大数，在10⁷量级；ε是人工参数，取0.001。

不可压缩流的求解以一类压力修正法最为普遍，例如SIMPLE算法就是其中一种。本案例使用同位网格法离散方程，同时使用网格边界上的动量插值以避免速度压力失偶。具体的求解过程是公知的，其结果，即水膜内部速度分布，带入能量方程，可获得整体内部区域的温度分布。结冰界面依靠控制方程中的液相率来判定。

11.按照公式(13)将剩余的水膜折算到外部计算域的边界中的超冷水滴的体积分数中；

12.重新划分外部流场计算域和内部结冰计算域；

13.回到步骤2，进行下一个时刻的模拟计算。

图7给出计算结束时的结冰状态，包括冰层的形状及冰层内部温度分布的等温线。来流温度和壁面温度都是243K，不考虑机翼壁面的热流。结果表明结冰沿着机翼前缘厚度约为0.02倍机翼弦长，冰层内部的温度等温线表明温度分布均匀，距离壁面的温度略低于壁面温度和外界温度。

Claims (4)

1.一种用于模拟飞行器的飞行结冰的数值方法，来模拟飞行器在空中飞行期间遭遇结冰时的状态。其主要特征是包括用于模拟空气-超冷水滴的运动的两相流流动的单流体模型中壁面水膜表面的速度分解和水膜厚度的算法、在计算冰层形状和内部的温度分布的水膜结冰状态模型中采用网格加密方法追踪结冰界面的算法、基于固定计算网格利用上述模型和算法进行飞行结冰数值模拟计算的流程。

2.根据权利要求1所述的一种用于模拟飞行器的飞行结冰的数值方法，其特征在于，所述的用于模拟空气-超冷水滴的运动的两相流流动的单流体模型中壁面水膜表面的速度分解和水膜厚度的算法，具体步骤是：

(1)设定一个壁面上超冷液态水的运动粘度v为初始值；

(2)按照不可压缩边界层速度分布公式求得假想水膜表面的速度u_2f；

(3)按照质量平均求液态水边界层的速度u₂；

(4)求空气的速度u₁；

(5)空气-水膜边界层的速度分布积分S与两相流单流体模型的边界层速度分布积分S_m进行比较，求二者的差值；

(6)比较步骤(5)的结果S和S_m，如果在误差范围内则结束，否则调整运动粘度，回到步骤(2)，直到求得精确的假想水膜表面速度；

(7)按照不可压缩流边界层理论求假想水膜沿着壁面的法向速度，而该速度也是超冷水滴在壁面法向速度分量v₂；

(8)用积分时间长度Δt乘以超冷水滴在壁面法向速度分量v₂，获得真实的水膜厚度；

(9)真实的水膜表面流动速度

由边界层内部流动速度的线性分布的假设获得。

3.根据权利要求1所述的一种用于模拟飞行器的飞行结冰的数值方法，其特征在于，所述的在计算冰层形状和内部的温度分布的水膜结冰状态模型中采用网格加密方法追踪结冰界面的算法，需要在水膜中生成至少三层计算网格，在其下方的冰层中至少生成三层网格，这六层网格形成相变区，水凝结成冰的相变过程将发生在其中，并形成新的结冰界面。

4.根据权利要求1所述的一种用于模拟飞行器的飞行结冰的数值方法，其特征在于，所述的基于固定计算网格利用上述模型和算法进行飞行结冰数值模拟计算的流程，具体步骤是：

(1)未结冰的飞行器周围计算网格的生成；

(2)规定计算的开始时刻；

(3)进行飞行器外部空气-超冷水滴的运动的单流体两相流流动模拟；

(4)求飞行器壁面网格内假想水膜的厚度；

(5)进行空气-水膜的两相流速度分解，获得假想水膜表面速度；

(6)求假想水膜表面运动的法向速度；

(7)求正是水膜的厚度和表面速度；

(8)检测壁面是否已经有结冰。如果没有结冰，用Messinger结冰模型计算结冰量，否则进入下个步骤；

(9)将水膜和其下方的冰层进行网格加密，构成相变区，进入结冰状态模型；

(10)在冰层内部计算温度分布，同时求得飞行器壁面结冰量并构成新的结冰界面；

(11)将剩余的水膜折算到外部计算域的边界中的超冷水滴的体积分数中；

(12)重新划分外部流场计算域和内部结冰计算域；

(13)回到步骤(2)，进行下一个时刻的流体两相流单流体流动模拟计算。

Images