地面输油管道泄漏流散数值模拟*

史晓蒙，吕宇玲，2，3，杨玉婷

(1.中国石油大学(华东) 储运与建筑工程学院，山东 青岛 266580;2.山东省油气储运安全省级重点实验室，山东 青岛 266580;3.青岛市环海油气储运技术重点实验室，山东 青岛 266580)

0　引言

随着大型油库数量越来越多，储罐和输油管道受到腐蚀、高压及外力破坏等作用，库区泄漏事故时有发生[1-2]。由于这种泄漏事故泄漏量大，波及范围广及原油流散行为的不可预见性，易造成巨大的人员及财产损失。因此研究外泄原油的流散规律对于库区安全有重要意义。

国内外关于原油泄漏后果的研究有很多。刘瑞凯[3]利用数值模拟研究了海底输油管道泄漏原油在海泥中的分布。姚志强[4]利用数值模拟确定埋地原油管道发生泄漏时，不同时刻油品在土壤中的扩散范围。李晰睿[5]以“11·22”中石化东黄输油管道泄漏爆炸事故为案例，研究了暗渠内输油管道泄漏、油品蒸发和爆炸过程。国内外学者关于液池扩展面积的计算方法及扩展状态进行了研究。William[6]和Cronin[7]利用大型实验研究了储罐发生破裂时液体对防火堤的冲击载荷及液池的扩展过程。Brambilla[8]提出了一种应用几何计算获得不同阻挡形式下液体扩散最大距离的方法。项小强等[9]采用经验公式计算液池扩展半径。王超等[10]提出了一种计算连续泄漏、泄漏停止时液池半径的半经验关系式，但由于忽略了液体粘度的影响，该方程不适用于高粘度液体在地面上的蔓延。李大全等[11]建立了一种可以评价泄漏事故危害程度的理论分析模型。潘旭海等[12-13]根据质量守恒定律，建立了动态液池蒸发模型，并用其对苯泄漏事故进行了模拟分析。

上述计算模型大都假设液池为圆柱体，扩展是轴对称的，这与实际情况不符，因为液体的流散由于受到障碍物阻挡、地形起伏等原因会呈现不规则的形状。Hirt[14]提出了用VOF方法来模拟动态的自由液面流动，得到不同时刻的相分布图，并与实验相吻合。以下主要通过建立库区地面输油管道泄漏扩散过程的三维CFD模型，应用FLUENT软件中的VOF模型进行模拟，分析油品泄漏后扩散趋势，研究结果有助于分析泄漏事故后果影响范围，合理制定防护及救援措施。

1　泄漏扩散模型

1.1　数学模型

由于流动主要属于自由泄流，近似为不可压缩流体。在空间直角坐标系下采用标准κ-ε湍流模型进行三维数值模拟，控制方程如下。

连续性方程：

(1)

动量方程：

(2)

ε方程：

(3)

κ方程：

(4)

式中：ρ为密度；U为速度矢量；ui为流体沿xi方向的速度分量(xi在三维情况下为x，y，z方向)；η为流体的动力粘度；σε，σκ分别为湍动能和湍动能耗散率所对应的Prandlf数；C1，C2为经验常数；t为时间。

VOF方法的原理是：在欧拉网格系统上定义一个体积分数函数f，每个网格上的f值由网格内含指定项的体积量来确定，由下列方程确定：

(5)

1.2　模型验证

选取Cronin & Evans[7]2001年进行的液体在陆地上扩散实验作为模型验证的依据。实验条件为：1个下端有25 mm裂口，高为1.8 m的水槽来模拟直径为70 m的储罐破裂，在混凝土表面设置多个阻力探测器监测记录水波到达时间及传播速度。基于以上实验，建立模型，以水作为计算流体，使用VOF模型和标准κ-ε湍流模型，模拟不同地面粗糙度情况下水在地面上的扩展，水泥混凝土面层粗糙度不大于1.0 mm，因此设置粗糙度分别为0，0.1，0.2，0.8 mm，计算结果与实验结果进行比较，如图1、图2所示。

图1　不同粗糙度下流散时间随距离变化Fig.1　Dispersion time varies with distance under different roughness

图2　水波沿中轴线传播速度随距离变化Fig.2　Frontal speed along a radial line varies with distance

当地面粗糙度为0.1 mm时，水波前锋位置随时间变化以及距泄漏口不同距离的水波传播速度与实验数据吻合地较好。因此，采用上文提出的液体流散模型是适用的。

1.3　几何模型与网格划分

根据石油库区消防安全要求，相邻罐组防火堤脚线之间宽度不小于7 m，故选择计算区域为60 m×10 m×3 m，内径为712 mm，考虑到管道外的保温层厚度，管道离地面高度为0.2 m[15]，在中心处有直径不同的圆形开孔，模拟接近地面敷设的站内管道发生破裂，模型几何尺寸如图3所示。网格采用四面体网格、六面体网格和结构化网格非均匀划分方法，对泄漏孔口处和地面进行网格加密。

图3　模型几何结构尺寸Fig.3　Geometric dimensions of the model

设置管道两端为压力入口，油库内输油管道运行压力一般在0.2～0.4 MPa之间，模拟管道压力分别为0.2，0.25，0.3，0.35，0.4 MPa。设置大气出口为压力出口，地面粗糙度为0.000 1 m，油品密度为846.47 kg/m3，粘度为3.458×10-6m2/s。

为了说明网格密度对数值计算结构的影响，对模型

进行疏密2种形式的网格划分(网格单元数分别为1 560 000和3 300 000)，计算结果相近。计算表明，网格具有一定的无关性，本文使用网格单元数为1 560 000的网格进行计算。

2　油品泄漏规律分析

2.1　管道运行压力不同的情况下的油品扩散

当泄漏孔口直径d=0.5 m，管道运行压力在0.2～0.4 MPa时，油品从破裂口泄漏，即薄壁孔口自由出流，由薄壁圆形孔口泄流公式[16]计算的5种管道运行压力下泄流喷射速度分别为：21.08，23.57，25.82，28.22，29.82 m/s，本文CFD模拟的速度分别为21.90，23.04，26.76，28.82，29.74 m/s，二者误差均不超过4%，再次验证此模型可以用于模拟管道破裂油品泄漏。

模拟得到油品相分布如图4所示，图中空白区域表示油品体积分数小于0.05%，代表油品未扩散到的地面，灰色区域表示油品扩散到的地面。

图4　压力不同时管道泄漏油品扩散区域变化Fig.4　Affected area under different operating pressure of the pipeline

泄漏油品流散距离随时间变化关系如图5所示。流散距离随时间变化过程可以分为2个阶段，第一阶段是油品从孔口泄漏，形成紊动射流，油品的扩散受到孔口射流的冲击影响，流散距离随着时间呈指数增长，由于受到地面摩擦阻力的作用，流散速度逐渐减小。第二个阶段认为射流过程结束，不再考虑射流作用，同时横向扩散受到阻挡，液体在地面摩擦力和液体表面张力的作用下扩展，扩展速度保持不变。流散距离随时间变化呈现线性函数增长趋势，其曲线斜率表示油品在地面扩展的速度。随着压力的增大，扩展速度增大。可以将流散0～5 s作为第一阶段，5 s之后作为第二阶段。

2.2　管道泄漏口面积不同时油品扩散

不同泄漏孔直径油品扩散区域变化如图6所示，管道泄漏面积的变化对油品的横向扩散影响较为明显。当泄漏口直径d=0.3 m时，泄漏油品从泄漏口呈直线喷射出，横向扩展较小，液体集中在射流轴线上向前推进，在短时间内流散距离最远。当d=0.5 m时，泄漏口附近油品呈扇形喷射，在中心线上的流散前锋也在迅速向前推进的同时，液体也发生横向扩展。

油品扩散距离随不同管道破裂孔径的变化如图7所示，当d=0.5 m时，由于泄漏口附近液体横向速度较大，油品横向扩展最远，液体厚度变薄，向前扩展受到影响，流散距离较小。当流散不受到射流的影响时，油品扩展速度随泄漏孔直径增大而增大。

图5　不同管道压力下泄漏流散距离与时间的变化关系Fig.5　Spread distance vs. spread time under different operating pressure of the pipeline

2.3　扩展速率与泄漏量的关系

由2.2节和2.1节分析得到，当液体扩散不受射流影响扩展速度稳定时，扩展速度随压力增大而增大，随泄漏孔径增大而增大，扩展速度跟泄漏流量有一定关系，如图8所示。随着泄漏流量增大，扩展速度增大，扩展速度可以表示为关于泄漏流量的函数：

v=1.869 61+0.084 54e0.493 28Q

(6)

式中：Q为泄漏流量，m3/s。

图6　不同泄漏孔直径油品扩散区域变化Fig.6　Affected area under different leakage hole diameter

图7　不同泄漏口直径流散距离随泄漏时间的变化关系Fig.7　Spread distance vs. spread time under different leakage hole diameter

图8　扩展速度与泄漏流量的关系Fig.8　Spread velocity vs. leakage flow rate

3　流散面积模型

假定在此种泄漏流散环境下，油品流散面积随管道压力、泄漏孔径和流散时间变化的关系为S=f(p,d,t)。

3.1　流散面积随管道压力的变化

图9　8 s时流散面积随管道压力的变化Fig.9　Affected area vs. operating pressure at t=8 s

(7)

式中：d为泄漏孔直径，m。

图10　压力相关系数随泄漏孔径的变化Fig.10　Coefficient of operating pressure vs. the leakage hole diameter

3.2　流散面积随破裂孔径的变化

图11是各个管道压力下泄漏8 s时流散面积随泄漏孔径的变化关系，可以回归成流散面积S关于泄漏孔径d的一元三次方程：

y=A+Bx+Cx2+Dx3

(8)

(9)

式中:p为管道压力，MPa。

图11　8 s时流散面积随泄漏孔径的变化Fig.11　Affected area vs. the leakage hole diameter at t=8 s

图12　泄漏孔径系数与管道压力的变化Fig.12　Coefficient of leakage hole diameter vs. the operating pressure

3.3　流散面积随泄漏时间的变化

(10)

图13　泄漏孔径为0.5 m时流散面积随时间的变化Fig.13　The affected area vs. time at d=0.5 m

图14　时间相关系数与管道压力的关系Fig.14　Coefficient of time vs. the operating pressure

3.4　扩展速率方程式和流散面积模型的验证

假设工况1：管道压力0.32 MPa，泄漏孔直径为0.4 m；假设工况2：管道压力为0.3 MPa，泄漏孔直径为0.38 m。使用得到的偏微分方程和数值模拟对2种假设工况分别进行计算，并对结果进行比较。以式(6)得到油品稳定扩展速率分别为2.31，2.22 m/s，利用数值模拟得到油品稳定扩展速率分别为2.39，2.19 m/s；使用式(7)，(8)，(9)得到的油品扩散面积分别为189.46，201.67 m2，利用数值模拟得到油品扩散面积分别为187.35，209.58 m2。二者相差不大，说明上文扩展速率方程式和流散面积模型在一定范围内具有适用性。

4　结论

1)若管道发生泄漏，在宽度为10 m的地面上扩散，当扩散速度稳定时，稳定的扩展速度随泄漏流量增大而增大，扩展速度与泄漏量的具体关系为：

v=1.869 61+0.084 54e0.493 28Q

(11)

2)在宽度为10 m的流散区域内，管道压力为0.2～0.4 MPa，泄漏孔径为0.3～0.5 m时，表示流散面积S与管道压力p、泄漏孔径d、泄漏时间t关系的偏微分方程组为：

(12)

当泄漏孔直径相同时，流散面积随管道压力增大呈线性函数增大，压力相关系数与泄漏孔径呈指数函数关系；当管道压力相同时，流散面积随泄漏孔径呈三次函数变化，先减小后增大，泄漏孔径相关系数与管道压力呈指数函数关系；在泄漏的开始阶段，流散面积随泄漏时间的变化随管道压力呈线性增长。

[1]Chang J I, Lin C C. A study of storage tank accidents [J]. Journal of Loss Prevention in the Process Industries, 2006, 19(1):51-59.

[2]潘旭海, 蒋军成. 重(特)大泄露事故统计分析及事故模式研究[J]. 化学工业与工程, 2002, 19(3):248-264.

PAN Xuhai, JIANG Juncheng. Analysis on important release accidents and modes studying [J].Chemical Industry and Engineering, 2002, 19(3):248-264.

[3]刘瑞凯，吴明，王同秀，等. 海底埋地热油管道泄漏扩散的数值模拟[J]. 中国安全生产科学技术, 2012, 8(8):63-68.

LIU Ruikai, WU Ming, WANG Tongxiu, et al. Numerical simulation on leakage and diffusion of submarine buried hot oil pipeline[J]. Journal of Safety Science and Technology, 2012, 8(8):63-68.

[4]姚志强,埋地原油管道油品泄漏扩散三维数值模拟[J]. 中国安全生产科学技术, 2015, 11(6):132-136.

YAO Zhiqiang，Three-dimensional numerical simulation on leakage and diffusion of oil in buried crude oil pipeline[J]. Journal of Safety Science And Technology, 2015,11(6):132-136.

[5]李晰睿. 暗渠下输油管道泄漏爆炸分析及事故后果研究[D]. 成都：西南石油大学, 2015.

[6]An experimental investigation of bund wall overtopping and dynamic pressures on the bund wall following catastrophic failure of a storage vessel[M]. HSE Books, 2005:4-42.

[7]CLEAVER R P, CRONIN P S, EVANS J A, et al. An experimental study of spreading liquid pools[C]//Proceedings of the Hazards XVI. Manchester, UK, 2001：167-179.

[8]Brambilla S, Manca D. On pool spreading around tanks: Geometrical considerations [J]. Journal of Hazardous Materials, 2008, 158(1): 88-99.

[9]项小强, 戴联双, 曹涛，等. 输油管道泄漏液池的蒸气云扩散模型[J]. 油气储运, 2011, 30(5): 334-336.

XIANG Xiaoqiang, DAI Lianshuang, CAO Tao，et al. Vapor cloud diffusion model for leakage pool of oil pipeline [J]. Oil and Gas Storage and Transportation, 2011, 30(5): 334-336.

[10]王超. 液体泄漏形成液池扩展面积的计算方法综述[J]. 安全与环境工程, 2012, 19(6):125-128.

WANG Chao. Review of area expansion calculation method for liquid leakage to liquid pool [J]. Safety and Environmental Engineering, 2012，19(6)：125-128.

[11]李大全, 姚安林. 成品油管道泄漏扩散规律分析[J]. 油气储运, 2006, 25(8):18-24.

LI Daquan, YAO Anlin. Analysis of leakage and diffusion law of oil product pipeline [J]. Oil and Gas Storage and Transportation, 2006, 25(8):18-24.

[12]潘旭海, 蒋军成. 泄漏液池动态蒸发过程模型研究与分析[J]. 工业安全与环保, 2006, 32(4):43-46.

PAN Xuhai, JIANG Juncheng. Research and analysis on model for dynamic evaporation process of pool [J]. Industrial Safety and Environmental Protection,2006, 32(4):43-46.

[13]潘旭海, 蒋军成. 危险品罐区突发性泄漏事故实时环境风险模拟分析[J]. 环境科学学报, 2004, 24(3):539-544.

PAN Xuhai, JIANG Juncheng. Real-time environment risk analysis for accident release of hazardous materials around tank area [J]. Acta Scientiae Circumstantiae, 2004, 24(3):539-544.

[14]Hirt C W, Nichols B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201-225.

[15]樊宝德, 朱焕勤. 油库设计规范[M].北京：中国石化出版社, 2007: 109-115.

[16]倪玲英. 工程流体力学[M]. 东营：中国石油大学出版社，2012：212-215.