随着我国海上石油资源的不断开采，海底输油管道铺设越来越多，其安全运行问题备受关注。在复杂的海洋环境中，一旦管道破损发生溢油事故，将会对海洋环境造成极大的破坏^[1]。因此，研究海底管道溢油的迁移扩散行为，有助于制定及时有效的应急措施，最大限度地减少溢油对海洋环境的污染。

对海底管道溢油扩散的研究，国外起步较早，Fay等^[2]提出了油膜三阶段扩展理论，得到相应的扩展公式，但未考虑海水紊流作用。Hirst^[3]对二维和三维浮射流溢油进行了数值模拟，并与实验数值进行了对比验证。Fanneløp等^[4]进行了静水条件下的溢油浮射流实验，Fan^[5]研究了不分层流动环境下的溢油实验，大部分海底溢油数值模型将其作为模型验证的数据资料。Brandvik等^[6]和Zhao等^[7]等均开展了大型水槽的溢油实验，分析水下油滴粒径分布的特征。Raznahan等^[8]建立了三相瞬态流动的溢油数值模型，研究了不同破冰下的溢油扩散特性。我国对于海底管道泄漏溢油的研究起步较晚。王晶等^[9]建立了海底管道小孔溢油的数值模型，可以预测油膜出现的位置。高清军等^[10]基于Fluent分析比较了不同管内压力和不同水流速度条件下的溢油行为；管永义等^[11]考虑了风波流共同作用下的海底管道溢油，但只分析了溢油的上升过程，并未考虑其他海况。曹学文等^[12]建立了二维的波流环境下海底管道溢油扩散数值模型，分析了海洋环境、油密度对原油扩散漂移行为的影响。欧阳锡钰等^[13]建立了二维的海面水动力模型，研究溢油在海面的扩散漂移过程。宋超等^[14]建立了埋地输油管道泄漏模型，分析了溢油在多孔介质中运移的影响因素与内在机理。对于物理实验方面，钱国栋等^[15]构建了水下溢油模拟实验装置，使用蓬莱19-3原油开展了在静水环境下不同溢油条件的物理实验。Zhu等^[16]对不同流速和压差下的溢油水下扩散进行了实验研究，并根据溢油实验提出了两种形式的破裂和聚结。廖国祥等^[17-18]开展了静水与非均匀流环境下的溢油浮射流实验与数值模拟分析，其数值结果与实验观测结果总体上符合较好，后又开展了横流环境下的水下溢油实验，为后续的物理实验和数值模拟研究提供参考。林先传等^[19]建立了小型溢油实验水槽模型进行实验验证，并数值模拟了波、流及波流共同作用下原油的水下扩散过程，但未考虑风速的影响。

目前，大多数学者对海底管道溢油问题的研究仅限于二维模型，未考虑溢油在海面上的横向扩散宽度变化，考虑复杂海况的影响也是非常少。鉴于此，本文以海底输油管道泄漏为例，基于VOF方法对N-S方程、k-ɛ湍流模型进行了三维数值计算，探究了溢油到达的海面时间和溢油的水下运动过程，细致分析了海流、波浪、海面风等因素对溢油漂移扩散的影响，计算了溢油在海面的横向扩散宽度及水平漂移距离变化，分析了溢油的漂移扩散运动规律，可以为相关部门处理突发溢油事故和制定应急策略提供理论指导。

1   材料与方法

1.1   基本控制方程

在三维空间中，流体质点的运动受质量、动量方程控制，具体形式为：

连续性方程：

动量方程（N-S方程）：

式中：i，j=1、2、3，代表三个维度；ρ表示流体密度（kg/m³），$ \rho ={\displaystyle\sum }_{1}^{3}{\varphi }_{q}{\rho }_{q} $；$ q $表示第$ q $相，$ {\varphi }_{q} $表示第$ q $相的体积分数；$ \overline{u} $表示混合质量平均速度，m/s；$ \mu $为流体动力黏度；$ P $为流体压力；$ {g}_{i} $为$ i $方向上的重力加速度；$ -\rho \overline{{u}'_{i}{u}'_{j}} $表示雷诺应力，表达式为：

式中：$ {\delta }_{ij} $为Kronecker函数。

湍流部分所用模型为Realisable k-ɛ湍流模型。Tauseef等^[20]研究表明，Realisable k-ɛ湍流模型具有更高的精度；此外，Zhu等^[21]的数值研究采用了流体体积法（VOF）结合Realisable k-ɛ湍流模型，计算结果与实验数据吻合良好，因此，本文采用Realisable k-ɛ湍流模型捕捉湍流的主要特征。Realisable k-ɛ湍流模型包含湍动能和湍流耗散率两个方程^[22-23]，即

式中：

式中：$ k $、$ \varepsilon $分别为湍动能和湍流耗散率；$ {\mu }_{t} $为湍流黏度；$ \upsilon $为流体运动黏度；$ S $为应变速率；$ {G}_{k}\mathrm{、}{G}_{b} $为平均速度梯度和升力引起的湍动能产生项；$ {Pr}_{t} $为普朗特数，取值为0.85；${C}_{\mu }、\;{C}_{1\varepsilon }、 {C}_{2}、\;{C}_{3\varepsilon }$为经验系数，分别为0.09、1.44、1.9和0.9；$ {\sigma }_{k} $和$ {\sigma }_{\varepsilon } $分别为$ k $和$ \varepsilon $的普朗特数，取值为1.0、1.2。

1.2   VOF方法

本文研究海底管道溢油数值模拟中，计算域存在油、水、气三相流体，由于流体互不相容，相界面明显，对于油、水、气三相的相界面捕捉，适合VOF方法。VOF方法是一种欧拉多相流模型，通过引进相体积分数$ {\varphi }_{q} $这一变量，求解一套动量方程并且跟踪穿过计算区域的每一种流体的体积分数来模拟两种或多种不相混合的流体^[24]，实现对每一个计算单元相界面的追踪。流体体积分数$ {\varphi }_{q} $满足公式（7）、公式（8）。

式中：$ q $为各相编号，取值1、2、3代表油相、水相和气相。

1.3   数值模型的建立

本文研究的溢油场景根据实际海洋尺度设置环境水深，考虑小尺度水域复杂海况下的溢油行为，由于实际海洋溢油的输运空间较大，本文仅考虑溢油发生后至海面漂移扩散一段时间内的溢油输运行为，对溢油的分布区域即对局部海域简化为水槽模型，在局部海域，海面风速、海流速度视为均匀稳定。

本文构建的数值模型（图1）长度为260 m、宽度为10 m、高度为23 m。由于渤海海湾一些油气管道位于水深几米到几十米，因此，本文构建模型的水深为18 m，泄漏口距左端速度进口60 m。泄漏口位于管道正上方，根据IOGP发布的《风险评估数据指南》^[25]，海底管道等效泄漏孔径分为小孔（d<20 mm）、中孔（20 mm<d<80 mm）、大孔（d>80 mm）和完全破裂，其中，小孔和中孔所占比例达90%。因此，本文选择泄漏孔直径为0.06 m，以开展复杂海况下海底管道泄漏的数值模拟研究。

图  1  计算域示意

注：$ {u}_{c} $为海流速度；$ {L}_{s} $为计算域长度；$ W $为计算域宽度；$ H $为水深；$ {H}_{air} $为空气域高度；$ {h}_{r} $为消波区高度，设置为两倍波高；$ {L}_{r} $为消波区长度，为1.5倍波长

Fig.  1  Schematic of the computational domain

下载: 全尺寸图片 幻灯片

网格采用切割体网格，重点对泄漏口周边及气液交界面进行局部加密处理，溢油区域网格为0.015 m×0.015 m，其余部分网格尺寸适当增大。经网格无关性验证后确定的网格数为2871423，具体的网格划分及局部加密如图2所示。

图  2  网格划分

Fig.  2  Diagram of grid division

下载: 全尺寸图片 幻灯片

1.4   边界条件

本文选择用Stokes五阶波来构造波浪环境，使用斯托克斯波理论的五阶近似对五阶波建模，与通过一阶方法生成的波相比，此波更接近于真正的波^[26]。对于Stokes五阶波，其速度势方程为：

波面方程为：

式中：$ \omega $、$ k $分别为圆频率、波数，即$ \omega = 2\pi /T $，$ k=2\pi /L $，$ T $为波周期，$ L $为波长；$ c＝L/T $；$ {\lambda }_{n} $为速度势系数。

计算域左侧为速度进口边界条件，为水流、波浪及风的入口，速度进口可表示为：

其中：

式中：$ {u}_{c} $为均匀流水平流速，$ \mathrm{m}/\mathrm{s} $；$ {u}_{x} $为波浪水平方向的速度，$ \mathrm{m}/\mathrm{s} $；$ {u}_{y} $为波浪竖直方向的速度，$ \mathrm{m}/\mathrm{s} $；$ H $为水深，$ m $；$ {u}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}} $为海面风速。

计算域的前后边界为对称平面边界；底部设置为壁面，无滑移条件。油品泄漏口为速度进口边界条件；右侧和上方设置为压力出口，计算域上方边界定义压力为0 hPa，右侧出口边界采用线性静压曲线来满足静水压力分布$ {P}_{\mathrm{o}\mathrm{u}\mathrm{t}}=\rho g\left(H-y\right) $。

1.5   数值验证

1.5.1   波形验证

以波长为36 m为例，验证了波高为1 m时的波形。图3是波长为36 m时，X=36 m处数值波形变化与理论值进行对比，结果表明，数值波形与理论波形吻合情况较好，数值解的峰值相对于理论解峰值误差较小，但在合理范围内。图4为计算域波形图，可以看出右侧边界波形趋于平缓，消波效果良好，无波的反射，可用于研究不同波长对于溢油扩散过程的影响。

图  3  波形对比（L=36 m）

Fig.  3  Comparison of waveform

下载: 全尺寸图片 幻灯片

图  4  计算域波形（L=36 m）

Fig.  4  Computational domain waveform

下载: 全尺寸图片 幻灯片

1.5.2   模型验证

为验证本文所建立数值模型的准确性，将数值模拟结果与Zhu等^[16]的实验数据进行对比。Zhu等^[16]的实验是基于小型实验水槽来探究溢油的水下扩散特性，溢油实验环境为符合对数分布的剪切流环境，水槽为1.5 m×0.1 m×0.3 m，水深为18 cm，泄漏口径为2 mm，位于管道的上部，实验油品选用染色白油，密度为905.7 kg/m³，黏度为0.0606 Pa·s。

根据实验几何条件与油品物性建立相匹配的数值模型。图5为实验和模拟条件下不同时刻的溢油结果，图6为实验和模拟条件下溢油垂直上升高度与泄漏时间的关系。经对比可以发现，数值计算结果与实验结果吻合良好，故模拟所采用的多相流和湍流模型是可行的，可用于后续溢油扩散规律的研究。

图  5  实验与模拟溢油体积分数云图

Fig.  5  Comparison of oil spilled experiments and simulation results

下载: 全尺寸图片 幻灯片

图  6  溢油垂直上升高度与泄漏时间的关系

Fig.  6  Vertical rise height of oil spilled with leakage time

下载: 全尺寸图片 幻灯片

2   结果与讨论

2.1   工况参数选取

本文所选工况参数是根据我国渤海海洋条件情况选取的，根据陈红霞等^[27]对中国近海波浪季节特征的研究，渤海冬季波浪有效波高为1～1.2 m，夏季波高平均低于冬季。因此，本文在研究波浪对溢油的影响时选择波高为H_W=1 m（振幅为0.5 m），波长定义分别为36 m、54 m和72 m；海流速度根据渤海特征，选取0.1～0.3 m/s，结合渤海海域风特性统计，假设海上风场匀速稳定，选取海面风速为6 m/s、10 m/s、15 m/s进行研究。空气和水的密度分别为1.225 kg/m³和1025 kg/m³，黏度分别为1.8×10⁻⁵ Pa·s和1.003×10⁻³ Pa·s，原油密度为950 kg/m³。在渤海海域，漏油速度通常为0.1～10 m/s^[21]，本文研究考虑溢油泄漏速度为3 m/s，对渤海的海底管道泄漏问题开展不同波长、海面风速、海流速度下溢油扩散数值模拟研究。模拟工况参数见表1。

表  1  工况参数

Tab.  1  Parameters of working conditions

工况	海流速度（VW）/m∙s−1	波长（L）/m	风速（Va）/m∙s−1
A1	0.1	36	0
A2	0.1	54	0
A3	0.1	72	0
A4	0.1	36	6
A5	0.1	36	10
A6	0.1	36	15
A7	0.2	36	10
A8	0.3	36	10

下载: 导出CSV 

| 显示表格

2.2   海浪波长的影响

本文仅考虑溢油发生后至海面漂移扩散一段时间内的溢油输运行为，将溢油的分布区域即局部海域简化为水槽模型，波浪选为Stokes五阶波，波高选为1 m，未考虑波浪破碎、水滚等能量损失的影响。在海面无风时，海流速度为0.1 m/s，不同波长下溢油在X（水平偏移距离L_S）、Y（上升高度H）、Z（横向扩散宽度）方向的扩散变化如图7所示。可以看出在泄漏条件一定时，波长对溢油浮升过程影响较小，主要影响溢油至海面后的水平漂移和横向扩散过程。

图  7  不同波长下溢油X、Y、Z向的扩散随时间变化

Fig.  7  Diffusion of oil spilled in X, Y and Z directions at different wavelengths

下载: 全尺寸图片 幻灯片

溢油到达海面后的漂移扩散过程主要受波浪的传播作用影响，溢油的水平漂移距离整体呈线性变化关系。在流速一定、波长更长的情况下，溢油的水平漂移距离明显增加，当溢油时间为250 s、波长为72 m时溢油的水平漂移距离为54.52 m，是波长（36 m）的1.3倍。

溢油横向扩散宽度变化在溢油至海面以后变化较快，对于溢油上升至海面以后的扩散宽度变化进行数据拟合（拟合公式及拟合优度在图中给出），发现溢油在海面的横向扩散宽度变化整体呈对数分布，溢油扩散宽度随着波长的增加而减小。这是由于在波高一定、波长较长时海水表层水质点向下游传播的速度较快，进而带动溢油向下游传播，溢油在海面漂移距离较大，横向扩散宽度较短，使溢油在海面的漂移距离高于横向扩散。

图8为不同波长下的溢油轨迹，可以看到在溢油初始阶段，在波长较长的条件下，油滴率先发生脱落现象，溢油射流发生明显的弯曲，当溢油即将上升到水面时，在波浪的扰动下，首滴油滴跟随水质点的运动到达水面较快，与后续生成的油滴相隔较远。当溢油发生时间为250 s时，可以看到油品上升至海面后的分布情况，在波长较小时，溢油水平偏移较小，横向扩散最大；当波长增加时，溢油横向扩散宽度较小，水平偏移距离较大。

图  8  不同波长下的溢油运动轨迹

Fig.  8  Oil spilled trajectories at different wavelengths

下载: 全尺寸图片 幻灯片

2.3   海面风速的影响

考虑海面风速的影响，当泄漏条件、波流海况条件一定时，不同海面风速下溢油在X、Y、Z向的扩散变化如图9所示。可以看到，海面风对溢油的上升过程的影响并不明显，主要影响溢油在海面上的漂移扩散过程。当溢油达到海面时，由于海面风的存在进一步带动海面水质点的运动，风速越大，海面水质点的流动速度也越大，使溢油的漂移距离越远，当风速为15 m/s、原油泄漏250 s时，溢油已到达计算域消波区处，漂移距离相比于海面无风时扩大4倍。溢油的横向扩散宽度随海面风速的增大而减小，宽度变化整体呈对数分布，在扩散初期溢油横向扩散速度较快，而后趋于平稳；当风速较大时，相同泄漏时间内溢油横向扩散宽度相比无风时小。

图  9  不同海面风速下溢油X、Y、Z向的扩散随时间变化

Fig.  9  Diffusion of oil spilled in X, Y and Z directions at different wind speeds

下载: 全尺寸图片 幻灯片

图10为不同海面风速下的溢油运动轨迹，可以看出海面风速对溢油在水下的运动影响较小，当风速较大时，首先到达海面的油滴与后续形成的油滴相距较远，对溢油的射流长度无明显影响。溢油至海面后，油膜在风波流的共同作用下向下游扩散，在风速较大时，油膜具有较小的扩散宽度和较大的漂移距离。当溢油事故发生时，需充分考虑海面风速的影响，避免溢油污染范围的进一步扩大。

图  10  不同海面风速下的溢油运动轨迹

Fig.  10  Oil spilled trajectories at different wind speeds

下载: 全尺寸图片 幻灯片

2.4   海流速度的影响

当风速、波浪条件一定时，研究不同海流速度对溢油扩散过程的影响。图11为不同海流速度下溢油在X、Y、Z向的扩散变化，可以看到，在泄漏条件一定的情况下，溢油的上升速度随海流速度的升高而降低，达到海面所需的时间也会更长。溢油在水下的偏移距离主要受海水流速的影响，流速越大，溢油在水下的偏移距离越远。在流速为0.30 m/s，溢油至水面时，溢油在水下的偏移距离为14.44 m，水平偏移距离是海流速度（0.10 m/s）的3.1倍。

图  11  不同海流速度下溢油X、Y、Z向的扩散随时间变化

Fig.  11  Diffusion of oil spilled in X, Y and Z directions at different current speeds

下载: 全尺寸图片 幻灯片

当溢油到达海面后，溢油的漂移扩散过程受风波流的共同作用，在海面风速、波长一定时，溢油在海面的漂移速度随着海流速度的增加而增加，当溢油发生250 s时，溢油已到达消波区处，溢油的漂移距离整体呈线性关系变化，相同泄漏时间内，海流速度越大，溢油的漂移距离越远。溢油在海面的横向扩散宽度变化呈对数曲线关系分布，在海流速度较大时，溢油横向的扩散速度略低于水流速度较小的工况。

图12为不同海流速度下的溢油运动轨迹，可以看出海流速度主导着溢油在水下的偏移和海面上的漂移运动。海流速度越大，溢油在水下的运动轨迹发生的偏转角度越大，水下偏移距离越远，油滴分布较为分散。当溢油泄漏时间为250 s时，海流速度较小时，油膜的宽度较大，溢油水平漂移距离较短。随着海流速度的增加溢油在海面漂移的距离变远，横向扩散相比流速较小时扩散程度较小。

图  12  不同海流速度下的溢油运动轨迹

Fig.  12  Oil spilled trajectories at different current speeds

下载: 全尺寸图片 幻灯片

3   结 论

（1）海流速度显著影响溢油的浮升过程，在泄漏条件一定的情况下，溢油的上升速度随着海流速度的升高而降低，达到海面需要更长的时间，溢油在水下的偏移距离也显著增加，溢油轨迹发生的偏转角度较大。而海浪波长和海面风速对溢油的浮升过程影响并不明显。

（2）溢油在海面的漂移扩散过程受波长、风速、海流速度的影响，溢油上升至海面后，油品在横向的扩散宽度随时间变化呈对数分布，水平漂移距离基本满足线性变化。随着波长、风速、海流速度的增加，溢油在海面的横向扩散速度变小，扩散宽度变短，溢油的水平漂移速度较快，漂移距离较远。