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

说 背景技术

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

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

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

飞行器飞行结冰的数值模拟技术中所采用的流 程都是基于三个主模块的调用。三个主模 块及其各自主要功能是分别为- 空气流动模块：用于求解飞行器外部流场（包 括飞行器表面）信息，即求解流体流动（此 处为空气运动） 控制方程；

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

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

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

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

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

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

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

(^ + «2 = 1。 ( 1 ) 混合物的密度 p _m 和粘性系数 / _m 按照体积分数进行加权平均

/ _m = + « ₂ / ₂ 。 ( 3 ) 速度 ii _m 的体积分数进行加权平均分别是 ti _m ( 4 )

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

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

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

其中， ν是流体的运动粘度系数。 所以， 来流速度 t/ _∞ 和运动粘度 V确定后， 公式（5)和（6) 可以确定边界层中某点 (χ, )的速度 Μ的大小。 一种情况是，边界网格中的假想水膜被包含在 边界层内，假想水膜和空气共同形成的边 界层的速度分布。在假想水膜和空气接触面上 ， 由于两相物性差异， 存在滑移速度。 为了求 得假想水膜表面速度。 该值需要从两相流单流体模型中的混合物速度 中分解出来。 图 2(c) 给出按照两相流单流体模型的边界层速度分布 。 其中 H是网格高度， 是几何中心， 1^是 几何中心处的两相混合速度。 图 2(d)是假想水膜 -空气边界层的速度分布。 其中， _Ml 是空气 的速度、 t/ ₂ 是假想水膜的速度、 t/ ₂ 是假想水膜表面的速度,其他符号和图 2(a)—致。此外， 需要建立假想水膜-空气边界层的速度分布积� �与两相流单流体的边界层速度分布积分相等 的假设， 即图 2(c)中的 ABC的面积等于图 2(d)中的 ABCD的面积。 目的是为了将假想水膜表 面的速度 _M , _f 从混合流动中分解出来。 图 2(e)给出了水膜表面速度的分解算法流程图， 其具 体步骤是:

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

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

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

1

U -—U (7)

/

2

求空气的速度 t/i。 因为混合流体的速度由公式（4)表示， 所以单相空气速度为， p _m u _m

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

S = ~^u ₂ h _f +Η _Χ {Η -h ), (9) S _m =l _Um H。 (10)

2

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

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

0.8604 - _{/ 1 1 λ}

V ₂ =— "2/ °

U _∞ ,x 至此可以将假想水膜厚度调整至真实的水膜厚 度， 即

h _f ^v ₂ -At, (12) 式中， At代表积分时间长度。 真实的水膜表面流动速度 ϋ ₂ 也可由边界层内部流动速度的 线性分布的假设获得， 即 ϋ ₂ 与 _Μ2 之比等于 与/ ¾之比。

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

因为飞行结冰是由于超冷水滴积聚在飞行器表 面上形成的水膜的内部的相变产生的，所 以，将整体计算区域按照水膜的边界位置重新 划分为外部流场区和内部结冰计算区。结冰过 程的计算区域应该包括水膜和与水膜接触的飞 行器壁面或是以前已经生成的冰层。水膜和冰 层的接触面， 也称作结冰界面。本发明提出的模型需要在水 膜中生成至少三层计算网格， 在 其下方的冰层中至少生成三层网格。这六层网 格形成相变区，是在原计算网格中加密的计算 网格， 水凝结成冰的相变过程将发生在其中， 并形成新的结冰界面。 图 3给出了结冰模型使 用的网格划分方法和加密方法的原理图。 其中， 图 3 (a) 中表示在 ί时刻,冰层的上部是超 冷液态水滴形成了水膜， 在水膜中生成三层计算网格， 即在原计算网格内加密， 可见重新划 分网格后， 内部区覆盖水膜和已经结成的冰层。 图 3 (b) 中表示在 ί^ί时刻,水膜中的下 面两层已经结成了冰， 同时， 在冰层上方有新的水膜形成， 同样在其中生成三层加密网格， 与原来的网格共同形成了内部结冰状态计算的 计算域。如果水膜是在未结冰的边界上，则直 接使用传统的 Messinger结冰模型， 不对水膜进行网格加密， 直接求取水膜的结冰量， 即飞 行器壁面形成的冰层。 Messinger结冰模型是零维模型，不计算水膜和冰 层内部的温度分布， 该模型的形式和求解方法是公知的。 水膜内部划分网格后开始进行结冰过程的计算 。水膜的结冰过程可以视为带有相变的不 可压缩液-固两相流流动问题。 其控制方程包括连续方程、 动量方程和以温度变化表示的能 量方程， 求解后给出水膜和冰层内部的温度分布和结冰 量。外部速度、温度已由流场的解获 得， 成为已知的边界条件， 其求解数值方法为公知的。 图 4给出了该模型的输入、输出变量 关系。输入参数是水膜表面速度、温度， 输出变量是冰层的高度既是结冰界面的位置和 冰层 内部温度分布。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(13) 回到步骤 （2) ， 进行下一个时刻的流体两相流单流体流动模拟 计算。 本发明提出的数值方法的流程的特点是：

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

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

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

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

图 1 Fensap软件的计算流程图

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

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

图中， 1 : 飞行器壁面、 2 : 超冷液态水膜的几何中心/ ¾、 3 :空气的几何中 心 / _{1 ;}

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

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

图中， 1 : 水膜速度^、 2： 水膜表面/ ^处速度 M ₂ (点 D)、 3 :水膜表面空 气速度点（点 C)、 4： 在/ ¾处空气速度 _Ml 、 5： 壁面网格边界的空气速度（点

B)、 6： 点八、 7： 水膜速度为零处 （点 0);

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

图中， 1： 混合流体在 处的速度^、 2： 混合流体在 H处的速度 （点 B )、 3 :点八、 4： 混合流体速度为零处 （点 0); (e) 水膜表面速度的分解算法的流程图

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

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

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

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

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

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

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

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

图 7 计算结束时的结冰状态 具体实施方式

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

数值模拟的条件是：

来流速度： 57m/s

来流温度： 243K

大气压力： lOOkPa

攻角： 0°

液态水含量： 2. 58g/m ³

超冷水滴直径： 50μιη

冰的密度： 917kg/m ³

根据图 5给出的本发明提出的飞行器飞行结冰的数值� �拟方法计算流程的详细步骤是， 1. 未结冰的飞行器周围计算网格的生成。 首先在 NACA0012周围生成计算网格。 这一个 C- 型二维结构化网格， 如图 6所示， 围绕机翼 256个网格， 沿壁面垂直方向 96 个网格， 并采用正交处理。 规定计算的开始时刻 ίθ,此时机翼表面尚未结冰。 同时， 规定时间步长 Ζ\ί。

进行飞行器外部空气-超冷水滴的运动的单流� �两相流流动模拟。首先给出模拟空气-超 冷水滴运动的单流体两相流模型。 以下变量的下标 7表示空气， 下标 2表示超冷水滴， 下标 表示混合物， 下标 ί表示湍流。密度、速度、压力温度、 时间， 分别用 p, ,p,7 表示。 速度矢量在笛卡尔坐标 方向上的分量分别为 M,V。 物性参数粘性系数、 定压 比热分别用/ Αθ^表示。 二元混合物中， 空气的体积分数为《,、 超冷水滴的体积分数为

« ₂ 。在单流体两相流模型中，上述混合物参 数以及密度和速度均可以由从两种组分中的 对应参数按照体积分数进行加权平均获得。 空气-超冷水滴运动的二维单流体两相流模 组分扩散方禾 '王;

连续方禾 '主 + V · l ? _m V, 0； (15) 动量方程 d[p _m y _m

+ V · [p _m y _m y _m = -Vp _m + _m V[V ·Υ +{μ _Μ +μ _τ )[v ^z y _m +VlV«V _m + p ,； (16) dt 能量方:

- V · {p _m Y _m )+ [ {V•Υ)+{μ _ιη +μ _τ )\^ ^Ζ Υ _Μ + Y _m V _} ,VV _m -V«Q _m + ? _m -V _ffl 其中， E _m 是混合流体内能、 g _m 是热流量、 是重力、 由 Stokes定律给出。 上述方程组的求解采用有限体积法，变量存储 在计算网格中心，方程组的空间离散采用 ：阶 Roe方法， 时间离散采用 LU-SGS方法， 均是公知的这里不再叙述。 在当地坐标下进行空气 -水膜的两相流速度分解， 获得水膜表面速度和厚度。 根据图 2 关于假想水膜表面速度的分解算法流程图。 其步骤是：

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

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

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

( 4) 求空气的速度 _Ml 。

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

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

5. 求假想水膜表面运动的法向速度。 按照公式 （11 ) 求假想水膜沿着壁面的法向速度， 而 该速度也是超冷水滴在壁面法向速度分量 v ₂ 。 6. 求真实水膜的厚度和表面速度。按照公式（12) 将假想水膜厚度调整至真实的水膜厚度。

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

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

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

10. 在冰层内部计算温度分布， 同时求得飞行器壁面结冰量并构成新的结冰界 面。 结冰状态模型的方程组是在当地坐标系下建立 。 当地坐标系由 ^, ζ"两正交方向组 成，其中方向 是壁面外法线方向。结冰状态模型将水膜和冰 层视为一个整体计算域， 在其中求解不可压缩流的流动， 其中用液相率表示发生在冰层和水膜的交界面 处的结 冰，用特大阻力表示冰层中水膜的运动特征。 水膜在 ^, ζ"两正交方向流动速度分别用 表示、 压力用 ρ表示， 此外， 密度 ρ、 定压比热 C _p 、 导热系数 、 结冰潜热 表示。 下表 w代表水膜、 代表冰层。 具体地，

连续方程 + _{0 ;} ( 18) θξ dg 能量方程 其中， 液相率/定义为

式中，

(21) k, = fk _w + (\- (22)

动量方程

其中， 特大阻力

0, ^ < /≤ 1

D = (26) 式中， ^是一个大数， 在 10 ⁷ 量级； s是人工参数， 取 0.001。

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

11. 按照公式 （13) 将剩余的水膜折算到外部计算域的边界中的超 冷水滴的体积分数中； 12. 重新划分外部流场计算域和内部结冰计算域；

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

图 7给出计算结束时的结冰状态，包括冰层的形� �及冰层内部温度分布的等温线。来流 温度和壁面温度都是 243K， 不考虑机翼壁面的热流。 结果表明结冰沿着机翼前缘厚度约为

0.02 倍机翼弦长， 冰层内部的温度等温线表明温度分布均勾， 距离壁面的温度略低于壁面 温度和外界温度。