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)