曹 权，王 吉，张 行，余东亮，蒋 毅

(1.中国石油大学(北京) 机械与储运工程学院，北京 102249； 2.国家管网集团西南管道有限责任公司技术中心，四川 成都 610000)

0 引言

管道运输具有运量大、密封性好、能耗小、长期运行稳定等优点。截至2017年底，中国油气长输管道总里程位列全球第3名，其中原油管道总里程达到3.2万km[1]。我国地形复杂，西南地区地貌以高原与河流为主，长距离铺设输油管道常需要穿越河流。穿河管道长期受到河床泥沙冲蚀，河水冲刷与腐蚀，极易发生泄漏，泄漏后原油进入河流，在河水的携带下会快速顺流而下，短时间内造成大范围的环境污染[2]。2003年8月5日，停泊在吴泾电厂的中海“长阳”轮被1艘小船碰撞，油箱受损导致85 t燃油泄漏，在江面造成近15万m2的污染[3]。如何快速阻止泄漏原油顺河而下的漂移，设置围油栏回收泄漏原油，对于相关部门做出应急决策，减轻事故风险具有重要意义。而指导应急决策的核心环节就是对原油顺河漂移过程的准确把握。

关于水面溢油，国内外学者开展相关研究，大多数研究针对海洋表面展开。针对海洋溢油，已开发出一系列溢油模拟系统，最新1代溢油软件如GNOME[4]，OILMAP[5]，MIKE21/3SA[6]等，基于油膜扩展理论[7]结合油粒子模型[8]，采用随机游走法和拉格朗日粒子追踪法预测油粒子时空分布，同时加入海面溢油的蒸发、溶解、乳化、沉降等过程[9]。杜文强[10]采用泊松圆盘采样算法，在油膜达到最大面积时将油膜离散为油粒子，解决油膜分散问题。现有研究多集中于海洋，港口等区域，对浅窄河道溢油研究极少。部分学者利用MIKE软件中的水动力模型及溢油模块对河流溢油进行了二维数值模拟，得到河道溢油的漂移轨迹及浓度分布等信息[11-12]。

随着计算流体力学技术发展，众多学者利用计算流体力学方法建立一系列水下原油泄漏运移模型。Yapa等[13]、Johansen[14]基于拉格朗日积分法建立水下原油喷射三维模型。Tkalich[15]采用欧拉法对原油海底扩散过程进行数值模拟，并预测溢出点油膜厚度。曹学文等[16]建立水下二维水箱模型，采用多相流VOF(Volume of Fluid)模型对水下原油管道泄漏进行二维非稳态数值模拟，研究不同漏油速率、流速等因素对水下原油扩散的影响。Cloete等[17]基于离散相DPM模型(Discrete Phase Model)，结合VOF模型，对水下天然气管道泄漏扩散行为进行模拟。

本文选用国内某一原油管道穿越的真实河流模拟原油在河面的漂移过程，采用Mixture(混合)多相流模型，计算获得原油泄漏后的漂移范围，为管道泄漏后制定及时准确的拦截方案提供理论依据。

1 数值计算模型

研究对象选取中国西南地区某条河(记为S江)，有1段原油管道穿河而过，穿越段地处平原，河流落差小，河水较浅，可以假设穿越段河水深度不变。原油在河流表面漂移主要受到河水流动的影响，忽略蒸发、溶解、乳化和沉降等因素的作用。水下原油管道发生泄漏后，上浮至河流表面的过程占据较少的时间，本文忽略原油上浮过程，从原油到达水面形成油膜之后开始考虑。由于原油被河水裹挟流动，原油在河面的漂移轨迹可近似为泄漏点位置的河水表面流体微团的流动轨迹，因此本文将在泄漏源位置处添加被标记的水分子质量源项模拟原油泄漏到达河流表面的过程，通过记录这些被标记的水分子在河流表面的移动得到被主流河水裹挟的原油在河面扩散漂移的过程。

1.1 多相流模型

多相流模型Mixture允许相间渗透，与标记水分子进入主流过程原理一致。河水与泄漏源位置产生的被标记的水均视为不可压缩流体，混合模型的连续性方程如式(1)～(3)所示：

(1)

(2)

(3)

式中：t为时间，s；ρm为混合密度，kg/m3；为哈密顿算子；为质量平均流速，m/s；k为相数；αk为第k相的体积分数；ρk为第k相密度，为第k相滑移速度，m/s。

混合相的动量方程如式(4)～(5)所示：

(4)

(5)

第k相的体积分数计算公式如式(6)所示：

(6)

将单相的标准k-ε模型扩展至多相流模型，将单相流的压力速度耦合SIMPLEC算法扩展至多相流中，利用Fluent求解器对方程进行数值求解。

1.2 边界条件

河流入口流速计算公式如式(7)所示：

(7)

式中：v为河流入口流速，m/s；Q为河水流量，m3/s；A为河流入口截面积，m2。

管道泄漏质量流速计算公式如式(8)～(10)所示：

(8)

(9)

(10)

式中：G为泄漏质量流速，kg/(m2·s)；As为单个泄漏质量源面积，m2；vl为管道泄漏质量流量，kg/s；Al为管道裂口面积，m2；Cv为流体泄漏系数，为经验常数，Cv=0.60；g为重力加速度，取9.81 m/s2；Δh为泄漏点处流体压力水头，m；p为输送管道内介质压力，Pa；p0为环境压力，Pa；ρ为泄漏流体密度，kg/m3；v0为管内介质流速，m/s。

河岸为无滑移固体边界，河道入口、出口分别为速度入口边界、自由出流边界。泄漏口定义为泄漏质量源项。

1.3 物理模型

穿越段河流二维物理模型如图1所示。在得到原始河流地图后，对原始图像进行截取和裁剪，只保留河道信息，得到仅包括河道区域的目标图像，河道信息地图如图1(a)所示。为获得河流边缘坐标信息，对目标图像进行二值化处理，二值目标图像如图1(b)所示。二值目标图像边缘坐标利用Canny算子提取，为重建连续通道边界，将这些离散点连接起来形成NURBS曲线。由于原始图片分辨率较高，直接重建的河道细节过多，存在尖锐边缘。为平滑河道形状，在重建过程中，沿边界均匀间隔选取1/10的离散数据点，以较少的数据表示主要几何特征，平滑插值获得可用于后续数值计算的河道形状。河道二维模型如图1(c)所示。

图1 穿越段河流二维物理模型Fig.1 Two-dimensional physical model of crossing section river

1.4 网格生成

自西向东在管道穿越位置上取10个泄漏点，即在水面等距建立10个泄漏质量源项，分别编号为1～10，代表原油溢出区域。对物理模型进行结构化网格划分，对河岸附近及泄漏口附近进行网格加密处理，最终生成网格总量277 952。建立的泄漏口模型及附近的网格如图2所示。

图2 泄漏口模型Fig.2 Model of leakage holes

1.5 模拟工况

河流、原油及管道参数见表1。

穿越河流管道发生泄漏位置不同，导致河面溢油区域位置不同，考虑河流表面河水速度分布以及河岸对原油漂移的影响，将泄漏位置划分为河岸边缘、近岸处及河流中部3个区域。S江全年平均流量为250.1 m3/s，计算段河流入口平均流速为2 m/s，考虑S江在枯水期及丰水区河水流速不同，将河流入口流速分别设置为1.58，2.58，3.58 m/s。水下原油管道工作环境恶劣，本文泄漏量按照管道发生完全断裂计算得到河流表面泄漏质量流速范围。模拟工况设置见表2。

2 计算结果分析

2.1 泄漏位置对原油漂移过程的影响

《地表水环境质量标准》(GB 3838—2002)[18]规定：地表水中石油类污染物体积分数大于5.638×10-8为污染区。河道边缘位置发生泄漏后，原油在河流表面的漂移过程如图3所示。原油在河流边缘处发生泄漏时，河流变窄与河流曲折对漂移过程影响不大，原油会紧贴河岸向下游漂移。河流扩张段对原油漂移影响较大，在西侧(泄漏口1)发生泄漏后，原油漂移至河道扩张段，会在扩张段西侧形成漩涡。与西侧发生泄漏相比，在东侧(泄漏口10)发生泄漏后，原油漂移进入河流扩张段未出现旋涡，但会在扩张段东侧末端近岸处出现堆积。经过扩张段后，原油漂移横向扩散范围变宽，基本充斥整个河道继续向下游漂移。

表1 河流、原油及管道参数Table 1 Parameters of river, crude oil and pipeline

河流近岸处发生泄漏时原油在河流表面的漂移过程如图4所示。原油向下游漂移过程中会逐渐向河岸靠近，最终紧贴河岸。在西侧近岸处(泄漏口3)发生泄漏时，原油漂移进入河流扩张段，会出现旋涡，在河面大量堆积，几乎充满整个扩张段，东侧近岸处(泄漏口9)发生泄漏时，原油漂移进入河流扩张段的角度偏向河流中部，会在扩张段出现微小旋涡。经过扩张段后，原油会充斥整个河流表面继续向下游漂移。

表2 模拟工况Table 2 Simulation conditions

图3 边缘处泄漏漂移过程Fig.3 Drift process of leaking at edges

图4 近岸处泄漏漂移过程Fig.4 Drift process of leaking at nearshore

河流中部(泄漏口5，泄漏口7)发生泄漏时原油在河流表面的漂移过程如图5所示。分析可知：原油在河流中部发生泄漏后，河流变窄与河流曲折对漂移过程影响不大，原油基本沿着河流中部向下游漂移，漂移至河流扩张段，原油会充斥整个扩张段。此外泄漏位置越靠近河流中部，原油漂移进入扩张段形成的污染面积越大。

图5 中部泄漏漂移过程Fig.5 Drift process of leaking at middle

不同泄漏位置发生泄漏后原油在计算域中的污染面积随漂移时间变化如图6所示。泄漏初期，原油在河流表面污染面积与漂移时间变化曲线呈近似线性关系，在约0.85 h时，原油漂移至河流扩张段，河水流速与方向变化较大，原油在河流表面污染面积的增长速率变大。在约2.3 h时，原油漂移到计算域中的河流出口，此后计算域内污染面积基本稳定下来，随时间变化不大。此外，泄漏位置越靠近河流中部，河流污染面积增长速率越大，最终造成的污染面积越大。

图6 不同泄漏位置原油漂移污染面积变化Fig.6 Variation of pollution area by crude oil drift under different leakage locations

2.2 河流入口流速对原油漂移过程的影响

图7 不同河流入口流速下原油漂移过程(2.78 h)Fig.7 Drift process of crude oil under different river inlet flow velocities (2.78 h)

发生泄漏2.78 h后，不同河流入口流速下原油在河流表面漂移过程如图7所示。原油在河流表面漂移造成的污染面积与漂移时间变化如图8所示。泄漏位置越靠近中部，造成污染面积越大。对不同泄漏口，河流入口流速一致时，漂移至河流出口所需时间基本相同，随着河流入口流速增加，原油在河流表面的漂移速度会明显变大，在河流入口流速为1.58 m/s时，原油在约2.3 h时漂移至河流出口，在河流入口流速为2.58 m/s时，原油在约1.7 h时漂移至河流出口，在河流入口流速为3.58 m/s时，原油在约1.3 h时漂移至河流出口。对同一泄漏口，原油漂移至河流扩张段时，河水流速的变大导致油膜变薄，污染面积随漂移时间增大的速率变大，但原油漂移至河流出口后，最终在计算域内河流表面形成的污染面积会趋于一致，说明河流入口流速对原油漂移在计算域内最终形成的污染面积影响不大。

图8 不同河流入口流速下原油漂移污染面积变化Fig.8 Variation of pollution area by crude oil drift under different river inlet velocities

2.3 泄漏质量流速对原油漂移过程的影响

发生泄漏2.78 h后，不同泄漏质量流速下原油在河流表面漂移过程如图9所示。原油在河流表面漂移造成的污染面积与漂移时间变化如图10所示。原油均在1.3 h左右漂移至计算区域出口，泄漏位置越靠近河流中部，造成污染面积越大。对同一泄漏口，随着泄漏质量流速增加，原油污染面积略有增加，油膜变厚。

图9 不同泄漏质量流速下原油漂移过程(2.78 h)Fig.9 Drift process of crude oil under different leakage mass flow rates (2.78 h)

图10 不同泄漏质量流速下原油漂移污染面积变化Fig.10 Variation of pollution area by crude oil drift under different leakage mass flow rates

3 结论

1)采用CFD软件建立原油河面漂移预测模型，对不同泄漏位置、不同河水入口流速、不同泄漏质量流速情况下发生泄漏后原油在河流表面的扩散漂移过程进行模拟分析，从而得到泄漏原油的漂移范围及污染面积等重要信息。

2)穿越河流管道泄漏位置对原油在河流表面漂移过程影响较大。若河流形状未出现扩张，在河流边缘及近岸处发生泄漏时，漂移过程基本贴近河岸；在河流中部发生泄漏时，基本不会漂移至河道边缘；突遇河流形状出现扩张，原油会在扩张段出现漩涡，泄漏位置越靠近河流中部，在扩张段形成的污染面积越大。

3)原油在河流表面漂移速度主要受到河水入口流速影响，与泄漏位置与泄漏质量流量关系不大，河流入口流速越大，漂移速度越快。泄漏质量流速大小对原油在河流表面污染范围影响不大。