海底管道泄漏原油扩散漂移规律数值模拟*

曹学文，温家铭，孙 媛

(中国石油大学(华东) 储运与建筑工程学院，山东 青岛 266580)

0 引言

海底输油管道的老化以及冲蚀和腐蚀所造成的冲刷穿孔失效是海管破裂泄漏的重要原因[1-2]。研究水下原油管道泄漏扩散漂移过程，对于相关部门作出应急决策和减轻事故风险具有重要意义[3]。

国外学者已研究发展了大量的水下溢油输运扩散模型，其中，Mcdugall[4]，Fannelop等[5]和Milgram[6]考虑泄漏垂直浮射流的气体膨胀作用来描述水下油井井喷过程，但未考虑与周围水体的相互作用；Yapa等[7]考虑了射流卷吸和原油溶解等影响因素，提出了基于拉格朗日积分法的水下原油泄漏喷射三维数学模型；Johansen[8]和Zheng等[9]又提出了DeepBlow模型和CDOG模型，能够描述深水油气泄漏的分离过程和溢油的运移扩散过程。

我国对于海底管道泄漏溢油的研究起步较晚。朱红钧等[10]、卢盛灿等[11]、李志刚等[12]采用VOF模型对水下原油泄漏进行了二维非稳态数值模拟，研究了泄漏方向、溢油速度和环境水深等因素对泄漏轨迹形态和输移扩散行为的影响；刘瑞凯等[13]模拟了海底埋地输油管道泄漏原油自海泥至海面中的过程；李新宏等[14]模拟了深水水下分离器不同泄漏位置的油气泄漏过程。

管道在浅海环境下泄漏上浮的时间较少，救援时间更紧张，以上研究未考虑波浪与洋流结合条件下浅海管道原油泄漏扩散漂移过程。本文考虑洋流波浪环境，探讨不同工况下溢油的扩散漂移规律，评估原油溢出海面的时长、海面溢出位置和漂移扩散范围等参数，为指导浅海海底管道溢油事故的处理提供理论参考。

1 水下原油泄漏数学模型

本文采用二维数值模拟，X坐标代表水平方向，Y坐标代表竖直方向。VOF模型是欧拉多相流模型的一种，基于混合相动量方程和流体相函数方程，若计算区域单元内的某流体体积分数取为0，则该单元内没有该流体，反之体积分数取为1，代表只有该流体。下标取值1，2，3分别代表气、水、油三相，二维输送方程和动量方程分别为：

(1)

(2)

式中：α为体积分数；V为单元平均流速，m/s；t为时间项，s；ρ为密度，kg/m3；p为压力项，Pa；μ为动力粘度，Pa·s；F为外部作用力，N。

为预测溢油在海水环境下非稳态湍流运动，湍流模型采用时均雷诺的带旋流修正κ-ε模型[15]。

小振幅波是一种简谐起伏振动的推进波，是研究复杂海面波动问题的基础[16]，故本文采用二维小振幅推进波，即艾利线性波来描述海面的波浪环境。

波面方程η和速度势方程φ分别为：

η=asin(kx-ωt)

(3)

(4)

式中：η为波面波动高度，m；a为波动振幅，m；D为水深，m；g为重力加速度，取值为9.81 m/s2；k为波动波数，rad/m；ω为海面波动圆频率，rad/s。

2 水下原油泄漏数值建模

图1为水下泄漏原油计算区域范围和计算边界条件的示意图，整个计算区域为长120 m、高14 m的二维矩形，水深D为10 m，水面上方为外界大气，设置波长L为20 m，波高H为1 m，振幅为0.5 m。为了在有限空间距离内模拟无限外域的波浪情况，减轻波浪在右边界产生反射波对波浪造成的干扰，需在右侧出口邻近区域设置一个关于水面对称的矩形阻尼层消波区[17]，此处设置消波区长度为1个波长的距离(20 m)，高度为2倍波高的高度(2 m)，泄漏口在入口下游20 m处，取泄漏口尺寸0.01～0.1 m。

计算区域左侧边界为速度入口，则波浪速度为：

(5)

式中：uw，vw和uc分别为波浪水平速度、波浪竖直速度和洋流速度，m/s。此处洋流均一速度取为0.1 m/s，水面上方左侧入口速度取0 m/s。

图1 洋流波浪环境泄漏原油扩散漂移的计算区域示意Fig.1 Schematic diagram of calculation area for diffusion and drift of oil spills under the combined action of wave and current

计算区域右侧出口边界为压力出口，泄漏口边界为速度入口，管内外压差、管径以及管位置等会影响原油泄漏速度，取泄漏速度2～6 m/s，取值间隔1 m/s。在阻尼层消波区，带附加源项的动量方程为：

(6)

(7)

(8)

式中：u和v分别为水平速度和竖直速度，m/s；μr为消波区的阻尼系数；FL为模拟区域长度，m；Lr为消波段长度，m；β为经验约束系数，值取为8[18]。

波浪洋流入口速度边界、出口压力边界和消波区的阻尼源项，均通过C++语言编程的用户自定义宏函数(UDF)编译到CFD软件中。模拟得到计算区内的波浪洋流场，波浪轮廓的相分布如图2所示。

图2 计算区域内的相分布Fig.2 The phase distribution of the calculation area

计算区域水深小，不考虑海水密度分层，设置空气和海水的密度分别为1.225和1 025 kg/m3，取原油密度为750～970 kg/m3。设置空气、海水和原油的粘度分别为1.8×10-5，1×10-3和5×10-2Pa·s。

非稳态速度与压力耦合方法采用PISO算法，启用隐式体积力选项，部分平衡动量方程的压力梯度和体积力，提高重力场中原油扩散漂移模型的稳定性。当考虑波浪洋流的泄漏环境时，先模拟稳定的波浪洋流场，再将泄漏口与环境流场进行耦合计算。

本文考虑不同因素对泄漏原油扩散漂移行为的影响，表1列出了对应的18种计算工况，其中最小和最大泄漏量的工况分别对应案例15和案例18。

表1 模拟工况Table 1 Simulation conditions

表1(续)

3 模拟结果与讨论

3.1 海洋环境对泄漏原油扩散漂移范围的影响

对静水、洋流、波浪、波浪洋流结合的环境工况做了对比，图3和图4显示了在不同环境工况下原油从泄漏口到水面的扩散上升过程。

在静水环境下，泄漏原油上升到水面的时长为13 s，当泄漏原油上升高度为4 m以内时，其泄漏形态为连续垂直柱状油流；当泄漏原油上升垂直高度大于4 m时，柱状油流会分散成油滴或油团，且整个分散油滴关于泄漏中心垂线对称分布。对于洋流环境，泄漏原油受到水流剪切力、重力、浮力以及惯性力等作用力，泄漏原油初始2 s上升的高度基本与静水环境下一致，这是因为泄漏口处原油受惯性力的影响要大于洋流剪切力的影响，上升高度高于2 m后，洋流剪切作用影响开始显现，其对应的扩散高度随时间变化曲线的斜率变小，油滴整体会向下游扩散，导致原油到达水面时间更长，且到达水面时长比静水环境下多3 s。对于波浪环境，由于水粒子的振荡，原油泄漏上升有时向上游漂移，有时向下游漂移，但整体上看来仍然关于泄漏中心垂线呈对称分布，由于波浪水粒子的水平最大速度(1.2 m/s)远大于洋流速度，故波浪对原油的剪切作用要大于洋流对原油的剪切作用，即波浪环境下原油上升高度对时间变化曲线的斜率小于洋流环境下的曲线斜率，原油上升到水面的时长相对洋流环境会更长。对于波浪洋流环境，由于洋流作用和波浪水粒子振荡作用，分散油滴向下游扩散的水平距离大于前3种工况。

泄漏原油上升到水面后，会继续在水面漂移，图5和图6显示了不同环境下原油在水面的漂移过程。在静水环境下，泄漏原油会像喷泉一样向水面扩展，水平扩散距离随着时间的变化曲线基本呈现线性关系，漏油水面漂移速率为0.3 m/s；对于洋流环境，漏油水面漂移速率为0.45 m/s，该值近似为环境漏油水面漂移速率与洋流速度之和；对于波浪环境，漏油水面漂移速率达到了0.48 m/s；而对于波浪洋流环境，波浪表面的水粒子速度有利于漏油在水面的扩散，漏油在水下升高过程的最大水平扩散距离能够达到5.5 m，在水面的水平扩散率达到了0.58 m/s，波浪洋流环境下的溢油事故需及时响应处置。

图3 不同环境下原油的水下扩散过程Fig.3 Underwater diffusion process of oil under different environments

图4 不同环境下原油上升高度变化Fig.4 Variation of oil rising height under different environments

3.2 原油密度对泄漏原油扩散漂移范围的影响

图7为不同原油密度的漏油上升高度随时间的变化曲线。在波浪洋流条件下，原油上升到水面的时间随原油密度增大而增大。在海底以上6 m高度范围内，原油上升高度基本随时间呈线性变化；随着上升高度进一步增加，扩散高度曲线斜率会波动性地变小，尤其当原油将要扩散到水面时，扩散高度曲线斜率会明显减小，斜率减小的幅度会随着原油密度的增大而增大。

图5 不同环境下原油的水面漂移过程Fig.5 Surface drift process of oil under different environments

图6 不同环境下原油漂移水平距离变化Fig.6 Variation of oil surface drift distance under different environments

图7 不同密度的原油上升高度变化Fig.7 Variation of rising height of oil with different densities

图8为不同密度原油水平扩散距离随时间的变化曲线。由图8可见，高密度原油上升到水面时的扩散范围更大，原油到达消波区的时间随着原油密度的增大而增大，密度为750 kg/m3的原油扩散至消波区的时长为127 s，而密度为970 kg/m3的原油扩散时长为180 s。原油水平漂移距离基本随时间呈线性变化，且高密度原油漂移速度更慢，高密度和低密度原油的水平扩散速率差别很大，密度分别为970和750 kg/m3原油的水平扩散速率分别为0.51和0.70 m/s。

图8 不同密度的原油漂移水平距离变化Fig.8 Variation of surface drift distance of oil with different densities

3.3 泄漏量对泄漏原油扩散漂移范围的影响

水下管道原油的泄漏量主要取决于泄漏速率和泄漏口孔径。图9为不同泄漏速率的原油上升高度变化曲线。由图9可知，泄漏速率越大，原油到达水面时间越短；当泄漏速率为2 m/s时，原油扩散上升到水面时长为25 s；泄漏速率为6 m/s时，上升时长仅为14 s。在海底附近高度，原油上升高度曲线呈线性变化，在水面附近高度，由于受到海浪作用，上升高度曲线开始出现起伏，原油上升速度减慢。

图9 不同泄漏速率下原油上升高度变化Fig.9 Variation of oil rising height under different leakage rate

图10为不同泄漏速率条件下原油水平扩散距离随时间的变化曲线。由图10可知，从泄漏速率快的原油到达水面时的扩散速度更快、扩散水平范围更广，原油在水面漂移到消波区所消耗的时长相差不大，即水面原油水平扩散速率基本相同。

图10 不同泄漏速率下原油漂移水平距离变化Fig.10 Variation of oil surface drift distance under different leakage rate

图11显示了在波浪洋流环境、泄漏速率为4 m/s时，不同泄漏孔径对漏油上升扩散高度变化的影响。由图11可知，原油上升扩散速度随着泄漏孔径增大而增大，这是由于大孔径泄漏口所溢出的原油初始动量更大，原油从0.1 m孔径泄漏口溢出上升到达水面的时长为14 s，0.01 m孔径泄漏口溢出原油到水面的时长是其2.5倍(35 s)。在海底4 m高度内，原油上升主要受自身惯性力影响，上升高度随时间线性变化，0.1 m孔径所对应的曲线变化斜率最大；随着原油上升高度的增大，波浪洋流成为影响原油扩散上升的主要因素，且原油上升高度曲线会有明显波动变化；当原油升高到接近水面时，其上升速度会显著降低，波浪水粒子的振动导致了高度曲线斜率的减小。

图11 不同泄漏孔径下原油上升高度变化Fig.11 Variation of oil rising height under different diameter of leak hole

图12显示了不同泄漏孔径溢出的原油水平扩散距离随时间的变化。由图12可知，从大孔径溢出的原油到达水面时的扩散速度更快、扩散水平范围更广，对海洋环境污染更严重；水面溢油扩散到消波区的时长随泄漏孔的增大而减小，对于小泄漏孔0.01 m的工况，溢油水平扩散速率仅为0.4 m/s；对于大泄漏孔0.1 m的工况，水平扩散速率高达0.7 m/s。

图12 不同泄漏孔径下原油漂移水平距离变化Fig.12 Variation of oil surface drift distance under different diameter of leak hole

4 结论

1)采用CFD软件建立了水下输油管道泄漏扩散漂移预测模型，可对不同水下环境、原油密度和泄漏量的工况进行泄漏扩散漂移模拟分析，从而得到不同工况下泄漏原油扩散上升至水面的时间和水面溢油的区域范围等重要信息。

2)在静水环境中，水下漏油的柱状油流以喷泉的形态呈现，整体关于漏油孔中心垂线呈现对称分布，到达水面时长最短。在洋流环境下，柱状油流在漏油孔2 m高度处分散成油滴，原油上升的同时向下游迁移，增加了原油的上升距离。在波浪环境下，波浪水粒子增加了原油在水面漂移速率。在波浪洋流环境下，泄漏原油在水下扩散范围更大，原油上升速率最慢但水面漂移速率最快。

3)在接近海底高度的范围内，泄漏原油的惯性主导其上升速率，随着原油不断升高远离海底，原油在水下的扩散则主要受波浪洋流影响，波浪在水面附近的作用更明显，波浪水粒子的振动会减小原油上升速率。

4)原油密度对原油水下上升扩散过程影响较大，对水面原油水平漂移影响较小，低密度原油在水下上升速度快，需尽快进行响应处置。高密度原油上升到水面时的水平扩散范围较大，需考虑海底污染治理。

5)海洋环境对水面原油水平扩散速率影响较大，而泄漏速率和泄漏孔径则主要会影响原油水下上升扩散速率，其中泄漏孔径对水面原油水平扩散速率也有较大影响，大泄漏孔径、高泄漏速率的原油在水下扩散速度更快、扩散范围更广，会导致更严重的事故后果。