摘 要 本文給出了我?！逗Ｑ笥蜌饧敗氛n程中計算機編程項目的詳細求解方法及參照答案，這些計算機編程項目在教學(xué)中應(yīng)用效果良好，能夠有效的幫助學(xué)生掌握油氣集輸課程的重要的理論，同時發(fā)展學(xué)生解決實際工程問題的能力?，F(xiàn)將實踐教學(xué)中使用的編程項目的解法及參考答案整理成文，方便相關(guān)專業(yè)的廣大師生在本科教學(xué)過程中使用。

關(guān)鍵詞 油氣集輸 編程項目 算例 閃蒸計算 管網(wǎng)計算

中圖分類號：G642文獻標(biāo)識碼：A

在《石油與天然氣工程挑戰(zhàn)性課程教學(xué)實踐——海洋油氣集輸計算機編程項目設(shè)計》一文中，介紹了中國石油大學(xué)（北京）在海洋油氣集輸全英文教學(xué)改革過程中使用的計算機編程項目的內(nèi)容設(shè)計和教學(xué)效果，本文詳細給出了三個計算機編程項目的求解方法和算例，以便開設(shè)類似課程的大專院校使用或在此基礎(chǔ)上對教學(xué)內(nèi)容做進一步改進。

1計算機編程項目1（Computer Project 1， CP1）：采用DAK方程進行天然氣壓縮因子計算

采用DAK方程進行天然氣壓縮因子計算的基本原理請參見文獻[1][2]?，F(xiàn)給出算例如下：

算例1：設(shè)氣體溫度為103.73F（313K），摩爾分?jǐn)?shù)如下表（馮叔初等，2006;例2-3， p.101）。求：壓力為1015.264psia（7.0MPa）時的氣體壓縮因子？

注：表中氣體摩爾分?jǐn)?shù)取自馮叔初等;臨界溫度、臨界壓力取自Ahmed;R = F + 459.67。

【解】：

第1步：計算擬臨界溫度（）、擬臨界壓力（）：;;注意在所有的計算過程中，溫度均需采用絕對溫標(biāo)：熱力學(xué)溫標(biāo)（K）或朗肯溫標(biāo)（R），本文使用朗肯溫標(biāo)。

第2步：如果天然氣組分中不含CO2和H2S，則直接跳至第3步。本例中，CO2和H2S組分的存在會使得天然氣壓縮因子偏離Standing-Katz壓縮因子圖（即DAK方程原始參考數(shù)據(jù)），Wichert和Aziz提出了一種針對酸氣擬臨界參數(shù)（和）的修正方法，使得經(jīng)過修正后的酸氣壓縮因子符合Standing-Katz壓縮因子圖。對于朗肯溫標(biāo)，Wichert和Aziz對擬臨界溫度的修正量可表示為：

（1）

此表達式與溫標(biāo)的選取有關(guān)（對于開氏溫標(biāo)：）。式（1）中代表CO2和H2S在天然氣中總的摩爾分?jǐn)?shù);B代表H2S在天然氣中的摩爾分?jǐn)?shù)。本例中，

第3步：計算擬對比溫度（）、擬對比壓力（）：;。將所得的與代入目標(biāo)函數(shù)（參文獻[1]中式（6）），形成以為唯一未知量的方程。采用牛頓法對其進行數(shù)值求解（參文獻[1]中式（7）），求得（表2給出了詳細的迭代計算結(jié)果）：

表2：DAK狀態(tài)方程求解壓縮因子計算迭代結(jié)果

意味著初始迭代時需給出的初始值，由參文獻[1]中式（4）可得， （取壓縮因子;臨界壓縮因子）。注：表2在計算時，收斂標(biāo)準(zhǔn)設(shè)為

;

第4步：參文獻[1]中式（4），

。

在進行CP1要求的天然氣壓縮因子圖（參文獻[1]，圖1、圖2）計算時，在給定的值下，取的取值范圍為[0，15]即可做出不同下的天然氣壓縮因子曲線;此外，將換算為即可得到壓縮因子與的關(guān)系曲線。

2計算機編程項目2（Computer Project 2， CP2）：二級地面油氣分離系統(tǒng)閃蒸計算

算例2：CP2的主要教學(xué)任務(wù)和教學(xué)目標(biāo)參見文獻[1]。

【解】求解CP2時，需要對每個分離器分別實施閃蒸計算。在完成閃蒸計算程序時，需要采用模塊化的編程思想將任務(wù)分解：（1）編寫一個求解立方形狀態(tài)方程的子函數(shù);（2）編寫一個求解Rachford-Rice目標(biāo)函數(shù)的子函數(shù);（3）編寫一段代碼更新平衡常數(shù)，調(diào)用前兩個子函數(shù)完成相平衡計算。具體求解過程如下：

第1步：求解立方形狀態(tài)方程（編寫Zfactor子函數(shù)）。

流體系統(tǒng)的PVT關(guān)系由狀態(tài)方程表述。Peng-Robinson狀態(tài)方程（PR EOS）和Soave-Redlich-Kwong狀態(tài)方程（SRK EOS）是目前工業(yè)界使用最廣泛的立方形狀態(tài)方程，其表達式分別如下：

（2）

（3）

式（2、3）中，為壓力（psia），為溫度（R），為摩爾容積（ft3/lbmol），為通用氣體常數(shù)10.73 ;和分別為狀態(tài)方程相對于混合物的引力項和斥力項參數(shù)，可由純組分參數(shù)（）和相應(yīng)的混合規(guī)則計算得到：

（4）

（5）

（6）

（7）

式（4-7中），分別為狀態(tài)方程對于第種組分對應(yīng)的純組分物質(zhì)的引力項和斥力項參數(shù);代表混合物中組分的數(shù)量;和是純組分參數(shù)，對于，;對于，;為混合物體系中第種組分的摩爾分?jǐn)?shù);和為混合物體系中第種組分的臨界溫度和臨界壓力;為混合物中第種組分和第種組分間的二元相互作用參數(shù)（Binary interaction parameter，BIP）;為第種組分的對比溫度（）;為與純組分偏心系數(shù)（）有關(guān)的表達式，對于PR EOS和SRK EOS分別如式（8.1）和式（8.2）所示：

（8.1）

（8.2）

在石油和天然氣行業(yè)的應(yīng)用中，立方形狀態(tài)方程常常表示為壓縮因子（z）的形式，根據(jù)實際氣體狀態(tài)方程，式（2、3）可分別改寫為：

（9）

（10）

其中，

;? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（11）

Coats進一步指出式（9）和（10）可以寫成統(tǒng)一的形式：

（12）

式（12）中，對于，;對于。求解式（12）的過程就是求三次方程的根。在教學(xué)實踐中，理解根的選擇是個教學(xué)難點。三次方程必然有三個根，但對于給定的溫度、壓力、組分的體系，只對應(yīng)一個具有物理意義的壓縮因子，需要在求解的三次方程時進行根的選擇（參圖1）。

圖1中，在溫度下，立方型狀態(tài)方程反映的P-v變化關(guān)系如圖中較粗的黑色曲線所示，該曲線在氣液兩相區(qū)內(nèi)分別與泡點線交于A點并與露點線交于E點。AE水平連線對應(yīng)的壓力（即為溫度下的飽和蒸汽壓。A、E點分別對應(yīng)“可能的”液相或氣相摩爾容積和，“可能的”意味著在求解時，我們通常并不預(yù)先知道所給定的混合物是以氣相還是液相存在。圖中C點位于曲線BCD段上，該段曲線上，有，說明在該段曲線上對于閉口系統(tǒng)容積增大壓力也增大，不符合實際情況;因此C點處對應(yīng)的摩爾容積僅僅是一個數(shù)學(xué)意義上的解，沒有實際的物理意義，在求解立方形狀態(tài)方程時應(yīng)當(dāng)舍去。

此時，在剩下的兩個具有物理意義的解和中，必須再舍棄一個，因為在計算過程中針對液相（）或氣相（）組分進行計算最終只能有一個具有物理意義的實根。通常有如下兩種處理方式：

如果預(yù)先知道計算過程針對的是氣相組分（），此時式（4）和式（5）中：，圖1中和兩個根選擇較大的根;對應(yīng)式（12）的解也選擇較大的根，。如果預(yù)先知道計算過程針對液相組分（），此時式（4）和式（5）中：，對應(yīng)式（12）的解也選擇較小的根，即。

在兩個具有物理意義的根（和）中，直接比較二者得到的系統(tǒng)Gibbs自由能大小，因為系統(tǒng)總是傾向于處于最小Gibbs自由能的狀態(tài)。Danesh指出：

（13）

式（13）大于零時，說明使得系統(tǒng)Gibbs自由能最小，取舍;反之，取。

驗證壓縮因子子程序計算的正確性。采用文獻[1]中表3提供的原油組分及物性參數(shù)數(shù)據(jù)，將自行編寫的壓縮因子計算子程序的計算結(jié)果與下表驗證后再進行下一步：

第2步：求解Rachford-Rice目標(biāo)函數(shù)（編寫VLE子函數(shù)）。

由氣液相物料守恒可以推導(dǎo)得出Rachford-Rice目標(biāo)函數(shù)具有如下形式：

（14）

式（14）中，為體系中第種組分在體系中總的摩爾分?jǐn)?shù);為體系達到平衡狀態(tài)時的氣化分率（氣相中物質(zhì)的量與體系中總物質(zhì)的量的比值）;為第種組分的平衡常數(shù)，可由Wilson經(jīng)驗公式近似計算得到：

（15）

式中，，。由式（14-15）可見，在各組分平衡常數(shù)已知的情況下，式（14）中唯一的未知量便是氣化分率，可采用數(shù)值解法（如：Newton-Raphson法）求解方程（14）得到氣化分率：

（16）

當(dāng)時停止迭代（為預(yù)先擬定的收斂標(biāo)準(zhǔn)）。在迭代計算過程中，需要給賦初值使得迭代計算開始（通常?。?，氣化分率的計算值雖然可能超出實際取值區(qū)間[0，1]，只要其最終的計算結(jié)果在區(qū)間[0，1]之間即可，為了在氣化分率超出實際取值區(qū)間時迭代計算仍然能夠繼續(xù)進行，在每一步迭代后都進行如下判斷：

（17）

式（17）中，和分別為組分平衡常數(shù)的最大值和最小值。當(dāng)?shù)諗亢?，氣液相的平衡組分可由式（18-19）求出：

（18）

（19）

此時，我們得到了平衡氣相組分（）以及平衡液相組分（）。通過采用氣相或液相組分替代式（4-5）中的（），我們可以得到氣相或液相的壓縮因子?，F(xiàn)采用文獻[1]中表3提供的原油組分及物性參數(shù)數(shù)據(jù)，讀者可將自行編寫的程序運行結(jié)果與下表對照，驗證求解Rachford-Rice目標(biāo)函數(shù)的正確性;學(xué)生經(jīng)與此表4內(nèi)的提供的數(shù)據(jù)驗證正確后，再進行下一步計算。

第3步：獲得準(zhǔn)確的平衡常數(shù)（編寫MAIN程序）。

命名MAIN主程序，實施平衡常數(shù)更新。第2步計算的前提是平衡常數(shù)已知（由Wilson公式計算）。但是在實際的工程應(yīng)用中，在特定的壓力、溫度和組分條件下，平衡常數(shù)的值往往并不能在一開始計算時就準(zhǔn)確獲得，Wilson經(jīng)驗公式只是較粗略的估算。準(zhǔn)確的平衡常數(shù)需要根據(jù)相平衡理論反復(fù)迭代更新后才能得到。

當(dāng)體系中氣液兩相處于平衡狀態(tài)時，相間各組分逸度相等，可表示為：

（20）

式（20）表示含有種組分的混合物體系處于氣液兩相平衡時，任意第種組分在氣相和液相中的逸度相等（上標(biāo)L表示液相;V表示氣相）。通過逸度的定義（）以及平衡常數(shù)的定義，可得：

（21）

當(dāng)且僅當(dāng)時，式（21）表示體系的平衡常數(shù)，因此可用連續(xù)迭代法（式22）對計算中求得的平衡常數(shù)進行更新：

（22）

在迭代開始時，需要給平衡常數(shù)賦初值，該初值往往由Wilson公式（式15）估算得到。在迭代過程中，當(dāng)滿足如下收斂標(biāo)準(zhǔn)時（式23）停止迭代，認(rèn)為平衡常數(shù)已不再變化。

（23）

由式（22）可以看出，平衡常數(shù)由逸度系數(shù)計算更新得到，而氣液相的逸度系數(shù)（）經(jīng)過嚴(yán)格的熱力學(xué)推導(dǎo)可由狀態(tài)方程得到（式24）：

（24）

式中，

（24.1）

（24.2）

（24.3）

（24.4）

（24.5）

注意在應(yīng)用式（24，24.3，24.5）時，要注意對于氣、液相來說，需分別帶入氣相組分（）或液相組分（）進行計算。圖2給出了CP2閃蒸計算主程序MAIN的算法流程圖。

表5給出了在一定的溫度和壓力條件下，采用文獻[1]中表3提供的原油組分及物性參數(shù)數(shù)據(jù)、依照圖2流程計算獲得的計算結(jié)果，如果學(xué)生自行編程計算的結(jié)果與下表各項數(shù)據(jù)吻合，則說明閃蒸計算流程及程序運行均正確無誤。

第4步：將步驟1至3中獲得的閃蒸計算程序應(yīng)用至2級分離計算（參文獻[1]CP2主要內(nèi)容）中，計算儲罐中原油收率、產(chǎn)出液氣油比（GOR）、產(chǎn)出原油API重度等數(shù)據(jù)。

（1）儲罐中原油收率：

假設(shè)1磅摩爾（1lbmol = 453.5924mol）井流進入高壓分離器，由計算知經(jīng)高壓分離器后得到0.5882磅摩爾液相原油;該部分原油再進入中壓分離器，閃蒸分離后得到0.4931磅摩爾原油;該部分原油進入儲罐后，在標(biāo)況下閃蒸分離后最終得到0.4821磅摩爾原油。由此可見，1磅摩爾井流最終在儲罐中得到0.4821磅摩爾原油，則儲罐中原油收率為48.2%。各級分離器閃蒸分離的計算結(jié)果如圖3所示。

（2）產(chǎn)出液氣油比（GOR）：

氣油比可以定義為產(chǎn)出氣的總體積（SCF， standard cubic feet）與儲罐中得到的原油體積（STB， stock tank barrel）之比。即：

（25）

式（25）中表示氣體在標(biāo)況下的體積;為原油在儲罐中的容積，由實際氣體狀態(tài)方程可得：

（26）

其中為原油在儲罐中的摩爾數(shù)（lb-mol）;為原油在儲罐條件下的壓縮因子;和分別為儲罐的溫度（R）和壓力（psia）。式（26）計算出的原油體積單位為SCF，轉(zhuǎn)化為STB需要除以5.615（1STB = 5.615SCF）。此外，假設(shè)產(chǎn)出氣在標(biāo)況（14.7psia， 60F）下可用理想氣體近似，1磅摩爾的氣體具有379.4SCF的體積，得：

（27）

式（27）中，代表經(jīng)各級分離器分離出氣體得總摩爾數(shù)。將式（26-27）帶入式（25）可得：

（28）

結(jié)合圖3中計算結(jié)果，，，帶入式（28），得到（SCF/STB）。

（3）儲罐原油API重度：

從實際氣體狀態(tài)方程，可得：

（29）

（30）

式中，為油水比重;為儲罐中原油的平均分子量;水密度（）取為62.4lbm/SCF，液相壓縮因子由閃蒸計算得到，則本例得到的結(jié)果為，為輕質(zhì)原油。

3計算機編程項目3（Computer Project 3， CP3）：天然氣管網(wǎng)管段流量及節(jié)點壓力計算

本例中管網(wǎng)系統(tǒng)求解時，要求學(xué)生分別使用兩種解法：（1）Q-formulation：以管段流量為未知數(shù)求解;（2）P-formulation：以節(jié)點壓力為未知數(shù)求解。如圖4所示（假設(shè)節(jié)點壓力已知）：

圖4：管網(wǎng)樣例（N：節(jié)點;B：管段;L：閉環(huán);S：氣源;D：用氣）

對于任意給定的天然氣管網(wǎng)穩(wěn)態(tài)流動過程，（1）其每個節(jié)點處流量守恒（流入節(jié)點流量等于流出節(jié)點流量）;（2）如氣流形成封閉環(huán)路，規(guī)定環(huán)路方向后（順時針或逆時針）則環(huán)路中各管段的壓力差之和為零。數(shù)學(xué)表達式如下：

（對所有節(jié)點） ? ? ? ? ? ? ? ? ? ? ? ?（31）

（對管流環(huán)路） ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （32）

Q-formulation法：圖4中，當(dāng)采用管段流量（）為未知數(shù)求解管網(wǎng)問題，將管流方程寫成如下形式：

（33）

對于管網(wǎng)中每個節(jié)點依照式（31）可列方程如下：

節(jié)點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （34.1）

節(jié)點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（34.2）

節(jié)點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（34.3）

節(jié)點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （34.4）

考慮到，式（34.4）其實是式（34.1至34.3）的線性組合（式（34.4）=式（34.1）+式（34.2）+式（34.3）），說明式（34.1至34.4）中，只有3個獨立方程，而未知量卻有5個，因此仍需要兩個獨立方程求解方程組。觀察圖4中的管路形成了2個閉環(huán)（假設(shè)氣流在閉環(huán)內(nèi)以順時針方向流動）：L1和L2;對每個閉環(huán)依照式（32）可列方程如下：

閉環(huán)：? ? ? ? ? ? ? ? ? ?（35）

閉環(huán)：? ? ? ? ? ? ? ? ? （36）

式（35、36）中采用式（33）將各管段的壓降轉(zhuǎn)換為管段阻力與流量的表達式;當(dāng)假設(shè)的管段流動方向與閉環(huán)的流動方向相反時，表示壓降需添加負(fù)號。聯(lián)立式（34.1至34.3）和（35-36）得到如下方程組（未知量）：

（37.1）

（37.2）

（37.3）

（37.4）

（37.5）

若記，由多變量牛頓法可得：

（38）

式（38）中，為Jacobian矩陣：

（39）

在用式（38）進行計算時，需要提供的初始值。各管段流量的初始值一般根據(jù)節(jié)點流量和管段阻力加以估算，例如：圖4中計算初始值可估算如下：;;;;。在進行迭代時，可將循環(huán)收斂標(biāo)準(zhǔn)定為當(dāng)向量（）中的最大元素小于。管段流量算出后，結(jié)合已知管網(wǎng)中某節(jié)點的壓力，可由式（33）解出官網(wǎng)中各節(jié)點的壓力。

P-formulation法：直接以管道中的節(jié)點壓力作為未知量進行求解。將管流方程寫成如下形式：

（40）

由式（31）列出任意三個節(jié)點的方程如下：

節(jié)點

（41.1）

節(jié)點：? （41.2）

節(jié)點：

（41.3）

由于節(jié)點壓力已知，方程組（41）中含有三個方程、三個未知量（，，），可直接用牛頓法求解。該法看似較Q-formulation法簡單，但對于氣體管網(wǎng)，方程的非線性很強，且對初值的選取非常敏感，直接求解可能會遇到不能收斂的問題。對此法的改進本文不做闡述，感興趣的讀者可參閱文獻[15]。

在使用通用性管流方程（Generalized）時，需要采用Colebrook公式求解Moody管道摩擦系數(shù)：

（42）

式（42）中，為管道絕對粗糙度，in;為管徑，ft;為雷諾數(shù)：

（43）

式（43）中，為氣體比重;為流量，MMSCFD;為氣體粘度，cp。當(dāng)采用牛頓法求解管道摩擦系數(shù)時，初始值可以令式（42）中的為無限大計算得到。結(jié)合上述求解方法，針對CP3給出的幾個穩(wěn)態(tài)工況下的氣體管網(wǎng)案例，現(xiàn)給出節(jié)點壓力及管段流量參考解。

算例3：管網(wǎng)系統(tǒng)及假設(shè)的管段流動方向如圖5所示（管道及氣體參數(shù)見文獻[1]）

【解】：

算例4：管網(wǎng)系統(tǒng)及假設(shè)的管段流動方向如圖6所示（管道及氣體參數(shù)見文獻[1]）

【解】：

算例5：管網(wǎng)系統(tǒng)及假設(shè)的管段流動方向如圖7所示（管道及氣體參數(shù)見文獻[1]）

【解】：

算例6：管網(wǎng)系統(tǒng)及假設(shè)的管段流動方向如圖8所示（管道及氣體參數(shù)見文獻[1]）

【解】：

4總結(jié)

本文給出了我校海洋油氣工程專業(yè)的本科生《海洋油氣集輸》挑戰(zhàn)性課程計算機編程項目的詳細求解方法和參照答案。該套編程項目對于中國石油大學(xué)（北京）大三的同學(xué)獨立完成存在一定困難，但在教學(xué)中取得了大多數(shù)學(xué)生的理解和認(rèn)可。本文所述的求解方法能夠幫助學(xué)生更好的完成相關(guān)教學(xué)環(huán)節(jié)，在獨立完成這些教學(xué)環(huán)節(jié)之后，學(xué)生不但能夠感受到自己在計算、分析及應(yīng)用能力上的進步，同時也對所學(xué)的內(nèi)容理解更加透徹。

基金項目：中國石油大學(xué)（北京）教育教學(xué)改革項目“海洋油氣集輸全英文核心課程建設(shè)”。

參考文獻

[1] 趙海寧.石油與天然氣工程挑戰(zhàn)性課程教學(xué)實踐——海洋油氣集輸計算機編程項目設(shè)計[J].化工高等教育，2019，36（06）：40-50.

[2] Dranchuk，P.M.and Abou-Kassem，J.H.Calculation of Z factors for natural gas using equations of state[J]. Journal of Canadian Petroleum Technology，1975，14（03）：34-36.

[3] 馮叔初，郭揆常等.油氣集輸與礦場加工（第二版）[M].中國石油大學(xué)出版社，2006.

[4] Ahmed，T.Equations of state and PVT analysis applications for improved reservoir modeling， 2nd edition[M].Gulf Professional Publishing，2016.

[5] Standing，M.B.&Katz，D.L.Density of natural gases[J].Pet. Trans.AIME，1942（146）：140-149.

[6] Wichert， E.&Aziz，K.Calculation of Z's for sour gases[J].Hydrocarbon Processing，1972，51（05）：119-122.

[7] Rachford Jr. H.H.&Rice，J.D.Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium[J].Pet. Trans. AIME，1952（195）：327-328.

[8] Peng， D.-Y.&Robinson，D.B.A new two-constant equation of state[J].Industrial Eng. Chem. Fundam，1976，15（01）：59-64.

[9] Soave，G.Equilibrium constants from a modified Redlich-Kwong equation of state[J].Chemical Engineering Science，1972（27）：1197-1203.

[10] Coats，K.H.Simulation of gas condensate reservoir performance[J].Journal of Petroleum Technology，1985，37（10）：1870-1886.

[11] Danesh，A.PVT and phase behaviour of petroleum reservoir fluids[M].Elsevier，1998.

[12] Wilson，G.A modified Redlich-Kwong Equation of state[A].Application to general physical data calculations [C].American Institute of Chemical Engineers National Meeting，1968.

[13] Whitson，C.H.&Michelsen，M.L.The negative flash[J].Fluid Phase Equilibria，1989（53）：51-71.

[14] Essenfeld，D.Network modeling and tracking of liquid preferential routes in natural gas pipeline systems [D].Master Thesis，The Pennsylvania State University，2010.

[15] Ayala， L.F.&Leong，C.Y.A robust linear-pressure analog for the analysis of natural gas transportation networks[J].Journal of Natural Gas Science and Engineering，2013（14）：174-184.

[16] Colebrook，C.F.Turbulent flow in pipes， with particular reference to the transition region between smooth and rough pipe laws[J].Journal of the Institution of Civil Engineers，1939，11（04）：133-156.

[17] Moody，L.F.Friction factors for pipe flow[J].Trans.ASME，1944（66）：671-684.