侯 慶 民
(哈爾濱商業(yè)大學(xué) 能源與建筑工程學(xué)院,哈爾濱 150028)
當(dāng)今世界,經(jīng)濟和社會飛速發(fā)展,能源需求高速增長,溫室氣體排放量數(shù)目驚人,環(huán)境污染日益嚴(yán)重,對經(jīng)濟的可持續(xù)發(fā)展和人類的長久生存構(gòu)成了巨大威脅,嚴(yán)峻的形勢使得世界各國高度注重節(jié)能和環(huán)保.天然氣具有熱值高、環(huán)境污染小、儲量豐富和經(jīng)濟效益高等方面的優(yōu)點,對改善能源結(jié)構(gòu)、緩解能源供需矛盾、提高環(huán)境質(zhì)量起到重要作用.天然氣已經(jīng)成為繼煤炭、石油之后的第三大能源,被譽為21世紀(jì)的環(huán)保能源[1-2].目前,公路、鐵路、水運和管道是目前天然氣最主要的運輸方式.其中,管道運輸與其他幾種方式相比,具有以下幾點優(yōu)勢:1)運量大,可以連續(xù)運輸,且自動化管理水平高;2)建設(shè)周期短、費用少、占地少;3)安全可靠,能源消耗少;4)成本低,效益好;5)易于適應(yīng)復(fù)雜地形、地貌和自然條件.因此,天然氣的管道運輸已經(jīng)成為天然氣最主要的運輸方式,在國民經(jīng)濟中起到了越來越重要的作用[3-6].雖然管道運輸天然氣的安全性優(yōu)于其他方式.但由于天然氣屬于危險物,一旦天然氣管道發(fā)生泄漏或斷裂,就會對其周圍的環(huán)境和人員造成極大的危害.因此研究正常工況和泄漏工況下的管內(nèi)天然氣流動,對于預(yù)防事故發(fā)生以及事故發(fā)生后的應(yīng)急處理都有著重要的實際意義.
天然氣流動的基本控制方程主要包括:連續(xù)性方程、動量方程和能量方程.另外,對于氣體,還應(yīng)該加上氣體狀態(tài)方程.天然氣流動的基本控制方程適用于天然氣管道的正常工況和泄漏工況.
在建立管道內(nèi)天然氣流動基本控制方程時,進行合理的簡化,提出如下三點假設(shè)[7-8]:
1)管內(nèi)天然氣流動為一維流動;
2)管道為剛體;
3)對一段天然氣管線,其傾斜度為定值.
在以上的假設(shè)條件下,可以得到天然氣流動的基本控制方程:
1)連續(xù)性方程
根據(jù)質(zhì)量守恒定律,可以建立天然氣流動的連續(xù)性方程:
(1)
其中:t為時間(s);x為距管道起點的距離(m);ρ為天然氣密度(kg/m3);u為天然氣流速(m/s).
2)運動方程
根據(jù)動量守恒定律,可以建立天然氣流動的動量方程:
(2)
其中:θ為管道傾斜角(rad);g為重力加速度(m/s2);D為管道內(nèi)徑(m);P為天然氣壓力(Pa);λ為水力摩擦系數(shù).
3)能量方程
根據(jù)能量守恒定律,可以建立天然氣流動的能量方程:
(3)
其中:q為單位流量的天然氣所放出的熱量(J/kg);e為天然氣的內(nèi)能(J/kg);z為天然氣管道的位置高度(m);h為天然氣的焓(J/kg).
4)氣體狀態(tài)方程
(4)
其中:T為天然氣溫度(K);Rm為氣體常數(shù)(J/(kg.K));Z為氣體壓縮因子.
管內(nèi)天然氣的流動可分為瞬態(tài)流動和穩(wěn)態(tài)流動兩大類,其中瞬態(tài)流動是指管道沿程上任一點的流動參數(shù)不僅與該點的位置有關(guān),還與時間有關(guān).因此以上天然氣流動基本控制方程(1)~(4)也可稱為天然氣管道瞬態(tài)流動方程.而穩(wěn)態(tài)流動是指管道沿程上任一點的流動參數(shù)只與該點的位置有關(guān),而與時間無關(guān).當(dāng)為穩(wěn)態(tài)流動時,方程(1)~(4)中對時間的偏導(dǎo)項為零.
方程(1)~(4)較精確地描述了天然氣管道的瞬態(tài)流動,對于工程實際問題,還可以進一步簡化.對于天然氣長直管道,管道外壁與接觸土壤之間的溫差很小,因而管內(nèi)天然氣與管壁的熱傳導(dǎo)很小,可近似認(rèn)為管內(nèi)天然氣的流動是等溫流動,在這種假定下計算得出的結(jié)果是可以接受的.因此當(dāng)不考慮管道內(nèi)天然氣溫度變化時,可以忽略能量方程(3);天然氣管道中天然氣流速和聲速相比很小,所以運動方程中的對流項可以忽略;當(dāng)標(biāo)高的差值不大時,重力項也可以忽略;同時,控制方程(1)、(2)是以單位體積(m3)流體為研究對象得出的,以單位質(zhì)量流體的控制方程在實際應(yīng)用中的方便性,引入質(zhì)量流量Q=ρuA.因此,可得到以下的簡化方程[9]:
(5)
(6)
其中:Q為管道天然氣的平均質(zhì)量流量(kg/s);A為管道的截面積(m2);c為聲速(m/s).
方程(5)、(6)是一組雙曲型偏微分方程,求解較復(fù)雜,往往采用數(shù)值解法來求解,而最成熟且最常用的數(shù)值解法就是特征線法.應(yīng)用特征線法求解偏微分方程(5)、(6)的基本思想為:先將方程(5)、(6)變換為沿特征線的常微分方程組,將常微分方程組化為差分方程組,再用數(shù)學(xué)方法化為非線性方程組,結(jié)合初始條件和邊界條件進行數(shù)值求解.
將方程(5)、(6)用E1、E2表示:
(7)
兩個方程用一任意的未知乘數(shù)進行線性組合如下:
(8)
設(shè)法選擇兩個δ值,使方程(8)變?yōu)橛米兞縋和Q表示的特殊的常微分方程.設(shè)P=P(x,t)和Q=Q(x,t)是方程(2-8)的解,則他們的全導(dǎo)數(shù)為:
(9)
(10)
將式(9)與式(8)的第一項、式(10)與式(8)的第二項括號內(nèi)進行比較,如果滿足以下條件:
(11)
(12)
則式(8)可以變?yōu)槌N⒎址匠蹋矗?/p>
(13)
式(11)、(12)是相等的,即:
(14)
由上式解得δ±A/c,說明δ是兩個不同的實數(shù),把δ代入式(11)和(12)得:
(15)
式(15)是式(13)必須滿足的兩個條件.
把δ的值分別代入式(13),并與式(15)對應(yīng)組合,就得到兩個常微分方程組:
(16)
(17)
以空間步長Δx和時間步長Δt將管道的時空平面離散化,劃分特征線網(wǎng)格,如圖1所示.
圖1 正常工況下的特征線網(wǎng)格
圖1中,MO、NO是兩條特征線,分別對C+、C-沿兩條特征線積分,引入一階近似,可以得到差分方程:
(18)
(19)
方程(18)、(19)是O點在jΔt時刻滿足的兩個方程,聯(lián)立方程(18)和方程(19),可以解得O點在jΔt時刻的各分段點上的壓力和流量.管道上除兩個端點外,每個分段點上都可以建立類似的差分方程.所以可以從t=0時刻一直到計算到所要求的時間,其中t=0時刻的初始條件已知.
如圖1所示,管道的兩端只存在一個特征線方程,無法求解端點處的壓力和流量.因此需要在管道的兩端輔助以關(guān)于壓力和流量的邊界條件,就可以求得邊界上的各個時刻的壓力和流量[10].
2.2.1 天然氣管道沿線壓力分布
天然氣流動的運動方程(2)可以變?yōu)椋?/p>
(20)
(21)
對于水平管道,式(21)可簡化為:
(22)
將流速公式代入式(22),并對其積分,可得:
(23)
其中:P(0)為管道首端壓力;P(L)為管道末端壓力.
整理式(23)可得天然氣管道的質(zhì)量流量為:
(24)
(25)
一段長度為L的水平天然氣管道AC,管道上任意一點B至管道首端A的距離為x,B處的壓力為P(x),管道首端壓力為P(0),末端壓力為P(L).
對AB段有:
(26)
對BC段有:
(27)
由質(zhì)量守恒定律可知,AB段與BC段的流量相等,所以聯(lián)立式(26)、(27)可以求得管道沿線壓力分布為:
(28)
2.2.2 天然氣管道沿線溫度分布
上文為了求解天然氣流動基本控制方程,近似認(rèn)為管道內(nèi)的天然氣流動是等溫流動.現(xiàn)實中,由于天然氣管道壁與接觸土壤間的熱傳導(dǎo),管道內(nèi)天然氣的溫度沿程是有變化的,因此與等溫流動會有一些偏差,對長直天然氣管線必須進行必要的修正[12].
對于穩(wěn)態(tài)流動,能量方程(3)可簡化為:
(29)
即:
(30)
由熱力學(xué)關(guān)系式可知:
(31)
將式(31)代入式(30)可得:
(32)
忽略式(32)中有關(guān)流速變化和高差的影響,方程可寫成:
(33)
由熱力學(xué)循環(huán)關(guān)系式有:
(34)
由氣體的定壓比熱與焦耳-湯姆遜系數(shù)關(guān)系式有:
(35)
(36)
其中:CP為氣體定壓比熱容(J/(kg·K));Di為焦耳-湯姆遜系數(shù)(℃/MPa).
所以式(33)可寫為:
-dq=CPdT-CPDidP.
(37)
由傳熱學(xué)知識可得:
(38)
其中:Kt為管道的總傳熱系數(shù)(W/(m2·K));Ts為管道周圍環(huán)境溫度(K).
式(37)可變?yōu)椋?/p>
(39)
(40)
式(40)為一階線性非齊次微分方程,其通解為:
(41)
記管道首端,即x=0處的溫度為T=T(0),利用式(41)可求得積分常數(shù)C為:
C=T(0)-Ts
(42)
將式(42)代入(41)就可求得管道沿線的溫度分布:
(43)
T(x)=Ts+(T(0)-Ts)e-βx+
(44)
在天然氣管道中,式(44)中的最后一項一般為3~5 ℃/MPa,如果忽略其影響,則管道沿線的溫度分布公式為:
T(x)=Ts+(T(0)-Ts)e-βx
(45)
與正常工況下的計算方法一樣,利用特征線法求解簡化方程(5)、(6),同時將管道上某點發(fā)生泄漏包含進去.如圖2所示,假設(shè)j時刻節(jié)點i處突發(fā)泄漏.
圖2 泄漏工況下的特征線網(wǎng)格
同樣引入一階近似,可得j時刻的節(jié)點i沿特征線C+和C-的差分方程:
(46)
(47)
方程(46)、(47)是在jΔt時刻滿足的兩個方程,因為方程(46)與(47)中有四個待求量(泄漏點前后的壓力和流量),所以聯(lián)立這兩個方程并不能直接解得jΔt時刻泄漏點前后的壓力和流量.所以在泄漏發(fā)生后,要在泄漏點處輔助關(guān)于泄漏點前后的壓力和流量的兩個邊界條件,即兩個邊界條件方程.一般情況下,假設(shè)泄漏在水平方向引起的動量變化可以忽略,即泄漏點前后的壓力不變;并且如果能求得泄漏量,根據(jù)質(zhì)量守恒原理,可以給出泄漏點前后的流量關(guān)系,此時邊界條件為:
(48)
(49)
將泄漏點處的邊界條件(48)和(49)與差分方程(46)和(47)聯(lián)立,四個方程中有四個待求量,即可求得任意時刻泄漏點前后的壓力和流量[13].
在管道端點,與正常工況的計算方法一樣,只要在邊界處輔助關(guān)于壓力、流量的邊界條件.即可以求得邊界上的各個時刻的壓力和流量.泄漏點和管道端點以外的各個節(jié)點的壓力和流量的計算方法與正常工況時的計算方法一樣.
瞬態(tài)下的天然氣管道泄漏量是隨時間變化的,最大泄漏量一般是在泄漏剛發(fā)生時.而在工程實際中,一般只要知道最大泄漏量即可,此時的最大泄漏量即為以穩(wěn)態(tài)下的泄漏量.因此本文主要研究穩(wěn)態(tài)下天然氣管道泄漏量的計算方法.
天然氣管道大多是埋地的,但由于埋地管道的外部環(huán)境比較復(fù)雜,所以很難建立其泄漏模型并計算泄漏量.因此,本文所研究的管道泄漏模型是指天然氣直接向大氣中泄漏的情況,具體為架空天然氣管道泄漏或者第三方破壞導(dǎo)致的埋地天然氣管道暴露在大氣環(huán)境下的泄漏.
圖3是典型的天然氣管道泄漏的示意圖,假設(shè)距離天然氣管道起點X處發(fā)生了泄漏.點1為天然氣管道軸線的起點位置,點4為天然氣管道軸線的終點位置,點3為泄漏點,點2為泄漏點所對應(yīng)的管道軸線上的位置.圖中P、T、u和ρ分別表示各點處的壓力、溫度、流速和密度,Q-和Q+表示泄漏點前后的質(zhì)量流量,K表示泄漏點處的泄漏量.當(dāng)泄漏發(fā)生時,點1處、點4處、大氣環(huán)境的參數(shù)值以及泄漏點前的質(zhì)量流量是已知的,而點2與點3處的參數(shù)值是未知的,需要從已知的參數(shù)入手,結(jié)合工程熱力學(xué)和流體力學(xué)的知識,推算出點2及點3處的參數(shù)值,從而計算得到天然氣泄漏量.即建立天然氣管道泄漏模型,得到泄漏量的計算方法[14].
圖3 天然氣管道泄漏示意圖
一般情況下,將天然氣管道泄漏視為孔口泄漏,因此泄漏過程也可看作可壓縮氣體的孔口出流,并且將不同形狀的泄漏孔都轉(zhuǎn)化為圓孔口來計算.為了便于計算,這里做如下的假定:(a)管道內(nèi)的天然氣流動是一維的;(b)管道內(nèi)為絕熱流動;(c)泄漏點處為等熵流動;(d)遵守理想氣體的運動規(guī)律.由于假設(shè)管內(nèi)天然氣為絕熱流動,因此,利用能量守恒和動量守恒方程可以得到:
(50)
其中:G為管道內(nèi)單位面積的質(zhì)量流量(kg/(m2.s));k為氣體絕熱系數(shù),對于天然氣一般取k=1.334;M為天然氣的摩爾質(zhì)量(kg/mol);Xe為點1到點2的管道當(dāng)量長度(m);Pi為點i處的壓力,i=1,2,3(Pa);Ti為點處的溫度,i=1,2,3(K);R為通用氣體常數(shù)(J/(mol.K)).
其中:水力摩擦系數(shù)λ是管道粗糙度δ和雷諾數(shù)Re的函數(shù).當(dāng)Re≤2 000時,λ=16/Re.
在過渡區(qū)和湍流情況下,λ和Re關(guān)系由柯爾布魯克公式給出:
(51)
當(dāng)Re較大,達到完全湍流時,λ和Re無關(guān).
把連續(xù)性方程(1)及氣體狀態(tài)方程(4)代入式(50),在泄漏孔處等熵流動的情況下,可求得孔口處天然氣泄漏量K的表達式為:
(52)
其中:K為天然氣泄漏量(kg/s);Aor為泄漏孔的面積(m2);CD為孔口流量修正系數(shù),亦指泄漏系數(shù).
對于孔口流量修正系數(shù),CD范圍為0.6~1.0.如果按照孔口形狀來分:圓形孔,CD=1.0;三角形孔,CD=0.95;長方形孔,CD=0.90;漸縮孔,CD=0.90~1.0;漸擴孔,CD=0.60~0.90.
泄漏點處的天然氣流動屬于音速流動還是亞音速流動對泄漏量的計算有很大的影響.考慮外界大氣壓已知,本文設(shè)臨界壓力比(Critical Pressure Ratio,CPR)為[15-16]:
(53)
其中:P2c代表點2處的臨界壓力,Pa.
如果點2處的壓力P2逐漸增加,天然氣泄漏速度也會一直增加,直到等于當(dāng)?shù)芈曀?這時如果P2繼續(xù)增加,天然氣泄漏速度則保持不變并一直等于當(dāng)?shù)芈曀?,這就是臨界狀態(tài).
(54)
式(52)與(54)是計算天然氣穩(wěn)態(tài)泄漏量的通用公式.
本文根據(jù)天然氣流動的基本控制方程,針對天然氣管道的正常和泄漏工況,利用特征線法分別推導(dǎo)了基本控制方程對應(yīng)的差分方程,給出了天然氣管道沿線瞬態(tài)壓力和流量的計算方法.根據(jù)天然氣流動的運動方程和能量方程推導(dǎo)出正常工況穩(wěn)態(tài)下的管道沿線壓力和溫度分布.推導(dǎo)研究了在泄漏工況穩(wěn)態(tài)下天然氣管道泄漏量的計算方法.