石油與天然氣工程挑戰性課程教學實踐

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

關鍵詞 油氣集輸 編程項目 算例 閃蒸計算 管網計算

中圖分類號：G642文獻標識碼：A

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

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

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

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

注：表中氣體摩爾分數取自馮叔初等;臨界溫度、臨界壓力取自Ahmed;R = F + 459.67。

【解】：

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

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

（1）

此表達式與溫標的選取有關（對于開氏溫標：）。式（1）中代表CO2和H2S在天然氣中總的摩爾分數;B代表H2S在天然氣中的摩爾分數。本例中，

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

表2：DAK狀態方程求解壓縮因子計算迭代結果

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

;

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

。

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

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

算例2：CP2的主要教學任務和教學目標參見文獻[1]。

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

第1步：求解立方形狀態方程（編寫Zfactor子函數）。

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

（2）

（3）

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

（4）

（5）

（6）

（7）

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

（8.1）

（8.2）

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

（9）

（10）

其中，

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

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

（12）

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

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

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

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

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

（13）

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

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

第2步：求解Rachford-Rice目標函數（編寫VLE子函數）。

由氣液相物料守恒可以推導得出Rachford-Rice目標函數具有如下形式：

（14）

式（14）中，為體系中第種組分在體系中總的摩爾分數;為體系達到平衡狀態時的氣化分率（氣相中物質的量與體系中總物質的量的比值）;為第種組分的平衡常數，可由Wilson經驗公式近似計算得到：

（15）

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

（16）

當時停止迭代（為預先擬定的收斂標準）。在迭代計算過程中，需要給賦初值使得迭代計算開始（通常?。瑲饣致实挠嬎阒惦m然可能超出實際取值區間[0，1]，只要其最終的計算結果在區間[0，1]之間即可，為了在氣化分率超出實際取值區間時迭代計算仍然能夠繼續進行，在每一步迭代后都進行如下判斷：

（17）

式（17）中，和分別為組分平衡常數的最大值和最小值。當迭代收斂后，氣液相的平衡組分可由式（18-19）求出：

（18）

（19）

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

第3步：獲得準確的平衡常數（編寫MAIN程序）。

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

當體系中氣液兩相處于平衡狀態時，相間各組分逸度相等，可表示為：

（20）

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

（21）

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

（22）

在迭代開始時，需要給平衡常數賦初值，該初值往往由Wilson公式（式15）估算得到。在迭代過程中，當滿足如下收斂標準時（式23）停止迭代，認為平衡常數已不再變化。

（23）

由式（22）可以看出，平衡常數由逸度系數計算更新得到，而氣液相的逸度系數（）經過嚴格的熱力學推導可由狀態方程得到（式24）：

（24）

式中，

（24.1）

（24.2）

（24.3）

（24.4）

（24.5）

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

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

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

（1）儲罐中原油收率：

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

（2）產出液氣油比（GOR）：

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

（25）

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

（26）

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

（27）

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

（28）

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

（3）儲罐原油API重度：

從實際氣體狀態方程，可得：

（29）

（30）

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

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

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

圖4：管網樣例（N：節點;B：管段;L：閉環;S：氣源;D：用氣）

對于任意給定的天然氣管網穩態流動過程，（1）其每個節點處流量守恒（流入節點流量等于流出節點流量）;（2）如氣流形成封閉環路，規定環路方向后（順時針或逆時針）則環路中各管段的壓力差之和為零。數學表達式如下：

（對所有節點） ? ? ? ? ? ? ? ? ? ? ? ?（31）

（對管流環路） ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （32）

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

（33）

對于管網中每個節點依照式（31）可列方程如下：

節點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （34.1）

節點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（34.2）

節點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（34.3）

節點：? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （34.4）

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

閉環：? ? ? ? ? ? ? ? ? ?（35）

閉環：? ? ? ? ? ? ? ? ? （36）

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

（37.1）

（37.2）

（37.3）

（37.4）

（37.5）

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

（38）

式（38）中，為Jacobian矩陣：

（39）

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

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

（40）

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

節點

（41.1）

節點：? （41.2）

節點：

（41.3）

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

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

（42）

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

（43）

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

算例3：管網系統及假設的管段流動方向如圖5所示（管道及氣體參數見文獻[1]）

【解】：

算例4：管網系統及假設的管段流動方向如圖6所示（管道及氣體參數見文獻[1]）

【解】：

算例5：管網系統及假設的管段流動方向如圖7所示（管道及氣體參數見文獻[1]）

【解】：

算例6：管網系統及假設的管段流動方向如圖8所示（管道及氣體參數見文獻[1]）

【解】：

4總結

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

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

參考文獻

[1] 趙海寧.石油與天然氣工程挑戰性課程教學實踐——海洋油氣集輸計算機編程項目設計[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].中國石油大學出版社，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.