煤矿采空区碳封存CO2泄漏地表扩散规律研究

丁 洋,陈文彬,林海飞,朱 冰,谭军红

(1.西安科技大学 安全科学与工程学院,陕西 西安 710054;2.西安科技大学 西部矿井开采及灾害防治教育部重点实验室,陕西 西安 710054)

0 引 言

近年来,全球气候变化加剧与生态环境不断恶化引起各国对碳减排的高度重视[1-3]。中国作为碳排放大国,2020年9月22日,习近平总书记在第七十五届联合国大会一般性辩论上郑重宣布中国CO2排放力争2030年前达到峰值,努力争取2060年前实现碳中和[4]。实现碳中和的根本路径是减排增汇,碳捕获、利用和封存(Carbon Capture,Utilization and Storage,CCUS)被认为是脱碳的关键战略,而CO2封存是CCUS链的最后一环,地质封存是目前发展最为成熟的一种封存方式[5-9]。

煤炭长期以来是中国的主体能源,在煤炭资源开采过程中,形成了大量的地下采空区[10-11],由于各种因素的限制,在煤矿采空区依然残留大量煤炭,煤基质中的CO2吸附优势要强于CH4,CO2注入采空区会与残留煤炭中的CH4产生竞争吸附[12-15],从而实现CO2以吸附态封存于采空区。但采空区原始应力遭到破坏,裂隙演化和渗流规律变得复杂[16-18],这导致了采空区封存CO2相较于其他封存方式可能具有更高的泄漏风险,而一旦泄漏,则可能对植物、动物以及人类造成严重威胁。因此,明确采空区封存CO2泄漏至地表后的扩散规律和浓度分布具有重要的安全和环保意义。

对于CO2泄漏至大气后的扩散规律,前人利用数值模拟进行了大量的研究。MAZZOLDI等将Kit Fox试验和高压管道泄漏CO2的扩散模拟进行比较,研究表明CFD模型的准确性比高斯模型更高[19]。XING等利用CFD软件模拟了CO2井喷事故的扩散情况,并将试验数据进行了对比验证,得到标准k-ε和SSTk-ω模型的计算结果与CHANG标准下的试验结果吻合较好[20]。WEN等在开源CFD代码OpenFOAM上开发了一个专门用于CO2在大气中扩散的求解器,模拟了在无地形影响下CO2的泄漏情况[21-22]。WANG等研究了地形和释放源温度对复杂地形条件下埋地管道CO2泄漏扩散的影响[23]。众多学者在CO2地表泄漏扩散方面取得了丰硕的研究成果,对于CO2的扩散模拟,目前大部分研究是高压储罐和运输管道的泄漏,对于常压泄漏扩散机制的研究较少,并且在研究地形对CO2扩散的影响方面,多是用单一的地形起伏来进行表征,对于真实地形条件下的复杂影响研究较少。

文中根据试验矿井采空区地表真实地形提取地形特征,建立CFD模型,利用Fluent模拟研究不同影响因素下常温常压CO2在封存区上方地表浓度分布规律和危险的扩散范围,为采空区封存CO2项目安全评价提供参考。

1 地形特征提取和数值建模

1.1 试验矿井位置

陕西省神府矿区历经多年的开采,形成了大量采空区,以浅埋煤层为主,并且有低透气性的红土层保证了低压储存CO2的可行性。同时,榆林能源化工基地每年会产生大量的纯CO2,基地与矿区处于同一区域,便于项目的开展,降低成本。因此,项目选择陕西省侏罗纪煤田神府矿区某煤矿作为典型浅埋煤层采空区吸储CO2实施地点。项目的中试位置拟选定在试验矿井首采25201工作面,试验平面范围初步确定为300 m×500 m,其封存位置如图1所示。

图1 CO2封存区域Fig.1 CO2 storage area

1.2 地形特征提取和数值建模

为了观测CO2泄漏后封存区域及周边的CO2分布,选择封存地点周边1 500 m×1 500 m的区域作为模拟区域。地形提取流程如图2所示,将所选区域的瓦片数据下载拼接后导入ArcGIS软件中,结合封存地点的Dem数据提取地形特征,再根据中央经线进行投影,将投影后的数据导入Global Mapper软件得到了研究区域的高程图(图3),该区域最高海拔1 225 m,最低海拔1 091 m,海拔高差134 m。最后在Sketch UP软件中完成建模。

图2 建模流程Fig.2 Modeling process

图3 研究区域高程图Fig.3 Elevation in study area

2 数值模拟

2.1 模拟方案

文中利用Fluent数值模拟软件,分析了不同泄漏速度、不同风速下CO2在目标区域的扩散规律,并且在平坦无起伏的地形模型上进行了对照模拟,以便观察分析地形影响。泄漏点设在图3红色圆圈区域内,该区域位于计算域的中部位置,更有利于观测CO2泄漏扩散规律,并且注入井也位于封存区域的中部位置,泄漏点设在此处能模拟出井筒泄漏时的CO2扩散情况,同时该区域四周海拔都相对较高,便于研究地形的影响。设定泄漏口为直径1 m圆孔,泄漏时间3 600 s,泄漏速度用CO2泄漏的质量流量来表征。模拟方案见表1。

表1 模拟方案Table 1 Simulation scheme

采空区封存CO2泄漏特点与管道或储罐泄漏特点不同。管道或储罐泄漏是由于高压突然释放,因此泄漏速度快、但并不持续;而采空区泄漏属于常压泄漏,泄漏速度缓慢、但持续泄漏。因此为了充分模拟这种泄漏特点,将泄漏速度根据梯度由小到大设置,并且对于每种情况,都模拟了从泄漏开始3 600 s后的CO2泄漏扩散结果。为了充分了解CO2泄漏的危害范围,选择较高的泄漏速度即10 kg/s作为初始条件模拟其他因素变化对CO2扩散的影响。根据调查,封存区域所在地风速常年小于5级(5级风速为8～10.7 m/s),因此选择0～10 m/s的风速进行数值模拟。由图3可知,研究区域总体海拔是沿着y轴正方向逐渐降低,因此为了研究地形和风速的相互影响,研究方案中将风速变化时的风向设为沿y轴负方向。

2.2 基本假设与流动方程

2.2.1 基本假设

对于采空区封存CO2泄漏释放到地表的气体,在抵达地表后,其压力温度已经和环境值基本一致,所以假设其泄漏压力温度与大气压和环境温度相同,即常温常压泄漏,因此用于扩散模拟的源流体为气态CO2。模拟研究基于以下假设：①气体为连续介质;②气体流动过程为不可压缩;③气体流动不考虑传热与化学反应;④地表组分传输模型只考虑空气和CO2。

2.2.2 流动方程[24]

1)质量守恒方程。

(1)

2)动量守恒方程。惯性(非加速)参考系中的动量守恒由下式描述

(2)

(3)

式中μ为分子粘度,Pa·s;I为单位张量,右侧的第二项是体积膨胀的影响。

3)能量守恒方程。

(4)

上述3个方程流体运动的基本控制方程,主要是将流体运动时必须遵循的物理定律用方程的形式表达出来,从而解出流体流动的未知量。

4)组分输运模型。组分输运模型是模拟混合物各组分之间或与其他相之间的相互作用,对于CO2泄漏扩散,主要是模拟CO2和空气的相互作用。选择求解组分的守恒方程时,Icepak通过解第i个组分的对流扩散方程来预测每个组分的局部质量分数Yi。这个守恒方程的一般形式如下。

(5)

式中Yi为物质i的质量分数,即质量扩散矢量,而Si为物质i的净生产速率。

湍流中的质量扩散

(6)

该方程主要表征CO2在地表扩散的行为。其中Di,m是组分i的质量扩散系数,Sct是湍流施密特数,默认设置为0.7。

5)标准k-ε模型的传输方程。湍流是自然界普遍的流动形式,湍流运动的特征是在流动过程中流体各质点具有的互相混合的特性,速度和压力等物理量在空间、时间上具有随机性的脉动值,而湍流模型就是表征湍流运动特点的方程。湍流动能k及其耗散率ε由以下输运方程获得

(7)

(8)

在这些方程中,Gk为由于平均速度梯度产生的湍流动能;Gb为由于浮力产生的湍流动能,根据k-ε模型中浮力对湍流的影响计算。YM为可压缩湍流中的脉动膨胀对总耗散率的贡献。C1ε,C2ε,C3ε为常数;σk,σε分别为k和ε的湍流普朗特数;Sk,Sε为用户定义的源项。

2.3 计算域和边界条件

根据图3高程数据建立采空区地表模型,通过红圈所在位置确定泄漏口中心点坐标为(715,725,47.22)。为了适应地表134 m的高程差异,以模型在z轴上的点即(0,0,114.48)作为基准点,建立高度为300 m的立方体模型。模型计算域如图4(a)所示。

图4 计算域和局部网格特征Fig.4 Computational domain and local mesh features

图4(b)显示了扩散模型的面网格特征,由于几何形状比较复杂,计算区域主要离散为六面体单元,其中分别对地表和泄漏源附近的网格进行了不同程度的加密。在最终仿真之前,进行了网格独立性研究,以确保网格对结果的影响可以忽略不计,最终确定的网格大约由400万个单元组成。

扩散模型边界条件设置如下：①CO2源：质量流量入口,环境压力和温度,质量流率根据模拟方案所述设定;②进风口：速度入口,环境压力和温度,风速根据模拟方案所述设定;③出口：具有环境压力和温度的压力出口,模拟无风条件时,进风口也设为压力出口;④地面和顶部：温度等于环境温度的防滑等温墙。

2.4 模型验证

利用上述k-ε模型和数值方法模拟了平地上泄漏速度10 kg/s,风速6 m/s的CO2扩散。WEN等利用SSTk-ω模型模拟了CO2和空气的混合物从泄漏源以50 kg/s的泄漏速度、2 m/s的风速泄漏至大气的扩散情况[21]。图5比较了距离泄漏源下风向200 m处CO2浓度的试验和两者的模拟结果,k-ε和SSTk-ω模型的模拟结果比较表明,两者的变化趋势具有高度的一致性,k-ε模拟的较大风速造成了CO2到达200 m处的时间更早,SSTk-ω模拟的较大泄漏速度造成了同一时刻更高的CO2浓度,符合实际情况,因此本研究所用的k-ε模型具有较高的可靠性。

图5 距离泄漏源下风向200 m处CO2浓度的试验和模拟结果的比较Fig.5 Comparison of experimental and simulated results of CO2 concentration at 200 m downwind from the leakage source

3 结果与讨论

研究的主要关注点是CO2浓度在10 000 ppm(即1%)水平的范围。根据参考文献[20]的研究结果,当CO2浓度超过1%时,可能引发人员出现气喘、头痛、眩晕等不适症状。为了分析CO2在浓度范围为0%～1%时的扩散规律,绘制了CO2浓度的云图,并将1%浓度水平的等值面作为对人体可能产生危害的范围进行描述。这一范围的描述包括在纵向(y轴方向)和横向(x轴方向)的影响距离以及等值面的面积。需要明确的是,上述危害范围仅表示对人体有害的范围,并非CO2扩散的范围(同样适用于后续描述)。在竖直方向(即z轴方向)上,由于重力的作用,CO2的扩散距离较小,主要分布在地表空间内,并且相对稳定。因此,本研究不考虑竖直方向上的扩散,而将地表危害区域作为立体空间内的危害区域加以考虑。

3.1 平坦无特征地形的泄漏

为了评估真实地形对危害范围的影响,提前模拟了平坦无特征地形上的CO2扩散作为对照组。图6显示了无风条件下泄漏速度为10 kg/s时,泄漏不同时间的CO2浓度分布云图。无风时,CO2泄漏后往四周均匀扩散,扩散范围逐渐增大,但高浓度区域相对于低浓度区扩散较为缓慢。

有风的条件下,整个模拟分为2个步骤进行：①稳态模拟：根据模拟方案提供的风速在地形上建立风场,为下一步的模拟提供初始条件;②瞬态模拟：在风场下模拟CO2在计算域内的泄漏扩散。图7显示了在平坦地形上,CO2在泄漏速度10 kg/s、风速6 m/s、风向沿y轴负方向时,泄漏不同时间的浓度分布云图。此时CO2主要沿着风向扩散,在190 s左右,其计算域内的浓度分布云图就趋于稳定,高浓度区域扩散相对于低浓度区域扩散也较为缓慢。

图7 无风时不同泄漏速度下CO2泄漏不同时间的浓度分布(0%～1%)云图Fig.7 Concentration distribution(0%～1%)of CO2 leakage at different time under different leakage velocity without wind

表2显示了以上2种情况泄漏3 600 s后的1%浓度等值面数据,蓝色区域代表可能对人体产生危害的区域。由表可见,无风和有风条件下的影响面积存在显著差异,这是因为无风时CO2均匀扩散,有风时CO2沿风向扩散,导致纵向影响距离增加,横向影响距离缩小,并且风力加速了CO2的扩散,使CO2累积变得困难,所以相比无风条件高浓度区域面积急剧缩小。

3.2 真实地形条件下泄漏速度的影响

图7显示了无风时不同泄漏速度下CO2浓度为0%～1%的扩散云图。在无风条件下,CO2泄漏后呈不规则分布,泄漏速度越大,CO2的横向和纵向扩散范围越大。泄漏速度较小时,CO2主要往地势较低的区域扩散,当CO2泄漏速度足够大时,其扩散也能克服部分地形影响,往地势高的区域扩散,但相对于地势低的方向,扩散有限。总体而言,CO2扩散的趋势是一致的,泄漏速度的变化只会影响扩散范围的大小。无风条件下,CO2扩散主要受2个因素的影响,一是地形因素,由于CO2是重气,在重力的作用下会下沉,因此会往地势低的区域扩散;二是浓度因素,CO2会由浓度高的区域往浓度低的区域扩散。

在相同的泄漏速度和扩散时间内,平坦地形的CO2扩散范围大于真实地形,因为CO2在平坦地形是从泄漏点往四周均匀扩散,而在真实地形CO2主要沿着y轴正方向运动,地势的升高限制了其他方向的扩散。然而,在y轴负方向,真实地形的扩散距离大于平坦地形。

表3显示了无风条件下CO2泄漏3 600 s后浓度为1%的等值面数据。由表可见,随着泄漏速度的增加,CO2扩散的横纵向影响距离以及影响面积都逐渐增加,因此泄漏速度是影响CO2扩散的关键因素,直接决定了CO2泄漏之后的危害范围。当泄漏速度为10 kg/s时,与平地相比,真实地形下的影响面积更大,说明在无风条件下封存区域地表整体地形更有利于CO2的局部累积,从而造成CO2浓度1%以上的区域面积较大,这一点在预防真实地形发生CO2泄漏方面需要引起重视。

表3 无风情况下真实地形泄漏3 600 s后的1%CO2浓度等值面数据Table 3 Isosurface data of 1% CO2 concentration after 3 600 s of real terrain leakage under no wind condition

3.3 真实地形条件下风速的影响

图8显示了在泄漏速度为10 kg/s,风向沿y轴负方向时,不同风速下的CO2浓度分布云图。在CO2云图前段,结合图3高程图所示,CO2有明显的沿山谷扩散的情况,这是由于扩散云图前段风向与山谷走向几乎一致,但在后段山谷走向发生了偏移,此时,风为CO2提供的力克服了地形影响,CO2不再沿山谷扩散。在可见区域(计算域)内,CO2在风力作用浓度分布趋于稳定,即随时间变化幅度很小或不变,并且风速越大,其达到稳定越快。

图8 泄漏速度10 kg/s时不同风速下CO2泄漏不同时间的浓度分布(0%～1%)云图Fig.8 Concentration distribution(0%～1%)of CO2 leakage at different time under different wind speeds at the leakage rate of 10kg/s

风速2 m/s时,在近泄漏源处,高浓度区域呈现沿x轴负方向偏移的倾向,这是由于地形的影响,从图3中可以看出,在该区域右侧是陡峭的山坡,上下落差超过50 m,所以阻挡了大部分CO2沿x轴正方向的横向扩散。当风速逐渐增大后,以上情形不再出现。风速相同时,CO2在真实地形的扩散达到稳定状态所需的时间比平坦地形更长,因为当风沿y轴负方向吹时真实地形下的CO2扩散需要克服地形对它的阻碍作用,而平坦地形则不需要。

表4显示了不同风速下CO2泄漏3 600 s后浓度为1%的等值面数据。由表可知,风速增大,CO2扩散的危害范围并没有随之扩大,其横向影响距离随风速的增大而减小,但其纵向影响距离是先增大后减小,这可能是由于风速到达一定值时,会加速CO2与空气的混合,从而造成了泄漏CO2的快速稀释,使高浓度区域的覆盖范围减小。

表4 有风情况下真实地形泄漏3 600 s后的1%CO2浓度等值面数据Table 4 Isosurface data of 1% CO2 concentration after 3 600 s of real terrain leakage under windy conditions

4 结 论

1)试验矿井采空区封存CO2的地表泄漏速度越大,扩散范围越大,危害范围也越大;风速越大,扩散速度越快,横向扩散越小。

2)在无风和风速较小(小于2 m/s)时,CO2会主要沿着地势低的方向扩散。

3)当风速大于等于2 m/s时,风和地形对CO2扩散都具有重要作用,并且随着风速逐渐增加,CO2的扩散会慢慢受风的主导,因此一旦发生泄漏需要注意泄漏源下风向的CO2浓度变化。