基于歐拉-拉格朗日方法的水下氣體泄漏擴散行為研究

李新宏, 陳國明, 朱紅衛(wèi), 暢元江

(中國石油大學(xué)(華東)海洋油氣裝備與安全技術(shù)研究中心,山東青島 266580)

1 歐拉-拉格朗日數(shù)值建模方法

歐拉-歐拉和歐拉-拉格朗日方法可用于求解多相流動過程。歐拉-歐拉方法認為不同相是相互貫穿的介質(zhì),一相占有的體積不能被另一相占有,通過體積率來對計算域內(nèi)不同相的分布情況進行描述。體積率是關(guān)于時間和空間的連續(xù)函數(shù),不同相的體積率之和為1。VOF模型是一種在固定歐拉網(wǎng)格下的表面追蹤方法,可以實現(xiàn)對多種不相容流體之間交界面的追蹤。采用VOF模型對水和空氣兩相自由表面進行追蹤,其連續(xù)性方程、動量方程和能量方程[13]分別為

(1)

ρg+F.

(2)

(3)

式中,aq為第q相體積分數(shù);ρq為第q相密度,kg/m3;ρ為混合相密度,kg/m3;vq為第q相速度,m/s;Sh為源項;μ為湍流黏度和分子混合黏度總和,N·s/m2;g為自由落體加速度,m/s2;F為外部作用力,N;keff為熱導(dǎo)率系數(shù);T為溫度,K。

將泄漏氣體作為離散相處理,認為其是分布在連續(xù)流場中離散的氣泡粒子。DPM是一種基于歐拉-拉格朗日方法的數(shù)值模型,采用DPM模型對海水中氣泡粒子的運動軌道進行追蹤,假設(shè)作為離散相的氣泡粒子體積比率很低,分散粒子之間的相互作用以及粒子體積比率對連續(xù)相的影響均可忽略,在實際應(yīng)用中要求離散相體積比率要小于10%～12%。

The old woman put the gingerbread man in the oven.Then she went to sleep in the chair.

DPM模型通過積分拉氏坐標(biāo)系下的氣泡粒子作用力的微分方程來求解離散相氣泡粒子的運移軌跡,為每一個氣泡粒子施加一個平衡力,使得離散相氣泡粒子在擴散過程中,其運動慣性與所受其他外力達到平衡,氣泡粒子受力平衡方程[2,13]為

(4)

(5)

Re=ρdp|up-u|/μ.

(6)

式中,u為流體相對速度,m/s;up為氣泡粒子速度,m/s;FD為拖曳力,N;μ為流體動力黏度,Pa·s;ρ為連續(xù)相密度,kg/m3;ρp為氣泡粒子密度,kg/m3;dp為氣泡粒子直徑,m;Re為相對雷諾數(shù);CD為拖曳力系數(shù),氣泡受到拖曳力的大小取決于氣泡的形狀。計算拖曳力系數(shù)[14]為

(7)

式中,E0為表征氣泡特征形狀的無量綱Eotvos數(shù)。

E0表示為

(8)

式中,ρl和ρg分別為水和氣體密度,kg/m3;σl為水的黏度,N·s/m2。

此外,泄漏氣體的擴散過程屬于復(fù)雜的非穩(wěn)態(tài)湍流運動,氣體從泄漏口噴射而出具有較高的雷諾數(shù),選用Realizableκ-ε模型對氣體擴散過程中的湍流特性進行描述,并使質(zhì)量守恒、動量守恒和能量守恒的基本控制方程封閉[12]。

2 水下氣體擴散仿真模型

2.1 網(wǎng)格模型

水下氣體泄漏以后,氣體從孔口噴射進入水中,在初始動能和浮力的作用下向水面方向擴散,整個泄漏擴散過程如圖1所示。假定泄漏口所處水深為80 m,水面上部空氣域高度為20 m,參考國際油氣生產(chǎn)者協(xié)會(IOGP)發(fā)布的風(fēng)險評估數(shù)據(jù)導(dǎo)則[12],取60 mm的泄漏孔徑,建立三維幾何模型,采用分塊映射網(wǎng)格劃分方法,對泄漏口以上中心區(qū)域的網(wǎng)格進行加密處理,以適應(yīng)流場變化和保證求解精度。整個計算空間采用六面體網(wǎng)格,通過上述方法得到整個計算域的結(jié)構(gòu)體網(wǎng)格模型如圖2所示,網(wǎng)格總數(shù)為642 468個。

圖2 網(wǎng)格模型Fig.2 Mesh model

2.2 邊界條件

計算域頂部即空氣域上部邊界流體變量梯度為零,采用對稱邊界。計算域側(cè)面和底部采用無滑移邊界。泄漏口的初值條件在DPM模型中完成。氣泡粒子從面射流源釋放進入海流場中,氣體密度變化服從理想狀態(tài)方程,據(jù)Deepspill大尺度水下油氣泄漏實驗[15],氣體在水下擴散過程中氣泡尺寸分布如圖3所示。通過UDF函數(shù)實現(xiàn)對氣泡粒子密度與尺寸變化的控制。

選用非穩(wěn)態(tài)壓力基求解器,VOF模型啟用隱式體力公式,部分平衡壓力梯度和動量方程中的體積力,提高重力場中泄漏擴散模型的穩(wěn)定性。PISO算法主要適用于非穩(wěn)態(tài)問題求解,壓力速度耦合采用PISO算法[12]。耦合求解過程中,采用隨機軌道模型計入氣泡粒子在湍流流場中擴散的附加拖曳力。

圖3 氣泡尺寸分布Fig.3 Bubbles size distribution

3 水下氣體擴散模擬結(jié)果

3.1 數(shù)值模型驗證

Engebretsen等[16]在挪威Statoil研究中心進行了小尺度的水下氣體泄漏擴散實驗,所用水池長9 m,寬6 m,高7 m。實驗過程中使用空氣作為介質(zhì),采用壓力調(diào)節(jié)器維持恒定的氣體泄漏速率,分別對泄漏速率為0.083、0.17和0.75 m3/s三種條件下的氣體擴散行為進行實驗研究。建立與實驗同尺度的仿真模型,采用實驗介質(zhì)和泄漏速率,重復(fù)實驗條件下的氣體擴散過程,從而實現(xiàn)對數(shù)值模型的驗證。圖4為數(shù)值模擬與實驗結(jié)果對比。從圖4可以看出,3種泄漏速率條件下的氣體上浮時間模擬值略高于實驗值,但其差值較小。中心線速度的模擬值大于實驗值,且隨著距離泄漏源的高度增加,模擬值與實驗值的差值減小,但模擬值與實驗值的總體變化基本一致。因此,采用VOF和DPM耦合方法開展水下氣體擴散行為研究是可行的。數(shù)值模型采用雷諾時均湍流模擬方法,其假設(shè)流場是完全發(fā)展的湍流,忽略分子間的黏性作用力,導(dǎo)致射流和羽流階段的中心線速度的預(yù)測值高于實驗值。此外,由于數(shù)值仿真是理想化的環(huán)境條件,難以和實驗做到完全相符,模擬環(huán)境中海水對氣體的阻力小于實驗環(huán)境,從而導(dǎo)致模擬值與實驗值存在一定偏差。

圖4 仿真與實驗結(jié)果對比Fig.4 Comparison of simulation results with experiment

圖5為0.170 m3/s泄漏速率條件下氣體水下運動軌跡和速度分布。據(jù)圖5可知,水下氣體泄漏后,氣體向水面運動過程中,由于環(huán)境壓力降低,其體積逐漸膨脹,擴散至海面形成一個倒立錐形的羽流結(jié)構(gòu)。氣體擴散過程中引起周圍海水運動,計算域內(nèi)速度分布形狀與氣體擴散軌跡較為一致,速度隨著距離泄漏源高度的增加而減小,由羽流內(nèi)層到外層逐漸減小。氣體運動至水面時,帶動表層水體運動,從中心迅速向周圍擴展,在水面形成一個遠大于水下羽流直徑的水體波動區(qū)域。

圖5 0.170 m3/s泄漏速率條件下氣體擴散模擬結(jié)果Fig.5 Simulation result of underwater gas dispersion under release rate of 0.170 m3/s

3.2 水下氣體羽流

泄漏口距離水面的距離為80 m,氣體泄漏速率為50 kg/s,以此作為初值條件研究實例場景條件下的水下氣體擴散行為。與大氣環(huán)境相比,水中的壓力和溫度均隨水深變化。因此,氣體在水中運移擴散是一個復(fù)雜的過程。氣體在水環(huán)境運移過程中,其密度變化服從理想氣體狀態(tài)方程,上升過程中密度逐漸減小。圖6為水下氣體羽流的發(fā)展過程。

圖6 水下氣體羽流發(fā)展過程Fig.6 Development of underwater gas plume

由圖6可知,泄漏發(fā)生時,由于管道內(nèi)的高壓作用,泄漏氣體以噴射狀涌入水中,快速向上部空間發(fā)展。羽流上升過程中,環(huán)境壓力逐漸降低,體積逐漸膨脹增大。當(dāng)氣體運動至水面時,氣體與表層水相互作用,以溢出點為中心向四周拓展,形成倒立的錐形羽流結(jié)構(gòu)。

圖7為水下氣體羽流垂向運動速度分布。由圖7可知,泄漏初期,由于受到水的阻力速度急劇降低。泄漏后期,氣體在水的浮力作用下,以近似恒定的速度,繼續(xù)向水面擴散。圖8為水下氣體羽流幾何尺寸變化過程。由圖8可知,氣體羽流垂向尺寸基本呈線性增長,泄漏時間為32 s時,氣體擴散至水面。由于自由膨脹作用,氣體向水面擴散的同時,逐漸向水平方向運動。據(jù)圖8可知,由于受到相同的自由膨脹作用,氣體羽流橫向尺寸變化過程基本一致。

圖7 氣體運移速度分布Fig.7 Movement velocity distribution of underwater bubbles

圖8 水下氣體羽流幾何尺寸Fig.8 Geometric sizes of underwater gas plume

3.3 水面涌流效應(yīng)

氣體運動至水面時,仍具有一定的速度,能夠帶動表層水運動,在水面溢出點附近形成沸騰效應(yīng),導(dǎo)致溢出點附近區(qū)域海水向上凸起,形成噴射水柱或涌流。溢出水面的氣體如被引燃甚至可能造成嚴重的火災(zāi)或爆炸事故,該類事件已被工程實際所證實。2011年12月中海油珠海橫琴天然氣處理終端附近的海底輸氣管道發(fā)生泄漏[17],在該海域作業(yè)的漁民發(fā)現(xiàn)泄漏氣體在海面形成的噴射水柱,后期運營方將泄漏天然氣在海面做燃燒處理,形成10 m高的火焰高度,如圖9所示。

圖10為氣體擴散至水面時空氣-水相界面變化與相界面區(qū)域速度矢量分布。由圖10可知,氣體擴散至水面時引起溢出點處相界面凸起,即圖9(a)中的涌流現(xiàn)象。氣體運動至海面時,在溢出點區(qū)域形成沸騰效應(yīng),引起表層水回流,在溢出點區(qū)域兩側(cè)形成漩渦流動,并向溢出點兩側(cè)傳遞,該區(qū)域的速度矢量分布如圖10(b)所示。

圖9 珠海海底輸氣管線泄漏事故Fig.9 Leakage accident of subsea gas pipeline in Zhuhai

圖11為水面氣池中心涌流高度隨時間的變化過程。氣體擴散至水面時,沸騰效應(yīng)引起的噴射水柱高度逐漸呈非線性和波動狀增加。由圖11可知,32～ 39 s之間噴射水柱高度迅速增加至3 m,由于表層水回流波動等因素,噴射水柱呈小幅度減小后,又逐漸增大。55 s時噴射水柱高度達到最大值4.4 m,55 s以后噴射水柱高度呈波動狀發(fā)展。由圖6可知,氣體擴散至水面時,引起表層水運動,形成一個以氣體溢出點為中心的圓形氣池。圖12為水面氣池尺寸隨時間的變化過程。由圖12可知,隨著氣體的持續(xù)泄漏,氣池半徑逐漸增大,75 s時氣池半徑增大至42 m,75 s以后氣池半徑趨于穩(wěn)定。

圖11 氣池中心涌流高度變化Fig.11 Variation of fountain height

圖12 氣體溢出區(qū)域尺寸變化Fig.12 Variation of radius of surface gas pool

4 結(jié) 論

(1)建立的三維水下氣體擴散數(shù)值模型,實現(xiàn)了對水下氣體擴散行為和水面涌流效應(yīng)的預(yù)測和定量評估,仿真結(jié)果與小尺度實驗數(shù)據(jù)對比論證了數(shù)值模型的可行性。

(2) 水下泄漏發(fā)生時,氣體以噴射狀涌入水中,在初始動能和浮力作用下向水面運動。由于環(huán)境壓力降低,氣體密度逐漸減小,體積膨脹增大,形成倒立錐形的羽流結(jié)構(gòu)。氣體擴散至海面時,帶動表層水運動,引起羽流頂部兩側(cè)水回流,產(chǎn)生漩渦流動,在水面形成涌流效應(yīng)。氣體從水面溢出,引起表層水波動,以氣體溢出點為中心向四周發(fā)展,在水面形成一個圓形的氣池。水面涌流高度隨時間呈非線性增長,達到峰值后呈小幅波動狀變化,水面氣池半徑隨時間逐漸增大后穩(wěn)定。

(3)離散相與連續(xù)相耦合的方法能夠有效用于水下氣體泄漏擴散精細分析,確定水下氣體運移軌跡、上浮時間、水面涌流高度和氣池尺寸等重要參數(shù),對水下氣體泄漏事件風(fēng)險評估和應(yīng)急決策具有較好的參考價值,從而提高風(fēng)險防控與事故處置的科學(xué)性與合理性。