基于群体平衡模型的冰浆流动与传热特性数值研究

邹良旭，马 非，孟昭男，张 鹏

(上海交通大学 制冷与低温工程研究所，上海 200240)

由于全球化石能源的不断消耗以及温室气体的排放使得各种环境问题日益凸显，人们愈加重视节能减排和提高能源利用率.蓄能技术可以解决能量供需在时间和空间上不均衡的问题，可以有效地提高能源利用率.冰浆是一种固-液两相可流动的相变材料，其固体颗粒直径通常不大于1.0 mm[1]，并具有较好的储能、流动和传热特性，广泛应用于建筑空调、区域供冷和食品存储等领域[2-4].在实际应用中，冰浆固相颗粒直径在流动和传热过程中不断变化，且冰浆的流动特性和固液两相之间的传热传质过程十分复杂，因此深入理解冰浆的流动与传热特性对其应用十分重要.

目前所进行的研究大多忽略了冰浆固相颗粒直径在流动与传热过程中的变化，采用数值计算获得冰浆的流动与传热特性.Wang等[5]通过对比实验与数值模拟结果发现，使用基于固相颗粒动力学理论的 Euler-Euler方法可以更准确地模拟出冰浆的流动过程，以及弯管处的阻力对冰浆流动产生的影响.Kumano等[6]研究了冰浆存储时间对冰浆流动特性的影响，作者利用显微镜对冰浆颗粒进行观察，发现流动工况下其直径会随时间的增长而变大，当储存时间超过12小时，圆管内因冰颗粒的聚集而出现堵塞现象.

传统的两相流数值模拟计算模型忽略了冰浆固相颗粒直径的变化.群体平衡模型是一种针对颗粒直径变化的数值计算模型，可较为准确地描述冰浆内固相颗粒的破碎、融化和聚并等过程.Xu等[7]采用基于Euler-Euler方法的计算流体动力学-群体平衡模型(CFD-PBM)对冰浆在水平直管内的流动进行了数值计算，计算获得的压降与实验数据较为吻合，验证了CFD-PBM模型的可靠性及准确性，但在该文献中尚未涉及冰浆的传热特性.Jin等[8]利用两相流模型和群体平衡模型对氮浆在水平管道中的流动与传热特性进行了模拟计算，研究结果表明加入群体平衡模型后的计算模型更准确，在不同工况下压降和传热系数的计算结果与实验结果的偏差在5%以内，但该文献未对颗粒直径变化对流动与传热特性的影响进行讨论.

针对上述问题，本文采用基于Euler-Euler法的CFD-PBM模型，对冰浆在水平圆管中的流动与传热特性及其固相颗粒的直径变化进行研究.在流动特性的研究中，分析固相体积分数与流速对颗粒直径在流动过程中的影响，为传热特性的研究提供依据；在传热特性的研究中，分析不同热流密度对颗粒直径变化的影响.同时在上述两类工况中，将CFD模型与CFD-PBM模型的计算结果进行对比，因CFD-PBM模型考虑了冰浆颗粒的直径变化，使得CFD-PBM模型对冰浆流动与传热过程的数值模拟更接近实际情况.

1 CFD模型

采用基于Euler-Euler法的CFD模型对冰浆在水平圆管中的流动与传热特性进行模拟计算.假设固-液两相均为连续性流体，且均具有质量、动量和热量的交换；冰浆在管内流动为不可压缩流体湍流动，固相为表面光滑、无弹性的球体[9].采用标准k-ε湍流模型描述冰浆在水平圆管中的湍流过程[10].因冰浆的两相密度差异较小，故相间作用力只考虑曳力、升力和湍流扩散力，颗粒动力学模型详见文献[9].

1.1 连续性方程

(1)

(2)

1.2 动量守恒方程

固相动量方程如下：

(3)

式中：p为局部压力；g为重力加速度；F为两相之间的总作用力；ps为固相压力；τs为固相应力张量.

液相动量方程如下：

(4)

式中：τl为液相应力张量.

1.3 能量守恒方程

在传热的工况下需要考虑能量守恒方程，基于Euler-Euler模型的假设，冰浆中的固-液两相都遵循能量守恒方程：

(5)

(6)

式中：H为焓值；λeff为等效热导率；hV为单位体积内固-液两相间的传热系数；T为温度.

2 群体平衡模型

2.1 模型描述

在冰浆的固-液两相流中，固相颗粒在实际物理过程中存在破碎、融化和聚并等现象，这些现象会导致固相颗粒直径分布和尺寸大小发生变化.群体平衡模型可用来描述固体颗粒在各种现象下的分布状态，使用正交矩方法对群体平衡方程进行求解并计算其Sauter平均直径，即固相颗粒直径[11-12].关于固相颗粒的群体平衡方程如下：

(7)

聚并是多个颗粒之间的团聚现象.聚并会导致颗粒尺寸增大，其颗粒聚并频率可用布朗运动和颗粒碰撞作用下的模型描述[13]：

(8)

式中：kB为Boltzmann常数；μ为液体黏度；ε为湍流扩散率；η为动力黏度.

破碎的作用与聚并相反，它表示颗粒之间由于碰撞导致的颗粒尺寸减小的过程.破碎现象会受到冰浆固相颗粒直径尺寸的影响，颗粒破碎频率由下式表示[7,14]：

δ(V)=c1ηc2εc3Vc4

其中：c1～c4均为系数，其值大小分别为 0.000 6，-1.25，0.75，1.在破碎模型中，固相颗粒破碎的概率密度函数描述了由破碎导致颗粒的体积由V′变为V的概率，并且可以描述颗粒数量以及由破碎导致的尺寸分布.本文的固相颗粒为球状颗粒，不同颗粒的大小不同，因此采用抛物型概率密度函数，其表达式为

β(V|V′)=

(9)

2.2 CFD与PBM模型的耦合过程

CFD-PBM模型首先利用CFD两相流模型，通过求解质量守恒、动量守恒、能量守恒等方程来获得固-液两相的流速、温度以及固相体积分数；再采用PBM模型通过正交矩方法求解固相颗粒的聚并、破碎以及融化过程，并且得到颗粒的尺寸及其分布；最后通过获得的固相颗粒直径，重新计算固相的流速、温度以及体积分数，如此循环以实现两种模型的耦合[15].

3 数值计算

考虑两类计算工况，一种仅计算冰浆在非均质状态下的等温流动工况，称为流动工况；另一种计算管道壁面有热流时冰浆的非均质、非等温流动的情况，称为传热工况.上述两类工况均采用质量分数为10.3%的乙醇-水溶液作为冰浆的载流体.在数值计算中，对于壁面边界条件，液相采用无滑移壁面条件，固相采用Johnson-Jackson壁面滑移条件[9].在传热工况中，还需在管壁添加热流密度边界，假定圆管入口处的固-液两相为相同流速和温度.

在流动工况中，圆管的长度为 1.3 m，直径为23 mm.工质的物性考虑为常数[9].固体颗粒的初始直径为125 μm.在传热工况中，为了观察到明显的冰浆融化过程，采用长度为 2.0 m的圆管，直径为16 mm，固-液两相的入口流速为1.0 m/s，固相体积分数为10%.在传热工况下，工质的物性随温度发生变化[16]，固相颗粒的初始直径设置为270 μm.除了考虑工况，还需对比CFD模型与CFD-PMB模型在相同工况下的结果，数值模拟计算的工况及条件如表1所示.其中：CFD表示不加PBM模型的模拟计算；CFD-PBM表示加入PBM模型的模拟计算；工况名后的数字为编号.

利用FLUENT 14.5软件进行数值模拟计算，计算时间步长为0.01 s.圆管网格划分采用六面体结构化网格，在近壁面处加密处理，网格均进行网格无关性分析以保证计算的准确性.流动工况下的网 格数量为 1 347 312，传热工况下的网格数量为2 739 380.在数值计算过程中，所有的控制方程均采用一阶迎风离散格式进行计算.每个时间步长的最大迭代次数为50次，当每一步迭代计算的速度项残差低于1.0×10-4、能量项残差低于1.0×10-6时，视为计算收敛并进行下一步计算.

表1 不同条件下的流动与传热工况

Tab.1 Flow and heat transfer of ice slurry under different conditions

工况入口流速/(m·s-1)固相体积分数/%壁面热流密度/(kW·m-2)CFD11.010-CFD-PBM11.010-CFD-PBM22.010-CFD-PBM31.015-CFD21.01050CFD-PBM41.01015CFD-PBM51.01025CFD-PBM61.01050

4 模型验证及数值计算结果

4.1 模型验证

CFD-PBM模型对冰浆流动的数值模拟计算的准确性验证见文献[7].由于冰浆固相颗粒直径的测量在实验中难以捕捉测量，所以采用压降作为模型准确性的衡量手段.在其实验中，冰浆载流体为质量分数6.5%的乙二醇-水溶液，在不同流动工况下，实验得到的冰浆在水平圆管中流动的压降与CFD-PBM模型的计算结果最大偏差为10.32%，验证了CFD-PBM计算冰浆流动工况的准确性.

4.2 CFD-PBM流动工况结果

在流动工况中不涉及到热量的传递，冰浆中固相颗粒之间只有聚并和破碎作用，不存在融化作用.固相颗粒直径D沿管道轴线的分布如图1所示，其中Lpip为管道轴线距离.由图1可知，CFD模型固相颗粒直径保持不变，而CFD-PBM模型的固相颗粒直径沿流动方向逐渐增大.分析CFD-PBM模型的计算结果可知，在Lpip=0～0.1 m区间内，固相颗粒的聚并和破碎效果很弱，因此颗粒直径几乎没有任何变化；Lpip=0.1～0.6 m区间是颗粒直径增长的发展段，聚并作用比破碎作用强，颗粒直径开始增长；Lpip=0.6～1.3 m区间颗粒直径以稳定的增长速率逐渐增大.在CFD-PBM1、CFD-PBM2、CFD-PBM3工况中，沿管道流动方向，颗粒直径分别从初始状态的125.0 μm增长到127.6、127.8和128.9 μm，说明增大固相体积分数及流速均可以促进颗粒直径的增长.

图1 颗粒直径沿管道轴线的分布Fig.1 Mean diameter distributions along central axial direction

各工况下管道出口截面对称轴上的颗粒直径及CFD-PBM1工况出口截面的颗粒分布如图2所示，其中R为对称轴上任意位置到轴线的距离与半径之比.由图2可知，CFD1模型中的固体颗粒直径无任何变化，而CFD-PBM模型中近壁面处的颗粒直径比主流区域的颗粒直径大.这是由于近壁面处的湍流扩散率较大，导致了聚并效果增强，最终使近壁面处的颗粒直径比主流区的颗粒直径大.对比图2(a)中CFD-PBM1工况和CFD-PBM2工况的结果发现，入口流速越大，颗粒直径就越较大，这是由于较大的流速有利于颗粒的碰撞，从而强化了颗粒的聚并作用.同时，流速越大，冰浆颗粒的分布越均匀，进而导致管内顶部与底部颗粒直径大小差异减小.通过上述分析可知，管道内颗粒的最大直径出现在出口顶部，在CFD-PBM1、CFD-PBM2和CFD-PBM33种工况的最大颗粒直径分别为139，142，146 μm.

图2 管道出口截面颗粒直径的分布Fig.2 Mean diameter distributions of ice slurry on the cross-section at the pipe outlet

图3 CFD-PBM1工况下出口截面固相体积分数分布Fig.3 Solid volume fraction distributions of outlet for CFD-PBM1

4.3 CFD-PBM传热工况结果

由于传热工况中热量的传递会导致冰浆固相颗粒融化，而融化会使冰浆固相体积分数降低、颗粒直径减小，因此冰浆中固体颗粒之间同时具有聚并、破碎和融化作用.不同工况下冰浆沿管道轴线的固相体积分数如图4所示.在冰浆的传热中，由于管壁的热量先传递到液相，使液相温度升高，随后液相将热量再传递给固相，使固相颗粒随着流动不断融化，固相体积分数不断降低.在CFD2、CFD-PBM4、CFD-PBM5及CFD-PBM6工况中，固相体积分数分别从初始状态的10%下降至出口截面处的1.21%、6.37%、4.33%及0.99%.

图4 沿流动方向轴线的固相体积分数变化Fig.4 Solid volume fraction along the axis of the flow direction

不同工况下，冰浆沿流动方向轴线的固相颗粒直径变化如图5所示.CFD模型下固相颗粒的直径保持不变，而CFD-PBM模型的固相颗粒直径沿着流动方向的轴线距离逐渐减小.在CFD-PBM4、CFD-PBM5和CFD-PBM6工况中，沿管道流动的轴线方向，颗粒直径分别从初始状态的270 μm融化至255、249和154 μm.结合图4和图5可知，在热流密度较大的条件下，固相颗粒直径减小的速度更快.颗粒直径的减小有利于冰浆的融化，使得冰浆固相体积分数减小得更多，因此随着颗粒直径的减小，CFD-PBM6工况的固相体积分数比同条件下CFD2工况的固相体积分数略低.

图5 沿流动方向轴线的固相颗粒直径变化Fig.5 Mean diameter of ice slurry along the axis of the flow direction

CFD-PBM6工况下的固相颗粒平均直径及其分布如图6所示.由图6可知，颗粒直径沿管道流动方向逐渐减小.对于管道不同位置的截面对称轴，颗粒直径分布如图7所示.由图7可知，由于固相体积分数与颗粒直径的相互影响，管道顶部固相颗粒直径比底部大.随着沿程距离的增大，管道近壁面处的颗粒完全融化，管道顶部与底部颗粒直径的差异逐渐减小，故出口附近的颗粒直径基本呈对称分布.

综上所述，在传热工况下，CFD模型并未考虑管道中的固相颗粒直径的变化及其分布规律，与实际情况不符.CFD-PBM模型同时考虑了颗粒的聚并、破碎和融化作用(以融化模型为主导作用)，最终导致了颗粒直径近壁面处的直径小于主流区，沿流动方向中轴线的固相颗粒直径逐渐减小.因此，CFD-PBM模型比CFD模型更能合理地反映出冰浆的实际传热特征与变化规律.

图6 CFD-PBM6工况下颗粒直径在管道轴向截面的分布Fig.6 Mean diameter distributions on the vertical symmetrical cross-section of the pipe for CFD-PBM6

图7 CFD-PBM6工况下颗粒直径在不同位置的截面对称轴上的分布Fig.7 Mean diameter distributions at different cross-sections along the flow direction for CFD-PBM6

5 结论

本文采用基于Euler-Euler法的CFD-PBM模型，对冰浆在水平圆管中的流动与传热特性及其固相颗粒直径的变化规律进行了数值模拟计算.通过对比不同条件下的CFD模型与CFD-PBM模型的计算结果可以获得以下结论：

(1) 使用CFD-PBM模型能够更准确地反映出冰浆在水平圆管中的流动与传热特性.研究发现，冰浆的固相体积分数与颗粒直径大小呈正相关，高固相体积分数区域的颗粒直径比低固相体积分数区域的颗粒直径大，同时颗粒直径的增长也导致了固相体积分数的增加.而流速越快会增加聚并作用，最终加速颗粒直径的增长.

(2) 在流动工况中，随着在圆管中的流动，固相颗粒直径逐渐增大.在入口流速为1.0 m/s与固相体积分数为10%的情况下，固相颗粒的直径最大可从初始状态的125 μm增长到139 μm.其中，管道顶部的颗粒直径大于底部，近壁面处的颗粒直径大于主流区.

(3) 在传热工况中，固相颗粒直径随着固相的融化不断减小.在入口流速为1.0 m/s、固相体积分数为10%与壁面热流密度为50 kW/m2的工况下，固相颗粒直径从初始状态的270 μm直至完全融化.其中，管道顶部的颗粒直径大于底部，近壁面处的颗粒直径小于主流区.出口处附近颗粒直径呈对称分布.热流密度的增大会加速固相的融化，进而加快颗粒直径的减小.