徐愛祥，劉志強，王小倩，趙騰磊

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

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

1 冰晶群體平衡模型

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

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

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

冰晶的團聚項為

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

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

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

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

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

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

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

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

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

1.2 求解方法

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

1.3 邊界條件與初始條件

1.3.1 邊界條件

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

1.3.2 初始條件

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

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

2 模型驗證

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

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

圖2 存儲22 h后粒徑概率密度的模擬值與實驗值Fig.2 Probability of ice crystal size distributions after 22 h storage obtained from model and experiment

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

3.1 純破碎作用

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

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

3.2 純團聚作用

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

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

3.3 破碎與團聚共同作用

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

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

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

4 結(jié)論

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

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

(3) 在冰晶存儲過程中，冰晶平均粒徑不斷長大，即從初始時刻的300 μm增長到存儲22 h后的392 μm，且模擬值與實驗值相對誤差在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)動力學(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.