渤海湾海上瞬时溢油三维模拟

张明 任鸿翔 朱天惠 周毅 唐海娜

摘要：为实时模拟及预测渤海湾三维海上溢油的动态特征，采用Mike 21 Flow模型计算渤海水动力潮汐，针对瞬时溢油的扩散和漂移过程，建立海上溢油运动数学模型；
采用泊松盘采样算法创建不同直径的油粒子，通过Marching Cube算法重建等值面，模拟实现海上溢油过程，并与真实案例的数据对比。溢油模拟结果表明：所建潮汐流数值模型评估溢油事故船舶位置的海面高程、海流速度和海流方向与观测站数据接近；
海上溢油运动数学模型原油运动轨迹与与实际轨迹相似；
相较于传统溢油模拟可视化工作局限于二维非交互效果，新方法从三维层面提高溢油模型的物理真实感，为海上搜救模拟器生成溢油事故的三维场景、降低人力和财物成本、为海上溢油应急预案的制定提供直观的信息和依据。

关键词：潮汐建模；
溢油模拟；
溢油扩散与漂移；
溢油可视化

中图分类号：U491；
U698.7；
X52文献标志码：A文章编号：1672-0032（2024）02-0123-08

引用格式：张明，任鸿翔，朱天惠，等.渤海湾海上瞬时溢油三维模拟［J］.山东交通学院学报，2024，32（2）：123-130.

ZHANG Ming， REN Hongxiang， ZHU Tianhui， et al. Three-dimensional simulation on instantaneous oil spill in the Bohai Sea Bay［J］.Journal of Shandong Jiaotong University，2024，32（2）：123-130.

0 引言

海上溢油是船舶交通运输行业面临的严峻环境问题。随全球化贸易的增长，由船舶事故、恶劣天气条件、技术故障或非法排放等原因造成的海上溢油事故频发，在海洋中释放大量有害物质，包括石油化学品、重金属和有毒物质，对附近的海岸线和沿海生存环境造成恶劣影响，对旅游业、渔业和当地经济带来重大损失，给海洋生态带来严重破坏。溢油分为持续溢油和瞬时溢油2种，持续溢油量可能为几百t甚至几十万t，发生海上溢油事故时，在重力、摩擦力和风、波浪场、潮流作用下，溢油立即向四周扩散。海上溢油的善后处理难度较大，周期较长，成本较高[1]。因此，有必要开发溢油三维模拟系统模拟特定的海上溢油事故，重点研究海水动力潮汐模型、海上溢油数学模型及溢油可视化方法。

采用流体动力模型进行潮汐预报的结果准确、快速，是该领域研究的主要方法之一。Guo等[2]采用普林斯顿海洋模型和近岸模拟波浪模型进行混合处理，模拟大连地区的潮汐。Badri 等[3]采用波浪和人工风场理论改进溢油模型，并应用于霍尔木兹海峡的溢油事故。Barros等[4]提出数值浅水流模型的干湿算法，对复杂拓扑区域的计算精度较好。Fay[5]忽略次要作用项，在不同膨胀阶段保留主要作用项，得到平静海面油膜膨胀的半理论计算公式。Lehr等[6]在Fay理论的基础上，考虑风的影响，提出油膜椭圆模型，正确反映油膜在风作用下的展开过程。Liu等[7]基于Fay的油膜膨胀理论建立胶州湾水动力模型，耦合环境场影响建立胶州湾溢油模型，研究风、浪、流各要素对溢油运动的影响。Li等[8]研究动态溢油覆盖海面区域的电磁散射机制，采用油膜膨胀理论与溢油体积结合，建立动态溢油覆盖海面区域的三维几何模型。Erdoan[9]采用尼霍尔方程和高斯随机游走技术，以油膜扩散长度步进粒子，通过维诺（Voronoi）图展示油面污染情况。Liu等[10]提出基于散列函数和空间四叉树算法的耦合技术，将溢油模型与海洋模型结合。Toz等[11]采用PISCES 2轨迹模型模拟溢油的主要风化过程，据此确定伊兹米尔湾船舶溢油的风险区域。Liu等[12]提出基于动态贝叶斯网络的新型动态评估方法，将平流、蒸发、扩散等模型转化为贝叶斯网络模型，量化极端风况下的溢油风险。Fraga等[13]提出拉格朗日粒子模型，预测平静海面上石油扩散第一阶段的浮油直径。溢油模型有回归模型和对流扩散模型2种，回归模型是建立在试验结果与试验观测数据拟合基础上的模型，对流扩散模型主要根据对流扩散方程模拟溢油整体的轨迹。

大多数海上溢油的可视化研究聚焦于海上油膜运动的二维可视化，海上溢油的三维可视化研究较少。Ren等[14]采用平面折射技术将溢油轨迹的模拟结果转化为溢油发展变化的纹理图，在图形处理器（graphics processing unit，GPU）的片元着色器中通过投影纹理坐标获取纹理，模拟再现海上溢油三维场景。余枫[15]采用软件OpenSceneGraph的纹理烘焙和平面折射技术，实现海上溢油模型的实时仿真和三维可视化。魏国[16]在文献[15]的基础上，对海上溢油应急工作流程进行三维可视化处理，实现溢油围油栏的三维可视化，但未解决物理真实性提升的问题。邹长军[17]基于Navier-Stock方程实现基于物理的流体仿真，并提出窄带流体隐式粒子方法解决流体与移动障碍物的交互问题，通过 GPU并行加速实现。李海江[18]采用速度无散的光滑粒子流体动力学（divergence-free smoothed particle hydrodynamics， DFSPH）方法提出改进的粒子数密度模型，基于等静态密度假设，可统一处理包括非均匀采样的边界粒子在内的所有领域粒子，解决流体交界面附近密度不连续的问题，有效模拟海上溢油扩散、漂移及破碎等典型场景。田睿等[19]提出基于正弦平移法线纹理的溢油可视化方法，实现基于二维网格的溢油交互可视化，提高油膜与海浪融合的仿真效率，但物理真实感不强。

本文研究海水动力潮汐模型、海上溢油数学模型和溢油可视化技术，采用Mike 21 Flow模型计算渤海潮汐流的方向和速度，将回归模型和对流扩散模型结合构建海上溢油混合模型，采用泊松盘采样算法创建油膜粒子，通过Marching Cube算法重建等值面，模拟实现海上溢油的三维可视化效果，以期提高海上溢油事故的处理效率。

1 潮汐流的数值模拟

发生溢油事故后须尽快确定油膜的位置和运动方向以便尽快收集溢油，阻止溢油扩散。预测风速、风向、潮汐速度和潮流方向有助于提高海上溢油应急与搜救时溢油模型的有效性。随计算机和信息系统技术的快速发展，潮汐的计算和预测已越来越精确。采用Mike Zero水动力模拟计算程序中的MIKE 21 Flow Model二维自由表面流模型系统，对潮汐水力学、风和波浪产生的潮流、风暴潮通过水动力模块建模。

通过MIKE 21 Flow Model中质量守恒和垂直方向上的动量积分表示流量和水位的变化，流体系统整体的质量守恒方程为：

-ξ-t+-p-x+-q-y=-d-t，

式中：ξ为表面高程，m；
p、q分别为x、y轴方向通量密度，为单位时间通过单位面积流体的体积，与对应方向的深度平均速度有关，m3/（s·m2）；
d为时变水深，m；
x、y为空间位置坐标，m；
t为时间，s。

x方向的动量方程为

-p-t+--xp2h+--ypqh+gh-ξ-x+gp p2+q2C2h2－1ρw--xhτxx+--yhτxy－Ωq－fvx+hρw-pa-x=0，

式中：h为水深，h=ξ－d，m；
g为重力加速度，m/s2；
C为Chezy数，m1/2/s；
ρw为水的密度，kg/m3；
τxx、τxy为有效剪应力在不同方向的分量；
Ω为科里奥利参数，与纬度有关，s-1；
f为风摩擦因数；
vx、vy分别为风速v在x、y方向的分量，m/s；
pa为大气压，Pa。

y方向的动量方程为

-q-t+--yq2h+--xpqh+gh-ξ-y+gq p2+q2C2h2－1ρw--y（hτyy）+--x（hτxy）+Ωp－fvy+hρw--y（pa）=0，

式中：τyy为有效剪应力的分量；

u、v 为x、y方向的深度平均速度。

Chezy数一般根据曼宁数M经验公式计算得到，计算公式为C=Mh1/6。如果可获得波浪引起的床阻力，根据实际情况计算得到C=u g/ufc ，其中u为波浪边界层速度，ufc为波浪边界层摩擦速度。

模拟渤海潮汐的第1步是创建大小适中的三角形网格，网格太小易忽略边界条件对研究区域的影响，网格太大增加模型的计算时间。边界条件须与观测或输入的计算数据适应，在研究区域内计算网格须足够详细地表达海区信箱，动态因素对地形的影响较大。计算网格如图1所示。计算网格为非构造网格，网格最大面积约为4 km2，黄海区域约为40 km2，仿真区域生成39 856个节点和68 623个元素。研究区域划分的网格与不同地理位置深浅不一的水深区域相结合，包含计算网格、水位、边界等信息，叠加水深数据的渤海计算网格水深图，如图2所示。

采用本文提出的模型评估2011年6月“蓬莱19-3号”钻井平台溢油事故中平台所处位置的水位、流速和流向，并与观测站的观测数据对比验证，如图3所示。由图3可知，“蓬莱19-3号”钻井平台溢油位置的海面高程、流速和流向的模拟数据和观测数据结果比较接近。

a） 海面高程

b） 流速

c） 流向

2 海上溢油运动数学模型

海上溢油会发生扩散、漂移、蒸发、溶解、乳化、溢油融合、油和障碍物相互作用等过程，通常在溢油较长时间后才发生油膜蒸发和乳化过程。本文主要考虑瞬时溢油后较短时间内以扩散和漂移过程为主的溢油建模与仿真。

在溢油的早期阶段，扩散是影响油膜轨迹的重要过程。海上溢油事故往往发生在大风浪的恶劣天气条件下，因此，海上溢油应考虑风的作用。在模拟油膜扩散时，将油膜形状设置为椭圆，反映油膜在风作用下的拉伸现象，以多边形模拟油膜的椭圆边界，椭圆的半长轴、半短轴分别为a、b，油膜扩散的面积

A=πab。（1）

基于Lehr的试验数据推导经验公式[9]，式（1）可改写为：

{A}=2 270Δρ/ρ023{V}23{t}12+40Δρ/ρ013{V}13{v}43{t} ， {Amax}=105{v}43，（2）

式中：{A}为以m2为单位的A的数值，Δρ为油水密度差，ρ0为溢油的密度，{V}为以barrel为单位的溢油体积的数值，{v}为以knot为单位的海平面风速的数值，{t}为以s为单位的时间的数值，{Amax}为以m2为单位的油品可扩散的最大面积Amax的数值。

溢油扩散到一定阶段时，以油膜的漂移为主要运动形式。风、浪和表面流是影响溢油漂移速度和预测油膜轨迹的因素，表面流对油膜影响较大，大油膜被分解为若干细小油膜。可将油膜整体分解为若干油粒子，通过计算每个油粒子的运动预测油膜的运动。单个油粒子的漂移速度

vsingle=dXi/dt=vdrift+vdiff ，

式中：Xi为第i个油粒子在世界空间坐标系下的位置；
vdrift为漂移速度，vdrift=awv10+acvc，其中，v10为海面上方10 m处的风速，vc为流速，aw为风速影响因子，aw=0.03，ac为流速影响因子，ac= 1.0；
vdiff为扩散速度，vdiff=vdriftRnexp（iαn），其中，Rn为随机数，Rn∈[-1，1]，αn为随机角度，αn∈[0，π]。

3 海上溢油可视化关键技术

采用泊松盘算法对油膜区域的油粒子进行随机取样。泊松盘取样是将点分布在一定空间的算法，2点间的距离不小于泊松盘半径。将油粒子抽象为无数个直径不同的球形粒子，对椭圆油膜区域的贡献随机。将椭圆的原点设为O（x， y， z），油膜平面为Oxz，海平面轴的高度为y。假定溢油初始阶段y=0，确定油粒子在椭圆油膜内的方法为：

（x-h）2/a2+（y-k）2/b2≤1，

式中（h， k）为椭圆中心位置坐标。

油膜粒子数是影响海上溢油可视化效果的重要因素。生成的粒子越少，对海上溢油精细尺度的特征模拟越不准确；
粒子越多，溢油可视化模型越灵活逼真。但粒子过多易导致计算成本较高，计算速度大幅减小，须根据计算设备的配置确定可视化溢油模型中的粒子总数。

通过泊松盘采样油膜粒子后，对基于粒子的流体模拟结果进行可视化研究。通过Marching Cube（MC）算法构建多边形重构流体表面，渲染多边形。先定义平滑核密度函数，将模拟区域划分为网格，计算各网格点周围粒子密度的总和，采用MC算法提取等值面。

通过7个步骤实现流体表面的可视化。

1）将三维空间划分为空间网格。在MC算法的初始化阶段，粒子较多，需在整个区域中划分网格。流体域是1组网格单元，在每步计算中遍历区域内的所有粒子，找到网格单元与立方体的表面交点。

2）计算每个空间网格顶点的密度。采用划分的网格作为油膜场景的采样点，得到油膜密度场数据。

3）定义并计算油膜表面的密度等值面的阈值。

4）比较每个顶点密度与等值面的阈值。等值面的阈值不大于顶点密度时，设顶点为1，顶点在等值面的内部；
反之，设顶点为0，顶点在等值面的外部。根据比较结果得到空间网格的结构。

5）计算空间网格密度等值面的顶点坐标。三角形网格的1条边将与等值面相交，顶点在等值面的内部时，另1个顶点在等值面的外部。确定与等值面相交的线段后，

通过线性插值方法定义网格边与等值面间的切割点。

6）计算各顶点的单位法线，对法线内插，输出三角形顶点和顶点法线。

7）根据顶点位置、顶点法线、照明类型等输入数据，采用软件Unity 3D着色器将流体域转化并可视化到屏幕上。

基于MC算法的油膜表面可视化结果如图4所示。黑色为油膜，褐色为海水，海水上面为大气。

4 海上溢油的模拟实现

海上溢油模型可视化实现过程如图5所示。

先建立溢油的初始数据，如溢油类型、释放量、油水密度等。通过Mike Zero和Mike 21 FM计算得到潮汐数据，水动力潮汐模型的计算输出数据可作为溢油扩散和漂移模型的输入数据计算油膜的演变与运动，达到最大溢油扩散面积时，采用泊松盘算法将油膜渲染成油粒子，再计算浮油漂移的轨迹和特征，采用MC算法重建等值面，实现海上油膜的可视化。

采用Unity 3D实现海上溢油模拟，通过C#语言编程，模拟渤海湾某海域溢油事件，模拟信息为：溢油事件发生日期为2022-07-04，溢油位置为北纬38.4°、东经120.1°附近，溢油体积为1 000 m3，ρ0 =806 kg/m3，ρw = 1 025 kg/m3，v10=0.5 m/s。

在溢油模型的开始扩散阶段，浮油的扩散形式为椭圆形增大面积，根据式（1）（2）计算椭圆的长、短轴。溢油扩散模拟结果如图6所示。达到最大油膜扩散面积时，采用泊松盘算法将油膜分成粒子。修改该算法，采用优化飞镖算法生成直径不同的油粒子，满足溢油模拟系统的要求。

本文模拟实现2种海上溢油场景。第1种是在广阔海域无障碍物时，油膜在海上自由扩散与漂移，溢油可视化结果如图7所示。在模型中，泊松圆盘采样算法产生约2 500个油粒子，通过添加粒子的浮力、颜色和其他属性可连续观察油膜的发展。溢油的初始扩散阶段，油膜足够厚，在波浪、风和表面张力的影响下，油膜不破裂。溢油不断扩散到足够薄时，波浪和风等因素抵消溢油的黏度影响，油膜逐渐分解为许多小部分。第2种场景是受岛屿等障碍物影响下的油膜扩散和漂移，油膜在海上扩散与漂移到岛屿附近后，受岛屿的阻碍产生一定堆积，溢油可视化结果如图8所示。

海上溢油后，大部分原油迅速扩散为1层薄膜，在更广阔的海洋范围内破碎成小的浮油面。受物理海洋学中朗缪尔环流的作用，溢油薄膜随时间的推移逐渐撕裂合并成条状溢油风积丘。

将可视化溢油模型模拟结果与2010年7月“深水地平线”钻井平台溢油事故现场图片对比，如图9所示。由图9可知：溢油事故现场的原油形成沿风向的溢油风积丘，在溢油模型中，将油膜抽象为油粒子，油粒子逐渐向图中红线所示的条状风积丘聚集，与真实溢油现场相似。本文的三维可视化成果基于计算机图形学领域相关研究基础，对成果的判断依据为与真实情况的接近程度。

5 结束语

本文采用Mike 21 Flow模型对渤海地区潮流进行数值模拟，将溢油运动过程分为油膜的扩散、漂移两阶段进行计算，建立基于风和潮流条件的溢油运动数学模型，实现溢油轨迹的实时模拟和三维可视化。此模型简单高效，满足应急溢油响应的基本需求，比传统纹理贴图的方法极大地提升可视化模拟的物理真实感。采用泊松盘采样方法动态实现油粒子的融合和分离，比基于网格的物理方法的真实性更好。

参考文献：

[1] LINDSEY B. Biodegradation andbioremediation[M].[S.l.]：
Scientific e-Resource， 2018.

[2] GUO W J， WANGY X. A numerical oil spill model based on a hybrid method[J].Marine Pollution Bulletin， 2009，58（5）：726-734.

[3] BADRI M A， AZIMIAN A R. Oil spill model based on the kelvin wave theory and artificial wind field for the Persian Gulf[J].Indian Journal of Geo-Marine Sciences， 2010，39（2）：
165-181.

[4] BARROS M L C， ROSMAN P， TELLES J C F. An effective wetting and drying algorithm for numerical shallow water flow models[J].Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2015，37（3）：
803-819.

[5] FAY J A. The spread of oil slicks on a calm sea[M].Boston，USA：
Springer， 1969：
53-63.

[6] LEHR W J， FRAGA R J， BELEN M S， et al. A new technique to estimate initial spill size using a modified Fay-type spreading formula[J].Marine pollution Bulletin， 1984，15（9）：326-329.

[7] LIU Z K， CHEN Q， ZHANG Y W， et al. Research on transport and weathering of oil spills in Jiaozhou Bight， China[J].Regional Studies in Marine Science， 2022， 51：102197.

[8] LI D， ZHAO Z， MA W， et al. Investigation of electromagnetic scattering mechanisms from dynamic oil spill-covered sea surface[J].Remote sensing， 2023， 15（7）：1777.

[9] ERDOAN M. Oil spill surface fate model through random walk method with drifted particles in voronoi diagram[D].[S.l.]：
Middle East Technical University， 2021.

[10] LIU D， MU L， HA S， et al. An innovation coupling technique for integrating oil spill prediction model with finite volume method-based ocean model[J].Marine Pollution Bulletin， 2022， 185：
114242.

[11] TOZ A C， KOSEOGLU B. Trajectory prediction of oil spill with Pisces 2 around Bay of Izmir， Turkey[J].MarinePollution Bulletin， 2018， 126：
215-227.

[12] LIU Z， HAN Z， CHEN Q， et al. Risk assessment of marine oil spills using dynamic Bayesian network analyses[J].Environmental Pollution， 2023， 317：
120716.

[13] FRAGA F， CARLOS A D. A Lagrangian analysis of the gravity-inertial oil spreading on the calm sea using the reflective oil-water interface treatment[J].Environmental Science and Pollution Research， 2021， 28（14）：
17170-17180.

[14] REN H X， ZHANG X， YIN Y. 3D real-time rendering of continuous oil spill on simulated ocean[J].Journal of System Simulation， 2008， 20（19）：
5117-5120.

[15] 余枫.海上搜救模拟器中近海溢油的实时仿真与可视化[D].大连：大连海事大学， 2010.

[16] 魏国.关于海上溢油应急处理的仿真和三维可视化的研究[J].电子技术与软件工程，2016，79（5）：100.

[17] 邹长军.航海模拟中船舶溢油应急仿真及可视化[D].大连：大连海事大学， 2018.

[18] 李海江.航海模拟器中基于物理模型的海洋场景建模[D].大连：大连海事大学， 2020.

[19] 田睿，刘秀文.海面溢油油膜交互可视化技术研究[J].中国航海，2022，45（1）：107-111.

Three-dimensional simulation on instantaneous oil spill in

the Bohai Sea Bay

ZHANG Ming1， REN Hongxiang2*， ZHU Tianhui2， ZHOU Yi3， TANG Haina4

1. Maritime College， Fujian Chuanzheng Communications College， Fuzhou 350007， China;

2. Navigation College， Dalian Maritime University， Dalian 116026， China;

3. CNOOC Energy Technology and Service-Oil Production Services Co.， Ltd.， Tianjin 300452， China;

4. School of Artificial Intelligence， University of Chinese Academy of Sciences ， Beijing 100049， China

Abstract：To simulate and predict the dynamic characteristics of three-dimensional offshore oil spills in the Bohai Sea Bay in real-time， the Mike 21 Flow model is used to calculate the hydrodynamic tides in the Bohai Sea. A mathematical model of offshore oil spill movement is established to exhibit the diffusion and drift processes of instantaneous oil spills. The Poisson disk sampling algorithm is used to create oil particles of different diameters， and the Marching Cube algorithm is employed to reconstruct iso-surfaces， simulating the offshore oil spill process and comparing it with data from real cases. The results of the oil spill simulation show that the sea surface elevation， current velocity， and direction at the location of the oil spill accident vessel assessed by the tidal flow numerical model are close to the data from observation stations. The trajectory of crude oil movement from the mathematical model of offshore oil spill movement is similar to the actual trajectory. Compared to traditional oil spill simulation visualization limited to two-dimensional non-interactive effects， the new method improves the physical realism of the oil spill model from a three-dimensional perspective. This method helps in generating three-dimensional scenarios of oil spill accidents for offshore search and rescue simulations， reducing manpower and material costs， and providing intuitive information and basis for the formulation of offshore oil spill emergency response plans.

Keywords：tidal modelling; oil spill simulation; the diffusion and drift processe of oil spill; oil spill visualization

（责任编辑：王惠）

收稿日期：2023-03-24

基金项目：国家自然科学基金资助项目（52071312）；
交通运输行业重点科技项目（2022-ZD3-035）；
辽宁省应用基础研究计划；
大连市科技创新基金项目（2021JJ12GX031）

第一作者简介：张明（1971—），男，福建南平人，副教授，工学硕士，主要研究方向为海运技术与海运安全，E-mail：zhangming_1971@126.com。

*通信作者简介：任鸿翔（1974—），男，黑龙江绥化人，教授，工学博士，主要研究方向为虚拟现实与视景仿真技术，E-mail：dmu_rhx@163.com。

DOI：10.3969/j.issn.1672-0032.2024.02.017