自然对流影响结冰的数值模拟及实验研究

刘仍通，潘 阳

(华东交通大学土木建筑学院，江西南昌330013)

冰蓄冷系统就是利用了低谷期的电在夜间满负荷开启制冷机组，并由蓄冷设备将冷量以冰的形式蓄存起来，在白天用电高峰期将冷量释放出来承担空调负荷。冰蓄冷系统包括动态蓄冰和静态蓄冰，其中静态蓄冰包括盘管式和冰球式蓄冰。许多学者对盘管式蓄冰进行了实验和理论研究。由于有自然对流和水的过冷现象的存在，再加上蓄冰过程是一个相变过程，随着时间的推移相变界面一直在变化。因此蓄冰过程是一个很复杂的传热问题。通常在处理浮力驱动时引入Boussinesq假设。然而，对于非Boussinesq流体(如水、镓、铋、液氮等)在特定的温度下有最大的密度值，如果在数值计算中，仍然采用线性的Boussinesq假设，则不会反映流体的实际流动过程。

N H Saeid[1]和Masaaki[2]研究了水的密度反转对自然对流的影响。刘道平[3-5]等人通过实验研究了不同管束的蓄冷特性，通过实验测得的水温和冰层轮廓，确定水的自然对流对管外蓄冰过程有影响。而Motoi Yamaha[6]对盘管式蓄冰桶性能进行了实验研究，指出了当底部水温降到4℃后，蓄冰桶上部的水温迅速下降，也是受水的自然对流的影响。白博峰等[7-8]也对矩形腔内的水的凝固进行了实验研究，指出在液相区也出现了流动反转现象并对此现象进行了解释。本文对竖直管外结冰过程进行了实验研究，并利用数值模拟软件对其结冰过程进行了数值模拟，对比了实验测得的温度和数值模拟计算的温度值，了解自然对流对结冰过程的影响。

1 理论模型建立

1.1 蓄冰物理模型

如图1的蓄冰桶，铜管中通入低温冷媒(-5℃的乙二醇)溶液，冷媒通过管壁与水进行热交换，使得靠近管壁的水温度降低，直到温度低于水的凝固点，水开始结冰。因为在实际中，由于水存在过冷度，所以即使当水的温度降到凝固点也不会结冰。因此，进行数值模拟需对实际问题进行假设，使得问题简单化。现做如下假设:

(1)乙二醇进出口温度相同恒为-5℃;(2)蓄冰桶的保温很好，认为其向周围环境散热忽略不计;(3)蓄冰桶内初始水温均匀一致;(4)除了浮力项中的密度外，其余的物性参数都为常数;(5)在结冰过程中，不存在过冷度，即水温达到0℃水开始结冰。

1.2 数学描述

如图2，假设r方向流动速度为u，z方向流动速度为v，控制方程[9]如下

式中:

β—液相分数;ε—不为零的常数;Amush—模糊区常数;

μl—动力粘度，P a◦s;g—重力加速度，9.8m◦s-2;

能量方程:

式中:L—潜热值;S—源项;H—焓值;Keff—有效导热系数;kl，ks分别为液相和固相的导热系数，w◦(m◦k)-1;cl，cs分别为液相和固相的比热容，kJ◦(kg◦k)-1

由于水不是Boussinesq流体，其在4℃时有最大的密度值，因此使用Boussinesq假设就不能反映实际的流动情况。为此，选择许多学者模拟水自然对流换热时使用的密度随温度的变化关系式:

为了观察结冰过程中自然对流对温度分布的影响，在蓄冰桶中设置9个点，距管壁5 mm设置3个点1～3;距管壁10mm的3个点4～6;距管壁30mm设置的3个点7～9。如图3。

图1 物理模型示意图

图2 简化模型坐标示意图

图3 温度测点示意图

2 数值模拟及实验结果及分析

使用流体计算软件Fluent对结冰过程进行数值模拟，该软件没有很明显的追踪固液界面，而是利用“焓-多孔度”理论来计算凝固过程，将流体和固体并存的糊状区域看作多孔介质来处理，流体凝固过程中，多孔度从1降低到0。Fluent软件采用有限体积法将非线性偏微分方程离散化为网格单元上的线性代数方程，然后通过求线性方程组得出流场的解，并用SIMPLE算法进行计算。

在Gambit中建模并对控制体进行网格的划分，网格划分好后，导入到Fluent软件中进行相关计算方法及参数的设置，并将自定义密度进行编译，即可开始计算。在计算中设置时间步长为0.1 s，每个时间步长的最大迭代次数为50次。

蓄冰桶内的水初始温度设置为288 K，为了验证数值模拟结果的正确性和真实性，对竖直铜管蓄冰进行了实验研究，铜管内通入低温(-5℃)的乙二醇溶液，使得铜管周围的水结冰，实验和数值模拟都进行2个小时。实验中通过T型热电偶和Fluke巡检仪进行温度的测量和记录。图4～9为纵向和径向温度分布。

图4 距管壁5mm纵向温度分布

图5 距管壁10mm纵向温度分布

图6 距管壁30mm纵向温度分布

图7 蓄冰桶上部径向温度分布

图8 蓄冰桶中部径向温度分布

从以上图可看出，数值计算出来的温度分布和实验测出的温度分布的发展趋势基本上是一致的，但是在实验中，加入蓄冰桶中的水的温度不均匀，且水在结冰过程中存在过冷度，因此使得实验数据和理论数据有一定的偏差。从径向温度分布图得出，点1，4，7实验和理论温度相差较大，而点2，5，8，6，9实验和理论值相差较小。观察到结冰的轮廓也是倒锥形，和数值计算得出的结冰轮廓一致。数值模拟得到的冰层厚度随时间的变化见图10，不同时间的流线变化见图11。

图9 蓄冰桶下部径向温度分布

图10 不同时间的结冰厚度

图11 不同时间的流线

刚开始冷却，靠近管壁的水由于重力的作用向下流动，而远离管壁的水向上流动，形成一个大的漩涡(见图11:120 s)。冷却进行到240 s时，靠近管壁的底角出现了相反的小漩涡。这是因为随着冷却的进行，水的温度不断降低，而温度首先降到4℃左右的水就出现在靠近管壁的底角处，由于水在4℃时密度最大，因此，使得温度高于4℃水向下流动时不能推动密度高的水流动，从而出现绕流。

随着时间的推移，靠近管壁的底角出现了新的小漩涡。这是因为，随着冷却的进行，底角处的水温温度降到低于4℃，其密度开始降低，受浮力的作用，这部分水开始向上流动。随着冷却的进行，水温的逐渐下降，密度最大的水层逐渐向上移动。而受到浮力作用的水逐渐向上流动，周围的水随之补充过来，形成了新的漩涡。随着水温的不断降低，漩涡逐渐扩大并且向上移动。直到1 200 s时，从图11可以看到，只有一个漩涡，而这个漩涡和开始时的漩涡的流动方向完全相反。刚开始漩涡是靠近管壁的水向下流动，而远离壁面的水则向上流动，形成的漩涡;而后面的漩涡是靠近管壁的水向上流动，而远离壁面的水则向下流动，形成的漩涡。

受到水的自然对流的影响，蓄冰桶内的水温出现了上下翻转。从距管壁10 mm处的纵向温度分布(见图5)，就可看出，刚开始从上到下的温度分布线为6点、5点和4点。而到后来，温度分布线从上到下变为4点、5点和6点。而从温度分布可看出，开始时温度一直均匀下降，当水温降到4℃左右，温度出现了突降现象。在冷却过程中，径向的温差很小，几乎没有温差(见图7，8，9)。

3 结论

由于水在4℃时有最大密度，Boussinesq线性假设不能真实的反映水的实际流动情况。因此采用了作为密度随时间的变化来考虑自然对流对蓄冰过程的影响。从结果中得出，纵向温度存在温差，径向各点温差很小。当温度达到4℃左右，出现温度突降和密度翻转现象。受水的密度翻转现象的影响，理论上计算得出的结冰的轮廓也不是圆柱状的，而是出现倒锥状，即上边结冰厚，下部结冰薄的现象。并且从数值模拟和实验测得的温度值对比可以看出数值模拟得出的温度分布与实验测量结果基本一致。

[1] SAEIDNH，SELANGOR，MALAYSIA.Maximum density effects on natural convetion in a porous cavity under thermal non-equilibrium condition[J].ActaMechanica，2007，188(1):55-68.

[2] MASAAKI ISHIKAWA，TETSUO HIRATA.Numerical simulation of natural convetion with density inversion in a square cavity[J].Numerical Heat Transfer.part A.2000，37(4):395-406.

[3] 刘道平，周文涛.伴随自然对流叠置圆管外蓄冰特性的试验研究[J].暖通空调，2002，32(4):14-16.

[4] 刘道平，周文涛.水平置恒壁温单圆管外蓄冰规律的实验研究[J].能源研究与信息，2005，21(1):37-42.

[5] 刘道平，周文涛.顺排管束外蓄冰过程的环水温度变化特性[J].上海理工大学学报，2003，25(3):205-208.

[6] MOTOIYAMAHA，NOBUONAKAHARA.Studies on thermal characteristics of ice thermal storage tank and amethodology for estimation of tank efficiency[J].Int T Energy Res，2008，32(3):226-241.

[7] 李恒，白博峰，陆军.圆柱腔体内水凝固过程热对流实验研究[J].工程热物理学报，2006，27(6):977-980.

[8] 白博峰，李恒，苏燕兵.矩形腔内纯水凝固过程试验研究[J].过程热物理学报，2007，28(3):435-438.

[9] 潘阳.热管解决土壤冻胀问题的数值分析及应用研究[J].华东交通大学学报，1998，15(4):44-50.

[10] NEMSTEELM W，MEDJANIK.Natural convection ofwater near density maximum in a rectangu larenclosure:Low Rayleighmumber calculation[J].Phys Fluids，1987，30(2):312-318.