內循環氣升式環流反應器生物降解苯酚廢水過程的計算傳質學模擬研究

張超，劉有智，焦緯洲，張巧玲

（中北大學超重力化工過程山西省重點實驗室，山西太原030051）

引 言

氣升式環流反應器(ALR)是在鼓泡塔反應器基礎上發展起來的一種多相流反應器，由于其具有出色的混合和傳熱性能，目前已在能源化工、環保和生物化工等領域[1?3]得到了廣泛的應用。當前ALR主要有內循環式(ILALR)和外循環式兩種類型，相比于鼓泡塔反應器，ALR 中部設置了導流筒，使液相流動更加規整，氣液相間動量傳遞更加均勻，更適用于對剪應力敏感的生化反應過程，可以為脆弱的微生物提供適宜生長的流體力學環境，避免因流體劇烈湍動和機械攪拌所帶來的物理損傷[4?5]。

苯酚作為一種重要的化工原料，對人體和環境會造成巨大危害，因而我國已將其列入須優先處理的污染物之一。苯酚廢水來源廣泛，焦化、食品、農藥、塑料樹脂和石油化工等行業都是苯酚廢水的主要來源，當水體中苯酚濃度超過0.015 mg/L后，會引發魚類的大量死亡，作為農田灌溉水則會造成農作物的減產。當前微生物降解法已被證明是處理低濃度有機廢水的一種有效方法，具有投資成本低、處理效率高、設備簡單等優勢，并且可將有機污染物降解為對環境和人體無害的CO2和水等，避免了二次污染。而ALR 由于具有較大的氣液接觸面積、良好的傳熱和傳質性能以及溫和的流體力學環境，在生物法處理有機廢水中受到了廣泛的關注和應用[6]。Loh 等[7]設計了一種新型外循環式ALR，并采用固定化惡臭假單胞菌對含酚廢水進行處理，實驗結果表明該方法可處理廢水含酚濃度高達3000 mg/L。Ryan 等[8]采用白根霉菌在ALR 內對苯酚廢水的生物降解過程進行了研究，發現反應器內流體流動環境非常適于菌體的生長，該體系內菌體生長速率和苯酚去除率都較高。Feng 等[9]針對一ALR 內生物降解苯酚廢水過程建立了相應的數學模型，考察了氣液兩相流動特性對生物降解過程的影響。

ALR 中氣相含率、液相速度、氣液相界面積等基本參數對于反應器的設計和放大至關重要[10?11]，而近年來隨著計算機技術的快速發展，計算流體力學(CFD)已經成為研究ALR 等多相流反應器中流體流動特性的一種主要方法。目前針對ALR 等多相流反應器中氣液兩相的模擬主要有兩種模型：歐拉?歐拉(E?E)模型和歐拉?拉格朗日(E?L)模型。在E?E 模型中，氣液兩相被看作是兩種相互滲透的連續流體，流場中的任一位置都同時被兩相所占據，氣液兩相按照自身的規律運動，相間存在著動量、能量以及質量交換，通過求解Navier?Stokes 方程來描述氣液兩相的流動狀態[12?13]。在E?L 模型中，氣相被看作是一群離散的顆粒，液相被視為連續相，離散顆粒與連續相之間存在動量、能量和質量交換，通過對大量顆粒的運動進行追蹤，最終得到反應器內流場信息以及離散顆粒的運動軌跡[14]。E?L模型可清晰地描繪出氣泡之間的相互作用，但由于需要追蹤每一個氣泡的運動軌跡，所以計算量巨大，目前僅適用于低氣含率情況下小尺度反應器的模擬計算[15]。E?E 模型所需計算量較小并且具有較高的精度，是目前多相流反應器模擬中使用最為廣泛的模型。

CFD 已廣泛應用于ALR 內化工和生物降解過程的模擬研究[16]，Huang等[17]采用E?E 模型對一中試規模的ALR 中煤直接液化過程進行了模擬，對反應過程中流體流動、質量傳遞和化學反應之間的相互作用規律進行了研究，模擬得到的氣相含率和液相速度的分布與實驗測量結果吻合較好。Wang 等[18]針對一實驗規模的ALR 中生物降解甲苯廢氣過程建立了三維瞬態CFD 模型，研究了反應器內甲苯和菌體濃度的動態分布。在ALR 內傳質過程的模擬計算中，為了對組分湍流傳質微分方程進行封閉，需要獲取湍流擴散系數Dt的數值。在已有研究中，研究者們普遍采用了經驗的湍流Schmidt 數(Sct)模型對Dt進行估算，而該方法的準確性嚴重依賴Sct的估值，并且也未能考慮組分濃度自身脈動特性對傳質過程的影響，因此該方法計算精度并不能達到要求[19?20]，需采用一個更為準確且嚴格的理論模型對ALR 內湍流傳質過程進行計算。計算傳質學-c2?εc模型是近年來提出的一種用于對化工傳質設備進行準確模擬的湍流擴散模型，該模型通過引入脈動濃度方差-c2和脈動濃度方差耗散率εc兩個模型方程實現了對Dt的理論計算，從而擺脫了傳統需要預估Sct的經驗模型，實現了對傳質過程的嚴格模擬。-c2?εc模型自提出以來已成功應用于精餾塔和氣固流化床中傳質過程的模擬研究[21?22]，然而其能否適用于氣液兩相相互作用強烈的ALR，目前仍有待進一步研究。

本文將建立三維瞬態CFD 模型對ILALR 中生物降解苯酚廢水過程進行模擬研究，采用RNG k-ε模型對流體湍流流動進行計算，引入氣泡群平衡模型(PBM)對反應器內氣泡尺寸分布進行預測，采用計算傳質學模型對湍流傳質過程進行模擬。最后將模擬結果與文獻中的實驗測量值進行對比，以驗證模型的有效性。

1 模擬對象

2007 年，Feng 等[9]報道了利用ILALR 作為生物反應器降解苯酚廢水的過程，反應器結構如圖1 所示。反應器長和寬均為0.2 m，高為1.0 m。導流筒長和寬均為0.14 m，高為0.45 m，安裝在距離反應器底部5 cm 高度處。反應器底部中心處的布氣板尺寸為50 mm×50 mm，布氣孔孔徑為1 mm。反應器內液相為含游離菌體的苯酚無機鹽水溶液，菌體為熱帶假絲酵母菌，該微生物可將苯酚作為主要碳源和能量來源，將其完全分解為二氧化碳和水，液面高度為0.6 m；氣相為空氣。反應器內溫度控制在25℃。在反應器中部設取樣點，用于取樣分析液相中的菌體和苯酚濃度。

圖1 ILALR結構示意圖Fig.1 Configuration of the ILALR

2 數值模型

本文采用E?E 模型來模擬ILALR 中氣液兩相流動，空氣是氣相，含游離菌體的苯酚無機鹽水溶液為液相。建立模型過程中做出如下假設：

(1) ILALR 內初始菌體濃度和初始苯酚濃度均一；

(2)生物降解過程和菌體生長過程僅與液相中菌體和苯酚濃度有關，其他底物，如無機鹽和氧氣等，均為非限制性因素。

2.1 連續性方程

式中，下角標k代表相應的相(k=G 代表氣相，k=L代表液相)，u為速度，α 代表相含率，ρ代表密度。

2.2 動量守恒方程

式中，p為壓力，Mk代表氣相與液相之間的動量交換即相間作用力，μk,eff為有效黏度。

2.3 湍流模型

動量方程中氣液兩相的有效黏度由分子黏度μlam和湍流黏度μt組成，其中μt可表示為：

采用RNG k-ε 模型來計算流體的湍流效應，湍動能k和湍動能耗散率ε方程如下[23]：

式中，Gk是由于流體黏性力而引起的湍流產生項，其計算式為：

2.4 相間作用力模型

氣液相間的動量交換可用相間作用力表示，而作用于氣相或液相上的力大小相等，方向相反：

式中，d32為Sauter氣泡尺寸，通過氣泡群平衡模型進行計算。CD為曳力系數，本文采用Ishii 等[25]提出的模型進行計算：

升力是促使氣泡徑向運動的關鍵作用力，是由于氣泡在液相中受浮力作用上升過程中表面受到不平衡的壓力作用而產生的一種與氣泡運動方向垂直的力：

式中，CL為升力系數，取0.5[26]。

2.5 氣泡群平衡模型

氣泡尺寸是決定氣液相間質量傳遞速率的關鍵參數之一，氣泡群平衡模型(PBM)目前已廣泛應用于氣液兩相流體系中氣泡尺寸分布的預測，其通過系統地對氣泡之間的聚并和氣泡破碎行為進行考察，可從不同尺度對氣液兩相的流動行為進行深入研究[27]。本文采用離散區間法[28?29]對PBM 模型方程進行求解，將直徑在1.2～12.5 mm 范圍內的氣泡分為十組，每組氣泡的守恒方程如下：

式中，fj為第j組氣泡的體積分數，Sj為由于氣泡的聚并和破碎而產生的源項，其表達式為：

其中

式中，vj是第j 組氣泡的氣泡體積，c(vm,vn)和b(vj)分別代表氣泡的聚并和破碎速率，β(v,vm)為子氣泡的分布函數。本文中氣泡的聚并模型采用Luo[30]提出的模型，氣泡的破碎模型采用Lehr 等[31]提出的模型。式(8)中的Sauter 氣泡尺寸d32可由式（15）計算得到：

式中，dj是第j組氣泡的氣泡尺寸。

2.6 局部最大剪應力

對于生物反應器，局部最大剪應力是一個重要參數，因為大多數微生物對剪應力的大小十分敏感，剪應力過高會對微生物細胞形態造成損傷，從而影響生化反應的效率。局部最大剪應力計算公式為[32]：

2.7 組分守恒方程

模型中相關參數數值為[34]：Cc0=0.11，Cc1=1.80，Cc2=2.20，Cc3=0.80，σc=1.0，σεc=1.0。

2.8 生化反應動力學模型

本文采用Jiang 等[35]得到的熱帶假絲酵母菌降解苯酚的本征動力學模型來計算菌體的生長速率rx和苯酚的降解速率rs：

式中，上角標x和s分別代表菌體和苯酚，μmax為最大比生長速率，Ks為苯酚的飽和常數，KI為苯酚的抑制常數，各動力學參數的數值如表1所示。

表1 生化反應動力學參數Table 1 Microbial kinetic constants

2.9 初始條件和邊界條件

初始條件和邊界條件對于模型方程的求解至關重要。ILALR 底部氣相入口設為速度入口條件，入口速度uin根據實驗條件確定，湍動能k 初始值和湍動能耗散率ε 初始值根據Liu 等[36]采用的方法進行計算，脈動濃度方差和脈動濃度方差耗散率εc的初始值采用Li等[21]提出的關系式進行計算：

氣泡在入口處尺寸均一，入口氣泡直徑db采用Chen等[37]提出的關聯式進行計算：

式中，do為氣體分布器孔徑。

液相苯酚初始濃度為1200 mg/L，壁面對于氣液兩相均設為非滑移條件，出口采用脫氣邊界條件[38]。

2.10 網格劃分和數值方法

ILALR 模型網格如圖2 所示，采用結構化網格對計算域進行劃分，通過對網格獨立性進行測試，最終采用網格數目為19665 個的網格進行模擬計算。

采用商業軟件Fluent 15.0 的三維非穩態求解器，微分方程組采用有限體積法進行求解，壓力速度耦合問題采用SIMPLE 算法，為了保證計算精度而采用了二階數值方法。為了提升計算過程的穩定性，使用Krishna 等[39]推薦的逐級遞增時間步長策略：計算第一個100 時間步長采用0.001 s；第二個100 時間步長為0.002 s；第三個100 時間步長為0.005 s；第四個100時間步長為0.01 s；第五個100時間步長為0.02 s；其余時間步長采用0.05 s。

圖2 模型網格示意圖Fig.2 Schematic view of the grid used in simulation

3 模擬結果與討論

3.1 氣相含率和液相速度的分布

圖3 和圖4 給出了不同表觀氣速下ILALR 中高度為400 mm 處模型計算得到的局部時均氣相含率和軸向液相速度分布與Feng 等[9]實驗測量值的對比。由圖可以看出，ILALR 內局部時均氣相含率和軸向液相速度隨著表觀氣速的增大而增大，并且計算得到結果與實驗測量值吻合較好。

圖3 不同表觀氣速下局部時均氣相含率的分布Fig.3 The profiles of local time?averaged gas holdup at different superficial gas velocities

圖4 不同表觀氣速下局部時均軸向液相速度的分布Fig.4 The profiles of local time?averaged axial liquid velocity at different superficial gas velocities

3.2 Sauter 氣泡尺寸和液相局部最大剪應力的分布

圖5 為采用氣泡群平衡模型計算得到的ILALR內Sauter 氣泡尺寸的瞬態分布云圖，可以看出隨著反應器高度的增加，氣泡尺寸逐步增大，并且由于表觀氣速較低，氣泡未能在反應器內形成環流，當氣泡隨液體流出導流筒進入下降區上部后，氣泡聚并現象更為顯著，所以在下降區上部氣泡尺寸最大。

圖5 Sauter氣泡尺寸分布云圖/mFig.5 Contour of Sauter bubble size distribution/m

現有研究表明液相剪應力對微生物形態會產生顯著影響，大部分微生物對剪應力大小十分敏感，甚至在高剪應力情況下會發生破裂[40?41]。所以剪應力分布是生物反應器設計和操作的一個重要指標。圖6 給出了表觀氣速為0.01 m/s 和0.02 m/s情況下，ILALR 中局部最大剪應力的分布。從圖中可看出反應器中剪應力隨著表觀氣速的增大而增大，剪應力在反應器中的分布較為均勻，但在導流筒上部數值較大。在反應器內靠近壁面的下降區和底部剪應力很小，這是由于下降區和底部內氣含率很低，而氣液相間的相互作用會對剪應力產生較大影響[42]。

圖6 不同表觀氣速下最大局部剪應力分布/(kg/(m·s2))Fig.6 Maximal shear stress distributions in the ILALR at different superficial gas velocities/(kg/(m·s2))

3.3 組分濃度變化

圖7 和圖8 給出了采用計算傳質學-c2?εc模型計算得到的ILALR 中局部液相苯酚濃度和菌體濃度隨時間的變化曲線，并將其與Feng 等[9]的實驗測量值進行了對比。

圖7 液相中菌體濃度的模擬值與實驗值的比較Fig.7 Comparison of simulated cell concentration with experimental measurement

從圖7 和圖8 可以看出，在反應初期菌體濃度增長速率較低，這是由于高濃度的苯酚對于菌體的生長有抑制作用，隨著液相中苯酚的不斷降解，菌體的增長速率逐漸增大。而隨著液相中菌體濃度的升高，苯酚的降解速率也在不斷提升。提升初始菌體的濃度，可以有效克服高濃度苯酚的抑制作用，增強菌體對于苯酚的降解能力，提升苯酚降解速率。從圖中還可看出，與Feng等采用的Sct經驗模型相比模型的計算結果與實驗值更為接近，具有更高的精度。

圖8 液相中苯酚濃度的模擬值與實驗值的比較Fig.8 Comparison of simulated dissolved phenol concentration with experimental measurement

3.4 湍流擴散系數分布

液相苯酚的湍流擴散系數Dt和液相湍流黏度μt的分布云圖如圖9 和圖10 所示。圖11 和圖12 給出了在不同反應器高度Dt和μt在徑向上的分布。

從圖9 和圖10 可以看出，Dt和μt的分布形式相似，但并非完全相同。從圖11 可知，計算得到的Dt隨著反應器高度的提升而增大，在固定高度上，升流區內Dt較大，靠近導流筒壁面區域Dt迅速降低，這是由于壁面處采用了非滑移邊界條件，在近壁區域液相的湍動受到了限制。

圖9 液相苯酚的湍流擴散系數云圖/(m2/s)Fig.9 Contour of turbulent mass diffusivity of dissolved phenol/(m2/s)

圖10 液相湍流黏度云圖/(kg/(m·s))Fig.10 Contour of turbulent viscosity/(kg/(m·s))

在ILALR 內傳質過程的模擬計算中，Dt常采用基于傳遞相似性假設的Sct模型進行估算，假定湍流傳質過程中Dt與μt呈正比。然而Sct的取值卻是經驗性的，針對不同情況，Sct的取值也各不相同[43]。Talvy 等[44]采用Sct=1 模擬了氣升式環流反應器中氧氣的吸收過程，Huang 等[17]則采用了Sct=0.75 模擬了氣升式環流反應器中煤直接液化過程。對比圖11和圖12 可以看出，雖然Dt和μt的分布形式十分相似，但并不是簡單的正比例關系，即計算得到的Sct在反應器中并不是一個常數。因此采用Sct模型來模擬ILALR 中的傳質過程是不嚴謹的。本文采用的計算傳質學-c2?εc模型優勢之一即是湍流擴散系數Dt可通過理論模型計算得到，而無須采用Sct等經驗參數。

圖11 不同ILALR高度液相苯酚湍流擴散系數的分布Fig.11 Profiles of turbulent mass diffusivity of dissolved phenol at different heights

圖12 不同ILALR高度液相湍流黏度分布Fig.12 Profiles of turbulent viscosity at different heights

4 結 論

本文采用計算傳質學方法對ILALR 中生物降解苯酚廢水過程進行了模擬研究，利用RNG k-ε 模型模擬了流體的湍流流動，利用氣泡群平衡模型對反應器內氣泡尺寸分布進行了預測，最后采用近年來提出的計算傳質學模型對湍流擴散系數Dt進行計算，主要結論如下。

(1)ILALR 中液相剪應力分布較為均勻，最大值出現在導流筒的上部，剪應力隨著表觀氣速的增大而增大。

(2) 模擬得到ILALR 中液相苯酚濃度和菌體濃度與實驗測量值符合較好，說明了本文采用的計算傳質學-c2?εc模型能夠適用于ILALR 中傳質過程的模擬。

(3) 計算得到ILALR 中湍流擴散系數Dt和湍流黏度μt之間并不是簡單的正比例關系，因而采用傳統Sct經驗模型計算得到的Dt是不準確的。

符 號 說 明

b(v)——氣泡破碎速率，s?1

c(vm,vn)——氣泡聚并速率，m3/s

Cc0,Cc1,Cc2,Cc3——計算傳質學雙方程模型常數

CD——曳力系數

CL——升力系數

-c2——脈動濃度方差

Dt——湍流擴散系數，m2/s

do——氣體分布器孔徑，m

d32——Sauter氣泡直徑，m

Eo——Eotvos數

fi——氣泡組i的體積分數

g——重力加速度，m/s2

KI——苯酚抑制常數，g/m3

Ks——苯酚飽和常數，g/m3

k——湍動能，m2/s2

p——壓力，Pa

Re——Reynolds數

Sct——湍流Schmidt數

u——速度矢量，m/s

v——氣泡體積，m3

α——相含率

β(v,vm)——子氣泡分布函數

ε——湍動能耗散率，m2/s3

εc——濃度方差耗散率，s?1

μeff——有效黏度，kg/(m·s)

μlam——分子黏度，kg/(m·s)

μmax——最大比生長速率，h?1

μt——湍流黏度，kg/(m·s)

ρ——密度，kg/m3

σ——表面張力，N/m

τ——剪應力，kg/(m·s2)