3D格子Boltzmann傳質模型模擬生物膜降解有機污水

楊艷霞，李 靜

（1. 太原理工大學熱能工程系，太原 030024；2. 低品位能源利用技術及系統教育部重點實驗室，重慶大學，重慶 400044；3. 太原理工大學建筑與能源應用工程系，太原 030024）

0 引 言

人類活動的多樣性和工業化程度的提高導致各種生活及工業污水大量排放，嚴重破壞了人類健康及生態環境平衡，這種嚴峻的形勢促進了污水處理工藝的新發展[1-3]。膜生物污水處理法是一種新型且高效的處理技術，已廣泛應用于處理各種有機污水。生物膜法是一種最常用的細胞固定技術，生物膜對污水中有害物質具有抗毒性，能夠增加反應器內的生物量，提高生化反應的穩定性和反應器性能[4]。其中，光合細菌能夠通過吸收太陽能降解污水中有機成分供自身新陳代謝生長，同時生成產物氫，因此，光合細菌生物膜（光生物膜）污水處理法是一種極具發展前途的處理技術[1,5-6]。生物膜是反應器內生化反應的關鍵組成部分[7]，它是一種具有微小孔隙的多孔介質，其結構對反應器內主流區及生物膜內流動傳質及生化反應過程有很大的影響，因此對生物膜結構特性的研究尤為重要。試驗研究中，主要是借助掃描電鏡、光學儀器等設備的高分辨率來獲取多孔介質的平面圖像[8-10]，這種方法主要受限于可視化、設備精度等條件，且費用昂貴。因此，本文將采用數值重構方法獲得生物膜多孔介質結構，進而對其結構特性進行模擬研究。

格子 Boltzmann方法是一種介觀尺度的數值計算方法，它是基于動力學理論，將流體離散為大量的流體粒子，通過跟蹤每個粒子的分布函數，利用統計學的方法研究流體運動及傳輸規律[11-12]，已成功地應用于流動傳熱傳質[13-15]、多相流[16-17]、多孔介質流動[18-20]、化學反應[19,21]等領域。與傳統算法不同，該方法無需求解非線性偏微分方程，是一種離散模型[22-23]，且算法簡單、容易處理各種復雜邊界條件[22-24]。本文將利用3D格子Botlzmann傳質模型對膜生物反應器內的生化降解過程進行模擬計算。生物膜多孔介質結構將通過四參數隨機生成法重構獲得[25]。由于三維計算量較大，將網格細化模型與格子Boltzmann模型[26-27]耦合來提高計算效率，同時保證計算精度。文中通過改變各參數重構得到不同結構的生物膜，研究分析其孔隙率、孔隙分布對膜生物反應器內流動傳質、生化反應過程及反應器性能的影響，進而對試驗研究進行預測和指導。

1 3D格子Boltzmann模型

文中采用3D格子Boltzmann傳質模型，離散方向i上的分布函數 fi和 gi,σ分別用于描述流場和σ組分的濃度場，表達式為[9]：

式中x為坐標矢量；ei為粒子離散速度矢量；t為時刻；δt，δx為時間和空間步長；格子速度 c=δx/δt；Rσ為 σ 組分的無量綱反應源項；τν，τσ分別為對應于 fi和 gi,σ的無量綱松弛時間。

對應于流場和濃度場的平衡態分布函數，的表達式為[9]：

描述流場時采用三維十五速（D3Q15）模型[28]，如圖 1所示為各離散方向。在 D3Q15模型中，權系數 wi為：w0=2/9；w1-6=1/9；w7-14=1/72。粒子離散速度矢量為：e0=0；e1-6=(±1,0,0)c, (0,±1,0)c, (0,0,±1)c；e7-14=(±1,±1, ±1)c。在保證精度要求的情況下，描述濃度場時粒子離散方向由15點降為7點[9]。模型中，常數Ki=1/2，常數Ji,σ=(1- J0)/6,i=1-6，其中 0≤J0≤1。

根據上述分布函數，利用統計學方法可計算得到宏觀參數密度ρ、速度u和σ組分的濃度cσ[29]：

圖1 D3Q15模型離散節點粒子分布Fig.1 Discretized velocity space of D3Q15 model

2 計算模型

計算中，為了便于比較，對各參數均進行了無量綱處理：X=x/H，U=u/u0，Cσ=cσ/c0。

2.1 模型驗證

采用一個對流-擴散-反應問題來驗證模型的正確性。溶液以速度U0、濃度Cin進入矩形區域反應，反應源項為R=kC (k=2)。給定各邊界條件為：進口邊界為定速度和定濃度邊界條件；四周邊界均為周期性邊界條件；出口處速度梯度為0，濃度為Cout。計算中，采用了網格細化技術[26]，將計算區域劃分為兩部分，0≤X≤L/2為粗網格，L/2≤X≤L為細網格，粗細網格比例 m =δx,c/δx,f=2。圖 2為計算區域內濃度場分布及沿中心線上濃度分布曲線。從圖中可看到，粗細網格界面處濃度的連續性很好，且LB模擬結果與理論解有很好的吻合度，從而證明了模型及程序的正確性。

圖2 濃度場分布和沿中心線上的濃度分布曲線Fig.2 Concentration contour and concentration profile along centerline

2.2 物理模型

圖 3所示為平板式膜生物反應器結構示意圖，在反應器下表面均勻生長一層穩定的光合生物膜，有機廢水溶液以一定的流速和濃度均勻流入反應器內，有機物成分（葡萄糖）通過生物膜界面向生物膜內傳輸，被生物膜內光合菌降解，同時生成產物（氫）。底物和產物反應源項r1、r2的表達式為[30]：

式 中 μmax= 0 .26exp(- 1 .2(I0/(I0)opt- 1 )2)， m=0.76 exp(-2.8(I /(I ) -1 )2)，α = 0 .019 2exp(-9 .5(I0/ (I0)opt- 1 )2)，0 0 opt

其中Cx，Yx/s，ks，β分別為細胞密度，細胞得率，飽和常數和產氫動力學常數，其值分別為：0.76 kg/m3，0.85，5.204 kg/m3和0.41 h-1。c1為局部底物濃度，I0為光照強度，最佳光照強度(I0)opt=6 000 lx。溶液進口流速和底物濃度分別為：60 mL/h和60 mmol/L，溶液pH值為7.0，溫度為30 ℃。

圖3 三維膜生物反應器物理模型Fig.3 Geometric configuration of 3D membrane bioreactor

反應器內生物膜是一種多孔介質，其孔隙結構通過四參數隨機生成法數值重構獲得。該方法可以通過調整各參數獲得不同結構的多孔介質。為了與 3D格子Boltzmann模型（D3Q15）相一致，生物膜數值重構將在14個離散方向上（除靜止點0方向）按一定的生長概率pi重構生長，見圖1。然后將獲得的多孔介質骨架參數導入格子 Boltzmann模型進行模擬計算其內部的流動傳質及生化反應過程。

設定各邊界條件為：反應器四周均為壁面，速度為U=0，濃度梯度為 0；進口邊界為流速 U0，濃度 C0；出口邊界為自由出流，速度梯度和濃度梯度均為0。由于三維模型的計算量較大，因此模型中耦合了網格局部細化處理技術來提高計算效率且不影響計算精度和穩定性。生物膜區域采用細網格，主流區采用粗網格，粗細網格步長比m=2。

以孔隙率ε=0.5，各方向生長概率pi相同的多孔介質（結構0）的模擬計算來進行網格無關性驗證。網格密度分別為H/δx,c=15，20，25，30時，計算得到底物降解效率η分別為51.09%，51.01%，50.97%，50.93%。當H/δx,c≥20已滿足精度要求。本文考慮了計算時間成本和計算精度，故選取H/δx,c=25的網格密度。粗細網格模塊的網格數分別為 51×26×11 和 101×51×29。

為了評估反應器性能，對底物降解效率進行了計算:

式中min、mout分別為進出反應器內的底物質量。

3 結果與分析

膜生物反應器內，生物膜結構會直接影響主流區與生物膜內流動傳質及生化反應過程，進而影響反應器性能。首先研究分析生物膜多孔介質孔隙率的影響，如圖4a所示為各方向生長概率相同（pi=0.005）條件下，數值重構獲得不同孔隙率（ε=0.3，0.5，0.7）的生物膜多孔結構及其反應器內的流線分布。生物膜孔隙率較小時（如ε=0.3），多孔介質結構比較緊湊密實，孔隙較小且分布較均勻；孔隙率較大時（如ε=0.7），生物膜結構變得松散，除均勻分散的小孔隙，還出現一些較大的孔隙通道。由于入口效應（右側為流體進口），入口段流線有一定的曲率；主流區流線逐漸變得平滑，而生物膜內流線不連續且形變較大。當孔隙率較小時，由于流體流通通道較小，生物膜內流線較短且有較大的彎曲度；隨著孔隙率的增大，生物膜內流線逐漸增長且曲率減小。這主要是由于較大孔隙率使得流體的流通區域增大，流動阻力減小，較多的流體更容易進入生物膜內部，流動更為順暢。

圖4b～4c為圖4a流場條件下，具有不同孔隙率的生物膜內底物和產物濃度等值面圖。從圖4b中可看到，左側進口處底物濃度較高，沿著 X方向上底物濃度逐漸降低，這是由于沿著流動方向底物隨流體向生物膜內輸運且沿程被微生物降解，使得底物濃度在 X方向上逐漸減小。從生物膜界面至生物膜內部（沿Z反方向上），底物濃度逐漸減小，這是由于生物膜內較大傳質阻力，流體速度降低使得生物膜內部的底物負載較少，同時生物膜內微生物降解消耗底物，導致越深入生物膜內部底物濃度越低，且下游區域更為明顯。隨著孔隙率的增大，生物膜內底物濃度明顯增加，這是由于孔隙率增大流動阻力減小，流體滲透能力增強，更多底物負載隨流體向生物膜內部輸運供微生物降解。從圖4c可以看到，生物膜內產物濃度沿X方向上逐漸增加。如對圖4b所分析，沿流動方向上底物被生物膜內微生物降解，同時生成的產物隨流體的對流擴散作用被攜帶至下游區域，所以沿 X方向上產物濃度逐漸增加。沿 Z方向（生物膜厚度）上產物濃度逐漸降低，這說明生物膜界面處流體較強的對流擴散作用使得產物被流體快速地攜帶進入主流區及下游區域，而生物膜內部傳質阻力較大，降解生成的產物不能較快地向主流區傳遞，使得生物膜內部產物濃度梯度較小，這將不利于生化降解過程。隨著孔隙率的增大，生物膜內產物濃度逐漸減小，且較大孔隙率時更為明顯（如 ε=0.7）。

圖 5給出了孔隙率對反應器內底物降解效率的影響規律。隨著孔隙率的增加，底物降解效率逐漸增大，且在孔隙率ε=0.5時達到最大，50.97%；隨孔隙率的繼續增大，底物降解效率則呈下降趨勢。這是由于孔隙率較小時（ε＜0.5），孔隙率的增加擴大了生物膜內流體流通截面，流體在生物膜內的滲透能力增強，傳質過程隨之增強，從而更多底物負載隨流體向生物膜內輸運，為微生物生化降解提供充足的原料；同時由于傳質能力的增強，生物膜內生成的產物也較快地隨流體向生物膜外傳遞，從而增大了生物膜內產物濃度梯度，促進了生物膜內微生物新陳代謝能力，有利于生化降解反應，因此反應器內底物降解效率會隨生物膜孔隙率的增大而增加。但是隨著孔隙率的繼續增大（ε＞0.5），生物膜內流通截面增大，且有較大尺寸的孔隙出現，更多流體將從流動阻力較小的大孔隙流過（見圖4a中ε=0.7）；孔隙率的增大也會使得流體滲透能力明顯增強，生物膜內底物濃度增大，而過高的底物濃度會導致底物抑制現象；同時流速增加縮短了水動力停留時間，底物沒有充足的時間被微生物降解而被流體帶入下游區域，因此反應器內底物降解效率隨孔隙率的進一步增大而降低。

此外，孔隙率一定時，各方向生長概率pi也會影響生物膜結構，進而影響反應器內流動傳質及底物降解性能。因此，ε=0.5時，數值重構了5種不同的生物膜結構，結構0：p1-14=0.005；結構1：p3-4=0.01，p1,2,5-14=0.005；結構 2：p5-6=0.01，p1-4,7-14=0.005；結構 3：p7-8,11-12=0.01，p1-6,9-10,13-14=0.005；結構 4：p9-10,13-14=0.01，p1-8,11-12=0.005。圖6a所示為孔隙率ε=0.5時，數值重構獲得的生物膜結構1～結構4，結構0如圖4中ε=0.5如所示。圖6b和6c給出了圖6a生物膜結構條件下，生物膜內底物和產物濃度等值面圖。比較5種結構的生物膜，可以看到底物濃度變化較小，而產物變化較為明顯。結構0和結構 2生物膜底部區域內濃度梯度較小，產物濃度較高將導致產物抑制現象，這將不利于生物膜內生化降解過程。

圖4 不同孔隙率的生物膜內流線分布、底物和產物濃度等值面圖（pi=0.005）Fig.4 Streamlines distribution, substrate and product concentration contours in biofilm with various porosities for pi=0.005

圖5 孔隙率對生物膜反應器內降解效率的影響Fig.5 Effect of porosity on substrate consumption efficiency in membrane bioreactor

表1所示為孔隙率ε=0.5時，5種生物膜結構對反應器內底物降解性能的影響。可以看到，生物膜結構不同使得底物降解效率有明顯的變化。結構 2的底物降解效率最低，47.45%；結構0次之，50.97%；而結構1底物降解效率最高，52.54%。這是由于生物膜各方向生長概率 pi不同使得生物膜孔隙分布及結構形態不同，這將直接影響生物膜內底物和產物流動及傳輸特性和生化反應過程。結果表明：結構 1的生物膜多孔結構特性有利于底物及產物在主流區和生物膜內的傳輸，傳質過程增強，促進了生物膜內微生物的新陳代謝，加快了生化反應過程，因此結構 1的生物膜反應器內底物降解效率最高。結構4、3、0的底物降解效率逐漸降低，結構 2時底物降解效率最低。這說明結構 2生物膜孔隙結構分布不利于物質的傳輸，傳質阻力較大使得生物膜內生成的產物不能及時向生物膜外傳遞，不斷積聚在生物膜內使得生物膜內部濃度較高而出現產物抑制現象，下游區域產物濃度較低（見圖6c），這也不利于生物膜內生化反應過程，因此結構 2條件下膜生物反應器內底物降解效率最低。此外，還將LB模擬結果和文獻[31]中相同條件下（各參數見2.2節）的試驗值和模擬值進行了比較，與試驗值相比誤差在0.1%～9%，從而證實了格子Boltzmann模型的可行性。與文獻[31]中的數值模型相比，LB模擬結果精度較高，且LB模型可以從孔隙尺度上研究分析生物膜多孔介質的結構特性，且易于處理多孔介質的復雜邊界條件，具有較大的優勢。

圖6 各方向生長概率不同的生物膜結構及生物膜內底物和產物濃度等值面圖（ε=0.5）Fig.6 Biofilm structures with various growth probabilities, substrate and product concentration contours in biofilm at ε=0.5

表1 不同結構生物膜對反應器內降解效率的影響Table 1 Substrate consumption efficiency in bioreactor with different biofilm structures %

4 結 論

將3D格子Boltzmann傳質模型耦合多孔介質四參數隨機生成，重構得到不同結構的生物膜多孔介質，進而對膜生物反應器內的流動傳質及生化反應過程進行模擬計算。研究分析了生物膜孔隙率、生長概率對反應器內流場、濃度場及反應器性能的影響規律。此外，將LB模擬結果與試驗結果進行比較，證明了該模型的可行性。模擬結果表明：隨著生物膜孔隙率的增加，底物降解效率逐漸增大，且在孔隙率 ε=0.5時達到最大值 50.97%，隨孔隙率的繼續增大則呈下降趨勢；各方向生長概率 pi不同所獲得的生物膜孔隙分布、結構形態、比表面積不同，進而影響底物降解效率，生物膜為結構1（p3-4=0.01，p1,2,5-14= 0.005）時，膜生物反應器內底物降解效率最高，52.54%，而結構 2（p5-6=0.01，p1-4,7-14=0.005）時則底物降解效率最低，47.45%，研究結果將對反應器的優化具有一定的預測和指導作用。

[1] Mostafa A, Elsamadony M, El-Dissouky A, et al. Biological H2potential harvested from complex gelatinaceous wastewater via attached versus suspended growth culture anaerobes [J]. Bioresour Technol, 2017, 231: 9－18

[2] Liu R, Mao Y, Shen C, et al. Can biofilm affect alum sludge adsorption: An engineering scope in a novel biofilm reactor for wastewater treatment[J]. Chemical Engineering Journal,2017, 328: 683－690.

[3] Asadi N, Alavijeh M K, Zilouei H. Development of a mathematical methodology to investigate biohydrogen production from regional and national agricultural crop residues: A case study of iran[J]. Int J Hydrogen Energy,2017, 42: 1989－2007.

[4] Kars G, Gündüz U. Towards a super h2 producer:Improvements in photofermentative biohydrogen production by genetic manipulations[J]. Int J Hydrogen Energy, 2010, 35:6646－6656.

[5] Chen C Y, Chang J S. Enhancing phototropic hydrogen production by solid-carrier assisted fermentation and internal optical-fiber illumination[J]. Process Biochem, 2006, 41:2041－2049.

[6] Kapdan I K, Kargi F. Bio-hydrogen production from waste materials[J]. Enzyme Microb Technol, 2006, 38: 569－582.

[7] Sheng G P, Yu H Q, Li X Y. Extracellular polymeric substances (eps) of microbial aggregates in biological wastewater treatment systems: A review[J]. Biotechnol Adv,2010, 28: 882－894.

[8] Khan F, Enzmann F, Kersten M, et al. 3d simulation of the permeability tensor in a soil aggregate on basis of nanotomographic imaging and lbe solver[J]. J Soil Sediment,2012, 12: 86－96.

[9] Sullivan S P, Gladden L F, Johns M L. 3d chemical reactor lb simulations[J]. Math Comput Simulat, 2006, 72: 206－211.

[10] Sullivan S P, Sani F M, Johns M L, et al. Simulation of packed bed reactors using lattice Boltzmann methods[J].Chem Eng Sci, 2005, 60: 3405－3418.

[11] Yan G. W, Chen Y. S, Hu S. X. Simple lattice Boltzmann model for simulating flows with shock wave[J]. Phys Rev E,1999, 59: 454－459.

[12] Raabe D. Overview of the lattice Boltzmann method for nano-and microscale fluid dynamics in materials science and engineering[J]. Modell Simul Mater Sci Eng, 2004, 12: R13.

[13] 段智英，郭玉明，王福貴. 基于格子Boltzmann方法分析果蔬真空冷凍干燥凍干速率[J]. 農業工程學報，2016，32:258－264.Duan Zhiying, Guo Yuming, Wang Fugui. Vacuum freeze-drying rate of fruits and vegetables based on lattice Boltzmann method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016,32: 258－264. (in Chinese with English abstract)

[14] Lu J H, Lei H Y, Dai C S. A simple difference method for lattice Boltzmann algorithm to simulate conjugate heat transfer[J]. Int J Heat Mass Transfer, 2017, 114: 268－276.

[15] Paradis H, Andersson M, Sunden B. Modeling of mass and charge transport in a solid oxide fuel cell anode structure by a 3d lattice Boltzmann approach[J]. Heat and Mass Transfer,2016, 52: 1529－1540.

[16] Xie C, Zhang J, Wang M. Lattice Boltzmann modeling of non-newtonian multiphase fluid displacement[J]. Chinese Journal of Computational Physics, 2016, 33: 147－155.

[17] Parker R R, Klausner J F, Mei R W. Supersonic two-phase impinging jet heat transfer[J]. Journal of Heat Transfer-Transactions of the Asme, 2013, 135:

[18] Wang T, Gao Q, Chen J, et al. Lattice Boltzmann simulation of mixed convection in an enclosure filled with porous medium[J]. Chinese Journal of Computational Physics, 2017,34: 39－46.

[19] Lei T, Meng X, Guo Z. Lattice Boltzmann study on influence of chemical reaction on mixing of miscible fluids with viscous instability in porous media[J]. Chinese Journal of Computational Physics, 2016, 33: 399－409.

[20] 趙凱，宣益民，李強. 基于格子Boltzmann方法的復雜多孔介質內雙擴散效應的對流傳熱傳質機理研究[J]. 科學通報，2010，55: 94－102.

[21] Xie C Y, Wang J K, Wang D, et al. Lattice Boltzmann modeling of thermal conduction in composites with thermal contact resistance[J]. Commun Comput Phys, 2015, 17: 1037－1055.

[22] 郭照立，鄭楚光. 格子Boltzmann 方法的原理及應用[M].北京：科學出版社，2008.

[23] 何雅玲，王勇，李慶. 格子Boltzmann方法的理論及應用[M]. 北京：科學出版社，2008.

[24] Chang C, Liu C H, Lin C A. Boundary conditions for lattice Boltzmann simulations with complex geometry flows[J].Comput Math Appl, 2009, 58: 940－949.

[25] Wang M, Wang J K, Pan N, et al. Mesoscopic predictions of the effective thermal conductivity for microscale random porous media[J]. Phys Rev E, 2007, 75: 036702.

[26] Yu D Z, Girimaji S S. Multi-block lattice Boltzmann method:Extension to 3d and validation in turbulence[J]. Physica A,2006, 362: 118－124.

[27] Yu D Z, Mei R W, Shyy W. A multi‐block lattice boltzmann method for viscous fluid flows[J]. Int J Numer Methods Fluids, 2002, 39: 99－120.

[28] Qian Y H, d'Humieres D, Lallemand P. Lattice bgk models for navier-stokes equation[J]. Europhys Lett, 1992, 17: 479－484.

[29] Sukop M C, Thorne D T. Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers[M]. Springer Verlag, 2006.

[30] Liao Q, Liu D M, Ye D D, et al. Mathematical modeling of two-phase flow and transport in an immobilized-cell photobioreactor[J]. Int J Hydrogen Energy, 2011, 36: 13939－13948.

[31] 郭成龍. 光合細菌生物膜反應器內傳輸特性及產氫性能強化[D]. 重慶：重慶大學，2013.Guo Chenglong. Transport Characteristics in Biofilm Photobioreactor with Photosynthetic Bacteria and Enhancement of Hydrogen Production Performance[D].Chongqing: Chongqing University, 2013. (in Chinese with English abstract)