張文龍，寧尚雷，靳海波，馬磊，何廣湘，楊索和，郭曉燕，張榮月

（1 北京石油化工學院新材料與化工學院，北京 102617；2燃料清潔化及高效催化減排技術北京市重點實驗室，北京 102617）

引 言

鼓泡塔反應器因操作方便、結構簡單、傳質和傳熱性好等優(yōu)點，被廣泛地應用于石油化工、煤化工、環(huán)境工程和食品工程等領域[1]，如費-托合成、環(huán)己烷氧化、甲醇合成和對苯二甲酸合成等[2]。近年來，我國已經(jīng)成為精對苯二甲酸（PTA）的生產(chǎn)大國，其中對二甲苯（PX）氧化反應器是PTA 生產(chǎn)過程中的關鍵部分，其決定了PTA 的產(chǎn)品質量和生產(chǎn)速率[3-4]。反應物性與操作條件是影響鼓泡塔內流體力學的主要因素，同時其氣含率和氣泡尺寸分布是設計和分析鼓泡塔反應器的重要參數(shù)，因此，研究有機體系下鼓泡塔流體力學特性對反應器的設計與放大具有重要意義。

目前，關于氣液鼓泡塔數(shù)值模擬的研究仍處于發(fā)展階段[5]。同時，雙流體模型中破碎和聚并模型較為復雜、參數(shù)較多，使得鼓泡塔的數(shù)值模擬成為一個充滿機遇與挑戰(zhàn)的研究領域。研究者采用計算流體力學（CFD）方法模擬鼓泡塔反應器內的多尺度流動問題，模擬結果的準確性取決于模型對塔內多相湍流描述是否合理，如相間作用力的設置、氣泡聚并和破碎的描述以及能量耗散機制等。Zhang等[6-7]考察了氣泡聚并、聚并效率、破裂所需的臨界能量、氣泡大小等因素對CFD-BPM 耦合模型數(shù)值計算的影響；Yu等[8]采用CFD-PBM 模擬來預測分散相液滴尺寸分布，并研究脈沖塔內局部液-液流動行為；Yang 等[9-10]考察了壓力對氣泡破碎的影響，并在氣泡破碎模型中加入了密度修正項；Gong 等[11]提出了一種改進的模擬湍流引起氣泡聚并的理論模型用來預測臨界聚并速度以及氣泡大小分布；王鐵鋒等[12-14]的研究認為氣泡群曳力系數(shù)與單個氣泡及氣泡群的運動速度有關，考慮了能量和壓力約束條件，并通過湍流渦體動能和氣泡表面能增量的大小來判斷氣泡破碎的原因；Prince 等[15]的聚并效率模型考慮了液膜排水模型的機理；Luo等[16]的聚并效率模型考慮了能量守恒；Lehr 等[17]的聚并效率模型考慮了臨界速度；Luo 等[18]認為，當湍流渦旋的湍動能大于氣泡表面能量增量時，氣泡發(fā)生破碎，這是一種能量約束的現(xiàn)象；Shi 等[19]結合這些結果，通過引入等效直徑進一步修正了氣泡破碎模型；Hinze[20]研究發(fā)現(xiàn)破碎函數(shù)可以通過慣性力和表面張力相比來定義；Kolmogorov 在Weber 數(shù)的基礎上提出的破碎函數(shù)，認為氣泡破碎行為和Weber 數(shù)有關[21]。因此，上述研究推動了氣泡破碎模型與聚并模型在氣液鼓泡塔數(shù)學模擬的發(fā)展。

本文在課題組前期工作的基礎上[22-26]，在空氣-醋酸體系的鼓泡塔內對曳力模型和聚并模型進行修正，通過CFD-PBM 耦合模型進行了二維和三維數(shù)值模擬，研究了醋酸濃度對鼓泡塔內流體力學參數(shù)的影響，并與實驗數(shù)據(jù)進行對比，驗證該模型的可行性。

1 數(shù)學模型

1.1 雙流體模型

氣液系統(tǒng)采用歐拉-歐拉模型，控制方程如下：質量守恒方程

動量守恒方程

式中，i代表氣相g 或液相l(xiāng)；α、ρ、u、t分別表示相含率、密度、速度矢量和時間；P表示壓力；μeff表示流體有效黏度。

1.2 湍流方程

選擇雷諾時均法（RANS）中的標準k-ε模型[27]，方程具體描述如下：

式中，Cε1=1.44,Cε2=1.92,Cμ=0.09,σk=1.0,σε=1.3。

1.3 相間作用力模型

歐拉-歐拉模型需要相間作用力來封閉，相間作用力的選擇對于模擬結果的精度十分重要。本文模擬選擇了曳力、升力、壁面潤滑力和湍流擴散力，忽略了虛擬質量力的作用。具體表達式為：

1.3.1 曳力 本文基于Li 等[28]的能量最小多尺度（EMMS）模型，在不同質量分數(shù)的醋酸-空氣體系中，通過表面張力函數(shù)f(σl/σ0)修正曳力模型：

簡化并優(yōu)化后的關系式為：

1.3.2 升力 氣泡在液相剪切流中上升時，會產(chǎn)生與氣泡上升方向垂直的一種徑向力，稱為升力。徑向升力對于氣含率的徑向分布至關重要，本文采用文獻[29]的升力，具體表達式為：

1.3.3 壁面潤滑力 壁面潤滑力是由靠近壁面處氣液速度梯度引起的，是使氣泡遠離鼓泡塔壁面的一種力。本文采用文獻[30]的壁面潤滑力，具體表達式如下：

式中，CWD、CWC為無量綱常數(shù)；m=1.5～2。CW為壁面潤滑力系數(shù)，是關于Eo的函數(shù)，具體表達式如下：

1.3.4 湍流擴散力 湍流擴散力是由液相湍流旋渦引起的，此力使得徑向氣含率分布更均勻。本文采用文獻[31]的湍流擴散力，具體表達式為：

式中，kl為湍流強度；CTD為湍流擴散力系數(shù)，取值范圍0.1～1，本文中取值為1。

1.4 群體平衡模型

群體平衡模型（PBM）主要是描述反應器中顆?；驓馀荽笮》植嫉囊环N方法。在鼓泡塔中，氣泡主要發(fā)生破碎和聚并，因此該狀況下群體平衡模型表示如下：

式中，V為母氣泡體積；V′為子氣泡體積；n(V′,t)為體積為V′的氣泡數(shù)密度函數(shù)；c(V-V′,V′)為氣泡聚并速率；β（V,V′）為體積V的氣泡破裂成體積V′的子氣泡分布函數(shù)。

1.4.1 聚并模型 氣泡聚并速率可表示為：

氣泡間碰撞頻率為：

本文基于文獻[16]的聚并模型，在空氣-醋酸體系中，通過引入表觀氣速和表面張力修正項對聚并系數(shù)進行了修正，如表1所示。

表1 聚并系數(shù)的修正Table 1 Correction of the coalescence coefficient

修正的聚并系數(shù)Ce具體表達式如下：

式中，σ0為25℃下標準濃度醋酸的表面張力;σ為實際醋酸溶液的表面張力。

1.4.2 破碎模型 本文主要采用Luo 等[18]的氣泡破碎模型進行數(shù)值模擬，該模型具體表達式為：

式中，β的取值為2.047；ξmin為湍流渦和母氣泡最小尺寸比，即ξmin=λmin/d，其中λmin為Kolmogorov最小湍流渦尺寸的11.4倍；cf為氣泡表面能增加量。

2 實驗裝置和模擬設置

2.1 實驗部分

實驗在玻璃鼓泡塔內進行，實驗裝置和測量方法詳見文獻[32]，塔內直徑D=0.15 m，塔高H=2.2 m，采用差壓變送器、光纖探針、ERT 技術手段測量實驗數(shù)據(jù)。同樣塔底部用以儲存氣體的氣室以及多孔氣體分布板均采用不銹鋼材料制作完成，氣體分布板上共有?2 mm×19的開孔，開孔率約為0.338%。

2.2 物理特性

實驗采用符合標準的工業(yè)冰醋酸（100%的醋酸），并根據(jù)其密度配制不同濃度的醋酸溶液，對所配溶液進行密度、表面張力、黏度的測量，不同濃度醋酸溶液的物性如圖1所示。

圖1 不同質量分數(shù)的醋酸物理性質Fig.1 Physical properties of acetic acid with different mass fractions

2.3 數(shù)值模擬條件的設置

數(shù)值模擬使用ANSYS FLUENT 軟件作為計算平臺，氣相為常溫常壓下的空氣，其作為離散相；液相為不同質量分數(shù)的醋酸溶液，其作為連續(xù)相。對于軟件程序，設置速度為鼓泡塔塔底入口邊界條件，壓力為塔頂出口邊界條件。鼓泡塔塔高為2.2 m，初始靜止液面高度為1.1 m，實驗值在塔高為0.86 m 處進行采集。模擬的時間步長固定為0.002 s，并認為在80 s內達到了準穩(wěn)態(tài)。

3 空氣-醋酸體系二維數(shù)值模擬

3.1 網(wǎng)格無關性驗證

圖2 描述了表觀氣速為0.094 m/s、系統(tǒng)壓力為101325 Pa下，采用網(wǎng)格數(shù)為2200、3080和4590個三種網(wǎng)格驗證徑向氣含率分布的計算精度。從結果看出，徑向氣含率隨著網(wǎng)格數(shù)的增加基本不發(fā)生變化，因此，在綜合考慮模擬的計算精度和計算時間的情況下，本模擬中最終選擇網(wǎng)格數(shù)為3080 個的Grid網(wǎng)格。

圖2 不同網(wǎng)格質量對徑向氣含率的影響Fig.2 Influence of different mesh quality on radial gas holdup

3.2 不同曳力模型對徑向氣含率的影響

在80%醋酸的條件下，將S-N、Tomiyama 單氣泡曳力模型以及本文修正的曳力模型等模擬的徑向氣含率結果進行對比分析，如圖3 所示。單氣泡直徑的設定采用光纖探針實驗測得的平均值，在0.071 和0.094 m/s 的表觀氣速下，測得的氣泡直徑分別為7.0 和6.5 mm。從圖中可以看出，徑向氣含率隨表觀氣速的增加而增加，從塔中心到塔壁處以拋物線形式分布，越接近塔壁處，徑向氣含率越小。而且可以發(fā)現(xiàn)壁面處的徑向氣含率模擬結果明顯低于光纖探針測得的實驗值，這是因為當鼓泡塔的塔徑較小時，具有較強的邊壁效應，導致塔壁面處預測值偏低。

圖3 不同曳力模型和表觀氣速下徑向氣含率對比Fig.3 Comparison of radial gas holdup under different drag force models and superficial gas velocities

通過對三種曳力模型的模擬結果進行對比，結果表明，修正后的曳力模型，直接以CD/db的形式計算曳力的大小，因此不會受到氣泡尺寸大小的影響，與光纖探針測得的實驗值吻合較好；而S-N單氣泡曳力模型的預測結果偏低，Tomiyama 單氣泡曳力模型的預測結果明顯偏高，與王鈺等[33]通過基于EMMS 方法的鼓泡塔反應器CFD 模擬結果相一致。

3.3 軸向高度對氣含率的影響

圖4 為醋酸濃度為100%，塔軸向高度為0.76 m和0.86 m 處的徑向氣含率分布。從圖中可以看出，徑向氣含率和軸向液速呈現(xiàn)出塔中心大，邊壁減小的趨勢。在0.76 m和0.86 m處徑向氣含率分布基本重合，說明在鼓泡塔的這兩個軸向高度處已經(jīng)為充分發(fā)展階段，與文獻[34]中列管型鼓泡塔中流動發(fā)展規(guī)律相一致。

圖4 不同軸向高度下的徑向氣含率分布Fig.4 Radial gas holdup and axial liquid velocity distribution at different column heights

3.4 空氣-醋酸體系徑向氣含率

圖5 為空氣-醋酸體系下徑向氣含率分布。從圖中可以看出，在塔高0.86 m 處，模型的模擬結果與光纖探針測得的實驗值基本吻合。在相同濃度下，徑向氣含率隨表觀氣速的增大而增大。且徑向氣含率從塔中心到塔壁處呈現(xiàn)減小趨勢，在塔壁處徑向氣含率減小幅度較大。

圖5 空氣-醋酸體系下徑向氣含率分布Fig.5 Radial gas holdup distribution with air-acetic acid system

3.5 空氣-醋酸體系平均氣含率

圖6為不同濃度和表觀氣速下醋酸平均氣含率分布。從圖中可以看出，實驗值和模擬結果在10%的誤差范圍內基本吻合。醋酸濃度在50%～100%的范圍內時，平均氣含率先增加后減小，濃度在70%～80%之間時，平均氣含率存在最大值。醋酸濃度在50%～80%范圍內時，表面張力呈現(xiàn)減小的趨勢，表面張力減小使得塔內氣泡穩(wěn)定性減弱，大氣泡破碎成小氣泡，氣含率增加。醋酸濃度在80%～100%范圍內時，表面張力減小，氣泡破碎的概率增大。但同時發(fā)現(xiàn)，在此醋酸濃度范圍內，黏度在1～3 mPa·s的低黏度范圍內是逐漸減小的。Ruzicka等[35]研究認為黏度在0～3 mPa·s 時，氣含率隨黏度的增加而增加。因此，在此濃度范圍內，黏度減小導致氣含率減小。黃娟等[36]通過實驗研究發(fā)現(xiàn)，氣含率在醋酸濃度為60%～80%的范圍內存在最大值，也驗證了本文的模擬結果。

圖6 不同醋酸濃度和表觀氣速下平均氣含率分布Fig.6 Distribution of average gas holdup under different acetic acid concentrations and superficial gas velocities

3.6 空氣-醋酸體系軸向液速

圖7 為不同濃度和表觀氣速下軸向液速分布。由圖可以看出，塔內軸向液速沿塔中心向上運動，且軸向液速從塔中心到塔壁處逐漸減小，在塔壁處軸向液速為負值，說明塔壁面處軸向液速沿壁面向下運動，而且塔中心軸向液速都在0.2～0.4 m/s 范圍內。在相同的醋酸濃度下，軸向液速隨著塔內表觀氣速的增加而增加，這是因為表觀氣速越大，向上運動的氣體對塔中心的液體施加的軸向向上的作用力越大，因此軸向液速越大。在相同的表觀氣速下，軸向液速隨著醋酸濃度的增大基本不變。Yan等[24]研究認為，在不同的表面張力下，塔內的軸向液速基本不變。Krishna 等[37]通過實驗研究發(fā)現(xiàn)，液體黏度對塔內軸向液速的影響作用可以忽略。

圖7 空氣-醋酸體系下軸向液速分布Fig.7 Axial liquid velocity distribution with air-acetic acid system

3.7 空氣-醋酸體系氣泡直徑分布

圖8為不同濃度和表觀氣速下醋酸的氣泡直徑分布。由圖8可看出，在相同的醋酸濃度下，隨著表觀氣速的增加，徑向氣泡直徑分布也隨著增大。這是因為隨著塔內表觀氣速增大，增加了塔內氣泡的聚并，小氣泡數(shù)量減少，大氣泡數(shù)量增多，徑向氣泡直徑分布增大。在相同的表觀氣速下，徑向氣泡直徑分布受到液相醋酸黏度和表面張力的影響，黏度和表面張力的變化導致氣泡發(fā)生聚并和破碎。醋酸濃度在50%～100%的范圍內，隨著黏度增加，徑向氣泡直徑分布相應增加，當醋酸濃度達到80%時，黏度減小，徑向氣泡直徑分布相應減?。欢诖舜姿釢舛确秶鷥?，表面張力在一直減小，徑向氣泡直徑分布也相應減小。因此，醋酸濃度為100%時，其黏度和表面張力都為最小值，此時氣泡直徑分布最小。

圖8 不同醋酸濃度和表觀氣速下徑向氣泡直徑分布Fig.8 Radial bubble diameter distribution under different acetic acid concentrations and superficial gas velocities

3.8 空氣-醋酸體系氣泡數(shù)密度分布

圖9 為不同濃度和表觀氣速下氣泡數(shù)密度分布。由圖9 可看出，氣泡在直徑為7 mm 左右所占有的比例最大，這也與圖8 光纖探針測得的氣泡直徑值分布相符合。在相同表觀氣速下，醋酸濃度為50%～100%范圍內，氣泡數(shù)密度變化程度不大。在醋酸濃度為50%時，氣泡尺寸分布較窄，而醋酸濃度為70%和80%時，氣泡尺寸分布較寬，此時黏度較大。而在醋酸濃度為100%時，表面張力和黏度等物理性質達到最小值，此時氣泡尺寸分布最窄。因為隨著醋酸濃度的增加，液相醋酸黏度先增大后減小，表面張力一直減小。液相黏度增加，氣泡發(fā)生聚并概率增加，小氣泡聚并成大氣泡，氣泡尺寸分布變寬；表面張力減小，大氣泡變得不穩(wěn)定，破碎成小氣泡，塔內小氣泡數(shù)量增多增加了氣泡在塔內的停留時間，氣含率增加，氣泡分布范圍變窄。在相同的醋酸濃度下，表觀氣速的增加使得進入塔內的氣體流量增加，氣泡數(shù)密度增加。

圖9 氣泡數(shù)密度分布（z=0.86 m）Fig.9 Bubble number density distribution(z=0.86 m)

4 醋酸體系三維數(shù)值模擬

4.1 三維網(wǎng)格劃分

本部分為不同濃度醋酸三維數(shù)值模擬，采用三維非結構性網(wǎng)格的劃分，三維網(wǎng)格的網(wǎng)格數(shù)為93856 個。與二維數(shù)值模擬相比，三維數(shù)值模擬可以得到塔截面的氣含率分布狀況、塔內瞬時氣含率分布以及塔內流型的變化情況。因此，進行三維數(shù)值模擬對了解塔內氣液相的瞬時變化情況具有一定的指導意義。

4.2 三維徑向截面氣含率分布

通過ERT 實驗測得的云圖和三維數(shù)值模擬的云圖進行了對比分析，圖10 為50%～80%醋酸的平均氣含率云圖分布。從圖中可以看出，三維數(shù)值模擬的云圖和ERT 測得的云圖變化規(guī)律基本一致。為了可視化不同濃度醋酸在不同表觀氣速下的變化規(guī)律，使用統(tǒng)一的圖例（0～0.40）進行說明。圖中中心處黃色區(qū)域代表較高的氣含率，而邊壁藍色區(qū)域代表較低的氣含率。隨著表觀氣速的增加，圖像中的黃色區(qū)域逐漸增加，藍色區(qū)域逐漸減小。說明在相同醋酸濃度下，氣含率隨著表觀氣速的增加而增加。同時，從圖中可以看出，當醋酸濃度為70%時，塔中心氣泡量相對較大，與黃娟等[36]實驗結果一致。

圖10 不同表觀氣速下不同濃度醋酸橫截面氣含率云圖分布（z=0.86 m）Fig.10 Cross-section contours of gas holdup distribution of different mass fractions of acetic acid under different superficial gas velocities（z=0.86 m）

4.3 三維柱體及軸向截面氣含率分布

從三維柱體圖11可以看出，氣含率隨表觀氣速的增加而增加，且塔中心氣含率較高，塔壁面處氣含率較低。三維柱體圖可以更全面地看出塔中心和塔壁面氣含率沿塔軸向高度的變化趨勢。也可以看出氣泡在上升過程中的流動形態(tài)，即以S 形螺旋式上升。

圖11 80%醋酸三維柱體氣含率分布Fig.11 Three-dimensional cylinder gas holdup distribution map of 80%acetic acid solution

4.4 不同時間三維軸向截面氣含率分布

在表觀氣速為0.094 m/s 的情況下，研究了1～30 s 時間軸向截面的氣含率變化，圖12 為濃度80%的醋酸在不同時間三維軸向截面氣含率云圖分布。從圖中可以看出，隨著模擬時間的增加，初始液面從1.1 m 上升到1 s 時的1.2 m，最后穩(wěn)定在1.4 m。氣泡流動從一開始沿塔高呈現(xiàn)出向上緩沖狀態(tài)，到5 s時垂直貫穿整個液面高度，且呈現(xiàn)出中心對稱分布，達到25 s 后呈現(xiàn)為S 形上升的趨勢。而且從整體上來看，軸向瞬時氣含率分布呈現(xiàn)出塔中心高、塔壁面處低的分布情況，這也與塔內實際的流動情況相符。

圖12 80%醋酸軸向截面氣含率云圖分布Fig.12 Axial-section contours of gas holdup distribution of 80%acetic acid solution

4.5 不同軸向高度處三維徑向截面氣含率分布

圖13 為濃度80%的醋酸三維軸向截面氣含率充分發(fā)展狀態(tài)下的云圖分布。從圖中可以看出，在表觀氣速為0.094 m/s 的條件下，不同軸向高度的氣含率呈現(xiàn)中心高、邊壁低的趨勢，而且較高的氣含率基本分布在塔中心，軸向高度的變化對徑向截面氣含率的大小分布影響不大。

圖13 不同濃度醋酸不同軸向高度處徑向截面氣含率云圖分布Fig.13 Radial-section contour distribution of gas holdup at different axial heights of different mass fractions of acetic acid solution

5 結 論

本文在冷態(tài)空氣-醋酸體系中進行了CFDPBM 耦合模型數(shù)值模擬研究，通過二維和三維數(shù)值模擬方法探究了鼓泡塔內流體力學參數(shù)的變化規(guī)律，具體結論如下。

（1）通過冷態(tài)空氣-醋酸體系下S-N 單氣泡曳力模型、Tomiyama 單氣泡曳力模型和修正的曳力模型等模擬結果的比較，S-N 單氣泡曳力模型徑向氣含率預測值偏低，Tomiyama 單氣泡曳力模型徑向氣含率預測值明顯偏高，而修正的曳力模型具有較好的預測性。

（2）冷態(tài)空氣-醋酸體系二維CFD-PBM 耦合模型數(shù)值模擬，探究了不同濃度醋酸對鼓泡塔流體力學參數(shù)的影響。通過模擬分析和實驗測得的徑向氣含率與徑向氣泡直徑等參數(shù)，發(fā)現(xiàn)模擬結果和實驗值吻合較好，說明修正后的模型具有較好的預測性。通過液體黏度和表面張力對平均氣含率的影響發(fā)現(xiàn)，在50%～100%醋酸濃度范圍內，平均氣含率先增加后減小，當醋酸濃度為70%～80%時，平均氣含率達到最大值。

（3）冷態(tài)空氣-醋酸體系三維CFD-PBM 耦合模型數(shù)值模擬結果表明，模擬測得的徑向截面氣含率云圖分布和ERT 測得的云圖分布結果基本相符，說明修正后的模型具有良好的預測性。不同時間的軸向截面氣含率云圖分布展示出了初始時間段塔內氣含率的變化規(guī)律，而塔高的變化對徑向截面氣含率分布的影響作用很小，且氣泡在塔內呈現(xiàn)S 形螺旋上升狀態(tài)。

符 號 說 明

CD——曳力系數(shù)

Ce——聚并系數(shù)

CL——升力系數(shù)

CTD——湍流擴散力系數(shù)

CW——壁面潤滑力系數(shù)

D——鼓泡塔內徑，m

db——氣泡直徑，mm

Eo——E?tv?s數(shù)

FD——曳力，N

FL——升力，N

FT——湍流擴散力，N

FWL——壁面潤滑力，N

fv——氣泡破碎比

H——鼓泡塔高度，m

H0——靜液面高度，m

Pc(di,dj)——尺寸為di和dj的氣泡間的聚并效率

r/R——徑向位置

T——液體溫度，℃

Ug——表觀氣速，m/s

ul——軸向液速，m/s

α——氣含率

?——湍流耗散率，m2/s3

ζ——氣泡相對直徑

ζmin——氣泡最小相對直徑

μl——液體黏度，Pa?s

μt——湍流黏度，Pa·s

ρl——液體密度，kg/m3

σl——液體表面張力，N/m

?c(di,dj)——尺寸為di和dj的氣泡間的碰撞頻率，m3/s