內河船舶溢油建模與仿真

熊正華，陳 茜

(1. 四川交通職業技術學院 航運工程系，四川 成都 611130； 2. 武漢理工大學 能源與動力工程學院，湖北 武漢 430063)

0 引 言

近年來溢油事故頻發，對經濟發展以及生態環境都造成了嚴重危害，研究船舶溢油的影響因素對于預防溢油事故的發生有著重要意義[1]。目前對于溢油行為的研究多集中于海洋，針對內河環境的溢油研究相對較少。內河水流流速較快且徑流量變化迅速，航道相對狹窄，并且內河作為人們生活工作的主要場所和生活用水的主要來源，內河溢油事故一旦發生，造成的危害更大、顯現的速度更快[2]。筆者在船舶溢油影響因素分析的基礎上，針對內河航道的特殊性選取江西撫河部分河段進行了水動力學模擬，并以流場模擬結果為基礎對溢油行為進行了模擬。

1 水動力模型

根據內河水文特征選取MIKE21的HD模塊進行水動力學模擬。MIKE21HD模塊采用交替方向隱格式法(ADI)離散水流的控制方程，以中心差分格式表示方程中的各項參數，采用追趕法對離散得到的矩陣方程進行求解，從而避免方程因為離散而產生計算結果失真的情況，計算過程中Taylor級數展開的截斷誤差在二階至三階以內[3]。

二維淺水連續性方程為

(1)

式中：ζ為水位；p為x方向的流量通量；q為y方向的流量通量；h為水深。x方向動量方程為

(2)

y方向動量方程為：

(3)

式中：H=h+ζ；C為謝才系數；ρ為流體密度；τxx、τxy、τyy為有效剪切力分量；f為科氏力系數；fω為風阻力系數；Wx、Wy分別為風速在x、y方向上的分量。

2 溢油模型

2.1 輸移過程

2.1.1 擴散過程

在溢油進入水體的初始階段，油膜較厚，油層在進行漂移的同時，油膜面積迅速擴大，厚度減少直至油膜到臨界值最終破裂。MIKE21SA模塊采用修正過的Fay公式來描述溢油的擴展運動，修正過的Fay公式為[4]

(4)

2.1.2 漂移擴散過程

在風力以及水流的作用下，油粒子出現相應的對流位移，油粒子的漂移速度表示為

UP=US+Cω×Uω×sin(θ-π+θω)

(5)

VP=VS+Cω×Uω×cos(θ-π+θω)

(6)

式中：UP和VP分別為油粒子的對流速度在x、y方向上的分量；US和VS分別為水流的表面流速在x、y方向上的分量；Cω為風漂移系數；Uω為距離計算點10 m處的水面風速；θ為風向角；θω為風偏轉角。

MIKE21SA模塊應用雙線內插法對油粒子的對流速度進行插值處理，解決油粒子無法恰好處于離散網格流速計算點的問題，原理如圖1。

圖1 流速內插法Fig. 1 Velocity interpolation method

油粒子的插值速度可以表示為

V=V1+(V4-V1)x+(V2-V1)y+(V1-V2+V3-V4)xy

(7)

式中：V1、V2、V3以及V4為離散網格點的已知流速；x和y為計算點到某網格點的距離。

水流在水平方向上的擴散距離可為

(8)

2.2 風化過程

2.2.1 蒸發過程

蒸發是溢油中的各項組分由液態氣化進入大氣的過程。溢油的蒸發率可表示為[5]

(9)

(10)

式中：k為蒸發系數；Sci為溢油組分i的Schmidts數。

2.2.2 溶解過程

蒸發是溢油中的各項組分由液態氣化進入大氣的過程。溢油的蒸發率可表示為[6]

(11)

2.2.3 吸附和沉降過程

從概率角度出發，油粒子被河岸吸附的概率可表示為

(12)

式中：Amax為河岸的最大吸附能力，Amax的大小與河岸類型有關；A為該區域目前已經吸附的溢油量。

3 溢油事故數值模擬實驗

3.1 仿真模型

選取江西撫河部分河段為計算區域應用MIKE21HD模塊建立二維水動力學模型[7]，并分析3種不同水流條件下的水位和流場模擬結果，用到的水流數據如表1。

表1 計算區域進出口流量及水位Table 1 Inlet and outlet discharge and water level at calculated region

時間步長進行計算分析，選取Δt=30 s。3種水流條件下的模擬情況分別如圖2～圖4。

圖2 50年一遇洪水模擬結果Fig. 2 Simulation results of a 50-year-once flood

圖3 20年一遇洪水模擬結果Fig. 3 Simulation results of a 20-year-once flood

圖4 10年一遇洪水模擬結果Fig. 4 Simulation results of a 10-year-once flood

3.2 模型驗證

為保證撫河二維水動力學模型的有效性，給溢油模型的模擬計算提供正確基礎，針對以上3種不同的水文情況，選取計算區域上的10個斷面，將實測水位與模擬水位進行比較驗證，考察上述水流模擬結果的有效性，斷面位置選擇如圖5。

圖5 驗證斷面Fig. 5 Validation section

斷面4～斷面13的實測水位以及由MIKE21HD模塊所建立的二維水動力學模型模擬出的模擬水位比較如圖6～圖8。

圖6 50年一遇洪水模擬驗證情況Fig. 6 Simulation and verification of 50-year-once flood

由圖6～圖8的比較情況來看，水位模擬結果與實測水位誤差在5%以內，且3種水文條件下各斷面水位均呈升高趨勢，以往文獻通常根據模擬水位與實測水位的誤差大小判斷模擬結果的可靠性，但是該河段水位變化較小，故僅依據模擬值與實測值之間的誤差大小難以說明模擬結果的可靠性。

圖7 20年一遇洪水模擬驗證情況Fig. 7 Simulation and verification of 20-year-once flood

圖8 10年一遇洪水模擬驗證情況Fig. 8 Simulation and verification of 10-year-once flood

在數學統計中較常應用于描述模擬結果與實際值擬合情況的參數為可決系數R2，可決系數R2是描述模型模擬結果與實測值的擬合優度的參數，被廣泛應用于水利工程計算的檢驗，R2的取值范圍為0～1，R2越接近1，證明擬合效果越好，根據SL104—2015《水利工程水力計算規范》，通常認為R2>0.8為模擬結果有效。

對于一個包含k個解釋變量的計算樣本和測量樣本，有：

(13)

(14)

式(14)中：ESS為可解釋的平方和，即：

(15)

TSS為總離差平方和，即：

(16)

RSS為剩余平方和，即：

(17)

(18)

由圖5～圖7計算得到的可決系數以及相關擬合優度曲線圖如圖9～圖11。

由圖9～圖11可知，3種水文條件的模擬結果均滿足SL104—2015《水利工程水力計算規范》中對于R2>0.8的要求，因此該撫河河段水流模擬具備有效性[8]。

圖9 50年一遇洪水擬合優度曲線Fig. 9 The goodness-of-fit curve of 50-year-once flood

圖10 20年一遇洪水擬合優度曲線Fig. 10 The goodness-of-fit curve of 20-year-once flood

圖11 10年一遇洪水擬合優度曲線Fig. 11 The goodness-of-fit curve of 10-year-once flood

本節介紹了應用MIKE21HD模塊進行水動力模擬的基本方程，以及求解基本方程的數值離散方法，以撫河部分河段為例，建立了二維水動力學模型，分別模擬了3種水流條件下的水位情況和流場情況，并根據實測數據驗證了該水動力學模型的正確性。

4 內河溢油建模與仿真

4.1 溢油模擬

本節對溢油模型的各項參數進行反復率定，各項參數的最終確定結果如表2。

表2 溢油模型參數取值Table 2 Parameter value of oil spill model

本節采用油粒子模型對溢油進行概化，即將溢油離散為有限數量的油滴粒子，每個粒子均代表一定量的溢油，并具有相應的油膜厚度、坐標、運動速度等物理量[8]。

溢油模擬的計算區域與二維動力學模型計算區域相同。選取水流條件以及風場兩個主要影響因素對撫河部分河段溢油行為進行模擬。根據撫州市氣象局的歷史風力統計結果，選取無風、東南風4.5 m/s以及西北風7 m/s，3種風況情況進行模擬。因此，本節中溢油模型針對如下所述9種工況展開模擬，工況條件整理如表3。

表3 計算工況Table 3 Calculation condition

4.2 溢油行為模擬

本次溢油事件選取溢油點網格坐標為(66，66)的位置進行模擬，溢油速度為2 m3/s，溢油持續時間為30 min，即溢油量為3 600 m3。

1)工況1模擬結果如圖12。由圖可知，溢油事故發生到13 h 20 min時，溢油開始到達計算區域出水口斷面，少量溢油往主流入水口右岸岸線漂移，無溢油到達支流入水口。因此，在此種風況下計算區域出水口遭受污染嚴重，主流入水口右岸岸線受到少量污染，支流入水口幾乎不受溢油影響。

圖12 工況1模擬情況Fig. 12 Simulation of case 1

2)工況2模擬結果如圖13。由圖可知，溢油事故發生到1 h 45 min時，溢油由于東南向4.5 m/s風逐漸向左漂移，開始到達主流入水口左岸，隨著時間的推移，溢油逐漸往主流入水口左岸岸線聚集，隨后沿主流入水口左岸岸線向計算區域出水口溢出。因此，在此種風況下主流入水口左岸岸線遭受污染嚴重，支流入水口幾乎不受溢油影響。

圖13 工況2模擬情況Fig. 13 Simulation of case 2

3)工況3模擬結果如圖14。由圖可知，溢油事故發生55 min時，溢油受到來自西北向7 m/s風而向東南方向漂移，開始到達主流入水口右岸，隨著時間的推移，溢油逐漸往東北方向漂移，即主流入水口右岸岸線聚集，隨后沿主流入水口右岸岸線向支流入水口溢出。因此，在此種風況下主流入水口右岸岸線、支流入水口遭受污染嚴重，計算區域出水口幾乎不受溢油影響。

圖14 工況3模擬情況Fig. 14 Simulation of case 3

4)工況4模擬結果如圖15：由圖可知，溢油事故發生到13 h 30 min時，溢油開始到達計算區域出水口斷面，對比工況1而言溢油更集中，少量溢油往主流入水口右岸岸線漂移，無溢油到達支流入水口。因此，在此種風況下計算區域出水口遭受污染嚴重，主流入水口右岸岸線受到少量污染，支流入水口幾乎不受溢油影響。

5)工況5模擬結果如圖16：由圖可知，溢油事故發生到1 h 45 min時，溢油由于東南向4.5 m/s風逐漸向左漂移，對比工況2而言溢油更集中，開始到達主流入水口左岸，隨著時間的推移，溢油逐漸往主流入水口左岸岸線聚集，隨后沿主流入水口左岸岸線向計算區域出水口溢出。因此，在此種風況下主流入水口左岸岸線遭受污染嚴重，支流入水口幾乎不受溢油影響。

圖15 工況4模擬情況Fig. 15 Simulation of case 4

圖16 工況5模擬情況Fig. 16 Simulation of case 5

6)工況6模擬結果如圖17：由圖可知，溢油事故發生到55 min時，溢油受到來自西北向7 m/s風而向東南方向漂移，較為集中的開始到達主流入水口右岸；隨著時間的推移，在溢油事故發生2 h 30 min時，溢油逐漸往東北方向漂移，即主流入水口右岸岸線聚集，隨后沿主流入水口右岸岸線沿岸向支流入水口溢出。因此，在此種風況下主流入水口右岸岸線、支流入水口遭受污染嚴重，計算區域出水口幾乎不受溢油影響。

圖17 工況6模擬情況Fig. 17 Simulation of case 6

7)工況7模擬結果如圖18: 由圖可知，溢油事故發生5 h時，溢油在事故點逐漸散開，但對比工況1而言溢油更集中；溢油事故發生13 h 20 min時，溢油開始到達計算區域出水口斷面，少量溢油往主流入水口右岸岸線漂移，無溢油到達支流入水口。因此，在此種風況下計算區域出水口遭受污染嚴重，主流入水口右岸岸線受到少量污染，支流入水口幾乎不受溢油影響。

圖18 工況7模擬情況Fig. 18 Simulation of case 7

8)工況8模擬結果如圖19:由圖可知，溢油事故發生到1 h 45 min時，溢油由于東南向4.5 m/s風逐漸向左漂移，開始到達主流入水口左岸；溢油事故發生到3 h 30 min時，溢油逐漸往主流入水口左岸岸線聚集，隨后沿主流入水口左岸岸線向計算區域出水口溢出。因此，在此種風況下主流入水口左岸岸線遭受污染嚴重，支流入水口幾乎不受溢油影響。

圖19 工況8模擬情況Fig. 19 Simulation of case 8

9)工況9模擬結果如圖20:由圖可知，溢油事故發生到55 min時，溢油受到來自西北向7 m/s風而向東南方向漂移，開始到達主流入水口右岸，溢油事故發生2 h 30 min時，溢油逐漸往東北方向漂移，即主流入水口右岸岸線聚集，隨后沿主流入水口右岸岸線向支流入水口溢出。因此，在此種風況下主流入水口右岸岸線、支流入水口遭受污染嚴重，計算區域出水口幾乎不受溢油影響。

圖20 工況9模擬情況Fig. 20 Simulation of case 9

4.3 溢油模擬結果分析

風況對油膜漂移擴散軌跡影響顯著。在無風條件下，油膜跟隨水流運動向計算區域出水口斷面溢出；在東南向4.5 m/s風況下，油膜向主流左岸岸線移動，最終從計算區域出水口溢出；在西北向7 m/s風況下，油膜向主流右岸岸線移動，最終從支流入水口溢出。

與風場的影響相比，流場對油膜漂移擴散的影響相對較小。流速越大，油膜漂移擴散速度越快。在局部流場特殊的區域，油膜的漂移擴散也受到了一定影響。

5 結 語

針對內河溢油事故的特殊性，分析了船舶溢油事故的影響因素，選取江西撫河部分河段為研究對象。利用MIKE21HD模塊對江西撫河部分河段的流場進行了模擬，并根據實測數據驗證了水動力學模型的正確性，為溢油模擬提供流場基礎。并在水動力模型的基礎上，利用油粒子模型對3種水流條件、3種風場條件共9種計算工況下的溢油行為進行了模擬，并得出水流條件和風場條件對溢油行為的影響規律，風場對于溢油行為的影響最為顯著，因此提高風況的測量精度，并對研究區域的風場進行高精度模擬對于溢油行為的研究具有重要意義。