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

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

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

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

1 實(shí)驗(yàn)部分

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

可應(yīng)用的狀態(tài)方程有很多種，但對于接近露點(diǎn)、氣體和液體都具良好適用性的真實(shí)流體方程主要有Soave-Redlich-Kwong（SRK），Peng-Robinson（PR），BWRS 等狀態(tài)方程。其中，BWRS 狀態(tài)方程應(yīng)用范圍較廣、精度較高，用于氣態(tài)和液態(tài)模擬的同時亦可用于氣液混合狀態(tài)下的模擬計(jì)算狀態(tài)方程可應(yīng)用于很寬的溫度和密度范圍，首選該方程為天然氣輕烴回收計(jì)算的狀態(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ù)，與各組分的密度、摩爾分?jǐn)?shù)、偏心因子及二元交互系數(shù)相關(guān)［15-16］。

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

1.2.1 等熵膨脹過程

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

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

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

1.2.2.1 密度的計(jì)算

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

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

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

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

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

1.2.2.2 壓縮因子的計(jì)算

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

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

1.2.2.3 氣化率的計(jì)算

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

相平衡方程：

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

物料平衡方程：

Fzi= Vyi+ Lxi（7）

F = V + L （8）

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

氣化率（e）按式（9）計(jì)算。

e = V/F （9）

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

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

1.2.3 等熵膨脹模型

對于實(shí)際混合氣體的熵可按式（13）計(jì)算。

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

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

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

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

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

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

熱平衡方程：

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

Ws= ΔHs= HF- HI（16）

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

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

膨脹機(jī)實(shí)際出口焓：

HF= HI- ΔH （18）

膨脹機(jī)絕熱效率（ηS）：

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

1.2.4 計(jì)算流程框圖

圖1 為等熵膨脹制冷計(jì)算流程。

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

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

2 實(shí)例驗(yàn)證及優(yōu)化

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

為了檢驗(yàn)軟件的可靠性，引用實(shí)例驗(yàn)證。表1為文獻(xiàn)數(shù)據(jù)和軟件計(jì)算數(shù)據(jù)對比。由表1 可知，與文獻(xiàn)數(shù)據(jù)對比，三組數(shù)據(jù)的絕對誤差分別為1.12，0.84，3.11，相對誤差分別為0.005 43，0.003 98，0.019 27，均小于0.02，可認(rèn)為計(jì)算結(jié)果真實(shí)有效，本軟件可用于燃?xì)獾褥嘏蛎涍^程計(jì)算模擬。

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

3 結(jié)論

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

2）與文獻(xiàn)數(shù)據(jù)對比，三組數(shù)據(jù)的絕對誤差分別為1.12，0.84，3.11，相對誤差分別為0.005 43，0.003 98，0.019 27，均小于0.02，可認(rèn)為計(jì)算結(jié)果真實(shí)有效，軟件可用于天燃?xì)獾褥嘏蛎涍^程計(jì)算模擬。

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

符 號 說 明

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

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

e 氣化率

F 進(jìn)料量，kmol/h

f 逸度系數(shù)

H 焓值，kJ/kmol

ΔH 膨脹機(jī)前后氣體的實(shí)際焓降，kJ/kmol

ΔHS膨脹機(jī)絕熱焓降，kJ/kmol

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

k 迭代次數(shù)

L 液相的量，kmol/h

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

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

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

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

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

V 氣相的量，kmol/h

W 理想功，kJ

WS實(shí)際功，kJ

xi氣相組成，摩爾分?jǐn)?shù)

yi液相組成，摩爾分?jǐn)?shù)

Z 壓縮因子

zi進(jìn)料組成，摩爾分?jǐn)?shù)

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

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

εp允許誤差

ηS膨脹機(jī)絕熱效率

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

上標(biāo)

L 平衡液相

V 平衡氣相

下標(biāo)

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

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