基于天然氣放空系統的水力計算軟件開發

婁 晨，林 棋

中國石油大學（北京）油氣管道輸送安全國家工程實驗室 （北京 昌平 102249）

高壓天然氣放空系統由干線內待放空管存氣體、放空管路和外界大氣組成（圖1），其中放空管路包括相應的管件及附件（包括管段、閥門、彎頭、三通及穩固定支架等）。管段內氣體流動屬于典型的非定常流動，放空過程將經歷超臨界流、臨界流及亞音速流3種狀態，管段水力計算不可簡單套用達西公式。相關的研究文獻內容較少，目前關于天然氣放空系統的放空量計算多采用一個沒有依據的公式 （即管容、放空時間），計算結果與現場結果差異極大，導致工程設計與工程實際嚴重不符，此種情況的存在將帶來嚴重的安全隱患 （比如，放空點火的熱輻射傷害）[1-2]，故有必要開發出可適用于工程實際的放空管路水力計算軟件。基于Fano方程，通過C語言編寫開發出水力計算軟件，并利用圖解法與該軟件進行瞬時放空量計算對比。同時結合東河伴生氣管線下沉工程放空現場數據，進一步驗證了計算軟件的準確性。此軟件的開發將為今后輸氣干線現場放空作業提供切合實際、安全、高效的放空方法，同時也為相關的科研研究提供便利的參考依據。

圖1 高壓天然氣放空系統示意圖

1 天然氣放空管路水力系統理論

1.1 放空過程狀態描述

高壓天然氣放空系統的實現是在短距離高壓差的條件下，其放空過程非常劇烈，放空管段沿線的氣體基本參數（溫度T、密度ρ、速度V、壓力p、馬赫數Mach等）變化快、差異大，該過程為非定常流動，需要經歷3個流動狀態：超臨界流（即：雍塞流）、臨界流、亞音速流。

1）超臨界流：放空初始時刻，放空干線內壓力p0大，從而使得入口截面處壓力p1很高，此時出口截面壓力p2遠大于外界大氣壓，放空管路出口處的氣體流速達到當地聲速 （即：臨界流速，此時馬赫數Mach2為1），對于給定的放空設施，此時將達到最大瞬時放空量（即：質量流量）。出口截面的氣體依靠剩余壓力差（p2-pa）進行膨脹，直至壓力降至外界大氣壓。隨著超臨界流的時間推移，放空管路進出口壓力均逐漸降低，但在出口壓力p2降至大氣壓之前，此放空過程仍處于超臨界流狀態，在新的壓力條件下放空過程達到新的最大瞬時放空量。

2）臨界流：隨著超臨界流的進行，p0、p1及 p2均逐步減小，當出口截面壓力p2下降到外界大氣壓力p0時，超臨界流結束，此刻便為臨界流狀態，此時放空管路出口處的氣體流速為臨界流速（Mach2為1）。

3）亞音速流：當放空過程越過臨界流繼續泄壓時，氣體進入亞音速流狀態。入口截面壓力p1繼續降低，出口截面壓力p2保持不變（等于pa），出口氣體流速逐漸減小（Mach2<1），質量流量也隨之下降，在放空干線內壓力p0下降至pa之前，該流動狀態將持續進行，直至天然氣放空系統內各壓力值相等（pa=p0=p1=p2），放空過程結束。

1.2 天然氣放空系統基本方程推導

設放空過程中在微元時間?t內，氣源滯止壓力p0保持不變，在等截面水平放空管路中取?x微元段作為控制體[3]。在此條件下的控制體流動基本方程（質量守恒方程、能量守恒方程及動量守恒方程）：

式中：m 為單位面積質量流率，kg/m2·s；ρ為氣體密度，kg/m3；i為單位質量的氣體熱焓，J/kg；vi為放空系統不同截面的氣體速度，m/s；d為放空管內徑，m；f為水力摩阻系數；V為氣體流速，m/s；式中下標1表示放空管路入口截面，下標2表示放空管路出口截面。

其他所需基礎方程 （氣體狀態方程及相關熱力方程）：

式中：Ti為放空系統不同截面的氣體溫度，K；γ為氣體的等熵指數；pi為放空系統不同截面的氣體壓力，Pa；Mach為氣體馬赫數；M為氣體摩爾質量，kg/mol；R為氣體常數；a為放空系統所處的當地聲速，m/s；Cp為氣體的定壓比熱，J/kg·K；Z 為截面處氣體壓縮因子。

將方程組（2）代入方程組（1）中的質量守恒方程可得到關聯方程組（3）：

同理將方程組（2）代入方程組（1）中的能量方程式可得關系式：

將關系式（4）代入關聯方程組（3）可將其轉化為馬赫數與壓力值的關系式：

同理將方程組（2）代入方程組（1）中的動量守恒方程可得關系式：

將方程組（3）及關系式（5）進行積分（考慮到氣體等熵指數為常數），并將結果代入方程（6）可得：

由科爾布魯克摩阻系數計算式（8）所示，管路中摩阻系數取決于雷諾數Re及管內壁粗糙度，在放空管路中，高速泄壓的高流速使得管路雷諾數Re非常大，進而其對管路沿線的摩阻系數影響很小，在此為簡化方程，忽略雷諾數Re的影響（即：管路摩阻系數僅由選定管路的內壁粗糙度決定）[4]：

式中：μ為放空管路內壁絕對粗糙度，m。

由此可對上述方程式（7）沿放空管路長度L進行積分，可得：

氣體從輸氣干線通過圓形擴孔進入放空管路（p0T0?p1T1）入口截面的熱力過程，可以視為一個典型的等熵熱膨脹過程，故由熱力學基本公式可得如下關系式[5]：

上述方程 （9）為有摩擦絕熱一維流動方程式（即：Fano方程），將上述10個方程式（組）作為天然氣放空系統水力求解計算關聯方程組。

2 放空管路水力計算軟件開發

2.1 程序計算流程簡述

1）放空管路基礎數據。方程式（9）中左邊項 fL/d（簡記為U）為阻力項，該項的求解可以分為3部分，如式（11）所示：

其中入口及出口的阻力項可以通過查閱資料獲取，且兩者占總阻力項的比例也較小。放空管段中L為其等效總長度（包括：彎頭、三通、閥門的等效長度以及放空管段本身管長），再結合選取管段的管內壁絕對粗糙度，即可求得Fano方程阻力項。

2）超臨界流初始時刻質量流量（瞬時放空量）。將方程式（7）中的氣體等熵指數γ視為定值（Cp/Cv=1.3），放空初始時刻放空管路出口處的氣體流速為臨界流速（Mach2為 1），結合方程式（11）可以將方程式（7）轉化為僅有Mach1的迭代計算式：

通過程序的迭代計算可以求得入口截面氣體馬赫數 Mach1，將其代入方程式（5）以及方程組（10）中的壓力關系式，即可求得放空管路進出口壓力p1、p2。將p1及Mach1代入方程組（10）中的溫度關系式，可求得入口截面氣體溫度T1，再將T1、Mach1及Mach2代入方程式（4）中，可求得出口截面氣體溫度T2。經過上述求解后，可由式（13）計算求得超臨界流的初始質量流量：

式中：W為瞬時放空量，kg/s；Z1為截面處天然氣壓縮因子。

3）臨界流時刻質量流量（瞬時放空量）。由上述放空過程描述可知，當放空氣體處于臨界流時，出口截面壓力p2始終等于大氣壓力pa，出口氣體馬赫數Mach2為1，由迭代計算式（12）可知，入口截面氣體馬赫數Mach1求解值不變，代入方程式（5）可得入口截面壓力p1，再代入方程組（10）壓力關系式，可求得臨界流狀態時干線管存壓力值p0，由方程組（10）溫度關系式可知入口截面氣體溫度求解值不變。將Mach1、Mach2及T1代入方程式（4）可求得出口截面氣體溫度T2，再同理利用方程式（13）可求得臨界流時刻質量流量。

4）亞音速流質量流量（瞬時放空量）。當放空過程處于亞音速流時，出口截面壓力p2保持不變（等于pa）。針對臨界流狀態計算所得的干線管存壓力值，采用壓力等分遞減求解計算方法。取一個合適的壓力遞減值Δp0，由此可得下一時刻的p0，代入方程組（10）壓力關系式，求得p1，由于此時迭代計算式（12）存在Mach1及Mach2兩個未知量，故先假設Mach2=1，通過迭代計算求得Mach1，之后將Mach1作為已知值，將方程式（9）轉化為求解Mach2的迭代計算式進行迭代求解，計算出Mach2。重復上述過程，直至上述兩個迭代式所求得的Mach1及Mach2均達到穩定收斂狀態，將最終時刻迭代值賦予Mach1及Mach2。此時同理由方程組 （10）溫度關系式求得入口截面氣體溫度、由方程式（4）可求得出口截面氣體溫度T2、有方程式（13）可求得亞音速流任一壓力遞減值所對應的質量流量。

5）放空過程累計放空時間。對于天然氣放空系統，最為關注的因素就是其完成整個放空過程的累計放空時間。上述介紹了各個流態下的質量流量計算思路，但僅由上述方法及關聯方程組是無法求得累計放空時間，在此引入一個經典基本數值積分法——“梯形法”。將初始干線管存壓力以某一微小壓力遞減值進行等分：

由2）～4）計算推導過程，可對上述每一個壓力等分遞減值計算出其所對應的瞬時放空量Wi，再結合方程組（2）及放空干線標準管存體積計算式，通過式（15）求得各個計算時刻的管存儲氣量Qi，由于微小壓力遞減值非常小，在每個計算區間可以看成勻速放空過程，由此可以求得各個壓力區間段的等效放空時間ti，將各區間的放空時間求和即可求得放空過程的累計放空時間T。

式中：T為累計放空時間，h；Q為輸氣干線管段儲氣量，m3。

6）放空管路管徑選取。上述推導了放空系統的水力計算詳細計算步驟，由規定的總累計放空時間可以反推計算求得所需要選取的放空管路管徑，由此可更好地運用在現場工程實際。反算流程如圖2所示。

圖2 部分C語言程序計算流程圖

2.2 算例計算分析

為檢驗程序的準確性，利用所開發的軟件 （圖3），建立放空系統算例模型（表1），并同時采用圖解法進行求解，API 521提供了以可壓縮流體有摩擦絕熱一維流動的Fano方程為基礎的圖解近似計算法（Lapple圖）[6]。將兩者結果進行對比驗證。

圖3 軟件計算界面圖

軟件計算可求得放空系統各時刻氣源壓力值所對應的瞬時放空量、累計放空時間及單位質量流率、進出口截面的壓力及溫度值、各時刻管段儲氣量部分數據提取見表2。從中獲知：

1）累計放空時間隨管存壓力值的降低逐漸增大，放空時間的增長速率也逐漸增大，且在亞音速流區域呈現出急速增長。

2）瞬時放空量及單位質量流率在泄壓放空過程中呈現出線性遞減的變化趨勢，與放空過程所經歷的3種流態形式沒有關系。

3）放空管路入口截面壓力p1及輸氣干線管存儲氣量隨干線壓力的減少呈現出與累計放空時間相反的變化規律，在雍塞流區域急速下降，而在亞音速流區域下降幅度平緩。

4）放空管路出口截面壓力p2在雍塞流區域急速下降至大氣壓力pa，然后在亞音速流區域維持該壓力狀態不變。

5）放空管路入口截面溫度T1略小于輸氣干線氣體溫度T0（由于開孔處微弱的節流），在雍塞流區域維持溫度不變，在亞音速流區域緩慢升溫至干線氣體溫度值，從而達到平衡（圖4）。

表1 放空系統程序計算算例基本參數

6）放空管路出口截面溫度T2受初始放空時刻的強烈節流效應，溫度急速下降至某一最低溫度，之后隨放空過程的進行，逐漸升溫，最終達到T2=T1=T0；此升溫過程在雍塞流區域上升幅度較慢，而在亞音速流區域呈現出急速回溫的變化趨勢（圖4）。

表2 放空系統部分計算結果匯總表

圖4 放空系統內各溫度變化曲線

7）在雍塞流狀態下壓降速度快，壓降變化值占到總放空過程的90%左右，且持續時間較短（小于50%），瞬時放空量大，會使得放空管路出口處急劇降溫（初始瞬間可達-40℃左右）；亞音速流狀態下壓降速度緩慢，壓降變化值僅占到總放空過程的10%左右，但持續時間卻大于雍塞流狀態，放空完成時刻系統各壓力值與大氣壓pa一致。

8）同時利用圖解法對算例進行了求解，發現程序計算瞬時放空量與圖解法計算結果誤差基本控制在2%以內（圖5）。

2.3 東河伴生氣管線下沉工程放空現場

2013年9月底在南疆鐵路二線建設中，輪臺縣至庫車縣A標志段B涵洞處，因通過涵洞下的公路路面基線與鐵路橋距離太小，當地車輛通行不便，需要降低路面，使此處東河伴生氣管線埋深不符合設計要求，因此需要進行伴生氣管線沉降施工。需將東河伴生氣管線出站至3號閥室進行放空作業。本工程的施工程序如下：原管線開挖—管線放空—氮氣置換—新管線預制—焊縫無損檢測—割管—封堵—新管線安裝—焊縫無損檢測—補口—電火花檢漏—補傷—管線吹掃—氮氣置換—管線復產—管線回填。

圖5 圖解法與計算軟件結果對比

東河氣站放空管線主閥為一電磁閥，通過系統操作全開時開度也不能達到100%，為提高放空速度，經協調通過打開旁通閥進行放空。放空管線長度34.4km，管徑168×6/7mm。放空前，管線壓力為6.15MPa，管段內天然氣量標準狀況下為40 408m3。放空管線規格DN100，放空閥門開度100%。在管線放空過程中，現場作業中工作人員對放空時間做了詳細記錄，相關數據見表3。

利用開發的水力計算軟件，根據現場進行天然氣放空系統建模，將程序計算結果與現場放空記錄數據進行對比（圖6），進一步驗證本計算軟件的準確性。

由對比結果可知：兩者數據對比基本一致，放空過程壓降變化曲線重合度較高，軟件計算所得壓降曲線較為圓滑規整，其超臨界流持續時間稍短于實際放空現場，軟件計算瞬時放空量稍大；現場記錄數據顯示，當旁通處于關閉時，放空速度較慢，打開旁通后，放空速度加快，壓降變化曲線與軟件計算基本一致，軟件計算累計放空時間略小于實際放空現場。上述存在的差異主要為：①軟件計算模型暫未考慮放空系統的等效長度（閥門、彎頭及三通等）；②現場放空時，旁通打開的時間是在放空進行了一個多小時之后。由于這兩點的差異性，導致了軟件計算天然氣放空速度略快。但從總體上看，軟件計算精度較高，可以應用到實際現場。

表3 東河伴生氣管線放空現場數據記錄表

圖6 放空現場與計算軟件結果對比

3 結 論

1）基于C語言編寫的天然氣放空系統水力計算軟件，操作方便、計算速度快，利于工程實際應用。

2）由軟件計算算例分析可得：超臨界流狀態下壓降速度快，壓降變化值大，持續時間較短，瞬時放空量大，使得放空管路出口處急劇降溫；亞音速流狀態下壓降速度緩慢,壓降變化值小,但持續時間較長。

3）計算軟件與圖解法計算結果進行對比，精度較高（誤差在2%以內），并彌補了圖解法僅可求得瞬時放空量而無法得到其他參數的缺陷。

4）通過東河伴生氣管線現場天然氣放空所記錄的數據，利用軟件建模計算，對比結果顯示軟件計算結果與實際現場放空吻合度高，可將軟件應用于工程實際，為現場放空作業提供準確、安全、高效的放空方案。

[1]周雪漪.計算水力學[M].北京：清華大學出版社,1995.

[2]API R P521 Guide for pressure relieving and de-pressuring system[S].

[3]葉學禮.天然氣放空管路水力計算[J].天然氣工業，1999,19（3）:90-94.

[4]李玉星,姚光鎮.輸氣管道設計與管理[M].東營：中國石油大學出版社,2009.

[5]葉學禮.圖解法求天然氣瞬時放空量[J].天然氣與石油,1999,17(2):1-3.

[6]K.K.BOTROS,W.M.JUNGOWSKI,M.H.WEISS.Models and MethodsofSimulating GasPipeline Blowdown[J].The CANADIAN JOURNAL OF CHEMICAL ENERGINEERING,1989，67（8）:529-539.