基于两相流的水封油库油气泄漏运移规律及控制措施

2024-02-21简回香王存利蒋中明

长江科学院院报 2024年2期

唐栋,简回香,王存利,李毅,蒋中明

(1.长沙理工大学水利与环境工程学院,长沙 410114; 2.长沙理工大学水沙科学与水灾害防治湖南省重点实验室,长沙 410114)

0 引言

地下水封洞库因其经济、安全等优势成为了国家战略石油储备的首选方式[1]。然而,因国外相关技术的封锁,国内地下水封洞库技术在应用理论和工程实践等方面存在着许多问题[2]。地下水封石油洞库面临的问题之一是洞库容积一般较大,所储存原油一般是低凝轻质原油,在洞库进出油作业过程中会有大量的油气产生和挥发,进而导致洞库的储存压力增加,长期运行情况下积聚在洞库洞顶位置处的油气就有可能泄漏到周围岩体乃至进入大气,导致潜在的健康、安全和环境问题[3]。

目前,国内外已开展的地下洞库气体泄漏研究工作主要集中在液化石油气(Liquefied Petroleum Gas,LPG)洞库和储气库方面。Yamamoto和Pruess[4]建立了日本某地下LPG洞库的二维三相流模型,考虑了随机渗透性等因素的影响,模拟了日本LPG洞库气体泄漏现象,讨论了洞库的储存安全性以及气体泄漏的潜在原因。Javadi和Sayadi[5]采用基于运移追踪算法和路径分析的数值技术,通过离散裂缝网络(Discrete Fracture Network,DFN),计算了地下无衬砌洞库周围裂缝中气体运移的发生情况。Liu等[6]采用数值模拟方法,分析了顶板盐岩坍塌和顶板石膏坍塌后的气体泄漏情况。Chen等[7]基于低倾角特征和层状岩盐夹层结构建立了盐岩洞库气体泄漏分析模型,揭示了层状盐洞室地下储库周围深部地层的气体泄漏运移规律。我国地下水封石油洞库多建造在结晶岩体中,所储存的一般是低凝轻质原油,油品的泄漏主要有液体渗漏和气体渗漏两种。时洪斌[8]在室内平行板模型试验基础上,采用数值模拟方法进行了不同工况下储库油气泄漏情况的模拟,定性地探讨了黄岛国家石油洞库储油过程气态油品泄漏的微观机理。可见以往研究多集中于气体泄漏机理和成因方面,对于气体泄漏后整体运移规律的研究则相对较少。

在地下洞库气体泄漏控制方面,常用的方法有围岩灌浆、洞库衬砌、冷冻法以及利用天然地下水和人工水幕等[9]。其中,利用人工水幕进行水动力控制是目前最为广泛采用的措施[10]。人工水幕之所以能密闭油气,是因为人工水幕可使储存油气的洞库围岩缝隙中形成指向洞库的渗流,当这些渗流的水力坡降大于某一临界值时,就可阻止气体进入岩石缝隙或阻止已进入缝隙的气泡向外运动。在这一临界水力梯度方面,国内外学者做了一系列相关的理论和实验研究。Åberg[11]研究了注水压力与洞库储存压力的关系,发现当垂直水力梯度>1时,可以保证洞库的水封性。Goodall 等[12]在Åberg的基础上扩展了这一标准,认为只要在所有可能的渗漏路径上某段距离内水压力不断增大,则可以保证不会发生气体泄漏。如今,国内外已有不少洞库发生气体泄漏后通过增设水幕的方式来补救的工程案例,如韩国济州岛地下水封油库的试验库和浙江象山地下油库等[13]。“水封幕”这个概念最初由Lindblom等[14]通过不衬砌地下储气洞库相关研究,发现地下饱和水在裂隙岩体中可以防止气体的逸出进而提出。谭忠盛等[15]采用数值模拟方法对汕头地下水封式LPG 储库的渗流场进行分析,证明了设置水幕是控制气体泄漏的有效方法。Li等[10]设计了物理模型来模拟水幕的水封情况,讨论了气体泄漏与水幕压力的关系。Li等[16]利用离散元法分析水幕的存在是否会影响储藏库的气密性。杨荣等[17]为了解水幕系统对地下水封洞库的储藏安全的作用,采用数值模拟方法,进行藏品水-油-气三相流运移模拟。郭得福等[18]基于地质统计学原理进行研究区随机渗透率计算,建立考虑研究区渗透性强烈非均质条件的丙烷洞库二维剖面多相流,进行不同工况下油气储品泄漏的运移模拟。Ueda等[19]为了防止LPG 储库气体泄漏,根据设计标准建造了水幕系统,并向洞室周围的岩体中连续注入大量的去离子海水,以保持施工和运行期间的水位恒定。一般说来,水幕系统压力越大,水封厚度也就越大,水封效果也越好,然而水幕系统的经济合理性常常被忽略。

本文以我国某大型地下水封石油储备库建设项目为依托,基于连续介质气液两相流理论,开展大型地下水封洞库油气泄漏运移数值模拟,模拟地下水封油库储油运行期间油气泄漏运移的过程,通过对比不同水幕压力下的渗流场和油气泄漏运移情况,对比了不同水封厚度下的油气泄漏控制效果。研究成果可为水封油库工程的设计及油气泄漏控制提供一定理论参考。

1 基本原理

1.1 控制方程

由于油气和水是不混溶的,地下水封石油洞库油气的泄漏运移可以基于气液两相流理论模型进行描述。其控制方程包括地下水的流动方程和不混溶相的传递方程,二者的耦合可以在有限元软件中通过多物理场耦合方法实现,各相的质量守恒方程由多孔介质相传递接口给出,其控制方程为

(1)

式中:t为时间;εp为孔隙率;ρi是i相流体密度(kg/m3);ui为i相的体积通量(m/s);si为饱和度(无量纲);∇为哈密顿算子。体积通量由扩展的达西定律确定,即

(2)

式中:μi为i相流体动力黏度(Pa·s);kri是相对渗透率(无量纲);k为多孔介质的渗透率(m2);pi为流体压力(Pa);g为重力加速度(m/s2)。i可取值w和g,分别代表为液相、气相,下同。满足sw+sg=1。

除了求解包含各相的质量守恒方程(1)外,多物理场耦合接口还需要求解包含其中一相的流体流动方程,可以采用满足连续方程的达西定律来描述,其表达式为

(3)

式中:ρ为流体密度(kg/m3);μ是流体动力黏度(Pa·s);p为流体压力(Pa);Qm为流体质量源项(kg/(m3·s))。

考虑到材料的非线性,将密度与孔隙率定义为压力的函数,应用链式法则,得到

(4)

总体控制方程最终表述为

(5)

式中S为储水系数(Pa-1)。

毛细压力是各相之间的压力差,并且是润湿相饱和度的函数,国内外学者基于试验数据推导出了相对渗透率和毛细压力的函数关系,研究发现它们都是润湿相饱和度的函数。Brooks和Corey[20]通过入口毛细压力和润湿相饱和度得到了一个描述排水过程的经验关系,其表达式为:

pc=psg-psw,

(6)

(7)

在 Brooks Corey 模型中,两相的相对渗透率由式(8)和式(9)推导得出:

(8)

(9)

(10)

(11)

在COMSOL Multiphysics中,达西定律只适用于计算多孔介质中的单相流体流动,为模拟两相流,可采用多孔介质多相流耦合接口,它将达西定律与多孔介质相传递物理场进行耦合求解,对于两相流体的密度和动力黏度可表示为:

(12)

(13)

1.2 油气泄漏运移控制原理

在地下水封石油洞库储油运行期间,原油蒸发等作用导致大量的油气积聚在洞库拱顶部处,引起洞库内部储油压力的增加。如果水封厚度不足,随着洞库内部油气压力的增加,油气可能会驱动洞库顶部周围岩体中的地下水产生流动。这一过程可以采用基于气液两相流理论的数值方法对油气的水封效果进行模拟分析[3]。对于油气泄漏来说,油气从地下无衬砌洞库向周围饱和围岩泄漏是由洞库中的气体压力与洞库周围围岩的水压之间的压力差驱动引起的。此外,由于油气密度低于水的密度,所以油气以向上泄漏运移为主,这可能导致大量油气在地表附近聚集以及进入大气,造成十分严重的后果。洞库在无水幕系统下开挖,会造成地下水位下降,当水位下降到一定位置时,该水封厚度下的洞库周围水压会低于洞库储油时的油气压力,从而造成油气泄漏。洞库在有水幕系统下开挖,通过水幕系统提高地下水位,对应水封厚度下的洞库周围水压会高于洞库储油时的油气压力,从而抑制油气泄漏。具体地下水封石油洞库建设应根据洞库具体条件确定合适的水幕系统,以对油气泄漏进行控制。

2 水封油库油气泄漏数值模拟

2.1 计算模型

采用文献[21] 中的洞库模型进行二维气-液两相流数值模拟。模型由3个主洞库和一系列水幕孔构成,模型长328 m,高210 m,左洞库与中洞库间距为30 m,中洞库与右洞库间距为58 m。模型下部边界为坐标系原点,模型上边界模拟至平均地下水位,标高设为210 m,在标高为90 m处设有人工水幕。为研究不同水封厚度对油气泄漏运移过程的影响,本文将水幕孔位置设在洞顶上方10 m处,水幕钻孔平行于洞库轴线,钻孔直径0.1 m,相邻钻孔间距10 m,左右各延伸至洞壁外侧16 m。模型左右边界各从洞壁外侧延伸90 m,底部边界从洞底向下延伸60 m(标高0 m)。网格剖分为四边形和三角形单元,洞库和水幕孔附近适当加密,模型计算完整网格单元数为20 068。数值计算采用的网格、初始及边界条件如图1所示。

图1 计算模型网格及边界条件Fig.1 Meshes and boundary conditions of computation model

2.2 初始条件及边界条件

初始条件:在洞库还未开挖之前,渗流场中初始地下水压力按静水压力计算。将洞库开挖3 a后得到的渗流场和饱和度计算结果作为后续运行期数值模拟的初始条件。油气泄漏运移模型分达西渗流场和多孔介质相传递场2个物理场,因此初始及边界条件也有2套。

对于地下水渗流的边界条件,模型左右边界为压pw=ρwg(210-y)。模型顶部和底部边界为零流量边界。施工工况下,洞库设计地下水位距洞库顶部120 m,水幕孔压力按表1设置,洞库边界按0 m压力水头处理;在储油工况下,液态油品浮在水垫层上,液态油品上部为饱和蒸气,本文算例中忽略水垫层的厚度,油品密度取850 kg/m3,饱和蒸气压取0.2 MPa,洞库拱顶、边墙和底板边界根据油气压力、储油压力分别施加。

表1 计算参数取值Table 1 Values of computation parameters

对于两相相传递场的边界条件,模型左右边界、水幕孔为液相边界,洞库顶拱设为气相边界,考虑本文主要研究洞顶处油气的泄漏,因此洞库边墙、底板按液相边界设置。模型底部边界为无通量边界,考虑油气有可能进入大气,模型顶部边界设为流出边界。

2.3 参数取值

库址区岩性属于花岗质片麻岩,参考文献[21] 、文献[22] ,洞库岩体渗透系数取kr=1.0×10-9m/s,孔隙率εp取为0.05;流体参数来自于文献[23] 。计算参数取值见表1。

通过调整水幕孔压力大小,可以获得不同的水封厚度,从而研究不同水封厚度对油气泄漏运移规律的影响。如表2所示,建立了无水幕系统模型(工况1)和有水幕系统不同水幕压力模型(工况2—工况6),对洞库经3 a施工期后再储油运行50 a的情况进行分析,其中洞库假设一次性开挖完成。

表2 数值模型工况Table 2 Cases of numerical simulation

3 结果与分析

3.1 施工期渗流场分析

地下洞室的开挖会改变洞库区的渗流场,在研究洞库储油期油气泄漏运移过程中,需要关注洞库开挖后周围的地下水位和孔隙水压力的分布情况。图2为施工期无水幕情况下(工况1)洞库一次性开挖3 a后的地下水位和孔隙水压力分布云图。由图2可知,在洞库无水幕开挖3 a后,由于洞库的施工开挖排水作用,原来的地下水位下降幅度明显,在洞库区形成了明显的水位降落漏斗,洞室顶部出现了大面积的疏干区,洞室内空气与大气直接接触,造成后期无法储油。可见,为确保洞库水封性,在洞库上方设置水幕系统非常必要。

图2 施工期末孔隙水压力分布云图(无水幕)Fig.2 Contours of pore water pressure at the end of construction in the absence of water curtain

图3为有水幕情况下施工期结束时不同工况对应的地下水位线。总体看来,由于洞室的开挖,洞库上方岩体中的地下水位都会形成较大幅度的下降,但由于水幕孔注水作用,水位降幅较无水幕开挖工况下要小,洞库上方岩体中地下水位没有出现降落到洞顶的情况,即洞库上方没有形成与大气相通的疏干区,洞库上方保留了一定的水封厚度。此外,工况2—工况6下地下水位线最低点到洞库的垂直距离(即最小水封厚度)分别约为14.51、22.29、29.12、35.23、40.74 m,这表明水幕孔注水压力越大,地下水位线的下降幅度越小,最小水封厚度越大。

图3 施工期末不同工况对应地下水位线(有水幕)Fig.3 Underground water levels in different cases at the end of construction in the presence of water curtain

3.2 储油期油气泄漏运移规律分析

通过对有水幕系统情况下不同工况油气泄漏运移分析发现,5个工况的油气泄漏运移规律类似,下面以工况2为例进行储油期油气泄漏运移规律分析。

图4为洞库不同储油期的油气饱和度分布云图,油气饱和度的分布能够表明油气向洞库外逃逸泄漏的行为和规律。由图4可知,储油运行1 a后,便在洞顶上方一定范围处积累大面积高饱和度的油气泄漏圈,并且随着储油时间的增加,油气不断向周边低压区扩散,泄漏圈的面积也随之增加,左洞库和中洞库由于距离近在储油运行10 a后出现了泄漏圈连通的情况。

图4 工况2下不同时间气相饱和度分布云图Fig.4 Contours of gas-phase saturation at different moments in case 2

以中间洞库竖直方向油气泄漏距离进行为例,储油运行1 a后,竖直方向上最远油气泄漏距离约为15 m;储油运行10 a后,竖直方向上最远油气泄漏距离约为27 m;储油运行30 a后,竖直方向上最远油气泄漏距离约为51 m;储油运行50 a后,竖直方向上最远油气泄漏距离约为72 m,这表明油气泄漏范围随着储油时间的增长而增加。

图5为工况2油气泄漏范围和泄漏量与储油运行时间的关系。本文定义油气泄漏范围A为油气饱和度sg≥0.01时所占据的面积大小,计算公式见式(14)、式(15),其中0.01是为考虑数值计算的稳定性而取的比0稍大的数值;油气泄漏量Q为油气所占据孔隙空间的体积(因本文为二维计算,单位为m2),计算公式见式(16)。

图5 工况2下油气泄漏范围和泄漏量的拟合曲线Fig.5 Fitting curves of range versus volume of oil vapor leakage in case 2

(14)

A=∬Df(sg)dS,

(15)

Q=∬DsgεpdS。

(16)

式中:dS表示单元的面积;D表示整域。

由图5可知,油气泄漏范围和油气泄漏量与储油运行时间整体呈正幂函数关系,拟合公式分别为A=1 418.18t0.28和Q=23.47t0.13,由决定系数R2可知,以上拟合公式能够准确描述储油期油气泄漏范围和泄漏量的变化趋势,具体表现为储油0～10 a时,油气泄漏范围和泄漏量迅速增加,10 a后的趋势则逐渐平缓。需要说明的是,这里的拟合公式是采用本项目的数据拟合得到的,仅适用于本工程,其它工程可以采用类似方程形式得到相应的拟合公式。

图6为不同水幕孔压力工况下油气泄漏范围演化过程。可以看出,各工况下油气泄漏范围分布规律大体相同,即中间洞库上方的油气泄漏面积均要大于左右两洞库,这是由于中间洞库水封厚度均要低于左右两洞库。但是,油气泄漏范围则随着水幕孔注水压力的增加而逐渐变小。相比于工况2和工况3,水幕压力从150 kPa开始,油气泄漏范围明显缩减,左洞库和中间洞库的油气泄漏范围不再出现相连通的情况。在工况6下,储油运行10 a后油气泄漏范围基本没有变化,保持在一个非常小的范围。

图6 不同工况下油气泄漏范围Fig.6 Ranges of oil vapor leakage in different cases

3.3 油气泄漏控制效果分析

图7为不同工况下油气泄漏范围和泄漏量随着储油运行时间变化曲线。由图7可知,在储油运行50 a期间,各储油时间段油气泄漏范围和泄漏量大小均为:工况2>工况3>工况4>工况5>工况6,这表明水幕孔压力越大,对油气泄漏控制效果越好。从图7还可知,工况2至工况6水幕孔压力逐级增量为50 kPa,油气的泄漏范围和泄漏量减小最明显阶段是工况3至工况4阶段,即水幕孔注水压力从100 kPa增加至150 kPa期间。此外,随着水幕压力从工况4的150 kPa再继续增加,对油气泄漏范围和泄漏量控制效果增加不明显,这表明水幕压力对油气泄漏的控制存在经济合理的范围。

图7 不同水幕压力下油气泄漏参数变化曲线Fig.7 Curves of oil vapor leakage parameters under different water curtain pressures

图8为工况2—工况6在储油50 a后油气泄漏范围(面积)、竖向最远泄漏距离和油气泄漏量与最小水封厚度之间的关系。可见,随着最小水封厚度的增加,油气泄漏范围和泄漏量逐渐减小,但是减小的幅度并不是线性变化的,最小水封厚度在30 m以内时,油气的泄漏范围和油气泄漏量减小非常显著,如水封厚度从约为15 m增加到30 m时,油气泄漏范围减少了约58%,油气泄漏量减少了约73%。当最小水封厚度>30 m时,油气泄漏范围和泄漏量减小趋势减缓。这表明虽然水封厚度越大,油气泄漏范围和泄漏量越小,但过大的水封厚度会大大增加工程成本,而油气泄漏控制效果并不能对应线性增加。对于本工程而言,水封厚度为30 m时对油气泄漏控制最为显著,所以建议水封厚度设置在30 m即可。

图8 不同水封厚度油气泄漏参数变化曲线Fig.8 Curves of oil vapor leakage parameters under different water sealing thickness

3.4 计算结果讨论

地下油库的水封实质是利用洞库围岩中水压大于油气压力从而将油气密封在一定封闭空间内,油气是否泄漏的关键在于是否满足一定的水体厚度。在工程当中,通常利用水幕系统人为增大洞库周围水压力来封存洞库中聚集的油气。根据《地下水封石洞油库设计标准》[24],洞室拱顶上距设计稳定地下水位垂直距离不应小于Hw,按式(19)计算。

Hw=100P+20 。

(19)

式中:Hw为地下水位至主洞库拱顶的最小垂直距离(m),即最小水封厚度;P为洞库内的气相设计压力(MPa)。本模型中主洞库运行期气体压力P=0.2 MPa,则由式(19)计算得到的最小水封厚度Hw=40 m。本文基于气液两相流理论对地下水封石油洞库在运行期油气泄漏运移规律进行了研究,发现当水封厚度设置为30 m时最为经济且能保证水封效果。而水封厚度为40 m对应本文中工况6的情况,从图7可以看出,工况6的水封效果也是非常显著的。需要说明的是,本文的模型采用的是均质模型,考虑到实际工程中岩体中裂隙的存在,《地下水封石洞油库设计标准》建议的40 m水封厚度是合理且具有一定的安全裕度的,因此本文计算结果也是对于《地下水封石洞油库设计标准》推荐公式的佐证。