李文欣 黄启玉 孙旭 张靖 孙超超

1中国石油大学（北京）油气管道输送安全国家工程实验室·石油工程教育部重点实验室

2中国石化工程建设有限公司

近年来，随着外浮顶储油罐不断朝大型化发展，由密封圈密封不严或失效导致的油气泄漏、起火爆炸事故逐年增加，不仅造成能源浪费，还会引起环境污染甚至人员伤亡。CHANG 和LIN[1]对国内外529 起油罐火灾事故进行统计表明，在雷击引起的浮顶油罐火灾事故中，因雷击密封圈引起的约占83%。因此，有必要对大型外浮顶储油罐密封圈油气扩散规律进行研究，以便正确预测油气分布规律、预估起火爆炸风险及范围，为外浮顶罐火灾事故的预防和救援提供理论支持。

在过去几十年里，人们开展了多种工业场景下的气体扩散规律研究。MAROTZKE、HUBER[2]、刘国梁[3]、吴晋湘[4]和施志荣[5]等开展了大量的风洞实验，研究可燃气体扩散规律和燃爆范围。其中MAROTZKE 等搭建平地、斜坡等实验场景，测取重气绕各种不同类型障碍物的扩散结果；HUBER等考虑了长方体建筑物对流场结构及气体扩散的影响，得出建筑物尾流区气体浓度分布的数学模型；刘国梁等建立直流式风洞实验系统，研究围墙和树对重气扩散特征的影响；吴晋湘等研究了液化石油气在有限空间内的扩散过程，实验结果表明，可燃区域随泄漏时间推移而逐渐向上方发展；施志荣运用示踪技术，在油气储运安全综合实验平台上模拟事故现场的有害气体的扩散规律。宋贤生[6]、王建[7]、赵刚[8]等采用雷诺时均（RANS）方程与湍流模型相结合的方式模拟流场，研究气体浓度分布和障碍物对扩散过程的影响。其中宋贤生等利用CFX软件模拟罐区油气扩散过程，研究群罐和防火堤对油气浓度分布的影响；王建等通过建立大型球罐区可燃气体泄漏扩散和燃爆的数值模型，确定可燃范围并对燃爆强度进行评估；赵刚等建立海上油气钻采平台的三维模型，分析平台结构、设备布局、风速对天然气浓度分布的影响。竺柏康[9]、赵晨露[10]、郝庆芳[11]等采用小型试验、数值模拟方法研究浮顶罐的油气泄漏问题。其中竺柏康等建立浮顶罐试验模型，发现风速比温度对密封圈油气浓度的不均匀分布影响更大；赵晨露、郝庆芳等采用RANS 方法与湍流模型相结合的方式模拟流场，研究不同风速、不同液位、不同浮盘泄漏孔隙位置的油气扩散规律。

已有储罐油气扩散数值研究中，均采用了RANS 与湍流模型相结合的求解方法。由于RANS方法仅能求解流场时均量，不能准确模拟流场旋涡结构和油气非定常耗散过程。为克服RANS 方法的缺点，采用大涡模拟（LES）方法，该方法可直接求解流场大尺度旋涡结构，更加准确描述油气随流场的输运和耗散过程。GOUSSEAU 等[12]对比了RANS 和LES 的计算精度，发现随着障碍物回流区对泄漏源影响的增大，RANS 模型对流通量的预测结果相较于LES 模型会产生较大误差。本文采用LES 方法，能更加准确描述大型外浮顶罐绕流流场及密封圈泄漏油气耗散过程。

1 控制方程及数值解法

根据大涡模拟的基本思想[13]，对不可压缩流体的N-S 方程作滤波处理，可得大涡模拟控制方程，即

式中：上划线为滤波后的场变量；下标i、j=1～3，为x、y、z三个方向的分量；ui、uj为气流速度分量，m/s；xi、xj为坐标分量，m；t为时间，s；ρ为气流密度，kg/m3；p为压强，Pa；ν为气流运动黏度，m2/s；τij为亚格子应力，m2/s2。

在控制方程中，和pˉ为待求量。由于式（3）中项不是和的显式表示，τij也是未知量，因此方程组（1）～（3）不封闭，需要通过亚格子尺度模型建立τij与、的联系。基于BOUSSINESQ 假设，亚格子应力可采用下式计算

式中：νSGS为亚格子涡黏系数，m2/s；为可解尺度应变率张量，s-1；τkk为亚格子应力的各向同性部分，m2/s2；δij为Kronecker 符号。

式（1）、（2）、（4）、（5）即为外浮顶罐绕流风场的控制方程。

在流过储罐的气流影响下，从密封圈泄漏出的油气将向周围环境扩散。由于油气组分与空气的浓度差较小，且环境温差不大，忽略浮力的影响，对不可压缩湍流中的组分输运方程进行过滤，可得油气扩散控制方程，即

式中：c为油气浓度，kg/m3；λ为质量耗散系数，m2/s；为亚格子组分输运过程，可通过建立与式（4）、（5）类似的亚格子模拟描述。

选用FLUENT 提供的DSM 模型（动态Smargorinsky-Lilly 模型）求解外浮顶罐油气扩散问题。DSM 模型以SM 模型（Smargorinsky-Lilly 模型）为基础，通过多次过滤，把湍流局部结构信息引入到亚格子应力中，可在计算过程中调整模型系数。该模型确定的亚格子涡扩散系数随流场和标量场的不同而发生变化，可以克服SM 模型中将湍流普朗特数视作常数的缺点。大涡模拟求解过程具体设置如下：控制方程离散采用三维双精度基于压力的分离式求解器；压力-速度耦合采用PISO 算法；空间离散中梯度项离散采用Least Squares Cell Based 格式，压力项离散采用Standard 格式，动量方程离散采用Bounded Central Differencing 格式，油气对流扩散项离散采用Second Order Upwind 格式，瞬态项求解采用Bounded Second Order Implicit 格式；压力项求解亚松弛因子设置为0.01；动量求解亚松弛因子设置为0.05；甲烷求解亚松弛因子设置为1；无量纲时间步长设置为Δt*=Δt×u∞/(1.5H)=0.062（Δt为时间步长，s；u∞为来流速度，m/s；H为储罐高度，m）；每个时间步长内的最大迭代数设置为200；连续性方程、油气浓度和x、y、z三个方向速度的残差收敛准则均为10-5。

2 大涡模拟解法验证

大型外浮顶罐的罐容大且储存介质易燃易爆，难以通过现场试验获取强风作用下的密封圈油气扩散数据。浮顶油罐的油气扩散研究归根结底是确定风场绕流障碍物时对油气的裹挟、运移规律，因此基于HUBER 试验，验证本文大涡模拟解法的准确性。根据HUBER 等的风洞实验建立数值计算模型，数值计算模型与实验模型的尺寸比为1∶1（图1）。在计算域为26H×14H×7H（H为长方体建筑物的高度，0.25 m）的流场中放置尺寸为H×H×2H的长方体建筑物（长边垂直于来流方向）。建筑物迎风面距流场上游进口5H，建筑物背风侧距流场下游出口20H，中心距流场两侧边界7H，建筑物顶部距流场顶部6H。紧贴建筑物后墙面有1 个高为1.5H、直径为0.042H的烟囱，烟囱顶部垂直向上排放泄漏气体（甲烷和空气的混合物，甲烷体积分数为1%），气体泄漏速度为We=1.5u∞。流场进口自由来流对应的雷诺数为Re=6.4×104（基于H）。边界条件设置如下：流场进口采用速度进口，流场出口采用自由出流边界条件，流场两侧采用对称边界条件，流场顶部采用滑移边界条件，烟囱顶部设置速度入口边界条件，流场底部、建筑物表面和烟囱侧面均采用无滑移边界条件。如图1 所示，在建筑物下游5H和10H处共设置6 条垂直于来流直线并检测其上浓度变化。其中，LINE1、LINE4 垂直于地面，LINE2、LINE3、LINE5 和LINE6 平行于地面。LINE2 和LINE5 距地面1.5H，LINE3 和LINE6 在地面上。

图1 气体扩散基准模型Fig.1 Benchmark model of gas diffusion

采用高质量六面体结构化网格对流场进行网格剖分，所得网格包括4 565 210 个节点，如图2 所示。根据HUBER 等的实验，各计算参数设置如下：参考压力为101 325 Pa；参考温度为300 K；启动组分输运模型，混合物为空气和甲烷（混合物之间不发生化学反应），空气的密度为1.225 kg/m3、动力黏度为1.789 4×10-5Pa·s，甲烷的密度为0.667 9 kg/m3、动力黏度为1.087×10-5Pa·s。浓度场求解稳定后，将t*（计算时间，无量纲）=600～900 的时均甲烷浓度无量纲化，与风洞实验数据进行对比，检验大涡模拟计算结果的准确性。定义浓度系数K为

式中：χ为甲烷浓度，kg/m3；A为建筑物的截面积，m2；u∞为对应参考高度处的自由来流速度，m/s；Qe为甲烷的泄漏速率，kg/s。

图2 气体扩散基准模型的计算网格Fis.2 Computational grid for gas diffusion benchmark model

如图3 所示，大涡模拟计算所得6 条检测线上的浓度系数分布与HUBER 等的实验结果吻合良好，验证了本文油气扩散模型和大涡模拟解法的准确性。

图3 大涡模拟计算的浓度系数与风洞实验对比Fig.3 Comparison of concentration coefficients calculated by LES and wind tunnel experiments

3 大型外浮顶罐密封圈油气扩散

3.1 计算模型

采用验证后的气体扩散大涡模拟解法模拟10×104m3大型外浮顶罐的油气扩散过程。大型外浮顶罐模型的流场计算域和密封圈泄漏区域划分如图4 所示，其中D为储罐直径，H为储罐的罐体高度（H=0.275D），H*为储罐最大安全液位（H*=0.25D），h为液位高度（即浮盘所在位置）。流场计算域范围为15D×5D×4H，采用结构化网格对计算域进行空间离散。流场进口距储罐中心5D，采用速度入口边界条件；出口位于储罐下游10D处，选用自由出流边界条件；两侧边界分别距离储罐中心2.5D，施加对称边界条件；顶部边界距离罐顶3H，采用滑移壁面边界条件；流场底部和储罐壁面选用无滑移边界条件。

图4 大型外浮顶罐模型的流场计算域和密封圈泄漏区域划分Fig.4 Flow field calculation domain and seal leakage area division of large external floating-roof tank model

外浮顶罐的密封圈在日常操作、风吹日晒等影响因素下，会出现不同程度的磨损和破坏，密封圈的油气泄漏位置和泄漏组分因罐而异，难以通过一个模型反映出所有储罐的泄漏场景。为反映不同的泄漏情况，本文的数值模型将环形密封圈均匀分为4 份，选用最易挥发的轻烃组分CH4作为泄漏物，对应设置4 种泄漏位置。如图4c、图4d 所示，上风侧为 -45°＜θ＜45° 的红色区域；侧风向为45°＜θ＜135° 的绿色区域和225°＜θ＜315° 的黄色区域；下风侧为135°＜θ＜225°的蓝色区域；整个环向为0°＜θ＜360°的整个环形密封圈区域（θ为泄漏点与迎风子午线的夹角）。储罐蒸发损耗与油品性质、风速、所在地的大气压及密封形式等因素有关，现场数据不易获取。因此，本文将某油田5×104m3原油外浮顶罐的年油气损耗量折算为泄漏速度[14]，设置泄漏区域的速度入口边界条件，以甲烷作为泄漏介质，假设其以u=0.03u∞的速度沿竖直方向从泄漏口逸出。对于现场实测数据已知的情况，可将实测泄漏速度与油气组分代入，采用与本文相同计算流程求解。来流雷诺数固定在Re=1.64 × 106（基于储罐直径D），计算液位h=25%H*、50%H*、75%H*、100%H*和泄漏位置为上风侧、侧风向、下风侧和整个环向等共16种工况。无量纲时间步长Δt*=Δt×u∞/D=0.037 5，待CH4浓度波动稳定后，将t*=30～70 的CH4体积分数进行时均统计。

3.2 环向油气浓度分布

图5 给出了不同情况下密封圈上方2.5%H*高度处的CH4体积分数分布。如图5a 所示，当油气从上风侧泄漏时，甲烷浓度沿环向缓慢升高，达到峰值后迅速减小；随液位升高，浓度峰值沿环向后移，但泄漏扩散范围变化不大。如图5b 所示，当油气从侧风向泄漏时，甲烷浓度沿两端逐渐减小，且45°＜θ＜135°范围内的上风侧甲烷浓度略高于下风侧。液位高度为h=25%H*和h=50%H*时，在泄漏口上风侧0°＜θ＜45°范围内充满了较低浓度的甲烷，而液位高度为h=75%H*和h=100%H*时上风侧的甲烷浓度几乎为0。说明低液位时泄漏气体在浮盘上方空间的扩散更显著。如图5c 所示，当油气从下风侧泄漏时，泄漏的甲烷聚集在泄漏口附近，不会向上游扩散；随液位升高，浓度零值沿环向略有前移。如图5d 所示，当油气从整个环向泄漏时，液位较低时，上风侧甲烷浓度高于下风侧，浓度最高值位于θ=45°附近；满液位时，甲烷气体沿整个环向分布均匀。

为理解泄漏气体在环向的分布规律，图6 给出了不同液位高度的浮盘表面瞬态流线图。如图6a、图6b 和图6c 所示，储罐液位较低时，浮盘上方的大部分气流撞击储罐内壁的后侧，沿反方向运动，使下风侧的油气向上风侧输运，上风侧泄漏的油气仍聚集于此，而不进入其他区域。当油气泄漏出现在侧风向和下风侧时，油气则会沿流线向上风侧扩散，小部分气流撞击储罐内壁前侧再次产生回流，沿环向向两侧流动，将油气向两侧输运，因此，迎风子午线处（θ=0°）的油气浓度不是最高值。这两股气流交汇于θ=45° 附近并产生旋涡，因此，当油气从整个环向均匀泄漏时，油气浓度最高值约位于θ=45° 。如图6d 所示，储罐接近满液位时，气流产生多股分流，在浮盘上方形成多处旋涡，因而此时的环向油气浓度分布较均匀。

图5 不同液位高度下泄漏口上方环向油气浓度分布（ t*=30～70）Fig.5 Annular oil vapor concentration distribution above the leakage port at different liquid level heights（t*=30～70）

图6 不同液位高度下浮盘表面瞬态流线（ t*=70）Fig.6 Transient flow line of floating roof surface at different liquid level heights（t*=70）

3.3 流向油气浓度分布

图7、图8 分别给出了环向泄漏时z=0 剖面上的CH4体积分数分布和瞬态流场。如图7a、图8a 所示，当储罐位于低液位h=25%H*时，在浮盘上方形成一个大尺度旋涡，涡量最大。在该旋涡的输运作用下，从密封圈泄漏出的油气可以充满整个空腔。常用涡量来描述旋涡运动的大小和方向，流场的中央量为流体速度的旋度，ω*为无量纲涡量。如图7b、图7c、图8b、图8c 所示，随着储罐液位从h=25%H*升高到h=50%H*和h=75%H*，空腔大尺度旋涡分裂成多个小尺度旋涡，涡量显著降低，旋涡对密封圈泄漏油气的输运作用减弱，因此，空腔内的油气扩散范围逐渐减小，油气主要聚集在密封圈附近。如图7d、图8d 所示，当储罐液位为h=100%H*时，浮盘上方空腔不再被旋涡充满，仅有浮盘上风侧被旋涡覆盖，由密封圈泄漏的油气极易被输运到主流中，因此，在浮盘表面只分布有少量油气。

3.4 可燃区分析

甲烷的爆炸极限为5%～15%（体积分数），低于爆炸下限时不燃烧，高于爆炸上限后安静燃烧。基于浓度计算结果，绘制爆炸下限CL=0.05 和爆炸上限CH=0.15 对应的时均浓度等值面图，进一步分析储罐浮盘上方空腔内的三维油气分布。如图9 所示，离泄漏口较近的红色曲面为爆炸上限等值面，包裹红色曲面的蓝色曲面为爆炸下限等值面。对爆炸下限的等值面进行面积积分，并除以相应的密封圈泄漏面积，可得单位泄漏面积对应的可燃区，用SL（无量纲）表示。如图9、图10 所示，低液位时（h=25%H*），浮盘上方大尺度旋涡的输运作用使可燃区的范围最大；当泄漏口位于下风侧时，SL最高可达55。随着液位增加，空腔内旋涡输运强度减弱，可燃区不断减小。满液位时（h=100%H*）的可燃区面积远小于其他三个液位，当泄漏口位于下风侧时SL达到最低值15。从图9、图10 还可以看出，当液位较低时，可燃区的面积差异较大，从下风侧泄漏出来的油气对应的SL最大，侧风向次之，上风侧最小。储罐接近满液位时，可燃区的面积大致相同，从侧风向泄漏出来的油气对应的SL最大，上风侧次之，下风侧最小。

图7 不同液位高度下油气瞬时浓度场（ t*=70）Fig.7 Instantaneous concentration fields of oil vapor at different liquid level heights（t*=70）

图8 不同液位高度下储罐绕流的瞬时涡量场和流线（t*=70）Fig.8 Instantaneous vorticity field and flow line of the flow around the rank at different liquid level heights（t*=70）

图9 密封圈泄漏的甲烷时均浓度等值面（CL=0.05，CH=0.15）Fig.9 Isosurface of the mean concentration of methane leaked from the seal（CL=0.05，CH=0.15）

图10 单位泄漏面积对应的可燃区Fig.10 Combustible area per unit leakage area

4 结论

利用大涡模拟高精度计算方法模拟10×104m3大型外浮顶罐密封圈泄漏油气的扩散过程，得到不同液位储罐瞬态绕流风场的旋涡分布和演变过程，重点分析泄漏位置和液位对罐顶油气分布和可燃区的影响。模拟结果如下：

（1）当液位较低时，气流在储罐浮盘上方空腔内形成大尺度旋涡，在旋涡的输运作用下，泄漏油气从下风侧向上风侧扩散，上风侧油气浓度高于下风侧；随着液位升高，浮盘上方空腔内大尺度旋涡分裂为小尺度旋涡，输运作用降低，油气聚集在泄漏口附近。

（2）当油气沿整个环向密封圈均匀泄漏时，油气浓度最高值位于迎风面左右两侧45°附近，应及时维护和检修这一区域的可燃液体监测装置。

（3）随着液位升高，泄漏油气在浮盘上的环向扩散作用减弱，主要聚集在泄漏口附近，可燃区减小；满液位时，浮盘只有在靠近迎风面的部分区域被旋涡覆盖，同时气流将部分泄漏油气带离浮盘，可燃区达到最低值。

（4）当油气分别从上风侧、侧风向和下风侧泄漏时，低液位储罐的可燃区范围由大到小依次为下风侧、侧风向、上风侧，满液位储罐的可燃区范围由大到小依次为侧风向、上风侧、下风侧，低液位储罐的可燃区范围比高液位储罐高1.3～1.7 倍。