金 戈，董仁杰，蔡 睿，陈红刚，金 生

(1.大连海事大学，辽宁 大连 116026；2.大连市水务集团水资源有限公司，辽宁 大连 116021；3.大连理工大学，辽宁 大连 116024)

海上溢油对海洋环境带来巨大破坏。能够准确快速地对溢油轨迹进行模拟预测对溢油事故的处理至关重要。一个简单、准确、有效的溢油模拟系统可以模拟溢油的运动轨迹、行为归宿与成分变化，从而对敏感地区溢油预防和灾情发生后即时处置提供重要支持。因此对溢油模型进行研究，开发完善的溢油模拟系统具有重要的理论与实际意义。

随着水环境模型与计算机技术的飞速发展，数值模拟计算软件的适用范围、计算能力、便利程度都在不断提高。很多成熟的水环境数值模拟软件随之涌现出来，如丹麦水力学研究所(DHI)的MIKE系列模拟软件、荷兰三角洲研究院(Deltares)的Delft 3D等。同样，溢油模拟也有许多较为成熟的软件，如美国国家海洋与大气研发中心(NOAA)研究开发的GNOME平台、ASA公司研发的OILMAP平台、挪威科技工业研究院的OSCAR、丹麦水环境研究院的MIKE/SA模块等。但它们大都需要通过其它渠道获得流场数据才能对溢油进行模拟，并没有水动力、溢油集成计算的功能。国内对于模拟软件研发工作的成果相对较少。本文搭建了包括水动力与溢油模拟的计算程序平台Hydroinfo，并使用该平台对近期发生的“桑吉轮”溢油事故进行了预测。

1 模型理论

1.1 Hydroinfo水利信息系统

Hydroinfo水利信息系统是一款计算复杂水流与输运问题的数值模拟软件。它能真实准确地模拟环境水流的多种形态，从河流海洋的一维、二维、三维近似，到非恒定的水流、波浪、泥沙与输运等各相关领域。其中包括三维自由水面流模拟、二维水流泥沙波浪问题、流域河网与管网、多维度耦合模块等，本文的潮流场则是基于三维自由水面流模拟模块进行计算。

1.2 水动力模型理论

1.2.1基本方程

如果直接对NS方程模拟则需要很高的时间与空间分辨率来获得不同尺度的湍流信息，导致计算量过大而无法计算。雷诺平均模拟和大涡模拟是常用模拟湍流运动的方法，其中雷诺平均模拟只计算平均运动，无需计算各尺度的湍流运动就可以得到湍流平均信息，计算量小，是目前工程中最广泛使用的模拟方法，也是本文所采用的方法。

连续性方程：

(1)

运动方程：

(2)

式中，ρ—流体密度；ν—运动粘度系数；fi—质量力强度；ui—速度，其中i=1、2、3。

Boussinesq假设认为湍流脉动所造成的附加应力，也类似层流运动应力，可以由时均的应变率与Boussinesq湍动粘度乘积表示。通过该假设得到雷诺应力与平均速度梯度的关系为：

(3)

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

湍动能k表达式：

(4)

在混合长度理论中Prandtl假定湍动粘度μt与时均速度ui的梯度和混合长度lm的乘积成正比，在二维问题中可以表示为：

(5)

混合长度lm是由实验确定或经验公式计算得出的值。

海洋中水平方向尺度远大于水深方向尺度，根据这一特点可以将方程进行简化，将水深方向应力与加速度忽略得到三维浅水Reynolds时均方程，将连续性方程沿水深方向积分，并通过Leibniz法则[7]可以得到水位方程。

(6)

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

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

自由表面的晃动往往会给数值模拟问题的求解带来困难，传统的描述方法通常使用拉格朗日描述法或欧拉描述法。纯欧拉描述与纯拉格朗日描述各有优缺点，而任意拉格朗日——欧拉描述法(ALE)则是结合两者特点产生的方法，该方法中生成的计算网格可以任何形式在空间运动，既可以独立于空间坐标系，也可以独立于物体坐标系。

使用ALE描述法对动边界处理有很大帮助，但同时会在方程中引入对流项使计算量增加，所以如何规定合适的网格运动方法是ALE描述中的关键问题。本文在自由表面法向方向使用拉格朗日描述以准确描述自由表面运动，在自由表面切线方向使用欧拉描述以防止网格发生扭曲。在这样的设定下网格只能沿自由表面法向做上下移动，既可以很好地描述动边界，也不会增加过多的计算量。

1.3 溢油模型理论

海上溢油的归宿是一个包含多种物理、化学变化的复杂过程，每个过程之间相互影响，共同组成了溢油的整个行为归宿。本文在“油粒子”法的基础上进行改进，将溢油视为由大量油粒子组成，并对每一个油粒子行为归宿进行模拟，所有油粒子的行为合在一起便模拟出整个溢油过程。

1.3.1扩展过程

本文是将Fay三段过程中的重力——粘性力阶段公式的微分形式应用于每一个油粒子上，模拟每一个油粒子单独的扩展过程，并不是先将整个溢油视作一个整体进行扩展过程模拟的传统方法，再拆分成多个油粒子进行后续过程的处理，而是直接对每一个油粒子在每个时间步长中的扩展情况进行模拟，从而得到每一时刻每个油粒子的扩展面积，这样可以针对持续性溢油进行模拟，同时也为后续过程模拟提供了必要的数据基础。通过量纲分析的方式得到扩展过程的计算公式，溢油表面积随时间变化率与溢油初始体积和溢油表面积有关，引入溢油扩展系数ks，计算表达式如下：

(7)

式中，A—油粒子表面积；t—时间；ks—溢油扩展系数；V—油粒子体积。

1.3.2迁移过程

本文从受力分析的角度出发，能够从原理上更加精确地对溢油迁移过程进行模拟。

油与水之间的摩擦力与海水密度、海水与溢油速度差的平方有关，所以通过量纲分析引入油与水摩擦系数fw。同理油与风之间的摩擦力与空气密度、风与溢油速度差的平方有关，引入油与风摩擦系数fa。溢油的受力与速度变化可表示为：

(8)

(9)

(10)

(11)

式中，τx、τy—溢油在x、y方向所受切应力；fw与fa—油与水摩擦系数和油与风摩擦系数；ρw与ρa—海水密度与空气密度；vw，x、vo，x、vw，y、vo，y—水与风x方向的速度和水与风y方向的速度；ax与ay—溢油在x、y方向的加速度；m—油粒子质量；A—油粒子表面积。

式(8)、(9)表示溢油所受水平切向力由溢油与海水的摩擦力和溢油与风的摩擦力构成，式(10)、(11)为牛顿定律，即加速度大小正比于合外力大小，与物体的惯性质量成反比。

1.3.3蒸发过程

本文通过量纲分析方式得到蒸发过程的计算公式，溢油蒸发体积量随时间变化率与溢油初始体积和溢油表面积有关，引入溢油蒸发系数ke，具体方程如下：

(12)

式中，A—油粒子表面积；t—时间；ke—溢油蒸发系数；V—油粒子体积。

1.3.4溶解过程

本文计算表面溶解量是通过量纲分析方式得到溶解过程的计算公式，溢油溶解量随时间变化率与溢油初始体积和溢油表面积有关，引入溶解扩展系数kd，计算表达式如下：

(13)

式中，A—油粒子表面积；t—时间；kd—溢油溶解系数；V—油粒子体积。

溶解物溶于水中后会随着水流运动，其运动满足物质运输的对流扩散方程，表达式为：

(14)

式中，S—溶解物浓度；t—时间；u、v、w—x、y、z三个方向的速度；εh、εv—水平方向与垂直方向的扩散系数。

2 “桑吉轮”溢油事故模拟

2.1 事故简介

2018年1月6日19∶50，“桑吉轮”与香港籍散货船“长峰水晶”轮发生剧烈碰撞，碰撞地点位于长江口以东约160n mile处，“桑吉轮”在碰撞后发生火灾，燃烧了整整8d后发生爆炸，随后沉入海底，伴随其沉入海底的还有13.6万t凝析油，船上32名船员全部遇难。“长峰水晶”号船员在碰撞后弃船登艇，之后被附近的浙岱渔03187船救起，该船船长郑磊等人目睹了整个碰撞过程。万幸的是长峰水晶轮最终在救助拖轮的护航下靠泊舟山港卸货。

2.2 模型参数

Hydroinfo水利信息系统与Google Earth有对接接口，可直接从Hydroinfo页面打开Google Earth并下载地图图片与地形数据。本文使用该功能下载地图图片作为底图，并获得其地形数据作为高程。由于底图尺寸与地形关系，还需在底图的基础上自行画出边界，圈出计算域，计算域如图1所示，计算时长为10d，从事故发生之日起。边界处的陆地部分可以设置为固壁，没有潮流流通也不允许油粒子穿越，而海面部分则可通过Hydroinfo中的潮位预报工具给出。

图1 计算域与高程二维展示

潮位预报通过对中国近海潮位数据进行傅里叶分析，从而得到该海域附近各点在各时间的潮位。海面边界潮位过程通过边界两端点的潮位过程插值给出，即使用两点潮位过程设置边界，十分快速便捷，避免了人工输入大量数据的繁琐。

本文的风场数据是从美国国家海洋和大气管理局(NOAA)地球系统研究实验室(ESRL)的物理科学部(PSD)网站上得到，风场0°为正东方向，90°为正北方向，风场数据见表1。

表1 风场数据

由于关心溢油点附近的变化情况，所以在该区域画了矩形加密线对其进行更密集的网格[12]。外边界的网格尺度为20000m，即平均每20000m生成一个网格，而加密线内区域的网格尺度被设置成5000m。最终生成16271个网格与8240个节点，生成的网格如图2所示，中间的矩形区域为加密区域。

图2 网格生成情况

参数设定包括物理参数与计算参数。其中物理参数包括：水位初值为0；海底糙率0.025；单个油粒子最大初始质量10000kg；溢油油品密度800kg/m3；空气密度1.293kg/m3；海水密度1025kg/m3；溢油挥发组分0.95；溢油挥发系数0.001/h；溢油溶解组分0.01；溢油溶解系数0.001/h；溢油面积扩展系数0.01/s。

控制参数包括：时间步长每步0.01h；计算步数24000步；垂向分层2层。

2.3 结果与分析

经过计算后得到溢油模拟的计算结果，由于篇幅有限，在此只展示10d的结果，如图3所示。中国大陆沿海的潮位变化最大，潮位最高处为高于海平面2m左右，潮位最低处为低于海平面2m左右。图中有部分区域始终显示为红色，如温州与福州沿海、台湾北部和琉球群岛，这是由于地形明显高于海平面导致的，并不是真实水位。

图3 10d潮流场

油粒子的运动轨迹是本次模拟最关心的内容。可以看出油粒子大体是向北方漂移并呈条带状，在10d溢油主要有两条轨迹。1月16日，国家海洋局中国海警在事故现场进行了监视监测，并在事故点附近发现油污带。9∶00，有长约9km，宽50～500m的油污带在北侧2km处被发现，10∶00，有长约6km，宽约1km的油污带在西北侧19km处被发现。这与模拟结果相类似。

油膜厚度图例范围为0～5mm，10d溢油轨迹与油膜厚度情况如图4所示，从图4中可以看出溢油油膜范围基本与溢油轨迹重叠，溢油点处油膜厚度较厚呈红色，距离溢油点较远处的油膜厚度较薄呈绿色。从原理上分析，由于溢油点附近新冒出的溢油时间短，所以蒸发量少，扩展不充分，导致油膜厚度较厚；距溢油点较远处的溢油由于溢出时间长，所以蒸发量大，扩展充分，导致油膜厚度较薄。这与模拟得到的结果相吻合。

由于每个油粒子的蒸发遵循相同的规则，所以选取第二个溢出的油粒子作为参考，根据数据显示得到初始质量为10000kg的油粒子在溢出后的第40小时质量降到5%以下，为437.8kg。

图4 10d溢油轨迹与油膜厚度情况

10d溢油轨迹与溶解物浓度情况如图5所示，从图5中可以发现溶解物浓度的分布与溢油轨迹并不重叠，其主要是在溢油点附近随时间向四周扩散。从原理上分析，这是由于溢油点附近的油量大，所以溢油溶解物多，而远离溢油点的地方油量少，所以溢油溶解物少，且风场对溶解物的漂移影响较小。

图5 10d溢油轨迹与溶解物浓度情况

3 结语

本文简述了流场模拟所使用的Hydroinfo水利信息平台与所涉及的基本方法，其中包括基本方程、混合长度紊流模型、ALE描述法以及溢油各风化过程的模拟计算方法。对“桑吉轮”溢油事故模拟的整个过程进行了阐述，对油粒子漂移、扩展、溶解等进行了细致的展示，溢油的运动轨迹朝向北方并在扩散后期大体由两条轨迹组成，这样的计算结果与观察到的现象也较吻合。由于新方法涉及参数的设定，在今后的工作中仍需在实例中不断调整以选取更加合理准确的溢油参数。