天然氣膨脹制冷工藝仿真優(yōu)化

唐 慧，劉 云，王泰人，李耀龍

（1. 長江大學(xué) 石油工程學(xué)院 油氣鉆采工程湖北省重點實驗室，湖北 武漢 430100；2. 湖北省天然氣發(fā)展有限公司，湖北 武漢 430223）

在油氣田開發(fā)過程中，為了達到天然氣管輸?shù)臒N露點標準和減少油氣損耗，通常采用輕烴回收工藝來回收乙烷及以上重?zé)N組分。常見的輕烴回收工藝是以冷凝回收法為主的液化加工技術(shù)［1-2］。輕烴回收不僅能提高企業(yè)效益，還能促進天然氣產(chǎn)品的充分利用，通過輕烴回收可對天然氣的熱值進行調(diào)整［3］。天然氣輕烴回收多采用以膨脹機制冷法為主的混合工藝［4-5］。通常采用深冷工藝，經(jīng)透平式膨脹機膨脹制冷，物料降溫程度大，輕烴回收率高［6-7］。目前，我國已經(jīng)在部分油氣田建設(shè)了膨脹制冷輕烴回收裝置，但只有少數(shù)學(xué)者研究了膨脹制冷仿真模型［8-11］。由于等熵膨脹模型需調(diào)用等溫閃蒸、等焓節(jié)流模型，包括很多迭代過程，如果初值不合適，很容易計算發(fā)散，因此需要改進模型的計算方法，解決模型計算的穩(wěn)定性問題，提高模型計算效率。

本工作基于Benedict-Webb-Rubin-Starling（BWRS）狀態(tài)方程，建立了天然氣等熵膨脹制冷仿真模型，主要采用牛頓法、二分法、快速迭代法、弦截法或幾種方法相結(jié)合的方法進行計算。通過實例計算表明，計算過程穩(wěn)定，結(jié)果誤差小，并通過優(yōu)化計算研究了壓力、溫度和膨脹機效率等對輕烴回收率的影響。建立的模型和求解方法可用于天然氣等熵膨脹計算，為輕烴回收系統(tǒng)流程仿真模擬和優(yōu)化奠定了基礎(chǔ)。

1 實驗部分

1.1 狀態(tài)方程的選擇

可應(yīng)用的狀態(tài)方程有很多種，但對于接近露點、氣體和液體都具良好適用性的真實流體方程主要有Soave-Redlich-Kwong（SRK），Peng-Robinson（PR），BWRS 等狀態(tài)方程。其中，BWRS 狀態(tài)方程應(yīng)用范圍較廣、精度較高，用于氣態(tài)和液態(tài)模擬的同時亦可用于氣液混合狀態(tài)下的模擬計算狀態(tài)方程可應(yīng)用于很寬的溫度和密度范圍，首選該方程為天然氣輕烴回收計算的狀態(tài)方程［13-14］，見式（1）。

式中，p 為系統(tǒng)壓力，kPa；ρ 為氣相或液相的密度，kg/m3；T 為系統(tǒng)溫度，K；R 為氣體常數(shù)，8.314 3 kJ/（kmol·K）；A0，B0，C0，D0，E0，γ，a，b，c，d，α 為非固定值參數(shù)，與各組分的密度、摩爾分數(shù)、偏心因子及二元交互系數(shù)相關(guān)［15-16］。

1.2 等熵膨脹數(shù)學(xué)模型的建立

1.2.1 等熵膨脹過程

等熵膨脹過程是天然氣以自身為介質(zhì)通過絕熱膨脹機，在氣體減壓的同時絕熱降溫，自身溫度下降。由于降溫，一部分露點低的組分凝結(jié)成液，以達到純化天然氣的目的。首先選擇物料組分，設(shè)定壓力、溫度、組分摩爾分數(shù)、熵、焓。經(jīng)過等熵膨脹制冷設(shè)備后，根據(jù)膨脹后物料的壓力計算膨脹終溫。不同溫度下可得到不同壓力、熵、焓、氣液比及氣液相各組分含量。

天然氣等熵膨脹過程計算較復(fù)雜，不僅需計算狀態(tài)方程中的參數(shù)，還需計算天然氣物性參數(shù)，如混合密度、逸度、壓縮系數(shù)、比熱容、焓、熵、相平衡常數(shù)、氣化率、泡點露點、偏心因子、臨界偏心因子等。

1.2.2 等熵膨脹過程的參數(shù)計算

1.2.2.1 密度的計算

密度是物性計算中的基本參數(shù)，根據(jù)指定的p，T 和混合物組成（xi，yi）用BWRS 狀態(tài)方程求解氣相密度（ρV）或液相密度（ρL）。二分法、弦截法、牛頓法均可用于天然氣混合密度的計算。牛頓法和弦截法的迭代次數(shù)少、收斂速度快，優(yōu)于二分法，然而運用牛頓法解決問題的過程中，需要同時計算密度函數(shù)和密度函數(shù)的導(dǎo)數(shù)，牛頓法計算函數(shù)值的次數(shù)高于弦截法，所以，無需求導(dǎo)數(shù)的弦截法是計算天然氣混合物密度的一個必然選擇［17］。

以ρ 為變量的BWRS 方程見式（2）。

指定T，p，xi，yi，當(dāng)F（ρ）=0 時，弦截法迭代求ρ，見式（3）。

式中，下標k 表示迭代次數(shù)。

按所設(shè)初值ρV和ρL分別由BWRS 方程求出F（ρ1）和F（ρ2）后，便可由弦截法的迭代公式依次求定下一次迭代用ρ。迭代計算進行至︱ρk+1-ρk︱≤ερ（ερ為允許誤差）為止。當(dāng)取ερ=10-4時，對ρV一般僅需迭代3 ～6 次，而求解ρL時一般則需迭代6 ～18 次。

1.2.2.2 壓縮因子的計算

天然氣處理過程中需考慮氣體的壓縮性，結(jié)合實際情況使計算結(jié)果滿足計算精度的要求?；贐WRS 狀態(tài)方程，壓縮因子（Z）的計算見式（4）。

密度由式（2）迭代計算所得?；旌衔锩芏扔嬎阆喈?dāng)重要，后續(xù)參數(shù)（如逸度、比熱容、泡點、露點等）計算仍需調(diào)用密度值。

1.2.2.3 氣化率的計算

部分氣化或部分冷凝過程的基本方程式見式（5）～式（9）。

相平衡方程：

yi= Kixi1≤i ≤n （5）氣液相組成和歸一化方程：

物料平衡方程：

Fzi= Vyi+ Lxi（7）

F = V + L （8）

式中，F(xiàn) 為進料量，kmol/h；zi為液相或氣相進料組成（摩爾分數(shù)）；V，L 分別為進料經(jīng)加熱或冷卻至溫度T，部分氣化或部分冷凝后進入分離器（壓力為p），分離成平衡的氣、液兩相的量，kmol/h；xi為液相組成（摩爾分數(shù)）；yi為氣相組成（摩爾分數(shù)）；Ki為i 組分的平衡常數(shù)。

氣化率（e）按式（9）計算。

e = V/F （9）

聯(lián)立式（4）～式（7）可得等熵膨脹過程的基本方程式（10）。

為了防止運行故障，本工作采用牛頓法和二分法相結(jié)合的方法求解，在求解非線性方程的同時又提高了運行速度。根據(jù)已求得的ek先用牛頓迭代法求解ek+1，若迭代速度慢或ek+1不在預(yù)定解區(qū)間，則用二分法繼續(xù)運算，直到得到滿足要求的解為止。

1.2.3 等熵膨脹模型

對于實際混合氣體的熵可按式（13）計算。

部分氣化或部分冷凝的計算一般是已知F，zi，p，T，求e，yi，xi。采用一個經(jīng)?？色@得快速收斂而對e 的初值無特殊要求的目標函數(shù)F（e）。令

一般F（e）=0 是非線性方程，可用牛頓迭代法求解，迭代公式見式（12）。

式中，迭代收斂條件為︱F（ek+1）︱≤ε；k 為迭代次數(shù)；ε 為預(yù)先給定好的足夠小的正數(shù)。

在計算過程中應(yīng)先判斷原料混合物在給定T，p 下的狀態(tài)是否處于兩相區(qū)，由于平衡常數(shù)Ki是未知的，需要假設(shè)初值，此計算方法不做贅述，具體參照Antoine 方程來計算。

式中，S 為氣相或液相混合物的實際熵，kJ/（kmol·K）；S0為理想混合氣體的熵，kJ/（kmol·K）。

等熵膨脹過程［18］基本計算方程見式（14）～式（19）。

熱平衡方程：

膨脹過程（等熵條件下）所產(chǎn)生的理想功（WS）：

Ws= ΔHs= HF- HI（16）

膨脹過程所產(chǎn)生的實際功（W）：

W = ΔH = ηsWs= ηsΔHs（17）

膨脹機實際出口焓：

HF= HI- ΔH （18）

膨脹機絕熱效率（ηS）：

式中，SI，SF 分別為膨脹機膨脹前后的熵，kJ/（kmol·K）；SV，SL 分別為膨脹機膨脹后平衡氣相和液相的熵，kJ/（kmol·K）；HI，HF 分別為膨脹機膨脹前后的焓，kJ/kmol；pI，pF 分別為膨脹機膨脹前后的壓力，kPa；TI，TF 分別為膨脹機膨脹前后的溫度，K；β 為多變指數(shù)，對天然氣取1.35 ～ 1.4；ΔH 為膨脹機前后氣體的實際焓降；ΔHS 為膨脹機絕熱焓降。實際情況下，任何膨脹設(shè)備都不可能達到完全絕熱，空分裝置的透平膨脹機的等熵效率一般在80% 左右。

1.2.4 計算流程框圖

圖1 為等熵膨脹制冷計算流程。

圖1 等熵膨脹計算流程Fig.1 Isentropic expansion calculation process.

由圖1 可知，輸入基礎(chǔ)數(shù)據(jù)時，首先應(yīng)選擇物料組分，其中包括選擇烴類物質(zhì)和非烴類物質(zhì)（H2O，H2S，CO2），輸入獨立變量參數(shù)，出口溫度需設(shè)定初值。在以上參數(shù)的基礎(chǔ)上計算初始焓和熵。等熵膨脹算法模擬是在等溫閃蒸算法的基礎(chǔ)上建立起來的，需調(diào)用等溫平衡閃蒸程序［19］計算出口理想的物性參數(shù)。當(dāng)絕對誤差小于0.000 1，滿足等熵膨脹的要求時，軟件繼續(xù)運行計算；否則，以步長為0.5ΔS 的溫變進行試算，直到誤差滿足要求為止。此次循環(huán)計算得出經(jīng)膨脹機后做的實際功、實際出口總焓、實際出口溫度。計算出實際出口溫度后，再次調(diào)用等溫平衡閃蒸程序計算出口實際溫度下的物性參數(shù)。在輸出實際出口總焓前，采用牛頓迭代法對函數(shù)G（T）進行核算，確保實際出口溫度與該溫度下焓的反推溫度一致。在此計算過程中如果f（T）數(shù)值較大，則以f（T）/2 為迭代步長，能減少計算步驟，迅速收斂［20］。

2 實例驗證及優(yōu)化

本工作在上述等熵膨脹的理論和方法基礎(chǔ)上，結(jié)合等溫閃蒸計算程序，以VC 為前臺應(yīng)用程序，Microsoft Access 為后臺數(shù)據(jù)庫，采用ADO 編程接口訪問屬性數(shù)據(jù)，研制了等熵膨脹過程模擬軟件。軟件的功能模塊除了含有物質(zhì)物性參數(shù)模塊和等溫閃蒸模塊外，還補充了等熵膨脹模塊。在等熵膨脹計算模塊輸入物料及組分、初始溫度、初始壓力、膨脹終壓，即可計算膨脹終溫。計算結(jié)果可顯示在表格中，清楚明了。等熵膨脹模型包括多個子程序，如求BWRS 方程11 個參數(shù)的子程序、求氣相或液相密度根的子程序、求焓/熵的子程序、相平衡常數(shù)子程序等。各個子程序的獨立存在使等熵膨脹模型功能模塊化，易于開發(fā)人員理解和改善，通過主程序的管理和子程序的運行，使等熵膨脹模型準確高效地運行。

為了檢驗軟件的可靠性，引用實例驗證。表1為文獻數(shù)據(jù)和軟件計算數(shù)據(jù)對比。由表1 可知，與文獻數(shù)據(jù)對比，三組數(shù)據(jù)的絕對誤差分別為1.12，0.84，3.11，相對誤差分別為0.005 43，0.003 98，0.019 27，均小于0.02，可認為計算結(jié)果真實有效，本軟件可用于燃氣等熵膨脹過程計算模擬。

表1 文獻數(shù)據(jù)和軟件計算數(shù)據(jù)對比Table 1 Comparison between literature data and software calculation data

3 結(jié)論

1）在等溫閃蒸模擬計算的基礎(chǔ)上研發(fā)了等熵膨脹計算軟件。軟件基于BWRS 狀態(tài)方程，利用弦截迭代法計算實際混合氣體密度，以牛頓迭代法計算實際出口溫度，利用二分法快速收斂。多種迭代方法結(jié)合計算可簡化計算過程、提高計算速度。

2）與文獻數(shù)據(jù)對比，三組數(shù)據(jù)的絕對誤差分別為1.12，0.84，3.11，相對誤差分別為0.005 43，0.003 98，0.019 27，均小于0.02，可認為計算結(jié)果真實有效，軟件可用于天燃氣等熵膨脹過程計算模擬。

3）等熵膨脹計算過程繁雜，手工計算或借用Excel 表格計算都不是最好的選擇。軟件隱藏計算程序，操作界面簡明、方法簡單，便于初學(xué)者學(xué)習(xí)，同時可為后續(xù)輕烴回收軟件開發(fā)提供技術(shù)參考。

符 號 說 明

A0，B0，C0，D0，E0，γ，a，b，c，d，α

非固定值參數(shù)，與各組分的密度、摩爾分數(shù)、偏心因子及二元交互系數(shù)相關(guān)

e 氣化率

F 進料量，kmol/h

f 逸度系數(shù)

H 焓值，kJ/kmol

ΔH 膨脹機前后氣體的實際焓降，kJ/kmol

ΔHS膨脹機絕熱焓降，kJ/kmol

Kii 組分的相平衡常數(shù)

k 迭代次數(shù)

L 液相的量，kmol/h

p 系統(tǒng)壓力，kPa

R 氣體常數(shù)， 8.314 3 kJ/(kmol·K)

S 氣相或液相混合物的實際熵，kJ/(kmol·K)

S0理想混合氣體的熵，kJ/(kmol·K)

T 系統(tǒng)溫度，K

V 氣相的量，kmol/h

W 理想功，kJ

WS實際功，kJ

xi氣相組成，摩爾分數(shù)

yi液相組成，摩爾分數(shù)

Z 壓縮因子

zi進料組成，摩爾分數(shù)

β 多變指數(shù)，對天然氣取1.35 ～1.4

ε 預(yù)先給定好的足夠小的正數(shù)

εp允許誤差

ηS膨脹機絕熱效率

ρ 氣相或液相的密度，kg/m3

上標

L 平衡液相

V 平衡氣相

下標

F 經(jīng)膨脹機膨脹后

I 經(jīng)膨脹機膨脹前