長(zhǎng)江溢油事故中數(shù)值模型的研究與應(yīng)用

韓龍喜張防修

(1.揚(yáng)州市職業(yè)大學(xué) 資源與環(huán)境工程學(xué)院，江蘇 揚(yáng)州 225000； 2.河海大學(xué) 環(huán)境科學(xué)學(xué)院，南京 210000)

1 研究背景

長(zhǎng)江航道是我國(guó)重要的水運(yùn)物資通道，近年來(lái)由于航運(yùn)繁忙，已發(fā)生多起因事故引發(fā)的溢油污染事故[1]。事故發(fā)生后，如何科學(xué)預(yù)測(cè)溢油輸移擴(kuò)散特征，對(duì)于提高應(yīng)急處置水平具有重要意義。溢油數(shù)值模擬是當(dāng)前該領(lǐng)域的研究熱點(diǎn)，現(xiàn)有的研究更多集中于海洋、河口區(qū)域[2-6]，而針對(duì)狹長(zhǎng)型內(nèi)陸河道的研究相對(duì)有限。Sayed等[7]、Goeury[8]在河道溢油事故處理中構(gòu)建的數(shù)值模型簡(jiǎn)便實(shí)用。Toz等[9]、Gautama等[10]建立的河道溢油數(shù)值模型綜合考慮了水流、隨機(jī)擴(kuò)展以及溢油風(fēng)化等因素的影響，并將遙感技術(shù)應(yīng)用于模型的驗(yàn)證中，取得了較好的應(yīng)用效果。在國(guó)內(nèi)，劉棟等[11]基于機(jī)理性實(shí)驗(yàn)分析了潮汐作用下河道中溢油的擴(kuò)展、漂移特性。趙琰鑫等[12]在長(zhǎng)江武漢段采用傳統(tǒng)經(jīng)驗(yàn)?zāi)Ｐ秃陀土Ｗ幽Ｐ拖嘟Y(jié)合的方法分析了溢油遷移特征。

溢油在擴(kuò)展過(guò)程中不接觸水陸邊界是傳統(tǒng)經(jīng)驗(yàn)?zāi)Ｐ蛻?yīng)用于溢油自身擴(kuò)展分析的前提條件。這種假設(shè)在事發(fā)水域較寬敞且溢油以中心排放時(shí)較易實(shí)現(xiàn)，但是如果事發(fā)水域?yàn)閮?nèi)陸河道或溢油以岸邊排放形式進(jìn)入水體時(shí)則很難實(shí)現(xiàn)。針對(duì)狹窄水域中發(fā)生的溢油事故，韓龍喜等[13]提出的基于油粒子模型分析溢油漂移擴(kuò)散和擴(kuò)展過(guò)程的研究方法能較準(zhǔn)確地反映溢油擴(kuò)散過(guò)程，但是該方法單純通過(guò)水平擴(kuò)散方式模擬油膜自身擴(kuò)展，導(dǎo)致油粒子擴(kuò)散速度響應(yīng)過(guò)慢，且未能有效反映溢油量對(duì)擴(kuò)散面積影響。

本文采用水深平均二維潮流模型模擬評(píng)價(jià)區(qū)域水流流場(chǎng)，采用油粒子模型分析溢油輸移擴(kuò)展。為了能準(zhǔn)確反映溢油自身擴(kuò)展的變化過(guò)程，采用隨機(jī)運(yùn)動(dòng)模擬了溢油的慣性擴(kuò)展、黏性擴(kuò)展和表面張力擴(kuò)展過(guò)程，使油粒子自身擴(kuò)展更符合Fay提出的三階段擴(kuò)展理論[14]。數(shù)值模型同時(shí)考慮了溢油風(fēng)化及觸岸吸附等過(guò)程。將數(shù)值模型應(yīng)用于假設(shè)工況的溢油事故分析中，取得了理想預(yù)測(cè)效果。該方法可為內(nèi)陸河道突發(fā)性溢油事故預(yù)報(bào)和應(yīng)急處理提供技術(shù)支持。

2 水動(dòng)力學(xué)模型

根據(jù)長(zhǎng)江河道特點(diǎn)，本文采用水深平均二維潮流模型[15]模擬研究區(qū)水流場(chǎng)，水動(dòng)力數(shù)據(jù)獲得準(zhǔn)確性驗(yàn)證后，可以作為溢油水動(dòng)力場(chǎng)的驅(qū)動(dòng)條件。水動(dòng)力模型通過(guò)連續(xù)方程和動(dòng)量守恒方程來(lái)表征流場(chǎng)和潮位的變化特征。

連續(xù)方程為

(1)

動(dòng)量守恒方程為

(2)

式中：Z，H分別表示水位和水深(m)；u，v分別為x，y向的流速分量(m/s)；νt為紊動(dòng)粘性系數(shù)(m2/s)；g為重力加速度；c為謝才系數(shù)；f為柯氏系數(shù)，f=2ωsinφ，其中ω為地球自轉(zhuǎn)角速，φ為計(jì)算地水域的地理緯度。

上述二維潮流模型基本方程中含有非線性混合算子，采用剖開(kāi)算子法[16]進(jìn)行離散求解。

3 油粒子預(yù)測(cè)模型

油粒子模型[17-18]最早由Johansen和Andunson(1982)提出，該模型假設(shè)溢油為大量離散的油粒子聚合體，每個(gè)油粒子代表一定油量。模型通過(guò)綜合分析油粒子在水流、風(fēng)向以及自身隨機(jī)游走作用下的運(yùn)動(dòng)輸移，分析溢油隨時(shí)間遷移及其空間分布特征。油粒子模型將溢油在Δt時(shí)間內(nèi)的運(yùn)動(dòng)分成對(duì)流、風(fēng)導(dǎo)漂移和隨機(jī)游走3個(gè)過(guò)程。單個(gè)油粒子的運(yùn)動(dòng)方程可表示為

Xn+1=Xn+ΔXC+ΔXW+ΔXD。

(3)

式中：Xn+1，Xn分別表示單個(gè)油粒子在(n+1)Δt，nΔt時(shí)刻的空間位置列向量；ΔXC，ΔXW，ΔXD分別表示該油粒子受表層水流對(duì)流運(yùn)動(dòng)影響、風(fēng)應(yīng)力影響以及水體紊動(dòng)擴(kuò)散影響而產(chǎn)生的空間位置變化的列向量。

3.1 溢油隨流運(yùn)動(dòng)模擬

在Δt時(shí)段內(nèi)，受水流對(duì)流作用產(chǎn)生的油粒子空間位移可以用確定性方法模擬，即

(4)

式中Un，Un+1分別表示油粒子在nΔt，(n+1)Δt時(shí)刻的水流速。

3.2 溢油風(fēng)導(dǎo)(應(yīng)力)漂移模擬

在Δt時(shí)段內(nèi)，受表層風(fēng)應(yīng)力直接作用，油粒子產(chǎn)生的漂移距離ΔXW可以表示為

ΔXW=αDW10Δt。

(5)

式中：α為風(fēng)漂移因子，取值范圍為0.03～0.04；D為考慮風(fēng)向偏轉(zhuǎn)角的轉(zhuǎn)換矩陣；W10為水體表面10 m高處的風(fēng)速向量。

3.3 溢油隨機(jī)游走模擬

受隨機(jī)游走影響，油粒子云團(tuán)的尺度和形狀隨時(shí)間不斷變化。本文基于隨機(jī)步長(zhǎng)法[19]計(jì)算油粒子在水平方向上產(chǎn)生的隨機(jī)運(yùn)動(dòng)距離ΔXD，即

(6)

式中：Rn表示[-1，1]之間的均勻分布隨機(jī)數(shù)；θ′表示位于[0，π]之間的均勻分布的隨機(jī)角度；De是由溢油自身擴(kuò)散引起的擴(kuò)散系數(shù)；DT是油粒子作布朗運(yùn)動(dòng)的擴(kuò)散系數(shù)。

DT的計(jì)算式可表示為

(7)

式中：nb是曼寧系數(shù)；Vc為水流速；h為水深。

采用油粒子模型模擬溢油機(jī)械擴(kuò)展時(shí)，存在溢油運(yùn)動(dòng)響應(yīng)慢，且未能反映溢油量對(duì)溢油擴(kuò)展影響的不足。鑒于Fay提出的溢油擴(kuò)展三階段理論能較真實(shí)地反映溢油機(jī)械擴(kuò)展過(guò)程[18]。本文采用Fay公式推導(dǎo)溢油自身擴(kuò)散系數(shù)De，如式(8)所示。該方法的實(shí)質(zhì)是通過(guò)溢油隨機(jī)運(yùn)動(dòng)反映溢油慣性擴(kuò)展、黏性擴(kuò)展和表面張力擴(kuò)展過(guò)程。

式中：常數(shù)K2i=1.14，K2v=1.45，K2t=2.30；Vp是油粒子的體積(m3)；σn為張力系數(shù)；ρw，γw分別為水的密度和運(yùn)動(dòng)黏性系數(shù)。

3.4 溢油風(fēng)化過(guò)程

溢油風(fēng)化包括蒸發(fā)、乳化及溶解等過(guò)程。風(fēng)化過(guò)程對(duì)溢油組成產(chǎn)生影響，但對(duì)溢油水平位置不發(fā)生影響[20]。

3.4.1 蒸 發(fā)

在溢油事故初期，溢油質(zhì)量傳輸主要通過(guò)蒸發(fā)完成。溢油蒸發(fā)率與溢油本身物化特征有關(guān)，通常通過(guò)溢油蒸發(fā)率表示，即

(9)

3.4.2 乳 化

溢油乳化表現(xiàn)為油包水的過(guò)程。溢油乳化進(jìn)程受水域條件影響，但最終乳化量?jī)H取決于溢油中乳化劑的含量。本文采用含水率來(lái)表征溢油乳化程度，即

(10)

3.4.3 溶 解

溶解是指油粒子在外界能量的擾動(dòng)下，均勻進(jìn)入水體的過(guò)程，溶解量和溶解速度主要受溢油自身物化性質(zhì)、風(fēng)速和水況條件的影響。通常由溢油溶解率反映溢油溶解程度，即

(11)

3.5 河岸吸附

溢油沿水流進(jìn)入長(zhǎng)江岸邊區(qū)域時(shí)，由于水流速往往與水邊線存在一定夾角。溢油粒子極易觸岸，根據(jù)Shen等[21]的研究結(jié)論，觸岸油粒子可能會(huì)出現(xiàn)被岸邊完全吸收、部分吸收、完全反射3種可能。本文作簡(jiǎn)化分析，只考慮完全吸收和完全反射2種情況。設(shè)定一個(gè)黏附概率，某油粒子觸岸時(shí)會(huì)獲得一個(gè)介于0～1之間的隨機(jī)數(shù)，當(dāng)該隨機(jī)數(shù)小于設(shè)定黏附概率時(shí)，油粒子被完全吸收，該油粒子坐標(biāo)被修正為計(jì)算網(wǎng)格中距離靠岸點(diǎn)最近的節(jié)點(diǎn)坐標(biāo)。

3.6 油膜厚度計(jì)算

油膜厚度變化是溢油事故應(yīng)急處理中一個(gè)重要分析指標(biāo)。在掌握油粒子空間分布規(guī)律后，可通過(guò)分析一定水體內(nèi)油粒子數(shù)量、質(zhì)量、體積等參數(shù)計(jì)算得到油膜厚度分布。

假設(shè)水中油粒子個(gè)體大小均等，A表示水面面積，N表示水面油膜中包含的油粒子個(gè)數(shù)，m表示受風(fēng)化作用后單個(gè)油粒子質(zhì)量。按液體厚度計(jì)算公式，油膜在t時(shí)刻的厚度ht可表示為

(12)

式中ρ為溢油密度。

4 溢油事故數(shù)值模擬

選擇鎮(zhèn)江揚(yáng)州段河道中儀征碼頭上游12 km至下游13 km范圍，作為溢油事故模擬的研究區(qū)域(119°13′E—119°08′E，32°22′N—32°13′N)。該江段包含了儀征市水廠的重要取水口，周邊分布有碼頭、潤(rùn)揚(yáng)大橋、江心洲、邊灘等，航道條件較為復(fù)雜(如圖1)。近年來(lái)長(zhǎng)江航運(yùn)業(yè)日益繁忙，過(guò)往船舶在水上運(yùn)輸或碼頭作業(yè)過(guò)程中一旦發(fā)生溢油事故，將對(duì)該地區(qū)工農(nóng)業(yè)生產(chǎn)和群眾飲水安全產(chǎn)生嚴(yán)重影響。

圖1 模擬江段位置示意圖Fig.1 Position of the simulated reach of the Yangtze River

溢油事故產(chǎn)生的影響受水文、氣象等多條件綜合作用。本文根據(jù)水環(huán)境敏感目標(biāo)位置、水文水動(dòng)力條件以及污染源位置的代表性確定6種計(jì)算方案，具體信息見(jiàn)表1。這里按照考慮最不利影響原則，選擇方案3(溢油時(shí)刻為大潮漲潮時(shí)，風(fēng)向取ESE方向，風(fēng)速取平均風(fēng)速2.6 m/s，化學(xué)品本底濃度取0 mg/L)具體分析說(shuō)明。

表1 溢油事故風(fēng)險(xiǎn)預(yù)測(cè)方案Table 1 Schedules of risk prediction for oil spill accidents

4.1 模擬事故工況

根據(jù)以往研究成果[22]確定長(zhǎng)江鎮(zhèn)揚(yáng)段河道的糙率系數(shù)，長(zhǎng)江主槽糙率一般取值0.018～0.022，河道灘地糙率一般取值0.024～0.028。分散系數(shù)選用Ex=αxhu確定，αx取為4.0，αy取為0.5。基于GIS平臺(tái)，采用矩形網(wǎng)格劃分水動(dòng)力模擬計(jì)算區(qū)域，平面共布置360(縱向)×300(橫向)個(gè)節(jié)點(diǎn)(網(wǎng)格)，網(wǎng)格邊長(zhǎng)50 m。依據(jù)1∶10 000水下地形等值線圖，計(jì)算提取各節(jié)點(diǎn)河底高程?？紤]到待評(píng)價(jià)江段與大通站間支流入流量相對(duì)較小，這里以大通站最小月平均流量作為一維水流模擬的上邊界條件；用同期的下游潮位站潮位過(guò)程作為下邊界條件[23]。經(jīng)一維水動(dòng)力學(xué)數(shù)學(xué)模型模擬后得到評(píng)價(jià)區(qū)域二維水動(dòng)力學(xué)模擬的上、下游邊界水文要素變化過(guò)程，并以此作為設(shè)計(jì)潮流量、潮位邊界條件，模擬設(shè)計(jì)潮流過(guò)程的水動(dòng)力特征。

模擬事故情景為某油船在儀征擬建碼頭處發(fā)生原油泄露，溢油源強(qiáng)取為140 t，泄漏時(shí)間為10 min，溢油油品為輕質(zhì)中原原油，其理化性質(zhì)如表2所示。

表2 油樣主要理化性質(zhì)Table 2 Main physical and chemical properties of oil samples

注：ρ0為油品密度；μoil為原油黏度

4.2 結(jié)果分析

采用二維非穩(wěn)態(tài)水動(dòng)力模型模擬分析研究區(qū)水動(dòng)力流場(chǎng)特征。水動(dòng)力模型通過(guò)多個(gè)測(cè)站的潮位潮流資料進(jìn)行驗(yàn)證。大潮漲急時(shí)刻流場(chǎng)見(jiàn)圖2。采用Willmott[24]的統(tǒng)計(jì)學(xué)方法對(duì)水動(dòng)力模型計(jì)算結(jié)果進(jìn)行效率評(píng)價(jià)，結(jié)果顯示效率評(píng)價(jià)值為極好。水動(dòng)力模型能較好地呈現(xiàn)研究區(qū)復(fù)雜水流特征。

圖2 大潮漲急時(shí)刻流場(chǎng)Fig.2 Flow field at spring tide

如圖3所示，根據(jù)模擬計(jì)算結(jié)果，溢油在大潮漲潮開(kāi)始時(shí)進(jìn)入長(zhǎng)江，在潮流和風(fēng)力的共同作用下，油粒子大致沿西方向沿上岸向上游漂移。儀征取水口離事發(fā)位置較近，溢油事故發(fā)生后對(duì)取水口的影響如表3所示，約20 min后油膜到達(dá)儀征水源地準(zhǔn)保護(hù)區(qū)下邊界，油膜在事故發(fā)生后150 min離開(kāi)水源地保護(hù)區(qū)，對(duì)儀征水源地保護(hù)區(qū)的持續(xù)影響時(shí)間為130 min，對(duì)取水口的持續(xù)影響時(shí)間為5 min。水源地保護(hù)區(qū)內(nèi)油膜最大厚度范圍為2.67～42.36 mm，其中取水口附近油膜最大厚度為6.77 mm。在未采取有效應(yīng)急措施的情況下，溢油事故對(duì)該水源地水質(zhì)造成直接影響。因此，應(yīng)及時(shí)啟動(dòng)水源地應(yīng)急監(jiān)測(cè)計(jì)劃，確保供水安全。

圖3 油粒子在事故發(fā)生后各時(shí)刻運(yùn)動(dòng)軌跡Fig.3 Trajectories of oil particles at differentinstances after an accident

位置影響目標(biāo)油粒子中心到達(dá)時(shí)間/min溢油持續(xù)影響時(shí)間/min折算油膜最大厚度/mm儀征取水口準(zhǔn)保護(hù)區(qū)下邊界202042.36二級(jí)保護(hù)區(qū)下邊界401014.60一級(jí)保護(hù)區(qū)下邊界50108.45儀征取水口6056.77一級(jí)保護(hù)區(qū)上邊界80156.46二級(jí)保護(hù)區(qū)上邊界130505.21準(zhǔn)保護(hù)區(qū)上邊界150202.67

5 結(jié) 論

溢油事故是長(zhǎng)江航運(yùn)過(guò)程中的易發(fā)事故類型。本文針對(duì)長(zhǎng)江河道溢油事故特點(diǎn)，采用二維潮流模型模擬評(píng)價(jià)區(qū)域水流流場(chǎng)。在溢油遷移擴(kuò)散分析中，采用改進(jìn)型油粒子模型分析溢油遷移擴(kuò)散特征，使得模型較準(zhǔn)確地反映了溢油遷移和擴(kuò)展過(guò)程。該模型同時(shí)考慮了溢油風(fēng)化及岸邊吸附對(duì)溢油遷移的影響。將建立的溢油事故分析模型應(yīng)用于鎮(zhèn)江—揚(yáng)州江段中虛擬事故的分析，結(jié)果表明在風(fēng)向取ESE方向、大潮漲潮時(shí)發(fā)生溢油污染事故，油膜大致沿西方向沿上岸向上游漂移。該模型不僅可以科學(xué)預(yù)測(cè)油膜對(duì)儀征水源地保護(hù)區(qū)和取水口的持續(xù)影響時(shí)間，以及油膜厚度變化，對(duì)于分析溢油遷移擴(kuò)散和厚度變化也具有較好的指導(dǎo)意義。

目前該主題的研究更多關(guān)注溢油在二維平面的遷移變化，在后續(xù)研究中，可以嘗試開(kāi)展溢油在縱向立面的三維分析。在研究方法上，可以嘗試將水動(dòng)力模型、水質(zhì)模型及GIS動(dòng)態(tài)可視化技術(shù)集成，以提高分析結(jié)果的呈現(xiàn)效果。