陳宏舉 王靖怡 康琦 吳海浩 宮敬

（1.中國石油大學（北京）北京 102249；2.中海油研究總院有限責任公司 北京 100028；3.中國石油天然氣管道工程有限公司 河北廊坊 065000；4.清華大學 北京 100084；5.深圳清華大學研究院 廣東深圳 518057）

油氣水多相混輸技術能夠降低分輸管道的建設費用，但多相流動的不穩(wěn)定性和氣液流動的交替性可能導致管道內壓力、持液率和氣液相流速發(fā)生周期性波動［1-2］，對精確預測多相流動現(xiàn)象提出了更高要求。目前國內多相流仿真模擬大多采用多相管流模擬軟件OLGA和LedaFlow，二者都是基于流體動力學機理的瞬態(tài)模擬工具［3-6］。中國使用多相管流模擬軟件一直受國外的制約，每年要繳納高額的使用授權費用。

為了減少對國外商業(yè)軟件的依賴，推動多相管流模擬軟件國產(chǎn)化發(fā)展，中國石油大學（北京）宮敬教授團隊自主研發(fā)了動態(tài)多相管流模擬軟件MPF［7-9］，目前處于測試階段。本文對MPF、OLGA和LedaFlow軟件的理論基礎進行介紹，以實際工程管線穩(wěn)態(tài)工況為例，比較3個軟件在預測壓降、溫降、持液率和氣液相速度等方面的差異；并研究瞬態(tài)工況下3個軟件對嚴重段塞流的處理情況，從而對比分析MPF與成熟商業(yè)軟件OLGA和LedaFlow在穩(wěn)態(tài)、瞬態(tài)工況方面的預測能力，為國產(chǎn)化軟件MPF的進一步完善提供參考。

1 OLGA、LedaFlow與MPF軟件概述

1.1 OLGA

OLGA采用基于6個守恒方程進行求解的拓展雙流體模型，包括氣相、液相或液膜、液滴3個質量守恒方程，液滴的氣相、液相或液膜2個動量方程及1個混合體系能量方程［10］。該模型采用歐拉法和歐拉形式的拉格朗日邊界跟蹤形式，通過人工假定條件或依據(jù)實驗關聯(lián)式反映本構關系。質量守恒方程中考慮了控制液滴夾帶和沉積的質量傳遞項，由此得出氣相、液相或液膜、液滴的質量守恒方程，計算氣相、液相之間的質量分布情況。此外，結合流體流動狀態(tài)的影響及氣液兩相之間因相變引起的質量交換，計算出多相流各項參數(shù)［11］。拓展的雙流體模型包括質量守恒方程、動量守恒方程、能量守恒方程［4］。

1.2 LedaFlow

LedaFlow的核心計算模型可分別求解油氣水三相的質量、動量及能量的守恒方程［5］。通過研究9個獨立物理場的特點，建立每個獨立的物理場所對應的質量方程以及每個連續(xù)區(qū)域的多混合模型所對應的動量守恒方程，LedaFlow內核模型會求解每一相的能量守恒方程，計算同一截面上三相的溫度，故可以獲得更多更準確的關于流體溫度的信息。LedaFlow的每個計算單元擁有16個方程，包括9個獨立場的質量守恒方程，連續(xù)相的3個動量方程，連續(xù)相的3個能量方程及1個體積守恒方程［5］。與OLGA相比，LedaFlow擁有更多的質量守恒方程，包含3個連續(xù)相方程，6個分散相方程，以及針對每個場的質量守恒方程、針對連續(xù)相中每一相的動量方程和針對連續(xù)相中每一相的能量方程。

關于動量方程的處理，LedaFlow認為在多相流混輸管道中移動介質的界面張力相對較小，因此可以忽略不同界面的壓差。采用交錯網(wǎng)格貢獻邊際法進行求解，把速度分配給網(wǎng)格邊界，把其余參數(shù)分配到網(wǎng)格內部。方程組內函數(shù)解析使用逆向差分法，從而提高計算穩(wěn)定性。通過修正后的SIMPLE算法進行求解［12-13］。

1.3 MPF

MPF采用雙流體模型，以變量在管道截面上的平均值為基礎，分別將兩相流視為連續(xù)介質而建立控制方程，分別描述各相的質量、動量、能量等性質［14］。MPF采用2個質量方程、2個動量方程、2個能量方程及1個混合能量方程組成的復雜7方程模型來研究兩相流動［15］。

各相質量方程：

式（1）～（5）中：k為g、l時分別表示氣相與液相；αk表示截面相含率；ρk表示密度；uk表示速度；t表示時間；x表示沿線管道的坐標；˙mk表示質量源變化率；p表示壓力；Fk表示除壓力外的其他表面力與質量力的合力；Ek表示總能；Hk表示焓；Wk表示能量源（包括外力做功和流體與外界換熱）。

關于能量方程的處理，MPF采用各相能量方程相加的方式來計算混合溫度，保證了流體總焓守恒。與OLGA、LedaFlow相比，MPF對摩阻系數(shù)的選?。ū?）更適用于中高持液率工況。

表1 MPF氣-液兩相流動摩阻封閉關系式Table 1 MPF closures for friction in gas-liquid flow

2 地形起伏濕氣管線穩(wěn)態(tài)模擬對比

2.1 管線運行基礎參數(shù)

某實際氣田新建一條管道輸送天然氣至處理廠，輸氣規(guī)模64.8×104m3／d，設計壓力4.0 MPa，出站壓力3.5 MPa，進站壓力2.5 MPa，出站溫度45℃，進站溫度3℃，總傳熱系數(shù)5.57 W／（m2·℃），全長約20.96 km。管道沿線縱斷面如圖1所示，流體組成如表2所示。

圖1 某氣田濕氣管線縱斷面圖Fig.1 Profile diagram of a wet gas pipeline in a gas field

表2 某氣田濕氣管線流體組成數(shù)據(jù)Table 2 Fluid composition data of a wet gas pipeline in a gas field

2.2 3個軟件穩(wěn)態(tài)模擬計算結果對比

2.2.1 壓降

3個軟件水力計算的內嵌模型較為相似，且都是給定終點壓力、起點溫度和起點輸量反算出起點壓力。為了驗證軟件計算結果的準確性，篩選出濕氣管線的10組運行數(shù)據(jù)（表3），模擬得到不同工況下3個軟件起點壓力及壓降計算結果（圖2）：①設計工況10中，OLGA和LedaFlow起點壓力計算值稍低于實際設計工況的起點壓力，相對誤差約為-3%；MPF起點壓力計算值稍高于實際設計工況的起點壓力，相對誤差為3%，3個軟件計算結果均在誤差允許范圍內。②工況1～9，3個軟件的起點壓力計算值均高于實際運行的起點壓力，但3個軟件之間的計算結果差距很?。籓LGA與Leda-Flow起點壓力和壓降的計算結果更接近，MPF起點壓力計算結果誤差高于另外2個軟件。結合設計工況10的計算結果來看，可能是現(xiàn)場氣質組分與給定設計工況下的組分存在差異，現(xiàn)場輸送的原料氣重組分占比更低，管道持液率偏小，因此壓降??；所以3個軟件反算出的起點壓力均大于實際運行情況。③與OLGA、LedaFlow相比，MPF的摩阻計算模型為BeggsBrill，水力算法為BeggsBrill-Modified，對摩阻系數(shù)的選取更適用于中高持液率工況，MPF壓力計算值偏大可能與軟件內置水力模型和摩阻系數(shù)的選取不同有關。但整體上，3個軟件的壓降計算偏差均在可接受范圍內，可用于工程模擬計算。

圖2 某氣田濕氣管線不同工況下3個軟件起點壓力及壓降計算結果Fig.2 Starting point pressure and pressure drop of a wet gas pipeline in a gas field under different conditions caculated by three software

表3 某氣田濕氣管線運行情況Table 3 Operating conditions of a wet gas pipeline in a gas field

2.2.2 溫降

基于現(xiàn)場運行數(shù)據(jù)，得到不同工況下3個軟件溫降計算結果（圖3）：MPF溫降變化與LedaFlow較為吻合，OLGA溫降值偏低，但3個軟件計算結果偏差不大。MPF熱力算法為Deng，傳熱模型為K＿U，采用3個能量方程分別計算各相的溫度和混合溫度。LedaFlow能夠計算同一截面上三相的溫度，從而達到求解每一相的能量守恒方程的要求。可見，MPF和LedaFlow對于傳熱模型的處理更加完善，認為其具備更細致的傳熱計算功能。

圖3 某氣田濕氣管線不同工況下3個軟件溫降計算結果Fig.3 Temperature drop of a wet gas pipeline in a gas field under different conditions caculated by three software

2.2.3 其他參數(shù)

以設計工況為例，比較穩(wěn)態(tài)運行下3個軟件的全線持液率計算結果（圖4）：3個軟件的持液率預測存在微小差異，OLGA在管線低點存在少量積液，LedaFlow持液率保持在0.001左右，MPF全線基本沒有積液產(chǎn)生。

圖4 某氣田濕氣管線設計工況下3個軟件全線持液率計算結果Fig.4 Holdup of a wet gas pipeline in a gas field under design condition caculated by three software

3個軟件對全線氣液相速度的計算結果（圖5）表明：MPF計算出的液相速度振蕩幅值與氣相速度基本重合，相較OLGA和LedaFlow的計算結果，MPF計算出的液相速度振蕩幅值范圍較大，但整體速度的變化趨勢與另2個軟件的結果相同。從水力段塞形成的角度分析，在水力段塞形成過程中，部分時刻會出現(xiàn)等效于均相的流動效果，段塞流動即液塞、分層流間歇出現(xiàn)，段塞體中會攜帶一些微小氣泡，其速度與段塞體基本一致，但OLGA與LedaFlow假設段塞內部為純液相，因此圖5a、b中氣液相流速變化波動較小，且相差較大，而MPF充分考慮段塞體內氣相流動的影響，更加貼合實際。

圖5 某氣田濕氣管線設計工況下3個軟件預測的全線氣液相速度Fig.5 Gas-liquid velocity of a wet gas pipeline in a gas field under design condition caculated by three software

3 臥底-懸鏈線立管系統(tǒng)瞬態(tài)模擬對比

3.1 臥底-懸鏈線立管運行基礎參數(shù)

以國外某深海油田一條實際運行的臥底-懸鏈線立管為例［16］對3款軟件的瞬態(tài)計算功能進行對比。立管高度約1 720 m，入口溫度90℃，出口壓力1.5 MPa，氣液比50 Sm3／Sm3；管徑203.2 mm，壁厚0.009 m；保溫材料為聚丙烯，厚度0.03 m。對該管線模型進行簡化，得到其高程和環(huán)境溫度變化情況（圖6）。

圖6 某臥底-懸鏈線立管系統(tǒng)高程和環(huán)境溫度變化Fig.6 Variation of elevation and ambient temperature of a undercover-catenary riser system

3.2 3個軟件瞬態(tài)模擬計算結果對比

3.2.1 沿線持液率

3個軟件均采用黑油模型，模擬得到該系統(tǒng)在同一時刻不同流量對應的沿線持液率變化（圖7）：3個軟件的沿線持液率變化總體趨勢較為一致；管線下傾管內液體積聚，持液率迅速升高直至為1，上傾立管內持液率基本穩(wěn)定在1左右，出口處持液率迅速下降；隨著入口流量的增大，沿線持液率走向更加穩(wěn)定。

圖7 3個軟件預測的某臥底-懸鏈線立管系統(tǒng)沿線持液率Fig.7 Holdup of a undercover-catenary riser system caculated by three software

3.2.2 立管底部壓力

3個軟件計算得到的臥底-懸鏈線立管系統(tǒng)立管底部壓力（圖8）均表現(xiàn)出良好的周期性變化，隨著入口流量的增加，立管底部壓力的波動周期逐漸縮短，符合真實物理規(guī)律［17-18］。MPF預測立管底部壓力的變化呈現(xiàn)出周期性波動趨勢，符合有回落的嚴重段塞流循環(huán)，具有穩(wěn)定的平臺期。LedaFlow預測結果為無回落的嚴重段塞流循環(huán)，與有回落循環(huán)類似，但液塞噴發(fā)過程表現(xiàn)不是很劇烈，有回落循環(huán)過程的液塞長度會高于無回落循環(huán)。OLGA預測結果為有回落循環(huán)過渡至無回落循環(huán)，立管底部壓力值也更接近LedaFlow。由于各軟件的嚴重段塞流預測模型存在差異，判別下傾管內氣體進入立管的劃分方式不同，導致三者的預測結果存在差異。MPF立管底部最大壓力值高于另外2個軟件。這是因為：MPF嚴重段塞流模型在液體回落過程中進行了簡化處理，MPF預測結果均為有回落循環(huán)，模型默認氣液噴發(fā)后的回落階段僅有液相，而立管中氣相全部噴出，在下一個液塞積累過程，認定立管內液位以下部分為純液相，而實際管內流動非常復雜。邱偉偉等［19-20］通過實驗觀察也證實了在氣液噴發(fā)過程中多相流動的復雜性，液體回落階段立管內也并非純液體。因此，MPF對于立管底部最大壓力的計算值高于OLGA和LedaFlow的計算結果，且具有穩(wěn)定的平臺期。整體上看，在嚴重段塞流的模擬中，MPF的表現(xiàn)稍遜于OLGA和LedaFlow。

圖8 3個軟件預測的某臥底-懸鏈線立管系統(tǒng)立管底部壓力Fig.8 Pressure at the bottom of the riser in a undercover-catenary riser system caculated by three software

3.2.3 段塞周期和液塞長度

3個軟件模擬得到不同流量下臥底-懸鏈線立管系統(tǒng)段塞周期和液塞長度的變化情況（表4）：①3個軟件的液塞長度計算結果均超過了立管的高度，這是嚴重段塞流形成的顯著標志；②隨著入口流量的增加，液塞長度逐漸縮短，段塞周期也不斷縮短。

表4 3個軟件預測的某臥底-懸鏈線系統(tǒng)液塞長度和段塞周期Table 4 Liquid plug length and slug fluctuation period of a undercover-catenary riser system caculated by three software

以OLGA的計算結果為基準，對比不同流量下LedaFlow和MPF計算結果的偏差（表5）?？梢姡孩費PF的段塞周期與OLGA的段塞周期偏差較??；②隨著入口流量的增大，LedaFlow的段塞周期變化幅度較為保守，不同流量下均為無回落的循環(huán)過程，因此最大液塞長度小于另2個軟件的計算結果。

表5 LedaFlow和MPF在某臥底-懸鏈線系統(tǒng)液塞長度和段塞周期的計算偏差（以OLGA為基準）Table 5 Error calculation of liquid plug length and slug fluctuation period for a undercover-catenary riser system caculated by LedaFlow and MPF software（based OLGA）

3個軟件對液塞長度的預測存在差異。MPF對于最大液塞長度的計算值高于LedaFlow和OLGA，主要原因在于：液體積累階段立管內充滿液體，最大液塞長度的計算差異主要是下傾管內液塞長度的不同導致的；MPF立管底部壓力計算值高于另2個軟件，因此基于公式計算下傾管內液塞長度相應也會偏大。此外，由于MPF對回落過程進行簡化處理，默認立管內全部為液相，但事實上油氣水三相流體在管內的流動過程十分復雜，實際液塞積累不一定全部充滿液相直至氣液噴發(fā)，且嚴重段塞流又會引發(fā)強烈的周期性波動，從而使MPF對回落過程的簡化處理導致了最大液塞長度的計算結果高于OLGA和LedaFlow的計算結果。因此，在后續(xù)開發(fā)中建議進一步研究完善MPF內嵌段塞流的計算模型。

3.2.4 立管底部氣相速度

3個軟件模擬得到的不同流量下臥底-懸鏈線立管系統(tǒng)立管底部氣相速度隨時間的變化規(guī)律（圖9）：隨著入口流量的增加，立管底部氣相速度總體呈現(xiàn)上升趨勢，氣相速度的波動周期在減小。在液塞積累過程中，OLGA和LedaFlow的氣相速度計算結果并不完全為0，說明兩款軟件在計算嚴重段塞流時，上升立管內并非純液相，這2個軟件在計算液塞積累時將立管內氣體流動情況考慮在內。MPF軟件在此方面需要進一步考慮液塞形成中氣速的精確描述。

圖9 3個軟件預測的某臥底-懸鏈線系統(tǒng)立管底部氣相速度Fig.9 Gas velocity at the bottom of the riser in a undercover-catenary riser system caculated by three software

4 結論及建議

MPF與OLGA、LedaFlow在計算管線的壓降、溫降、沿線持液率、立管底部壓力、段塞周期與液塞長度及立管底部氣相速度等方面存在差異，表現(xiàn)在：①穩(wěn)態(tài)工況下，3個軟件水力計算與現(xiàn)場實際數(shù)據(jù)相對誤差絕對值在3%以內，MPF和LedaFlow熱力計算模型能夠精確到每一相，二者計算結果更為接近；氣液相速度方面，MPF對水力段塞流形成過程的處理更符合實際物理過程。②嚴重段塞流瞬態(tài)工況下，3個軟件的持液率波動位置及頻率變化十分吻合，OLGA與MPF的波動周期相差很小，呈現(xiàn)了良好的預測結果，但由于MPF對嚴重段塞流模型進行了簡化處理，故其對立管底部壓力、最大液塞長度的計算結果均高于OLGA和LedaFlow的計算結果，立管底部氣相速度波動預測不夠準確。建議在后續(xù)研究中，進一步完善MPF中的段塞流等模型，并結合更豐富、完整的實驗數(shù)據(jù)和現(xiàn)場運行數(shù)據(jù)進行測試更新，全方位提升MPF的功能及可靠性，從而使其能準確應用于實際工程的各類工況計算。