基于FLUENT 的水下热滑翔机相变过程仿真分析

刘鸿瑨，胡欲立，杨 威，郝泽花

(西北工业大学 航海学院，陕西 西安，710072)

0 引言

我国海洋资源十分丰富，人们对海洋资源的探索开发离不开各种水下航行器[1]。水下热滑翔机是一种新型无人自主水下航行器，其利用相变材料在固液相变过程中产生的体积变化，将海洋温差能转化为机械能，并完成水下的之字形运动。

国内学者很早就对水下热滑翔机展开了研究，孔巧玲[2]、倪园芳[3]、任龙飞[4]等建立一维的水下热滑翔机热管换热模型，采用Matlab 软件，利用有限容积法对滑翔机上相变材料的固液相变过程换热特性展开了研究，计算过程中利用等效导热系数的方法来表现凝固融化过程中自然对流作用的影响，计算过程并未扩展到二维或三维情况。王延辉[5]和Yang 等[6]对温差能驱动的水下滑翔机展开了动力学分析与设计，并设计出了一种带内柔囊的换热热管，但是并未对热管内相变换热过程温度场及液相率的分布进行深入研究。梁泽德等[7]使用正十六烷和石墨为9：1 组成的复合相变材料，研究对比了不同尺寸圆筒形和圆柱形竖直换热容器在不同温度下的相变换热特性，但是并未考虑到压力对相变过程的影响。李国道[8]以石蜡为研究对象，采用添加膨胀石墨的方法来提高石蜡换热性能，并仿真计算了水平圆管内的相变过程。杨锐[9]以海洋浮标上低温固液相变材料为研究对象，采用实验的方法研究了添加改进剂对正十六烷相变温度的影响。从前人研究内容来看，对水下热滑翔机相变过程研究主要集中在数值仿真、热管结构以及相变材料几方面。

在此基础上，文中从融化温度的角度，引入了压力对相变过程的影响。建立了二维水下热滑翔机换热热管中的相变换热过程仿真模型。为方便模型建立以及网格绘制，同时也为得到自然对流作用对相界面的影响，使用FLUENT 软件对水下热滑翔机热管内相变过程温度及液相率的变化进行了仿真分析，对水下热滑翔机温差能动力系统设计改进具有积极意义。

1 相变换热模型

1.1 物理模型

图1 为热管截面相变过程示意图，管内部存放相变材料，整个圆筒处于海水环境中。材料初始状态为液态，当管外海水低于材料相变温度时，相变材料开始凝固；当管外海水温度高于材料相变温度时，管内相变材料开始融化。热滑翔机所用相变材料为正十六烷，其物性参数受温度影响较小。

为简化计算模型进行如下处理：

图1 热管结构示意图Fig.1 Schematic diagram of heat pipe structure

1) 假设固液相变材料正十六烷各向同性，固液两相的物性参数均为常数，不随温度发生变化，比如热导热系数、粘度等；

2) 相较于相变材料，管壁材料的导热系数很大，管壁厚度可忽略；

3) 液态工质温度分布不均会使其密度分布不均匀，进而产生自然对流作用，在仿真计算过程中密度选项选择布辛奈斯克(Boussinesq)近似，其作用是在动量方程中浮力项中，将原本为常数的密度改为与温度有关的函数。

由于水下热滑翔机热管为细长圆管，温度随管长变化不明显，故截取横截面分析，使用二维结构网格。网格用ICEM CFD 绘制，圆面直径为24 mm。经网格无关性检验，取网格单元数量9 396，节点数量9 289 时，计算能保证较高的精度和较快的计算速度。将网格导入FLUENT 18.0，选择二维分离隐式双精度求解器，并将圆截面中心点作为观测节点，监测该点温度变化。热滑翔机在水下以缓慢速度滑翔，其热管与周围海水之间为对流换热，假设滑翔机运行速度为0.28 m/s[10]，经计算取对流换热系数为324 W/(m2·K)，仿真计算时间步长取0.1 s。

正十六烷在固液相变过程体积变化大，且固态密度大于液态密度，其融点刚好在海洋温跃层温差范围内(4～25℃)，比较适合作为温差能热机的相变工质。文中采用正十六烷为相变材料，其物性参数如表1 所示。

1.2 数学模型

FLUENT 软件中的Solidification/Melting 模型[11]可以用来仿真材料相变过程，模型中能量方程采用焓法[12]求解，在固态、液态和固液糊状区建立统一的能量方程，不需要跟踪相界面的移动，以简化计算过程。

表1 正十六烷物性参数Table 1 Physical parameters of n-hexadecane

式中：H为单位体积相变材料的总焓；ρ为材料密度；υ为流体速度；S为动量源项；k为导热系数；ν为运动粘度；p为压强。

焓孔隙率技术将糊状区域(部分凝固区域)作为多孔介质，每个微元内的孔隙度等于该孔内的液体分数。在动量方程添加动量源项

式中：f为液相率，在液体区f=0，在固体区f=1，在糊状区域f在0～1 之间变化；ε为常数，作用是防止仿真计算过程中式(4)分母变成为0，通常取一个极小数值0.001；Amush为糊状区常数，一般取1×105。

FLUENT Solidification/Melting 模型中液相率f的定义如下

对于纯物质，其融化温度lT等于其凝固温度sT，但在计算过程中为得到较好的收敛性，2 个温度之间通常取一个极小的差值，文中在后续计算过程中差值取0.5 K。

焓法求解能量方程的过程中，液相率是用温度定义的，故须建立液相率和工质焓的关系。在计算过程中，工质焓定义为

式中：L为材料相变潜热；href为参考焓值；Tref为参考温度；cp为相变材料的比定压热容。

根据以上公式建立焓和温度的关系。对于相变过程中材料的导热系数采用如下表达形式

式中，ks、kl分别为材料固相和液相的导热系数。

材料比热容示意见图2，具体定义如下

图2 不同区域比热容示意图Fig.2 Schematic diagram of specific heat capacity of different regions

1.3 压力对相变过程的影响

热滑翔机相变材料在融化时体积膨胀会压缩蓄能器中的气体，其动力系统相当于一个压力系统，而压力变化又会改变相变材料的融化温度，因此研究压力对相变过程的影响很有必要。

Ma 等[13]通过分析压力对相变材料和液压油中气体溶解度的影响来研究不同压力条件下动力系统的存储效率；Xia 等[14]通过在能量方程加入压力做功的源项，建立了一维相变换热模型，发现由于压力做功相对于相变潜热非常小，压力对相变过程的影响也非常小。

材料发生固液相变时在融化温度处是一个两相平衡的状态，依据克拉配龙方程[15]，物质在固液两相共存时，其熔化温度和压力存在如下关系

式中，ρl和ρs分别为材料液相和固相的密度，进一步积分可得

式中，下角标1、2 分别代表相变材料的不同温度压力状态。

依据上式即可计算得到材料在不同压力条件下的融化温度，如图3 所示。

图3 融化温度随压力变化曲线Fig.3 Curve of melting temperature versus pressure

2 仿真结果及分析

2.1 融化过程

计算中设置管内的相变材料初始温度为16℃，此时低于正十六烷的融化温度，正十六烷处于凝固状态。管外设置为恒壁温边界条件，海水温度为25℃，热管从周围海水中吸收热量，温度升高，当温度高于正十六烷融化温度时，正十六烷开始融化过程。

从图4 和图5 可以看出，随着时间的推进，热管外层相变材料最先开始融化，并缓慢地向热管中心推进。随着压力的上升，热管内相变材料融化速度明显变慢。在4 000 s 时5 Mpa 压力条件下，已有大约85%的相变材料完成融化；而相同时间30 Mpa 压力条件下，只有不到30%相变材料变成液态。这是因为，随着压力升高，材料融点升高，相变材料需要吸收更多的热量来达到相变温度，在换热条件相同的条件下就需要花费更长的时间。融化过程中，相界面开始能保持为圆形，后面由于液体部分的增加，自然对流作用使得相界面形状发生了改变。

图5 4 000 s 时融化过程液相率分布Fig.5 Liquid fraction distribution of melting process at 4 000 s

图6 为融化过程中不同压力条件下液相率的变化，开始时刻材料全部为固态，此时液相率的值为0。当此值变为1 时，表明热管内材料已全部变为液态。从图中可以看出，压力越高，曲线越靠近横轴，液相率变化梯度也越小。在0.1，5，10，15，25，30 MPa 压力下，材料完全融化所需时间分别约为64，72，86，101，179，271 min，其中5，15 MPa 条件下的融化时间比0.1 MPa 条件下融化时间分别长了12%，34%，而35 MPa 条件下，在300 min 时融化过程仍未全部完成。

图6 不同压力下液相率变化曲线Fig.6 Curves of liquid fraction at different pressure

图7 为不同压力下中心点温度的变化。初始时刻，中心点温度上升很快，此时热量主要以显热的形式存储在材料中，当达到融化温度时，材料温度不再上升，热量以相变潜热的形式存储，材料开始进行融化过程。当材料融化为液态后，由于外界温度仍然高于此时液态材料，温度继续上升。

图7 中心节点温度变化曲线Fig.7 Temperature change curves at the center node

对比图6、图7 及液相率分布图可知，当中心点完成融化时，热管内的整个融化过程并未全部完成，这正好与云图结果相吻合。出现这种情况，是因为固态正十六烷比液态正十六烷密度大，在融化过程中，热管中心的固态正十六烷会下沉，同时由于液体部分自然对流的作用，热管内相变材料最后融化的位置是在中心节点的下方。

2.2 凝固过程

为研究压力对正十六烷凝固过程的影响，仿真计算时设置热管内相变材料初始温度为25℃，管外设置为恒壁温边界条件，海水温度设置为16℃。图8、图9 分别为凝固过程1 000 s 和4 000 s时热管内液相率分布云图。

图8 1 000 s 凝固过程液相率分布Fig.8 Liquid fraction distribution of solidification process at 1 000 s

对比图8 和图9 可知，材料最先从外部开始凝固并逐渐向内部发展，压力越高凝固速率越快。同一压力条件下，周向上凝固速率相同，整个过程相界面都保持为圆形。

图9 4 000 s 凝固过程液相率分布Fig.9 Liquid fraction distribution of solidification process at 4 000 s

图10 为凝固过程中液相率的变化。初始时刻材料全部处于融化状态，此时液相率的值为1。当此值变为0 时，表明材料已全部凝固。从图中可知，随着压力升高，曲线越靠近纵轴，下降梯度也越来越大，材料凝固速率加快。在0.1，5，10，15，25，30，35 MPa 压力下，正十六烷完全凝固所需时间为333，252，203，168，124，113，101 min，其中5，15 MPa 压力下凝固速率分别比0.1 MPa条件下凝固速率提高了24%和50%。

图10 凝固过程液相率变化曲线Fig.10 Change curves of liquid fraction during solidification process

图11 为凝固过程中心节点的温度变化。对比图10、图11 可以看出，在中心节点凝固完成时，整个热管内相变材料也都完成了凝固。

图11 凝固过程中心节点温度变化曲线Fig.11 Temperature change curves at the center node during solidification process

3 结束语

计算分析表明，压力会改变材料的融化温度进而影响材料相变过程。压力升高使得材料融化温度上升，在相同温度条件下完全融化所需时间会变长，但是凝固所需时间反而变短。如果融点高于了海水暖水层温度，相变材料将无法完成相变过程，热滑翔机将不能稳定工作。文中每个算例仿真计算过程中压力都是取一定值进行计算，实际相变过程中压力是受体积变化的影响，并处于一个连续变化的状态，如何实现压力变化和相变过程耦合仿真是一个值得研究的方向。下一步的工作重点将是建立压力和相变过程的耦合仿真，使模型更接近于实际的滑翔机上的相变过程；搭建滑翔机相变过程仿真平台，实测相变过程的压力变化和相变时间，并与仿真结果对比分析。