超重力下相变材料熔化过程的数值模拟

赖柄竹，潘文彬，李 玲，邸志惠，陈 萱，许 玉

（南京航空航天大学航空学院飞行器环境控制与生命保障工业和信息化部重点实验室，南京210016）

随着机载电子设备的不断增多和功率的不断增大以及高能武器的应用，机载设备的热负荷不断增加，传统的风冷和液冷方式已较难满足散热需求[1]。近年来，利用相变材料（PCM）冷却发热元件的方法，得到了广泛关注[2-4]。作为一种被动散热技术，PCM的吸热与放热过程均在特定温度下进行，可使发热元件的温度保持恒定或近似恒定，从而确保其正常工作[5]。由于利用了PCM的大量潜热，因而该技术具有能量密度高、温度恒定和质量轻、体积小等优点，对于机载设备的热管理具有很大的应用潜力，不仅可以有效散热，而且能在脉动热流冲击等特殊情况下对设备进行短时过热保护[6]。然而，现代高性能飞机机动飞行时，通常处于超重力环境下[7-8]，例如三代机上的超重力可达6倍地球重力，四代机可达9倍，此时PCM的熔化过程很可能发生显著变化，继而导致换热性能发生改变，致使地球重力下所得到的设计方法不再适用。因此，有必要研究超重力对PCM 熔化过程的影响。

国内外相关研究[9-15]大都集中在提升PCM的换热性能上，这是因为PCM 的导热系数普遍较小，不利于设备的高效散热。一种有效的解决方式是在PCM 中添加导热增强器，如翅片、泡沫金属和金属颗粒等。Pakrouh 等[9]对含有不同翅片的PCM 热沉进行了数值模拟研究，并考虑了自然对流的影响，发现最优结构与PCM和翅片的体积分数存在复杂关系，且强烈依赖于翅片的数量、高度和厚度等。Arshad 等[10]实验研究了PCM体积分数对散热器换热性能的影响，发现当翅片体积占总容积的9%左右时，换热性能最好。Mahmoud等[11]实验研究了PCM种类和热沉结构等对电子设备冷却装置性能的影响，发现随着翅片数量的增加，热沉峰值温度降低。Baby 和Balaji 等[12]实验研究了含有33、72 和120 个翅片的PCM 散热器的换热性能，发现72 个翅片的散热器具有最高的换热性能，表明在一定条件下，翅片数量的增加并不一定能提高散热性能。盛强和邢玉明[13]对泡沫金属复合相变材料储放热过程进行了实验和数值模拟研究，发现泡沫铜的填充缩短了PCM 熔化和凝固的时间。林佳等[14]对某高功率电子设备模块进行了数值模拟，发现复合相变材料（膨胀石墨/高碳醇）热沉的散热性能显著优于铝热沉。何智航[15]模拟研究了某弹载天线热管PCM 复合热控装置的温控性能，发现使用较高导热系数的石墨基石蜡与较多的铜翅片均有利于提高装置的热控性能，且效益互补。

部分学者[16-21]研究了航天微重力环境下PCM的换热性能。Humphries和Griggs[16]研究了PCM在航空航天方面的应用，检查了PCM在地球重力和微重力环境下的换热性能，提出了相变蓄热装置的优化设计方法。Veilleux 等[17]模拟研究了电磁场对低重力环境下PCM（镓）熔化时换热性能的影响，发现电磁场可以增强换热，加快PCM的熔化速度。麻才新等[18]对一种空间相变换热器进行了模拟仿真，并对比了该换热器在地球重力和微重力条件下的性能差异。他们发现，在地球重力下，PCM 液体有明显的无规则流动扰动（不稳定Rayleigh-Benard 对流），可强化PCM 的换热性能。而在微重力下，这种流动扰动较弱，PCM 的传热阻力增大。Chen等[19]模拟研究了PCM 在地球重力和微重力条件下的熔化过程，结果表明，在微重力下，热传导起主导作用，而在地球重力下，自然对流支配换热过程，因而微重力下的换热性能弱于地球重力下的换热过程。此外，随着融化过程的进行，PCM 在2 种重力下的焓值之差逐渐增大。张靖驰等[20]通过数值方法研究了微重力下填充泡沫铜和泡沫碳对PCM 换热性能的影响，结果表明，泡沫复合相变材料可将热源热量有效分散到各个区域，减缓加热面温度上升的速度，并且能够降低重力变化对换热的影响。阮世庭等[21]通过数值方法探究了微重力下PCM 熔化过程的换热特性，结果表明，当受微重力影响时，PCM熔化速度明显下降，热量主要通过热传导传递。

综上所述，目前关于PCM换热性能的研究主要集中在地球重力和微重力条件下，且PCM在微重力下的性能与地球重力下有所不同，预示着在超重力环境下PCM可能存在更为明显不同。然而，目前对超重力条件下PCM换热性能的研究仍很少。鉴于此，本文通过数值模拟的方法研究超重力下PCM 熔化时的换热特性，并与地球重力下的换热特性进行对比，揭示超重力的影响。

1 物理模型

本文模拟所用的物理模型如图1 所示，其外形尺寸为0.1 m×0.022 m，底部为加热面。依据机载电子设备常见热流密度（q）范围[21]，选取q=50 kW/m2。重力（ah）大小分别为地球重力（1g）、5 倍地球重力（5g）和10倍地球重力（10g），重力方向垂直向下，即与热流的进入方向相反。

图1 物理模型Fig.1 Physical model

石蜡类PCM 具有相变温度范围宽、潜热高、可逆性好、无毒无腐蚀、价格低廉、在使用温度范围内发生燃烧的概率小等优点，是较为理想的机载电子设备散热用PCM[2]。石蜡RT56 的相变温度为329.15～331.15 K，符合一般机载电子设备正常工作的温度要求[4，10]，因而本文模拟时选用该材料，其具体物性参数见表1。

表1 PCM物性参数Tab.1 Physical parameters of the PCM

2 数学模型

采用基于焓法的凝固/融化模型进行模拟，主要控制方程如下：

式（1）～（3）中：ρ 为密度；t 为时间；u 为速度；μ 为动力粘度；k 为导热系数；P*为压力；T 为温度；H 为焓；Sm为动量源项；Se为能量源项。

式中：h 为显焓；ΔH 为潜焓。

式中：Tref为参考温度（PCM 初始温度）；href为Tref对应的焓；cp为定压比热容。

式中：β 为液相率；L 为相变潜热。

式中：Tsol为固相温度；Tliq为液相温度。

式中：ε 为一个很小的计算常数，避免公式除0，取0.000 1；Amush为固液模糊区常数，取105；Sb为浮力源项。

Sb满足Boussinesq 假设，即除浮力项外所有项中的密度不变，浮力项中密度随温度呈线性变化。

式中：ρref为Tref对应的密度；αV为热膨胀系数；ah为重力加速度。

从式（9）中可以看出，ah变化会改变浮力源项Sb，继而改变动量方程（2），从而对PCM熔化过程产生影响。

3 网格划分和求解设置

3.1 网格划分

如图2 所示，采用Pointwise 18.0 将模型划分成四边形结构网格。

图2 局部网格Fig.2 Partial meshes

3.2 边界和初始条件

加热面的热流密度设置为q=50 kW/m2，而容器左、右侧壁和顶面均设置为绝热条件，即q=0 kW/m2。ah大小设置为ah=1g、5g 和10g，ah方向设置成垂直向下。初始时刻，PCM 的温度设定为329.15 K，此时PCM处于固态。

3.3 求解设置

采用Fluent 17.0进行非稳态计算，选择凝固/熔化模型，开启能量方程。根据瑞利数（Ra）判断自然对流的强度，即Ra ＜108为层流，Ra ＞1010为湍流，其定义[22]如下：

式（10）中：ΔT 为温差；l 为特征尺寸；α 为热扩散系数。

压力-速度耦合选择SIMPLE 算法，压力离散采用标准算法，动量和能量的离散格式设置为二阶迎风格式。收敛标准为所有的残差均为10-5。时间步长取0.005 s。

4 模拟结果

4.1 地球重力（1g）时熔化特性

图3 1g时加热面温度和液相率动态变化Fig.3 Dynamic change of bottom temperature and liquid fraction under 1g

图4展示了1g情况下，t=22 s 时的固液相分布和速度场。从中可以看出，在液相区域内，底部率先熔化的液体温度较高，密度较小，而上方液体温度较低，密度较大，在浮力作用下，底部液体逐渐上升，上方液体逐渐下降，如此形成对流，此时固液交界面出现一定程度的波动。

图4 1g时固液相分布和速度场Fig.4 Distribution of solid and liquid phases and velocity field under 1g

4.2 5倍地球重力（5g）时熔化特性

图5 展示了5 倍地球重力下（5g），模型的Tw和β的动态变化趋势。从中可以看出，5g 时Tw和β 的变化趋势与1g时的类似，但是PCM熔化速度明显加快，不同时刻的Tw均显著降低，在t=114.775 s 时，PCM即完全熔化。此时，Tw=484.51 K ，β=1。此外，在t=2.985 s ，第一个过热温度达到峰值，此时Tw=427.31 K ，β=0.024 1 ；在 t=5.15～9.31 s 时，β=0.049 5～0.087 1，Tw的波动幅度相对于1g时增大。

图5 5g时加热面温度和液相率动态变化Fig.5 Dynamic change of bottom temperature and liquid fraction under 5g

图6展示了5g条件下，t=22 s 时的固液相分布和速度场。从中可以看出，此时液体的运动速度相对于1g 时明显增大，对流作用明显增强，固液交界面波动幅度增大。

图6 5g时固液相分布和速度场Fig.6 Distribution of solid and liquid phases and velocity field under 5g

4.3 10倍地球重力（10g）时熔化特性

图7 展示了10 倍地球重力下（10g），模型的Tw和β 的动态变化趋势。从中可以看出：10g 时，Tw和β的变化趋势亦类似于1g 和5g 时的，但PCM 熔化速度相对于5g 时又进一步加快，不同时刻的Tw进一步降低，在t=100.955 s 时，PCM 已完全熔化。此时，Tw=449.88 K，β=1。值得注意的是：10g 时，PCM 熔化速度加快的幅度没有从1g 到5g 时的幅度大；与此同时，Tw下降的幅度亦比从1g到5g时的小。此外，在t=2.595 s ，第一个过热温度达到峰值，此时Tw=415.28 K ，β=0.021 3 ；在t=4.66～8.485 s 时，β=0.048 7～0.085 8，Tw亦有波动，其波动幅度比5g时略大。

图8 展示了10g 情况下，t=22 s 时的固液相分布和速度场。从中可以看出，此时液体的运动速度相对于5g 时进一步增大，对流作用进一步增强，固液交界面波动幅度继续增大。

图7 10g时加热面温度和液相率动态变化Fig.7 Dynamic change of bottom temperature and liquid fraction under 10g

图8 10g时固液相分布和速度场Fig.8 Distribution of solid and liquid phases and velocity field under 10g

5 结论

本文通过数值模拟的方法，对PCM在地球重力和超重力下的熔化换热过程进行了研究，主要结论如下：

1）由于PCM导热系数过低，导致熔化初期加热面温度出现过热。在1g、5g 和10g 时，初期最大过热温度分别为480.02 K、427.31 K 和415.28 K，十分不利于设备散热。须通过耦合其他强化换热措施，如翅片或高导热系数的纳米尺度添加物，来减弱甚至消除该种过热。

2）超重力对于PCM 熔化过程的影响十分明显。无论是熔化速度，还是加热面温度，它们的变化趋势均随着超重力的变化而出现显著变化。在1g、5g 和10g 时，PCM 完全熔化所需时间依次为158.37 s 、114.775 s 和100.955 s，完全熔化时的加热面温度依次为585.48 K、484.51 K 和449.88 K。

3）超重力对PCM 熔化过程的影响可归结为随着超重力的增加，动量方程中浮力项的作用不断增大，使得液相PCM的自然对流不断增强，继而导致熔化界面的波动幅度增大，熔化速度加快，加热面的温度下降，十分有利于设备散热。

4）超重力对PCM熔化过程的增强幅度，随着超重力的增加而减小。加热面温度的下降幅度由1～5 g 时的100.97 K 变为5～10 g 时的34.63 K，PCM 完全熔化所需时间的缩短幅度由1～5g 时的43.595 s 变为5～10 g 时的13.82 s，表明超重力的强化换热作用有一定限度，当超重力超过一定值，其作用相对弱化。

5）在设计机载环境下使用的电子设备PCM 散热系统时，可充分利用超重力对PCM熔化过程的强化作用，实现同等散热性能下，系统体积和重量的减小。