基于介尺度PBM模型的生物反应器放大模拟及实验研究

1、基于介尺度PBM模型的生物反应器放大模拟及实验研究进入21世纪，生物技术因资源可再生、环境友好等特点受到人们越来越多的关注。工业生物技术已成为世界各国争相发展的重点项目，其中生物反应器更是占据工业生物技术的中心位置。生物反应周期长、连续操作、效率低的特点导致了生物反应器具有高能耗、低设备利用率的缺点，因此高效节能的反应器成为科技与工程研发核心。由于生物反应器的设计和放大过程中存在着诸多困难1-2，大型生物反应器在设计初期很难通过实验的方式提供可靠的设计依据。而计算流体力学（CFD）技术可有效提供反应器内的流场、传质、传热等信息，因此利用CFD相关软件进行数值模拟已成为大型生物反应器设计的重要途

2、径3-4。耦合CFD与群体平衡模型(PBM)，涵盖了实际存在的气泡聚并、破碎等气泡动力学行为，从而可准确预测生物反应器内流场分布以及传质能力。随着对气泡动力学行为的深入认识，人们开发了多种气泡聚并与破裂的模型5-6。在气泡聚并模型方面，Luo模型7被广泛应用于鼓泡塔模拟，但研究表明该模型会过高估计聚并频率，因此杨宁等8-9将气液体系内介尺度能量的描述和PBM模型相结合，在气液体系稳定的状态下，介尺度能量趋于最小，与此同时，不同尺寸气泡的破碎和聚并状态达到平衡，故修正之后的模型能够更加准确地估计聚并频率。在气泡破碎模型方面，Luo等10将能谱与涡能量相结合，他们认为涡不仅具有尺度，还具有能量，由

3、此提出了湍流涡诱导碰撞气泡破碎的模型。Han等11基于表面能量密度的增加提出了一种新型破碎依据，并将表面振荡和大尺度涡引起的破碎影响纳入模型；此外，Han等12-13还改进了数学模型，使其能够适用于湍流的全能谱和气泡的多次破碎。Solsvik等14的模型考虑了黏流剪切对气泡破碎的影响。Shi等15提出了气泡诱导湍流的修正破碎模型，考虑了气泡诱导湍流影响下涡旋的平均湍流速度和对应的特征波数/长度尺度。以上模型为模拟气泡聚并破裂的动力学行为提供了充分的理论基础，但是生物反应器内存在着湍流涡和气泡群介尺度结构，如何选择合适的模型阐明这种介尺度结构仍然面临巨大的挑战。本课题组在5 L通气搅拌式工业生物

4、反应器中通过组合不同聚并、破碎模型的方式，对八种模型组进行考察，分别从体积传质系数kLa和气泡尺寸两个方面与实验值对比，这八种组合模型包括以下四种诱导气泡碰撞的因素，它们分别为湍流涡随机运动诱导碰撞、浮升力差异诱导碰撞、大气泡的尾流夹带诱导碰撞和速度梯度差诱导碰撞。结果显示，与实验数据最为吻合的是基于介尺度理念的修正聚并模型与考虑黏流剪切的破碎模型组合。整个计算过程在商业软件Ansys Fluent 17.0平台上进行。由于小型生物反应器的工业应用前景不大，故本文在选定最优气泡聚并破碎模型的基础上，通过叶轮末端剪切力相等的放大原则将5 L生物反应器放大到400 m3。搅拌是化工生产过程中的重要

5、环节，开发新型搅拌桨一直都是化学工程领域的热点16-18，本文通过软件模拟新开发了四种组合桨，考察不同桨型对生物反应器内气泡分散的效果，在分析了各个组合桨的流场特性、混合死区分布等特征之后，选出400 m3生物反应器的最优组合桨，为大型通气搅拌式工业生物反应器的设计与研究提供参考。本文组合气泡不同的聚并效率和破碎频率，得到了贴合实验值的气泡聚并破碎模型；并通过放大准则得到400 m3通气搅拌式工业生物反应器相关的结构和工况参数，进而使用最优聚并破碎模型模拟不同桨型的气泡分散情况。1 实验装置与方法1.1 实验装置实验装置如图1所示，该生物反应器外壁是具有椭圆封头的圆柱形玻璃材质，内壁有四块不锈

6、钢挡板，挡板高232 mm，宽20 mm，厚度为3 mm。生物反应器内径T=150 mm，液面高度为298 mm，封头曲面高度30 mm。采用三层涡轮刚性桨，桨叶为不锈钢材质，桨叶直径D=50 mm,叶片长度17 mm，叶片宽度17 mm，叶片厚度2 mm，桨叶间距H0=70 mm，下层桨离底高度C=35 mm，气体分布器是内径为46 mm、宽度为12 mm的不锈钢圆环，圆环上方等距离分布10个直径为1 mm圆孔，气体从中空的圆环通过圆孔向生物反应器喷出，生物反应器内设有溶氧电极，可间隔36 s记录氧气浓度。图1图1生物反应器实验装置Fig.1Bioreactor experimental s

7、etup1.2 实验方法由于反应器内气相氧气浓度恒定，且假设混合过程是理想的，所以体积传质系数kLa的计算式19为CS-CICS-CI0=e-kLat-kLape-t/p1-kLap(1)式中，CI0为初始溶解氧浓度；p为DO电极时间常数；CS为饱和溶解氧浓度；CI为实时溶解氧浓度。如果反应器的时间常数R=1/kLa，远大于p。则式（1）可简化为CS-CICS-CI0=e-kLat(2)kLat=-lnCS-CI/CS-CI0(3)因为饱和溶解氧浓度CS和初始溶解氧浓度CI0在整个实验过程中不会变化，故只需记录不同时刻t的实时溶解氧浓度CI，便可绘制函数-ln(CS-CI)/(CS-CI0)的

8、图像。由式(3)可知，kLa可通过函数-ln(CS-CI)/(CS-CI0)对时间t求导得到。实验开始前，先注纯水到生物反应器指定液位，校正溶氧值后开启搅拌（转速可调节200、300、400 r/min），打开连接氮气钢瓶的气体流量计，调节流量为4 L/min，通过溶氧电极记录体系中溶氧数据的变化；体系内的溶氧量会逐渐降低，待平衡后，停止通入氮气，打开氧气阀通入氧气，继续监测溶氧值的变化，每组转速重复3次实验，以此减少实验设备对结果的影响。通过高速摄像机(Phanton v2640)记录气泡信息，采集速度为1000帧/秒，采集时间为30 s；得到气泡图像后，进一步处理得到气泡尺寸分布。2 数学

9、模型和数值模拟2.1 控制方程采用Eulerian-Eulerian双流体方法进行三维数值模拟20-21。该方法将气相和液相分别视为分散相和连续相。两相滞留量满足相容条件g+l=1(4)各相的连续性平衡方程为ggt+ggug=0llt+llul=0(5)各相的动量平衡方程为ggugt+ggugug=-gp+geff,gug+ugT+ggg-Ml,lgsllult+llulul=-lp+leff,lul+ulT+llg+Ml,lg(6)2.2 相间作用力根据 Sanyal等22-23的研究结果，本文模拟只考虑曳力，忽略其他相间作用力，曳力系数选择 Schiller-Naumann模型24模拟气液

10、体系相间作用力。FgasD=-FliquidD=34gCDdblugas-uliquidugas-uliquid(7)其中，曳力系数为CD,0=24Re1+0.15Re0.687Re10000.44Re1000(8)式中，CD是气泡群曳力系数;CD,0为单气泡曳力系数，CD=CD,0(1-g)p，本模拟中p=1。2.3 湍流模型本文采用RNGk-湍流模型25来预测气泡塔内流体的流动特性。方程为llklt+llulkl=llkkl+lGk,l-lll+llk,llllt+llull=lll+llklC10Gk,l-C20ll+lle,l(9)其中k,l=34gCDg,lkg,l-2k1+ug,l

11、d,e,l=C30lklk,l,C10=1.42,C20=1.68,C30=1.3式中，k,l和e,l表示考虑气泡诱导湍流的源项。2.4 群平衡模型2.4.1 聚并模型气泡聚并速率通常采用Luo模型7，为cVi,Vj=wijTPcdi,djMc(10)其中，wijT和Pc(di,dj)分别为碰撞频率和聚并效率；Mc为聚并速率的修正系数。湍流波动引起的气泡碰撞频率26(Prince and Bl)为wijT=4di+dj2213di23+dj231/2(11)Luo7提出的聚并效率模型为Pcdi,dj=exp-0.751+xij21+xij30.5gl+0.50.51+xij3Weij0.5We

12、ij0.5=ldiAij2,Aij=1.4313di23+dj2312(12)杨宁等8-9将气液体系内介尺度能量的描述和PBM模型相结合，气液体系在稳定状态下，介尺度能量趋于最小，同时不同直径之间的气泡破碎和聚并处于动态平衡，从而提出了对聚并模型的修正方法。其中Luo聚并模型7对应修正因子关系式见式（13）。Mc=0.13491+1.13403e-Ug0.01148Ug0.101m/s0.10164+0.53737e-Ug0.04961Ug0.101m/s(13)本文所研究的聚并模型如表1所示。表1本文研究的聚并模型详细信息Table 1Details of the coalescence m

13、odel studied in this paper聚并模型简称聚并模型C1LuoC2Mc-Luo新窗口打开|下载CSV2.4.2 破碎模型体积为V0的母气泡破碎成fV和(1-fV)V0两个子气泡的破碎率27为bV0,fVV0=minmaxPd0,fV13d0,(,d0)dmin=11.4(14)式中，为涡旋大小；d0为母气泡直径；fV为破碎体积分数。Luo等10提出的碰撞频率模型为1(,d0)=0.9231-g131+x2x-113d0-530,x=/d0,d0(15)由Han等11提出的模型中气泡与小涡(d0)和大涡(d0)的碰撞频率为2(,d0)=0.9231-g131+C0+x2x-1

14、13d0-53,d0(16)3(,d0)=0.431-g13x143d053min1,1/22xsin24x,d0(17)式中，C0为无量纲振荡比。黏性剪切引起的碰撞频率由Solsvik等14提出的模型计算。4(,d0)=4+d020.5d0+l0.8221-g4,d0(18)由Shi等15提出的破碎模型考虑了在气泡诱导湍流影响下涡流的平均湍流速度和对应于BIT的特征波数/长度。5(,d0)=C4(1-)ni(H)2d03CbguslipV,d0(19)C4=4C3C12,C3=12C(2)2(20)破碎概率为Pd0,fV13d0,=1-0cexp-d=exp-c=ecre,e=L63A22(

15、21)式中，c为破碎的临界无量纲能量；e()为涡的平均动能；ecr()为破碎所需临界涡的动能。本文所研究的破碎模型的详细信息见表2。表2本文研究的破碎模型详细信息Table 2Details of the crushing model studied in this paper破碎模型简称碰撞频率B11B22+3B32+4B45新窗口打开|下载CSV2.5 气泡直径和体积传质系数在采用PBM方法预测气泡尺寸分布时，利用class method(CM)对气泡进行离散，根据直径大小可将气泡分成若干组。界面面积a由气含率g和Sauter平均直径 (d32)求得a=6gd32(22)d32=inidi

16、3inidi2(23)文献提出了很多计算体积传质系数kLa的方法，其中渗透理论28和表面更新理论29最常用。根据渗透理论，传质系数可描述为kL=2DLvl14(24)式中，DL为氧分子扩散率，设为2.01 10-5cm2/s30；为湍流耗散率；vl为液相运动黏度。2.6 数值模拟为了求解上述方程，本文采用耦合多网格求解器的有限体积格式对方程进行了“高分辨率”离散。搅拌槽的计算域划分为旋转叶轮域和静止域，在网格划分时采用较大的全局网格尺寸，并对叶轮旋转区域进行局部加密，保证计算精度。本文采用四个网格尺寸，分别为501244、757201、1001872、1249354个，来进行网格无关性验证。体

17、积传质系数kLa及平均气泡直径d32是表征气液流动混合的重要参数。采用4种密度的网格得到的体积传质系数kLa及平均气泡直径d32数值模拟结果如图2所示，由图可知，当网格数量增加到100万个左右时，kLa和d32的数值模拟结果基本不变。因此，在随后的模拟中，均采用100万个网格尺寸作为计算基准。图2图2网格数量对模拟kLa和d32的影响Fig.2Effects of mesh number on simulatedkLaandd32液面设置为Degassing边界条件，搅拌桨设置为无滑移壁面，搅拌釜的外壁和挡板设置为静止壁面，初始气泡尺寸采用文献31中的方法计算。采用多重参考系法（MRF）进行稳

18、态计算，气泡的聚并和破碎模型采用udf编写，动量方程、体积分数以及湍流方程均采用一阶迎风格式离散，平衡状态为收敛残差低于10-5且进出口物料守恒。3 实验结果与讨论3.1 聚并和破碎模型对气泡尺寸模拟的影响图3（a）分别为0、200、300、400 r/min转速下拍摄的气泡图片。在同一位置进行多次拍摄，然后通过软件处理得到气泡尺寸分布。由图可知，随着搅拌转速的增加，气泡的尺寸明显减小且形状趋于球形。图3（b）为不同模型组合的平均气泡尺寸和实验测量值以及搅拌桨的剪切速率。由图可知，加入修正因子后，气泡尺寸会在一定程度上减小，相比其他破碎模型，B2和B3破碎模型预测的气泡尺寸更小，原因是这两种模

19、型不仅考虑了传统的小涡引起的气泡碰撞，还考虑了大尺度湍流涡以及黏流剪切引起的气泡碰撞。在八组模型组合中，只有C2-B3（基于介尺度修正并考虑流体剪切的模型）最为接近实验值，而其他组合预测值的误差相对较大；在标准Luo模型中，随着转速的增加，搅拌桨剪切速率上升，对气泡的切割作用越显著，使得气泡直径越小，这与实际情况相吻合。四种诱导气泡碰撞的因素中，湍流涡随机运动诱导碰撞的影响最为显著，其他气泡组合模型无法准确地描述这种气泡行为，使得气泡聚并效率要高于实际值。本文C2-B3将气液体系内介尺度能量的描述和PBM模型相结合，在气液体系稳定的状态下，介尺度能量趋于最小，与此同时，不同尺寸气泡的破碎和聚并

20、状态达到平衡，故该模型能准确描述湍流涡随机运动诱导碰撞的影响，从而与实验结果吻合。图3图30、200、300、400 r/min的气泡图(a); 不同模型不同转速下的平均直径(b)Fig.3Bubble diagrams at 0, 200, 300, and 400 r/min (a); Mean diameter of different models at different speeds (b)图4给出了最优模型组合下200、300、400 r/min转速的气泡平均尺寸分布云图。由图可得，转速为200 r/min时，搅拌桨能够在一定程度上使其周围气泡破碎，但无法切割底部的大气泡；转速为

21、300 r/min时，大气泡数量明显减少，部分区域气泡尺寸更加均匀；转速为400 r/min时，整个生物反应器内的气泡大小都比较相近，且气泡分布最为均匀。图4图4最优模型组合下不同转速的气泡直径分布Fig.4Distribution of bubble diameters at different rotational speeds under the optimal model combination3.2 聚并和破碎模型对模拟体积传质系数的影响图5（a）为生物反应器的溶氧率（DO值）随时间的变化曲线，体系的通气量为1 m3/(m3min)，搅拌转速为200、300、400 r/min。由图

22、可知，随着氮气的通入，溶氧不断被吹出，导致体系内DO值减小；当DO值到达最低点后改通空气，体系的溶氧逐渐增加直至饱和。此外，随着转速增加，体系内氧浓度的下降和上升都更加迅速。图5图5DO值随时间的变化曲线(a)和拟合曲线(b)，不同模型组合在不同转速下的kLa(c)Fig.5Variation curve (a) and fitting curve (b) of DO value with time，kLaof different model combinations at different speeds (c)实验数据如表3所示。在转速相同的情况下，三次实验得到的kLa都很接近，经过计算可

23、得误差范围都控制在3%左右，这说明了实验测量的体积传质系数是可靠的。表3不同实验次数下kLa随转速的变化情况Table 3kLachanges with speed under different number of tests试验次数传质系数kLa200 r/min300 r/min400 r/min第一次0.01550.02420.0384第二次0.01530.02370.0378第三次0.01580.02440.0387新窗口打开|下载CSV图5（c）为不同模型组合模拟的体积传质系数。由图可知，八种模型组合预测的体积传质系数均低于实验值，其中组合C1-B1即文献中常用的Luo聚并和破碎模

24、型7与实验值偏差最大，组合C2-B3模型与实验值最为接近，这与图3（b）气泡尺寸分布结果一致，因为根据式（22）、式（24），体积传质系数与气泡尺寸成反比。值得注意的是在相同破碎模型下，加入基于介尺度的修正因子能够降低聚并效率从而提高体积传质系数，这种修正效果在5 L生物反应器中并不明显，这是因为5 L生物反应器的表观气速较低，由式（13）可知较低的表观气速对应接近于1的修正因子。但是当该模型应用于大型通气搅拌式工业生物反应器时，其表观气速较大，这种修正效果将变得极为明显。此外，考虑黏流剪切的破碎模型将更适合于生物反应器的模拟，这是因为高速搅拌会造成极强的黏流剪切作用。图6为最优模型组合下体积

25、传质系数分布云图，由图可知，kLa值大致呈现由中心向外壁逐渐降低的趋势，除此之外，一方面提高转速造成体系中气含率增大，kL增大；另一方面破碎气泡的作用也使得气泡尺寸变小，相际接触面积a增大，两方面作用下使得kLa值与转速成正相关。图6图6最优模型组合下不同转速的体积传质系数分布Fig. 6Distribution of volumetric mass transfer coefficients at different rotational speeds under the optimal model combination3.3 400 m3通气搅拌式工业生物反应器桨型的优化在工业化生产中，

26、为减少建筑成本和人工成本，工厂通常使用大型通气搅拌式生物反应器进行生物发酵，其中搅拌桨型对生物反应器效能而言至关重要，合适的搅拌桨型组合能使反应器内氧气的分布更加均匀，搅拌死区体积分数更小，从而缩短好氧生物的发酵时间，提高发酵效率。本文在设计大型生物反应器的罐体时，保持其与5 L反应器几何相似，并微调高径比进行系列放大。400 m3生物反应器模型内壁有四块挡板，挡板高11.45 m，宽0.45 m，厚度为6 cm。生物反应器内径T=5.54 m，液面高度为12.5 m，封头曲面高度1.4 m。采用四层涡轮桨，桨叶直径D=2 m,叶片长度41 cm，叶片厚度2 cm，桨叶间距2.88 m，下层桨

27、离底高度1.6 m，气体分布器是内径为1.56 m、宽度为33 cm的圆环，圆环上方均匀分布着5圈气孔，每圈有72个直径为2 cm的圆孔。生物反应器的放大准则分为单位体积功率相等、末端剪切力相等、单位体积传质系数相等。本文通过末端剪切力，即末端线速度相等的放大原则确定400 m3生物反应器的转速为50 r/min。本文考察的轴流桨有六斜叶圆盘搅拌桨，径流桨有非对称式抛物线搅拌桨、布鲁马金式搅拌桨以及六直叶圆盘搅拌桨。将这些桨分为图7所示的四种组合组合(A), 非对称抛物线式-六斜叶圆盘式-六斜叶圆盘式-非对称抛物线式搅拌桨; 组合(B), 非对称抛物线式-六斜叶圆盘式-六斜叶圆盘式-布鲁马金式

28、搅拌桨; 组合(C), 六直叶圆盘式-六斜叶圆盘式-六斜叶圆盘式-布鲁马金式搅拌桨; 组合(D), 布鲁马金式-六斜叶圆盘式-六斜叶圆盘式-布鲁马金式搅拌桨，综合对比生物反应器组合桨相应的气含率及kLa，可选择最优桨型。图7图7400 m3生物反应器桨型组合Fig.7400 m3bioreactor paddle type combination如图8所示，在相同工况下(50 r/min, 0. 3 m3/(m3min)组合(A)、(D)对应的生物反应器底部存在近1/3的死区，导致有效发酵空间减少；组合(C)对应的生物反应器虽然死区最小，但传质效果差，气泡分散效果不佳且生物反应器顶部气含率不高；组合(B)对应生物反应器的死区不仅较小，而且对气泡分散效果要优于其他三种组合。如图9所示，组合（A）、（D）对应生物反应器的传质主要集中在搅拌桨周围，原因是该搅拌桨组合并不能