低温模型表面水蒸气/CO2凝华数值模拟研究★

2022-08-24彭福胜郭佳全黄春杰

山西建筑 2022年17期

王杰，彭福胜，郭佳全，黄春杰

(陆军工程大学，江苏南京 210000)

低温风洞[1]以氮气为运行介质，通常工作在110 K左右的低温工况下。当模型进行气动特性实验时，若风洞内空气露点、CO2含量不满足要求，则水蒸气和CO2会在模型表面结露、结霜，进而影响机翼表面粗糙度，这将大幅影响试验数据精准度。目前世界上仅有三座生产型低温风洞NTF[2]，KKK和ETW[3]能够模拟大型飞机模型试验所需的雷诺数，均采用相应措施对空气露点和CO2含量进行严格控制。然而，常温常压下水蒸气含量随露点升高呈非线性升高趋势，-60 ℃附近的露点稍有偏差，水蒸气含量可能呈量级增加，这对低温风洞参数控制提出更高要求。因此，获得干燥空气中水蒸气/CO2的凝华特性及随时间的变化规律，可以大幅降低风洞参数设计风险，确保风洞性能。

结霜[4]是生活中常见的一种现象，广泛存在于制冷空调、低温技术及其他相关领域。目前已有大量文献和著作对水蒸气凝华现象进行研究，云雾室研究[5]、航空航天的冷凝结霜[6]、天然气液化及分离装置[7]等，这些都是自然界和工程中发生的冷凝结霜现象。孟繁炯等[8]提出一种将霜晶看作由毛细管组成的多孔介质物理模型，建立了一维霜层弯曲度函数关系式方程。基于平板气流结霜物理模型[9]，可以测算出霜粗糙度和求解空气当量导热系数与霜内层导热系数的组合方程。CUI等[10]提出一种新模型来预测霜的形成和生长，借助成核理论研究潮湿空气中流动冷平板的霜特性，得到霜结构及其性质的空间和时间变化规律。QI等[11]研究氮气流动中低温表面上的痕量水蒸气结霜机制，讨论了氮气密度对霜层性质的影响。夏斌等[12]进行了快速来流条件下一维干模态结霜模拟分析，掌握单一因素对结霜的影响规律。CHOI[13]提出了一种初始壁冷却的霜形成晶核的数值模型，较好地解释低温表面的传热传质过程以及填充和滞留过程。YOONG提出了霜冻条件下ASHP系统的仿真模型，找出了在结霜过程中引起热泵系统性能变化的因素，指出性能的下降与霜的沉积没有直接关联。

考虑到目前国内低露点环境凝华研究较少，缺乏相关凝华数据，霜层测量难度大等问题，为解决上述问题，因此开展了干燥空气中的水蒸气/CO2凝华数值模拟研究。基于分散的多相模型，分析并得到霜层的生长过程和流场对结霜厚度的影响[14]。

1 物理模型

为更好复现低温风洞真实模型表面结霜情况，选取二元翼型作为冷却模型，模型尺寸为309 mm×300 mm×35.6 mm，如图1所示。

设定物理模型时，为了简化复杂的计算过程，考虑翼型整体温度分布均匀，将二元翼型表面设置为固定的工况温度。模拟的环境空间尺寸1 200 mm×1 000 mm×900 mm，如图2所示。考虑翼型延展方向(Y向)上的一致性，翼型处处相同且没有扭转，仅在二元翼型宽度方向上保留25.4 mm的厚度，形成如图3所示的模型。

选取六面体网格进行划分，在小面及细节处进行网格加密处理，为了准确捕捉翼型的几何特征，保留了翼型的特征线，在翼型附近进行了网格加密；进行体网格划分，方便网格变形表征翼型表面的结霜效果，对翼型表面划分了三层边界层网格，划分的体网格局部放大图如图4所示。

前期对网格进行无关性验证，当网格单元数量超过14万时，网格数量对结果的影响基本保持稳定，并在模拟几何中使用148 003个元素获得网格独立性。

2 数值方法

结霜过程并非仅仅是简单的相变过程，是一种复杂多变的过程，将气相转换为液滴，基于分散的多相模型模拟机翼结霜。机理假设如下：

1)空气流速低且压力变化小，空气被视为不可压缩理想气体。

2)水蒸气和二氧化碳使用微小液滴代替，且体积分数小，对气流无影响。

3)不考虑液滴运动过程中的热传递或蒸发，且物理参数不改变，不受湍流的影响。

2.1 控制方程

1)通过控制体积的质量平衡的连续性方程：

(1)

其中，ρ为密度；v为连续介质的速度。

2)线性动量的时间变化率等于作用在连续体上的合力：

(2)

其中，⊗为外积;fb为作用在连续体上的每单位体积力(如重力和离心力)的合力；σ为应力张量。

(3)

其中，τ为表面应力分量，其下标的第一个字母代表应力作用表面，第二个字母代表应力的作用方向。对于流体，应力张量通常被写成法向应力和剪切应力之和：

σ=-pI+T

(4)

其中，p为压力；T为黏性应力张量，给出：

(5)

3)能量守恒方程可以写成：

(6)

其中，E为单位质量的总能量；q为热流；SE为单位体积的能源。

2.2 凝固模型

凝固模型可以使用焓公式来确定液固相凝固部分的分布。

(7)

(8)

其中，T*为标准化温度，定义为：

(9)

f(T*)=1-T*

(10)

3 计算设置

初始条件设置流场速度沿X轴正方向，大小为0.5 m/s，设置重力方向沿Z轴负方向，数值为-9.81 m/s2。模型所有壁面都为滑移壁面。Z轴负方向壁面特别设置允许离散相逃逸，防止假设的液滴在重力的作用下在Z轴负方向壁面堆积，从而结霜。在换热方面除了机翼本身外，周边环境、模拟域进出口的温度均为常温300 K。物质物性参数如表1所示。

表1 物性参数表

计算时间设置为120 min，进行瞬态求解分析，由于流速较低，其流态很快趋于稳定，因此设置较大的时间步长，在保证计算不发散的情况下设置时间步长为10 s，每个时间步内部迭代10次，即总计算步数为7 200步。

4 结果分析

4.1 霜层生长分析

为了解低温模型表面霜层生长过程，通过比较不同工质的霜层生长，在不同的离散相浓度下，其霜层生长变化过程随时间的变化是相似的，因此以冷表面温度110 K，空气流速0.5 m/s，水蒸气含量(质量浓度)10×10-6为例，观察霜层生长过程，不同时刻的霜层厚度云图如图5所示。

从图5可以看出，随着时间的增加，霜层变化较明显，霜层高度不断增加，霜逐步累积，霜层变得越来越紧密。模型的前端和后端的结霜量较小，霜层主要集中在中部区域，这主要受流场分布的影响。从总体来看，模拟后的机翼表面霜厚度增加变化明显，霜层生长趋势较为一致。图6为10×10-6含量下霜层厚度随时间变化的模拟结果。由图6可得，以20 min为时间基准，随着时间的不断推移，霜层厚度增加变化明显，霜厚度分别为3.6 μm,7.6 μm,11 μm,14 μm,16.1 μm,18.3 μm，但霜的增长率较小，每20 min仅增加2 μm～4 μm，时间对霜厚度贡献率较低。

图7为不同时刻的速度云图。由图7可知，在初始阶段,翼型上表面流动分离不明显，霜层生长主要受重力影响，过冷液滴与模型前部撞击较多，模型前端捕捉更多液滴形成霜层，60 min后气流在翼型尾部逐步形成周期性扰动，更多的小液滴逐步在机翼中部及后部撞击，霜层迅速在中部发展，霜层不断加厚。另一方面，受重力作用以及气体流动拖曳力共同作用的影响，气流未充分发展前，前期离散相高浓度分布比较靠近模型前缘，后期随着气流的不断发展，气流将环绕模型中部，离散相高浓度在此处不断集结，最终造成霜层中部厚，模型两端较薄的现象。

4.2 霜层厚度分析

通过模拟得到了不同含量水蒸气的霜层云图，以冷表面温度110 K，空气流速0.5 m/s，凝华时间120 min为例，观察其结霜厚度变化情况，不同含量下的霜层厚度云图如图8所示。

通过霜层云图发现，不同浓度下的霜层生长不一，霜层中部有明显三处增厚，这与离散相浓度分布有关，气流一直波动，造成三个区域浓度高，同样，速度分布也会加剧这种浓度分布，水蒸气含量越多，现象越明显。从霜层增长变化看，增加水蒸气含量，模型表面结霜量增加，霜的厚度也有所增加。当水蒸气含量(质量浓度)为5×10-6时，霜层厚度仅为9 μm，平均厚度7 μm，表面结霜量较少，霜层极薄。结果表明，成倍增加水蒸气含量，厚度变化与含量呈线性相关，水蒸气含量增加至100×10-6，霜层厚度是 5×10-6环境下的20倍，霜层厚度为0.18 mm。此外，如此低水蒸气含量形成的霜层厚度和密度都比较小，对实验研究的测量造成较大的困难。图9为不同水蒸气含量下的霜层厚度变化图。

仅改变结霜工质，设置为二氧化碳，保持其他条件不变，图10为不同二氧化碳含量下的霜层厚度变化图。

从图10中可以得出，二氧化碳厚度增长变化和水蒸气相似，呈正比例增加，比较厚度变化值，二氧化碳在10×10-6含量下，厚度增加16 μm，与同等条件下的水蒸气相比，霜层厚度略厚些。在30×10-6，50×10-6，80×10-6，100×10-6的环境下，对应霜层厚度分别为52 μm，86 μm，137 μm，169 μm，每增加10×10-6含量，霜层厚度平均变化幅度值15 μm～20 μm。总体来看，随着含量的增加，二氧化碳凝华结霜几率也就会越大，单位面积内有更多的液滴凝结成霜，霜层覆盖率也就越大，结霜厚度变化也越明显。二氧化碳浓度越大，表面的霜层变化越明显，表面形态破坏程度也越严重。因此，需严格控制风洞实验室中二氧化碳含量，以便后续实验的正常开展。