复杂海湾水质预测中油粒子模型的寻址优化法

周晨琦，季智灵，孔 俊

(河海大学水文水资源与水利工程科学国家重点实验室，南京 210098)

随着社会经济的迅速发展，石油作为重要的发展能源，其需求与日俱增，20世纪60年代以来，海上石油的勘探与开采不断发展，但突飞猛进的海洋石油运输业与石油开采业，也导致海洋溢油事故频频发生，这不仅对人们造成了巨大的经济损失，也对附近的海洋生态环境造成了严重的破坏，发生在近海的溢油事故甚至会危害附近居民的健康[1]。鉴于此，国内外的学者们对溢油事故的成因与海上溢油的过程进行了一系列的研究，以完善溢油应急[2]体系。

溢油在进入海洋之后，会经过一系列复杂的物理化学变化[3]，主要包括三大类。①扩展过程：油膜在海面上面积扩大；②输移过程：溢油在海洋动力作用下的迁移运动，包括水平漂移和垂直掺混等；③风化过程：蒸发、溶解、 乳化、光氧化和生物降解等。

随着研究的深入，人们发现数学模型能够有效模拟、预测溢油的迁徙规律，为后续的科学决策提供理论依据。目前，根据不同的数学模型计算理论，溢油运用的数值模拟方法可分为油膜扩展模型、对流扩展模型和油粒子模型三类[4]。

在有关溢油扩展的研究中，文献[5]只考虑重力和溢油体积，建立了扩展直径公式。Fay[6]全面考虑了重力、表面张力、惯性力和黏性力的作用，认为油膜始终保持圆形，并提出了三阶段油膜扩展理论，在溢油扩展模型方面取得了阶段性的突破。后续众多学者[7-8]对Fay模型进行了修正，提出了众多扩展模型。对流扩散方程在溢油运动轨迹的研究中也得到了应用[9]，但是对流扩散模型的稳定性不高，容易失真。王琰等[10]率先提出了油粒子的概念，建立了油粒子模型，将油膜视为众多油粒子的组合，其中每一个粒子表征一定数量的示踪物质，这些模型粒子的平流过程具有拉格朗日性质，可用拉格朗日方法模拟，而剪切流和湍流引起的粒子扩散属于随机运动，可以用随机走动法进行模拟。油粒子模型创新性地使用粒子的宏观运动计算溢油的运动轨迹，目前在中外得到了广泛的应用。但模型中大量的油粒子可能会使模拟时间过长，对此，提供了一种新的油粒子定位方式，以实现快速定位。

1 理论依据

1.1 溢油过程

当海上石油泄漏时，油膜在重力、黏性力、惯性力以及表面张力的作用下在水平方向不断扩张，这种自身扩展是溢油初期油膜运动最主要的形式，在整个扩展过程中，油膜厚度不断减小，最后在水流、风场因素等作用下发生破碎，此时，油膜主要进行紊动扩散。使用油粒子模型，能够准确模拟出溢油在第二阶段的运动情况，但也忽略了初期油膜的自身扩展，为了更好地还原整个溢油扩散过程，应该对具体的溢油事故进行分析，对于大规模瞬时泄漏，油膜具有一个快速扩展的过程，此时若忽略自身扩展，直接采用油粒子模型对紊动扩散过程进行模拟，必然导致溢油模拟产生较大误差，显然是不合理的，而对于一些小规模的溢油事故，油膜扩展面积小、时间短，可以考虑忽略此过程，直接进行紊动扩散的模拟，不会对最终结果造成太大的影响。

1.1.1 油膜扩展过程

Fay模型[6]考虑了重力、表面张力、惯性力和黏性力的共同作用对油膜扩展的影响，但其理论建立在静水假定基础上，认为油膜始终保持圆形，而在实际的海上溢油事故中，油膜扩展呈现明显的各向异性，因此采用修正的Fay模型对扩展直径进行计算。

(1)

式(1)中：r为短轴方向的扩展直径；c1为经验常数；Δg为约化重力加速度，Δg=(1-ρ0/ρw)g，其中ρ0为油膜密度，ρw为海水密度；V为油膜体积；t为时间；l为长轴方向上的扩展直径；Uw为海面10 m处风速；c、δ、ε均为随油的种类及性质变化的经验常数。

油膜自身扩展的持续时间tf计算公式为

(2)

式(2)中：c1、c2均为经验系数；ν为运动黏性系数。

油膜在扩展的同时，由于受到风场和流场的作用，会产生漂移，油膜漂移通常以等效圆油膜中心处的位移来判断，假设油膜中心初始时间t0位置为S0，经过Δt时间后，油膜的位置S计算式为

(3)

式(3)中:VL为拉格朗日速度。

油膜中心的漂移速度和方向由表面水流和风引起的流速矢量求和所得，表达式为

v0=Uf+CwUw

(4)

式(4)中：v0为油膜中心漂移速度；Uf为表面水流速；Cw为风拽力系数，一般取0.03～0.04；Uw为海面10 m处风速。

1.1.2 紊动扩散

油膜经过一段时间的自身扩展运动后，油膜变薄、破碎，此时油膜运动以紊动扩散为主，将经过自身扩展阶段的油膜转化为一系列的粒子，就可以采用油粒子方法经行模拟，油粒子方法将油粒子定义为很小的圆球，直径为10～1 000 μm，由于每个油粒子都很微小，要做到精确地表示一个油膜需要大量的实际粒子数，而在计算机中同步堆栈太多的粒子特性参数是不可能的。因此，在实际应用中，最大可能的粒子数目需要通过计算机的容量、运行时间等因素进行确定，采用附加体积参数的方法来实现对油粒子的模拟，每个油粒子代表溢油体积的一部分。将某个油粒子体积参数定义为Vi，其所占的油膜总体积的百分比为fi，则每个油粒子的特征体积Vi=fiV0，其中V0表示初始油膜的体积。

当油膜被分割为N个油粒子后，以油膜质心的位置为直角坐标系坐标原点，可以得到每个油粒子的初始坐标(xj,yj,j=1,2,…,N)，随后可以进行单个油粒子运动的计算。油粒子在海中的漂移过程主要分为两个部分，即对流过程和紊动扩散。这两个过程综合作用出油粒子在每个时间步长的具体分布。假定粒子ti时刻所处的位置为(xi,yi)，经过Δt时间的运动后，粒子在ti+1时刻新的位置坐标(xi+1,yi+1)可以表示为

(5)

式(5)中：ΔxCi、ΔyCi分别表示油粒子由于对流运动在x、y方向上产生的位移；ΔxDi、ΔyDi分别表示油粒子由于湍流弥散在x、y方向上产生的位移。

其中对流运动产生的x、y方向上的位移分量(ΔxCi,ΔyCi)可表示为

(6)

式(6)中：ui、vi分别表示粒子在ti时刻沿x、y方向的流速，Δt=ti+1-ti。

油粒子的对流位移由水流力和风曳力引起，因此，油粒子的漂移速度ui、vi计算公式为

(7)

式(7)中：Uwxi、Uwyi分别表示ti时刻水面以上10 m的风速在x、y方向上的分量；Ufxi、Ufyi分别表示ti时刻水面流速在x、y方向上的分量。

湍流弥散具有随机性, 而海上的溢油扩散过程实际上正是湍流的一个弥散过程，因此紊动扩散即油粒子由水流的随机性脉动导致的空间位移，可采用抽取随机数的方法来计算得到一个时间步长内粒子的可能扩散距离。由于水体紊动扩散引起的随机扩散位移在x、y方向上的位移分量(ΔxDi,ΔyDi)可表示为

(8)

式(8)中：Rx、Ry为分布在(-1，1)区间内的标准正态分布随机数；Kx、Ky分别表示x、y方向上的紊动扩散系数。

1.2 风化过程及模型

溢油在海面经历扩展、漂移等物理输移过程的同时，也经历着蒸发、溶解及乳化等风化过程，风化过程导致溢油的溢油量、油粒子的组成发生变化，不影响油粒子的水平位置。初期的油粒子模型忽略了溢油的风化过程，对于风化造成的溢油量损失等情况无法准确模拟，经过中外众多学者的相关研究，溢油的风化过程被逐步完善。

1.2.1 蒸发

溢油的蒸发过程较复杂，油粒子组成、气温水温、溢油面积、油膜厚度、风速和太阳辐射等因素均影响油膜蒸发。其蒸发速率首先取决于油膜的化学成分，其次包括溢油的物化属相、环境分度及风的作用等。

采用SA模型假定：①油膜内部扩散无限制(适用于气温高于0 ℃以及油膜厚度小于10 m)；②油膜完全混合均匀；③油膜组分在大气中的分压与蒸汽压比忽略不计。如此假定后，油膜中某一组分的蒸发速率可以表示为

(9)

(10)

式(10)中：k为蒸发系数；Aoil为油膜面积；Sci为组分油i的蒸汽Schmidts数。

1.2.2 乳化

乳化是指海上溢油在风化过程中油与海水混合形成油水乳化物的过程，因此乳化过程可看作两个部分，分别是溢油向水体扩散后形成水包油乳化物过程和油滴吸收水分以形成油包水乳化物过程。

(1)形成水包油乳化物过程。溢油向水体中的运动机理包括溶解、垂向扩散、沉淀等。垂向扩散是溢油发生后最初的主要过程，水流的紊动能量使油膜破碎成油滴，形成水包油的乳化物。这些乳化物被表面活性剂稳定，阻止油滴返回到油膜。恶劣的天气状况下，油膜的波浪破碎是最主要的扩散作用力，在平静的天气的状况下伸展压缩运动则是最主要的扩散作用力。油膜扩散到水体中形成乳化物的油分损失量计算公式为

D=DaDb

(11)

式(11)中：Da为油膜进入水体的分量；Db为进入水体后没有返回油膜的分量，计算公式为

(12)

式(12)中：μoil为油膜的黏度；γow为油-水界面张力；hs为油膜厚度。

另外，油滴返回油膜的速率为

(13)

式(13)中：Voil表示油膜体积。

(2)形成油包水乳化物过程。可用平衡方程[式(14)]表示油中的含水率yw变化为

(14)

式(14)中：R1为水的吸收速率；R2为水的释出速率，表示为

(15)

1.2.3 溶解

溢油的溶解度较小，在全部溢油过程中溢油的溶解量一般不超过总量的1%，占主要溶解量的是低碳的轻组分油。溢油某一组分i的溶解率Vdsi用式(16)表示：

(16)

实际海洋预警预报活动中，往往要求能在最短时间内(1～2个潮周期)对油膜的运动范围做出判断，在如此短的时段内，乳化、蒸发等化学过程减少的表面原油仅占总量的1%～5%，因此在开展短期快速预警过程中，偏于保守考虑，暂不考虑油膜的化学消耗作用。

2 油粒子模型优化

2.1 优化寻址法

针对非结构的三角形网格，水位、流速、流向等信息可由潮流数值模型预测后布置于三角形的节点或单元中心位置，空间布局没有规律。基于此动力输出文件，构建油粒子模型时，不可避免地存在寻址难度，跟踪某一个油粒子上一时刻和当前时刻的位置，都需要在众多的三角形中去一一判别，可采用的方法是面积法，即

S1+S2+S3=S总

(17)

式(17)中：S1、S2和S3分别表示油粒子与某个三角形网格3个网格节点组成的三个三角形的面积;S总表示该三角形网格的面积，当等式成立时，即代表油粒子处于当前三角形网格。如图1所示，该方法尽管能够准确判别出油粒子的位置，但较为耗时，若模型中油粒子数量数以万计，则在每次模拟过程中寻址等待的时间就非常漫长。

针对这一不足，提出了栅格法技术，即在真实的非结构网格框架上，另构建一套虚拟的规则网格系统，该套网格的精度根据需要在横向和纵向上可疏密优化，如图2所示，虚线网格为虚拟栅格格网，实线为真实的水动力模型计算网格，前者的密度可以高于后者。虚拟的规则网坐标系统有助于快速判别油粒子的位置，即油粒子的空间坐标位置为

(18)

式(18)中：x、y代表油粒子在以油膜初始质点为坐标原点的直角坐标系上的坐标；dx、dy分别代表虚拟网格单元在x轴与y轴上的长度。

图1传统方法油粒子判别示意Fig.1 Conventional oil particle discrimination scheme

图2 真实计算网格和虚拟栅格Fig.2 Real compute grid and virtual grid

当上述位置确定后，则根据规则网格四个节点的流速信息，采用双线性插值法，可获得当前点位的速度和方向，为下一次油粒子运动提供动力信息，表达式为

F=F1+(F2-F1)y+(F4-F1)x+(F1-F2+F3-F4)xy

(19)

式(19)中：F1、F2、F3、F4分别为网格点上的已知流速；x、y表示与网格节点的距离，流速内插方法示意图如图3所示。

图3 流速内插法示意图Fig.3 Schematic diagram of flow rate interpolation

值得注意的是，虚拟网格系统下，矩形网格每个节点的流速信息也是通过插值获得，由真实非结构网格节点或中心点的信息提供。矩形单元每个节点的空间位置，决定了其所处的三角形网格信息，这就需要构建虚拟网格和真实网格的拓扑关系，可在溢油模型运行前，实现构建，而不影响油粒子模型的运行效率。

2.2 复杂结构海湾处的改良方法

通过虚拟栅格法，可以快速判别任一油粒子的空间位置，但也存在不足。如在海湾内存在防波堤、复杂结构建筑物时，采用传统的真实网格，可以在模型中设置通过固壁边界的法向流速为0。而采用虚拟网格构建时，很容易将这类占海用地区域也视为可通行水域，导致在模拟过程中出现油粒子穿越建筑物的现象，为此，提出了寻址优化方法，实现油粒子的绕行。

如图4所示，A点为n时刻油粒子位置，C点为n+1时刻油粒子位置，B、D分别为建筑物(如防波堤)的两个头部端点，在计算过程中当油粒子穿越结构物时，满足以下判别指标：

SABC+SACD=SABD+SBCD

(20)

式(20)中：SABC、SACD、SABD、SBCD分别代表四个三角形的面积。

图4 特殊建筑物阻隔油粒子运动处理示意图Fig.4 Schematic diagram of motion processing of special building blocking oil particles

此时，在程序中，通过插值方法，将油粒子修正至建筑物与油粒子路径的交点处。

3 改进后的油粒子模型应用分析

3.1 案例构造

为了验证新方法下油粒子的运动，这里以一个具有复杂防波堤形式的海湾为例。在此案例中，防波堤在港区外侧共布置了4条，以满足消浪作用。外海潮位进入港区后，沿防波堤之间的水道往复运动，涨、落急时刻流场如图5所示。

图5 港区的防波堤布置及流场Fig.5 Breakwater layout and flow field of the port area

3.2 模拟结果分析

油粒子运动的模型预测结果如图6所示，溢油后每隔3 h给出油膜的整体分布图，从图6中可以看出，溢油点位于海湾北岸处，受溢油发生后前6 h的涨潮流影响，油膜绕过4条防波堤堤坝，进入港区，至3 h后油膜完全进入港区，并沿海岸发生漂移，经过6 h后主体部分运动到南岸一侧。此后，油膜跟随落潮流流出港区，从9～12 h之后的油膜分布图中，可以看出油膜在设有多个防波堤的复杂地形处，均沿着防波堤之间的通航水道运动，并且油粒子的运动不会出现穿越防波堤的迹象，比较真实地反映了油膜的漂移过程。

图6 防波堤港区溢油后油膜运动过程Fig.6 Process of oil film movement after oil spill in breakwater port area

4 结论

针对常规油粒子模型在网格中的定位法对于大数量级油粒子模型计算时间较长的缺陷，提出了一种适用于复杂海湾地形的优化寻址法，得出以下结论。

(1)虚拟栅格法可以有效提高油粒子定位的速度，目前油粒子在网格中的定位多采用面积法判别每个油粒子在实际非规则网格中的具体位置，这一步骤计算烦琐，会大大提高模型的运算时间，虚拟栅格法构建了一套虚拟的规则网格系统，其横纵精度可自行调整，不会受到复杂的实际网格的约束，可以有效实现油粒子位置的快速判别。

(2)通过增设一个判别指标，优化寻址，可以使本文方法适用于复杂的海湾地形。由于虚拟栅格法并非传统的真实网格，不能设置固壁边界条件，因此通过一个判别指标限制油粒子的运动，避免油粒子穿越结构物。根据本文复杂海湾地形的实际案例计算结果，可以看出采用新计算方法的油粒子模型，既大大提高了模型的运算效率，同时又具有很强的适应性，可以应用于不同的地形条件。