二項式分布在種群平衡模型模擬粒度分布中的應用

李振亮，盧培利，張代鈞，周志恩，張晟，何強

?

二項式分布在種群平衡模型模擬粒度分布中的應用

李振亮1,2，盧培利3,4，張代鈞3,4，周志恩2，張晟2，何強1,5

（1重慶大學城市建筑與環(huán)境工程學院，重慶400044；2重慶市環(huán)境科學研究院，重慶401147；3重慶大學煤礦災害動力學與控制國家重點實驗室，重慶400044；4重慶大學環(huán)境科學系，重慶400044；5重慶大學三峽庫區(qū)生態(tài)環(huán)境教育部重點實驗室，重慶400044）

提出了一種適用于幾何網(wǎng)格的子粒子二項式分布函數(shù)，并應用于種群平衡模型模擬活性污泥絮凝后的粒度分布。結果表明：與二元分布相比，該二項式分布可以得到更準確的粒度分布和平均粒度模擬值；通過校核二項式分布參數(shù)C的取值，可以提高粒度分布和平均粒度的模擬精度。相比于二元分布或正態(tài)分布只能描述一種類型的子粒子分布，該二項式分布具有較強的適應性，調整參數(shù)C的取值，可以得到更多可能的子粒子分布；參數(shù)C還可以表征粒子的破碎方式——較小的C值表征粒子具有較強的穩(wěn)定性，易破碎生成較大的子粒子；較大的C值表征粒子具有較弱的穩(wěn)定性，易破碎生成較小的子粒子。

種群平衡；二項式分布；粒度分布；數(shù)值模擬；破碎；子粒子

引 言

種群平衡模型（population balance model，PBM）是描述多相流體系中分散相大小分布的重要方法，在結晶體系、絮凝體系、粒子或氣泡形成體系中得到了廣泛應用[1-5]。應用種群平衡模型描述絮凝動力學，分析粒子聚并和破碎行為對粒度分布的影響，有助于更好地認識絮體粒子的形成、演變及運動規(guī)律[6-11]。

粒子可通過侵蝕（小粒子從大粒子表面侵蝕脫落）和破裂（大粒子分裂成兩個或更多較小的粒子）等方式發(fā)生破碎[12-13]，進而形成不同的子粒子分布[14]。由于難于實測子粒子分布規(guī)律，以往的種群平衡模型中一般采用假設的子粒子分布函數(shù)，如二元分布（binary distribution）、三元分布（ternary distribution）和高斯分布（Gaussian distribution）[15-16]。其中，假設粒子破裂生成兩個等大子粒子的二元分布應用最為廣泛，而高斯分布認為粒子破裂生成子粒子的概率密度符合高斯分布。值得注意的是，這些子粒子分布函數(shù)只能描述一種類型的子粒子分布，并不能描述其他可能的子粒子分布。

另外，如果應用離散區(qū)間法求解種群平衡模型，上述子粒子分布的應用還需結合計算網(wǎng)格，如二元分布多應用于幾何網(wǎng)格v=2v-1，高斯分布多應用于均勻網(wǎng)格（v=，為常數(shù)），且可近似為二項分布[17-18]

式中，v為粒子的體積，γ為粒子破碎為子粒子的分布函數(shù)。根據(jù)式（1），粒子（體積為，以下簡稱）破碎后生成兩個等大子粒子（2）的發(fā)生概率最大。種群平衡模型求解中，劃分網(wǎng)格的數(shù)量直接決定了計算量[19-20]；式（1）雖然簡化了高斯分布的計算，但是均勻網(wǎng)格對應的計算量較大（尤其是粒子粒度范圍較大時），從而限制了式（1）或高斯分布的應用。

本文提出了一種適用于幾何網(wǎng)格、能夠描述粒子破碎后子粒子多種可能分布的二項式分布函數(shù)，并將其應用于種群平衡模型對活性污泥絮凝動力學的模擬，通過對絮體粒子粒度分布和平均粒度模擬結果的對比，分析其適用性和優(yōu)越性。

1 模型與方法

1.1 種群平衡模型

描述粒子在聚并或破碎后的數(shù)量濃度變化率的種群平衡模型為[21]

式（2）右邊第1、2項分別為由于粒子聚并產(chǎn)生的生成項和消亡項，第3、4項分別為由于粒子破碎產(chǎn)生的生成項和消亡項；其中，為粒子（）的數(shù)量濃度,為聚并效率，()為粒子（）和粒子（）間的碰撞頻率，()為粒子（）破碎為粒子（）的概率密度函數(shù)，()為粒子（）的破碎速率。

應用離散區(qū)間法，式（2）可轉化為離散的種群平衡模型[22]

式中，N為粒子數(shù)量濃度，δ,k為Dirac Delta函數(shù)，η為分配系數(shù)

(4)

聚并效率通常被視作0～1的常數(shù)或者與粒度相關的函數(shù)[8,23-25]，可通過模型校核來確定。由于粒子的聚并和破碎是同時發(fā)生的，所以聚并效率和破碎速率系數(shù)往往同時進行校核確定[24,26]。碰撞頻率和破碎速率與粒度相關，考慮粒子間作用更依賴于粒子的最大粒度[27]。因此，碰撞頻率和破碎速率可分別描述為[28]

() =c()(6)

式（5）中，為速度梯度（本文只考慮了流體剪切引起的碰撞頻率），c為最大粒度；式（6）中為破碎速率系數(shù)。最大粒度c與等效粒度L的關系可近似為[28]

式中，P為基本粒度，2為二維分形維數(shù)。

1.2 子粒子分布函數(shù)

假設粒子均由基本粒子構成，所含基本粒子的數(shù)量可以代表粒子的類別[29]。以幾何網(wǎng)格（v=2v-1）為例，粒子類別如圖1所示。

只考慮粒子破碎生成兩個子粒子的情況，以粒子（P= 8）的破碎為例，存在4種破碎方式（圖2）。

圖2中，粒子（P= 3，5，6，7）被按比例分配至相鄰的粒子類別中；如1個P= 3的粒子可分配為1/2個P= 2的粒子和1/2個P= 4的粒子。假設4種破碎方式的概率相同，可計算出粒子破碎為子粒子的數(shù)量分布（表1）。

表1 粒子j破碎為子粒子i的數(shù)量分布

由表1可見，子粒子的生成概率符合一定的分布規(guī)律，即體積為原粒子體積1/2的子粒子的生成概率最大，如粒子（= 4）破碎后生成子粒子（= 3）的概率最大。將表1的子粒子數(shù)量分布歸一化處理后，子粒子概率分布可以較好地被二項式分布函數(shù)來描述（圖3）；該二項式分布函數(shù)為

式中，P為粒子破碎為子粒子的生成概率，為參數(shù)。由二項式分布定義可知，如果子粒子-1的生成概率最大，則×=-1。

圖3 二項式分布函數(shù)的子粒子概率分布模擬結果

Fig.3 Simulated probability distribution of daughter particles by binomial distribution function

○ analysis results; —— simulation results

已知生成概率P，則子粒子分布函數(shù)為

以上分析是以特定幾何網(wǎng)格（v=2v-1）劃分粒子類別為例，對于一般情形（v=kv1，1<<2），可以通過定義參數(shù)C來決定發(fā)生概率的分布

= (-C)/(10)

如果粒子（）破碎后子粒子（/2）的生成概率最大，則

C= lg2/lg(11)

1.3 參數(shù)校核

種群平衡模型的最終狀態(tài)是粒子的聚并與破碎達到平衡，即式（3）中dN/d= 0。將式（6）代入式（3），得

由式（12）可見，當粒度分布（數(shù)量濃度）一定時，聚并效率與破碎速率系數(shù)之比是唯一且可確定的值。

參數(shù)校核采用最小誤差法，定義

式中，為粒子類別總數(shù)，()和′() 分別為粒子體積分數(shù)的實測值和模擬值。

定義

式中，L和L分別為平均粒度的實測值和模擬值

(15)

式中，V為粒子總體積。

1.4 實測數(shù)據(jù)

采用活性污泥絮凝后的粒度分布實測數(shù)據(jù)（= 28.2 s-1）進行模型驗證與參數(shù)校核?；钚晕勰嗳∽猿鞘形鬯幚韽S，粒度分布采用激光粒度分析儀S3500（Microtrac, Enhanced型號）測量，粒子二維分形維數(shù)采用顯微圖像分析法測量[28,30]。通過S3500測量得到的累計體積分數(shù)可插值得到不同粒度粒子的體積分數(shù)，即某類粒子總體積（V）與所有類別粒子的總體積（V）之比。種群平衡模型描述的是粒子數(shù)量濃度（N），因此需要建立V與N之間的關系[7]

其中

(17)

式中，T為樣本總體積；為樣本濃度；f為絮體密度；為絮體中液相與固相之比

式中，ρ和s分別為液相和固相的密度。

本研究中，活性污泥樣本濃度為0.1 kg·m-3，ρ、s和f分別為1000、1040、1700 kg·m-3。

2 結果與討論

圖4為二項式分布和二元分布所對應的粒度分布模擬結果，參數(shù)校核結果見表2。當≠ 2時，二元分布的處理思路是：假設破碎生成兩個等大的子粒子，參考式（4）的分配方法計算出相應類別粒子的分配比例值，并將該分配比例值視作生成概率。由圖4可見，二項式分布與二元分布對應的模擬曲線形狀有差異，說明子粒子分布對粒度分布模擬結果有影響。

表2 二項式分布和二元分布對應的參數(shù)校核結果

表2中，二項式分布與二元分布對應的參數(shù)校核結果差異不顯著，與文獻[3]的結論一致。原因可能在于采用式（10）、式（11）計算的二項式分布結果中粒子破碎生成兩個等大子粒子的發(fā)生概率最大，與二元分布結果類似。無論是二項式分布，還是二元分布，離散區(qū)間劃分網(wǎng)格數(shù)量越多（值越?。?，粒度分布模擬誤差越小。相比于二元分布，二項式分布對粒度分布和平均粒度的模擬誤差更?。? 2.0的粒度分布模擬例外）。

二項式分布中，按照式（11）計算參數(shù)C的前提假設是粒子破碎方式（圖2）的發(fā)生概率相同。如果破碎方式的發(fā)生概率不同（如更傾向于侵蝕破碎或者完全破裂），C的值可能會發(fā)生變化。

分別選取3種不同的活性污泥（取自不同的培養(yǎng)階段，有機質含量不同），測得粒度分布差異明顯，“中心位置”和“形狀”均不同（圖5）。本文采用兩種不同的二項式分布處理方法對粒度分布進行了模擬，圖5（a）中，方法1[C為定值，由式（11）計算]較好地模擬了粒度分布“中心位置”的變化，但是對于“形狀”變化的模擬結果不理想；圖5（b）顯示，方法2（C為模型校核參數(shù)）較好地模擬了粒度分布“中心位置”和“形狀”的變化，粒度分布和平均粒度的模擬結果誤差更?。ū?）。由此可見，二項式分布具有較強的適應性，通過調整二項式分布參數(shù)C的取值，可以提高粒度分布和平均粒度的模擬精度。

表3 不同活性污泥粒度分布對應的模型參數(shù)校核結果

二項式分布參數(shù)C的取值還可以反映粒子的破碎方式。表3中，方法1中由式（11）計算得到的C= 1.47，表示粒子（）破碎后生成兩個等大子粒子（）的發(fā)生概率最大。C< 1.47表示粒子（）破碎生成子粒子（介于/2和之間）的發(fā)生概率最大；C>1.47表示粒子（）破碎生成子粒子（介于0和/2）的發(fā)生概率最大（圖6）。較小的C值意味著粒子具有較強的穩(wěn)定性，破碎后易生成較大的子粒子；反之，較大的C值意味著粒子具有較弱的穩(wěn)定性，破碎后易生成較小的子粒子。表3中，方法2校核得到的C值為0.97～1.78，活性污泥絮凝機理[31-32]與該值代表意義相符：當活性污泥EPS含量較低時，不易絮凝形成較大的絮體，而且絮體穩(wěn)定性較強，對應C的值較??；當活性污泥粒子EPS含量較高時，易絮凝形成較大的絮體，但是形成的絮體穩(wěn)定性較差，破碎易生成為較小的子絮體，對應C的值較大。

3 結 論

（1）提出的二項式分布函數(shù)[式（8）、式（10）、式（11）]適合于幾何網(wǎng)格，能夠描述粒子破碎后的子粒子分布；且具有較強的適應性，通過調整參數(shù)C的取值，得到可能更符合特定對象實際的子粒子分布。

（2）按式（11）計算參數(shù)C，可得到等大子粒子生成概率最大的二項式分布；將其應用于種群平衡模型模擬活性污泥絮凝動力學，與傳統(tǒng)的二元分布相比，可以得到更準確的粒度分布和平均粒度模擬結果。

（3）針對不同活性污泥粒度分布差異明顯的情況，通過校核二項式分布參數(shù)C的取值，可以提高粒度分布和平均粒度的模擬精度；校核所得的參數(shù)C的取值還可以反映粒子的破碎方式，較小的C值表明粒子具有較強的穩(wěn)定性，破碎易生成較大的子粒子；較大的C值表明粒子具有較弱的穩(wěn)定性，破碎易生成較小的子粒子。

符 號 說 明

b(v,w)——粒子破碎的概率密度函數(shù) Cp——二項式分布參數(shù) D2——二維分形維數(shù) E——破碎速率系數(shù)，m-1·s-1 Err1——粒度分布的平均誤差 Err2——平均粒度的平均誤差 G——平均速度梯度，s-1 L——粒子等效粒度，μm Lc——粒子最大粒度，μm Lm——粒子平均粒度，μm LP——基本粒度，μm m(i)——粒子體積分數(shù) Ni——粒子i的數(shù)量濃度，m-3 n——粒子種類數(shù)量 Pj,i——粒子j破碎為子粒子i的生成概率 p——參數(shù) S(i)——粒子i的破碎速率，s-1 Vf——所有類別粒子的總體積，m3 Vi——粒子 i的總體積，m3 VT——樣本總體積，m3 vi——粒子i的體積，m3 X——樣本濃度，kg·m-3 α——聚并效率 β(i,j)——粒子i和粒子j間的碰撞頻率，m3·s-1 γj,i——粒子j破碎為子粒子i的分布函數(shù) δj,k——Dirac Delta函數(shù) ηi——分配系數(shù) ρf——絮體密度，kg·m-3 ρl——液相密度，kg·m-3 ρs——固相密度，kg·m-3

References

[1] 鄭建祥, 許帥, 王京陽,等. 基于微分代數(shù)積分矩量法的聚并器超細粒子聚團研究[J]. 化工學報, 2017, 68(1): 119-128. ZHENG J X, XU S,WANG J Y,Simulation of ultrafine particle aggregation in aggregation device by differential-algebraic quadrature method of moments[J]. CIESC Journal, 2017, 68(1): 119-128.

[2] 王鐵峰. 氣液-漿-反應器流體力學行為的實驗研究和數(shù)值模擬[D]. 北京: 清華大學, 2004. WANG T F. Experimental study and numerical simulation on the hydrodynamics in gas-liquid (slurry) reactors [D]. Beijing: Tsinghua University, 2004.

[3] BRIESEN H. Simulation of crystal size and shape by means of a reduced two-dimensional population balance model[J]. Chemical Engineering Science, 2006, 61(1): 104-112.

[4] LEE K F, DOSTA M, MCGUIRE A D,Development of a multi-compartment population balance model for high-shear wet granulation with discrete element method[J]. Computers & Chemical Engineering, 2017, 99: 171-184.

[5] WANG B, MOSBACH S, SCHMUTZHARD S,Modelling soot formation from wall films in a gasoline direct injection engine using a detailed population balance model[J]. Applied Energy, 2016, 163: 154-166.

[6] DUCOSTE J. A two-scale PBM for modeling turbulent flocculation in water treatment processes[J]. Chemical Engineering Science, 2002, 57(12): 2157-2168.

[7] NOPENS I, KOEGST T, MAHIEU K,PBM and activated sludge flocculation: from experimental data to calibrated model[J]. AIChE Journal, 2005, 51(5): 1548-1557.

[8] DING A, HOUNSLOW M J, BIGGS C A. Population balance modelling of activated sludge flocculation: investigating the size dependence of aggregation, breakage and collision efficiency[J]. Chemical Engineering Science, 2006, 61(1): 63-74.

[9] YEOW Y L, LIOW J, LEONG Y. A general procedure for obtaining the evolving particle-size distribution of flocculating suspensions[J]. AIChE Journal, 2012, 58(10): 3043-3053.

[10] JELDRES R I, CONCHA F, TOLEDO P G. Population balance modelling of particle flocculation with attention to aggregate restructuring and permeability[J]. Advances in Colloid & Interface Science, 2015, 224: 62.

[11] LI Z L, LU P L, ZHANG D J,Population balance modeling of activated sludge flocculation: investigating the influence of extracellular polymeric substances (EPS) content and zeta potential on flocculation dynamics[J]. Separation & Purification Technology, 2016, 162: 91-100.

[12] JARVIS P, JEFFERSON B, GREGORY J,A review of floc strength and breakage[J]. Water Research, 2005, 39(14): 3121–3137.

[13] YUAN Y, FARNOOD R R. Strength and breakage of activated sludge flocs[J]. Powder Technology, 2010, 199(2): 111-119.

[14] SCHUETZ S, PIESCHE M. A model of the coagulation process with solid particles and flocs in a turbulent flow[J]. Chemical Engineering Science, 2002, 57(20): 4357-4368.

[15] SPICER P T, PRATSINIS S E. Coagulation and fragmentation: universal steady-state particle-size distribution[J]. AIChE Journal, 1996, 42(6): 1612-1620.

[16] ZHANG J J, LI X Y. Modeling particle-size distribution dynamics in a flocculation system[J]. AIChE Journal, 2003, 49(7): 1870-1882.

[17] SERRA T, CASAMITJANA X. Modelling the aggregation and break-up of fractal aggregates in a shear flow[J]. Flow, Turbulence and Combustion, 1997, 59(2): 255-268.

[18] MAGGI F, MIETTA F, WINTERWERP J C. Effect of variable fractal dimension on the floc size distribution of suspended cohesive sediment[J]. Journal of Hydrology, 2007, 343(1): 43-55.

[19] 蘇軍偉, 顧兆林, XU X Y. 離散相系統(tǒng)種群平衡模型的求解算法[J]. 中國科學: 化學, 2010, 40(2): 144-160. SU J W, GU Z L, XU X Y. Advances of solution methods of population balance equation for disperse phase system[J]. Scientia Sinica Chimica, 2010, 40(2): 144-160.

[20] LI Z L, ZHOU Z E, ZHANG S,Comparison of the accuracy and performance of different numbers of classes in discretised solution method for population balance model[J]. International Journal of Chemical Engineering, 2016, 2: 1-6.

[21] RAMKRISHNA D. Population Balances: Theory and Applications to Particulate Systems in Engineering[M]. London: Academic Press, 2000.

[22] KUMAR S, RAMKRISHNA D. On the solution of population balance equations by discretization(Ⅱ): A fixed pivot technique[J]. Chemical Engineering Science, 1996, 51(8): 1311-1332.

[23] THOMAS D N, JUDD S J, FAWCETT N. Flocculation modelling: a review[J]. Water Research, 1999, 33(7): 1579-1592.

[24] LI X Y, LOGAN B. Collision frequencies between fractal aggregates and small particles in a turbulently sheared fluid[J]. Environmental Science & Technology, 1997, 31(4): 1237-1242.

[25] BIGGS C A, LANT P A. Modelling activated sludge flocculation using population balances[J]. Powder Technology, 2002, 124(3): 201-211.

[26] MIETTA F, CHASSAGNE C, VERNEY R,On the behavior of mud floc size distribution: model calibration and model behavior[J]. Ocean Dynamics, 2011, 61(2): 257-271.

[27] FLESCH J C, SPICER P T, PRATSINIS S E. Laminar and turbulent shear‐induced flocculation of fractal aggregates[J]. AIChE Journal, 1999, 45(5): 1114-1124.

[28] LI Z L , ZHANG D J, LU P L,Population balance model and calibration method for simulating the time evolution of floc size distribution of activated sludge flocculation[J]. Desalination and Water Treatment, 2017, 67: 41-50.

[29] MAGGI F. Flocculation dynamics of cohesive sediment[D]. Delft : Delft University of Technology, 2005.

[30] 李振亮, 張代鈞, 盧培利,等. 活性污泥粒度分布與分形維數(shù)的影響因素[J]. 環(huán)境科學, 2013, 34(10): 3975-3980. LI Z L, ZHANG D J, LU P L,Influencing factors of floc size distribution and fractal dimension of activated sludge[J]. Environmental Science, 2013, 34(10): 3975-3980.

[31] WILEN B M, JIN B, LANT P. The influence of key chemical constituents in activated sludge on surface and flocculating properties[J]. Water Research, 2003, 37(9): 2127-39.

[32] SHENG G P, YU H Q, LI X Y. Extracellular polymeric substances (EPS) of microbial aggregates in biological wastewater treatment systems: a review[J]. Biotechnology Advances, 2010, 28(6): 882-894.

Application of binomial distribution for daughter particles in simulation of particle size distribution by population balance model

LI Zhenliang1,2, LU Peili3,4, ZHANG Daijun3,4,ZHOU Zhi’en2, ZHANG Sheng2, HE Qiang1,5

(1College of Urban Construction and Environmental Engineering, Chongqing University, Chongqing 400044, China;2Chongqing Research Academy of Environmental Sciences, Chongqing 401147, China;3State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400044, China;4Department of Environmental Science, Chongqing University, Chongqing 400044, China;5Key Laboratory of the Three Gorges Reservoir Region’s Eco-environment, Ministry of Education, Chongqing University, Chongqing 400045, China)

A binomial distribution function suitable for geometric grid was proposed and applied in population balance model (PBM) to simulate the particle size distribution (PSD) of activated sludge after flocculation. The results showed that the PSD and mean size simulated by using binomial distribution give better agreement with the experimental data than those simulated by using binary distribution. The accuracy of simulation of PSD and mean size can be improved by calibrating the value of parameterCof binomial distribution function. Different from other daughter particle distribution which can only describe one-type daughter particle distribution, the binomial distribution function shows a strong adaptability and can present more probable daughter particle distributions through adjusting the parameterC. Moreover, the value of parameterCmight also characterize breakage behavior of particle: a smallerCvalue might imply that the particle display strong stability leading to large numbers of large daughter particle, whereas a largerCvalue might imply that the particle display weak stability and the large-scale fragmentation resulting in large numbers of small daughter particle.

population balance; binomial distribution; particle size distribution; numerical simulation; breakage; daughter particle

10.11949/j.issn.0438-1157.20170675

TQ 028.5

A

0438—1157（2017）09—3397—07

2017-05-25收到初稿，2017-06-29收到修改稿。

盧培利。

李振亮（1981—），男，博士，副教授。

國家自然科學基金項目（5160091049）；重慶市科委專項基金項目（2015CSTC-JBKY-01605）；重慶市基礎科學與前沿技術研究專項項目（CSTC2016jcyjA0506）。

2017-05-25.

Prof.LU Peili, lupl@cqu.edu.cn

supported by the National Natural Science Foundation of China (5160091049), the Natural Science Foundation of Chongqing (2015CSTC-JBKY-01605) and the Chongqing Research Program of Basic Research and Frontier Technology(CSTC2016jcyjA0506).