楊大慎，陳 強，謝 成，陳曉華，姜紅濤，趙彥琳

(1.國家石油天然氣管網(wǎng)集團有限公司華南分公司，廣東 廣州 510620；2.中國石油大學(xué)(北京)機械與儲運工程學(xué)院，北京 102249)

近年來，隨著社會經(jīng)濟的發(fā)展，對能源的需求不斷增長，管道運輸在成品油遠距離運輸中具有不可替代的作用。運輸管道的油品泄漏會對企業(yè)經(jīng)濟和生態(tài)環(huán)境帶來巨大的損失和破壞[1]。2009—2019年海底輸油管道事故統(tǒng)計信息顯示，我國共發(fā)生100余起泄漏事件，僅2001年發(fā)生的東海油田管道泄漏事故造成的直接經(jīng)濟損失就超過2 000萬元。因此，油氣運輸管道的安全運行越來越受到人們的重視，而開展泄漏油品的軌跡和擴散范圍的研究可為事故的緊急處理提供及時有效的參考。

針對油品泄漏擴散特征，國內(nèi)外學(xué)者進行了一些相關(guān)研究，例如BEMPORD[2]考慮水流的影響，對分層流工況下圓孔的浮射流軌跡進行了研究，建立了穩(wěn)定分層流動液體中水下浮力圓形射流的數(shù)值模擬模型，用于分析穩(wěn)定或中性分層情況下不同流動條件時潛水浮力射流的行為。JOHVANSEN[3]提出了一種基于拉格朗日概念的海底井噴模型，研究了分層水柱中上升氣泡在水中的溶解，得到模擬結(jié)果與野外試驗觀測結(jié)果較吻合的結(jié)論。UDITHA等[4]研究了氣泡尺寸分布在深水天然氣或石油運輸中的作用機制，通過數(shù)值模擬和實驗對比證明了氣泡破裂和聚結(jié)在泄漏模型中的重要性。KANDEL等[5]通過實驗合成水溶性陰離子分散劑，研究了不同分散劑對抑制成品油污染擴散的影響。高龍驤等[6]改進了經(jīng)典溢油模型，建立了內(nèi)河溢油動態(tài)模型，實現(xiàn)了對溢油過程的仿真與預(yù)測。劉冰等[7]針對不同的河流水體和河流環(huán)境溢流的特點，結(jié)合基本控制方程與水動力數(shù)學(xué)模型，改進Fay模型、簡化Navy漂移模型，建立了新的河流溢油環(huán)境影響評價模型。?；勰鹊萚8]結(jié)合水環(huán)境的脆弱性和對人體的危害，建立了河流水環(huán)污染風險模糊綜合評價模型。

綜上可以看出，輸油管道成品油泄漏研究大多數(shù)是從實驗研究和數(shù)值模擬兩個方面開展，而且主要針對海水管道、井下泄漏等油品污染的擴散情況進行研究，較少涉及河流跨越段輸油管道泄漏后成品油在河流流域擴散速度與擴散范圍的研究。因此，筆者基于北盤江流域成品油管道跨越段，在不同工況下對不同泄漏位置、不同泄漏時間和不同泄漏方式下成品油在江中的擴散速度與擴散范圍進行數(shù)值模擬，獲得不同影響因素作用下油品擴散的軌跡、速度和污染范圍，以期為管道運輸安全、工程搶險救災(zāi)、水環(huán)境保護提供技術(shù)參考。

1 研究原理及方法

1.1 研究方法概述

研究區(qū)選取某一河段流域，通過Getdata軟件獲得河岸坐標點，再導(dǎo)入矢量作圖軟件SolidWorks中構(gòu)建幾何模型，并設(shè)置油品污染物的泄漏位置；通過CFD前處理軟件對構(gòu)建的幾何結(jié)構(gòu)模型進行網(wǎng)格劃分，生成結(jié)構(gòu)性網(wǎng)格；根據(jù)實際問題進行不同程度的歸一理想化，調(diào)整各物性參數(shù)、邊界條件，建立污染擴散模型，對計算結(jié)果進行對比分析，整理出不同條件下該流域的污染擴散規(guī)律。具體研究思路如圖1所示。

圖1 研究思路

1.2 網(wǎng)格劃分

采用ICEM 軟件來生成網(wǎng)格。考慮到結(jié)構(gòu)化網(wǎng)格的優(yōu)點，如容易實現(xiàn)區(qū)域的邊界擬合，適于流體和表面應(yīng)力集中等方面的計算，網(wǎng)格生成速度快、質(zhì)量好，數(shù)據(jù)結(jié)構(gòu)簡單等，故將計算域模型畫成結(jié)構(gòu)化網(wǎng)格，并在計算過程中不斷加密網(wǎng)格，找出網(wǎng)格無關(guān)化解。

1.3 離散方法

采用ANSYS R19.2進行數(shù)值模擬，選擇有限體積法求解，采用QUICK離散格式，壓力和速度耦合求解采用SIMPLE(semi-implicit method for pressure linked equations)算法[9]?；赟IMPLE算法，以壓力為基本求解變量，通過密度延遲修正方法將SIMPLE算法擴展到可壓縮流動的計算[10]。使用非交錯網(wǎng)格使程序更加簡潔通用，運用動量插值法消除由非交錯網(wǎng)格所產(chǎn)生的非物理性壓力振蕩現(xiàn)象[11]。使用標準κ-ε湍流模型對粘性流場進行模擬，對于固體壁面附近的區(qū)域采用壁面函數(shù)法處理[12]。SIMPLE算法是一種壓力修正法，采用預(yù)測-修正計算方法，加快了每個時間步長的收斂速度，效率較高。動量空間離散采用二階迎風格式，體積分數(shù)的空間離散應(yīng)用Geo-Reconstruct格式，湍動能和比耗散率的空間離散均采用一階迎風格式[13]。

1.4 湍流模型控制方程

按照假設(shè)，流體為不可壓縮的牛頓流體，其基本控制方程采用連續(xù)性方程與不可壓縮流體的雷諾N-S方程:

▽·u=0

(1)

(2)

式中：u為流體速度；ρ為流體密度；p為流體壓力；v為流體運動黏度；fi為作用在流體上的外力。

1.5 VOF多相流模型

VOF多相流模型基于流體體積函數(shù)進行計算[14]，即利用許多單個網(wǎng)格中一種流體體積與網(wǎng)格體積的比值f來進行界面追蹤[15]。定義流體體積分數(shù)fq為網(wǎng)格內(nèi)第q相流體所占體積與該網(wǎng)格總體積之比[16]。若fq=1，則表示網(wǎng)格內(nèi)均為第q相流體，不含其他流體；若fq=0，則表示網(wǎng)格內(nèi)不存在第q相流體；若0

(3)

∑fq=1

(4)

在二維計算域中，假設(shè)fo、fw分別為成品油和水在一個網(wǎng)格中所占的體積分數(shù)，則fo、fw的計算公式為：

fo=Vo/V

(5)

fw=Vw/V

(6)

其中，V、Vo、Vw分別為網(wǎng)格體積、網(wǎng)格中成品油的體積和網(wǎng)格中水的體積。

結(jié)合式(3)～式(6)，可推導(dǎo)出計算域內(nèi)成品油和水兩相混合流動過程的連續(xù)性方程：

(7)

(8)

此外，在VOF方法中，多相流的物性參數(shù)φ由計算區(qū)域內(nèi)各相的物性參數(shù)和體積分數(shù)共同決定，即：

φ=∑φqfq

(9)

1.6 模型假設(shè)

在數(shù)值模擬中，建立模型需要針對現(xiàn)實情況做出相應(yīng)的簡化和假設(shè)，忽略次要因素，考慮主要影響，以提高計算效率。

(1)所建立的模型只包含兩種流體，油品和水，且二者在流動過程中完全不相溶，二者均不可壓縮，密度、粘度等物理性質(zhì)均保持不變；泄漏油品以油膜的形式隨河流表面向下游擴散遷移。

(2)在管道發(fā)生完全斷裂時，上下游的閥門關(guān)閉后，管道中的殘余油品經(jīng)0.5 h就完全流入河流中，在山坡上無殘留油品，河流兩側(cè)泄漏口泄漏時完全沿著山坡流入河道，油品泄漏時風速保持恒定。

(3)在管道發(fā)生小孔泄漏時，從上下游閥門關(guān)閉到封堵成功，管道內(nèi)壓力和小孔泄漏流量基本保持不變。

(4)油品和水不發(fā)生任何化學(xué)反應(yīng)，油水混合物揮發(fā)量不計。

2 北盤江計算模型

2.1 北盤江概述

北盤江為典型的山區(qū)性河流，洪枯水季節(jié)分明。每年6月至次年5月為枯水期。洪水期徑流主要由降雨形成，水量約占全年的80%，水位暴漲；枯水期徑流主要由地下水補給，流量約占全年的20%，水位相對平穩(wěn)。北盤江多年平均流量為309 m3/s，多年平均徑流量為98.18億m3。由于北盤江跨越段所處地區(qū)，山區(qū)、河谷交錯，成品油管道一旦發(fā)生泄漏，水流彌散作用將造成泄漏油品迅速向下游遷移，短時間內(nèi)即對下游造成大面積污染，給當?shù)氐纳鷳B(tài)、漁業(yè)、旅游、農(nóng)業(yè)等造成嚴重影響。因此，對北盤江跨越段成品油管道的泄漏擴散特征進行研究是非常必要的。

研究區(qū)選取放馬坪閥室至小水井閥室的北盤江跨越管段，采用單跨懸索跨越方案，管道通過支架固定到橋上，管道總長203 m，共設(shè)置支架21個。第一個支架距離橋端3 m，支架之間距離大約為10 m?，F(xiàn)場調(diào)研發(fā)現(xiàn)，距離橋端33 m處的支架因雨水進入較多，腐蝕相對嚴重，故將其設(shè)置為典型泄漏計算點。

2.2 流域計算模型的建立

通過CAD建模，選取有限的區(qū)域段建立相應(yīng)的計算域。以北盤江跨越段為起點，以直線距離20.2 km的位置為終點。東西方向(x方向)距離跨度約為9.4 km，南北方向(y方向)距離跨度約為16.0 km。對北盤江跨越段下游河兩岸進行數(shù)據(jù)分析處理，進一步細化所需河道信息，順序提取河岸兩側(cè)輪廓，獲得河道坐標數(shù)據(jù)。將河道坐標數(shù)據(jù)點導(dǎo)入矢量建模軟件中，并將北盤江跨越段在河流中的起點設(shè)為坐標原點(0,0)。通過繪圖軟件將地理位置坐標進行離散化，獲得北盤江的下游河道計算模型的坐標參數(shù)。北盤江河岸兩側(cè)輪廓及河道坐標數(shù)據(jù)點分布如圖2所示。

圖2 北盤江河岸兩側(cè)輪廓及河道坐標數(shù)據(jù)點分布

以圖2的邊界數(shù)據(jù)為基礎(chǔ)，進行計算域的建模。將河道兩岸的散點坐標導(dǎo)入icem軟件，得到河流兩岸邊界的軌跡線，在起始位置和出口位置形成閉合區(qū)域，再對封閉曲線進行內(nèi)部填充，從而得到二維的計算域模型圖，計算域y方向長度為12 014 m，x方向長度為9 634 m。入口邊界條件設(shè)為velocity-inlet；出口邊界條件設(shè)為pressure-outlet；兩側(cè)河岸的邊界條件設(shè)為wall，無滑移邊界條件；成品油進入河道區(qū)域的邊界長度設(shè)為10 m，邊界條件設(shè)為mass-flow-inlet。

在模擬計算中，北盤江跨越段的河道信息及相關(guān)流域的水文資料均被考慮進模型中，由于北盤江是季節(jié)性河流，河水的流速和流量均受季節(jié)影響，故筆者將北盤江分為枯水期、平水期和豐水期3種，以探討不同時期泄漏溢油的污染狀況和擴散范圍。北盤江跨越段不同時期水流流量和流速的平均值如表1所示。

表1 北盤江跨越段不同時期水流流量與流速的平均值

同時，將泄漏工況按照破裂嚴重性分為兩大類：管道發(fā)生完全斷裂和局部小孔泄漏。管道發(fā)生完全斷裂可能是因為山區(qū)意外地質(zhì)災(zāi)害、人為破壞或者年久失修管道銹蝕，不同位置發(fā)生的完全斷裂所造成的溢油擴散范圍可能不同，故將泄漏點分為管道南側(cè)、管道中間和管道北側(cè)3處，以探討不同泄漏點位置對溢油擴散情況的影響。當管道破損不大，僅發(fā)生小孔泄漏時，成品油泄漏量相對于管道完全斷裂來說比較小，因此僅考慮成品油在跨越段管道中間泄漏的情況，成品油從小孔泄漏后直接落入北盤江中。

2.3 泄漏量計算

北盤江跨越段的管道長度L為203 m，管道直徑D為406.4 mm，壁厚為12.7 mm，設(shè)計壓力為10 MPa。管道材質(zhì)為X60鋼，制管方式為LSAW，輸送介質(zhì)為97#、93#汽油和0#柴油等成品油?？缭蕉螛蛄壕嚯x水面的垂直距離約為18 m。泄漏檢測系統(tǒng)可以在成品油泄漏發(fā)生20 s內(nèi)做出正確的泄漏報警，閥門由全開至全關(guān)需20 s。假設(shè)該段成品油最大輸量Q年為1 000萬t/年，成品油密度ρ為850 kg/m3。

2.3.1 管道完全斷裂的泄漏量

當管道發(fā)生完全斷裂事故時，從事故開始到閥門關(guān)閉這段時間溢油量Q1的估算公式為：

(20+20)≈14.92 m3

(10)

式中：T為一年的周期長度；t為泄漏時間。

緊急關(guān)閉閥門后，管道殘余量一般為上下游閥室之間所剩余的成品油量。估算時按照泄漏的最大情況來分析，所以留在管道內(nèi)的成品油泄漏量Q2為：

(11)

則最大成品油泄漏量Q為：

Q=Q1+Q2=41.24 m3

(12)

2.3.2 管道小孔泄漏的泄漏量

當管道發(fā)生小孔泄漏，設(shè)定小孔的孔徑d為10%D(按照整理的有關(guān)數(shù)據(jù)，管道發(fā)生泄漏的最大小孔泄漏尺寸為10%D)。小孔的通流長度即管道的壁面厚度l為12.7 mm，則l/d=12.7/40.6=0.31<0.5，符合流體力學(xué)中關(guān)于薄壁小孔的定義要求。

通過流體力學(xué)中的伯努利方程，可得小孔泄漏的連續(xù)性方程為：

(13)

式中：ρ為溢油的密度；v1、v2分別為泄漏小孔內(nèi)側(cè)截面和外側(cè)截面處溢油的速度；p1、p2分別為小孔內(nèi)側(cè)截面和外側(cè)截面處的壓力，p1取管道壓力，即10 MPa，p2取環(huán)境壓力，即101 325 Pa；ζ為動能修正系數(shù)；Cv為小孔流速系數(shù)。

則管道小孔泄漏的成品油流量q為：

(15)

式中：A為泄漏小孔外側(cè)截面的面積；Cc為小孔收縮系數(shù)，即最小通流截面的面積與泄漏小孔外側(cè)截面面積的比值；Cd為小孔流量系數(shù)，由實驗測得。完全收縮時，流體在泄漏孔處的雷諾數(shù)較大，Cd的取值范圍為0.61～0.62；不完全收縮時，Cd的取值范圍為0.7～0.8。本次計算為完全收縮，故Cd的取值為0.62。經(jīng)計算，小孔泄漏的成品油流量q為0.11 m3/s。

當管道突發(fā)小孔泄漏事故時，會經(jīng)歷泄漏警報和閥門關(guān)閉過程。雖然閥門關(guān)閉后泄漏的成品油流量相對比較小，但是泄漏持續(xù)時間較長，根據(jù)經(jīng)驗推算緊急堵漏的時間為2 h。則從泄漏發(fā)生至閥門關(guān)閉再到緊急封堵泄漏這段時間的溢油量Q′為：

Q′=q×(20+20+7 200)=796.4 m3

(16)

3 模擬結(jié)果分析

參考北盤江跨域段管道的數(shù)據(jù)與水文資料確定模擬計算中所需參數(shù)，根據(jù)所提供的參數(shù)條件設(shè)置計算條件，如北盤江江水的密度為998.2 kg/m3，動力黏度為1.003 mPa·s；成品油密度設(shè)為850.0 kg/m3，動力黏度為0.08 Pa·s。

3.1 管道完全斷裂擴散結(jié)果分析

工況①～工況③分別為北盤江枯水期管道南側(cè)、管道中間、管道北側(cè)發(fā)生完全斷裂時成品油的泄漏工況，枯水期3個工況不同時刻泄漏成品油擴散的數(shù)值模擬結(jié)果分別如圖3～圖5所示。計算結(jié)果顯示，泄漏0.5 h后，成品油完全流進入北盤江，泄漏成品油剛開始集中在河流中央，成塊狀聚集。泄漏發(fā)生2 h后，油品逐漸被沖散，塊狀主體主要集中在跨越段下游4.1～4.4 km；泄漏4 h后，成品油擴散范圍更大，油水混合物還附著在岸邊，主要集中在跨越段下游7.3～8.2 km。

圖3 工況①不同時刻的成品油擴散情況

圖4 工況②不同時刻的成品油擴散情況

圖5 工況③不同時刻的成品油擴散情況

工況④～工況⑥分別為北盤江平水期管道南側(cè)、管道中間、管道北側(cè)發(fā)生完全斷裂時成品油的泄漏工況，平水期3個工況不同時刻泄漏成品油擴散的數(shù)值模擬結(jié)果分別如圖6～圖8所示。在平水期3個工況，同樣是泄漏0.5 h后成品油全部流入北盤江，泄漏成品油最初集中在跨越段管道附近河道內(nèi)，呈大范圍塊狀聚集；泄漏2 h后，成品油主要集中在跨越段下游9.6～11.8 km處的河道中；泄漏4 h后，成品油大多已經(jīng)成為油水混合細帶狀，進一步污染到下游14.6～15.7 km處。

圖6 工況④不同時刻的成品油擴散情況

圖7 工況⑤不同時刻的成品油擴散情況

圖8 工況⑥不同時刻的成品油擴散情況

工況⑦～工況⑨分別為北盤江豐水期管道南側(cè)、管道中間、管道北側(cè)發(fā)生完全斷裂時成品油的泄漏工況，豐水期3個工況不同時刻泄漏成品油擴散的數(shù)值模擬結(jié)果如圖9～圖11所示。由于豐水期水流速度較快，泄露4 h后工況⑦～工況⑨的油品擴散范圍已超出計算域，因此選取了t=3 h的擴散情況進行對比分析。豐水期3個工況泄漏的成品油在北盤江中擴散得較為迅速，按照泄漏0.5 h后成品油全部流入北盤江，此時成品油富集于管道附近大部分河岸。泄漏2 h后，油膜集中在跨越段下游12.2～13.8 km的河道范圍內(nèi)；泄漏3 h后，成品油擴散成小碎塊狀，且有部分滯留于河流彎頭處，在下游15.6～17.8 km處有較大塊油膜存在。

圖9 工況⑦不同時刻的成品油擴散情況

圖10 工況⑧不同時刻的成品油擴散情況

圖11 工況⑨不同時刻的成品油擴散情況

3.2 管道小孔泄漏擴散結(jié)果分析

圖12 工況⑩不同時刻的成品油擴散情況

圖13 工況不同時刻的成品油擴散情況

4 結(jié)論

筆者基于有限元體積法和κ-ω模型，采用追蹤多相流界面的VOF方法，建立了一種跨越段管道泄漏污染擴散模型，得到以下結(jié)論：

(1)模型計算結(jié)果可得出不同工況下成品油擴散范圍與泄漏時間的關(guān)系，豐水期管道完全斷裂泄漏擴散范圍最廣，平水期次之，枯水期最小。

(2)管道泄漏點位置在泄漏前期對成品油擴散影響較大，但隨著擴散時間的增加，在泄漏后期成品油擴散范圍主要由水流速度和河岸的地形決定。在成品油污染擴散前期，泄漏的成品油大部分聚集在一起，尚未被水流沖散，適合收集處理。因此，在管道泄漏點附近，宜設(shè)置圍油欄等應(yīng)急處理裝置進行油品攔截回收。

(3)河道中的彎道容易產(chǎn)生漩渦，漩渦附近會造成油品滯留，也會使經(jīng)過的大塊油膜被沖散，部分油品會附著在此處的河岸上。應(yīng)急搶險時可在彎道處安裝攔油裝置，有利于節(jié)約時間和回收成本。

(4)小孔泄漏造成的成品油擴散泄漏流量小、泄漏持續(xù)時間長，不會產(chǎn)生大塊的油膜。豐水期水流速度快，小孔泄漏的油品基本被水流沖散，污染擴散快且不利于收集。但無論是枯水期還是豐水期，泄漏前期的應(yīng)急處理都很關(guān)鍵，對降低泄漏造成的流域污染損失起到重要作用。

(5)VOF方法追蹤互不相容且物性參數(shù)相差較大的兩種或多種流體的多相流流動界面時效果較好，且可根據(jù)離散方法改進計算精度。