， ， 建星， ，

(1. 天津大學 水利工程仿真與安全國家重點實驗室， 天津 300072；2. 中國艦船研究設計中心， 湖北 武漢 430064)

0 引 言

近年來中國快速上升的綜合國力與能源供應之間的矛盾十分突出，海洋的油氣資源儲量巨大，成為解決能源問題的重要手段。海底油氣輸送系統(tǒng)面臨極端環(huán)境和復雜地質條件，一旦發(fā)生失效，維修與更換困難，不但會導致正常管道的運輸功能停滯，使經濟收益受到大幅損失，而且將污染水質，破壞海洋生態(tài)環(huán)境，造成災難性后果。若水下輸氣管道破裂的位置處于海洋平臺附近，溢出的氣體有可能意外引燃，進而對上部平臺造成更為嚴重的事故影響。

對于水下氣體泄漏，國內外學者的研究包括試驗研究、理論研究和數(shù)值模擬等3部分。在試驗研究方面，BULSON較早用試驗的方法對氣泡羽流及其工程應用作了系統(tǒng)研究；MCDOUGALL[1]和MILGRAM[2]建立了氣泡羽流模型，并通過大量的試驗研究得出氣泡擴散半徑的經驗公式；20世紀80年代末期，吳鳳林等[3-4]專門對圓形氣泡羽流形成區(qū)運動形態(tài)作了細致分析。在理論研究方面，BUSCAGLIA等[5]分析了質量轉移效應對大尺度氣泡羽流的數(shù)值計算模型的影響；馬霞[6]通過編制程序，預測靜止環(huán)境中的圓形氣泡羽流流場并與試驗數(shù)據(jù)對比；徐婷婷[7]對氣泡羽流流型試驗進行研究，改變曝氣量和水深條件來研究氣泡羽流的流動形態(tài)。在數(shù)值模擬方面，DHOTRE等[8]采用計算流體動力學(Computational Fluid Dynamics， CFD)仿真對大尺寸氣泡羽流進行研究；尹群等[9]運用軟件NEPTUNE對不同破損情況下天然氣泄漏擴散情況進行模擬；范開峰等[10]采用VOF模型研究海底管道天然氣泄漏在水中上升過程的體積分數(shù)分布特性。

國內對水下油氣泄漏過程的研究主要是對國外研究的吸收和再現(xiàn)[11]，對氣泡羽流運動特性數(shù)學模型的完善，然而國內實際海況下氣體泄漏問題的分析很少，特別是在實際工程應用中對海底管道氣體泄漏的一些重要預測參數(shù)的確定還存在一定困難。本文基于一組水下氣體泄漏釋放試驗，在MATLAB中進行數(shù)學理論模型計算，結合渤海實際工況，將FLUENT軟件中的多相流模型(VOF模型)和離散相模型耦合，建立水下氣體泄漏有限元數(shù)值模擬模型。對潮流場進行二次擬合的方式加入海洋環(huán)境載荷，對比數(shù)學模型理論計算結果和有限元數(shù)值模擬結果，得到天然氣泄漏的偏移距離和擴散半徑等關鍵參數(shù)，探究水下氣體泄漏運動的規(guī)律，為海洋工程實踐的安全仿真系統(tǒng)提供研究基礎。

1 水下氣體泄漏運動理論

對水下油氣泄漏的情況來說，氣體占據(jù)控制體內部空間且不會形成水合物，可用理想氣體的狀態(tài)方程來表述氣體在整個過程中的變化，本文研究中有如下假定。

管道氣體泄漏假定：(1)氣體在管道內做一維穩(wěn)態(tài)流動；(2)泄漏口的流動為等熵流動；(3)管道內氣體服從理想氣體運動規(guī)律；(4)氣體在管道內的流動為絕熱流動。此假定用于天然氣泄漏速率計算。

圖1 水下油氣泄漏理想狀態(tài)

氣泡在水中運動基本方程假定：(1)初始半徑較小的氣泡在短距離上升過程中形狀始終保持為球形；(2)氣泡內氣體溫度在運動過程中保持不變。

當環(huán)境流體的速度不很大時，如圖1所示，為淺水情況下的典型油氣泄漏情況，此時氣體始終沒有從主控制體側面逸出，主控制體的運動軌跡幾乎不會發(fā)生彎折或彎折曲率的絕對值很小。由于氣泡與控制體內部的液體之間存在滑移速度ωb，從而使氣體從控制體上部逸出，并最終上浮至水平面滑移半徑為βb。βb的取值范圍為0.65～0.80(可取0.70)，ωb的取值范圍為0.25～0.35 m/s。

根據(jù)控制單元體內的具體情況，利用有限差分對控制方程進行離散，得到理論計算結果。

1.1 氣泡體積分數(shù)計算

控制單元體內部氣泡的體積分數(shù)ε可以通過下式[12]進行計算：

(1)

式中：ρl為控制單元體內部液體的密度；ρ為氣液混合態(tài)下的混合體的密度；ρb為氣體的密度。這些單位都是kg/m3。

1.2 質量控制方程

如果考慮氣體溶解的影響，假定每個截面的質量流率保持為定值，其質量控制方程[12]如下：

(2)

關于氣體溶解的計算，在通常情況下，氣泡的溶解率可以通過如下方程得到：

(3)

式中：n為氣泡所含氣體的體積摩爾數(shù),mol；K為質量遷移系數(shù)，m/s；Co為溶解氣體的摩爾濃度，mol/m3；Cs為氣體的飽和濃度，即氣體的溶解度。其中，質量遷移系數(shù)和氣體的溶解度為主要控制參數(shù)。

1.3 動量控制方程

在動量方程中，控制單元體是均勻的，忽略其表征參數(shù)(如溫度、濃度等)的分布特征，同時將氣泡的滑移速度考慮在內。此外，根據(jù)一組對氣體泄漏的試驗，將文獻[12]的系數(shù)β修正為β2，能更準確地描述氣體運動的動量控制。不考慮拖曳力影響，具體方程形式如下：

·u]=uaρaQε

(4)

(5)

(6)

式(2)～式(4)中：β為氣體截面面積與控制單元體截面面積之比；ml為控制單元體內部液體的質量，m1=ρlπb2(1-β2ε)h,kg；mb為控制單元體內氣體的質量,kg；wb為滑移速度；u、v、w分別為控制單元體在3個方向上的速度；ua、va、wa為環(huán)境流體在3個方向上的速度分量；g為重力加速度；ρb為控制單元體總密度；ρa為氣體在周圍環(huán)境流體中的密度；ρl為控制單元體內部液體的密度；Qe為環(huán)境流體夾帶率；Qε為氣泡的體積流量；h為控制單元體的高度；b為控制單元體的半徑。

2 有限元數(shù)值模型基本理論

FLUENT軟件流體力學計算有速度快、穩(wěn)定性強和精度高等優(yōu)點。VOF模型作為其提供的一種多相流模型，在整個計算域內對不互溶的流體求解同一個動量方程組并追蹤每種流體的體積分數(shù)來模擬多相流?；谒職怏w泄漏運動中氣泡球形顆粒構成的第二相分布在海水這一連續(xù)相中的特點，結合使用DPM模型可以更好地計算這些氣泡顆粒的運動以及由此引起的熱量/質量傳遞。

2.1 VOF模型

作為一種多相流模型，VOF模型能描繪連續(xù)相行為，通過求解單獨的動量方程和處理穿過區(qū)域每一流體的溶劑比來模擬2種或3種不能混合的液體。VOF模型基本控制方程[13]如下：

質量守恒方程：

ρi)+·(aiρivi)=0

(7)

動量守恒方程：

(8)

組分質量守恒方程：

(9)

式(7)～式(9)中：ρ為環(huán)境流體密度；V為速度矢量；Yj為組分j的質量分數(shù)；Jj為組分j的擴散通量；Rj為系統(tǒng)內部單位時間內單位體積通過化學反應消耗或生成該種組分的凈生成率；Sj為通過其他方式(如異相反應、相變)所生成該種組分的凈生產率以及用戶自己定義的其他質量源項。

VOF模型典型的應用包括流體噴射、流體中氣泡運動、流體在大壩壩口的流動、氣液界面的穩(wěn)態(tài)和瞬態(tài)處理等。

2.2 DPM模型

離散顆粒模型(Discrete Particle Model， DPM)能計算水下油氣泄漏中油滴和氣泡顆粒的運動軌跡，充分考慮相間耦合作用。DPM模型的優(yōu)點在于能準確追蹤模擬計算區(qū)域內的離散粒子，粒子受拖曳力、重力和其他外力的聯(lián)合作用達到受力平衡，具體如下：

(10)

3 淺水海底管道氣體泄漏實例

3.1 目標油田環(huán)境條件

選取渤海某目標油田氣體泄漏模擬分析案例。由海圖和附近海域的水深數(shù)據(jù)可推算目標油田海域水深(海圖水深基準面起算)為17 m，油田工程設計水位數(shù)據(jù)見表1。

表1 目標油田設計水位 m

圖2 浪玫瑰圖(年)

考慮目標油田海域的風、浪、流因素，具體的風浪分布圖和環(huán)境參數(shù)見表2～表4和圖2。

表2 目標油氣田海域氣溫 ℃

表3 目標油氣田海域相對濕度 %

表4 目標油氣田海域海水溫度 ℃

3.2 環(huán)境數(shù)據(jù)擬合

依據(jù)油田環(huán)境條件分別對各個方向上的潮流流速進行二次擬合，得到其平均流速和最大流速的函數(shù)擬合公式，主流向NE方向潮流擬合圖如圖3～圖4所示。

圖3 NE方向上的潮流流速擬合 圖4 WSW方向上的潮流流速擬合

3.3 管道泄漏口參數(shù)

管道所處水深為17.0 m，壓強P=0.17 MPa，管道起點處的壓力為6 MPa,天然氣溫度為323.15 K，管道摩擦系數(shù)取0.5，管內天然氣甲烷含量為94.72%，摩爾質量為17.097 kg/ kmol，假定距管道起點等效長度Le=1 km處管道發(fā)生泄漏，氣體常數(shù)取8.314 Pa·m3/(mol·K)，管道直徑為1.219 m，泄漏孔徑取圓孔。計算不同泄漏孔徑下氣體泄漏速度，具體情況如表5所示。天然氣管道的泄漏事故規(guī)模一般是按泄漏口徑的大小劃分的。根據(jù)工程經驗，小規(guī)模(小口徑)泄漏事故發(fā)生的可能性較大，大規(guī)模(大口徑)泄漏事故發(fā)生的可能性較小。即天然氣管道泄漏一般為孔口泄漏，泄漏過程可看作是可壓縮氣體的孔口出流。為方便建模，這里以圓孔口泄漏為基礎來推導計算公式。在該實例中定義泄漏口徑100 mm以下為小孔泄漏，泄漏口徑達100 mm及以上為大孔泄漏。

表 5 不同破損口徑下氣體泄漏速率

3.4 有限元數(shù)值模型

在FLUENT中建立的模型是一個長度L=80 m，高度H=28 m的矩形，其中水面在從底邊起17 m高處，水面上的介質是空氣，具體如圖5所示。A、B、D 3個邊界定義為速度入口類型，F(xiàn)與C兩條邊界定義為wall類型，E定義為出口中的溢出類型。其中水和天然氣的湍流黏度比設為17，而D邊上來流速度由UDF加載3.3中定義的各方向下潮流的擬合速度。

圖5 淺海管道氣體泄漏有限元模型

網格劃分為正方形，單位長度為0.1 m。由于海底管道為細長型構件，因此采用雙精度模式。湍流設置選用標準雙k-e方程，歐拉相數(shù)為3(其中水為主相，空氣、甲烷為第二相)。設置模型底邊為管道，甲烷氣體泄漏口在底邊中點F處，呈射流狀?？衫肍LUENT自帶的Injections功能，設置破損孔徑的大小、氣體泄漏速度和質量速度等參數(shù)。

為使流場內的初始速度為穩(wěn)定海流速度，使用自定義初始化程序進行初始化DEFINE_INIT，初始化過程采用相應的穩(wěn)定流場參數(shù)，計算時間為小孔泄漏(20 mm)氣體開始泄漏后25 s左右，大孔泄漏(100 mm、300 mm、500 mm)氣體開始泄漏后10～15 s。

4 計算結果對比

4.1 數(shù)學理論計算結果

在數(shù)學模型中加載目標油田具體的管道泄漏口參數(shù)，將潮流條件代入到模型初始條件中進行計算，選取主風向NE方向這一組氣體泄漏數(shù)學模型計算結果，如圖6和圖7所示。

圖6 泄漏口徑d=20 mm時，NE流向氣體泄漏軌跡 圖7 泄漏口徑d=100 mm時，NE流向氣體泄漏軌跡

4.2 有限元模擬結果及對比

4.2.1 小孔徑泄漏

隨機選取小孔徑泄漏的兩組FLUENT模擬圖如圖8所示。

圖8 兩個示例流向氣體泄漏密度分布

通過圖8，結合氣體泄漏動態(tài)過程，可以看出：FLUENT數(shù)值模擬與實際天然氣泄漏后的擴散運動過程類似，兩種模型的耦合一方面完成了氣液兩歐拉相在自由表面上的相互作用，使泄漏氣體離散成更加貼切的氣泡粒子，另一方面能真實地模擬泄漏的氣體羽流形態(tài)以及運動軌跡。

研究水下氣體泄漏運動的規(guī)律，關鍵參數(shù)是泄漏的偏移距離和擴散半徑，綜合統(tǒng)計16個潮流方向的兩種模型計算結果見表6和圖9。

表6 各潮流方向的兩種模型計算結果對比誤差率 %

圖9 小口徑泄漏的羽流對比圖

由表5～圖9可知：在淺水海底管道氣體的小口徑泄漏中，數(shù)學理論計算模型和有限元數(shù)值模型模擬氣體泄漏運動的偏移距離和擴散半徑這兩種關鍵參數(shù)的誤差均控制在10%以內，結果擬合度較高。

4.2.2 大孔徑泄漏

隨機選取大孔徑泄漏的兩組FLUENT模擬圖如圖10所示。

圖10 大口徑泄漏的羽流對比圖

變換大孔徑泄漏的口徑，按照3.2部分泄漏口參數(shù)，分別選取口徑為100 mm、300 mm和500 mm，綜合兩種模型的計算結果如圖11和圖12所示。

圖11 不同泄漏口徑時不同潮流流向的羽流半徑

圖12 不同泄漏口徑時不同潮流流向的偏移距離

由圖11和圖12可知：在淺水海底管道氣體的大口徑泄漏中，數(shù)學理論計算模型和有限元數(shù)值模型模擬氣體泄漏運動得到的偏移距離和擴散半徑這兩種關鍵參數(shù)的結果擬合度高，兩種模型相互驗證，結果合理正確。

4.2.3 結果對比

在不同管道泄漏口徑下，對比不同潮流流向氣體泄漏的氣泡羽流半徑、流方向偏移距離，如圖13和圖14所示。由圖13和圖14可知：在相同潮流作用、不同泄漏率(泄漏孔徑一樣)下，泄漏氣體到達海面的擴散半徑、流方向漂移距離是不同的。氣泡羽流到達水面的擴散半徑隨泄漏率的增大而增大，流方向漂移距離隨泄漏率的增大而減小。原因在于泄漏率越大，相對的運動速度越大，到達海面的運動時間則越小。

圖13 不同泄漏口徑的擴散羽流半徑對比 圖14 不同口徑泄漏的氣體運動偏移距離對比

5 結 論

運用數(shù)學理論計算模型與有限元數(shù)值模型相結合來研究淺水海底管道氣體泄漏問題，可以很好地考慮海洋環(huán)境載荷因素的影響，相互對比驗證，更加準確客觀。其中，將FLUENT軟件中的VOF模型和DPM模型耦合一方面能夠準確地仿真水下氣體泄漏的整體羽流運動形態(tài)，另一方面能夠模擬氣體泄漏過程中氣泡離散顆粒的運動軌跡。

對于淺水氣體泄漏得到以下結論：(1)海洋的潮流流向、泄漏口徑(即泄漏率)大小影響泄漏氣體的運動，在海面得到不同的擴散半徑、流方向偏移距離。(2)泄漏氣體到達水面的擴散半徑隨著泄漏率的增大而增大，流方向偏移距離隨著泄漏率的增大而減小。

上述結論對于進行管道氣體泄漏點的定位、預報海底管道氣體泄漏情況以及實施應對處理措施都有著很大的參考作用。