基于改進(jìn)有限容積法的數(shù)群平衡方程數(shù)值方法研究

趙騰磊，劉志強(qiáng)，徐愛(ài)祥，何娜萍

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

數(shù)群平衡模型(population balance model, PBM)主要用于研究離散系統(tǒng)中顆粒粒度分布函數(shù)的時(shí)間演變過(guò)程[1]。它用數(shù)群平衡方程(population balance equation,PBE)描述由于顆粒增長(zhǎng)、團(tuán)聚、破碎等作用引起顆粒粒度分布隨時(shí)間的變化。數(shù)群平衡模型首先是被用在化學(xué)工程學(xué)科中[2-3]。隨著數(shù)值計(jì)算和高速計(jì)算機(jī)的發(fā)展，數(shù)群平衡模型被應(yīng)用在多個(gè)科學(xué)與工程學(xué)學(xué)科中[4-6]。實(shí)際運(yùn)用過(guò)程中，數(shù)群平衡方程精確的數(shù)值模擬需要一個(gè)合適的數(shù)值計(jì)算方法。許多該領(lǐng)域的學(xué)者已經(jīng)發(fā)展了多個(gè)數(shù)值計(jì)算方法解數(shù)群平衡方程，如Lee 等[7]的Monte Carlo 法，Hill 等[8]的有限差分法，Daniele 等[9]的正交矩方法，Everson 等[10]的有限元法。經(jīng)過(guò)研究發(fā)現(xiàn)，運(yùn)用這些方法所得數(shù)值解具有擴(kuò)散性或者不穩(wěn)定性[11]?，F(xiàn)在運(yùn)用最多的數(shù)值方法是改進(jìn)有限容積法，即特征法(method of characteristics, MOC)和有限容積法(finite volume schemes, FVS)結(jié)合的數(shù)值方法，這種數(shù)值計(jì)算方法是根據(jù)方程的數(shù)學(xué)特征進(jìn)行數(shù)學(xué)變換，用有限容積法對(duì)變換后的方程求解。Qamar等[4]對(duì)改進(jìn)有限容積法進(jìn)行了描述，并分別運(yùn)用改進(jìn)有限容積法和有限容積法對(duì)數(shù)群平衡方程進(jìn)行求解，把用2 種數(shù)值方法計(jì)算結(jié)果分別和分析解對(duì)比，分析了2 種數(shù)值方法的效果。數(shù)學(xué)變換是根據(jù)生長(zhǎng)項(xiàng)和團(tuán)聚項(xiàng)的數(shù)學(xué)特征對(duì)兩項(xiàng)進(jìn)行不同的變換，Qamar 等[4]只對(duì)幾種動(dòng)力學(xué)事件下粒子的分布進(jìn)行求解，并沒(méi)有對(duì)單獨(dú)考慮1 種變換以及同時(shí)考慮2 種變換的情況進(jìn)行計(jì)算，分析每一種變換對(duì)計(jì)算顆粒粒度分布的影響。為此，本文作者采用非均勻?qū)?shù)網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行離散，用有限容積法解數(shù)學(xué)變換后的數(shù)群平衡方程，分別求解純生長(zhǎng)事件、純團(tuán)聚事件和同時(shí)考慮生長(zhǎng)和團(tuán)聚作用下的顆粒粒度分布，把數(shù)值解和分析解進(jìn)行對(duì)比，分析數(shù)值求解方法的有效性；用所得數(shù)值結(jié)果與Qamar等[4]和Kumar 等[12]直接用有限容積法解數(shù)群平衡方程計(jì)算所得結(jié)果進(jìn)行比較，分析每一種變換對(duì)顆粒粒度分布的影響以及本文所用數(shù)值方法的優(yōu)點(diǎn)。

1 數(shù)群平衡模型

離散系統(tǒng)的數(shù)群平衡方程是描述體積分布函數(shù)的動(dòng)態(tài)變化。這里主要是用常用的數(shù)學(xué)方法對(duì)一維數(shù)群平衡方程進(jìn)行變換。Perry 等[13]給出了離散系統(tǒng)的一維數(shù)群平衡方程：

式中：f(t,x)為單位體積載流體中分散相的數(shù)量密度；G(t,x) 為生長(zhǎng)速率，表示t 時(shí)刻單位時(shí)間體積為x 的顆粒粒度增長(zhǎng)的大小[14]；Aagg(t,x) 為團(tuán)聚項(xiàng)；Bbreak(t,x)為破碎項(xiàng)；S(t,x) 為成核項(xiàng)。本文沒(méi)有涉及到成核項(xiàng)的變換，破碎項(xiàng)和團(tuán)聚項(xiàng)的變換是相同的，為了便于計(jì)算與對(duì)比，只考慮生長(zhǎng)和團(tuán)聚作用下顆粒粒度分布。故方程(1)變換為：

初始分布為：

團(tuán)聚項(xiàng)定義為：

方程(3)右邊第1 項(xiàng)表示體積為y 和x-y 的顆粒聚合為體積為x 的顆粒數(shù)量；第2 項(xiàng)表示體積為x 的顆粒和其他顆粒團(tuán)聚而減少的數(shù)量；β(t,x,y)表示顆粒的團(tuán)聚率，即t 時(shí)刻體積為x 與y 的顆粒碰撞后發(fā)生團(tuán)聚事件的速率[1]。

生長(zhǎng)作用主要影響顆粒體積大小，對(duì)顆?？倲?shù)量沒(méi)有影響；團(tuán)聚作用主要影響顆粒數(shù)量，顆?？傮w積不變。理論上顆粒會(huì)發(fā)生生長(zhǎng)和團(tuán)聚直到體積無(wú)窮大，但是實(shí)際顆粒會(huì)受到其他作用影響，體積不會(huì)無(wú)限增長(zhǎng)。體積比較大的顆粒數(shù)量較少，在方程的計(jì)算過(guò)程中會(huì)用截?cái)嗾`差對(duì)方程進(jìn)行處理。假如直接把有限容積法應(yīng)用到數(shù)群平衡方程中，需要把團(tuán)聚項(xiàng)作為源項(xiàng)處理。當(dāng)團(tuán)聚作用大于生長(zhǎng)作用時(shí)，可能導(dǎo)致截?cái)嗾`差較大，顆粒質(zhì)量損失較大，方程不守恒性增大，得到不真實(shí)的解。為了解決這個(gè)問(wèn)題，把團(tuán)聚項(xiàng)作為網(wǎng)格通量處理，即把數(shù)量守恒形式的數(shù)群平衡方程變換為質(zhì)量守恒形式的方程。

把方程(2)和(3)左右兩邊同時(shí)乘以x：

令：xAagg(t , x)=-?Fagg(t , x )/?x ，根據(jù)方程(3)得出：

2 求解方法

2.1 區(qū)域離散

為了運(yùn)用數(shù)值方法對(duì)方程求解，首先要對(duì)計(jì)算區(qū)域進(jìn)行離散。區(qū)域離散包括建立均勻網(wǎng)格和非均勻網(wǎng)格。當(dāng)顆粒粒度范圍較小時(shí)，非均勻幾何網(wǎng)格更加合理。本文所取計(jì)算范圍較小，故對(duì)計(jì)算區(qū)域建立非均勻幾何網(wǎng)格。一維幾何網(wǎng)格為：

式中：q 表示在小顆粒粒度范圍內(nèi)網(wǎng)格加密程度。

2.2 數(shù)值方法

當(dāng)生長(zhǎng)速率較大時(shí)，直接用有限容積法解數(shù)群平衡方程質(zhì)量損失較大，具體表現(xiàn)為方程的截?cái)嗾`差較大。在數(shù)群平衡方程中生長(zhǎng)項(xiàng)為方程的對(duì)流項(xiàng)，為了減小方程生長(zhǎng)項(xiàng)所引起的誤差，主要針對(duì)方程的對(duì)流項(xiàng)進(jìn)行變換。根據(jù)生長(zhǎng)速率的數(shù)學(xué)特征對(duì)原始數(shù)群平衡方程進(jìn)行數(shù)學(xué)變換，用有限容積法對(duì)變換后的方程進(jìn)行計(jì)算。根據(jù)線性標(biāo)量守恒定律，沿著顆粒粒度分布曲線移動(dòng)方向有唯一的特征曲線。假如解沿著傳播路徑的方向移動(dòng)，利用生長(zhǎng)速率的特征使方程(6)的對(duì)流項(xiàng)消失。

用以下方程替代生長(zhǎng)速率：dx/dt=G(t,x)

將式(6)在控制體積?i=[xi-1/2(t)，xi+1/2(t)]上積分：

Filbet 等[15]給出了團(tuán)聚項(xiàng)的離散形式：

其中：αi,k是當(dāng)xi+1/2- xk∈Ωαi,k-1(t )時(shí)，網(wǎng)格節(jié)點(diǎn)的數(shù)字符號(hào)。

初始條件：

在計(jì)算區(qū)域以外的部分?jǐn)?shù)量密度為零，即方程的邊界上的數(shù)量密度設(shè)置為零。也就是說(shuō)，在設(shè)定計(jì)算區(qū)域以外沒(méi)有更大和更小的顆粒。

3 數(shù)值結(jié)果與討論

用以上數(shù)群平衡方程數(shù)值計(jì)算方法計(jì)算純生長(zhǎng)、純團(tuán)聚、生長(zhǎng)和團(tuán)聚同時(shí)作用下顆粒粒度分布。為了說(shuō)明改進(jìn)有限容積法計(jì)算數(shù)群平衡方程的有效性，把數(shù)值計(jì)算結(jié)果和分析結(jié)果進(jìn)行對(duì)比。本文還給出了Qamar 等[4]和Kumar 等[12]直接用有限容積法解數(shù)群平衡方程所得結(jié)果，通過(guò)2 種數(shù)值方法計(jì)算結(jié)果與分析解的對(duì)比，分析每一種變換對(duì)顆粒粒度分布的影響以及本文所用數(shù)值方法的優(yōu)點(diǎn)。本文不針對(duì)具體的物理模型計(jì)算，只考慮數(shù)值計(jì)算方法對(duì)數(shù)群平衡模型的計(jì)算效果。本文用Matlab 自編程對(duì)數(shù)群平衡方程的離散公式求解。

3.1 生長(zhǎng)

對(duì)比在常數(shù)生長(zhǎng)(constant growth)和線性生長(zhǎng)(linear growth)條件下計(jì)算數(shù)群平衡方程的數(shù)值解和分析解。

常數(shù)生長(zhǎng)速率G(t,x)=G0，初始分布為：

式中：N0和x0分別代表初始時(shí)刻顆?？倲?shù)量和顆粒平均體積，分析解為

令G0=1，N0=5，x0=0.01，q=3，取60 個(gè)網(wǎng)格節(jié)點(diǎn)。

圖1 所示為2 種數(shù)值計(jì)算方法計(jì)算常數(shù)生長(zhǎng)條件下顆粒粒度分布數(shù)值結(jié)果和分析解的對(duì)比。圖1(a)所示為用改進(jìn)有限容積法計(jì)算的數(shù)值解和分析解的對(duì)比，從圖1(a)可知數(shù)值解和分析解基本一致。圖1(b)[11]所示為用有限容積法計(jì)算的數(shù)值結(jié)果和分析解的對(duì)比，從圖1(b)可知數(shù)值解和分析解基本一致，但是在小顆粒處出現(xiàn)了解的失真。由圖1 可知：小顆粒生長(zhǎng)速度較快，而大顆粒粒度幾乎不變。

圖1 顆粒常數(shù)生長(zhǎng)作用下的粒徑分布Fig.1 Ice crystal size distributions for constant growth process

線性生長(zhǎng)速率G(t,x)=G0x，初始分布為方程(13)。在數(shù)值計(jì)算中，令G0=1，N0=5，x0=0.01，q=3，取60個(gè)網(wǎng)格點(diǎn)。Kumar 等[16]給出分析解為：

圖2 所示為2 種數(shù)值計(jì)算方法計(jì)算線性生長(zhǎng)條件下顆粒粒度分布數(shù)值解和分析解的對(duì)比。圖2(a)所示為用改進(jìn)有限容積法計(jì)算的數(shù)值解和分析解的對(duì)比，從圖2(a)可知數(shù)值解和分析解一致性較好，隨著時(shí)間增加計(jì)算誤差增大，但數(shù)值解沒(méi)有出現(xiàn)不真實(shí)、不穩(wěn)定等現(xiàn)象。圖2(b)[11]所示為用有限容積法計(jì)算的數(shù)值結(jié)果和分析解的對(duì)比，從圖2(b)可知數(shù)值解和分析解基本一致，但在小顆粒處出現(xiàn)了解的失真。由圖2 可知：線性增長(zhǎng)使顆粒粒度不斷生長(zhǎng)的同時(shí)使顆粒粒度分布概率改變，小顆粒的概率不斷減小，大顆粒概率增加，總的顆粒數(shù)目沒(méi)有改變。

圖2 顆粒線性生長(zhǎng)作用下的粒徑分布Fig.2 Ice crystal size distributions for linear growth process

圖1 和圖2 所示為在純生長(zhǎng)作用下，顆粒粒度不斷長(zhǎng)大，而顆粒數(shù)量幾乎不變。從圖1 和圖2 可知：粒度較小的顆粒生長(zhǎng)速率較大，使得質(zhì)量損失對(duì)小顆粒粒度的分布影響較大，故小顆粒粒度分布的數(shù)值解出現(xiàn)失真現(xiàn)象。數(shù)群平衡方程是一個(gè)非標(biāo)準(zhǔn)輸運(yùn)控制方程，對(duì)流項(xiàng)中增加了生長(zhǎng)速率的影響，直接用有限容積法對(duì)方程進(jìn)行計(jì)算，計(jì)算所得的數(shù)值結(jié)果會(huì)出現(xiàn)失真現(xiàn)象。與有限容積法對(duì)比，改進(jìn)有限容積法使得方程中對(duì)流項(xiàng)消失，減少了生長(zhǎng)速率對(duì)方程計(jì)算結(jié)果的影響。從以上2 種計(jì)算結(jié)果可知對(duì)流項(xiàng)消失使得方程的解沒(méi)有出現(xiàn)失真。

3.2 團(tuán)聚

對(duì)比用 2 種數(shù)值方法計(jì)算常數(shù)團(tuán)聚(constant aggregation)和合團(tuán)聚(sum aggregation)問(wèn)題所得數(shù)值解和分析解。

常數(shù)團(tuán)聚速率β(t,x,y)=β0，初始分布為方程(13)，令β0=1，N0=5，x0=0.01，取91 個(gè)網(wǎng)格節(jié)點(diǎn)。Kumar等[12]給出常數(shù)團(tuán)聚作用下數(shù)群平衡方程的分析解為：

由圖3(a)和圖3(b)[12]對(duì)比可知：用2 種數(shù)值計(jì)算方法計(jì)算數(shù)群平衡方程所得數(shù)值解和分析解誤差很小，在±5%以?xún)?nèi)。在常數(shù)團(tuán)聚作用下用2 種數(shù)值方法計(jì)算的數(shù)值解沒(méi)有出現(xiàn)不穩(wěn)定、擴(kuò)散以及失真的現(xiàn)象。隨著時(shí)間不斷增加，小顆粒數(shù)量不斷減少，大顆粒數(shù)量不斷增加，顆粒總的體積幾乎沒(méi)有變化。

圖3 顆粒常數(shù)團(tuán)聚作用下的粒徑分布Fig.3 Ice crystal size distributions for constant aggregation process

合團(tuán)聚速率β(t,x,y)=β0(x+x′)，初始分布為方程(13)。在數(shù)值計(jì)算中，令β0=1，N0=5，x0=0.01，取91個(gè)網(wǎng)格節(jié)點(diǎn)。Kumar 等[12]給出分析解為

其中：τ=1-exp(-Ta)；I1是一階貝塞爾函數(shù)。

由圖4(a)和圖4(b)[12]對(duì)比可知：用2 種數(shù)值計(jì)算方法計(jì)算數(shù)群平衡方程所得數(shù)值解和分析解基本一致，與分析解的誤差都很小。在線性團(tuán)聚作用下用兩種數(shù)值方法計(jì)算的數(shù)值解同樣沒(méi)有出現(xiàn)不穩(wěn)定、擴(kuò)散或者失真等現(xiàn)象。隨著時(shí)間不斷增加，小顆粒數(shù)量不斷減少，大顆粒數(shù)量不斷增加，小顆粒減少速率比大顆粒增加速率小，顆?？偟捏w積幾乎沒(méi)有變化。

圖4 顆粒合團(tuán)聚作用下的粒徑分布Fig.4 Ice crystal size distributions for sum aggregation process

通過(guò)以上2 種團(tuán)聚事件作用下計(jì)算的顆粒粒度分布數(shù)值解和分析解的對(duì)比，可知在純團(tuán)聚作用下，2種方法計(jì)算結(jié)果誤差不大。團(tuán)聚作用只對(duì)顆粒數(shù)目有影響，對(duì)顆粒總的體積幾乎沒(méi)有影響，也就是說(shuō)質(zhì)量是守恒的，質(zhì)量損失較小，故用2 種方法計(jì)算的數(shù)值解沒(méi)有出現(xiàn)失真。

3.3 生長(zhǎng)和團(tuán)聚

對(duì)比用2 種數(shù)值方法計(jì)算同時(shí)考慮生長(zhǎng)和團(tuán)聚的數(shù)群平衡方程的數(shù)值解和分析解。本文考慮不同生長(zhǎng)速率和團(tuán)聚速率結(jié)合的3 種顆粒粒度演化情況，3 種情況初始分布都為方程(13)。3 種不同的動(dòng)力學(xué)過(guò)程如表1 所示[12,16-17]。

表1 生長(zhǎng)和團(tuán)聚作用下不同動(dòng)力學(xué)參數(shù)組合的動(dòng)力學(xué)過(guò)程Table 1 Kinetic processes of different kinetic parameters under growth and breakage

圖5 顆粒常數(shù)增長(zhǎng)和常數(shù)團(tuán)聚作用下的粒徑分布Fig.5 Ice crystal size distributions for constant growth and constant aggregation process

圖5 所示為表1 中第1 種情況的數(shù)值計(jì)算結(jié)果與分析解的對(duì)比。取G0=1，β0=100，N0=5，x0=0.01，q=4，取91 個(gè)網(wǎng)格節(jié)點(diǎn)。圖5(a)和圖5(b)[11]所示為用2 種求解方法計(jì)算的數(shù)值結(jié)果和分析解基本一致。甚至要求誤差非常小的情況下，這2 種方法所得的數(shù)值解也是非常有效的。然而有限容積法計(jì)算結(jié)果在左邊界出現(xiàn)失真現(xiàn)象，在其他區(qū)域與分析解是一致的。

圖5(a)所示在左邊界數(shù)值解和分析解有一些差別，但是不明顯。這是由于這個(gè)區(qū)間范圍的分析解有誤差。Kumar 等[15]指出在生長(zhǎng)速率和團(tuán)聚速率相等時(shí)，常數(shù)增長(zhǎng)和常數(shù)團(tuán)聚的分析解只是在大顆粒處有效；當(dāng)團(tuán)聚作用明顯時(shí)，分析解的有效范圍增大；當(dāng)生長(zhǎng)速率明顯時(shí)，分析解會(huì)產(chǎn)生一定的誤差。

圖6 所示為表1 中第2 種情況的數(shù)值計(jì)算結(jié)果。取G0=1，β0=10，N0=5，x0=0.01，q=4，取91 個(gè)網(wǎng)格節(jié)點(diǎn)。圖6(a)和圖6(b)[11]所示為用2 種求解方法計(jì)算的數(shù)值結(jié)果和分析解基本一致。有限容積法計(jì)算結(jié)果在左邊界同樣出現(xiàn)失真現(xiàn)象，在其他區(qū)域與分析解是一致的。這種情況下整個(gè)計(jì)算區(qū)域的顆粒團(tuán)聚作用比較明顯，分析解在整個(gè)計(jì)算區(qū)域都有效，所以圖6(a)中數(shù)值解沒(méi)有和分析解產(chǎn)生較大誤差。

圖7 所示為表1 中第3 種情況的數(shù)值計(jì)算結(jié)果。取G0=1，β0=1，N0=5，x0=0.01，q=4，取91 個(gè)網(wǎng)格節(jié)點(diǎn)。圖7(a)和圖7(b)[11]表示用2 種求解方法計(jì)算的數(shù)值結(jié)果和分析解基本一致。有限容積法計(jì)算結(jié)果在左邊界同樣出現(xiàn)失真現(xiàn)象，在其他區(qū)域與分析解同樣是一致的。這種情況同樣由于在整個(gè)計(jì)算域上團(tuán)聚作用比較明顯，所以圖7(a)中的分析解比較有效，沒(méi)有與數(shù)值解產(chǎn)生較大誤差。

圖6 顆粒線性增長(zhǎng)和常數(shù)團(tuán)聚作用下的粒徑分布Fig.6 Ice crystal size distributions for linear growth and constant aggregation process

圖7 顆粒常數(shù)增長(zhǎng)和合團(tuán)聚作用下的粒徑分布Fig.7 Ice crystal size distributions for linear growth and sum aggregation process

直接運(yùn)用有限容積法對(duì)數(shù)群平衡方程進(jìn)行計(jì)算時(shí)，只能保證顆粒數(shù)量守恒，不能保證顆粒質(zhì)量守恒。假如同時(shí)考慮生長(zhǎng)和團(tuán)聚作用對(duì)顆粒粒度的分布，當(dāng)團(tuán)聚作用較小時(shí)，顆粒質(zhì)量損失較小，方程數(shù)值解不會(huì)出現(xiàn)失真；當(dāng)團(tuán)聚作用比生長(zhǎng)作用明顯時(shí)，顆粒質(zhì)量損失較大，方程質(zhì)量不守恒性也會(huì)增大，方程數(shù)值解會(huì)出現(xiàn)失真。表1[17]中3 種情況都是團(tuán)聚速率比生長(zhǎng)速率大的情況，而且團(tuán)聚速率第1 種情況最大，第3 種情況最小。根據(jù)3 種情況計(jì)算結(jié)果可知：直接運(yùn)用有限容積法解數(shù)群平衡方程會(huì)出現(xiàn)解的擴(kuò)散，而且團(tuán)聚速率越大失真現(xiàn)象越明顯；而用經(jīng)過(guò)數(shù)學(xué)變換后的有限容積法計(jì)算減小了顆粒質(zhì)量損失，使方程不守恒性減小，方程數(shù)值解沒(méi)有出現(xiàn)失真。

4 結(jié)論

(1) 純生長(zhǎng)條件下求解數(shù)群平衡方程，與有限容積法相比，本文所用求解方法消除了方程中對(duì)流項(xiàng)的影響，減小了方程解的失真。

(2) 團(tuán)聚作用只是影響顆粒數(shù)量分布，幾乎沒(méi)有質(zhì)量損失。純團(tuán)聚條件下求解數(shù)群平衡方程，兩種數(shù)值求解方法所得數(shù)值結(jié)果都和分析解基本一致。

(3) 有限容積法只可以保證顆粒分布數(shù)量的守恒，不能保證顆粒質(zhì)量守恒。在生長(zhǎng)和團(tuán)聚同時(shí)作用條件下對(duì)數(shù)群平衡方程求解，當(dāng)團(tuán)聚作用比較大時(shí)，運(yùn)用有限容積法所得數(shù)值解會(huì)出現(xiàn)失真；而且當(dāng)團(tuán)聚作用越大時(shí)，數(shù)值解的失真現(xiàn)象越明顯。本文對(duì)數(shù)群平衡方程進(jìn)行數(shù)學(xué)變換，即使方程團(tuán)聚項(xiàng)作為網(wǎng)格通量處理，使得方程質(zhì)量守恒，數(shù)值解的失真現(xiàn)象消失。

[1] 趙海波, 鄭楚光. 離散系統(tǒng)動(dòng)力學(xué)演變過(guò)程的顆粒群平衡模擬[M]. 北京: 科學(xué)出版社, 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.

[2] Hulbert H M, Katz S. Some problems in particle technology[J].Chemical Engineering Science, 1995, 19: 555-574.

[3] Randolph A, Larson M A. Theory of particulate processes[M].2nd ed. San Diego: Academic Press, 1988: 7-19.

[4] Qamar S, Warnecke G, Elsner M P. On the solution of population balances for nucleation, growth, aggregation and breakage processes[J]. Chemical Engineering Science, 2009, 64:2088-2095.

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

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

[7] Lee K, Matsoukas T. Simultaneous coagulation and breakage using constant—N Monte Carlo[J]. Powder Technology, 2000,110: 82-89.

[8] Hill P J, Ng K M. New discretization procedure for the breakage equation[J]. AIChE Journal, 1995, 41: 1204-1216.

[9] Daniele L M, Dennis R V, Rodney O F. Quadrature method of moments for aggregation-breakage processes[J]. Colloid Interface Science, 2003, 258: 322-334.

[10] Everson R C, Eyre D, Campbell Q P. Spline method for solving continuous batch grinding and similarity equations[J].Computers & Chemical Engineering, 1997, 21: 1433-1440.

[11] Qamar S, Warnecke G. Numerical solution of population balance equations for nucleation, growth and aggregation processes[J].Computers & Chemical Engineering, 2007, 31: 1576-1589.

[12] Kumar J, WarnecKe G, Peglow M, et al. Comparison of numerical methods for solving population balance equations incorporating aggregation and breakage[J]. Power Technology,2009, 189: 218-299.

[13] Perry R H, Green D W. Perry’s chemical engineers’handbook[M]. 7th ed. New York: McGraw-Hill, 1997: 107-116.

[14] Ramkrishna D. Population balances theory and applications to particulate systems in engineering[M]. San Diego: Academic Press, 2000: 47-59.

[15] Filbet F, Lauren?ot P. Numerical simulation of the Smoluchowski coagulation equation[J]. SIAM Journal Scientific Computing, 2004, 25(6): 2004-2048.

[16] Kumar S, Ramkrishna D. On the solution of population balance equations by discretization-III. Nucleation, growth and aggregation of particles[J]. Chemical Engineering Science, 1997,52: 4659-4679.

[17] Ramabhadran T E, Peterson T W, Seinfeld J H. Dynamics of aerosol coagulation and condensation[J]. AIChE Journal, 1976,22: 840-851.