基于天然氣放空系統(tǒng)的水力計算軟件開發(fā)

婁 晨，林 棋

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

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

圖1 高壓天然氣放空系統(tǒng)示意圖

1 天然氣放空管路水力系統(tǒng)理論

1.1 放空過程狀態(tài)描述

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

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

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

3）亞音速流：當(dāng)放空過程越過臨界流繼續(xù)泄壓時，氣體進(jìn)入亞音速流狀態(tài)。入口截面壓力p1繼續(xù)降低，出口截面壓力p2保持不變（等于pa），出口氣體流速逐漸減?。∕ach2<1），質(zhì)量流量也隨之下降，在放空干線內(nèi)壓力p0下降至pa之前，該流動狀態(tài)將持續(xù)進(jìn)行，直至天然氣放空系統(tǒng)內(nèi)各壓力值相等（pa=p0=p1=p2），放空過程結(jié)束。

1.2 天然氣放空系統(tǒng)基本方程推導(dǎo)

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

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

其他所需基礎(chǔ)方程 （氣體狀態(tài)方程及相關(guān)熱力方程）：

式中：Ti為放空系統(tǒng)不同截面的氣體溫度，K；γ為氣體的等熵指數(shù)；pi為放空系統(tǒng)不同截面的氣體壓力，Pa；Mach為氣體馬赫數(shù)；M為氣體摩爾質(zhì)量，kg/mol；R為氣體常數(shù)；a為放空系統(tǒng)所處的當(dāng)?shù)芈曀?，m/s；Cp為氣體的定壓比熱，J/kg·K；Z 為截面處氣體壓縮因子。

將方程組（2）代入方程組（1）中的質(zhì)量守恒方程可得到關(guān)聯(lián)方程組（3）：

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

將關(guān)系式（4）代入關(guān)聯(lián)方程組（3）可將其轉(zhuǎn)化為馬赫數(shù)與壓力值的關(guān)系式：

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

將方程組（3）及關(guān)系式（5）進(jìn)行積分（考慮到氣體等熵指數(shù)為常數(shù)），并將結(jié)果代入方程（6）可得：

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

式中：μ為放空管路內(nèi)壁絕對粗糙度，m。

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

氣體從輸氣干線通過圓形擴(kuò)孔進(jìn)入放空管路（p0T0?p1T1）入口截面的熱力過程，可以視為一個典型的等熵?zé)崤蛎涍^程，故由熱力學(xué)基本公式可得如下關(guān)系式[5]：

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

2 放空管路水力計算軟件開發(fā)

2.1 程序計算流程簡述

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

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

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

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

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

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

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

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

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

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

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

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

2.2 算例計算分析

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

圖3 軟件計算界面圖

軟件計算可求得放空系統(tǒng)各時刻氣源壓力值所對應(yīng)的瞬時放空量、累計放空時間及單位質(zhì)量流率、進(jìn)出口截面的壓力及溫度值、各時刻管段儲氣量部分?jǐn)?shù)據(jù)提取見表2。從中獲知：

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

2）瞬時放空量及單位質(zhì)量流率在泄壓放空過程中呈現(xiàn)出線性遞減的變化趨勢，與放空過程所經(jīng)歷的3種流態(tài)形式?jīng)]有關(guān)系。

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

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

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

表1 放空系統(tǒng)程序計算算例基本參數(shù)

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

表2 放空系統(tǒng)部分計算結(jié)果匯總表

圖4 放空系統(tǒng)內(nèi)各溫度變化曲線

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

8）同時利用圖解法對算例進(jìn)行了求解，發(fā)現(xiàn)程序計算瞬時放空量與圖解法計算結(jié)果誤差基本控制在2%以內(nèi)（圖5）。

2.3 東河伴生氣管線下沉工程放空現(xiàn)場

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

圖5 圖解法與計算軟件結(jié)果對比

東河氣站放空管線主閥為一電磁閥，通過系統(tǒng)操作全開時開度也不能達(dá)到100%，為提高放空速度，經(jīng)協(xié)調(diào)通過打開旁通閥進(jìn)行放空。放空管線長度34.4km，管徑168×6/7mm。放空前，管線壓力為6.15MPa，管段內(nèi)天然氣量標(biāo)準(zhǔn)狀況下為40 408m3。放空管線規(guī)格DN100，放空閥門開度100%。在管線放空過程中，現(xiàn)場作業(yè)中工作人員對放空時間做了詳細(xì)記錄，相關(guān)數(shù)據(jù)見表3。

利用開發(fā)的水力計算軟件，根據(jù)現(xiàn)場進(jìn)行天然氣放空系統(tǒng)建模，將程序計算結(jié)果與現(xiàn)場放空記錄數(shù)據(jù)進(jìn)行對比（圖6），進(jìn)一步驗證本計算軟件的準(zhǔn)確性。

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

表3 東河伴生氣管線放空現(xiàn)場數(shù)據(jù)記錄表

圖6 放空現(xiàn)場與計算軟件結(jié)果對比

3 結(jié) 論

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

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

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

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

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

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

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

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

[5]葉學(xué)禮.圖解法求天然氣瞬時放空量[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.