引言

油品在储运过程中不可避免地会由于设施建材的质量问题、施工隐患、油品对设备的腐蚀穿孔及外力损伤等因素，致使装置破裂而引发泄漏事故。泄漏的油品在地面形成液池并向周围蔓延，挥发具有毒性的油气混合气体，极易导致发生火灾、爆炸等次生事故。如2010年大连石化公司的储罐区发生多次油品泄漏事故，溢油在地面形成液池后发生燃烧、爆炸等次生灾害，严重危害周围人员的生命安全。因此研究油品泄漏后液池的形成及扩展特征具有重要意义。

在实验研究方面，Atherton^[1]和Cleaver等^[2]利用实验模型研究了储罐破裂后泄漏的液体在地面的覆盖距离以及对围堰产生的冲击压力；Brambilla等^[3]提出了一种应用几何计算的方法，得到在不同围堰阻挡形式下液体扩散的距离，与相关实验部分结果吻合较好。在理论研究方面，Fay^[4-5]在总结前人进展的基础上以平静水面为背景开创性地提出液池扩展三阶段模型，认为油品在自身重力、黏性力、惯性力和表面张力的共同作用下向周围扩展。之后众多学者基于该模型进行研究并改进，牟林等^[6]对模型中第二阶段的黏性扩展公式进行了修正；Lehr等^[7]分析了风场等因素对液池非对称性扩展的影响，使液池扩展模型更加符合实际情况。

综上所述，已有的实验方法局限于难以还原实际尺寸设施的泄漏过程，理论方法则通过多种假设条件推导，忽略了某些关键因素的影响。近年来随着计算机技术的快速发展，数值模拟成为众多学者的选择，特别是采用计算流体力学(CFD)方法构建实际尺寸模型，结合环境影响、流体流动特性、流体间作用力等因素，使计算结果更贴合实际。李泽等^[8]建立了埋地管道的三维模型，模拟管道发生泄漏后汽油、柴油及原油在地下及地表的扩散范围。马跃等^[9]构建土壤的多孔介质模型并进行研究，得出埋地管道发生泄漏后油品在土壤中的渗透、扩散范围主要受泄漏孔位置的影响。在多相流模型的选择方面，Hirt等^[10]通过实验获得不同时刻下的液相分布，验证了流体体积(volume of fluid，VOF)模型在用于研究动态自由液面流动时的良好应用性。詹翔等^[11]基于VOF模型模拟偏二甲肼泄漏，得到不同时刻偏二甲肼沿壁面扩展及在地面形成液池的分布特征。史晓蒙等^[12]采用VOF模型对管道泄漏进行数值模拟，得到了油品的扩展速度与泄漏流量之间的关系式，以及管道压力、泄漏孔径等因素对流散面积的影响结果。

现阶段利用CFD对油品泄漏的研究多侧重于液池在水面的扩展或地下土壤的渗透，而关于地面上油品泄漏的研究较少。本文基于CFD方法建立油品泄漏模型，模拟油品泄漏后在地面形成液池的过程，进一步探究地面粗糙程度、泄漏口尺寸、泄漏速度、油品物性等因素对液池面积的影响规律，并根据液池扩展理论和模拟结果拟合得到计算液池面积的关联式，为确定事故影响范围及采取相应防护措施提供理论基础和技术指导。

1 数值模拟方法 1.1 几何模型及网格划分

以某输油管道及周围环境建立三维几何模型，高1 m，底边长8 m，管道直径0.5 m；中心处的管壁右侧有一圆形泄漏口，直径0.1 m，假定油品从泄漏口流出并在地面形成液池。由于模型内含有开孔、管道等复杂结构，采用非结构性网格建模，并对泄漏口和地面区域进行局部网格加密，如图 1所示。划分不同数量级的网格进行计算，可知不同泄漏时间下液池的面积随网格数的增加而增大，当网格数超过250万后变化幅度较小，因此经网格无关性验证，本文模型的网格数量定为2 470 177。

1.2 数学模型

本文基于VOF模型，在计算油品的泄漏和扩展过程中涉及到气相和液相，其质量守恒方程表示为

$ \frac{\partial }{{\partial t}}\left( {{\alpha _q}{\rho _q}} \right) + \nabla \left( {{\alpha _q}{\rho _q}{\mathit{\boldsymbol{u}}_q}} \right) = 0 $

(1)

动量守恒方程表示为

$ \begin{array}{l} \frac{\partial }{{\partial t}}(\rho \mathit{\boldsymbol{u}}) + \nabla (\rho \mathit{\boldsymbol{uu}}) = - \nabla \mathit{\boldsymbol{p}} + \nabla \left[ {\mu \left( {\nabla \mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{u}}^{\rm{T}}}} \right)} \right] + \\ \rho \mathit{\boldsymbol{g}} + \mathit{\boldsymbol{F}} \end{array} $

(2)

对于气液两相，VOF模型中涉及到的物理参数由各分相决定

$ {\rho = \sum {{\alpha _q}} {\rho _q}} $

(3)

$ {\mu = \sum {{\alpha _q}} {\mu _q}} $

(4)

在每一控制体内需满足条件

$ \sum\limits_{q = 1}^n {{\alpha _q}} = 1 $

(5)

式中，ρ为密度，kg/m³；u为流体速度，m/s；μ为流体黏度，Pa·s；p为压力，Pa；F为表面张力，在Fluent中基于连续表面力(CSF)模型计算；α_q为q相的体积分数，q表示液相或气相。

由于本文模拟的油品流速较低且黏度较大，因此选用更适合模拟低雷诺数下流动的SST k-ω湍流模型，求解方程如下。

k方程

$ \begin{array}{l} \frac{\partial }{{\partial t}}(\rho k) + \frac{\partial }{{\partial {x_j}}}\left( {\rho k{\mathit{\boldsymbol{u}}_j}} \right) = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + \\ {P_k} - {\beta ^\prime }\rho k\omega + {P_{kb}} \end{array} $

(6)

ω方程

$ \begin{array}{l} \frac{\partial }{{\partial t}}(\rho \omega ) + \frac{\partial }{{\partial {x_j}}}\left( {\rho \omega {\mathit{\boldsymbol{u}}_j}} \right) = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\omega }}}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + \\ \alpha \frac{\omega }{k}{P_k} - \beta \rho {\omega ^2} + {P_{\omega b}} \end{array} $

(7)

式中，k为湍动能；ω为湍动能耗散率；P_k为湍动能生成项；P_kb为湍动能源项；P_ωb为湍动耗散率源项；其他参数采用默认值：β′=0.09，α=5/9，β=0.075，σ_k=1.176和σ_ω=2。

1.3 边界条件及求解方法

模型底面和管道壁面设置为wall，代表地面和管道的管壁，通过改变地面的粗糙常数和粗糙高度来模拟不同地面条件，粗糙常数取值为0.5~1(表示地面粗糙的均匀度)，粗糙高度取0.001~0.30 m(表示地面凸出物的平均高度)；泄漏口为速度入口，速度范围0.25~1 m/s；其余为压力出口。采用基于压力求解器的非稳态PISO算法进行求解。

2 结果分析与讨论 2.1 液池的形成特征及结果对比 2.1.1 液池的形成过程

地面粗糙高度取0.010 m，粗糙系数为0.5，速度入口为1 m/s，模拟汽油从泄漏口流到地面形成液池的过程，如图 2所示。1~5 s内由于油品具有初始动能，且泄漏口距地面有一定高度，泄漏的油品碰撞地面后呈放射状向周围扩展；10 s后泄漏的油品量随时间增加，各溢油点逐步形成一个液池整体，由于地面无障碍物阻挡，液池在自身重力等的作用下近似呈圆形向周围扩展，液池面积随时间线性增加，直至铺满整个计算区域。

2.1.2 模型的分析验证

对于第1节所述泄漏模型，采用理论分析计算法及工程经验推导法进行比较验证。理论分析基于质量守恒定律，一定时间内泄漏口流出的油品质量与地面上液池内油品的质量相等，即

$ m = \frac{{\rm{ \mathsf{ π} }}}{4}{d^2} \cdot u \cdot t \cdot \rho = \frac{{\rm{ \mathsf{ π} }}}{4}{D^2} \cdot h \cdot \rho $

(8)

式中，m为油品质量, kg；d为泄漏口直径, m；D为液池直径, m；u为泄漏速度, m/s；ρ为汽油密度，取720 kg/m³；h为液池厚度, m。由于地面性质难以精确表示，液池厚度采用式(9)计算。

$ h = \max (\sqrt {\frac{{2\sigma (1 - \cos \theta )}}{{\rho g}}} , \varepsilon ) $

(9)

式中，σ为油品表面张力, N/m；θ为接触角，°；ε为地面粗糙高度，m。由于泄漏区域多为混凝土或土壤地面，粗糙高度较大，经计算比较后取较大值，即h=ε=0.010 m。

经验推导法是基于相关工程试验，以泄漏点为中心认为油品在光滑平面上呈扁的圆柱形向周围扩展^[13]。对于连续泄漏过程，液池的直径为

$ D = 2{\left( {\frac{{32gm{t^3}}}{{{\rm{ \mathsf{ π} }}p}}} \right)^{\frac{1}{4}}} $

(10)

图 3为理论分析、经验推导和CFD模拟3种方法下液池直径随时间的变化规律。由图可知模拟值和理论值的增长趋势相同，且模拟值比理论值高约9%，这是由于CFD模型中泄漏的油品具有初速度和势能，泄漏到地面后液池的扩展趋势较大，因此导致模拟值稍大于理论值。此外经验推导值和CFD模拟值相比，在前20 s内两者相差在±0%以内，之后两者的差距逐渐增大。这是由于经验推导法假设油品在光滑平面上扩展，忽略了实际情况下油品黏度和地面粗糙度的影响，导致液池面积增长率较模拟值高。以上计算结果表明本文所建立的油品泄漏CFD模型考虑了地面粗糙因素及油品自身物性，更符合油品的实际流动特征，因此具有良好的可行性。

2.2 油品泄漏扩散过程的影响因素分析 2.2.1 地面粗糙高度和粗糙常数

图 4为地面粗糙高度对液池面积的影响。可以看出液池面积随粗糙高度的减小而增大，且地面粗糙高度小于0.005 m时，液池面积基本不再变化。这是由于地面粗糙高度增大，即相当于地面上砂粒的高度增大，在较大砂砾的阻碍下，油品向周围地面流动时所受到的摩擦阻力上升，导致液膜边缘处的扩展速率下降，液池厚度增加，液池面积减小。

图 5为地面粗糙常数对液池面积的影响。由图可知液池面积随粗糙常数增大而略有下降，幅度在±5%以内，影响较小。粗糙常数可理解为地面凸出物(如砂砾)的非均匀性分布，在砂砾高度一定的情况下，砂砾的非均匀性分布对流动阻力的影响小于砂砾高度的影响；且当大量油品溢流到地面后，由于流体自身的重力和流动性，油品会迅速覆盖地面，降低了砂砾非均匀性分布的影响，导致粗糙常数对液池面积的影响较小。

2.2.2 泄漏速度和泄漏口直径

从图 6和图 7可以看出，液池面积随油品泄漏速度、泄漏口直径的增大而增大。在管道、储罐等设施的破损程度范围内(如管道的最大破损程度即为管道完全断裂)，当不考虑地面围堰等阻碍时，随泄漏口直径、泄漏速度的增大，一定时间内泄漏的油品总量也增多，在地面形成的液池面积也将增大；同时根据式(8)也可得液池的面积随泄漏速度、泄漏口直径的增大呈线性增长的趋势。

2.2.3 油品的密度和黏度

模拟常见的油品(表 1)发生泄漏，由图 8可知，汽油、柴油与燃油发生泄漏后形成的液池面积基本相等，相差在±4%以内，差距较小。本文模拟较短时间内油品的连续泄漏，根据Fay液池扩展理论^[5]，此时液池主要处于重力-惯性力扩展阶段，黏性力作用较小，因此在油品密度相差较小的情况下所形成的液池面积也基本相等，即所模拟油品的密度和黏度对短时间内泄漏形成的液池面积的影响较小。

2.3 液池面积的经验关联式

通过对以上各影响因素进行模拟分析可得，不同黏度、不同密度的油品泄漏形成的液池面积相差在±4%以内，不同地面粗糙常数下液池的面积相差在±5%以内，因此推测油品的密度和黏度、地面粗糙常数对液池面积的影响较小，在本节计算关联式的建立中不作考虑。另外增大地面粗糙高度会导致液池的面积减小，增大泄漏口直径或泄漏速度将导致面积增大，因此推测液池面积S与地面粗糙高度ε、泄漏口直径d、泄漏速度u有关，即S=f(d, u, ε, t)。将模拟所得数据进行拟合，得到液池面积的计算关联式为

$ S = 0.701 \cdot {\left( {\frac{d}{{{d_0}}}} \right)^{1.45}} \cdot \left( {\frac{u}{{{u_0}}}} \right) \cdot \left( {\frac{{{t^{1.1}}}}{{1 + \varepsilon }}} \right) $

(11)

式中，S为液池面积, m²；d为泄漏口直径, m，d₀为0.1 m；u为泄漏速度, m/s，u₀为1 m/s；ε为地面粗糙高度, m。

基于Webber^[14]对液池扩展过程的推导，得到液池半径关于时间的函数为R∝t^a，a为常数^[15]，与本文关联式S∝t^1.1表达形式相近。将液池面积的关联式计算值与CFD模拟值作于图 9中比较，可以看出两者相差在±5%以内，差值较小，表明拟合的关联式可用于预估油品泄漏后在地面形成的液池面积。

3 结论

(1) 与理论模型和经验模型相比，考虑了地面粗糙因素与油品物性的CFD模型更符合油品的实际流动过程。

(2) 在连续泄漏过程中，液池面积随泄漏速度、泄漏口直径的增大而增大；随地面粗糙高度的增大而减小，当地面粗糙高度小于0.005 m时，液池面积基本不再变化。地面粗糙常数、油品的密度和黏度对液池面积的影响较小。

(3) 根据CFD模拟结果拟合得到液池面积的计算关联式，计算值与模拟值误差在±15%以内。