水升华器升华模式理论分析与数值仿真

廖俊元杨春信杨涵

(北京航空航天大学航空科学与工程学院，北京 100191)

1 引言

随着中国航天技术的发展及载人登月工程的逐步推进，航天员的舱外活动需求逐渐提高。舱外航天服的热控系统需要保证航天员在出舱活动时处于安全、舒适的热环境之中。水升华器是一种利用水的升华潜热进行散热的消耗型相变散热装置，相比于辐射器等热控器件具有体积小、质量轻、效率高、微重力下工作可靠等优点，是一种性能优良的空间热沉，十分适用于应对小散热面、短时大功耗任务的舱外航天服的热控系统。

国际对水升华器的研究开展较早，从上世纪60年代初已有关于水升华器作为热沉应用于航天器热控系统的相关研究，并制作了如LM-107、LM-209等一系列原理样机进行了初步的实验探究，如图1所示。在上世纪70年代初至上世纪末的期间，水升华器处于大量工程应用的阶段，在首次成功应用于美国的阿波罗探月工程中后，还在苏联(俄罗斯)的Orlan航天服等航天工程中得到应用。进入新世纪后，为了满足更高的航天工程需求，水升华器出现了追求防污染、良好动态性能的发展趋势，对此NASA的研究人员提出了X-38、CIS、SDC、ISDC 4种水升华器的新构型，如图2所示。

图1 早期水升华器Fig.1 Early application of sublimator

图2 新型水升华器Fig.2 New type of sublimator

中国对水升华器的相关研究起步较晚。吴志强等最早对水升华器进行了散热性能的理论分析及一系列实验研究；李森等建立了多孔板单孔内水相变过程的物理数学模型，并进行了周期模式的数值模拟；王玉莹等采用一维热阻网格模型，通过仿真得到了周期模式下界面位置、水升华器温度等参数的变化，并对不同加热边界的水升华器工作性能进行了实验研究。航天工程方面，中国关于水升华器的首次成功应用为2008年神舟七号飞船出舱活动的“飞天”舱外航天服。

由于水的蒸发、升华等相变过程发生在水升华器内部的给水腔、多孔板微孔等尺寸微小的区域，实验仅能测得外部温度等宏观数据，因此对于水升华器工作原理的研究主要通过理论分析、数值仿真来进行。但是，目前国内对水升华器的仿真分析较少且主要为一维数值仿真，没有考虑水升华器不同工作模式的转换。为了探究水升华器工作模式转换的规律，分析不同设计参数对水升华器工作特性的影响，本文在总结前人理论分析的基础上，对水升华器工作模式转换特性进行量纲分析，通过Fluent数值仿真的方法对水升华器在低热载荷下的稳态工作模式(升华模式)开展研究，分析不同多孔板结构参数、热载荷大小下水升华器工作特性的变化规律。

2 水升华器工作过程分析

2.1 水升华器工作原理

本文研究对象为平板型水升华器，其基本结构由供水入口、加热板、给水腔、多孔板组成，如图3所示。

图3 升华模式下的水升华器Fig.3 Water sublimator working in sublimation mode

系统工作过程为:工质水在一定的供水压力下被送往高真空度的给水腔，水的压力骤降至其三相点以下而迅速蒸发。蒸发潜热带走大量热量，若此时水温降至三相点温度以下，则会凝固形成固态冰。外表面的冰层暴露于高真空环境下，热载荷将通过升华潜热排向外太空。

Hamilton Standard公司根据工质水的不同消耗形式，将水升华器工作模式划分为4种:升华模式、蒸发模式、混合模式、周期模式。当热载荷较小时，水升华器工作在升华模式。由于水蒸气流量较小，升华界面对应的饱和蒸气压低于水的三相点压力。根据图4，此时水不会进入多孔板的毛细孔，而是首先在多孔板与给水腔的界面上形成一定厚度的冰层(图3)。此时，供水流量等于升华流量，升华带走的热量等于底面加热及供水热量之和，水升华器内部达到整体的热质动态平衡。升华模式是水升华器最期望的工作模式[12]。

图4 水的三相图Fig.4 Three-phase diagram of water

当热载荷继续增大，水将进入多孔板形成一定厚度的冰层。冰层在升华、熔化的共同作用下逐渐耗尽后，水继续流入并凝固成冰层，形成了周期性的循环工作过程，即周期模式。若多孔板温度较高，则冰层不会出现，工质将以液态水蒸发的形式进行散热，即蒸发模式。水升华器实际工作中可能会出现以上3种模式共存的情况，这种模式被称为混合模式。

2.2 升华模式分析

吴志强等提出升华模式下水升华器整体处于稳态工作的状态，可以通过理论分析计算升华热流、加热面温度等参数。假设如下:

1)多孔板、加热板及工质水的物性为常数；

2)不考虑结构漏热、辐射漏热等热损失；

3)多孔板毛细孔简化为众多同径的竖圆孔；

4)忽略给水腔内部水对流造成的换热；

5)加热板、给水腔、多孔板具有相同的面积。

升华模式下水升华器内部的热质平衡如图3中展示，图中红色虚线箭头代表热量的流动，蓝色实线箭头代表质量的流动。定义加热板热流密度如式(1):

式中，

q

为底面加热热流，W/m；

Q

为水升华器热载荷，W；

A

为加热板面积，m。进一步对液态水、冰层分别进行热平衡分析，可得式(2)、(3):

式中，

q

、

q

分别为水层、冰层的导热热流，W/m；

m

为供水流量，kg/s，在升华模式下等于总体升华流量；

t

、

t

分别为供水温度、水-冰相变温度，K；

c

为水的比热，J/(kg·K)；Δ

h

、Δ

h

分别为水-冰相变潜热及冰-水蒸气升华潜热，J/kg。同时，还可以分别写出水层与冰层的导热关系如式(4)、(5):

式中，

λ

、

λ

分别为冰和水的导热系数，W/(m·K)；

t

、

t

为升华温度、底面温度，K。需要说明的是，式(4)仅考虑了水层的稳态导热，而忽略了入口较高温度供水的影响。王玉莹在相关仿真中做了类似的处理，其仿真与实验结果较为一致，可以认为该处理方式的误差在可接受范围内。

此外根据式(6)所示，升华模式下水层、冰层的厚度之和与给水腔厚度相等:

式中，

d

、

d

、

d

分别为冰层厚度、水层厚度、给水腔厚度，m。根据式(3)，可推导出水升华器消耗工质水的质量流量如式(7):

式中，

ρ

为水的密度，kg/m；，

S

、

u

分别为供水入口的面积与供水水流速度。由于升华质量流量很小，水蒸气分子平均自由程远大于流动特征长度，具有较大的克努森数

Kn

，可以认为处在自由分子流的流动状态。自由分子流状态下的质量流量为式(8):

式中，

n

为多孔板单孔个数，

ϕ

为多孔板孔径，m；

d

为水蒸气流动长度，m；

m

为单个水分子质量，kg；

k

为玻尔兹曼常数，

p

为真空环境压力，Pa；饱和蒸气压

p

可根据式(9)给出:

此外，将式(8)带入式(3)中，可得水升华器的升华热流如式(10):

由于水不进入多孔板，此时水蒸气的流动长度

d

即多孔板厚度

d

。联立式(7)～(9)可以数值求解升华温度

t

，并根据式(4)～(6)计算冰层厚度

d

、水层厚度

d

及底面温度

t

。若热载荷大于一定值，使升华表面蒸气压

p

大于水的三相点压力

p

、升华温度

t

超过三相点温度

t

，水会进入多孔板转换为周期模式。存在一个水升华器以升华模式工作的临界热载荷，标志升华模式向周期模式的转换，基于式(7)～(8)可得式(11):

当水升华器的热载荷小于式(11)时，认为处于升华模式；而当热载荷大于该值时，则认为水升华器处于周期模式下。此外，升华模式下热载荷还需要满足式(12)的限制，避免给水腔被冰层完全填充可能造成的损坏。这个最低热载荷的数值同样可以通过数值方法计算得出。

式中，

d

｜代表热载荷

Q

下根据导热关系式及式(6)计算得出的冰层厚度。

综上，对于尺寸确定的水升华器，其升华模式稳定工作的最高与最低热载荷均是一个常数，可以通过式(11)～(12)确定。取升华面积0.04 m、多孔板厚度1 mm、给水腔厚度3 mm，可得图5所示的结果，多孔板孔度、孔径的增大会使热载荷范围变大、总体数值升高。

图5 升华模式工作范围Fig.5 Working scope of sublimation mode

2.3 热载荷公式量纲分析

将式(11)中的小孔个数

n

用孔度

ε

展开并移项，两边再同除质量流量

m

，可得式(13):

将上式进行量纲分析可以得到式(14):

式中，将式(13)的左边将临界热载荷

Q

替换为实际热载荷

Q

，可定义斯坦顿数

St

、雅各布数

Ja

倒数差为式(15):

式中，

ρ

为水蒸气密度。定义为一个新的无量纲数

Su

如式(17)，包含水升华器的孔度、孔径、面积等所有设计参数:

结合密度、速度、面积比项，定义无量纲质量流量如式(18):

根据2.2节的内容，当热载荷大于一定值时水升华器将转入周期模式。由于此前的分析保留了模式转换对应的三相点参数，因此热载荷

Q

取临界热载荷

Q

时，可将式(11)转化为如式(14)所示的无量纲形式，反映临界热载荷与水升华器设计参数之间的关系。对于孔度、孔径及热载荷分别在0.3～0.7、1～9μm、100～300 W间随机取值，可得到如图6所示的水升华器工作模式分布及模式转换的无量纲特性曲线。对于水升华器不同设计参数组合的无量纲数

Su

，若热载荷对应的

St

取值在图中的特性曲线以上，可以认为水升华器处于周期模式下工作，反之则认为进入了升华模式。此外，由式(12)决定的升华模式下边界曲线同样可以在图中给出。

图6 模式转换无量纲特性曲线Fig.6 Dimensionless curve of mode transition

3 仿真模型

为了进一步探究水升华器升华模式的具体特性，基于2.2～2.3节理论分析，进行水升华器升华模式的数值建模与仿真计算。

3.1 焓-多孔模型

本文采用焓-多孔模型求解固-液相变问题。该模型以温度和焓作为共同变量，在全计算区域内建立统一形式的能量守恒方程。温度和焓作为共同因变量的关系可以式(19)用来表示:

引入糊状区(Mushy zone)的概念，将温度处于固相线温度与液相线温度之间的区域视为糊状区，该区域内网格的液相分数视为孔度，通过不断更新所有网格单元内的孔度追踪固液界面的位置。液体分数的表达为式(20):

对于水-冰相变问题，其固相线与液相线温度相等，均为273.15 K。糊状区的动量表达为式(21):

式中，

β

为液体分数，

ε

为防止分母为0的小量，

A

为糊状区常数，

V

为速度矢量，

V

为固态离开原区域的漂移速度。基于式(20)～(21)，可得相变问题的控制方程组如式(22)～(24)所示。

连续性方程为:

由于水升华器工作在失重的空间环境中，可以忽略重力的影响，得到动量方程:

能量方程为:

式中，总焓

H

为显焓

h

与潜热焓Δ

H

之和，关系如式(25):

式中，

h

为参考焓值，J/kg；

T

为273.15 K的参考温度；

L

为水-冰相变潜热，J/kg。

3.2 网格划分及Fluent求解模型

使用商业软件ANSYS ICEM CFD及ANSYS Fluent 2020 R1版本作为网格划分工具及求解器。由于供水流量很小，流动模型选择层流Laminar。能量方程的求解调用Fluent中基于焓-多孔模型的Solidification/Melting模型，选择基于压力的隐式稳态单精度求解器进行求解，当整体热量达到平衡后计算结束。

仿真分析相同外形的水升华器在不同孔度、孔径时的升华模式特性，水升华器结构参数及相关常数如表1所示，其中3个工况已在图7中标出。

表1 水升华器参数及相关常数Table 1 Parameters of sublimator and relative constants

对给水腔进行二维矩形结构网格的划分，并且选取了数量为7150、8400、10500的网格进行仿真网格无关性验证计算。3套网格的工况2计算结果如图7所示，不同网格数量下的计算结果最大差异不超过0.15%。

图7 网格无关性验证Fig.7 Grid independence validation

综合考虑计算的准确性及计算资源的占用，最终选取600×14共8400数量的网格，网格划分及边界条件设置情况如图8所示。其中，底面热载荷所在的加热面设置为等热流密度边界，为根据图6确定的热载荷区间选取的合适的热载荷。顶面为升华界面，定义为由温度决定的热流边界，与孔度、孔径及升华温度等有关，如式(9)～(10)所示。计算域右侧加入一个由式(7)确定的等质量流量入口，同时在顶面施加一个质量源项UDF模拟升华带走的质量流。边界条件结果如表2所示。

图8 网格与边界条件Fig.8 Mesh and boundary condition

表2 不同工况边界条件取值Table 2 Boundary condition settings

4 结果与分析

由于Fluent中是基于式(26)计算液体分数的，松弛因子

γ

的存在会导致计算域中出现糊状区。为避免不确定性，采用

t

=273.15 K等温面表示水-冰界面如式(26):

三工况整体分布较为相似，取工况2的计算结果进行分析。整体及入口附近液体分数云图如图9所示，红色代表液态水，绿色、蓝色代表相变形成的糊状区及固态冰。可见受到温度较高的供水的影响，入口附近几乎没有形成冰层；而向后直至给水腔尽头，入口效应逐渐消失，形成了厚度较为一致且水-冰界面与上下边界大致平行的薄冰层。

图9 液体分数云图Fig.9 Contour of liquid fraction

工况2的流场矢量图如图10所示。受固液界面位置的影响，水流从入口进入后仅在下侧的液态区域中流动，上侧的固态冰层区域中不存在流动；随后水流继续向前流动远离入口，速度逐渐降低并最终停止流动。此外，还能观察到在冰面及底部壁面附近的流动边界层效应。

图10 速度矢量图Fig.10 Vector of velocity

工况2的温度分布云图如11所示。可以发现，温度分布与液体分数分布具有很强的关联性。同样将入口附近放大，可见入口的高温水流前进的同时很快地受到了冷却，等温线逐渐倾斜；入口效应逐渐消失后，受到上下热流边界的影响，等温线也逐渐平行于上下边界。此外，还能观察到右下角区域在底面加热作用下达到了全域最高温度，已超过了293.15 K的供水温度。

图11 温度云图Fig.11 Contour of temperature

统计各工况计算稳定后的底面温度、升华温度与其对应的冰层厚度的仿真结果与理论计算值对比，如表3所示。各工况的温度结果相对误差都在0.5%以内，底面温度仿真结果整体比理论值高1 K左右，而升华温度仿真结果比理论值稍低。造成这种误差的原因是计算域顶部升华热流、升华温度耦合的边界条件设置，其热流是由升华界面平均温度决定的，因此入口高温水流对升华界面的边界有着不可忽略的影响。

表3 温度、界面位置结果对比Table 3 Comparison of temperature and thickness of ice

各工况冰层厚度的计算结果与理论值相对误差的绝对值均小于1%。同样受到入口效应的影响，计算域右侧没有形成冰层，同时温度的误差也会影响冰层厚度，因此冰层平均厚度整体较理论值稍低，对应的相对误差偏向于负值。

5 结论

1)分析推导得到了水升华器工作模式转换的特性曲线，该曲线能给出不同热载荷及孔度、孔径等结构参数下水升华器的工作模式；

2)对不同孔度、孔径共3种工况进行仿真并与理论值比较，各温度结果的相对误差均在0.5%以内、冰层厚度的相对误差绝对值均在1%以内，仿真模型通过了理论验证；

3)入口附近较高温度的供水会造成等温线倾斜、无冰层形成等入口效应及冰层厚度仿真结果偏低等误差。

本文仅研究了低热载荷下的升华模式，可在后续研究中开展对热载荷较大时的周期模式的理论分析及数值仿真工作。