天然氣膨脹制冷工藝仿真優化

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

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

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

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

1 實驗部分

1.1 狀態方程的選擇

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

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

1.2 等熵膨脹數學模型的建立

1.2.1 等熵膨脹過程

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

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

1.2.2 等熵膨脹過程的參數計算

1.2.2.1 密度的計算

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

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

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

式中，下標k 表示迭代次數。

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

1.2.2.2 壓縮因子的計算

天然氣處理過程中需考慮氣體的壓縮性，結合實際情況使計算結果滿足計算精度的要求。基于BWRS 狀態方程，壓縮因子（Z）的計算見式（4）。

密度由式（2）迭代計算所得。混合物密度計算相當重要，后續參數（如逸度、比熱容、泡點、露點等）計算仍需調用密度值。

1.2.2.3 氣化率的計算

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

相平衡方程：

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

物料平衡方程：

Fzi= Vyi+ Lxi（7）

F = V + L （8）

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

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

e = V/F （9）

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

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

1.2.3 等熵膨脹模型

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

部分氣化或部分冷凝的計算一般是已知F，zi，p，T，求e，yi，xi。采用一個經常可獲得快速收斂而對e 的初值無特殊要求的目標函數F（e）。令

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

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

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

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

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

熱平衡方程：

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

Ws= ΔHs= HF- HI（16）

膨脹過程所產生的實際功（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；β 為多變指數，對天然氣取1.35 ～ 1.4；ΔH 為膨脹機前后氣體的實際焓降；ΔHS 為膨脹機絕熱焓降。實際情況下，任何膨脹設備都不可能達到完全絕熱，空分裝置的透平膨脹機的等熵效率一般在80% 左右。

1.2.4 計算流程框圖

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

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

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

2 實例驗證及優化

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

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

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

3 結論

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

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

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

符 號 說 明

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

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

e 氣化率

F 進料量，kmol/h

f 逸度系數

H 焓值，kJ/kmol

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

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

Kii 組分的相平衡常數

k 迭代次數

L 液相的量，kmol/h

p 系統壓力，kPa

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

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

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

T 系統溫度，K

V 氣相的量，kmol/h

W 理想功，kJ

WS實際功，kJ

xi氣相組成，摩爾分數

yi液相組成，摩爾分數

Z 壓縮因子

zi進料組成，摩爾分數

β 多變指數，對天然氣取1.35 ～1.4

ε 預先給定好的足夠小的正數

εp允許誤差

ηS膨脹機絕熱效率

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

上標

L 平衡液相

V 平衡氣相

下標

F 經膨脹機膨脹后

I 經膨脹機膨脹前