王少雄，李玉星，劉翠偉，梁杰，李安琪，薛源

（中國石油大學(xué)（華東）山東省油氣儲(chǔ)運(yùn)安全省級(jí)重點(diǎn)實(shí)驗(yàn)室，山東青島266580）

引 言

天然氣的管道輸送過程中不可避免地會(huì)遇到水體環(huán)境、河流沖刷、水體對(duì)管道的腐蝕，以及設(shè)計(jì)安裝不當(dāng)、第三方破壞等人為因素，使得水下輸氣管道面臨的風(fēng)險(xiǎn)增大，泄漏事故也逐年增加[1-3]。當(dāng)天然氣連續(xù)地從泄漏孔進(jìn)入水體中時(shí)，氣體由于自身與周圍水域的壓力差以及氣體與水體交界面表面張力的作用，將會(huì)破裂成一個(gè)個(gè)氣泡。氣泡受到浮力的作用而上升，上升的氣泡將會(huì)夾帶周圍的海水一起向上運(yùn)動(dòng)，從而形成氣泡羽流[4]。水下輸氣管道一旦發(fā)生泄漏，上升到水面的氣泡羽流可能會(huì)引起火災(zāi)、造成水域污染，甚至?xí)T發(fā)二次事故對(duì)人身財(cái)產(chǎn)安全造成威脅。因此，科學(xué)地認(rèn)識(shí)水下輸氣管道泄漏所形成的氣泡羽流擴(kuò)散規(guī)律，對(duì)于及時(shí)地采取防治措施和制定風(fēng)險(xiǎn)評(píng)價(jià)方案具有一定的意義。

目前對(duì)于輸氣管道泄漏的研究大多針對(duì)陸地的架空管道[5-7]，對(duì)于水下輸氣管道泄漏的模擬研究主要包括積分模型和計(jì)算流體動(dòng)力學(xué)模型。積分模型把氣泡羽流中的氣泡和卷吸進(jìn)來的水體看作一個(gè)多相控制體，依據(jù)相似性原理，假設(shè)羽流徑向的質(zhì)量和動(dòng)量守恒為某個(gè)概率分布從而得到特征參數(shù)隨羽流軌跡變化的常微分方程。Morton 等[8-10]提出了一種積分模型，該模型假設(shè)羽流速度和密度分布具有相似的形式，這一假設(shè)簡化了羽流軸向上的質(zhì)量和動(dòng)量守恒的計(jì)算方法。Ditmars等[11-12]進(jìn)一步完善了Morton 等的模型，在模型中考慮了氣泡之間的滑移速度，并且認(rèn)為氣泡羽流對(duì)水體產(chǎn)生了拖動(dòng)作用。Milgram 等[13-14]提出了自己的積分模型，并且在估計(jì)氣泡羽流的動(dòng)量時(shí)引入了動(dòng)量放大系數(shù)的概念。Billeter等[15]通過假設(shè)“反向停滯流”改進(jìn)了積分模型，并且提出了新的預(yù)測羽流輪廓的方法。計(jì)算流體動(dòng)力學(xué)模型是通過求解流體流動(dòng)的Navier-Stokes 方程來研究流體動(dòng)力學(xué)問題，它直接模擬氣泡夾帶和湍流運(yùn)輸產(chǎn)生的影響。Cloete 等[16]建立了CFD 數(shù)值模型來研究水下輸氣管道斷裂引起的泄漏擴(kuò)散。Olsen 等[17]提出了一種基于歐拉-拉格朗日準(zhǔn)則的CFD 模型來模擬開闊水域的氣泡羽流的釋放，并將模型與Milgram 等[13-14]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證，發(fā)現(xiàn)模型具有良好的適應(yīng)性。Tessarolo 等[18]提出了一種基于拉格朗日法的數(shù)值模型來模擬深水井噴的油氣擴(kuò)散，并通過模型驗(yàn)證了三種夾帶關(guān)系式，發(fā)現(xiàn)JETLAG 關(guān)系式與實(shí)驗(yàn)結(jié)果最為吻合。Wu 等[19]用CFD 模型來模擬海底氣體釋放，并且評(píng)估了四種數(shù)值方法（雷諾時(shí)均N-S、非穩(wěn)態(tài)雷諾時(shí)均N-S、LES 大渦模擬、SAS 尺寸自適應(yīng)模擬）在捕捉上升氣泡羽流的特征行為上的適用性。國內(nèi)學(xué)者也針對(duì)天然氣在水體中的泄漏擴(kuò)散進(jìn)行了相關(guān)的模擬研究。李長俊等[20]采用CFD建模的方法研究了穿越河道天然氣管道發(fā)生泄漏時(shí)不同的泄漏朝向和水流速度下氣體的水平遷移距離變化規(guī)律。李新宏等[21-22]利用Fluent 建立了海底管道泄漏二維擴(kuò)散模型，研究了不同海流流速和泄漏孔徑對(duì)氣體在水體中的擴(kuò)散形態(tài)和泄漏速率的影響。張煥好等[23]建立了二維軸對(duì)稱水下超聲速氣體射流的數(shù)值計(jì)算模型,重點(diǎn)研究了水下超聲速氣體射流的初始流動(dòng)結(jié)構(gòu)。相比于模擬研究，國內(nèi)外專家對(duì)于氣體水下淹沒射流進(jìn)行了大量實(shí)驗(yàn)研究，并取得了豐富成果。Topham[24]在一個(gè)水深為23 m 的水池中做了現(xiàn)場實(shí)驗(yàn)，測量了氣泡羽流在軸線上和橫斷面上的速度變化規(guī)律。Lees 等[25]進(jìn)行了現(xiàn)場和實(shí)驗(yàn)室實(shí)驗(yàn)來研究水下氣體釋放時(shí)的氣體濃度分布，發(fā)現(xiàn)氣體濃度呈現(xiàn)高斯分布。王志剛等[26]進(jìn)行了室內(nèi)氣體水下泄漏實(shí)驗(yàn)，研究了不同工況下的羽流直徑的變化規(guī)律。郭強(qiáng)等[27]進(jìn)行了水下超聲速氣體射流實(shí)驗(yàn)，研究了不同工況下的水下高速氣體垂直射流的演化過程。王超等[28-29]通過實(shí)驗(yàn)研究了超聲速氣體射流在水中的噴射所形成的流場結(jié)構(gòu)，并用數(shù)值計(jì)算的方法成功模擬了射流初期氣泡運(yùn)動(dòng)演化的復(fù)雜過程。施紅輝等[30]進(jìn)行了三維水下超聲速氣體射流實(shí)驗(yàn)，重點(diǎn)研究了射流初期氣泡運(yùn)動(dòng)演化過程。李婷婷等[31]進(jìn)行了水下豎直環(huán)形噴管出口氣體射流實(shí)驗(yàn)，并通過Matlab 對(duì)環(huán)形噴管射流的特性進(jìn)行了分析。

上述模擬研究大多側(cè)重于特定的管路過程參數(shù)對(duì)氣體的流動(dòng)特性的影響而對(duì)氣體上升到水面時(shí)的擴(kuò)散特性以及氣體在水體中具體的擴(kuò)散過程研究較少，并且目前還沒有關(guān)于兩種模型對(duì)水下輸氣管道泄漏擴(kuò)散特性的預(yù)測結(jié)果準(zhǔn)確性的對(duì)比評(píng)估。因此本文通過建立水下輸氣管道泄漏三維CFD 模型和積分?jǐn)?shù)學(xué)模型來研究水下輸氣管道泄漏后的擴(kuò)散特性，分析水下輸氣管道發(fā)生泄漏時(shí)形成的氣泡羽流在水體中擴(kuò)散的速度、羽流半徑和泉涌高度的變化規(guī)律，并利用實(shí)驗(yàn)數(shù)據(jù)對(duì)兩種模型的準(zhǔn)確性進(jìn)行了對(duì)比驗(yàn)證，為水下輸氣管道發(fā)生泄漏后的風(fēng)險(xiǎn)評(píng)價(jià)以及后續(xù)的模擬研究提供一定的參考。

1 水下氣體泄漏CFD模型理論基礎(chǔ)

歐拉-拉格朗日模型可以用來求解水下輸氣管道泄漏后的多相流動(dòng)過程，該模型將泄漏氣體作為拉格朗日離散相而周圍的水體作為連續(xù)的歐拉相處理，這些離散相粒子在連續(xù)相中運(yùn)動(dòng)，并通過交換質(zhì)量、動(dòng)量和能量與之相互作用。歐拉參考系中的方程組求解周圍水體的速度和湍流、可溶性氣體的濃度以及不同連續(xù)相（水和大氣）之間的界面位置，這些變量將用于對(duì)拉格朗日系中氣泡運(yùn)動(dòng)軌跡的追蹤。氣泡上升過程中在浮力、拖曳力和虛擬質(zhì)量力的作用下達(dá)到平衡，拖曳力和虛擬質(zhì)量力也會(huì)對(duì)歐拉相進(jìn)行耦合從而影響周圍水體的速度，因此歐拉-拉格朗日模型具有強(qiáng)大的雙向耦合功能，所以又稱為耦合的DPM 和VOF模型。

1.1 VOF模型

在VOF 模型中，水體上方的大氣和周圍水體都被視為連續(xù)相，二者相互貫穿，一相的體積不能被另外一相的體積所占有。模型中采用體積率來衡量氣體在水中所占的體積比，氣體和水的體積之和為1。使用VOF 方法，首先要定義一種流體體積函數(shù)，通過網(wǎng)格內(nèi)該流體的體積與網(wǎng)格體積的比值來確定自由面，繼而追蹤流體的變化。對(duì)于水下天然氣管道泄漏，VOF 模型通過求解連續(xù)性方程、動(dòng)量方程和能量方程來實(shí)現(xiàn)對(duì)水和大氣兩相自由表面的追蹤[2]。

連續(xù)性方程如式（1）所示。

式中，αi為第i相歐拉相的體積分?jǐn)?shù)；ρi為第i相歐拉相的密度，kg/m3；vi為第i相歐拉相的速度，m/s。

動(dòng)量方程取決于各相的體積分?jǐn)?shù)，如式（2）所示。

式中，ρ = ∑αiρi為混合相的密度，kg/m3；μ =μT+ μM為分子混合黏度和湍流黏度之和，(N·s)/m2；F為外部作用力，N。

能量方程如式（3）所示。

式中，keff表示流體的傳熱系數(shù)；Sh為源項(xiàng)。

1.2 DPM模型

DPM 模型將泄漏氣體作為離散相處理，假設(shè)氣體粒子在拉格朗日計(jì)算區(qū)域中的體積分?jǐn)?shù)低于10%～12%，通過對(duì)離散相粒子施加一個(gè)平衡力實(shí)現(xiàn)對(duì)粒子運(yùn)動(dòng)軌跡的追蹤。氣體受力平衡方程如式（4）所示。

式中，ub表示泄漏氣體的速度，m/s；ρb表示泄漏氣體的密度，kg/m3；FD表示氣泡顆粒受到水的拖曳力，N；Cd表示拖曳力系數(shù)；E0為無量綱數(shù)，用來表征氣泡上升過程中形狀的變化；σ表示環(huán)境水的黏度，(N·s)/m2；Fvm表示虛擬質(zhì)量力，N；Cvm表示虛擬質(zhì)量力系數(shù)，通常為0.5；db表示氣泡顆粒的直徑，m；Re表示相對(duì)Reynolds數(shù)。

1.3 氣泡密度和尺寸模型

隨著氣泡在水體中的上升，氣泡所處水深變淺，使得氣泡運(yùn)動(dòng)過程中周圍的壓力發(fā)生顯著變化，從而導(dǎo)致氣泡密度隨水深的降低而降低。因此，水下氣體的密度是水深的函數(shù)，可以通過推導(dǎo)理想氣體方程進(jìn)行求解，如式（9）所示。

式中，Mb表示泄漏氣體的相對(duì)分子質(zhì)量；R 為通用氣體常數(shù)，J/(kmol·K)；P 表示氣泡所承受的靜水壓力，Pa；P0表示大氣壓，Pa；Hp表示氣泡所處的水深，m。

由式（7）可知，氣泡的大小決定了氣泡的特征形狀，進(jìn)而影響了氣泡在水中所受拖曳力的大小。同時(shí)氣泡上升過程中與環(huán)境水存在質(zhì)量傳遞，而這取決于氣泡的表面積大小。現(xiàn)有的模型出于簡化計(jì)算過程大多將氣泡大小設(shè)定為常數(shù)，但氣泡上升時(shí)由于發(fā)生破裂和合并使得氣泡的大小時(shí)刻發(fā)生變化，因此現(xiàn)有的模型未能準(zhǔn)確地預(yù)測氣泡羽流的擴(kuò)散特性。Laux 等[32]通過求解輸運(yùn)方程建立了一個(gè)歐拉框架下的氣泡尺寸模型，如式（11）所示。

式中，τrel為松弛時(shí)間，s；為氣泡的平衡直徑，m；αb為氣泡的體積分?jǐn)?shù)；σ1為表面張力；C1= 4 為無量綱常數(shù)；C2=100 μm表示氣泡的最小尺寸。

1.4 網(wǎng)格模型和數(shù)值計(jì)算方法

由于本文模擬的泄漏工況為水下輸氣管道的小孔泄漏，F(xiàn)luent 中的二維網(wǎng)格會(huì)將泄漏孔默認(rèn)為狹縫泄漏，不能真實(shí)反映泄漏孔徑對(duì)天然氣泄漏擴(kuò)散的影響，而三維網(wǎng)格可將泄漏孔形狀畫為圓形，所以本文對(duì)水下輸氣管道泄漏擴(kuò)散模型進(jìn)行簡化，泄漏孔位于模型底部的中心位置。以水深h=1.1 m，泄漏孔徑d=3 mm 為例，模型所研究的區(qū)域大小為1 m×1 m×1.5 m，由于泄漏孔徑相對(duì)較小，對(duì)泄漏口附近的網(wǎng)格進(jìn)行加密處理，先對(duì)整體block 劃分OBlock 網(wǎng)格，再對(duì)中心的block 進(jìn)行O-Block 網(wǎng)格劃分，泄漏孔向模型邊界的網(wǎng)格線劃分方式為Exponential2，Spacing2=0.0002，Ratio2=1.135，網(wǎng) 格數(shù)為376000個(gè)，分別選擇Determinant2×2×2和Angle作為網(wǎng)格質(zhì)量的判定標(biāo)準(zhǔn)，所有網(wǎng)格的Determinant 2×2×2值大于0.7，大于0.85的占93.832%，所有網(wǎng)格的Angle 值大于40.5°，大于60°的占70.331%，生成的網(wǎng)格如圖1所示。

圖1 網(wǎng)格模型劃分Fig.1 Mesh model

計(jì)算域通過Patch 選項(xiàng)將其分為水域和空氣域，空氣域上部設(shè)為壓力出口邊界，計(jì)算域的側(cè)面和底部設(shè)為無滑移固壁邊界，底部中心半徑1.5 mm的泄漏孔為空氣入口，泄漏孔的初值條件在DPM 模型中設(shè)置，在DPM 模型中將入口邊界條件設(shè)置為流量入口，該值與實(shí)驗(yàn)中的泄漏速率相對(duì)應(yīng)。水下輸氣管道的泄漏是一種瞬態(tài)湍流過程，因此模型采用k-ε 湍流模型和非穩(wěn)態(tài)壓力基求解器以及適用于求解非穩(wěn)態(tài)流動(dòng)的PISO 算法。同時(shí)由于氣體的可壓縮性以及與水體之間的密度差，所以在VOF 模型中選用隱式體積力公式。對(duì)于氣泡上升過程中的破裂和合并通過編寫UDF 函數(shù)導(dǎo)入氣泡密度和尺寸模型。為了對(duì)網(wǎng)格進(jìn)行獨(dú)立性檢驗(yàn)，通過設(shè)置最大網(wǎng)格尺寸、網(wǎng)格層數(shù)、第一層網(wǎng)格點(diǎn)與端點(diǎn)之間的距離等參數(shù)，得到不同網(wǎng)格劃分方案，具體方案如表1所示。圖2是水深1.1 m、泄漏孔徑為3 mm 時(shí)不同網(wǎng)格數(shù)下羽流上升時(shí)間與實(shí)驗(yàn)值的對(duì)比，從圖中可以看出，網(wǎng)格數(shù)為37.6 萬個(gè)和102.9 萬個(gè)對(duì)羽流的上升時(shí)間影響不大，網(wǎng)格數(shù)越多越接近實(shí)驗(yàn)值，為了加快計(jì)算速度并考慮模型的準(zhǔn)確性，選擇網(wǎng)格數(shù)為37.6 萬個(gè)的模型來模擬氣體在水體中的泄漏擴(kuò)散過程。

表1 網(wǎng)格劃分方案Table 1 Mesh classification scheme

圖2 不同數(shù)目網(wǎng)格下的羽流上升時(shí)間Fig.2 Plume rise time under different numbers of grids

2 水下氣體泄漏積分模型理論基礎(chǔ)

2.1 積分模型推導(dǎo)過程

當(dāng)水下輸氣管道發(fā)生泄漏時(shí)，泄漏氣體向上連續(xù)地進(jìn)入水體中形成氣泡羽流。積分模型假設(shè)氣泡羽流截面上的速度和空隙率時(shí)均分布為高斯分布[31]，分別如式（13）和式（14）所示。

式中，v 為羽流的速度，m/s；下角標(biāo)c 表示軸線處的值；r 和z 分別為距離羽流軸心的水平距離和距離泄漏點(diǎn)處的垂直高度，m；b 為羽流軸心到邊緣的寬度，m；φ 為羽流的空隙率；λ 為空隙率分布和速度分布之比。

氣相的連續(xù)性方程如下

式中，g 為重力加速度；Qg為氣體的體積流量，m3/s；Hv為水深，m；Hp表示與大氣壓力相對(duì)應(yīng)的水深，在清水條件下，該值為10.33 m，在海水條件下，該值在標(biāo)準(zhǔn)大氣壓下稍微低于10 m；γ 是動(dòng)量放大系數(shù)；s～ =(1+ λ2)v～s為表征氣體上升過程中滑移速度影響的無量綱量；z～、b～、v～ 分別為無量綱軸向坐標(biāo)、無量綱羽流寬度和無量綱軸向速度，表達(dá)式如式（17）、式（18）、式（19）所示。

液相的連續(xù)性方程和氣液混合物的動(dòng)量方程如式（20）、式（21）所示。

在該簡化的積分模型中，上述兩個(gè)方程的近似解可由式（22）、式（23）表示

作用在羽流上的力包括水面上的大氣壓力Pa、底部的壓力Pb、靜止水面下水的重力Gw、泉涌的重力Gf和浮力B，如圖3所示。氣體上升到水面形成的泉涌高度可以通過動(dòng)量守恒方程來求解，積分模型假設(shè)控制體在橫向是無限大的，所以沒有垂直動(dòng)量通量通過外部邊界。平衡方程如式（24）所示

圖3 作用在羽流上的力Fig.3 Force acting on plume

令控制體的最低邊界位于泄漏孔下方，底部的靜水壓力Pb將不會(huì)影響羽流的流動(dòng)，那么Pb的值將是恒定的，因此式（24）中的項(xiàng)Pa、Pb和Gw可以相互抵消，簡化后的方程如下

假設(shè)泉涌的密度ρf和靜水面處羽流的密度ρ(Hv)相等以及泉涌底部的動(dòng)能完全轉(zhuǎn)化為勢(shì)能，那么式（26）可以簡化為式（27）。

式中，β是經(jīng)驗(yàn)常數(shù)；v = vc(Hv)是泉涌底部軸線處羽流的速度，m/s；hf是泉涌高度,m。

則羽流徑向泉涌高度分布的表達(dá)式如式（28）所示。

式中，hoffset是高斯分布基準(zhǔn)線在靜水面上的偏移量。

2.2 積分模型中的經(jīng)驗(yàn)系數(shù)

2.2.1 多變指數(shù)n 假設(shè)氣泡羽流在上升過程中氣體是等溫膨脹的，則多變指數(shù)n=1，因?yàn)槿绻窃诮^熱情況下膨脹將導(dǎo)致上升氣體的溫度大幅度下降，這與實(shí)際情況不符。

2.2.2 夾帶系數(shù)α 夾帶系數(shù)α 隨著氣體流速的增加而增大，Milgram等[14]在1983年提出了一個(gè)關(guān)于夾帶系數(shù)的半經(jīng)驗(yàn)公式

其中，k=0.165；A=7.598；氣泡Froude 數(shù)Fr 由式（30）表示

在實(shí)驗(yàn)室尺度下，夾帶系數(shù)的值通常較小，因?yàn)閷?shí)驗(yàn)室條件下氣體流量較小，而現(xiàn)場實(shí)驗(yàn)因氣體泄漏速度較大得出的夾帶系數(shù)的值也相對(duì)較大。

2.2.3 寬度比λ 寬度比λ 的范圍較小，通常處于0.6～1之間，并且對(duì)羽流的影響較小。實(shí)驗(yàn)室尺度下得出的寬度比λ 的值通常比現(xiàn)場實(shí)驗(yàn)得出的值要小。

2.2.4 滑移速度vs在氣泡羽流中，氣泡相對(duì)于液體的滑移速度是氣體濃度、氣泡直徑等多因素的函數(shù)，并且由于氣泡是以氣泡群的形式存在，所以氣泡間的受力較單個(gè)氣泡的受力復(fù)雜。但氣泡滑移速度對(duì)羽流的影響很小，Levich[33]在1962 年提出當(dāng)氣泡直徑處于0.2～1.5 cm 之間時(shí)滑移速度為28～30 cm/s,而對(duì)于直徑更大的氣泡滑移速度將增長到35～40 cm/s，但是在湍流羽流中大氣泡通常會(huì)破碎成較小尺寸的氣泡。

2.2.5 動(dòng)量放大系數(shù)γ 動(dòng)量放大系數(shù)是總的動(dòng)量通量與平均流動(dòng)所引起的動(dòng)量通量之比，當(dāng)羽流中的氣泡與羽流的尺寸相比變得非常小時(shí)，氣泡的動(dòng)力以及氣泡之間的相互作用就變得沒那么重要，羽流流動(dòng)表現(xiàn)為單相流，動(dòng)量放大系數(shù)將趨向統(tǒng)一。因此，在求解上述羽流控制方程時(shí)，可以忽略該參數(shù)的影響。

表2 推薦的經(jīng)驗(yàn)參數(shù)Table 2 Recommended empirical parameters

由上可知，羽流的形成與發(fā)展和夾帶系數(shù)密切相關(guān)，其他參數(shù)的變化對(duì)結(jié)果影響很小，各個(gè)參數(shù)的范圍和推薦值如表2所示。

對(duì)積分模型進(jìn)行數(shù)值求解，模型計(jì)算程序框圖如圖4所示。

圖4 模型算法Fig.4 Model algorithm

3 模型驗(yàn)證

為了驗(yàn)證CFD 模型和積分模型的準(zhǔn)確性，采用水下泄漏實(shí)驗(yàn)環(huán)道裝置，如圖5所示，以空氣作為實(shí)驗(yàn)介質(zhì)，進(jìn)行水下輸氣管道泄漏實(shí)驗(yàn)，研究不同泄漏速率下氣體在水體中以及水面上的擴(kuò)散特性。實(shí)驗(yàn)是在一個(gè)長、寬為1 m，高1.5 m 的水槽中進(jìn)行的，環(huán)道總長為251.5 m，內(nèi)徑為42 mm。泄漏孔前后裝有壓力表7、9（量程為0～1 MPa，測量精度為0.5%）以及質(zhì)量流量計(jì)5、10（量程為0.01～1 m3/h，測量精度為0.2%）。水槽前方布置高速攝像機(jī)（型號(hào)為FASTCAM SA-X2，拍攝速度為5000 幀/秒），可通過相機(jī)架桿調(diào)節(jié)高度，對(duì)氣體在水體中的擴(kuò)散過程進(jìn)行拍攝，高速攝像機(jī)與PC 端連接，實(shí)時(shí)保存泄漏擴(kuò)散過程。通過專業(yè)的視頻圖像處理軟件ProAnalyst 對(duì)羽流擴(kuò)散過程進(jìn)行分析，將實(shí)驗(yàn)數(shù)據(jù)與模型進(jìn)行對(duì)比，結(jié)果討論如下。

圖5 水下管道氣體泄漏實(shí)驗(yàn)系統(tǒng)Fig.5 Experimental system for gas leakage of underwater pipeline

3.1 水下擴(kuò)散過程分析

圖6 為水深30 cm、泄漏孔徑5 mm、泄漏壓力300 kPa、泄漏速率0.75 m3/h時(shí)CFD 模擬的羽流擴(kuò)散過程和實(shí)驗(yàn)的對(duì)比。由圖可知，水下輸氣管道的泄漏過程可以分為三個(gè)不同的階段，每個(gè)階段都對(duì)應(yīng)著一個(gè)流動(dòng)區(qū)域且在每個(gè)階段羽流的上升動(dòng)力以及所產(chǎn)生的物理現(xiàn)象均不相同。首先是初始階段，該階段的羽流為動(dòng)量射流，對(duì)應(yīng)為流動(dòng)形成區(qū)。氣體從孔口泄漏到水體環(huán)境中，起始會(huì)在孔口形成一個(gè)完整的小氣泡，氣泡由最初一個(gè)突起逐漸上升變大，形成半球。隨著泄漏過程的進(jìn)行，小氣泡逐漸上升拉長變成了橢球狀，具體形態(tài)如圖6(a)所示。其次是充分發(fā)展階段。氣體將從流動(dòng)形成區(qū)進(jìn)入定型流動(dòng)區(qū)，如圖6(b)所示。在這個(gè)區(qū)域內(nèi)氣體的自身浮力將取代初始動(dòng)量占據(jù)主導(dǎo)作用，初始動(dòng)量對(duì)羽流產(chǎn)生的影響變得微乎其微。氣體在以管內(nèi)壓力作為初始動(dòng)力沖出泄漏孔口進(jìn)入到水體環(huán)境之后，克服了周圍環(huán)境的阻力開始快速地上升和膨脹。氣體在上升過程中，氣泡羽流速度間斷面的不穩(wěn)定會(huì)引起湍動(dòng)，從而把周圍靜止的水體卷入到向上運(yùn)動(dòng)的羽流當(dāng)中。隨著湍動(dòng)的發(fā)展，被卷吸進(jìn)入羽流內(nèi)部的流體不斷增多并且隨著羽流一起向上運(yùn)動(dòng)，羽流的邊界不斷向兩側(cè)發(fā)展擴(kuò)張。最后是表面流動(dòng)階段，隨著氣泡羽流的上升，氣泡羽流所處的水深也就越淺，此時(shí)所受周圍的環(huán)境壓力也就越小。環(huán)境水被夾帶到羽流中并被攜帶到表面，當(dāng)羽流到達(dá)表面時(shí)，夾帶的水并不會(huì)隨著氣體進(jìn)入大氣而是從水平表面向外徑向流動(dòng)，帶動(dòng)氣泡遠(yuǎn)離羽流軸心，如圖6(c)所示。大部分氣體會(huì)繼續(xù)上升進(jìn)入到大氣中，形成表面紊動(dòng)或沸騰，稱為表面流動(dòng)區(qū)。

圖6 羽流擴(kuò)散過程Fig.6 Plume diffusion process

圖7 不同水深下的羽流擴(kuò)散過程Fig.7 Plume diffusion process at different water depths

圖8 不同水深下的羽流半徑分布Fig.8 Plume radius distribution under different water depths

當(dāng)水下輸氣管道發(fā)生泄漏時(shí)，水深（泄漏點(diǎn)距離水面的高度）是泄漏氣體在水下擴(kuò)散時(shí)最重要也是最敏感的影響參數(shù)之一。圖7 是泄漏孔徑為3 mm、泄漏壓力為200 kPa、水深分別是1.5 m 和1 m時(shí)CFD 模擬的羽流擴(kuò)散過程。由圖可知，1.5 m 水深下氣泡羽流的擴(kuò)散過程和水深1 m 時(shí)的擴(kuò)散過程基本一致，但相比于水深1.5 m 的泄漏工況水深1 m時(shí)初始動(dòng)能的作用更加顯著，受到水體的阻礙作用相對(duì)較小，氣泡羽流的平均上升速度更快，因此到達(dá)水面的時(shí)間也更短。圖8 是泄漏孔徑為3 mm、水深分別是1.5 m和1 m時(shí)CFD模擬的羽流半徑分布，從圖中可以看出水深1 m 時(shí)的羽流半徑稍大于水深1.5 m 時(shí)相同水位高度的羽流半徑。這是因?yàn)樵谙嗤男孤毫托孤┛讖较?，水深?huì)改變泄漏孔外的水壓。水深增加，水壓變大，此時(shí)氣體沖出孔口就受到更大的阻力，氣體在水體中的擴(kuò)張也受到了更大的阻礙，因此在水深較大時(shí)羽流半徑就會(huì)隨之減小。

隨著氣泡羽流的上升，氣泡羽流所處的水深也就越淺，此時(shí)所受周圍的環(huán)境壓力也就越小。周圍水被夾帶到羽流中并被攜帶到表面，當(dāng)羽流到達(dá)表面時(shí)，夾帶的水并不會(huì)隨著氣體進(jìn)入大氣而是從水平表面向外徑向流動(dòng)，帶動(dòng)氣泡遠(yuǎn)離羽流軸心，形成一個(gè)圓形氣池，如圖9 所示。圖10 是泄漏孔徑為3 mm、水深分別是1.5 m和1 m時(shí)氣池半徑隨時(shí)間的變化。當(dāng)羽流到達(dá)水面時(shí)，以氣體溢出點(diǎn)為中心繼續(xù)向外擴(kuò)散，擴(kuò)散半徑持續(xù)增大，當(dāng)?shù)竭_(dá)一定程度時(shí)將會(huì)趨于穩(wěn)定。從圖中還可以看出水深越深，氣池?cái)U(kuò)散半徑越大。

圖9 不同水深下的氣池?cái)U(kuò)散過程Fig.9 Gas pool diffusion process at different water depths

圖10 不同水深下的氣池半徑隨時(shí)間的變化曲線Fig.10 Curves of gas pool radius with time under different water depths

3.2 羽流速度變化規(guī)律

在實(shí)驗(yàn)中，通過高速攝像機(jī)記錄輸氣管道泄漏發(fā)生的過程，將拍攝的視頻分幀解析成圖片，用ProAnalyst 圖像軟件處理連續(xù)幀的圖片。測量出每幀圖片中羽流寬度和高度，用連續(xù)幀寬度和高度的差值除以每幀間隔時(shí)間即為單位時(shí)間內(nèi)羽流的速度，取多次實(shí)驗(yàn)均值進(jìn)行統(tǒng)計(jì)，將CFD 模型和積分模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，如圖11、圖12所示。

圖11 為氣泡羽流在水深1.1 m、泄漏速率為0.21 m3/h 時(shí)，水位高度在0.5、0.7 和0.9 m 處的徑向速度分布。為了避免分析和闡述上的混亂，本文引入無量綱水位高度hz的概念（水位高度與水深的比值），將不同水位高度與水深區(qū)分開來。由圖可知，在較低水位高度（hz為0.45 m 和0.64 m）處，CFD 模型和積分模型計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果基本相同，羽流的徑向速度由于夾帶作用隨著徑向距離的增加逐漸降低。但是，在水位高度為0.9 m 處，積分模型計(jì)算結(jié)果與實(shí)驗(yàn)存在偏差，在羽流中心處觀察到的模擬速度低于實(shí)驗(yàn)測量值。這是由于氣泡羽流接近水面時(shí)，羽流的流動(dòng)方向發(fā)生了變化，垂直方向的流動(dòng)變?yōu)槎喾较蛄鲃?dòng)，其中徑向流動(dòng)占據(jù)主導(dǎo)地位。

圖11 不同水位高度處的羽流徑向速度分布Fig.11 Radial velocity distribution of plume at different water depths

圖12 為水深1.1 m、泄漏孔徑3 mm、泄漏壓力100 kPa、泄漏速率為0.21 m3/h 和水深1.1 m、泄漏孔徑3 mm、泄漏壓力300 kPa、泄漏速率為0.37 m3/h 工況下的羽流軸向速度分布，由圖可知，積分模型預(yù)測的結(jié)果與實(shí)驗(yàn)結(jié)果基本相同，羽流的軸向速度隨著水深的增加逐漸降低且在接近泄漏孔口的位置，軸線上羽流的速度急速衰減。這是因?yàn)樵跉怏w剛剛進(jìn)入水體環(huán)境時(shí)，在泄漏口的上方氣流受到水體較大阻力，此時(shí)氣體因克服阻力消耗掉大量動(dòng)能，因此羽流速度會(huì)迅速減小。而在之后的上升過程中，氣泡羽流把動(dòng)量傳遞給環(huán)境水體，周圍的水與分散氣泡沿著相同的路徑一起向上運(yùn)動(dòng)，氣液流在羽流邊界的剪切作用下速度逐漸減弱，變化幅度較為平緩。而CFD 模擬計(jì)算的結(jié)果大于實(shí)驗(yàn)結(jié)果，但隨著水位的升高，二者的差值變小，這是因?yàn)镃FD 模型是基于雷諾時(shí)均湍流方法，忽略了氣體在水中的溶解以及氣泡顆粒之間的黏性作用力。對(duì)比圖12（a）、(b)，可以看出初始泄漏速率越大，羽流軸向上升速度越大。這是因?yàn)樾孤┧俾试酱螅孤怏w的初始動(dòng)能也就越大，相應(yīng)會(huì)產(chǎn)生較大的沖量，水壓對(duì)氣體的阻礙作用相對(duì)變?nèi)?，初始段?dòng)量射流對(duì)應(yīng)的長度將會(huì)增大，羽流的上升速度相應(yīng)變大。

3.3 羽流半徑變化規(guī)律

當(dāng)水下輸氣管道發(fā)生泄漏時(shí)，氣體在管內(nèi)壓力作用下射入水體之中，氣泡羽流在上升的過程中會(huì)形成一個(gè)倒立的圓錐體結(jié)構(gòu)。圖13 是泄漏速率分別為0.21 m3/h 和0.37 m3/h 工況下的羽流半徑分布，由圖可知，積分模型預(yù)測的結(jié)果與實(shí)驗(yàn)結(jié)果基本相同。隨著氣泡羽流紊動(dòng)的發(fā)展，周圍靜止的水體被卷吸進(jìn)氣泡羽流中一起向上運(yùn)動(dòng)，流量沿程增大，羽流半徑逐漸向兩側(cè)發(fā)展并且隨著水深近似呈線性增長。但是CFD 模擬計(jì)算的結(jié)果稍微大于實(shí)驗(yàn)結(jié)果，這是因?yàn)樵贑FD 模型中羽流邊緣的氣泡顆粒分布較為分散，無法準(zhǔn)確地定義羽流半徑，因此在模型中的羽流半徑是測量羽流中心到空隙率為2%的氣泡顆粒之間的距離。

3.4 羽流泉涌高度變化規(guī)律

圖12 不同泄漏速率下中線處的羽流軸向速度分布Fig.12 Axial velocity distribution at the midline at different leak rates

當(dāng)氣泡羽流接觸到水面時(shí)，仍具有一定的速度，上升的羽流會(huì)繼續(xù)帶動(dòng)自由水面向上運(yùn)動(dòng)，導(dǎo)致溢出點(diǎn)附近區(qū)域海水向上凸起形成泉涌，水面上升的高度即為泉涌高度，如圖14 所示。溢出水面的氣體在不同條件下的泄漏行為會(huì)對(duì)水體環(huán)境和水面上的作業(yè)和人員安危造成不同程度的影響，如被引燃甚至可能造成嚴(yán)重的火災(zāi)或爆炸事故，如圖15 所示，因此對(duì)羽流泉涌高度變化規(guī)律的研究對(duì)于水下輸氣管道泄漏事件的應(yīng)急決策具有重要價(jià)值。

圖16 為水深1.1 m 時(shí)積分模型預(yù)測的不同泄漏速率下水面處羽流泉涌高度的徑向分布，由圖可知，泄漏速率為0.21、0.37和0.84 m3/h 時(shí)積分模型預(yù)測的中線處泉涌高度hf分別為0.0508、0.0748 和0.1305 m。Friedl 等[34]在2000 年提出了最大泉涌高度hpmax以及平均初始泉涌高度與積分模型預(yù)測的泉涌高度hf之間的關(guān)系式，如式（32）所示，由此可以得出積分模型預(yù)測的最大泉涌高度hpmax以及平均初始泉涌高度的值，如表3所示。

圖13 不同泄漏速率下羽流半徑分布Fig.13 Plume radius distribution at different leakage rates

圖14 氣體在水面形成的泉涌現(xiàn)象Fig.14 Fountain phenomenon of gas formation on water surface

圖15 水面溢散燃燒現(xiàn)象Fig.15 Submarine gas pipeline leakage accident

圖16 不同泄漏速率下水面處羽流徑向泉涌高度分布Fig.16 Distribution of plume radial fountain height at different leakage rates

表3 積分模型計(jì)算的泉涌高度Table 3 Fountain height calculated by integral model

圖17為實(shí)驗(yàn)與兩種模型預(yù)測的泉涌高度對(duì)比。由圖可知，在泄漏速率較小時(shí)，CFD 模擬預(yù)測的結(jié)果略低于實(shí)驗(yàn)值，但泄漏速率較大時(shí)二者差值變大，平均相對(duì)誤差為8.9%。這是因?yàn)楫?dāng)泄漏速率較大時(shí)在實(shí)驗(yàn)中會(huì)觀察到氣泡羽流上升到水面時(shí)出現(xiàn)水花濺射現(xiàn)象，這會(huì)使得測量的泉涌高度值偏大，而在CFD 模擬中并未觀察到這種現(xiàn)象。從圖中還可以看出，積分模型預(yù)測的最大泉涌高度和初始泉涌高度與實(shí)驗(yàn)差值較大，平均相對(duì)誤差為17.2%。這是因?yàn)榉e分模型中經(jīng)驗(yàn)常數(shù)β對(duì)羽流的泉涌高度影響很大，通常對(duì)于羽流的自由上升損耗過程經(jīng)驗(yàn)常數(shù)β 等于0.5，但是Friedl 等[34]將實(shí)驗(yàn)與數(shù)值模擬的數(shù)據(jù)進(jìn)行擬合后認(rèn)為β 為0.39 預(yù)測的結(jié)果最好，因此在模型中選取的經(jīng)驗(yàn)常數(shù)β 為0.39。同時(shí)由于積分模型沒有考慮在水面處水的回流以及自由水面對(duì)氣泡羽流的阻礙作用，這也會(huì)使得積分模型預(yù)測的結(jié)果大于實(shí)驗(yàn)值。

4 結(jié) 論

通過建立水下輸氣管道泄漏三維CFD 模型和積分?jǐn)?shù)學(xué)模型來研究水下輸氣管道泄漏后的擴(kuò)散特性，并利用實(shí)驗(yàn)數(shù)據(jù)對(duì)兩種模型的準(zhǔn)確性進(jìn)行了對(duì)比驗(yàn)證，具體結(jié)論如下。

圖17 實(shí)驗(yàn)與模擬預(yù)測的泉涌高度對(duì)比Fig.17 Comparison of fountain height predicted by simulations and experiments

（1）水下輸氣管道泄漏后在水體中的擴(kuò)散大致經(jīng)歷了三個(gè)不同的階段：初始階段、充分發(fā)展階段和表面流動(dòng)階段，這三個(gè)階段對(duì)應(yīng)著三個(gè)流動(dòng)區(qū)域，并且在氣泡羽流上升過程中伴隨著卷吸現(xiàn)象、紊動(dòng)沸騰現(xiàn)象。當(dāng)羽流到達(dá)表面時(shí)，夾帶的水并不會(huì)隨著氣體進(jìn)入大氣而是從水平表面向外徑向流動(dòng)，帶動(dòng)氣泡遠(yuǎn)離羽流軸心。該過程將會(huì)由于夾帶水自身的動(dòng)量而使得水面升高，從而形成泉涌。

（2）水下輸氣管道泄漏后，羽流的徑向速度近似呈高斯分布并且隨著徑向距離的增加逐漸降低；羽流的軸向速度隨著水深的增加而降低，且在接近泄漏孔口的位置由于受到水體較大的阻力，軸線上羽流的速度急速衰減；隨著氣泡羽流紊動(dòng)的發(fā)展，周圍靜止的水體被卷吸進(jìn)氣泡羽流中一起向上運(yùn)動(dòng)，羽流半徑逐漸向兩側(cè)發(fā)展并且隨著水深近似呈線性增長。

（3）耦合的VOF和DPM模型能夠用于水下氣體泄漏擴(kuò)散運(yùn)動(dòng)規(guī)律的描述，其對(duì)氣泡羽流在水體中泄漏擴(kuò)散運(yùn)動(dòng)規(guī)律的預(yù)測與實(shí)驗(yàn)基本相符，但由于模型中未考慮氣體在水中的溶解以及氣泡顆粒之間的黏性作用力，所以CFD 模擬計(jì)算的羽流軸向速度結(jié)果要略大于實(shí)驗(yàn)結(jié)果；相比于CFD 模型，由于積分模型忽略了表層水的回流以及自由水面對(duì)氣泡羽流的阻礙作用，模型計(jì)算的工作量大大減少因而計(jì)算速度也更快，但計(jì)算精度比CFD 模型稍差。同時(shí)模型中的經(jīng)驗(yàn)系數(shù)如夾帶系數(shù)和動(dòng)量放大系數(shù)等對(duì)模型計(jì)算結(jié)果影響較大，因此仍需要進(jìn)行更多深水條件下高泄漏流速的實(shí)驗(yàn)以確定準(zhǔn)確的經(jīng)驗(yàn)系數(shù)來進(jìn)一步提高積分模型的計(jì)算精度。

符 號(hào) 說 明

g——重力加速度，m/s2

Hv——水深，m

hf——泉涌高度，m

n——多變指數(shù)

Pa——大氣壓力，Pa

Qg——?dú)怏w體積流量，m3/s

r——距離羽流中心的水平距離，m

v——羽流速度，m/s

vs——滑移速度，m/s

α——夾帶系數(shù)

γ——?jiǎng)恿糠糯笙禂?shù)

λ——寬度比

下角標(biāo)

c——羽流軸線處的值

f——泉涌