不同管道内冰浆流动换热特性及冰堵分析

2021-08-18张雪刘圣春宋丽莹徐智明

制冷技术 2021年3期

张雪，刘圣春，宋丽莹，徐智明

（天津商业大学机械工程学院，天津市制冷技术重点实验室，天津 300134）

0 引言

随着全球电力需求的增长，许多环境问题如全球变暖、臭氧层损耗等受到了广泛关注，并严重影响全球能源政策[1]。为了缓解能源紧张，提高能源利用效率，储能技术发挥着关键作用。冰浆是由冰晶、液态水和防冻剂组成的絮状混合物，因其具有良好的流动特性和较大的能量密度，在同一运行工况下相较于传统冷冻水冷却能力高4 倍至5 倍[2-3]，已成为一种新型相变材料在蓄冷系统中被广泛使用[4-8]。

目前，国内外学者对冰浆系统的研究较多。在换热特性方面，JEAN-PIERRE 等[9]实验研究了双管换热器中冰浆的换热情况，结果表明，随着含冰率的减小，传热系数先减小，在低含冰率时基本保持不变。赵腾磊等[10]以含冰率在20%以下的冰浆为模型进行实验研究，分析冰浆的密度、导热系数、比热和动力黏度，结果表明冰浆的传热系数与含冰率成正比。张曼等[11]建立了动态冰浆传热特性的数理模型，数值研究了流速、管径和冰浆体积分数对传热系数的影响。LEE 等[12]实验研究了以冰浆作为冷却介质的双管换热器，结果表明换热率随质量流量和含冰率的增加而增加，流速低时，冰粒对换热的影响更为明显。BELLAS 等[13]实验研究了带冰浆的板式换热器的性能，得出当含冰率在5%～20%时，传热系数基本保持不变。LIU 等[14]数值研究了水平直管直径、长度、含冰率和流速对冰浆流动压降的影响，验证了仿真模型的可靠性，进一步发展了黏度与流速的数学模型。简夕忠等[15]数值研究了含冰率、流速、管径和添加剂种类等因素对冰浆换热效果的影响。由于冰浆系统涉及到相对比较复杂的流动及传热问题，目前，对冰浆传热特性的研究在理论分析、数值计算和推广应用方面仍需完善。

本文利用颗粒动力学理论，通过数值模拟的方法得出冰浆在不同类型管道内的流动换热特性，并对直管内冰堵现象进行了分析。

1 数学模型

在Fluent 软件中采用双欧拉流体模型，结合颗粒动力学理论建立数值模型的控制方程。

1.1 质量守恒方程

冰浆本质上为一种液固两相混合流体，冰粒固体为直径较小的颗粒，内置颗粒相模型，设两相分别为q相和p相。

相位q的体积：

q相的有效密度按式(3)计算：

q相的连续性方程为：

由质量守恒方程，得：

式中，αq为q相的体积分数；ρq为q相的密度，kg/m3；mpq为p相到q相的质量传递速率，kg/(m3·s)；vq为q相的速度，m/s。

1.2 动量守恒方程

将冰浆看作固液两相流体，用液相动力参数描述动量方程[16]，可得：

式中，下标s 表示固相；l表示液相；vsl为固液两相间的相对速度，若msl＞0，则vsl=vs；若msl＜0，则vsl=vl，且vsl=vls；τl为液相切应力，Pa；Fl为外部体积力，N；Flift,l为升力，N；FVm,l为虚拟质量力，N；Ksl为固液相间的动量交换系数。

固体颗粒相动量方程为：

液相与固体颗粒相之间的作用力可采用GIDASPOW[17]模型确定。

当αl＞0.8 时，固-液动量交换系数Ksl为：

当αl≤0.8 时，Ksl有如下形式：

式中，Re为相对雷诺数；ds为颗粒直径，m。

q相和p相的相对雷诺数为：

1.3 能量方程

能量控制方程[18]如下：

式中，λeff为冰浆的有效导热系数，W/(m·K)；ST为源项，W/m3。

1.4 湍流方程

选择k-ε湍流模型，考虑固液两相间的相互作用[16]：

式中，C1ε=1.44；C2ε=1.92；σε=1.3；σk=1；μt,m=ρmCμk2/ε；Cμ=0.09；下标m为固液混合相；k为湍流脉动动能，J；ε为湍流脉动动能耗散率；μt,m为混合相湍流动力黏度，kg/(m·s)；σk为湍流脉动动能普朗特数。

1.5 黏度模型

冰浆是一种固液两相混合物，物理密度、动力黏度以及结冰点不仅与受添加物体积分数影响，还与含冰率有关，THOMAS[19]通过反复实验和理论推算，得出黏度关于含冰率的公式：

式中，φi为含冰率，%。

1.6 相变模型

在冰浆换热特性模拟过程中，开启能量方程，将冰-水换热模拟写入相变UDF 程序中，收敛精度设为10-4。

2 数值模拟模型图及边界条件设置

2.1 管道模型

90°弯管和T 型管模型分别如图1(a)和图1(b)所示。管道直径为0.02 m，管总长为1 m，弯管和T 型管均在中间位置处发生改变，弯管曲率半径为0.1 m，网格质量在0.7 以上，管道分流处采用结构化网格形式。溶液进出口位置及重力方向见图1。

图1 90°弯管和T 型管管道网格划分模型

2.2 边界条件设置

对冰浆在90°弯管、T 型管内的换热特性进行模拟。具体参数设置为：含冰率IPF 为5%～25%，流速为0.2～1.5 m/s。冰粒直径约为10-4m。边界条件设为速度入口，出口设为出流，壁面设为无滑移，采用定壁温条件，温度为1 ℃；入口处冰相和水相的温度均为0 ℃；湍流模型选用k-ε模型。冰浆两相流体的物性参数如表1所示。

表1 冰浆两相流体的物性参数[20]

3 数值模拟结果分析

3.1 90°弯管内冰浆换热特性结果分析

图2 和图3所示为含冰率为10%时，90°弯管截面上溶液平均温度分布和焓值分布。由图2 和图3 可知，纵截面上，溶液平均温度变化较小，保持在0～1 ℃，弯管处溶液焓值最大为2.46×104J/kg。

图2 90°弯管截面上溶液平均温度分布

图3 90°弯管截面上溶液的焓值分布

图4所示为90°弯管截面上流体流动的速度分布。由图4 可知，在90°弯管中，流体流动速度从入口处逐渐增大直至弯管处达到最大，约为1 m/s，经过转弯处，出现明显分层现象，在离心力的作用下，冰粒靠近管外侧速度明显大于内侧，逐渐在外侧堆积，这与文献[21]的结果相符。

图4 90°弯管截面上流动速度分布

综合以上分析，在90°弯管处，由于存在摩擦和二次流现象[22]，冰浆流体流动阻力大，扰动强烈，摩擦损失加大。而流速从入口处逐渐增大，在弯曲处流速最大且开始出现分层现象，在流速较大区域，冰-水换热效果缓慢，冰粒不易完全融化[23]。因此，溶液平均温度随流速增大而升高，在弯曲处温度高于入口处，且冰的密度略小于水，由于受到浮升力作用和弯管处冰浆流速较快的影响，冰粒在流体的带动下，逐渐沿着管壁一侧分散开来。经过弯曲处，平均流速减小，温度稍降低。溶液温度分布从大到小依次为弯曲处、出口处和入口处，但总体温度比较均匀。

3.2 T 型管内冰浆换热特性结果分析

含冰率为10%时，T 型管纵截面处冰粒焓值分布，如图5所示。

图5 T 型管纵截面冰粒焓值分布

由图5 可知，入口处溶液的焓值较高，随着管道流动呈梯形下降趋势，在分流处焓值减小。这是由于冰浆在管内流动中，摩擦损失和紊流效应使能量损失。在交叉处，扰动强烈，液态水与冰粒之间交换动量。分流之后固液两相流动趋于稳定，在出口处压力再次降低。

3.3 直管内冰浆流动冰堵分析

为了解决冰浆堵塞造成换热效果差这一问题，对直管中冰浆的分布进行数值模拟。管道直径为0.02 m；入口边界条件为速度型，出口为出流型。重力方向为Y轴负方向。当含冰率为20%，取管道中间横截面位置进行分析，得到冰粒的体积分数分布结果如图6所示。可以分为低体积分数区、高体积分数区和过渡区。图7所示为文献[24]中的3 种冰浆的流动流态。

图6 不同流速下冰粒在直管横截面上体积分数分布

图7 文献[24]中3 种冰浆流动流态

由图6(a)和图6(b)可知，当流速为0.2 m/s 和0.3 m/s 时，冰粒体积分布分层明显，上层冰粒体积分数最高约为0.4，底部冰粒体积分数最低，几乎为0，中间为过渡区，这与图7(a)对应，由于流速较小，在高体积分数区黏滞力较大，使冰层沉积，管道堵塞。过渡区流速较高、体积分数区大，小于临界沉积速度，仍有可能造成冰堵。当流速为1 m/s时，由图6(c)可知，冰粒体积分数分布整体较为均匀，上层为移动床流动，下层为非均质流动，即悬浮床流动，与图7(b)对应。当流速大于1.5 m/s 后，如图6(d)和图6(e)所示，此时，流速大于临界沉积速度，整个管道内呈现为悬浮床模型，冰粒速度相当，体积分数分布均匀，与图7(c)相对应，发生冰堵情况较少。

通过对不同流速下的冰粒体积分数模拟结果的比较，选择1.5 m/s 作为冰浆由移动床流动变为悬浮床流动的第二临界流速，冰粒分布趋于最均匀。

不同含冰率条件下，冰浆体积分数的最值与流速的关系如图8所示。由图8 可知，随着含冰率的增加，A 点的第一临界速度逐渐减小，在1～2 m/s之间，C 点的第二临界速度减小，在10 m/s 左右。其原因是含冰率增加对湍流效应的影响大于对黏性摩擦的影响，因此，含冰率高的冰浆更有可能在较低流速下转化为均匀流，减少冰堵发生。根据固体冰的分布轮廓，冰浆的总体积为：