冰浆存储过程冰晶粒径演化数值模拟
赵腾磊，刘志强，徐爱祥，何娜萍
(中南大学 能源科学与工程学院，湖南 长沙，410083)
摘要：为了探索冰浆存储过程冰晶粒径演化规律，运用数群平衡模型对冰浆存储过程冰晶粒径演化进行数值模拟，采用改进有限容积法对数群平衡方程进行离散求解，把数值结果与实验结果进行对比验证模型的可靠性，并研究冰晶团聚、破碎2个动力学事件对冰晶粒径分布演化的影响。研究结果表明：冰晶粒径演化数值模拟结果与实验结果的拟合曲线基本吻合，因此数群平衡方程模拟冰晶粒径演化可靠。对破碎与团聚作用下冰浆存储过程中冰晶粒径演化进行模拟，团聚作用使得大粒径冰晶颗粒增多，小粒径冰晶颗粒减少；破碎作用使得小粒径冰晶颗粒增多，大粒径冰晶颗粒减少；团聚和破碎共同作用下，团聚比破碎对冰晶粒径演化的影响明显，使得大粒径冰晶颗粒增多，小粒径冰晶颗粒减少。
关键词：储能；冰浆；粒径演化；数群平衡模型；破碎；团聚
中图分类号：TB61.1 文献标志码：A 文章编号：1672-7207(2014)10-3651-06
Numerical simulation of evolution of ice crystal size distribution during storage
ZHAO Tenglei, LIU Zhiqiang, XU Aixiang, HE Naping
(School of Energy Science and Engineering, Central South University, Changsha 410083, China)
Abstract: In order to study the principle of the evolution of ice crystal size distribution during the process of ice slurry storage, population balance model was introduced to make a numerical simulation on the evolution of ice crystal size distribution. The improved finite volume method is used to solve the population balance equation. Population balance model’s reliability was tested based on the comparison of numerical results and experiment results. Then it is time to research the effect of breakage and aggregation on the evolution of ice crystal size distribution. The results show that the changing curve for the evolution of ice crystal size distribution based on numerical results and experiment results is generally the same. So it is reliable to use population balance model to imitate the evolution of ice crystal size distribution. Aggregation makes the big size ice crystal increase and the small size ice crystal decrease. Breakage makes the small size ice crystal increase and the big size ice crystal decrease. When breakage and aggregation happen at the same time, the effect of aggregation on the ice crystal size distribution is more striking than that of breakage.
Key words: energy storage; ice slurry; particle size evolution; population balance model; breakage; aggregation
冰浆是由载流溶液和冰晶组成的固液两相流[1]，具有良好的流动与传热特性、高潜热、无污染等优势，因此受到化工、环保、燃烧等多个领域的广泛关注[2-3]。冰浆存储过程冰晶会受到成核、生长、破碎、团聚等动力学事件的作用[4-6]，这些动力学事件会使冰晶粒径分布随着时间发生变化，冰晶粒径演化会对冰浆的流动与传热特性产生重要影响[7]。国内外许多学者现今还主要集中于从宏观角度对冰浆存储进行研究，比如对冰浆存储过程中的富集与分层现象进行研究[8-9]、在富冰层模型基础上运用CFD模拟研究了冰浆存储过程中冰晶粒子的运动及变化轨迹[10]等，但是这些研究无法对冰浆存储过程冰晶动力学演化进行阐述。Pronk等[11]从微观角度对冰晶团聚、磨损以及熟化效应进行了研究，但只是对动力学事件的作用机理进行了阐述，并没有探讨动力学事件对冰浆粒径演化的影响。以上研究表明冰浆存储过程冰晶演化研究还较少，特别是冰晶动力学事件对冰晶粒径演化的影响还有待深入研究。数群平衡模型(population balance model, PBM)是模拟离散系统中离散单元动力学演化过程的有效方 法[5]。数群平衡模型在离散系统中离散相演化研究过程的应用对冰浆存储过程冰晶颗粒动力学演化的研究带来了启示。冰浆流体系统是由载流溶液和冰晶颗粒组成的离散系统，冰晶的成核、生长、团聚与破碎等动力学作用都可以运用数群平衡方程(population balance equation, PBE)进行表达。徐爱祥等[12]运用数群平衡模型对团聚和破碎作用下冰晶粒径演化规律进行了研究，但在研究过程中是采用分区法对群体平衡方程进行直接离散求解，而数群平衡方程中团聚项和破碎项是基于冰晶数量变化所形成的数学表达式，只能保证数量守恒而忽略了质量守恒。本文作者运用数群平衡模型对冰浆存储过程冰晶粒径尺度演化进行模拟，把方程中团聚项和破碎项变换成质量通量的形式进行处理，有效地避免了文献[12]基于数量守恒进行计算而引起质量损失过大的问题，其次本文直接把数值结果与实验结果进行对比，改善了文献[12]把实验结果变换成概率分布后与数值结果进行对比而带来的误差，最后着重讨论了破碎作用与团聚作用对冰晶粒径演化的影响。
1 数群平衡模型
1.1 数群平衡方程
Perry等[13]给出了离散系统中离散单元零维数群平衡方程：
(1)
式中：f(t, x)为离散颗粒尺度分布函数；G(t, x)为生长速率，表示t时刻单位时间尺度为x的离散单元尺度增长的大小[14]；Aagg(t, x)和Bbreak(t, x)分别为团聚项和破碎项；S(t, x)为成核项。
团聚项定义为
(2)
方程右边第1项表示尺度为y和x-y的离散单元聚合为尺度x的离散单元数量；第2项表示尺度为x的离散单元和其他离散单元团聚而减少的数量。β(t, x, y)在统计学上表示离散单元单位时间的团聚概率(计算区域内离散单元团聚速率)。
破碎项定义为
(3)
方程右边第1项表示尺度为y的离散单元破碎为尺度x的离散单元数量；第2项表示尺度为x的离散单元因破碎而减少的数量。p(t, x, y)表示单位时间内尺度为y的离散单元破碎为尺度x的离散单元概率，是与破碎离散单元尺度有关的函数；b(x) 在统计学上表示单位时间尺度为x离散单元破碎的概率(计算区域内离散单元破碎速率)。
令：
(4)
(5)
将方程(1)左右两边乘以x，且把方程(4)和(5)代入方程(1)，变换方程为质量守恒方程：
(6)
其中：
(7)
(8)
运用有限容积法对方程(6)进行离散：
(9)
其中：
(10)
(11)
(12)
1.2 计算区域离散
为了运用数值方法对数群平衡方程进行数值计算，首先要对计算区域进行离散。区域的离散方法分为均匀离散和非均匀离散2种。本文运用均匀离散对计算区域进行离散。
令N为计算区域总数，用xi-1/2表示区域[xmin, xmax]中的一个离散区域节点，xi表示第i个离散区域，xmin和xmax分别表示离散单元内部变量的最小值和最大值，i∈(1，…，N)。
对变量进行以下设定：
(13)
令计算区域为1~1 401 μm，减小时间步长与空间步长可以减小计算误差，但是结合计算时效性，选取离散区域总数N为100，离散时间步长为1 s。
2 冰晶粒径演化数学模型
2.1 模拟条件
为了运用数群平衡模型对冰浆存储过程冰晶粒径演化进行数值计算，进行如下假设：1) 不考虑冰晶的空间分布，认为冰晶颗粒在空间上具有同一性；2) 冰晶在存储过程中成核与生长已不明显，忽略冰晶成核与生长；3) 忽略载流体流动对冰晶演化的影响。因此，冰浆存储过程中冰晶粒径演化的零维空间数群平衡模型离散方程为
(14)
现今还没有运用数群平衡模型研究冰浆存储过程冰晶动力学演化，缺少冰晶颗粒动力学参数模型公式，因此，针对冰晶特征、实验条件以及假设的模拟条件选取其他与冰浆相似的多相流体系统离散相颗粒动力学参数模型公式作为研究冰晶动力学演化的动力学参数模型公式。
假设冰浆存储过程中冰晶粒子在惯性区域受流体扰动作用下保持各向同性，不同冰晶粒子碰撞时破碎速率[15]的表达如下：
(15)
式中：C1为常数，其值为0.6×10-3；ν为载流体运动黏度，m2/s；ε为耗损率，表示单位质量的冰浆在外界扰动或内部粒子间相互作用下所耗散的功率，W/kg。
粒径为x′的冰晶破碎为粒径为x的冰晶概率函数[16]如下：
(16)
式中：C为形状因子，本文假设冰晶为球形，这时形状因子C=1分布呈凹抛物线形。
描述载流体中固相冰晶粒子的团聚速率[16]表达式如下：
(17)
式中：μ为冰浆中载流溶液的动力黏度，Pa·s；KB为Boltzman常数，是有关于温度与能量的物性常数，J/K；T为冰浆存储时的平衡温度，K；μ为冰浆中载流溶液的动力黏度，Pa·s。
根据实验条件对模型物性参数进行设置，实验过程中是用浓度为4.9%的氯化钠水溶液进行过冷形成的冰浆溶液，冰浆溶液物性参数如表1[17]所示。
表1 冰浆溶液物性参数[17]
Table 1 Physical parameters of ice slurry solutions[17]
边界条件：冰浆存储过程中冰晶碰撞形成的最小颗粒范围是1~10 μm[18]，因此计算过程中取最小粒径冰晶颗粒为1 μm，1~1 401 μm区域的冰晶粒子数量占总粒子数量的96%以上，因此在模拟冰晶粒径分布演化的过程中取计算区域为1~1 401 μm，舍去计算区域以外的冰晶颗粒，也就是说假设在计算区域以外的冰晶粒径分布数量密度为零。
初始条件：对实验测定数据进行拟合，用数学方法对冰浆存储过程初始时刻实验数据拟合曲线进行处理，把实验过程中冰晶初始粒径分布用数学表达式进行表示如下：
(18)
其中：x为取值范围，1 μm＜x＜1 401 μm。
2.2 模型验证
为了验证数群平衡方程求解方法的有效性以及冰晶数群平衡模型的可靠性，将冰晶粒径数量密度分布曲线演化的模拟结果与Pronk等[11]的实验结果进行对比。Pronk等[11]较早对冰浆存储过程冰晶动力学演化进行了实验研究，实验运用双层器壁的搅拌结晶器作为冰浆制取以及存储简易实验装置，提出运用电子显微镜法观测冰晶演化过程。Pronk等[11]进行了多组冰浆存储实验，对冰晶尺度分布数量密度、尺度、形状的演化过程进行研究，因此实验结果可以作为本文模拟结果对比验证的对象。图1和图2所示分别为冰浆存储2 h和6 h后冰晶粒径数量密度分布演化曲线，如图1和图2所示实验拟合冰晶粒径分布演化曲线与数值模拟曲线基本吻合。分析图1和图2中冰晶粒径分布曲线变化可知：随着冰浆存储时间的增加，冰晶粒径分布演化实验拟合曲线与数值模拟曲线误差越来越小，也就是两条曲线重合度越来越高。从图1和图2中的曲线变化分析可知本文数群平衡方程求解方法有效且冰晶数群平衡模型模拟冰晶粒径分布演化准确可靠。实验结果与模拟结果存在着一定的偏差，主要原因是：1) 冰浆存储实验过程中所取样本冰晶颗粒形状不规则，因此实验数据是通过实验处理工具得到。本文在模拟冰晶粒径分布演化过程中是把冰晶作为圆形处理。2) 冰浆存储过程中冰晶粒径分布演化是成核、生长、团聚和破碎4个动力学事件共同作用的结果，而本文假设冰晶只受到团聚、破碎的作用。3) 在数值计算时，忽略载流体的流动及传热、传质过程对冰晶粒子动力学演化的影响。
图1 冰浆存储2 h后冰晶粒径数量密度分布的实验拟合曲线与数值模拟曲线
Fig. 1 Curve of ice crystal size distributions after 2 h storage obtained from model and experiment
图2 冰浆存储6 h后冰晶粒径数量密度分布的实验拟合曲线与数值模拟曲线
Fig. 2 Curve of ice crystal size distributions after 6 h storage obtained from model and experiment
3 模拟结果与分析
3.1 纯破碎作用
首先探讨破碎作用对冰浆存储过程冰晶粒径演化的影响。图3所示为冰晶只发生破碎作用时冰晶粒径分布随时间变化的曲线。由图3可知：冰晶粒径分布曲线随着时间的增加有不断向左移动的趋势，也就是说小粒径冰晶颗粒数量不断增加，大粒径冰晶颗粒数量不断减少。根据冰晶破碎机理可知大粒径冰晶颗粒不断破碎为多个小粒径冰晶颗粒，因此冰晶颗粒总数增加且冰晶颗粒主要集中在小粒径范围内，所以引起了图3中曲线向左移动的同时峰值也有不断增大的趋势。图3中冰晶粒径分布曲线变化幅度不大，且随着存储时间的推移，冰晶粒径分布曲线的变化幅度会越来越小。分析可知冰浆存储过程中破碎作用主要是使冰晶颗粒数量不断增加，随着时间推移冰晶颗粒数量浓度不断增加，当冰晶数量浓度增加到一定程度后，破碎作用效果就被削弱。
图3 冰晶在破碎作用下的粒径数量密度分布演化曲线
Fig. 3 Curve of ice crystal size distributions in time under breakage
3.2 纯团聚作用
其次探讨团聚作用对冰浆存储过程冰晶粒径分布演化的影响。图4所示为冰晶只发生团聚作用时冰晶粒径分布随时间变化的曲线。由图4可知：冰晶粒径分布曲线随着时间的增加有不断向右移动的趋势，也就是说大粒径冰晶颗粒数量不断增加，小粒径冰晶颗粒数量不断减少。根据冰晶团聚机理可知多个小粒径冰晶颗粒不断团聚为大粒径冰晶颗粒，因此冰晶颗粒总数减少且冰晶颗粒主要集中在大粒径范围内，所以引起了图4中曲线峰值微弱右移同时峰值不断增小。与破碎作用效果相比，团聚作用效果更加明显。图4中冰晶粒径分布曲线变化幅度较大，但是同样随着存储时间的推移，冰晶粒径分布曲线的变化幅度越来越小。总体分析可知团聚作用效果主要是使冰晶颗粒尺度不断增加且数量不断减小，随着时间推移冰晶颗粒数量不断减小，冰晶数量减小到一定程度后，团聚作用效果就不断被削弱。
图4 冰晶在团聚作用下的粒径数量密度分布演化曲线
Fig. 4 Curve of ice crystal size distributions in time under aggregation
3.3 团聚和破碎共同作用
最后探讨团聚和破碎同时作用下冰浆存储过程冰晶粒径分布演化情况。在冰浆存储过程中冰晶粒径分布演化实际是团聚和破碎共同作用的结果，因此研究团聚与破碎共同作用下冰晶粒径分布演化有重要意义。图5所示为冰晶团聚和破碎共同作用时冰晶粒径分布随时间变化的曲线。根据图5中冰晶粒径分布演化曲线所示可知：1~200 μm范围内的冰晶颗粒数量密度略有增加，超过200 μm的小粒径冰晶颗粒数量密度不断减小，大颗粒数量密度不断增加。分析可得1~200 μm范围内的冰晶颗粒破碎作用大于团聚作用，而大于200 μm的冰晶颗粒团聚作用大于破碎作用。总体分析可知：冰晶粒径分布曲线峰值不断减小，小颗粒数量不断减少，大颗粒数量不断增加，冰晶总数量不断减少；对比纯团聚和纯破碎工况，在团聚和破碎共同作用下，团聚作用占主导地位，使得冰晶粒径分布曲线总的变化趋势是向着团聚作用效果方向发展。
图5 冰晶在团聚和破碎共同作用下的粒径数量密度分布演化曲线
Fig. 5 Curve of ice crystal size distributions in time under breakage and aggregation
4 结论
1) 引入数群平衡模型对冰浆存储过程冰晶粒径演化进行模拟，运用改进有限容积法对冰晶数群平衡方程进行求解，获得了冰晶粒径数量密度分布演化状况。
2) 实验拟合冰晶粒径分布演化曲线与数值模拟曲线基本吻合，且随着存储时间的增加，实验拟合曲线与数值模拟曲线重合度越来越高，也就是说实验结果与模拟结果误差越来越小。
3) 团聚作用和破碎作用是影响冰浆存储过程冰晶粒径演化的2个重要因素。破碎作用使小粒径冰晶颗粒数量增加，大粒径冰晶颗粒数量减少，冰晶颗粒总数量增加；团聚作用使大粒径冰晶颗粒数量增加，小粒径冰晶颗粒数量减少，冰晶颗粒总数量减少；两者共同作用下团聚作用占主导地位，冰晶粒径分布曲线向团聚作用效果方向变化。
参考文献：
[1] Egolf P W, Kauffeld M. From physical properties of ice slurries to industrial ice slurry applications[J]. International Journal of Refrigeration, 2005, 28(1): 4-12.
[2] 简夕忠. 二元冰浆流动与传热特性研究[D]. 南京: 南京航空航天大学宇航学院, 2008: 2-6.
JIAN Xizhong. Study of the flow and heat transfer characteristics of binary ice slurry[D]. Nanjing: Nanjing University of Aeronautics and Astronautics. School of Astronautics, 2008: 2-6.
[3] 柳飞. 动态冰晶制取系统的实验研究[D]. 武汉:华中科技大学能源与动力工程学院, 2005: 2-8.
LIU Fei. The experimental study of the dynamic ice-generating system[D]. Wuhan: Huazhong University of Science and Technology. Shool of Energy and Power Engineering, 2005: 2-8.
[4] Michael K, Masahiro K. Handbook of ice slurry[M]. Paris: International Institute of Refrigeration Press, 2005: 19-30.
[5] 赵海波, 郑楚光. 离散系统动力学演变过程的颗粒群平衡模拟[M]. 北京:科学出版社, 2008: 13-18.
ZHAO Haibo, ZHENG Chuguang. Population balance modeling for the process of the dynamic evolution in dispersed systems[M]. Beijing: Science Press, 2008: 13-18.
[6] Pronk P. Fluidized bed heat exchangers to prevent fouling in ice slurry systems and industrial crystallizers[D]. DELFT: Delft University of Technology. Department of Process & Energy, 2006: 139-168.
[7] Kumano H, Hirata T, Hagiwara Y, et al. Effects of storage on flow and heat transfer characteristics of ice slurry[J]. International Journal of Refrigeration, 2012, 35(1): 122-129.
[8] Yoshiyuki N A, Masayuki T. Study on ice storing charcteristics in dynamic-type ice storage system by using supercooled water. Effects of the supplying conditions of ice-slurry at deployment todistrict heating and cooling system[J]. International Journal of Refrigeration, 2005, 28(2): 73-82.
[9] Masayuki Y K. Ice water two-phase flow behavior in ice heat storage system[J]. International Journal of Refrigeration, 2001, 24(1): 639-651.
[10] 青春耀, 肖睿, 宋文吉, 等. 冰浆在蓄冰槽内的动态蓄冰过程模拟[C]//中国工程热物理学会多相流2009年学术会议论文集. 北京, 2009: 1-10.
QING Chunyao,XIAO Ri,SONG Wenji, et al. Simulation of the dynamic storage Process of ice slurry in the ice storage tank[C]//The 2009 Chinese Society of Engineering Thermalphsics Conference on Multiphase Flow. Beijing, 2009: 1-10.
[11] Pronk P, Hansen T, Ferreira C A I, et al. Time-dependent behavior of different ice slurries during storage[J]. International Journal of Refrigeration, 2005, 28(1): 27-36.
[12] 徐爱祥, 刘志强, 王小倩, 等. 破碎与团聚对冰浆存储过程冰晶粒径演化的影响[J]. 中南大学学报(自然科学版), 2013, 44(11): 4720-4725.
XU Aixiang, LIU Zhiqiang, WANG Xiaoqian, et al. Effect of breakage and aggregation on evolution of ice crystals during storage[J]. Journal of Central South University (Science and Technology), 2013, 44(11): 4720-4725.
[13] Perry R H, Green D W. Perry’s chemical engineers’ handbook[M]. 7th ed. New York: McGraw-Hill, 1997: 206-219.
[14] 苏军伟, 顾兆林. 离散相系统群体平衡模型的求解算法[J]. 中国科学: 化学, 2010, 40(2): 144-160.
SU Junwei, GU Zhaolin. Advances of solution methods of population balance equation for disperse phase system[J]. Science in China: Chemistry, 2010, 40(2): 144-160.
[15] Marchisio D L, Dennis Vigil R, Fox R O. Implementation of the quadrature method of moments in CFD codes for aggregation–breakage problems[J]. Chemical Engineering Science, 2003, 58: 3337-3351.
[16] 唐恒. 二元冰蓄冷系统添加剂的研究[D]. 上海: 上海海事大学商船学院, 2005: 1-9.
TANG Heng. Study of additives of the binary ice storage system[D]. Shanghai: Shanghai Maritime University. School of Merchant Marine, 2005: 1-9.
[17] Pronk P, Ferreira C A I, Witkamp G J. Effect of long-term ice slurry storage on crystal size distribution[C]//5th Workshop on Ice Slurries of the IIR. Stockholm, Sweden, 2002: 82-93.
[18] Meewisse J W. Fluidized bed ice slurry generator for enhanced secondary cooling systems[D]. DELFT: Delft University of Technology. Department of Process & Energy, 2004: 34-38.
(编辑 杨幼平)
收稿日期：2013-10-28；修回日期：2013-12-29
基金项目(Foundation item)：国家自然科学基金项目(51376198)；湖南省自然科学基金资助项目(11JJ22029)；教育部第47批留学回国人员科研启动基金资助项目(2013)(Project (51376198) supported by National Natural Science Foundation of China; Project (11JJ22029) supported by Natural Science Foundation of Hunan Province of China; Project (2013) supported by Scientific Research Foundation for Returned Overseas Chinese Scholars, Ministry of Education of China)
通信作者：刘志强(1970-)，男，湖南岳阳人，博士，教授，从事储能系统的热动力学研究；电话：0731-88876039；E-mail：liuzq@csu.edu.cn
