內河溢油污染源逆時追蹤數值模擬

付 玉,童思陳,2,王 嘯,蔣聘鳳,黃國鮮

(1.重慶交通大學 河海學院,重慶 400074; 2.國家內河航道整治工程技術研究中心,重慶 400074; 3.河南省交通規劃設計研究院有限公司,鄭州 450000; 4.中國環境科學研究院,北京 100012)

0 引 言

隨著我國水路運輸蓬勃發展,內河船舶運輸成為推動經濟增長的重要力量[1]。由于內河航道相對狹窄,船舶海損事故時有發生,或導致石油、燃油泄漏,對水域環境造成污染。自20世紀60年代以來,溢油研究多針對海洋區域,且主要關于溢油正向擴展漂移模擬。孫蘭君等[2]基于蒙特卡羅方法及Mia散射理論建立了溢油海水雙分布函數模型,為海洋激光遙感測油提供了理論基礎;Wang等[3]基于“粒子法”開發了包含海洋環境預測的近海溢油預測系統;Simecek-Beatty等[4]結合空間驗證的技巧評分法進行了石油泄漏預測,研究了海上溢油預測的全新方式。近年來,通過對河道中溢油現象的研究,學者們發現與海洋環境不同,油粒子在擴散過程中更容易受河岸和河底的固壁作用影響[5],主要受水流對流作用控制。Gautama等[6]將合成孔徑雷達(Synthetic Aperture Radar,SAR)遙感技術應用于河道溢油模型的驗證;蔣文燕[7]在二維模型基礎上,考慮增加垂向尺度油滴擴散和上浮作用對垂向油滴濃度分布的影響,建立了河道三維溢油數學模型;黃瑞等[8]將改進油粒子模型應用于長江鎮揚段上,分析了溢油擴散特征和油膜厚度變化,并考慮了溢油風化及岸邊吸附等影響;Jacketti等[9]拓展了開源地下石油模擬器(Subsurface Oil Simulator,SOSim)用于預測河流底部緩慢移動的石油軌跡。

相對于溢油污染正向預測和預報,針對污染事件反向追蹤溢油源的研究較為缺乏。辛小康等[10]采用優化算法建立了河流水污染事故污染源識別模型;楊海東等[11]結合MCMC(Markov Chain Monte Carlo)和微分進化方法提出河流突發性水污染的溯源方法; Soong等[12]使用隨機行走粒子跟蹤算法模擬油粒子聚集物在河流中的輸移特性并開發了FluOil模型。大多數學者將污染追溯問題看作源項反演問題,采用概率法或模擬優化法等反演污染源的位置和歷史運動軌跡。國內外開發的具有污染源上游回溯功能的模型有RiverSpill溢油模型和OILMAP數值模擬系統,其中RiverSpill更適用于河流污染情況,但其使用的是一維流場計算方法[13],OILMAP具有二維及三維計算能力,多用于模擬溢油與岸線、海床等的相互作用[14]。除此以外,目前針對內河溢油污染源的逆時追蹤研究很少,關于多維度模擬狹長河道復雜地形的溯源追蹤有待進一步研究深化。本文在平面二維水流數值模擬基礎上,基于Fortran語言從底層自主開發了溢油污染源逆時追蹤模型,可實現溢油污染過程的逆向運算,反演溢油的漂移路徑和尋找溢油污染源的位置。通過在三峽庫區典型河段的溢油污染源逆時追蹤數值模擬,可較好地模擬尋找出溢油污染源,可對溢油污染事故的應急處置和污染源的確定提供支撐。

1 模型及求解

內河具有河寬較窄、水深較淺、平面形態復雜等特點,河道溢油主要受水流流速控制,對流作用突出。溢油污染數學模型是在水流數學模型基礎上建立的,本文采用的二維水流數學模型基于課題組前期研究成果[15]。計算采用貼體正交曲線網格,差分方程采用三對角矩陣算法求解,數值模擬采用SIMPLEC算法,并使用欠松弛技術促進非線性迭代的收斂。

在溢油應急時間尺度下[16],基于拉格朗日方法建立了溢油擴展漂移油粒子模型,將油膜離散為1 000顆(可根據需要自行設置)直徑為10～1 000 μm的球形油粒子,油膜的擴散和漂移被視為多個油粒子在外部驅動力作用下的疊加效應。油粒子模型可以綜合流場和河岸地形的影響,適用于模擬水文條件較為復雜的河流中的溢油事故。本文重點研究瞬時溢油導致的短期環境影響,在此期間溢油的運動主要為自身的擴展和漂移運動,重力和黏性力起主導作用[17]。

1.1 溢油正向擴展漂移模型

1.1.1 擴展過程

根據FAY理論,擴展主要分3個階段:重力-慣性力擴展階段、重力-黏性力擴展階段和表面張力-黏性力擴展階段[18]。本文主要考慮在重力和黏性力作用下的油膜擴展情況,即采用Mackay修正的FAY理論重力-黏性力公式計算油膜的擴展,即

(1)

V0=πR02h0。

(2)

式中:A為油膜表面積(m2);t為時間(s);K0為油膜面積增長率;V0為溢油體積(m3);R0為油膜半徑(cm);h0為油膜初始厚度(cm)。

1.1.2 漂移過程

溢油發生后,油膜自身會做擴展運動,還會在水流、波浪和風等外界動力的驅動下漂移擴散。漂移過程包括平移過程和隨機擴散過程。在平移過程中,對于庫區河道水流,影響漂移的主要動力因素是水流和風,油粒子的平移速度可由式(3)計算,即

vp=vf+avw。

(3)

式中:vp為油粒子平移速度(m/s);vf為水流速度(m/s);a為風力因子,一般取0.03～0.04;vw為水面以上10 m處風速(m/s)。擴散過程中,油粒子在水中的運動是具有隨機性的湍流彌散過程。本文采用生成隨機數的方法來計算油粒子可能擴散距離,表達式為

(4)

式中:S0為擴散距離(m);R為利用函數生成的[-1,1] 區間內的均勻分布的隨機數;D0為水平方向上的擴散系數;Δt為時間步長(s)。

1.1.3 油粒子的單步距離

一個時間步長內單顆油粒子的單步距離St計算方法為

St=vpΔt+S0=(vf+avw)Δt+S0。

(5)

1.1.4 油粒子的流速插值

水流數學模型計算過程中,采用了貼體正交曲線網格,首先可得到每個網格節點的流速。在二維歐拉法水流數學模型基礎上,采用拉格朗日法把每個油粒子視作獨立的模擬對象進行模擬反演,其原理如圖1所示。在水動力模塊計算中只能得到每個網格節點流速,而每顆油粒子的位置是隨機的,因此可先根據面積相等法編程,對油粒子進行搜索定位,找出每顆油粒子所在網格中的位置。比如,假定經搜索,確定某顆油粒子O(x,y)在某個網格A1A2A3A4之中(如圖1所示)。

圖1 油粒子流速插值原理Fig.1 Schematic of oil particle velocity interpolation

各網格節點的坐標如圖1,網格中心點的坐標為C(xc,yc)。設各網格節點坐標分別為A1(x1,y1)、A2(x2,y1)、A3(x2,y2)、A4(x1,y2),計算網格邊長分別為L1、L2,設坐標變換函數為

(6)

式中:x、y是原始坐標系中油粒子O的坐標;ξ、η是經過坐標變換后油粒子O在貼體正交曲線坐標系中的坐標。

式中ξi、ηi是網格上各頂點Ai的局部坐標。插值函數f(xi)為

(8)

于是,由二維水流數學模型計算可得出每個網格節點A1、A2、A3、A4的流速vx、vy,利用式(6)—式(8)即可插值得出該顆油粒子O(x,y)的流速vxO、vyO。

1.2 溢油污染源逆時追蹤模型

1.2.1 實現思路

(1)首先,獲取溢油污染發生時刻末的污染范圍和面積,并確定時刻末溢油油品和油粒子位置。

(2)進行平面二維水動力數值模擬并獲取每個網格節點的水動力參數。

(3)自時刻末起,對每顆油粒子進行逆時追蹤反演計算,模擬油粒子反方向的漂移軌跡,從而確定溢油污染源,具體實現步驟:①對油粒子進行逆時追蹤運算,每顆油粒子前一時刻位置坐標為

(9)

式中:X0和Y0為油粒子當前時刻的位置坐標;X和Y為油粒子前一時刻的位置坐標;dt為時間步長;vx0和vy0分別為x和y方向的水流流速;vwx和vwy分別為x和y方向的風速;S0x和S0y分別為油粒子在一個步長時間內x和y方向的擴散距離。

②計算所追蹤油粒子區域面積,若所得面積大小滿足污染源判別標準,則將此區域作為溢油污染源;若不滿足,則重新進行計算。追蹤過程如圖2所示。

圖2 溢油污染源逆時追蹤示意圖Fig.2 Inverse time tracking of oil spill pollution source

1.2.2 判別標準

1.2.2.1 油膜面積的計算

由油粒子模型可知,油膜可視為大量油粒子組成的云團,即油粒子群。其面積可采用任意多邊形面積公式進行計算,即

(10)

其中,每時刻最外層的油粒子坐標分別為A1(x1,y1)、A2(x2,y2)、…、An(xn,yn)。

1.2.2.2 污染源判定

根據溢油的正向擴展漂移運動特征可知,在反向追蹤溢油源頭的過程中,將會顯示溢油面積逐漸縮小且油粒子逐漸聚攏的過程,在模型中可預設一個油膜厚度閾值dmax。研究表明當油膜厚度為1 mm時,溢油在水中的擴展過程將停止[19],因此可初步設定dmax=1 mm。再根據油膜厚度反算油膜面積Smin,當某時刻溢油污染面積S≤Smin時,則可認為已尋找到溢油污染源。

2 模型驗證

三峽庫區變動回水區洛磧段船舶易發生海損事故,溢油風險較大,因此將洛磧段作為典型河段進行模擬。

2.1 水流模型驗證

利用2015年5月28日實測的水位和流速對二維水流模型進行了驗證(圖3),實測流量為6 180 m3/s?？梢?水位、流速計算值與實測值均符合良好。

圖3 水位、流速分布驗證Fig.3 Validation of water level and flow velocity

2.2 溢油正向擴展漂移模型驗證

課題組開展了溢油擴展漂移的概化水槽試驗[20],采用其中一組工況對溢油數學模型進行了驗證(圖4、圖5)。其中,水流速度v為0.25 m/s,溢油量m為31.0 g?？芍?油膜面積、厚度模擬值與試驗值均吻合良好,偏差在±15%之內。

圖4 溢油擴展漂移模擬與室內試驗過程比較Fig.4 Comparison of oil spill spread and drift between simulation and indoor test

圖5 溢油數學模型與室內試驗油膜面積與厚度的比較Fig.5 Comparison of area and thickness of oil film between mathematical modelling and indoor test

3 溢油逆時追蹤模擬

3.1 流場模擬

采用洛磧段為典型河段,模擬計算范圍長度約5.1 km,寬度約1.0 km。選取洪水期流量為22 000 m3/s,尾水位為170.22 m。經分析,該河段中洪水期與平水期流量為6 180 m3/s的河道條件并無明顯差別,因此借鑒了前文率定的糙率。

3.2 污染源追蹤模擬

設油品為潤滑油,擴散系數D0=3.0×10-4,溢油量為1 000 kg,溢油密度為940.5 kg/m3,風速為1.5 m/s。時刻末(溢油發現時刻)油膜面積為70 440 m2,油粒子顆粒數設為1 000,計算時間步長為10 s。模擬可得到溢油污染源逆時追蹤過程如圖6所示,相應各典型時刻的溢油追蹤過程特征參數如表1所示。

表1 不同追蹤時刻溢油油膜面積和厚度特征統計Table 1 Areas and thicknesses of oil spill film at different moments

注:數字表示沿邊界方向分布的網格線數量。圖6 模擬的溢油逆時追蹤過程軌跡及局部軌跡Fig.6 Local trajectories of inverse time tracking simulation

溢油源油膜厚度閾值dmax=1 mm,反算得到預設油膜面積Smin=1 063.26 m2。當模擬得到油膜面積S

3.3 溢油事故應急處置建議

溢油事故發生后,或將嚴重危害河流水體環境、水中動植物、河岸灘涂及沿岸的人類生產生活,若能及時尋找溢油事故源頭,可為溢油事故應急處置提供有效指導。一旦發現河道溢油事故,應立即啟動應急響應機制,根據溢油泄漏量、發生地點、河流流量等劃分事故風險等級,制定相應的溢油處理方案?？蛇x擇水流流速較小或支流交匯口等位置作為溢油圍控點[21],及時布設圍油欄或攔污浮筒,在取水口和自然保護區等敏感水域設置多層油柵[22],或采用浮標、阻隔膜和重力錨組成的軟隔離設施,便于隔離和收集漂浮污染物與可溶或沉積污染物[23],減小油膜的擴散范圍。當流速較大時,可使用高速撇油器回收油污,然后利用吸附劑或化學清潔劑對河面進行清理。此外,也可采用原位燃燒法[24]、表面清洗劑[25]、分散劑[26]和生物修復[27]等方法修復含油污水。

4 結 論

在溢油應急時間尺度范圍內,基于平面二維水流數學模型,采用Fortran語言從底層開發建立了溢油污染擴展漂移及污染源逆時追蹤模型。在水流數學模型和溢油擴展漂移數學模型驗證基礎上,考慮了內河水動力對流作用突出的特點,選取了三峽庫區洛磧河段作為典型河段進行了溢油污染源的模擬反演。結果表明,各油粒子之間距離逐漸減小,油粒子逐漸聚攏。能夠反演重現溢油污染物隨水流的漂移擴散過程,通過預設油膜面積閾值可判斷出污染源位置。

溢油模型具有逆推功能,可以把其應用在內河船舶溢油事件的航政調查決策。根據特定地點所觀察到的油污情況逆推溢油的可能來源,可為不明溢油源污染事故調查提供參考,以便及時提出事故處理預案,防止油污的進一步擴散,減小溢油污染對河流環境造成的不良影響。