水下输气管道泄漏扩散特性模拟研究

王少雄，李玉星，刘翠伟，梁杰，李安琪，薛源

（中国石油大学（华东）山东省油气储运安全省级重点实验室，山东青岛266580）

引 言

天然气的管道输送过程中不可避免地会遇到水体环境、河流冲刷、水体对管道的腐蚀，以及设计安装不当、第三方破坏等人为因素，使得水下输气管道面临的风险增大，泄漏事故也逐年增加[1-3]。当天然气连续地从泄漏孔进入水体中时，气体由于自身与周围水域的压力差以及气体与水体交界面表面张力的作用，将会破裂成一个个气泡。气泡受到浮力的作用而上升，上升的气泡将会夹带周围的海水一起向上运动，从而形成气泡羽流[4]。水下输气管道一旦发生泄漏，上升到水面的气泡羽流可能会引起火灾、造成水域污染，甚至会诱发二次事故对人身财产安全造成威胁。因此，科学地认识水下输气管道泄漏所形成的气泡羽流扩散规律，对于及时地采取防治措施和制定风险评价方案具有一定的意义。

目前对于输气管道泄漏的研究大多针对陆地的架空管道[5-7]，对于水下输气管道泄漏的模拟研究主要包括积分模型和计算流体动力学模型。积分模型把气泡羽流中的气泡和卷吸进来的水体看作一个多相控制体，依据相似性原理，假设羽流径向的质量和动量守恒为某个概率分布从而得到特征参数随羽流轨迹变化的常微分方程。Morton 等[8-10]提出了一种积分模型，该模型假设羽流速度和密度分布具有相似的形式，这一假设简化了羽流轴向上的质量和动量守恒的计算方法。Ditmars等[11-12]进一步完善了Morton 等的模型，在模型中考虑了气泡之间的滑移速度，并且认为气泡羽流对水体产生了拖动作用。Milgram 等[13-14]提出了自己的积分模型，并且在估计气泡羽流的动量时引入了动量放大系数的概念。Billeter等[15]通过假设“反向停滞流”改进了积分模型，并且提出了新的预测羽流轮廓的方法。计算流体动力学模型是通过求解流体流动的Navier-Stokes 方程来研究流体动力学问题，它直接模拟气泡夹带和湍流运输产生的影响。Cloete 等[16]建立了CFD 数值模型来研究水下输气管道断裂引起的泄漏扩散。Olsen 等[17]提出了一种基于欧拉-拉格朗日准则的CFD 模型来模拟开阔水域的气泡羽流的释放，并将模型与Milgram 等[13-14]的实验结果进行对比验证，发现模型具有良好的适应性。Tessarolo 等[18]提出了一种基于拉格朗日法的数值模型来模拟深水井喷的油气扩散，并通过模型验证了三种夹带关系式，发现JETLAG 关系式与实验结果最为吻合。Wu 等[19]用CFD 模型来模拟海底气体释放，并且评估了四种数值方法（雷诺时均N-S、非稳态雷诺时均N-S、LES 大涡模拟、SAS 尺寸自适应模拟）在捕捉上升气泡羽流的特征行为上的适用性。国内学者也针对天然气在水体中的泄漏扩散进行了相关的模拟研究。李长俊等[20]采用CFD建模的方法研究了穿越河道天然气管道发生泄漏时不同的泄漏朝向和水流速度下气体的水平迁移距离变化规律。李新宏等[21-22]利用Fluent 建立了海底管道泄漏二维扩散模型，研究了不同海流流速和泄漏孔径对气体在水体中的扩散形态和泄漏速率的影响。张焕好等[23]建立了二维轴对称水下超声速气体射流的数值计算模型,重点研究了水下超声速气体射流的初始流动结构。相比于模拟研究，国内外专家对于气体水下淹没射流进行了大量实验研究，并取得了丰富成果。Topham[24]在一个水深为23 m 的水池中做了现场实验，测量了气泡羽流在轴线上和横断面上的速度变化规律。Lees 等[25]进行了现场和实验室实验来研究水下气体释放时的气体浓度分布，发现气体浓度呈现高斯分布。王志刚等[26]进行了室内气体水下泄漏实验，研究了不同工况下的羽流直径的变化规律。郭强等[27]进行了水下超声速气体射流实验，研究了不同工况下的水下高速气体垂直射流的演化过程。王超等[28-29]通过实验研究了超声速气体射流在水中的喷射所形成的流场结构，并用数值计算的方法成功模拟了射流初期气泡运动演化的复杂过程。施红辉等[30]进行了三维水下超声速气体射流实验，重点研究了射流初期气泡运动演化过程。李婷婷等[31]进行了水下竖直环形喷管出口气体射流实验，并通过Matlab 对环形喷管射流的特性进行了分析。

上述模拟研究大多侧重于特定的管路过程参数对气体的流动特性的影响而对气体上升到水面时的扩散特性以及气体在水体中具体的扩散过程研究较少，并且目前还没有关于两种模型对水下输气管道泄漏扩散特性的预测结果准确性的对比评估。因此本文通过建立水下输气管道泄漏三维CFD 模型和积分数学模型来研究水下输气管道泄漏后的扩散特性，分析水下输气管道发生泄漏时形成的气泡羽流在水体中扩散的速度、羽流半径和泉涌高度的变化规律，并利用实验数据对两种模型的准确性进行了对比验证，为水下输气管道发生泄漏后的风险评价以及后续的模拟研究提供一定的参考。

1 水下气体泄漏CFD模型理论基础

欧拉-拉格朗日模型可以用来求解水下输气管道泄漏后的多相流动过程，该模型将泄漏气体作为拉格朗日离散相而周围的水体作为连续的欧拉相处理，这些离散相粒子在连续相中运动，并通过交换质量、动量和能量与之相互作用。欧拉参考系中的方程组求解周围水体的速度和湍流、可溶性气体的浓度以及不同连续相（水和大气）之间的界面位置，这些变量将用于对拉格朗日系中气泡运动轨迹的追踪。气泡上升过程中在浮力、拖曳力和虚拟质量力的作用下达到平衡，拖曳力和虚拟质量力也会对欧拉相进行耦合从而影响周围水体的速度，因此欧拉-拉格朗日模型具有强大的双向耦合功能，所以又称为耦合的DPM 和VOF模型。

1.1 VOF模型

在VOF 模型中，水体上方的大气和周围水体都被视为连续相，二者相互贯穿，一相的体积不能被另外一相的体积所占有。模型中采用体积率来衡量气体在水中所占的体积比，气体和水的体积之和为1。使用VOF 方法，首先要定义一种流体体积函数，通过网格内该流体的体积与网格体积的比值来确定自由面，继而追踪流体的变化。对于水下天然气管道泄漏，VOF 模型通过求解连续性方程、动量方程和能量方程来实现对水和大气两相自由表面的追踪[2]。

连续性方程如式（1）所示。

式中，αi为第i相欧拉相的体积分数；ρi为第i相欧拉相的密度，kg/m3；vi为第i相欧拉相的速度，m/s。

动量方程取决于各相的体积分数，如式（2）所示。

式中，ρ = ∑αiρi为混合相的密度，kg/m3；μ =μT+ μM为分子混合黏度和湍流黏度之和，(N·s)/m2；F为外部作用力，N。

能量方程如式（3）所示。

式中，keff表示流体的传热系数；Sh为源项。

1.2 DPM模型

DPM 模型将泄漏气体作为离散相处理，假设气体粒子在拉格朗日计算区域中的体积分数低于10%～12%，通过对离散相粒子施加一个平衡力实现对粒子运动轨迹的追踪。气体受力平衡方程如式（4）所示。

式中，ub表示泄漏气体的速度，m/s；ρb表示泄漏气体的密度，kg/m3；FD表示气泡颗粒受到水的拖曳力，N；Cd表示拖曳力系数；E0为无量纲数，用来表征气泡上升过程中形状的变化；σ表示环境水的黏度，(N·s)/m2；Fvm表示虚拟质量力，N；Cvm表示虚拟质量力系数，通常为0.5；db表示气泡颗粒的直径，m；Re表示相对Reynolds数。

1.3 气泡密度和尺寸模型

随着气泡在水体中的上升，气泡所处水深变浅，使得气泡运动过程中周围的压力发生显著变化，从而导致气泡密度随水深的降低而降低。因此，水下气体的密度是水深的函数，可以通过推导理想气体方程进行求解，如式（9）所示。

式中，Mb表示泄漏气体的相对分子质量；R 为通用气体常数，J/(kmol·K)；P 表示气泡所承受的静水压力，Pa；P0表示大气压，Pa；Hp表示气泡所处的水深，m。

由式（7）可知，气泡的大小决定了气泡的特征形状，进而影响了气泡在水中所受拖曳力的大小。同时气泡上升过程中与环境水存在质量传递，而这取决于气泡的表面积大小。现有的模型出于简化计算过程大多将气泡大小设定为常数，但气泡上升时由于发生破裂和合并使得气泡的大小时刻发生变化，因此现有的模型未能准确地预测气泡羽流的扩散特性。Laux 等[32]通过求解输运方程建立了一个欧拉框架下的气泡尺寸模型，如式（11）所示。

式中，τrel为松弛时间，s；为气泡的平衡直径，m；αb为气泡的体积分数；σ1为表面张力；C1= 4 为无量纲常数；C2=100 μm表示气泡的最小尺寸。

1.4 网格模型和数值计算方法

由于本文模拟的泄漏工况为水下输气管道的小孔泄漏，Fluent 中的二维网格会将泄漏孔默认为狭缝泄漏，不能真实反映泄漏孔径对天然气泄漏扩散的影响，而三维网格可将泄漏孔形状画为圆形，所以本文对水下输气管道泄漏扩散模型进行简化，泄漏孔位于模型底部的中心位置。以水深h=1.1 m，泄漏孔径d=3 mm 为例，模型所研究的区域大小为1 m×1 m×1.5 m，由于泄漏孔径相对较小，对泄漏口附近的网格进行加密处理，先对整体block 划分OBlock 网格，再对中心的block 进行O-Block 网格划分，泄漏孔向模型边界的网格线划分方式为Exponential2，Spacing2=0.0002，Ratio2=1.135，网 格数为376000个，分别选择Determinant2×2×2和Angle作为网格质量的判定标准，所有网格的Determinant 2×2×2值大于0.7，大于0.85的占93.832%，所有网格的Angle 值大于40.5°，大于60°的占70.331%，生成的网格如图1所示。

图1 网格模型划分Fig.1 Mesh model

计算域通过Patch 选项将其分为水域和空气域，空气域上部设为压力出口边界，计算域的侧面和底部设为无滑移固壁边界，底部中心半径1.5 mm的泄漏孔为空气入口，泄漏孔的初值条件在DPM 模型中设置，在DPM 模型中将入口边界条件设置为流量入口，该值与实验中的泄漏速率相对应。水下输气管道的泄漏是一种瞬态湍流过程，因此模型采用k-ε 湍流模型和非稳态压力基求解器以及适用于求解非稳态流动的PISO 算法。同时由于气体的可压缩性以及与水体之间的密度差，所以在VOF 模型中选用隐式体积力公式。对于气泡上升过程中的破裂和合并通过编写UDF 函数导入气泡密度和尺寸模型。为了对网格进行独立性检验，通过设置最大网格尺寸、网格层数、第一层网格点与端点之间的距离等参数，得到不同网格划分方案，具体方案如表1所示。图2是水深1.1 m、泄漏孔径为3 mm 时不同网格数下羽流上升时间与实验值的对比，从图中可以看出，网格数为37.6 万个和102.9 万个对羽流的上升时间影响不大，网格数越多越接近实验值，为了加快计算速度并考虑模型的准确性，选择网格数为37.6 万个的模型来模拟气体在水体中的泄漏扩散过程。

表1 网格划分方案Table 1 Mesh classification scheme

图2 不同数目网格下的羽流上升时间Fig.2 Plume rise time under different numbers of grids

2 水下气体泄漏积分模型理论基础

2.1 积分模型推导过程

当水下输气管道发生泄漏时，泄漏气体向上连续地进入水体中形成气泡羽流。积分模型假设气泡羽流截面上的速度和空隙率时均分布为高斯分布[31]，分别如式（13）和式（14）所示。

式中，v 为羽流的速度，m/s；下角标c 表示轴线处的值；r 和z 分别为距离羽流轴心的水平距离和距离泄漏点处的垂直高度，m；b 为羽流轴心到边缘的宽度，m；φ 为羽流的空隙率；λ 为空隙率分布和速度分布之比。

气相的连续性方程如下

式中，g 为重力加速度；Qg为气体的体积流量，m3/s；Hv为水深，m；Hp表示与大气压力相对应的水深，在清水条件下，该值为10.33 m，在海水条件下，该值在标准大气压下稍微低于10 m；γ 是动量放大系数；s～ =(1+ λ2)v～s为表征气体上升过程中滑移速度影响的无量纲量；z～、b～、v～ 分别为无量纲轴向坐标、无量纲羽流宽度和无量纲轴向速度，表达式如式（17）、式（18）、式（19）所示。

液相的连续性方程和气液混合物的动量方程如式（20）、式（21）所示。

在该简化的积分模型中，上述两个方程的近似解可由式（22）、式（23）表示

作用在羽流上的力包括水面上的大气压力Pa、底部的压力Pb、静止水面下水的重力Gw、泉涌的重力Gf和浮力B，如图3所示。气体上升到水面形成的泉涌高度可以通过动量守恒方程来求解，积分模型假设控制体在横向是无限大的，所以没有垂直动量通量通过外部边界。平衡方程如式（24）所示

图3 作用在羽流上的力Fig.3 Force acting on plume

令控制体的最低边界位于泄漏孔下方，底部的静水压力Pb将不会影响羽流的流动，那么Pb的值将是恒定的，因此式（24）中的项Pa、Pb和Gw可以相互抵消，简化后的方程如下

假设泉涌的密度ρf和静水面处羽流的密度ρ(Hv)相等以及泉涌底部的动能完全转化为势能，那么式（26）可以简化为式（27）。

式中，β是经验常数；v = vc(Hv)是泉涌底部轴线处羽流的速度，m/s；hf是泉涌高度,m。

则羽流径向泉涌高度分布的表达式如式（28）所示。

式中，hoffset是高斯分布基准线在静水面上的偏移量。

2.2 积分模型中的经验系数

2.2.1 多变指数n 假设气泡羽流在上升过程中气体是等温膨胀的，则多变指数n=1，因为如果是在绝热情况下膨胀将导致上升气体的温度大幅度下降，这与实际情况不符。

2.2.2 夹带系数α 夹带系数α 随着气体流速的增加而增大，Milgram等[14]在1983年提出了一个关于夹带系数的半经验公式

其中，k=0.165；A=7.598；气泡Froude 数Fr 由式（30）表示

在实验室尺度下，夹带系数的值通常较小，因为实验室条件下气体流量较小，而现场实验因气体泄漏速度较大得出的夹带系数的值也相对较大。

2.2.3 宽度比λ 宽度比λ 的范围较小，通常处于0.6～1之间，并且对羽流的影响较小。实验室尺度下得出的宽度比λ 的值通常比现场实验得出的值要小。

2.2.4 滑移速度vs在气泡羽流中，气泡相对于液体的滑移速度是气体浓度、气泡直径等多因素的函数，并且由于气泡是以气泡群的形式存在，所以气泡间的受力较单个气泡的受力复杂。但气泡滑移速度对羽流的影响很小，Levich[33]在1962 年提出当气泡直径处于0.2～1.5 cm 之间时滑移速度为28～30 cm/s,而对于直径更大的气泡滑移速度将增长到35～40 cm/s，但是在湍流羽流中大气泡通常会破碎成较小尺寸的气泡。

2.2.5 动量放大系数γ 动量放大系数是总的动量通量与平均流动所引起的动量通量之比，当羽流中的气泡与羽流的尺寸相比变得非常小时，气泡的动力以及气泡之间的相互作用就变得没那么重要，羽流流动表现为单相流，动量放大系数将趋向统一。因此，在求解上述羽流控制方程时，可以忽略该参数的影响。

表2 推荐的经验参数Table 2 Recommended empirical parameters

由上可知，羽流的形成与发展和夹带系数密切相关，其他参数的变化对结果影响很小，各个参数的范围和推荐值如表2所示。

对积分模型进行数值求解，模型计算程序框图如图4所示。

图4 模型算法Fig.4 Model algorithm

3 模型验证

为了验证CFD 模型和积分模型的准确性，采用水下泄漏实验环道装置，如图5所示，以空气作为实验介质，进行水下输气管道泄漏实验，研究不同泄漏速率下气体在水体中以及水面上的扩散特性。实验是在一个长、宽为1 m，高1.5 m 的水槽中进行的，环道总长为251.5 m，内径为42 mm。泄漏孔前后装有压力表7、9（量程为0～1 MPa，测量精度为0.5%）以及质量流量计5、10（量程为0.01～1 m3/h，测量精度为0.2%）。水槽前方布置高速摄像机（型号为FASTCAM SA-X2，拍摄速度为5000 帧/秒），可通过相机架杆调节高度，对气体在水体中的扩散过程进行拍摄，高速摄像机与PC 端连接，实时保存泄漏扩散过程。通过专业的视频图像处理软件ProAnalyst 对羽流扩散过程进行分析，将实验数据与模型进行对比，结果讨论如下。

图5 水下管道气体泄漏实验系统Fig.5 Experimental system for gas leakage of underwater pipeline

3.1 水下扩散过程分析

图6 为水深30 cm、泄漏孔径5 mm、泄漏压力300 kPa、泄漏速率0.75 m3/h时CFD 模拟的羽流扩散过程和实验的对比。由图可知，水下输气管道的泄漏过程可以分为三个不同的阶段，每个阶段都对应着一个流动区域且在每个阶段羽流的上升动力以及所产生的物理现象均不相同。首先是初始阶段，该阶段的羽流为动量射流，对应为流动形成区。气体从孔口泄漏到水体环境中，起始会在孔口形成一个完整的小气泡，气泡由最初一个突起逐渐上升变大，形成半球。随着泄漏过程的进行，小气泡逐渐上升拉长变成了椭球状，具体形态如图6(a)所示。其次是充分发展阶段。气体将从流动形成区进入定型流动区，如图6(b)所示。在这个区域内气体的自身浮力将取代初始动量占据主导作用，初始动量对羽流产生的影响变得微乎其微。气体在以管内压力作为初始动力冲出泄漏孔口进入到水体环境之后，克服了周围环境的阻力开始快速地上升和膨胀。气体在上升过程中，气泡羽流速度间断面的不稳定会引起湍动，从而把周围静止的水体卷入到向上运动的羽流当中。随着湍动的发展，被卷吸进入羽流内部的流体不断增多并且随着羽流一起向上运动，羽流的边界不断向两侧发展扩张。最后是表面流动阶段，随着气泡羽流的上升，气泡羽流所处的水深也就越浅，此时所受周围的环境压力也就越小。环境水被夹带到羽流中并被携带到表面，当羽流到达表面时，夹带的水并不会随着气体进入大气而是从水平表面向外径向流动，带动气泡远离羽流轴心，如图6(c)所示。大部分气体会继续上升进入到大气中，形成表面紊动或沸腾，称为表面流动区。

图6 羽流扩散过程Fig.6 Plume diffusion process

图7 不同水深下的羽流扩散过程Fig.7 Plume diffusion process at different water depths

图8 不同水深下的羽流半径分布Fig.8 Plume radius distribution under different water depths

当水下输气管道发生泄漏时，水深（泄漏点距离水面的高度）是泄漏气体在水下扩散时最重要也是最敏感的影响参数之一。图7 是泄漏孔径为3 mm、泄漏压力为200 kPa、水深分别是1.5 m 和1 m时CFD 模拟的羽流扩散过程。由图可知，1.5 m 水深下气泡羽流的扩散过程和水深1 m 时的扩散过程基本一致，但相比于水深1.5 m 的泄漏工况水深1 m时初始动能的作用更加显著，受到水体的阻碍作用相对较小，气泡羽流的平均上升速度更快，因此到达水面的时间也更短。图8 是泄漏孔径为3 mm、水深分别是1.5 m和1 m时CFD模拟的羽流半径分布，从图中可以看出水深1 m 时的羽流半径稍大于水深1.5 m 时相同水位高度的羽流半径。这是因为在相同的泄漏压力和泄漏孔径下，水深会改变泄漏孔外的水压。水深增加，水压变大，此时气体冲出孔口就受到更大的阻力，气体在水体中的扩张也受到了更大的阻碍，因此在水深较大时羽流半径就会随之减小。

随着气泡羽流的上升，气泡羽流所处的水深也就越浅，此时所受周围的环境压力也就越小。周围水被夹带到羽流中并被携带到表面，当羽流到达表面时，夹带的水并不会随着气体进入大气而是从水平表面向外径向流动，带动气泡远离羽流轴心，形成一个圆形气池，如图9 所示。图10 是泄漏孔径为3 mm、水深分别是1.5 m和1 m时气池半径随时间的变化。当羽流到达水面时，以气体溢出点为中心继续向外扩散，扩散半径持续增大，当到达一定程度时将会趋于稳定。从图中还可以看出水深越深，气池扩散半径越大。

图9 不同水深下的气池扩散过程Fig.9 Gas pool diffusion process at different water depths

图10 不同水深下的气池半径随时间的变化曲线Fig.10 Curves of gas pool radius with time under different water depths

3.2 羽流速度变化规律

在实验中，通过高速摄像机记录输气管道泄漏发生的过程，将拍摄的视频分帧解析成图片，用ProAnalyst 图像软件处理连续帧的图片。测量出每帧图片中羽流宽度和高度，用连续帧宽度和高度的差值除以每帧间隔时间即为单位时间内羽流的速度，取多次实验均值进行统计，将CFD 模型和积分模型计算结果与实验结果进行对比，如图11、图12所示。

图11 为气泡羽流在水深1.1 m、泄漏速率为0.21 m3/h 时，水位高度在0.5、0.7 和0.9 m 处的径向速度分布。为了避免分析和阐述上的混乱，本文引入无量纲水位高度hz的概念（水位高度与水深的比值），将不同水位高度与水深区分开来。由图可知，在较低水位高度（hz为0.45 m 和0.64 m）处，CFD 模型和积分模型计算结果和实验结果基本相同，羽流的径向速度由于夹带作用随着径向距离的增加逐渐降低。但是，在水位高度为0.9 m 处，积分模型计算结果与实验存在偏差，在羽流中心处观察到的模拟速度低于实验测量值。这是由于气泡羽流接近水面时，羽流的流动方向发生了变化，垂直方向的流动变为多方向流动，其中径向流动占据主导地位。

图11 不同水位高度处的羽流径向速度分布Fig.11 Radial velocity distribution of plume at different water depths

图12 为水深1.1 m、泄漏孔径3 mm、泄漏压力100 kPa、泄漏速率为0.21 m3/h 和水深1.1 m、泄漏孔径3 mm、泄漏压力300 kPa、泄漏速率为0.37 m3/h 工况下的羽流轴向速度分布，由图可知，积分模型预测的结果与实验结果基本相同，羽流的轴向速度随着水深的增加逐渐降低且在接近泄漏孔口的位置，轴线上羽流的速度急速衰减。这是因为在气体刚刚进入水体环境时，在泄漏口的上方气流受到水体较大阻力，此时气体因克服阻力消耗掉大量动能，因此羽流速度会迅速减小。而在之后的上升过程中，气泡羽流把动量传递给环境水体，周围的水与分散气泡沿着相同的路径一起向上运动，气液流在羽流边界的剪切作用下速度逐渐减弱，变化幅度较为平缓。而CFD 模拟计算的结果大于实验结果，但随着水位的升高，二者的差值变小，这是因为CFD 模型是基于雷诺时均湍流方法，忽略了气体在水中的溶解以及气泡颗粒之间的黏性作用力。对比图12（a）、(b)，可以看出初始泄漏速率越大，羽流轴向上升速度越大。这是因为泄漏速率越大，泄漏气体的初始动能也就越大，相应会产生较大的冲量，水压对气体的阻碍作用相对变弱，初始段动量射流对应的长度将会增大，羽流的上升速度相应变大。

3.3 羽流半径变化规律

当水下输气管道发生泄漏时，气体在管内压力作用下射入水体之中，气泡羽流在上升的过程中会形成一个倒立的圆锥体结构。图13 是泄漏速率分别为0.21 m3/h 和0.37 m3/h 工况下的羽流半径分布，由图可知，积分模型预测的结果与实验结果基本相同。随着气泡羽流紊动的发展，周围静止的水体被卷吸进气泡羽流中一起向上运动，流量沿程增大，羽流半径逐渐向两侧发展并且随着水深近似呈线性增长。但是CFD 模拟计算的结果稍微大于实验结果，这是因为在CFD 模型中羽流边缘的气泡颗粒分布较为分散，无法准确地定义羽流半径，因此在模型中的羽流半径是测量羽流中心到空隙率为2%的气泡颗粒之间的距离。

3.4 羽流泉涌高度变化规律

图12 不同泄漏速率下中线处的羽流轴向速度分布Fig.12 Axial velocity distribution at the midline at different leak rates

当气泡羽流接触到水面时，仍具有一定的速度，上升的羽流会继续带动自由水面向上运动，导致溢出点附近区域海水向上凸起形成泉涌，水面上升的高度即为泉涌高度，如图14 所示。溢出水面的气体在不同条件下的泄漏行为会对水体环境和水面上的作业和人员安危造成不同程度的影响，如被引燃甚至可能造成严重的火灾或爆炸事故，如图15 所示，因此对羽流泉涌高度变化规律的研究对于水下输气管道泄漏事件的应急决策具有重要价值。

图16 为水深1.1 m 时积分模型预测的不同泄漏速率下水面处羽流泉涌高度的径向分布，由图可知，泄漏速率为0.21、0.37和0.84 m3/h 时积分模型预测的中线处泉涌高度hf分别为0.0508、0.0748 和0.1305 m。Friedl 等[34]在2000 年提出了最大泉涌高度hpmax以及平均初始泉涌高度与积分模型预测的泉涌高度hf之间的关系式，如式（32）所示，由此可以得出积分模型预测的最大泉涌高度hpmax以及平均初始泉涌高度的值，如表3所示。

图13 不同泄漏速率下羽流半径分布Fig.13 Plume radius distribution at different leakage rates

图14 气体在水面形成的泉涌现象Fig.14 Fountain phenomenon of gas formation on water surface

图15 水面溢散燃烧现象Fig.15 Submarine gas pipeline leakage accident

图16 不同泄漏速率下水面处羽流径向泉涌高度分布Fig.16 Distribution of plume radial fountain height at different leakage rates

表3 积分模型计算的泉涌高度Table 3 Fountain height calculated by integral model

图17为实验与两种模型预测的泉涌高度对比。由图可知，在泄漏速率较小时，CFD 模拟预测的结果略低于实验值，但泄漏速率较大时二者差值变大，平均相对误差为8.9%。这是因为当泄漏速率较大时在实验中会观察到气泡羽流上升到水面时出现水花溅射现象，这会使得测量的泉涌高度值偏大，而在CFD 模拟中并未观察到这种现象。从图中还可以看出，积分模型预测的最大泉涌高度和初始泉涌高度与实验差值较大，平均相对误差为17.2%。这是因为积分模型中经验常数β对羽流的泉涌高度影响很大，通常对于羽流的自由上升损耗过程经验常数β 等于0.5，但是Friedl 等[34]将实验与数值模拟的数据进行拟合后认为β 为0.39 预测的结果最好，因此在模型中选取的经验常数β 为0.39。同时由于积分模型没有考虑在水面处水的回流以及自由水面对气泡羽流的阻碍作用，这也会使得积分模型预测的结果大于实验值。

4 结 论

通过建立水下输气管道泄漏三维CFD 模型和积分数学模型来研究水下输气管道泄漏后的扩散特性，并利用实验数据对两种模型的准确性进行了对比验证，具体结论如下。

图17 实验与模拟预测的泉涌高度对比Fig.17 Comparison of fountain height predicted by simulations and experiments

（1）水下输气管道泄漏后在水体中的扩散大致经历了三个不同的阶段：初始阶段、充分发展阶段和表面流动阶段，这三个阶段对应着三个流动区域，并且在气泡羽流上升过程中伴随着卷吸现象、紊动沸腾现象。当羽流到达表面时，夹带的水并不会随着气体进入大气而是从水平表面向外径向流动，带动气泡远离羽流轴心。该过程将会由于夹带水自身的动量而使得水面升高，从而形成泉涌。

（2）水下输气管道泄漏后，羽流的径向速度近似呈高斯分布并且随着径向距离的增加逐渐降低；羽流的轴向速度随着水深的增加而降低，且在接近泄漏孔口的位置由于受到水体较大的阻力，轴线上羽流的速度急速衰减；随着气泡羽流紊动的发展，周围静止的水体被卷吸进气泡羽流中一起向上运动，羽流半径逐渐向两侧发展并且随着水深近似呈线性增长。

（3）耦合的VOF和DPM模型能够用于水下气体泄漏扩散运动规律的描述，其对气泡羽流在水体中泄漏扩散运动规律的预测与实验基本相符，但由于模型中未考虑气体在水中的溶解以及气泡颗粒之间的黏性作用力，所以CFD 模拟计算的羽流轴向速度结果要略大于实验结果；相比于CFD 模型，由于积分模型忽略了表层水的回流以及自由水面对气泡羽流的阻碍作用，模型计算的工作量大大减少因而计算速度也更快，但计算精度比CFD 模型稍差。同时模型中的经验系数如夹带系数和动量放大系数等对模型计算结果影响较大，因此仍需要进行更多深水条件下高泄漏流速的实验以确定准确的经验系数来进一步提高积分模型的计算精度。

符 号 说 明

g——重力加速度，m/s2

Hv——水深，m

hf——泉涌高度，m

n——多变指数

Pa——大气压力，Pa

Qg——气体体积流量，m3/s

r——距离羽流中心的水平距离，m

v——羽流速度，m/s

vs——滑移速度，m/s

α——夹带系数

γ——动量放大系数

λ——宽度比

下角标

c——羽流轴线处的值

f——泉涌