基于介尺度PBM 模型的生物反應器放大模擬及實驗研究

萬景，張霖，樊亞超，劉勰民，駱培成，張鋒，張志炳

（1 南京大學化學化工學院，江蘇 南京 210023； 2 中國石化大連石油化工研究院，遼寧 大連 116045；3東南大學化學化工學院，江蘇 南京 211189）

引 言

進入21 世紀，生物技術因資源可再生、環境友好等特點受到人們越來越多的關注。工業生物技術已成為世界各國爭相發展的重點項目，其中生物反應器更是占據工業生物技術的中心位置。生物反應周期長、連續操作、效率低的特點導致了生物反應器具有高能耗、低設備利用率的缺點，因此高效節能的反應器成為科技與工程研發核心。由于生物反應器的設計和放大過程中存在著諸多困難[1-2]，大型生物反應器在設計初期很難通過實驗的方式提供可靠的設計依據。而計算流體力學（CFD）技術可有效提供反應器內的流場、傳質、傳熱等信息，因此利用CFD 相關軟件進行數值模擬已成為大型生物反應器設計的重要途徑[3-4]。

耦合CFD 與群體平衡模型(PBM)，涵蓋了實際存在的氣泡聚并、破碎等氣泡動力學行為，從而可準確預測生物反應器內流場分布以及傳質能力。隨著對氣泡動力學行為的深入認識，人們開發了多種氣泡聚并與破裂的模型[5-6]。在氣泡聚并模型方面，Luo 模型[7]被廣泛應用于鼓泡塔模擬，但研究表明該模型會過高估計聚并頻率，因此楊寧等[8-9]將氣液體系內介尺度能量的描述和PBM 模型相結合，在氣液體系穩定的狀態下，介尺度能量趨于最小，與此同時，不同尺寸氣泡的破碎和聚并狀態達到平衡，故修正之后的模型能夠更加準確地估計聚并頻率。在氣泡破碎模型方面，Luo 等[10]將能譜與渦能量相結合，他們認為渦不僅具有尺度，還具有能量，由此提出了湍流渦誘導碰撞氣泡破碎的模型。Han 等[11]基于表面能量密度的增加提出了一種新型破碎依據，并將表面振蕩和大尺度渦引起的破碎影響納入模型；此外，Han 等[12-13]還改進了數學模型，使其能夠適用于湍流的全能譜和氣泡的多次破碎。Solsvik 等[14]的模型考慮了黏流剪切對氣泡破碎的影響。Shi 等[15]提出了氣泡誘導湍流的修正破碎模型，考慮了氣泡誘導湍流影響下渦旋的平均湍流速度和對應的特征波數/長度尺度。

以上模型為模擬氣泡聚并破裂的動力學行為提供了充分的理論基礎，但是生物反應器內存在著湍流渦和氣泡群介尺度結構，如何選擇合適的模型闡明這種介尺度結構仍然面臨巨大的挑戰。本課題組在5 L 通氣攪拌式工業生物反應器中通過組合不同聚并、破碎模型的方式，對八種模型組進行考察，分別從體積傳質系數kLa和氣泡尺寸兩個方面與實驗值對比，這八種組合模型包括以下四種誘導氣泡碰撞的因素，它們分別為湍流渦隨機運動誘導碰撞、浮升力差異誘導碰撞、大氣泡的尾流夾帶誘導碰撞和速度梯度差誘導碰撞。結果顯示，與實驗數據最為吻合的是基于介尺度理念的修正聚并模型與考慮黏流剪切的破碎模型組合。整個計算過程在商業軟件Ansys Fluent 17.0 平臺上進行。由于小型生物反應器的工業應用前景不大，故本文在選定最優氣泡聚并破碎模型的基礎上，通過葉輪末端剪切力相等的放大原則將5 L 生物反應器放大到400 m3。攪拌是化工生產過程中的重要環節，開發新型攪拌槳一直都是化學工程領域的熱點[16-18]，本文通過軟件模擬新開發了四種組合槳，考察不同槳型對生物反應器內氣泡分散的效果，在分析了各個組合槳的流場特性、混合死區分布等特征之后，選出400 m3生物反應器的最優組合槳，為大型通氣攪拌式工業生物反應器的設計與研究提供參考。本文組合氣泡不同的聚并效率和破碎頻率，得到了貼合實驗值的氣泡聚并破碎模型；并通過放大準則得到400 m3通氣攪拌式工業生物反應器相關的結構和工況參數，進而使用最優聚并破碎模型模擬不同槳型的氣泡分散情況。

1 實驗裝置與方法

1.1 實驗裝置

實驗裝置如圖1 所示，該生物反應器外壁是具有橢圓封頭的圓柱形玻璃材質，內壁有四塊不銹鋼擋板，擋板高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 生物反應器實驗裝置Fig.1 Bioreactor experimental setup

1.2 實驗方法

由于反應器內氣相氧氣濃度恒定，且假設混合過程是理想的，所以體積傳質系數kLa的計算式[19]為

式中，CI0為初始溶解氧濃度；τp為DO 電極時間常數；CS為飽和溶解氧濃度；CI為實時溶解氧濃度。如果反應器的時間常數τR=1/kLa，遠大于τp。則式（1）可簡化為

因為飽和溶解氧濃度CS和初始溶解氧濃度CI0在整個實驗過程中不會變化，故只需記錄不同時刻t的實時溶解氧濃度CI，便可繪制函數-ln[(CS-CI)/(CSCI0)]的圖像。由式(3)可知，kLa可通過函數-ln[(CSCI)/(CS-CI0)]對時間t求導得到。

實驗開始前，先注純水到生物反應器指定液位，校正溶氧值后開啟攪拌（轉速可調節200、300、400 r/min），打開連接氮氣鋼瓶的氣體流量計，調節流量為4 L/min，通過溶氧電極記錄體系中溶氧數據的變化；體系內的溶氧量會逐漸降低，待平衡后，停止通入氮氣，打開氧氣閥通入氧氣，繼續監測溶氧值的變化，每組轉速重復3次實驗，以此減少實驗設備對結果的影響。通過高速攝像機(Phanton v2640)記錄氣泡信息，采集速度為1000 幀/秒，采集時間為30 s；得到氣泡圖像后，進一步處理得到氣泡尺寸分布。

2 數學模型和數值模擬

2.1 控制方程

采用Eulerian-Eulerian 雙流體方法進行三維數值模擬[20-21]。該方法將氣相和液相分別視為分散相和連續相。兩相滯留量滿足相容條件

各相的連續性平衡方程為

各相的動量平衡方程為

2.2 相間作用力

根據Sanyal 等[22-23]的研究結果，本文模擬只考慮曳力，忽略其他相間作用力，曳力系數選擇Schiller-Naumann模型[24]模擬氣液體系相間作用力。

其中，曳力系數為

式中，CD是氣泡群曳力系數;CD,0為單氣泡曳力系數，CD=CD,0(1-αg)p，本模擬中p=1。

2.3 湍流模型

本文采用RNGk-ε湍流模型[25]來預測氣泡塔內流體的流動特性。方程為

其中

式中，Пk,l和Пe,l表示考慮氣泡誘導湍流的源項。

2.4 群平衡模型

2.4.1 聚并模型 氣泡聚并速率通常采用Luo 模型[7]，為

湍流波動引起的氣泡碰撞頻率[26](Prince and Bl)為

Luo[7]提出的聚并效率模型為

楊寧等[8-9]將氣液體系內介尺度能量的描述和PBM 模型相結合，氣液體系在穩定狀態下，介尺度能量趨于最小，同時不同直徑之間的氣泡破碎和聚并處于動態平衡，從而提出了對聚并模型的修正方法。其中Luo 聚并模型[7]對應修正因子關系式見式（13）。

本文所研究的聚并模型如表1所示。

表1 本文研究的聚并模型詳細信息Table 1 Details of the coalescence model studied in this paper

2.4.2 破碎模型 體積為V0的母氣泡破碎成fV和(1-fV)V0兩個子氣泡的破碎率[27]為

式中，λ為渦旋大小；d0為母氣泡直徑；fV為破碎體積分數。

Luo等[10]提出的碰撞頻率模型為

由Han 等[11]提出的模型中氣泡與小渦(λ≤d0)和大渦(λ>d0)的碰撞頻率為

式中，C0為無量綱振蕩比。

黏性剪切引起的碰撞頻率由Solsvik等[14]提出的模型計算。

由Shi 等[15]提出的破碎模型考慮了在氣泡誘導湍流影響下渦流的平均湍流速度和對應于BIT 的特征波數/長度。

破碎概率為

式中，χc為破碎的臨界無量綱能量；eˉ(λ)為渦的平均動能；ecr(λ)為破碎所需臨界渦的動能。

本文所研究的破碎模型的詳細信息見表2。

表2 本文研究的破碎模型詳細信息Table 2 Details of the crushing model studied in this paper

2.5 氣泡直徑和體積傳質系數

在采用PBM 方法預測氣泡尺寸分布時，利用class method(CM)對氣泡進行離散，根據直徑大小可將氣泡分成若干組。界面面積a由氣含率αg和Sauter平均直徑(d32)求得

文獻提出了很多計算體積傳質系數kLa的方法，其中滲透理論[28]和表面更新理論[29]最常用。根據滲透理論，傳質系數可描述為

式中，DL為氧分子擴散率，設為2.01× 10-5cm2/s[30]；ε為湍流耗散率；vl為液相運動黏度。

2.6 數值模擬

為了求解上述方程，本文采用耦合多網格求解器的有限體積格式對方程進行了“高分辨率”離散。攪拌槽的計算域劃分為旋轉葉輪域和靜止域，在網格劃分時采用較大的全局網格尺寸，并對葉輪旋轉區域進行局部加密，保證計算精度。

本文采用四個網格尺寸，分別為501244、757201、1001872、1249354 個，來進行網格無關性驗證。體積傳質系數kLa及平均氣泡直徑d32是表征氣液流動混合的重要參數。采用4種密度的網格得到的體積傳質系數kLa及平均氣泡直徑d32數值模擬結果如圖2 所示，由圖可知，當網格數量增加到100 萬個左右時，kLa和d32的數值模擬結果基本不變。因此，在隨后的模擬中，均采用100萬個網格尺寸作為計算基準。

圖2 網格數量對模擬kLa和d32的影響Fig.2 Effects of mesh number on simulated kLa and d32

液面設置為Degassing 邊界條件，攪拌槳設置為無滑移壁面，攪拌釜的外壁和擋板設置為靜止壁面，初始氣泡尺寸采用文獻[31]中的方法計算。采用多重參考系法（MRF）進行穩態計算，氣泡的聚并和破碎模型采用udf 編寫，動量方程、體積分數以及湍流方程均采用一階迎風格式離散，平衡狀態為收斂殘差低于10-5且進出口物料守恒。

3 實驗結果與討論

3.1 聚并和破碎模型對氣泡尺寸模擬的影響

圖3（a）分別為0、200、300、400 r/min 轉速下拍攝的氣泡圖片。在同一位置進行多次拍攝，然后通過軟件處理得到氣泡尺寸分布。由圖可知，隨著攪拌轉速的增加，氣泡的尺寸明顯減小且形狀趨于球形。圖3（b）為不同模型組合的平均氣泡尺寸和實驗測量值以及攪拌槳的剪切速率。由圖可知，加入修正因子后，氣泡尺寸會在一定程度上減小，相比其他破碎模型，B2 和B3 破碎模型預測的氣泡尺寸更小，原因是這兩種模型不僅考慮了傳統的小渦引起的氣泡碰撞，還考慮了大尺度湍流渦以及黏流剪切引起的氣泡碰撞。在八組模型組合中，只有C2-B3（基于介尺度修正并考慮流體剪切的模型）最為接近實驗值，而其他組合預測值的誤差相對較大；在標準Luo模型中，隨著轉速的增加，攪拌槳剪切速率上升，對氣泡的切割作用越顯著，使得氣泡直徑越小，這與實際情況相吻合。四種誘導氣泡碰撞的因素中，湍流渦隨機運動誘導碰撞的影響最為顯著，其他氣泡組合模型無法準確地描述這種氣泡行為，使得氣泡聚并效率要高于實際值。本文C2-B3將氣液體系內介尺度能量的描述和PBM 模型相結合，在氣液體系穩定的狀態下，介尺度能量趨于最小，與此同時，不同尺寸氣泡的破碎和聚并狀態達到平衡，故該模型能準確描述湍流渦隨機運動誘導碰撞的影響，從而與實驗結果吻合。

圖3 0、200、300、400 r/min的氣泡圖(a);不同模型不同轉速下的平均直徑(b)Fig.3 Bubble 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 時，攪拌槳能夠在一定程度上使其周圍氣泡破碎，但無法切割底部的大氣泡；轉速為300 r/min 時，大氣泡數量明顯減少，部分區域氣泡尺寸更加均勻；轉速為400 r/min 時，整個生物反應器內的氣泡大小都比較相近，且氣泡分布最為均勻。

圖4 最優模型組合下不同轉速的氣泡直徑分布Fig.4 Distribution of bubble diameters at different rotational speeds under the optimal model combination

3.2 聚并和破碎模型對模擬體積傳質系數的影響

圖5（a）為生物反應器的溶氧率（DO 值）隨時間的變化曲線，體系的通氣量為1 m3/(m3·min)，攪拌轉速為200、300、400 r/min。由圖可知，隨著氮氣的通入，溶氧不斷被吹出，導致體系內DO 值減小；當DO值到達最低點后改通空氣，體系的溶氧逐漸增加直至飽和。此外，隨著轉速增加，體系內氧濃度的下降和上升都更加迅速。

圖5 DO值隨時間的變化曲線(a)和擬合曲線(b)，不同模型組合在不同轉速下的kLa(c)Fig.5 Variation curve(a)and fitting curve(b)of DO value with time，kLa of different model combinations at different speeds(c)

實驗數據如表3 所示。在轉速相同的情況下，三次實驗得到的kLa都很接近，經過計算可得誤差范圍都控制在3%左右，這說明了實驗測量的體積傳質系數是可靠的。

表3 不同實驗次數下kLa隨轉速的變化情況Table 3 kLa changes with speed under different number of tests

圖5（c）為不同模型組合模擬的體積傳質系數。由圖可知，八種模型組合預測的體積傳質系數均低于實驗值，其中組合C1-B1 即文獻中常用的Luo 聚并和破碎模型[7]與實驗值偏差最大，組合C2-B3 模型與實驗值最為接近，這與圖3（b）氣泡尺寸分布結果一致，因為根據式（22）、式（24），體積傳質系數與氣泡尺寸成反比。值得注意的是在相同破碎模型下，加入基于介尺度的修正因子能夠降低聚并效率從而提高體積傳質系數，這種修正效果在5 L 生物反應器中并不明顯，這是因為5 L 生物反應器的表觀氣速較低，由式（13）可知較低的表觀氣速對應接近于1的修正因子。但是當該模型應用于大型通氣攪拌式工業生物反應器時，其表觀氣速較大，這種修正效果將變得極為明顯。此外，考慮黏流剪切的破碎模型將更適合于生物反應器的模擬，這是因為高速攪拌會造成極強的黏流剪切作用。圖6為最優模型組合下體積傳質系數分布云圖，由圖可知，kLa值大致呈現由中心向外壁逐漸降低的趨勢，除此之外，一方面提高轉速造成體系中氣含率增大，kL增大；另一方面破碎氣泡的作用也使得氣泡尺寸變小，相際接觸面積a增大，兩方面作用下使得kLa值與轉速成正相關。

圖6 最優模型組合下不同轉速的體積傳質系數分布Fig.6 Distribution of volumetric mass transfer coefficients at different rotational speeds under the optimal model combination

3.3 400 m3 通氣攪拌式工業生物反應器槳型的優化

在工業化生產中，為減少建筑成本和人工成本，工廠通常使用大型通氣攪拌式生物反應器進行生物發酵，其中攪拌槳型對生物反應器效能而言至關重要，合適的攪拌槳型組合能使反應器內氧氣的分布更加均勻，攪拌死區體積分數更小，從而縮短好氧生物的發酵時間，提高發酵效率。本文在設計大型生物反應器的罐體時，保持其與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，下層槳離底高度1.6 m，氣體分布器是內徑為1.56 m、寬度為33 cm 的圓環，圓環上方均勻分布著5 圈氣孔，每圈有72 個直徑為2 cm 的圓孔。生物反應器的放大準則分為單位體積功率相等、末端剪切力相等、單位體積傳質系數相等。本文通過末端剪切力，即末端線速度相等的放大原則確定400 m3生物反應器的轉速為50 r/min。

本文考察的軸流槳有六斜葉圓盤攪拌槳，徑流槳有非對稱式拋物線攪拌槳、布魯馬金式攪拌槳以及六直葉圓盤攪拌槳。將這些槳分為圖7所示的四種組合[組合(A), 非對稱拋物線式-六斜葉圓盤式-六斜葉圓盤式-非對稱拋物線式攪拌槳;組合(B),非對稱拋物線式-六斜葉圓盤式-六斜葉圓盤式-布魯馬金式攪拌槳;組合(C),六直葉圓盤式-六斜葉圓盤式-六斜葉圓盤式-布魯馬金式攪拌槳;組合(D),布魯馬金式-六斜葉圓盤式-六斜葉圓盤式-布魯馬金式攪拌槳]，綜合對比生物反應器組合槳相應的氣含率及kLa，可選擇最優槳型。

圖7 400 m3生物反應器槳型組合Fig.7 400 m3 bioreactor paddle type combination

如圖8 所示，在相同工況下(50 r/min, 0. 3 m3/(m3·min))組合(A)、(D)對應的生物反應器底部存在近1/3 的死區，導致有效發酵空間減少；組合(C)對應的生物反應器雖然死區最小，但傳質效果差，氣泡分散效果不佳且生物反應器頂部氣含率不高；組合(B)對應生物反應器的死區不僅較小，而且對氣泡分散效果要優于其他三種組合。如圖9 所示，組合（A）、（D）對應生物反應器的傳質主要集中在攪拌槳周圍，原因是該攪拌槳組合并不能很好地徑向分散氣泡，軸流效果明顯大于徑流效果，進而導致氣泡主要集中在攪拌槳附近；組合(B)、(C)對應的生物反應器有較為明顯的環流效應，攪拌槳對氣泡的破碎效果較好，這導致平均體積傳質系數高，且分布較為均勻。綜上所述，組合(B)對應的生物反應器擁有攪拌死區體積分數小、平均體積傳質系數高等特點，故組合(B)對應的攪拌槳為最佳槳型。

圖8 400 m3生物反應器不同槳型氣含率分布Fig.8 Distribution of gas holdup with different paddle types in a 400 m3 bioreactor

圖9 400 m3生物反應器kLa分布Fig.9 kLa distribution of 400 m3 bioreactor

4 結 論

針對生物反應器內存在的復雜湍流渦和氣泡群介尺度結構，本文以5 L 通氣攪拌式工業生物反應器為對象，考察了基于介尺度提出的修正聚并模型以及四種氣泡破碎模型對模擬流體流動行為以及氣液傳質能力的影響；后續通過葉輪末端剪切力相等的放大原則將5 L生物反應器放大到400 m3，在采用最優氣泡聚并破碎模型的基礎上對400 m3通氣攪拌式生物反應器進行了攪拌槳型的優化，結論如下。

（1）氣泡的聚并破碎模型對預測氣泡尺寸分布有很大影響，最優氣泡模型無論從氣泡的數密度還是氣泡所占體積兩方面考察，所得模擬結果均與實驗測量結果最接近。在氣泡聚并模型方面，基于介尺度修正的聚并模型能夠降低體系內氣泡聚并效率，從而模擬的氣泡尺寸和體積傳質系數更加準確。聚并效率降低的程度因反應器尺寸的不同而有所差別，如在小型生物反應器中，加入修正對體系的影響并不顯著，但是在大型生物反應器中，由于表觀氣速較大，這種修正效果將更加明顯；在氣泡破碎模型方面，B3破碎模型除了考慮傳統的小湍流渦引起的氣泡碰撞外，還考慮了黏流剪切引起的氣泡碰撞，使模擬破碎效果與實際氣泡破碎效果十分接近。綜上所述，C2-B3（基于介尺度修正且考慮流體剪切的模型）預測的氣泡尺寸和體積傳質系數最貼合實驗值。

（2）對于400 m3工業生物反應器，非對稱拋物線式-六斜葉圓盤式-六斜葉圓盤式-布魯馬金式攪拌槳的組合形式具有氧傳質能力強、死區體積分數小的特點，適合應用于大型工業生物反應器。