徐愛祥，劉志強(qiáng)，王小倩，趙騰磊

(中南大學(xué) 能源科學(xué)與工程學(xué)院，湖南 長沙，410083)

冰漿(ice slurry)是指含有大量懸浮冰晶粒子的固液兩相溶液，其具有較高的儲(chǔ)能密度和良好的流動(dòng)與換熱特性而受到人們廣泛關(guān)注[1?2]，并且在眾多領(lǐng)域具有廣闊的應(yīng)用前景[3?4]。冰漿存儲(chǔ)期間冰晶的數(shù)量和粒徑存在著顯著的變化，這些變化會(huì)對冰漿的流動(dòng)與傳熱特性產(chǎn)生重要影響[5]。Kozawa等[6?7]從宏觀角度研究了冰漿在存儲(chǔ)槽內(nèi)富集與分層現(xiàn)象，而這無法精確闡述冰漿存儲(chǔ)中冰晶演化過程與機(jī)理。Pronk等[8]從微觀角度出發(fā)，認(rèn)為存儲(chǔ)過程冰晶演化受磨損、團(tuán)聚及Oswald熟化效應(yīng)共同作用，而熟化效應(yīng)的影響較為顯著，但計(jì)算結(jié)果與實(shí)驗(yàn)存在較大差異。最近，Kumano[5]等分析發(fā)現(xiàn)團(tuán)聚與破碎作用對冰晶粒徑演化也有重要影響。上述研究表明，存儲(chǔ)過程中冰晶演化的研究還不夠全面，特別是微尺度下冰晶演化過程與機(jī)理還有待深入研究。群體平衡模型(population balance model,PBM)是描述離散體系中分散相實(shí)體尺度分布的有效方法[9]。PBM在化工領(lǐng)域的廣泛應(yīng)用[10?11]為冰漿存儲(chǔ)過程離散相間相互作用的研究帶來很大啟示。冰漿系統(tǒng)是典型的離散系統(tǒng)，其基本特點(diǎn)是載流溶液攜帶大量的冰晶顆粒群。冰漿存儲(chǔ)過程中冰晶生長與磨損、團(tuán)聚與破碎等微觀現(xiàn)象都是群體平衡方程(population balance equation, PBE)所要表達(dá)的內(nèi)容。為此，本文作者針對冰漿存儲(chǔ)過程離散特點(diǎn)建立冰晶顆粒群體平衡模型，通過Matlab自編程序，采用分區(qū)法對群體平衡方程進(jìn)行離散求解，著重討論破碎作用與團(tuán)聚作用對冰晶粒徑演化的影響。

1 冰晶群體平衡模型

1.1 數(shù)學(xué)建模

在建立冰漿存儲(chǔ)過程冰晶粒徑演化群體平衡數(shù)學(xué)模型時(shí)，進(jìn)行如下假設(shè)：(1) 冰晶粒子在空間上均勻分布，且具有各向同一特性；(2) 存儲(chǔ)過程冰晶的核化現(xiàn)象已不明顯，過冷度隨著冰晶的形成而逐漸變小，忽略冰晶的核化與生長；(3) 忽略外部空間對冰晶分布與演化的影響。因此，冰漿存儲(chǔ)過程中冰晶粒徑演化的零維空間群體平衡模型為

式中：Bbr為由于破碎動(dòng)力學(xué)行為造成冰晶數(shù)值密度的增加量，表示粒徑為x′的冰晶破碎成粒徑為x子顆粒的生產(chǎn)率；Dbr為由于破碎動(dòng)力學(xué)行為造成冰晶數(shù)值密度的減小量，表示冰晶破碎成子顆粒后粒徑為x原顆粒的損失率；b(x)為破碎率(breakage frequency)，表示單位時(shí)間粒徑為x的冰晶碰撞后發(fā)生破碎事件的速率；p(x|x′)為粒子概率密度函數(shù)(particle distribution function, PDF)，表示粒徑為x′冰晶破碎后會(huì)形成粒徑為x冰晶的概率；v(x′)為粒徑為x′的冰晶破碎一次產(chǎn)生粒徑為x冰晶粒子的平均個(gè)數(shù)，本文取v(x′)=2。

冰晶的團(tuán)聚項(xiàng)為

式中：Bag為由于團(tuán)聚動(dòng)力學(xué)行為造成冰晶數(shù)值密度的增加量，表示粒徑為x?x′的冰晶與粒徑為x′冰晶團(tuán)聚形成粒徑為x冰晶的生產(chǎn)率；Dag為由于團(tuán)聚動(dòng)力學(xué)行為造成冰晶數(shù)值密度的減小量，表示粒徑為x的冰晶在參與下一次的團(tuán)聚事件后會(huì)產(chǎn)生更大粒徑的冰晶，從而造成粒徑為x原顆粒的損失率；a(x,x′)為團(tuán)聚率(aggregation frequency)，表示單位時(shí)間粒徑為x與x′冰晶碰撞后發(fā)生團(tuán)聚事件的速率。需要注意的是：式(2)與式(3)積分項(xiàng)的上限∞指的是在計(jì)算時(shí)冰晶粒子粒徑的最大值。

根據(jù)模型假設(shè)，冰漿中冰晶粒子在慣性區(qū)域受流體擾動(dòng)作用仍保持各向同性，故不同粒徑的冰晶碰撞時(shí)破碎率[12]的表達(dá)式如下：

式中：ε為耗損率(dissipation rate)，表示單位質(zhì)量的冰漿在外界擾動(dòng)或內(nèi)部粒子間相互作用下所耗散的功率，W/kg；γ為冰晶與載流溶液間的表面張力，mN/m；φice為含冰率，%；ρ為載流溶液的密度，kg/m3；C1與C2為破碎率中常數(shù)，其取值分別為0.004 81和0.08。

粒徑為x′的冰晶破碎為粒徑x的冰晶概率分布函數(shù)[13]為

式中：c為形狀因子，當(dāng)冰晶為球形時(shí)，c可取為1。

根據(jù)上述假設(shè)，模型只考慮布朗碰撞作用下的團(tuán)聚，因此，描述冰晶粒子在液相中的團(tuán)聚率[12]為

式中：KB為Boltzman常數(shù)，J/K；T為冰漿存儲(chǔ)時(shí)的平衡溫度，K；μ為冰漿中載流溶液的動(dòng)力黏度，mPa·s。

模型中物性參數(shù)設(shè)置參照實(shí)驗(yàn)條件，計(jì)算所需要的各物性參數(shù)如表1所示[12]。實(shí)驗(yàn)中采用輕柔的攪拌方式以保持冰晶空間上的均勻分布，數(shù)值計(jì)算過程將其換算成耗損率。

表1 載流溶液(9.2% NaCl)的物性參數(shù)Table 1 Physical parameters of ice slurry solutions(9.2% NaCl)

1.2 求解方法

群體平衡方程(PBE)是一個(gè)典型的混合積分與微分方程，其解呈非線性雙曲形式，求解過程需采取特殊的處理方式以封閉方程。PBE的求解方法主要有以下3種[14]：(1) 分區(qū)法；(2) 矩法；(3) Monte Carlo法。其中，分區(qū)法的離散方式思路直觀、簡潔, 且能夠得到顆粒尺度譜的動(dòng)力學(xué)演變過程細(xì)節(jié)信息，計(jì)算精度與計(jì)算代價(jià)較適中。本文采用分區(qū)法對PBE進(jìn)行離散求解：將冰晶尺度分布曲線進(jìn)行離散，劃分為有限數(shù)目個(gè)區(qū)間，在每個(gè)區(qū)間內(nèi)針對某個(gè)冰晶粒徑的分布函數(shù)，建立平衡方程，再聯(lián)立求解這些平衡方程。需要指出的是：求解過程中利用基于核函數(shù)的積分中值定理封閉PBE右邊的團(tuán)聚與破碎項(xiàng)積分形式。上述求解過程通過Matlab自編程序來實(shí)現(xiàn)。

1.3 邊界條件與初始條件

1.3.1 邊界條件

在 0～10?11m 范圍的冰晶粒徑概率密度在方程中無法表達(dá)與計(jì)算，而據(jù)成核理論可知[15]，0～10?11m的冰晶粒子會(huì)有消融。因此，在任意時(shí)刻計(jì)算，可添加 1個(gè)與初始時(shí)刻相同的概率密度作為左邊界條件，即。

1.3.2 初始條件

考慮實(shí)驗(yàn)中測定的初始狀態(tài)冰晶粒徑概率密度近似呈正態(tài)分布[16]，在數(shù)值計(jì)算中，初始條件參照實(shí)驗(yàn)結(jié)果中冰晶的初始分布形式設(shè)定，并通過實(shí)驗(yàn)結(jié)果擬合而得到，擬合關(guān)聯(lián)式如下：

其中：x的取值范圍為 10?11＜x＜vM+1。

2 模型驗(yàn)證

為驗(yàn)證本文建立的模型可靠性，將冰晶粒徑概率密度及冰晶平均粒徑的模擬結(jié)果與 Pronk的實(shí)驗(yàn)結(jié)果[17]進(jìn)行對比，如圖1和圖2所示。從圖1可以看出：冰晶粒徑實(shí)驗(yàn)結(jié)果要比模擬結(jié)果稍偏高。從圖2可以看出：冰漿存儲(chǔ)22 h后，粒徑分布實(shí)驗(yàn)結(jié)果與模擬結(jié)果在曲線谷峰有一定偏差，其余位置較吻合。實(shí)驗(yàn)結(jié)果與模擬結(jié)果存在一定的偏差，其主要原因是：(1) 由于實(shí)驗(yàn)過程中冰晶并非規(guī)則球狀，所獲取的冰晶圖片也非規(guī)則的圓形，因此，冰晶粒徑是經(jīng)過專門的粒徑圖像處理而得到，而模擬中冰晶形狀因子則按規(guī)則的球形來計(jì)算；(2) 冰漿實(shí)驗(yàn)存儲(chǔ)過程可能涉及多種復(fù)雜微觀過程，而模擬計(jì)算忽略了冰晶生長與磨損，這也會(huì)給計(jì)算結(jié)果帶來偏差。通過定量比較，冰晶平均粒徑的實(shí)驗(yàn)結(jié)果與模擬結(jié)果相對誤差在 5%以內(nèi)，冰漿存儲(chǔ)22 h后冰晶粒徑概率密度基本吻合，說明本文模型準(zhǔn)確、可靠。

圖1 存儲(chǔ)期間冰晶平均粒徑的變化Fig.1 Variations of average ice crystal size during storage

圖2 存儲(chǔ)22 h后粒徑概率密度的模擬值與實(shí)驗(yàn)值Fig.2 Probability of ice crystal size distributions after 22 h storage obtained from model and experiment

3 模擬結(jié)果與分析

3.1 純破碎作用

首先考察純破碎工況下冰漿存儲(chǔ)過程冰晶粒徑分布與演化。圖3所示為理想狀態(tài)下冰晶碰撞只發(fā)生破碎現(xiàn)象時(shí)的粒徑概率密度。由圖3可知：隨存儲(chǔ)時(shí)間的推移，冰晶分布曲線整體有向左移動(dòng)的趨勢，也就是說，小粒徑冰晶出現(xiàn)的概率增大，大粒徑冰晶出現(xiàn)的概率減??；同時(shí)，峰值的左側(cè)更遠(yuǎn)離初始分布，而峰值的右側(cè)更靠近初始分布，這表明破碎作用對小粒徑的冰晶影響更大，而對大粒徑的影響較小。這主要是由于冰漿在存儲(chǔ)過程中破碎作用使大粒徑冰晶減少及顆粒冰晶不斷增加。總體上，分布曲線移動(dòng)幅度不大，且變化的幅度由快到慢，即在0～2 h內(nèi)冰晶粒徑有較小變化，2～22 h則無明顯變化。這表明在存儲(chǔ)過程中，破碎作用對冰晶的影響較小，且在短時(shí)間內(nèi)，其作用有一定的影響；隨著存儲(chǔ)時(shí)間的增加，其作用逐漸削弱。

圖3 冰晶在破碎作用下的粒徑概率密度Fig.3 Probability of ice crystal size distributions in time under breakage

3.2 純團(tuán)聚作用

圖4 所示為理想狀態(tài)下冰晶碰撞只發(fā)生團(tuán)聚現(xiàn)象時(shí)的粒徑概率密度。由圖4可知：與破碎作用下相反，隨著存儲(chǔ)時(shí)間的推移，團(tuán)聚作用下冰晶分布曲線整體有向右移動(dòng)，即小粒徑冰晶出現(xiàn)的概率減少，大粒徑冰晶出現(xiàn)概率增加；此外，峰值的左側(cè)與右側(cè)都遠(yuǎn)離初始分布，這表明團(tuán)聚作用對所有粒徑的冰晶都有影響。究其原因，主要是在存儲(chǔ)過程冰晶粒子群在團(tuán)聚作用下，顆粒與顆粒之間不斷的黏附而導(dǎo)致冰晶整體呈現(xiàn)變大趨勢?？傮w上來看，分布曲線移動(dòng)明顯，且變化的幅度由慢到快，即在0～2 h內(nèi)冰晶粒徑有較小變化，2～22 h則變化明顯。這說明在存儲(chǔ)過程中，團(tuán)聚作用對冰晶的影響是顯著的，且隨著存儲(chǔ)時(shí)間的增加，其作用也更加明顯。

圖4 冰晶在團(tuán)聚作用下的粒徑概率密度Fig.4 Probability of ice crystal size distributions in time under aggregation

3.3 破碎與團(tuán)聚共同作用

在實(shí)際存儲(chǔ)過程中，破碎作用與團(tuán)聚作用同時(shí)發(fā)生，共同作用，因此，很有必要分析破碎與團(tuán)聚共同作用工況下冰漿存儲(chǔ)過程冰晶粒徑分布與演化。圖 5所示為冰晶在破碎與團(tuán)聚共同作用下的粒徑概率密度，如圖5所示。從圖5可見：隨著冰漿存儲(chǔ)過程的推進(jìn)，冰晶粒徑概率密度曲線向右移動(dòng)且分布曲線坡度逐漸變緩，小粒徑冰晶出現(xiàn)的概率減小，大粒徑冰晶出現(xiàn)概率增大，隨時(shí)間推移冰晶粒徑分布范圍變廣。圖6所示為共同作用下冰晶粒徑分布函數(shù)。從圖6可以看出：初始時(shí)刻冰晶粒徑主要分布在100～500 μm，隨著存儲(chǔ)時(shí)間的增加，大粒徑冰晶比例明顯增加。對比純破碎作用、純團(tuán)聚作用及破碎與團(tuán)聚共同作用對粒徑分布與演化的影響可看出：冰漿存儲(chǔ)經(jīng)歷了破碎與團(tuán)聚這2個(gè)不同作用的微觀過程，2個(gè)“作用方向”相反，破碎作用下大粒徑冰晶產(chǎn)生小粒徑冰晶，而團(tuán)聚作用下小粒徑冰晶形成大粒徑冰晶；此外，在存儲(chǔ)初期，破碎與團(tuán)聚單獨(dú)作用下(見圖3與圖4)粒徑概率密度曲線的變化比兩者共同作用下(見圖 5)粒徑概率密度曲線的變化大，表明在共同作用時(shí)破碎與團(tuán)聚對冰晶粒徑的影響相互抵消；隨著存儲(chǔ)時(shí)間的增加，破碎作用不斷削弱，而團(tuán)聚作用依然明顯，因此，二者的相互作用由最初的破碎與團(tuán)聚“均衡作用”發(fā)展成團(tuán)聚“主導(dǎo)”。這一現(xiàn)象的具體表現(xiàn)是冰晶在存儲(chǔ)期間平均粒徑會(huì)增大，圖1證明了這一結(jié)論，冰晶的平均粒徑從初始時(shí)刻的 300 μm 增長到存儲(chǔ) 22 h后的392 μm。

圖5 冰晶在破碎與團(tuán)聚共同作用下粒徑概率密度Fig.5 Probability of ice crystal size distributions in time under breakage and aggregation

圖6 冰晶在破碎與團(tuán)聚共同作用下粒徑分布函數(shù)Fig.6 Distribution function of ice crystal size distributions in time under breakage and aggregation

4 結(jié)論

(1) 引入群體平衡模型可較好地模擬冰漿存儲(chǔ)過程的粒徑分布及演化，獲得了冰晶粒徑概率密度隨時(shí)間的變化關(guān)系：純破碎作用下小粒徑冰晶出現(xiàn)的概率增大，大粒徑冰晶出現(xiàn)概率減小，純團(tuán)聚及兩者共同作用下小粒徑冰晶出現(xiàn)的概率減少，大粒徑冰晶出現(xiàn)概率增加。

(2) 破碎與團(tuán)聚是冰晶粒徑演化的2個(gè)重要作用，在存儲(chǔ)過程中，兩者“作用方向”相反，且破碎與團(tuán)聚經(jīng)歷從最初的破碎與團(tuán)聚“均衡作用”發(fā)展成團(tuán)聚“主導(dǎo)”的過程。

(3) 在冰晶存儲(chǔ)過程中，冰晶平均粒徑不斷長大，即從初始時(shí)刻的300 μm增長到存儲(chǔ)22 h后的392 μm，且模擬值與實(shí)驗(yàn)值相對誤差在5%以內(nèi)。

[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]Ayel V, Lottin O, Peerhossaini H. Rheology, flow behaviour and heat transfer of ice slurries: a review of the state of the art[J].International Journal of Refrigeration, 2003, 26(1): 95?107.

[3]Kauffeld M, Wang M J, Goldstein V, et al. Ice slurry applications[J]. International Journal of Refrigeration, 2010,33(8): 1491?1505.

[4]Bellas I, Tassou S A. Present and future applications of ice slurries[J]. International Journal of Refrigeration, 2005, 28(1):115?121.

[5]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.

[6]Kozawa Y, Aizawa N, Tanino M. Study on ice storing characteristics in dynamic-type ice storage system by using supercooled water: Effects of the supplying conditions of ice-slurry at deployment to district heating and cooling system[J].International Journal of Refrigeration, 2005, 28(1): 73?82.

[7]青春耀, 肖睿, 宋文吉, 等. 冰漿在蓄冰槽內(nèi)的蓄冰特性及其均勻度研究[J]. 低溫與超導(dǎo), 2009, 37(5): 41?46.QING Chunyao, XIAO Rui, SONG Wenji, et al. Study on ice slurry storing characteristics and its uniformity in ice storage tank[J]. Cryogenics and Superconductivity, 2009, 37(5): 41?46.

[8]Pronk P, Hansen T M, 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.

[9]Ramkrishna D. Population balances: Theory and applications to particulate systems in engineering[M]. San Diego: Academic Press, 2000: 30?45.

[10]顧兆林, 蘇軍偉, 李云, 等. 兩相及多相體系的離散相行為與群體平衡模型[J]. 化學(xué)反應(yīng)工程與工藝, 2007, 23(2): 162?167.GU Zhaolin, SU Junwei, LI Yun, et al. Behaviors of the dispersed phase in the multiphase system and population balance model[J]. Chemical Reaction Engineering and Technology, 2007,23(2): 162?167.

[11]Gerstlauer A, Motz S, Mitrovic A, et al. Development, analysis and validation of population models for continuous and batch crystallizers[J]. Chemical Engineering Science, 2002, 57(20):4311?4327.

[12]Marchisio D L, Vigil R D, Fox R O. Implementation of the quadrature method of moments in CFD codes for aggregation-breakage problems[J]. Chemical Engineering Science, 2003, 58(15): 3337?3351.

[13]Hill P J, Ng K M. New discretization procedure for the breakage equation[J]. AIChE Journal, 1995, 41(5): 1204?1216.

[14]趙海波, 鄭楚光. 離散系統(tǒng)動(dòng)力學(xué)演變過程的顆粒群平衡模擬[M]. 北京: 科學(xué)出版社, 2008: 70?79.ZHAO Haibo, ZHENG Chuguang. Population balance modeling for the process of the dynamic evolution in dispersed systems[M]. Beijing: Science Press, 2008: 70?79.

[15]Chégnimonhan V, Josset C, Peerhossaini H. Ice slurry crystallization based on kinetic phase-change modeling[J].International Journal of Refrigeration, 2010, 33(8): 1559?1568.

[16]Hansen T M, Radosevic M, Kauffeld M. Behaviour of ice slurry in thermal storage systems[R]. Atlanta: ASHRAE, 2002:93?107.

[17]Pronk P, Ferreira C A I, Witkamp G J. Effect of long-term ice slurry storage on crystal size distribution[C]// Melinder A.5th Workshop on Ice Slurries of the IIR. Stockholm: International Institute of Refrigeration, 2002: 151?160.