基于改進有限容積法的數群平衡方程數值方法研究

趙騰磊，劉志強，徐愛祥，何娜萍

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

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

1 數群平衡模型

離散系統的數群平衡方程是描述體積分布函數的動態變化。這里主要是用常用的數學方法對一維數群平衡方程進行變換。Perry 等[13]給出了離散系統的一維數群平衡方程：

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

初始分布為：

團聚項定義為：

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

生長作用主要影響顆粒體積大小，對顆?？倲盗繘]有影響；團聚作用主要影響顆粒數量，顆?？傮w積不變。理論上顆粒會發生生長和團聚直到體積無窮大，但是實際顆粒會受到其他作用影響，體積不會無限增長。體積比較大的顆粒數量較少，在方程的計算過程中會用截斷誤差對方程進行處理。假如直接把有限容積法應用到數群平衡方程中，需要把團聚項作為源項處理。當團聚作用大于生長作用時，可能導致截斷誤差較大，顆粒質量損失較大，方程不守恒性增大，得到不真實的解。為了解決這個問題，把團聚項作為網格通量處理，即把數量守恒形式的數群平衡方程變換為質量守恒形式的方程。

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

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

2 求解方法

2.1 區域離散

為了運用數值方法對方程求解，首先要對計算區域進行離散。區域離散包括建立均勻網格和非均勻網格。當顆粒粒度范圍較小時，非均勻幾何網格更加合理。本文所取計算范圍較小，故對計算區域建立非均勻幾何網格。一維幾何網格為：

式中：q 表示在小顆粒粒度范圍內網格加密程度。

2.2 數值方法

當生長速率較大時，直接用有限容積法解數群平衡方程質量損失較大，具體表現為方程的截斷誤差較大。在數群平衡方程中生長項為方程的對流項，為了減小方程生長項所引起的誤差，主要針對方程的對流項進行變換。根據生長速率的數學特征對原始數群平衡方程進行數學變換，用有限容積法對變換后的方程進行計算。根據線性標量守恒定律，沿著顆粒粒度分布曲線移動方向有唯一的特征曲線。假如解沿著傳播路徑的方向移動，利用生長速率的特征使方程(6)的對流項消失。

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

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

Filbet 等[15]給出了團聚項的離散形式：

其中：αi,k是當xi+1/2- xk∈Ωαi,k-1(t )時，網格節點的數字符號。

初始條件：

在計算區域以外的部分數量密度為零，即方程的邊界上的數量密度設置為零。也就是說，在設定計算區域以外沒有更大和更小的顆粒。

3 數值結果與討論

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

3.1 生長

對比在常數生長(constant growth)和線性生長(linear growth)條件下計算數群平衡方程的數值解和分析解。

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

式中：N0和x0分別代表初始時刻顆?？倲盗亢皖w粒平均體積，分析解為

令G0=1，N0=5，x0=0.01，q=3，取60 個網格節點。

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

圖1 顆粒常數生長作用下的粒徑分布Fig.1 Ice crystal size distributions for constant growth process

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

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

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

圖1 和圖2 所示為在純生長作用下，顆粒粒度不斷長大，而顆粒數量幾乎不變。從圖1 和圖2 可知：粒度較小的顆粒生長速率較大，使得質量損失對小顆粒粒度的分布影響較大，故小顆粒粒度分布的數值解出現失真現象。數群平衡方程是一個非標準輸運控制方程，對流項中增加了生長速率的影響，直接用有限容積法對方程進行計算，計算所得的數值結果會出現失真現象。與有限容積法對比，改進有限容積法使得方程中對流項消失，減少了生長速率對方程計算結果的影響。從以上2 種計算結果可知對流項消失使得方程的解沒有出現失真。

3.2 團聚

對比用 2 種數值方法計算常數團聚(constant aggregation)和合團聚(sum aggregation)問題所得數值解和分析解。

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

由圖3(a)和圖3(b)[12]對比可知：用2 種數值計算方法計算數群平衡方程所得數值解和分析解誤差很小，在±5%以內。在常數團聚作用下用2 種數值方法計算的數值解沒有出現不穩定、擴散以及失真的現象。隨著時間不斷增加，小顆粒數量不斷減少，大顆粒數量不斷增加，顆?？偟捏w積幾乎沒有變化。

圖3 顆粒常數團聚作用下的粒徑分布Fig.3 Ice crystal size distributions for constant aggregation process

合團聚速率β(t,x,y)=β0(x+x′)，初始分布為方程(13)。在數值計算中，令β0=1，N0=5，x0=0.01，取91個網格節點。Kumar 等[12]給出分析解為

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

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

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

通過以上2 種團聚事件作用下計算的顆粒粒度分布數值解和分析解的對比，可知在純團聚作用下，2種方法計算結果誤差不大。團聚作用只對顆粒數目有影響，對顆?？偟捏w積幾乎沒有影響，也就是說質量是守恒的，質量損失較小，故用2 種方法計算的數值解沒有出現失真。

3.3 生長和團聚

對比用2 種數值方法計算同時考慮生長和團聚的數群平衡方程的數值解和分析解。本文考慮不同生長速率和團聚速率結合的3 種顆粒粒度演化情況，3 種情況初始分布都為方程(13)。3 種不同的動力學過程如表1 所示[12,16-17]。

表1 生長和團聚作用下不同動力學參數組合的動力學過程Table 1 Kinetic processes of different kinetic parameters under growth and breakage

圖5 顆粒常數增長和常數團聚作用下的粒徑分布Fig.5 Ice crystal size distributions for constant growth and constant aggregation process

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

圖5(a)所示在左邊界數值解和分析解有一些差別，但是不明顯。這是由于這個區間范圍的分析解有誤差。Kumar 等[15]指出在生長速率和團聚速率相等時，常數增長和常數團聚的分析解只是在大顆粒處有效；當團聚作用明顯時，分析解的有效范圍增大；當生長速率明顯時，分析解會產生一定的誤差。

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

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

圖6 顆粒線性增長和常數團聚作用下的粒徑分布Fig.6 Ice crystal size distributions for linear growth and constant aggregation process

圖7 顆粒常數增長和合團聚作用下的粒徑分布Fig.7 Ice crystal size distributions for linear growth and sum aggregation process

直接運用有限容積法對數群平衡方程進行計算時，只能保證顆粒數量守恒，不能保證顆粒質量守恒。假如同時考慮生長和團聚作用對顆粒粒度的分布，當團聚作用較小時，顆粒質量損失較小，方程數值解不會出現失真；當團聚作用比生長作用明顯時，顆粒質量損失較大，方程質量不守恒性也會增大，方程數值解會出現失真。表1[17]中3 種情況都是團聚速率比生長速率大的情況，而且團聚速率第1 種情況最大，第3 種情況最小。根據3 種情況計算結果可知：直接運用有限容積法解數群平衡方程會出現解的擴散，而且團聚速率越大失真現象越明顯；而用經過數學變換后的有限容積法計算減小了顆粒質量損失，使方程不守恒性減小，方程數值解沒有出現失真。

4 結論

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

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

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

[1] 趙海波, 鄭楚光. 離散系統動力學演變過程的顆粒群平衡模擬[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.

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