内河船舶溢油建模与仿真

熊正华，陈 茜

(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种计算工况下的溢油行为进行了模拟，并得出水流条件和风场条件对溢油行为的影响规律，风场对于溢油行为的影响最为显著，因此提高风况的测量精度，并对研究区域的风场进行高精度模拟对于溢油行为的研究具有重要意义。