金 戈，董仁杰，蔡 睿，陳紅剛，金 生

(1.大連海事大學，遼寧 大連 116026；2.大連市水務集團水資源有限公司，遼寧 大連 116021；3.大連理工大學，遼寧 大連 116024)

海上溢油對海洋環(huán)境帶來巨大破壞。能夠準確快速地對溢油軌跡進行模擬預測對溢油事故的處理至關重要。一個簡單、準確、有效的溢油模擬系統(tǒng)可以模擬溢油的運動軌跡、行為歸宿與成分變化，從而對敏感地區(qū)溢油預防和災情發(fā)生后即時處置提供重要支持。因此對溢油模型進行研究，開發(fā)完善的溢油模擬系統(tǒng)具有重要的理論與實際意義。

隨著水環(huán)境模型與計算機技術的飛速發(fā)展，數(shù)值模擬計算軟件的適用范圍、計算能力、便利程度都在不斷提高。很多成熟的水環(huán)境數(shù)值模擬軟件隨之涌現(xiàn)出來，如丹麥水力學研究所(DHI)的MIKE系列模擬軟件、荷蘭三角洲研究院(Deltares)的Delft 3D等。同樣，溢油模擬也有許多較為成熟的軟件，如美國國家海洋與大氣研發(fā)中心(NOAA)研究開發(fā)的GNOME平臺、ASA公司研發(fā)的OILMAP平臺、挪威科技工業(yè)研究院的OSCAR、丹麥水環(huán)境研究院的MIKE/SA模塊等。但它們大都需要通過其它渠道獲得流場數(shù)據(jù)才能對溢油進行模擬，并沒有水動力、溢油集成計算的功能。國內(nèi)對于模擬軟件研發(fā)工作的成果相對較少。本文搭建了包括水動力與溢油模擬的計算程序平臺Hydroinfo，并使用該平臺對近期發(fā)生的“桑吉輪”溢油事故進行了預測。

1 模型理論

1.1 Hydroinfo水利信息系統(tǒng)

Hydroinfo水利信息系統(tǒng)是一款計算復雜水流與輸運問題的數(shù)值模擬軟件。它能真實準確地模擬環(huán)境水流的多種形態(tài)，從河流海洋的一維、二維、三維近似，到非恒定的水流、波浪、泥沙與輸運等各相關領域。其中包括三維自由水面流模擬、二維水流泥沙波浪問題、流域河網(wǎng)與管網(wǎng)、多維度耦合模塊等，本文的潮流場則是基于三維自由水面流模擬模塊進行計算。

1.2 水動力模型理論

1.2.1基本方程

如果直接對NS方程模擬則需要很高的時間與空間分辨率來獲得不同尺度的湍流信息，導致計算量過大而無法計算。雷諾平均模擬和大渦模擬是常用模擬湍流運動的方法，其中雷諾平均模擬只計算平均運動，無需計算各尺度的湍流運動就可以得到湍流平均信息，計算量小，是目前工程中最廣泛使用的模擬方法，也是本文所采用的方法。

連續(xù)性方程：

(1)

運動方程：

(2)

式中，ρ—流體密度；ν—運動粘度系數(shù)；fi—質(zhì)量力強度；ui—速度，其中i=1、2、3。

Boussinesq假設認為湍流脈動所造成的附加應力，也類似層流運動應力，可以由時均的應變率與Boussinesq湍動粘度乘積表示。通過該假設得到雷諾應力與平均速度梯度的關系為：

(3)

式中，μt—引入的湍動粘度；k—湍動能。

湍動能k表達式：

(4)

在混合長度理論中Prandtl假定湍動粘度μt與時均速度ui的梯度和混合長度lm的乘積成正比，在二維問題中可以表示為：

(5)

混合長度lm是由實驗確定或經(jīng)驗公式計算得出的值。

海洋中水平方向尺度遠大于水深方向尺度，根據(jù)這一特點可以將方程進行簡化，將水深方向應力與加速度忽略得到三維淺水Reynolds時均方程，將連續(xù)性方程沿水深方向積分，并通過Leibniz法則[7]可以得到水位方程。

(6)

式中，η—水位；zb—底高程。

1.2.2ALE任意拉格朗日——歐拉描述法

自由表面的晃動往往會給數(shù)值模擬問題的求解帶來困難，傳統(tǒng)的描述方法通常使用拉格朗日描述法或歐拉描述法。純歐拉描述與純拉格朗日描述各有優(yōu)缺點，而任意拉格朗日——歐拉描述法(ALE)則是結合兩者特點產(chǎn)生的方法，該方法中生成的計算網(wǎng)格可以任何形式在空間運動，既可以獨立于空間坐標系，也可以獨立于物體坐標系。

使用ALE描述法對動邊界處理有很大幫助，但同時會在方程中引入對流項使計算量增加，所以如何規(guī)定合適的網(wǎng)格運動方法是ALE描述中的關鍵問題。本文在自由表面法向方向使用拉格朗日描述以準確描述自由表面運動，在自由表面切線方向使用歐拉描述以防止網(wǎng)格發(fā)生扭曲。在這樣的設定下網(wǎng)格只能沿自由表面法向做上下移動，既可以很好地描述動邊界，也不會增加過多的計算量。

1.3 溢油模型理論

海上溢油的歸宿是一個包含多種物理、化學變化的復雜過程，每個過程之間相互影響，共同組成了溢油的整個行為歸宿。本文在“油粒子”法的基礎上進行改進，將溢油視為由大量油粒子組成，并對每一個油粒子行為歸宿進行模擬，所有油粒子的行為合在一起便模擬出整個溢油過程。

1.3.1擴展過程

本文是將Fay三段過程中的重力——粘性力階段公式的微分形式應用于每一個油粒子上，模擬每一個油粒子單獨的擴展過程，并不是先將整個溢油視作一個整體進行擴展過程模擬的傳統(tǒng)方法，再拆分成多個油粒子進行后續(xù)過程的處理，而是直接對每一個油粒子在每個時間步長中的擴展情況進行模擬，從而得到每一時刻每個油粒子的擴展面積，這樣可以針對持續(xù)性溢油進行模擬，同時也為后續(xù)過程模擬提供了必要的數(shù)據(jù)基礎。通過量綱分析的方式得到擴展過程的計算公式，溢油表面積隨時間變化率與溢油初始體積和溢油表面積有關，引入溢油擴展系數(shù)ks，計算表達式如下：

(7)

式中，A—油粒子表面積；t—時間；ks—溢油擴展系數(shù)；V—油粒子體積。

1.3.2遷移過程

本文從受力分析的角度出發(fā)，能夠從原理上更加精確地對溢油遷移過程進行模擬。

油與水之間的摩擦力與海水密度、海水與溢油速度差的平方有關，所以通過量綱分析引入油與水摩擦系數(shù)fw。同理油與風之間的摩擦力與空氣密度、風與溢油速度差的平方有關，引入油與風摩擦系數(shù)fa。溢油的受力與速度變化可表示為：

(8)

(9)

(10)

(11)

式中，τx、τy—溢油在x、y方向所受切應力；fw與fa—油與水摩擦系數(shù)和油與風摩擦系數(shù)；ρw與ρa—海水密度與空氣密度；vw，x、vo，x、vw，y、vo，y—水與風x方向的速度和水與風y方向的速度；ax與ay—溢油在x、y方向的加速度；m—油粒子質(zhì)量；A—油粒子表面積。

式(8)、(9)表示溢油所受水平切向力由溢油與海水的摩擦力和溢油與風的摩擦力構成，式(10)、(11)為牛頓定律，即加速度大小正比于合外力大小，與物體的慣性質(zhì)量成反比。

1.3.3蒸發(fā)過程

本文通過量綱分析方式得到蒸發(fā)過程的計算公式，溢油蒸發(fā)體積量隨時間變化率與溢油初始體積和溢油表面積有關，引入溢油蒸發(fā)系數(shù)ke，具體方程如下：

(12)

式中，A—油粒子表面積；t—時間；ke—溢油蒸發(fā)系數(shù)；V—油粒子體積。

1.3.4溶解過程

本文計算表面溶解量是通過量綱分析方式得到溶解過程的計算公式，溢油溶解量隨時間變化率與溢油初始體積和溢油表面積有關，引入溶解擴展系數(shù)kd，計算表達式如下：

(13)

式中，A—油粒子表面積；t—時間；kd—溢油溶解系數(shù)；V—油粒子體積。

溶解物溶于水中后會隨著水流運動，其運動滿足物質(zhì)運輸?shù)膶α鲾U散方程，表達式為：

(14)

式中，S—溶解物濃度；t—時間；u、v、w—x、y、z三個方向的速度；εh、εv—水平方向與垂直方向的擴散系數(shù)。

2 “桑吉輪”溢油事故模擬

2.1 事故簡介

2018年1月6日19∶50，“桑吉輪”與香港籍散貨船“長峰水晶”輪發(fā)生劇烈碰撞，碰撞地點位于長江口以東約160n mile處，“桑吉輪”在碰撞后發(fā)生火災，燃燒了整整8d后發(fā)生爆炸，隨后沉入海底，伴隨其沉入海底的還有13.6萬t凝析油，船上32名船員全部遇難?！伴L峰水晶”號船員在碰撞后棄船登艇，之后被附近的浙岱漁03187船救起，該船船長鄭磊等人目睹了整個碰撞過程。萬幸的是長峰水晶輪最終在救助拖輪的護航下靠泊舟山港卸貨。

2.2 模型參數(shù)

Hydroinfo水利信息系統(tǒng)與Google Earth有對接接口，可直接從Hydroinfo頁面打開Google Earth并下載地圖圖片與地形數(shù)據(jù)。本文使用該功能下載地圖圖片作為底圖，并獲得其地形數(shù)據(jù)作為高程。由于底圖尺寸與地形關系，還需在底圖的基礎上自行畫出邊界，圈出計算域，計算域如圖1所示，計算時長為10d，從事故發(fā)生之日起。邊界處的陸地部分可以設置為固壁，沒有潮流流通也不允許油粒子穿越，而海面部分則可通過Hydroinfo中的潮位預報工具給出。

圖1 計算域與高程二維展示

潮位預報通過對中國近海潮位數(shù)據(jù)進行傅里葉分析，從而得到該海域附近各點在各時間的潮位。海面邊界潮位過程通過邊界兩端點的潮位過程插值給出，即使用兩點潮位過程設置邊界，十分快速便捷，避免了人工輸入大量數(shù)據(jù)的繁瑣。

本文的風場數(shù)據(jù)是從美國國家海洋和大氣管理局(NOAA)地球系統(tǒng)研究實驗室(ESRL)的物理科學部(PSD)網(wǎng)站上得到，風場0°為正東方向，90°為正北方向，風場數(shù)據(jù)見表1。

表1 風場數(shù)據(jù)

由于關心溢油點附近的變化情況，所以在該區(qū)域畫了矩形加密線對其進行更密集的網(wǎng)格[12]。外邊界的網(wǎng)格尺度為20000m，即平均每20000m生成一個網(wǎng)格，而加密線內(nèi)區(qū)域的網(wǎng)格尺度被設置成5000m。最終生成16271個網(wǎng)格與8240個節(jié)點，生成的網(wǎng)格如圖2所示，中間的矩形區(qū)域為加密區(qū)域。

圖2 網(wǎng)格生成情況

參數(shù)設定包括物理參數(shù)與計算參數(shù)。其中物理參數(shù)包括：水位初值為0；海底糙率0.025；單個油粒子最大初始質(zhì)量10000kg；溢油油品密度800kg/m3；空氣密度1.293kg/m3；海水密度1025kg/m3；溢油揮發(fā)組分0.95；溢油揮發(fā)系數(shù)0.001/h；溢油溶解組分0.01；溢油溶解系數(shù)0.001/h；溢油面積擴展系數(shù)0.01/s。

控制參數(shù)包括：時間步長每步0.01h；計算步數(shù)24000步；垂向分層2層。

2.3 結果與分析

經(jīng)過計算后得到溢油模擬的計算結果，由于篇幅有限，在此只展示10d的結果，如圖3所示。中國大陸沿海的潮位變化最大，潮位最高處為高于海平面2m左右，潮位最低處為低于海平面2m左右。圖中有部分區(qū)域始終顯示為紅色，如溫州與福州沿海、臺灣北部和琉球群島，這是由于地形明顯高于海平面導致的，并不是真實水位。

圖3 10d潮流場

油粒子的運動軌跡是本次模擬最關心的內(nèi)容?？梢钥闯鲇土Ｗ哟篌w是向北方漂移并呈條帶狀，在10d溢油主要有兩條軌跡。1月16日，國家海洋局中國海警在事故現(xiàn)場進行了監(jiān)視監(jiān)測，并在事故點附近發(fā)現(xiàn)油污帶。9∶00，有長約9km，寬50～500m的油污帶在北側2km處被發(fā)現(xiàn)，10∶00，有長約6km，寬約1km的油污帶在西北側19km處被發(fā)現(xiàn)。這與模擬結果相類似。

油膜厚度圖例范圍為0～5mm，10d溢油軌跡與油膜厚度情況如圖4所示，從圖4中可以看出溢油油膜范圍基本與溢油軌跡重疊，溢油點處油膜厚度較厚呈紅色，距離溢油點較遠處的油膜厚度較薄呈綠色。從原理上分析，由于溢油點附近新冒出的溢油時間短，所以蒸發(fā)量少，擴展不充分，導致油膜厚度較厚；距溢油點較遠處的溢油由于溢出時間長，所以蒸發(fā)量大，擴展充分，導致油膜厚度較薄。這與模擬得到的結果相吻合。

由于每個油粒子的蒸發(fā)遵循相同的規(guī)則，所以選取第二個溢出的油粒子作為參考，根據(jù)數(shù)據(jù)顯示得到初始質(zhì)量為10000kg的油粒子在溢出后的第40小時質(zhì)量降到5%以下，為437.8kg。

圖4 10d溢油軌跡與油膜厚度情況

10d溢油軌跡與溶解物濃度情況如圖5所示，從圖5中可以發(fā)現(xiàn)溶解物濃度的分布與溢油軌跡并不重疊，其主要是在溢油點附近隨時間向四周擴散。從原理上分析，這是由于溢油點附近的油量大，所以溢油溶解物多，而遠離溢油點的地方油量少，所以溢油溶解物少，且風場對溶解物的漂移影響較小。

圖5 10d溢油軌跡與溶解物濃度情況

3 結語

本文簡述了流場模擬所使用的Hydroinfo水利信息平臺與所涉及的基本方法，其中包括基本方程、混合長度紊流模型、ALE描述法以及溢油各風化過程的模擬計算方法。對“桑吉輪”溢油事故模擬的整個過程進行了闡述，對油粒子漂移、擴展、溶解等進行了細致的展示，溢油的運動軌跡朝向北方并在擴散后期大體由兩條軌跡組成，這樣的計算結果與觀察到的現(xiàn)象也較吻合。由于新方法涉及參數(shù)的設定，在今后的工作中仍需在實例中不斷調(diào)整以選取更加合理準確的溢油參數(shù)。