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

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

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

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

(1)

ρg+F.

(2)

(3)

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

將泄漏氣體作為離散相處理,認(rèn)為其是分布在連續(xù)流場中離散的氣泡粒子。DPM是一種基于歐拉-拉格朗日方法的數(shù)值模型,采用DPM模型對海水中氣泡粒子的運(yùn)動(dòng)軌道進(jìn)行追蹤,假設(shè)作為離散相的氣泡粒子體積比率很低,分散粒子之間的相互作用以及粒子體積比率對連續(xù)相的影響均可忽略,在實(shí)際應(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)系下的氣泡粒子作用力的微分方程來求解離散相氣泡粒子的運(yùn)移軌跡,為每一個(gè)氣泡粒子施加一個(gè)平衡力,使得離散相氣泡粒子在擴(kuò)散過程中,其運(yùn)動(dòng)慣性與所受其他外力達(dá)到平衡,氣泡粒子受力平衡方程[2,13]為

(4)

(5)

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

(6)

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

(7)

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

E0表示為

(8)

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

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

2 水下氣體擴(kuò)散仿真模型

2.1 網(wǎng)格模型

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

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

2.2 邊界條件

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

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

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

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

3.1 數(shù)值模型驗(yàn)證

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

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

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

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

3.2 水下氣體羽流

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

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

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

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

圖7 氣體運(yùn)移速度分布Fig.7 Movement velocity distribution of underwater bubbles

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

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

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

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

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

圖11為水面氣池中心涌流高度隨時(shí)間的變化過程。氣體擴(kuò)散至水面時(shí),沸騰效應(yīng)引起的噴射水柱高度逐漸呈非線性和波動(dòng)狀增加。由圖11可知,32～ 39 s之間噴射水柱高度迅速增加至3 m,由于表層水回流波動(dòng)等因素,噴射水柱呈小幅度減小后,又逐漸增大。55 s時(shí)噴射水柱高度達(dá)到最大值4.4 m,55 s以后噴射水柱高度呈波動(dòng)狀發(fā)展。由圖6可知,氣體擴(kuò)散至水面時(shí),引起表層水運(yùn)動(dòng),形成一個(gè)以氣體溢出點(diǎn)為中心的圓形氣池。圖12為水面氣池尺寸隨時(shí)間的變化過程。由圖12可知,隨著氣體的持續(xù)泄漏,氣池半徑逐漸增大,75 s時(shí)氣池半徑增大至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)建立的三維水下氣體擴(kuò)散數(shù)值模型,實(shí)現(xiàn)了對水下氣體擴(kuò)散行為和水面涌流效應(yīng)的預(yù)測和定量評估,仿真結(jié)果與小尺度實(shí)驗(yàn)數(shù)據(jù)對比論證了數(shù)值模型的可行性。

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

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