曹學文，王春智，孫文娟,3，牟林升

(1.中國石油大學(華東)儲運與建筑工程學院，山東 青島 266580；2.濱州學院 化工與安全學院，山東 濱州 256600；3.中海油研究總院有限責任公司 規(guī)劃研究院，北京 100028)

能源結構低碳化是世界能源發(fā)展的趨勢，與石油和煤炭相比，天然氣更加清潔低碳，中國《能源發(fā)展“十三五”規(guī)劃》提出，到 2020 年，天然氣占一次能源消費量的比重力爭達到10%[1]，能源生產和消費革命將進一步激發(fā)天然氣需求。含硫化氫天然氣是天然氣資源的重要組成部分[2]，目前我國已探明天然氣儲量中有相當部分為含硫天然氣，且有些氣田中硫化氫含量較高。根據天然氣氣質特點選擇高效的脫硫工藝，是高含硫氣田開發(fā)需要解決的重要問題。目前，胺法仍是應用最廣泛的酸性天然氣凈化方法，但其存在設備系統復雜、成本和運行費用高等缺點，尤其是硫化氫含量較高時，脫硫溶劑循環(huán)量大、能耗高[3-5]。為降低含硫天然氣開發(fā)成本，應大力研究發(fā)展能耗低、污染小的新型天然氣脫硫技術。

超聲速旋流分離是一種新興的混合氣體分離技術，最初主要應用于空調空氣中水分的分離，后來荷蘭Shell石油公司和俄羅斯ENGO石油公司將其引入天然氣加工處理領域[6-7]。超聲速旋流分離裝置集膨脹降溫、旋流式氣-液分離、再壓縮等過程于一體，具有結構簡單可靠、無轉動部件、無需化學藥劑、支持無人值守等優(yōu)點[8]，國內外學者針對其在天然氣脫水及重烴方面的應用開展了大量理論及實驗研究[9-10]，在天然氣液化[11-12]、天然氣脫碳[13-14]方面也取得了一定的研究進展，若能將其用于天然氣中硫化氫的脫除，對于完善和發(fā)展天然氣凈化工藝、降低含硫天然氣開發(fā)成本意義重大。

超聲速旋流分離器主要由超聲速噴管、旋流裝置、擴壓段等組成，酸性天然氣進入超聲速旋流分離裝置后，在Laval噴管中高速膨脹至超聲速，溫度降低，當達到一定過飽和狀態(tài)時，天然氣中的硫化氫組分將發(fā)生凝結，形成液滴，同時在旋流場作用下實現氣-液分離。硫化氫氣體在噴管內的自發(fā)凝結是天然氣超聲速旋流分離脫硫技術的前提條件。目前關于高速膨脹過程中氣體自發(fā)凝結特性已開展較多理論、模擬及實驗研究[15-18]，對數學模型的建立、求解及實驗方法等均有探討，對于本文研究的開展具有較好的借鑒意義，但超聲速氣流凝結特性研究多以濕蒸汽、醇類、烴類氣體為介質，目前尚缺乏對于硫化氫氣體凝結特性的研究。因此，筆者結合凝結相變模型，利用計算流體力學軟件FLUENT對天然氣中硫化氫氣體在超聲速噴管內的凝結相變過程進行數值模擬研究，并分析入口壓力、溫度及背壓對凝結流動過程的影響。

1 噴管設計

基于流動特征，Laval噴管可劃分為亞聲速收縮段、達到聲速的喉部和超聲速擴張段[19]。計算喉部尺寸是噴管設計的首要任務，考慮真實氣體效應，選用BWRS方程進行喉部臨界參數的計算，進而確定喉部截面積[20]；采用雙三次曲線法設計噴管收縮段;采用圓弧加直線法設計噴管擴張段。同時考慮噴管邊界層黏性修正，認為邊界層厚度從喉部開始沿軸向線性發(fā)展，修正角取為0.5°[19]。

選取天然氣主要組分甲烷和硫化氫構成二元物系，硫化氫的摩爾分數設為10%。根據上述設計方法，利用MATLAB軟件編制噴管結構設計程序，針對入口壓力5 MPa、入口溫度273 K、入口流量5000 m3/h的工況，對噴管尺寸進行計算，設計所得噴管結構如圖1所示。喉部半徑為5.65 mm，位于x=128.27 mm處。

圖1 Laval噴管結構圖Fig.1 Diagram of Laval nozzle

2 自發(fā)凝結數值模型

2.1 流動控制方程

基于歐拉-歐拉模型，建立氣、液相流動控制方程組。氣相流動控制方程包括連續(xù)性方程、硫化氫氣體質量守恒方程、動量守恒方程及能量守恒方程，如公式(1)～(4)所示，諸方程中添加由于硫化氫氣體自發(fā)凝結而產生的源項。

(1)

(2)

(3)

(4)

由于凝結所形成的液滴粒徑較小(亞微米級)，故忽略氣、液相間速度滑移，認為氣、液兩相具有相同的流動速度。液相控制方程包括液滴數目及液相質量守恒方程，如式(5)～(6)所示，通過用戶自定義標量(User-defined scalar,UDS)輸運方程添加。

(5)

(6)

流動控制方程中的源相表達式如式(7)～(10)所示，公式中m按公式(11)計算。源相利用C語言編寫用戶自定義函數(User-defined function,UDF)程序嵌入到FLUENT中，計算凝結相變對質量、動量以及能量方程的影響。

Sm=-m

(7)

Su=-mu

(8)

Sh=m(hlv-h)

(9)

SY=m

(10)

(11)

2.2 凝結模型

成核是氣體達到一定過飽和度后繼而自發(fā)產生的結果。成核率定義為單位時間、單位體積過飽和氣體中自發(fā)產生的臨界液滴數目。由于凝結相變機理復雜，眾多學者開展了大量研究，提出了許多理論和半經驗模型。目前，Girshick等提出的內洽經典成核理論(Internally consistent classical theory，ICCT)[21-22]因其相對較高的準確性得到了較為廣泛的應用。Lamanna[23]在實驗研究的基礎上，提出利用經驗系數對ICCT進行修正，修正公式預測所得成核率數值與實驗更為相符。因此，筆者利用Lamanna修正后的ICCT模型計算成核率，如公式(12)所示。

(12)

凝結核心形成后，更多氣體分子將在液滴表面凝結，促使凝結液滴繼續(xù)生長?？紤]液滴與周圍氣體的傳熱、傳質，采用Gyarmathy模型[24]計算液滴生長速率，如公式(13)所示。

(13)

2.3 數值計算方法

針對噴管中氣體高速可壓縮流動特征，選用基于密度法的求解器。雷諾應力項的計算采用k-ω湍流模型，該模型適用于墻壁束縛流動計算，且包含了可壓縮性的影響。采用二階迎風格式離散流動控制方程、湍流動能方程和湍流耗散率方程。

噴管入口采用壓力入口，指定總壓、靜壓、總溫、氣體組成及湍流參數；噴管出口采用壓力出口，指定靜壓、回流總溫、湍流參數；壁面采用無滑移、無滲流、絕熱邊界。

噴管內氣流為甲烷和硫化氫混合氣體，氣體狀態(tài)方程選用石油與天然氣工程領域常用的Peng-Robinson (PR)[25]方程，氣體混合物黏度、導熱系數等熱力學參數的計算[26]采用C語言編寫UDF加入到FLUENT軟件中。

模擬過程中，同時監(jiān)測殘差變量及流量的變化情況，當殘差小于10-3，且入口和出口的質量流量(氣、液相質量流量之和)相對誤差小于0.05%時，才認為計算收斂。

2.4 網格無關性驗證

采用CFD專業(yè)軟件ANSYS ICEM CFD進行網格劃分，利用Laval噴管結構的對稱性，計算模型采用二維，考慮到噴管結構的不規(guī)則性，采用非結構化網格。為排除網格劃分情況對數值模擬結果的影響，逐步加密網格，對網格進行無關性驗證。網格數為2408、9216、13938及18541情況下噴管中心軸線處液相質量分數分布如圖2所示(入口壓力5 MPa，入口溫度273 K)。從圖2可以看出，當網格數大于13938時，計算結果趨于穩(wěn)定。因此，綜合考慮計算精度要求及計算機求解效率，在對噴管內酸性天然氣凝結流動特性的數值計算中，網格數量選定為13938，網格劃分情況如圖3所示。

圖2 不同網格數下液相質量分數(Y)沿噴管中心線分布Fig.2 Distributions of liquid-phase mass fraction (Y) along nozzle axis with different mesh cells Tin=273 K;pin=5 MPa

圖3 網格劃分Fig.3 Mesh display

3 結果與討論

3.1 模型驗證

為驗證前述模型的準確性，設計矩形截面可視化噴管及相應的測試系統，開展Laval噴管內氣體凝結相變實驗研究。室內實驗過程中，天然氣存在較大危險，硫化氫是劇毒氣體，實驗室條件下難以采用真實的天然氣組分進行實驗研究，因此只能采用其他實驗介質替代。結合實驗室實際狀況，采用濕空氣為介質進行超聲速流動凝結特性的研究。采用壓力分布測試系統測量噴管內壓力分布，多波長消光法顆粒測試系統測量噴管內凝結液滴數目和粒徑，設計噴管總長為177.3 mm，喉部尺寸為8 mm×8 mm，入口尺寸為40 mm×8 mm，出口尺寸為12.51 mm×8 mm。實驗工況為：入口壓力0.574 MPa，入口溫度292.95 K，入口濕度90.3%。

實驗測得噴管沿程壓力、液滴數量密度、液滴半徑數據與模擬所得結果對比如圖4所示。從圖4(a)中可以看出，實驗測得噴管沿程壓力變化趨勢與數值模擬結果基本一致，隨著氣體在噴管內的流動，壓力逐漸降低，由于凝結液量非常小，凝結時釋放的潛熱對周圍氣體的加熱作用被高速膨脹引起的溫降效應所掩蓋，并未捕捉到明顯的凝結沖波現象。從圖4(b)中可以看出，實驗測得凝結液滴數目分布與模擬結果呈現出相同規(guī)律，液滴數目達到一定數量級后基本不變。從圖4(c)中可以看出，實驗測得凝結液滴粒徑與數值模擬結果均在亞微米級。但由于噴管擴張段內氣體流速達到超聲速，測試裝置對于液滴的捕捉存在一定困難。另外，雖然實驗系統中采用了兩級過濾器進行固體、液體顆粒的過濾，但仍無法完全消除外界顆粒對氣體凝結的影響，使得實驗測得數據與模擬結果之間存在一定偏差。以lg(Nexp/Nsim)衡量液滴數目之間的偏差在-1.78～1.01個數量級，與目前國內外常用氣體凝結參數測試結果的偏差相當[17,27]，且模擬所得凝結起始位置與實驗結果較為吻合，進一步說明了所建氣體自發(fā)凝結數值模型能較為準確地描述噴管內氣體的超聲速凝結流動過程。

圖4 噴管內凝結流動參數實驗數據與模擬結果對比Fig.4 Comparison of flow and condensation parameters in nozzle between experimental data and simulation resultsTin=292.95 K;pin=0.574 MPa(a)p；(b)N；(c)r

3.2 凝結參數分布

基于上述凝結流動數學模型，模擬計算甲烷-硫化氫混合氣體在噴管中的凝結流動過程，入口壓力5 MPa，入口溫度273 K工況條件下，模擬所得凝結參數沿噴管軸線分布如圖5所示。

圖5 凝結參數沿噴管中心線分布Fig.5 Distributions of condensation parameters along nozzle axisTin=273 K;pin=5 MPa(a)ΔT，J；(b)N，r；(c)Y，ys

由圖5(a)看到，隨著氣流在噴管內的膨脹，硫化氫氣體由不飽和狀態(tài)逐漸進入過飽和狀態(tài)，過冷度不斷增加，最大可達40 K左右，氣體處于極度熱力學不平衡狀態(tài)。從圖5(a)可以看出，氣體達到飽和狀態(tài)后并不會立即發(fā)生凝結，而是達到一定過冷度后才出現大規(guī)模成核現象，凝結核在較短距離內急劇產生，成核率峰值可達1021數量級，并釋放凝結潛熱，過冷度迅速降低，成核環(huán)境遭到破壞，成核率由峰值劇降為0，成核過程結束。由于成核過程的急劇性，液滴數量密度也急劇升高至1015數量級，成核過程結束后，液滴數目保持不變，如圖5(b)所示，但由于氣體仍處于過冷狀態(tài)，硫化氫氣體在凝結核心表面繼續(xù)凝結，使得液滴繼續(xù)生長，液滴半徑和液相質量分數持續(xù)增加，混合氣體中硫化氫的摩爾分數相應降低，如圖5(c)所示。

在噴管出口處，氣相中硫化氫摩爾分數可由10%降至3.33%，顯示了較好的硫化氫液化脫出能力，但尚未達到商品天然氣質量指標，需考慮將超聲速旋流分離技術與其他脫硫工藝相結合。當天然氣中硫化氫摩爾分數較高時，可將超聲速旋流分離技術作為第一步脫硫化氫措施，將硫化氫摩爾分數降至較低水平，然后再進入胺法等常規(guī)天然氣脫硫裝置進一步處理，使凈化氣達到規(guī)定指標，在滿足產品處理指標的同時降低裝置能耗，使傳統技術與新技術的優(yōu)勢得以充分發(fā)揮。

3.3 入口壓力對硫化氫凝結的影響

保持入口溫度273 K不變，不同入口壓力情況下噴管內過冷度、成核率、液相質量分數分布如圖6 所示。從圖6可以看出，更高入口壓力情況下入口過冷度較大，硫化氫氣體能更早達到大規(guī)模成核所需要的過冷狀態(tài)，凝結成核過程發(fā)生在更靠近噴管喉部處，更多的硫化氫氣體從氣相中凝結出來，生成的液相質量分數增加。

凝結起始位置隨入口壓力變化情況如圖7所示。由圖7看到，隨著入口壓力的降低，凝結起始位置向后移動。在入口溫度273 K條件下，當入口壓力低至2 MPa時，在靠近噴管出口處才開始有凝結核心生成，硫化氫氣體不能在噴管內完成自發(fā)成核過程，噴管冷凝效果幾乎為0，此時超聲速旋流分離裝置將無法實現對天然氣中硫化氫氣體的脫除。

3.4 入口溫度對硫化氫凝結的影響

保持入口壓力5 MPa不變，不同入口溫度條件下噴管內過冷度、成核率及液相質量分數分布如圖8 所示。從圖8可以看出，隨著入口溫度的升高，入口過冷度降低，更晚達到凝結成核所需要的極限過冷狀態(tài)，成核起始點后移，且凝結生成的液相質量分數隨入口溫度的升高而減少。

圖6 不同入口壓力噴管內過冷度(ΔT)、成核率(J)及液相質量分數(Y)沿噴管中心線的分布Fig.6 Distributions of supercooled degree (ΔT),nucleation rate (J)and liquid-phase mass fraction (Y) in nozzle at different inlet pressures Tin=273 Kpin/MPa:(1)3；(2)4;(3)5(a)ΔT；(b)J；(c)Y

凝結起始點隨入口溫度變化情況如圖9所示。由圖9看到，隨著入口溫度的升高，凝結起始點后移，在入口壓力5 MPa條件下，當入口溫度高于313 K時，在靠近噴管出口處才開始有凝結核心生成，硫化氫氣體無法在噴管內完成自發(fā)成核過程，此時噴管冷凝效果幾乎為0，超聲速旋流分離裝置將無法實現對天然氣中硫化氫氣體的脫除。

圖7 凝結起始點位置(X)隨入口壓力(pin)變化Fig.7 Nucleation onset points (X)vs inlet pressures (pin)Tin=273 K

圖8 不同入口溫度(Tin)噴管內過冷度(ΔT)、成核率(J)及液相質量分數(Y)分布Fig.8 Distributions of supercooled degree (ΔT),nucleation rate (J)and liquid-phase mass fraction (Y)in nozzle at different inlet temperatures (Tin)pin=5 MPaTin/K：(1)293；(2)283;(3)273(a)ΔT；(b)J；(c)Y

圖9 凝結起始點位置(X)隨入口溫度(Tin)變化Fig.9 Nucleation onset points (X)vs inlet temperatures (Tin)pin=5 MPa

3.5 背壓對硫化氫凝結的影響

保持入口溫度(273 K)和壓力(5 MPa)不變，不同背壓情況下噴管內壓力、過冷度、成核率和液相質量分數如圖10所示。從圖10可以看出，隨著背壓的增大，噴管內產生激波，并逐步向噴管喉部方向移動，激波引起壓力突躍，破壞氣體的凝結環(huán)境。當背壓不是很大時，激波未到達成核區(qū)間，對成核過程不會產生影響，如圖10(c)中(2)所示，但會造成凝結液滴再次蒸發(fā)，液相質量分數變?yōu)?，如圖10(d)中(2)所示。當激波前移至成核區(qū)間時，成核過程受到影響，如圖10各圖中(1)所示，激波的產生導致壓力、溫度的突升，破壞噴管內的低溫環(huán)境，過冷度迅速降低，成核環(huán)境遭到破壞，使得成核條件不再具備，凝結核心無法繼續(xù)生成，噴管冷凝效果幾乎為0。因此，在超聲速旋流分離裝置運行過程中，應合理控制背壓，使激波不進入噴管，避免激波對凝結過程的破壞。

圖10 不同背壓條件(pb)下噴管內壓力(p)、過冷度(ΔT)、成核率(J)及液相質量分數(Y)分布Fig.10 Distributions of pressure (p),supercooled degree (ΔT),nucleation rate (J)and liquid phase mass fraction (Y)in nozzle at different backpressures (pb)Tin=273 K;pin=5 MPapb/MPa:(1)3.2；(2)2.0;(3)0.5(a)p；(b)ΔT；(c)J；(d)Y

4 結 論

(1)氣流進入噴管后高速膨脹，達到一定過冷度后硫化氫氣體將發(fā)生凝結成核現象，成核過程在較窄區(qū)域內完成，此后液滴繼續(xù)生長，液相質量分數不斷增大，氣相中硫化氫含量隨之降低，在入口壓力5 MPa、入口溫度273 K條件下，噴管出口氣相中硫化氫摩爾分數可由10%降至3.33%，但尚無法滿足商品天然氣要求，需與其他脫硫方法結合才能達到規(guī)定的處理指標。

(2)升高入口壓力或降低入口溫度將使硫化氫氣體更早達到凝結起始條件，從而使成核起始點前移，在噴管出口能達到更大的液相質量分數，因此可通過調節(jié)來流入口壓力和溫度來促進硫化氫氣體的凝結過程；過低的入口壓力或過高的入口溫度將使硫化氫氣體無法在噴管內完成自發(fā)凝結過程，超聲速旋流分離裝置將無法實現天然氣中硫化氫的脫除。

(3)隨著背壓的升高，激波在噴管內產生并逐漸前移，激波的產生會破壞凝結所需冷凝環(huán)境，造成凝結液滴的再蒸發(fā)，為避免激波對硫化氫氣體凝結過程的影響，應合理控制背壓，防止激波進入噴管。

符號說明：

dr/dt——液滴生長速率，m/s；

E——總能，J/kg；

h——氣體總焓，J/kg；

hlv——凝結潛熱，J/kg；

J——成核率，m-3·s-1；

kB——Boltzmann常數，1.3806505×10-23J/K；

keff——有效導熱系數，W/(m·K)；

Kn——Knudsen數；

m——單位時間、單位體積內凝結的液滴質量，kg/(m3·s)；

mv——單個硫化氫分子質量，kg；

N——液滴數目，kg-1；

p——壓力，Pa；

pb——背壓，Pa；

pin——入口壓力，Pa；

Prv——氣體Prandtl數；

r——液滴半徑，m；

rc——臨界液滴半徑，m；

S——過飽和度；

Sh——能量源相，J/(m3·s)；

Sm——質量源相，kg/(m3·s)；

Su——動量源相，kg/(m2·s2)；

SY——濕度源相，kg/(m3·s)；

t——時間，s；

T——氣體溫度，K；

Tin——入口溫度，K；

Ts——飽和溫度，K；

u——速度，m/s；

ui、uj——速度分量，m/s；

u′i、u′j——速度波動，m/s；

x——軸向坐標，m;

xi、xj——軸向與徑向位置坐標，m;

X——凝結起始點位置，m;

y——徑向坐標，m;

ys——硫化氫氣相摩爾分數，%;

Y——液相質量分數；

γ——氣體比熱容比；

δij——Kronecker delta數；

ΔT——過冷度，ΔT=Ts-T，K；

ε——Lamanna給出的修正系數，0.01；

θ——無因次表面張力；

λv——導熱系數，W/(m·K)；

μ——黏度，kg/(m·s)；

ρ——混合相密度，kg/m3；

ρs——硫化氫氣體密度，kg/m3；

ρl——液相密度，kg/m3；

ρv——氣相密度，kg/m3；

σ——表面張力，N/m；

τeff——有效應力張量。