內(nèi)河船舶溢油建模與仿真

熊正華，陳 茜

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

0 引 言

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

1 水動力模型

根據(jù)內(nèi)河水文特征選取MIKE21的HD模塊進(jìn)行水動力學(xué)模擬。MIKE21HD模塊采用交替方向隱格式法(ADI)離散水流的控制方程，以中心差分格式表示方程中的各項(xiàng)參數(shù)，采用追趕法對離散得到的矩陣方程進(jìn)行求解，從而避免方程因?yàn)殡x散而產(chǎn)生計(jì)算結(jié)果失真的情況，計(jì)算過程中Taylor級數(shù)展開的截?cái)嗾`差在二階至三階以內(nèi)[3]。

二維淺水連續(xù)性方程為

(1)

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

(2)

y方向動量方程為：

(3)

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

2 溢油模型

2.1 輸移過程

2.1.1 擴(kuò)散過程

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

(4)

2.1.2 漂移擴(kuò)散過程

在風(fēng)力以及水流的作用下，油粒子出現(xiàn)相應(yīng)的對流位移，油粒子的漂移速度表示為

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

(5)

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

(6)

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

MIKE21SA模塊應(yīng)用雙線內(nèi)插法對油粒子的對流速度進(jìn)行插值處理，解決油粒子無法恰好處于離散網(wǎng)格流速計(jì)算點(diǎn)的問題，原理如圖1。

圖1 流速內(nèi)插法Fig. 1 Velocity interpolation method

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

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

(7)

式中：V1、V2、V3以及V4為離散網(wǎng)格點(diǎn)的已知流速；x和y為計(jì)算點(diǎn)到某網(wǎng)格點(diǎn)的距離。

水流在水平方向上的擴(kuò)散距離可為

(8)

2.2 風(fēng)化過程

2.2.1 蒸發(fā)過程

蒸發(fā)是溢油中的各項(xiàng)組分由液態(tài)氣化進(jìn)入大氣的過程。溢油的蒸發(fā)率可表示為[5]

(9)

(10)

式中：k為蒸發(fā)系數(shù)；Sci為溢油組分i的Schmidts數(shù)。

2.2.2 溶解過程

蒸發(fā)是溢油中的各項(xiàng)組分由液態(tài)氣化進(jìn)入大氣的過程。溢油的蒸發(fā)率可表示為[6]

(11)

2.2.3 吸附和沉降過程

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

(12)

式中：Amax為河岸的最大吸附能力，Amax的大小與河岸類型有關(guān)；A為該區(qū)域目前已經(jīng)吸附的溢油量。

3 溢油事故數(shù)值模擬實(shí)驗(yàn)

3.1 仿真模型

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

表1 計(jì)算區(qū)域進(jìn)出口流量及水位Table 1 Inlet and outlet discharge and water level at calculated region

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

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

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

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

3.2 模型驗(yàn)證

為保證撫河二維水動力學(xué)模型的有效性，給溢油模型的模擬計(jì)算提供正確基礎(chǔ)，針對以上3種不同的水文情況，選取計(jì)算區(qū)域上的10個(gè)斷面，將實(shí)測水位與模擬水位進(jìn)行比較驗(yàn)證，考察上述水流模擬結(jié)果的有效性，斷面位置選擇如圖5。

圖5 驗(yàn)證斷面Fig. 5 Validation section

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

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

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

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

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

在數(shù)學(xué)統(tǒng)計(jì)中較常應(yīng)用于描述模擬結(jié)果與實(shí)際值擬合情況的參數(shù)為可決系數(shù)R2，可決系數(shù)R2是描述模型模擬結(jié)果與實(shí)測值的擬合優(yōu)度的參數(shù)，被廣泛應(yīng)用于水利工程計(jì)算的檢驗(yàn)，R2的取值范圍為0～1，R2越接近1，證明擬合效果越好，根據(jù)SL104—2015《水利工程水力計(jì)算規(guī)范》，通常認(rèn)為R2>0.8為模擬結(jié)果有效。

對于一個(gè)包含k個(gè)解釋變量的計(jì)算樣本和測量樣本，有：

(13)

(14)

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

(15)

TSS為總離差平方和，即：

(16)

RSS為剩余平方和，即：

(17)

(18)

由圖5～圖7計(jì)算得到的可決系數(shù)以及相關(guān)擬合優(yōu)度曲線圖如圖9～圖11。

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

圖9 50年一遇洪水?dāng)M合優(yōu)度曲線Fig. 9 The goodness-of-fit curve of 50-year-once flood

圖10 20年一遇洪水?dāng)M合優(yōu)度曲線Fig. 10 The goodness-of-fit curve of 20-year-once flood

圖11 10年一遇洪水?dāng)M合優(yōu)度曲線Fig. 11 The goodness-of-fit curve of 10-year-once flood

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

4 內(nèi)河溢油建模與仿真

4.1 溢油模擬

本節(jié)對溢油模型的各項(xiàng)參數(shù)進(jìn)行反復(fù)率定，各項(xiàng)參數(shù)的最終確定結(jié)果如表2。

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

本節(jié)采用油粒子模型對溢油進(jìn)行概化，即將溢油離散為有限數(shù)量的油滴粒子，每個(gè)粒子均代表一定量的溢油，并具有相應(yīng)的油膜厚度、坐標(biāo)、運(yùn)動速度等物理量[8]。

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

表3 計(jì)算工況Table 3 Calculation condition

4.2 溢油行為模擬

本次溢油事件選取溢油點(diǎn)網(wǎng)格坐標(biāo)為(66，66)的位置進(jìn)行模擬，溢油速度為2 m3/s，溢油持續(xù)時(shí)間為30 min，即溢油量為3 600 m3。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3 溢油模擬結(jié)果分析

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

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

5 結(jié) 語

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