付建民 趙洪祥 陳國(guó)明 鄭曉云 朱 淵 任 婷

1.中國(guó)石油大學(xué)（華東）機(jī)電工程學(xué)院 2.中國(guó)石化工程建設(shè)有限公司3.中國(guó)石化石油工程地球物理有限公司河南分公司

隨著天然氣管網(wǎng)規(guī)模的不斷擴(kuò)大，由于各種原因?qū)е碌男孤┦鹿蕰r(shí)有發(fā)生。泄漏速率計(jì)算是天然氣管網(wǎng)安全領(lǐng)域的重要問(wèn)題，也是泄漏后擴(kuò)散、風(fēng)險(xiǎn)評(píng)估的前提和基礎(chǔ)。

較早的泄漏速率模型研究主要基于理論推導(dǎo)和數(shù)值計(jì)算，假設(shè)氣體在管內(nèi)為絕熱流動(dòng)，在泄漏點(diǎn)為等熵流動(dòng)，遵循理想氣體狀態(tài)方程和泊松方程，根據(jù)能量守恒和動(dòng)量守恒方程定量描述管內(nèi)氣體的流動(dòng)過(guò)程。Montiel［1］提出了適用于不同泄漏孔徑的管孔綜合模型，討論了管內(nèi)為亞臨界流、孔口為臨界流或亞臨界流以及管內(nèi)和孔口均為臨界流的泄漏速率計(jì)算方法。國(guó)內(nèi)學(xué)者［2－4］通過(guò)對(duì) Montiel模型的研究發(fā)現(xiàn)，管內(nèi)和孔口均出現(xiàn)臨界流的情形在工程實(shí)際中很少見(jiàn)。因此在氣體泄漏速率計(jì)算過(guò)程中僅保留前兩種情形。這類模型本身存在一定的缺陷，如沒(méi)有考慮裂口幾何形態(tài)對(duì)泄漏速率的影響，對(duì)于理論推導(dǎo)的泄漏速率缺少相應(yīng)的驗(yàn)證等。

實(shí)驗(yàn)和計(jì)算流體力學(xué)（CFD）可以更詳盡地研究管道泄漏的流場(chǎng)特征。趙金輝等［5］提出了大孔亞臨界流泄漏速率的計(jì)算公式，并搭建氣體管道泄漏模擬實(shí)驗(yàn)臺(tái)驗(yàn)證公式的正確性和適用性。Kostowski等［6］利用實(shí)驗(yàn)和CFD方法比較絕熱模型和等溫模型的適用性。目前，基于實(shí)驗(yàn)和CFD方法研究輸氣管道泄漏主要針對(duì)圓孔，考慮到不同的泄漏部位和破壞原因可能造成不同的裂口形狀和方向［7］。因此無(wú)法全面描述裂口幾何形態(tài)對(duì)泄漏速率的影響。

筆者的主要研究對(duì)象為中低壓輸氣管道因腐蝕等原因引發(fā)的小孔泄漏，這主要是因?yàn)橄鄬?duì)而言，小孔泄漏在實(shí)際情況中更難發(fā)現(xiàn)與定位［8］。利用相似原理搭建氣相管道小孔泄漏實(shí)驗(yàn)系統(tǒng)，設(shè)計(jì)并加工不同的泄漏模塊以模擬不同的泄漏場(chǎng)景；以FLUENT為平臺(tái)建立仿真模型，并對(duì)模型的可行性展開(kāi)全面的理論分析和實(shí)驗(yàn)驗(yàn)證，在此基礎(chǔ)上研究裂口幾何形態(tài)對(duì)泄漏速率的影響機(jī)理并得到泄漏孔口附近的氣體動(dòng)力學(xué)特征量，包括速度分布、馬赫數(shù)分布等，使研究趨于精細(xì)化，為管道泄漏事故的應(yīng)急救援和事故調(diào)查提供理論依據(jù)。

1 實(shí)驗(yàn)方法

1.1 實(shí)驗(yàn)系統(tǒng)

利用相似原理搭建氣相管道小孔泄漏實(shí)驗(yàn)系統(tǒng)，實(shí)驗(yàn)過(guò)程中以空氣為介質(zhì)，主要是出于安全和成本方面的考慮，盡管天然氣與空氣性質(zhì)不同，但并不會(huì)影響研究裂口幾何形態(tài)對(duì)泄漏速率的影響。

圖1為氣相管道小孔泄漏實(shí)驗(yàn)系統(tǒng)流程圖，系統(tǒng)主要由動(dòng)力單元、泄漏模塊與測(cè)量?jī)x表單元、數(shù)據(jù)采集單元組成。其中：

動(dòng)力單元：包括空氣壓縮機(jī)、緩沖罐、減壓閥，為系統(tǒng)提供壓力恒定的氣源。

泄漏模塊與測(cè)量?jī)x表單元：包括一套可拆卸更換的泄漏模塊以及安裝在泄漏模塊兩端的壓力表、渦街流量計(jì)，模擬泄漏過(guò)程發(fā)生及監(jiān)測(cè)過(guò)程參數(shù)的變化。

數(shù)據(jù)采集單元：包括PLC控制器、計(jì)算機(jī)，實(shí)現(xiàn)過(guò)程參數(shù)的實(shí)時(shí)采集與處理。

實(shí)驗(yàn)過(guò)程中，空氣通過(guò)壓縮機(jī)增壓至緩沖罐，經(jīng)減壓閥減壓后，在泄漏模塊處發(fā)生泄漏，未泄漏的空氣對(duì)大氣進(jìn)行放空；通過(guò)安裝在泄漏模塊兩端的壓力表和渦街流量計(jì)監(jiān)測(cè)壓力、流量參數(shù)的實(shí)時(shí)變化；通過(guò)泄漏模塊兩端渦街流量計(jì)的差值換算得到泄漏速率［9］。

圖1 氣相管道小孔泄漏實(shí)驗(yàn)系統(tǒng)流程圖

1.2 技術(shù)思路

實(shí)驗(yàn)通過(guò)調(diào)節(jié)減壓閥、更換泄漏模塊等手段構(gòu)建不同的泄漏場(chǎng)景，研究管內(nèi)壓力、泄漏孔口形狀、裂口方向?qū)π孤┧俾实挠绊憴C(jī)理。綜合考慮實(shí)驗(yàn)的可操作性及相關(guān)標(biāo)準(zhǔn)規(guī)范，確定具有代表性的孔口尺度與管道壓力范圍。為盡可能降低讀數(shù)的觀測(cè)誤差及舍入誤差，采用多次實(shí)驗(yàn)取平均值的方法。

1.3 管內(nèi)壓力

實(shí)驗(yàn)所用管道為不銹鋼材質(zhì)，內(nèi)徑50mm，壁厚3.5mm，壓力的選取主要依據(jù)《城鎮(zhèn)燃?xì)庠O(shè)計(jì)規(guī)范》（GB 50028—2006），結(jié)合減壓閥的工作特性，確定本實(shí)驗(yàn)的管道壓力范圍，如表1所示。

表1 管內(nèi)壓力范圍表

1.4 泄漏模塊

設(shè)計(jì)并加工一套可拆卸更換的泄漏模塊，即3根并聯(lián)連接的管道，每根管道長(zhǎng)1.2m，在管道頂部開(kāi)圓孔和矩形方孔，圓孔和矩形方孔具有等水力半徑的幾何特征，矩形方孔的長(zhǎng)邊沿周向或軸向分布，通過(guò)閥門(mén)切換管道模擬不同裂口形狀及方向的泄漏過(guò)程。圓孔孔徑為5mm，主要參照《基于風(fēng)險(xiǎn)的檢驗(yàn)》（API 581—2008）關(guān)于小孔泄漏孔徑的推薦范圍。圖2為泄漏模塊流程圖。

圖2 泄漏模塊流程圖

1.5 實(shí)驗(yàn)結(jié)果

圖3為圓孔、周向矩形方孔、軸向矩形方孔的泄漏速率實(shí)驗(yàn)值對(duì)比。由圖3可以看出，對(duì)于圓孔、周向矩形方孔、軸向矩形方孔，管內(nèi)壓力增大時(shí)泄漏速率隨之增加，呈近似線性關(guān)系，因?yàn)闅怏w通過(guò)孔流出的過(guò)程可假定為等熵膨脹，孔口面積一定，隨著管內(nèi)壓力增大，更多的壓力能轉(zhuǎn)化為氣體的動(dòng)能，泄漏速率隨之增加。相同的管內(nèi)壓力下，矩形方孔的泄漏速率顯著高于圓孔，從實(shí)驗(yàn)數(shù)據(jù)來(lái)看，裂口方向?qū)π孤┧俾实挠绊懖伙@著，壓力為0.245MPa時(shí)，周向矩形方孔與軸向矩形方孔泄漏速率的差值僅為0.032 3×10－3kg／s，但總體上，相同的管內(nèi)壓力下，周向矩形方孔泄漏速率略高于軸向矩形方孔。

圖3 不同裂口泄漏速率實(shí)驗(yàn)值對(duì)比圖

2 數(shù)值仿真

2.1 控制方程及計(jì)算方法

流體的流動(dòng)應(yīng)遵循質(zhì)量守恒定律、動(dòng)量守恒定律和能量守恒定律，而控制方程是上述守恒定律的數(shù)學(xué)描述。建立氣體管道泄漏數(shù)值模擬的控制方程組，包括可壓縮理想流體的質(zhì)量守恒方程、動(dòng)量守恒方程和能量守恒方程，其通用形式是［10］：

從左至右，方程中各項(xiàng)依次為瞬態(tài)項(xiàng)、對(duì)流項(xiàng)、擴(kuò)散項(xiàng)、源項(xiàng)。φ是通用變量；Γ是廣義擴(kuò)散系數(shù)；S是廣義源項(xiàng)。對(duì)于特定的方程，φ、Γ、S具有特定的形式。此外，控制方程組還應(yīng)包括附加的湍流輸運(yùn)方程。u是速度失量。

湍流模型采用標(biāo)準(zhǔn)k－ε雙方程模型，具有較高的適用性和精度。標(biāo)準(zhǔn)k－ε雙方程模型假定湍動(dòng)黏度是各向同性的，基本輸運(yùn)方程為：

式中Gk是由平均速度梯度引起的湍動(dòng)能k的產(chǎn)生項(xiàng)；Gb是由于浮力引起湍動(dòng)能k的產(chǎn)生項(xiàng)；YM代表可壓湍流中脈動(dòng)擴(kuò)張的貢獻(xiàn)；C1ε、C2ε和C3ε為經(jīng)驗(yàn)常數(shù)；αk和αε分別為湍動(dòng)能k和湍流耗散率ε對(duì)應(yīng)的Prandtl數(shù)；Sk和Sε是用戶定義的源項(xiàng)。

求解控制方程組采用基于壓力的穩(wěn)態(tài)求解器，選擇該方法是因?yàn)楸疚闹饕P(guān)心近似穩(wěn)定狀態(tài)下的泄漏速率和氣體動(dòng)力學(xué)特征量分布，對(duì)于中低壓輸氣管道，選擇基于壓力的求解器可提高計(jì)算效率同時(shí)可以保證較高的精度。

2.2 模型與網(wǎng)格劃分

模擬區(qū)域?yàn)榘霃?.6m，高3m的圓柱形區(qū)域，管道長(zhǎng)1.2m，內(nèi)徑50mm，壁厚3.5mm，在管道中心處分別開(kāi)圓孔、周向矩形方孔和軸向矩形方孔，具體尺度與相似實(shí)驗(yàn)相同?？紤]到管道長(zhǎng)度較短，認(rèn)為在靜風(fēng)條件下管道內(nèi)外流場(chǎng)具有軸對(duì)稱的性質(zhì)，因此對(duì)管道及模擬區(qū)域進(jìn)行簡(jiǎn)化。網(wǎng)格采用六面體單元、結(jié)構(gòu)化非均勻劃分方法，為準(zhǔn)確描述泄漏孔口附近的流場(chǎng)特征，對(duì)泄漏孔口附近的網(wǎng)格進(jìn)行局部加密，網(wǎng)格數(shù)為60×104，模型及網(wǎng)格如圖4所示。

圖4 模型與網(wǎng)格劃分圖

為了說(shuō)明網(wǎng)格密度對(duì)數(shù)值計(jì)算結(jié)果的影響，對(duì)模型進(jìn)行疏密兩種形式的網(wǎng)格劃分（網(wǎng)格單元數(shù)分別為40×104和60×104），流場(chǎng)計(jì)算結(jié)果表明，網(wǎng)格疏密程度對(duì)計(jì)算結(jié)果的影響很小，網(wǎng)格具有一定的無(wú)關(guān)性。

2.3 結(jié)果分析

2.3.1 可行性驗(yàn)證

保持出口壓力為0.101MPa不變，通過(guò)調(diào)節(jié)入口壓力得到不同入口壓力下泄漏速率的模擬值，并與實(shí)驗(yàn)值進(jìn)行對(duì)比，圖5為不同裂口泄漏速率的實(shí)驗(yàn)值與模擬值對(duì)比。

從圖5可以看出，泄漏速率的模擬值略高于實(shí)驗(yàn)值，但總體上較為一致，基于CFD方法研究裂口幾何形態(tài)對(duì)輸氣管道小孔泄漏的影響是可行的。從建模的角度來(lái)看，產(chǎn)生偏差的原因主要包括：①模型簡(jiǎn)化處理，建模過(guò)程中為提高計(jì)算效率同時(shí)又能保證一定的模擬精度，假設(shè)一段短管道內(nèi)的初始?jí)毫Φ戎?、均勻分布，?dǎo)致泄漏孔口的上游壓力高于實(shí)際值，即泄漏速率的模擬值略高于實(shí)驗(yàn)值；②管道壁面模型選擇，建模過(guò)程中管道壁面的熱通量設(shè)為0，即假設(shè)壁面為絕熱壁，Levenspiel指出倘若源處的溫度和壓力相同時(shí)，絕熱模型通常會(huì)高估實(shí)際流動(dòng)［11］；③粗糙度的選擇難以做到與實(shí)驗(yàn)所用管道完全相符，導(dǎo)致模擬值與實(shí)驗(yàn)值存在一定的偏差。

圖5 不同裂口泄漏速率實(shí)驗(yàn)值與模擬值對(duì)比圖

2.3.2 速度分布

圖6、7分別為出口壓力0.101MPa、入口壓力0.205MPa下不同裂口孔口附近速度矢量分布及截面速度分布。從圖6可以看出，由于膨脹作用，氣體的壓力能轉(zhuǎn)化為動(dòng)能，擴(kuò)散速度急劇增加，在孔口附近可以觀察到由于噴出氣體的高速運(yùn)動(dòng)而形成的負(fù)壓效應(yīng)，即周圍大氣以較低的速度向孔口運(yùn)動(dòng)的現(xiàn)象［12］。從圖7可以看出，孔口截面速度呈環(huán)形分布，最大速度位于孔口中心區(qū)域，靠近壁面的過(guò)程中，速度不斷減小直至為零，這主要是因?yàn)楸诿鏋闊o(wú)滑移邊界，切向速度為零。圓孔、周向矩形方孔、軸向矩形方孔的最大速度分別為315.81m／s、329.05m／s、329.22m／s，即相同管內(nèi)壓力下，矩形方孔的最大速度明顯高于圓孔，而裂口方向?qū)λ俣确植嫉挠绊懖伙@著。

圖6 不同裂口孔口附近速度矢量分布圖

圖7 不同裂口孔口截面速度分布圖

2.3.3 馬赫數(shù)分布

圖8為出口壓力0.101MPa下不同裂口中心線馬赫數(shù)分布。從圖8可以看出，馬赫數(shù)在很短的距離內(nèi)由0急劇增加到最大值，然后呈單調(diào)遞減。最大馬赫數(shù)發(fā)生在孔口截面，即Y為0.028 5m處。隨著管內(nèi)壓力的增大，相同位置的馬赫數(shù)也隨之增加。

圖8 不同裂口中心線馬赫數(shù)分布圖

從圓孔中心線馬赫數(shù)分布可以看出，當(dāng)壓力為0.215MPa時(shí)，孔口截面馬赫數(shù)為1，氣體在孔口處的流速為聲速，即臨界流；對(duì)于周向矩形方孔和軸向矩形方孔，當(dāng)壓力為0.205MPa時(shí)，氣體在孔口處達(dá)到臨界流。此時(shí)，出口壓力與管內(nèi)壓力的比值定義為臨界壓力比（CPR），即圓孔、矩形方孔的臨界壓力比分別為0.471和0.494，說(shuō)明孔口面積一定，矩形方孔比圓孔更容易在孔口處達(dá)到臨界流。對(duì)于空氣，臨界壓力比理論值的計(jì)算依據(jù)［13］：

式中pa為出口壓力，MPa；p2c為管內(nèi)臨界壓力，MPa；κ為絕熱指數(shù)，對(duì)于空氣，κ＝1.4。

由此可見(jiàn)，圓孔和矩形方孔的臨界壓力比均小于理論值，這主要是因?yàn)橛?jì)算過(guò)程中考慮了孔口處的摩擦損失［14］。對(duì)比圖8－b、c可知，周向矩形方孔、軸向矩形方孔中心線馬赫數(shù)分布基本重合，說(shuō)明裂口方向?qū)χ行木€馬赫數(shù)分布的影響不顯著。

3 結(jié)論

1）對(duì)于小孔泄漏，泄漏速率與管內(nèi)壓力呈近似正相關(guān)；管內(nèi)壓力一定，矩形方孔的泄漏速率明顯高于圓孔，而裂口方向?qū)π孤┧俾实挠绊懖伙@著。

2）基于實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證仿真模型的可行性，在此基礎(chǔ)上，得到泄漏孔口附近的氣體動(dòng)力學(xué)特征量分布，包括速度分布、馬赫數(shù)分布等，仿真結(jié)果表明，最大速度發(fā)生泄漏孔口截面中心處，矩形方孔的最大速度明顯高于圓孔，而裂口方向?qū)ζ溆绊懖伙@著；從臨界壓力比來(lái)看，孔口面積一定，矩形方孔更容易在孔口處達(dá)到臨界流，圓孔、周向矩形方孔、軸向矩形方孔的臨界壓力比模擬值均低于理論計(jì)算值，這主要是因?yàn)榭卓谔幍哪芰亢纳ⅰ?/p>

［1］MONTIEL H，VILCHEZ J A，CASAL J，et al.Mathematical modelling of accidental gas releases［J］.Journal of Hazardous Materials，1998，59（2／3）：211－233.

［2］霍春勇，董玉華，余大濤，等.長(zhǎng)輸管線氣體泄漏率的計(jì)算方法研究［J］.石油學(xué)報(bào)，2004，25（1）：101－105.HUO Chunyong，DONG Yuhua，YU Datao，et al.Estimation of accidental gas release flow rate in long transmission pipelines［J］.Acta Petrolei Sinica，2004，25（1）：101－105.

［3］董玉華，周敬恩，高惠臨，等.長(zhǎng)輸管道穩(wěn)態(tài)氣體泄漏率的計(jì)算［J］.油氣儲(chǔ)運(yùn)，2002，21（8）：11－15.DONG Yuhua，ZHOU Jing＇en，GAO Huilin，et al.Estimation of steady state gas release flow rate in long distance pipeline［J］.Oil＆ Gas Storage and Transportation，2002，21（8）：11－15.

［4］王大慶，霍春勇，高惠臨.長(zhǎng)輸管線氣體泄漏率簡(jiǎn)化計(jì)算方法［J］.天然氣工業(yè)，2008，28（1）：116－118.WANG Daqing，HUO Chunyong，GAO Huilin.A simplified method for calculating long－distance pipeline leakage rate［J］.Natural Gas Industry，2008，28（1）：116－118.

［5］趙金輝，譚羽非.燃?xì)夤芫€大孔亞臨界流泄漏實(shí)驗(yàn)［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2011，43（2）：84－87.ZHAO Jinhui，TAN Yufei.Experiment study on subcritical flow nozzle model in the gas pipelines［J］.Journal of Harbin Institute of Technology，2011，43（2）：84－87.

［6］KOSTOWSKI W J，SKOREK J.Real gas flow simulation in damaged distribution pipelines［J］.Energy，2012，45（1）：481－488.

［7］張一先，張隆.管道煤氣泄漏事故評(píng)估的不確定性［J］.煤氣與熱力，2000，20（1）：9－13.ZHANG Yixian，ZHANG Long.Uncertainties on assessment of pipeline gas leakage accidents［J］.Gas ＆ Heat，2000，20（1）：9－13.

［8］謝昱姝，吳宗之，呂良海，等.城市管道天然氣在土壤中泄漏擴(kuò)散實(shí)驗(yàn)研究［J］.中國(guó)安全生產(chǎn)科學(xué)技術(shù)，2012，8（4）：13－17.XIE Yushu，WU Zongzhi，LYU Lianghai，et al.Experimental research on diffusion behavior of leaked gas from underground gas pipeline［J］.Journal of Safety Science and Technology，2012，8（4）：13－17.

［9］趙洪祥，付建民，陳國(guó)明，等.氣相管道小孔泄漏壓力響應(yīng)試驗(yàn)與仿真研究［J］.安全與環(huán)境學(xué)報(bào)，2014，14（1）：47－50.ZHAO Hongxiang，F(xiàn)U Jianmin，CHEN Guoming，et al.Experimental and simulated study on the pressure response to the gas－pipeline leaking［J］.Journal of Safety and Environment，2014，14（1）：47－50.

［10］劉振翼，張應(yīng)安，鄭遠(yuǎn)攀，等.含CO2天然氣井井噴流場(chǎng)特征［J］.石油勘探與開(kāi)發(fā)，2011，38（1）：90－96.LIU Zhenyi，ZHANG Ying′an，ZHENG Yuanpan，et al.Flow fields in the blowout of natural gas wells containing carbon dioxide［J］.Petroleum Exploration and Development，2011，38（1）：90－96.

［11］LEVENSPIEL O.Engineering flow and heat exchange［M］.New York：Plenum Press，1984.

［12］鄭遠(yuǎn)攀，蘇曉珂，劉新新，等.氣相井噴的泄漏源模型及其仿真［J］.高壓物理學(xué)報(bào)，2012，26（2）：189－198.ZHENG Yuanpan，SU Xiaoke，LIU Xinxin，et al.Study on leakage model of gas blowout and its numerical simulation［J］.Chinese Journal of High Pressure Physics，2012，26（2）：189－198.

［13］DANIEL A C，JOSEPH F L.Chemical process safetyfundamentals with application［M］.New Jersey：Prentice－Hall，1990.

［14］丁英濤，謝君堂，仲順安.微型噴管內(nèi)氣體流動(dòng)的流量壅塞現(xiàn)象［J］.北京理工大學(xué)學(xué)報(bào)，2006，26（12）：1086－1089.DING Yingtao，XIE Juntang，ZHONG Shun′an.Phenomenon of gas flow choking in micro－nozzles［J］.Transactions of Beijing Institute of Technology，2006，26（12）：1086－1089.