孔隙率对埋地超临界CO2管道泄漏扩散的影响

陈兵， 赵琼， 郭焕焕

(西安石油大学机械工程学院， 西安 710065)

全球变暖对人们生活和生产的影响与日俱增，将生活及工业产生的废气捕集净化后，加工提取其中的CO2并输送至目的地利用和封存，这一高效减排CO2技术称为CCUS(carbon capture、use and storage)。在CCUS过程中，碳运输是一个重要的中间环节，管道运输因其输送量大、可长期持续输送等众多优点，成为最经济可行的方式[1-2]。为获得更高的传输效率，管道输送的CO2一般处于超临界相态，超临界态的CO2密度与液相密度接近，黏度系数非常小，接近于气相且具有非常高的扩散性[3-4]。管道的施工问题或水合物形成而产生的局部腐蚀等原因易导致管道发生泄漏，泄漏的CO2会下沉到地面并聚集在低凹处，过高的CO2浓度会导致人和动物窒息[5]。

国外在CO2管道方面的模拟研究与实际应用已经十分成熟。Cortis等[6]采用在稳定状态下的Navier-Stokes方程，建立二维系统研究CO2封存点泄露扩散的密度效应对其在大气中扩散的影响。Wen等[7]使用FOAM(field operation and manipulation)求解器，采用均匀的弛豫模式，基于湍流动能k和耗散率ω的k-ω湍流模型对埋地低雷诺数流体管道穿孔时释放致密相CO2并在大气中的扩散进行模拟研究，并将模拟结果应用于实际工程中。针对埋地CO2管道泄漏这一课题的探讨，中国仍然处于发展阶段，实验及模拟的相关成果也逐渐呈现。李康[8]通过自主搭建的地面超临界CO2管道泄漏扩散实验平台，采用理论模拟与实验比对相结合的方法，研究不同影响因素对泄漏口附近参数的影响。刘正刚[9]通过自主搭建的三套CO2埋地管道输送装置，对发生穿孔和断裂的CO2管道在气相和液相这两种输送相态下，通过实验验证泄漏口径和泄漏口压力等对泄漏和泄放的影响。陈国龙等[10]认为现有的泄漏事故均为小孔泄漏(泄漏孔径小于20 mm)，搭建了泄漏孔径为1 mm和3 mm的实验台。郭晓璐等[11]对国内外泄漏扩散进行了整理，认为初期泄漏为小孔泄漏，持续发生泄漏则会逐渐变为大孔泄漏(泄漏孔径在20 mm至管径之间)。刘少荣[12]根据工程多发事故对管道进行了孔泄漏实验，得到了不同泄漏孔径下的泄漏规律。孔泄漏和横向断裂是长输管网的两种泄漏方式，但在实际中，管道一旦发生泄漏导致工况不稳就会被监测并抢修，所以工程上CO2管道泄漏事故多为小孔泄漏。

中国对CO2管道泄漏扩散的研究大多停留在气相和液相短程架空管道的实验及模拟。而输送量大、输送效率高且温度压力更适合CCUS工程要求的超临界CO2长输管道相关研究和应用还没有深入展开。由于超临界CO2管道基本为埋地管道，只有少数管段会因地理因素架空，所以针对埋地CO2管道泄漏过程的扩散规律展开研究在工程上具有更高的实际价值。根据国外学者的研究方法，基于标准k-ε高雷诺数模型，建立埋地管道三维模型，与中国相关的实验数据对比验证模型的正确性，对超临界CO2埋地输送管道因腐蚀等因素导致管道发生初期泄漏(小孔泄漏)后，模拟土壤孔隙率对埋地管道不同截面处CO2浓度扩散规律的影响。

1 理论基础与控制方程

以流体力学[13]为基本理论，根据多孔介质理论建立埋地管道发生小孔泄漏空间模型[14],在此基础上利用GAMBIT软件采用疏密结合的方式建立网格；结合流体力学基本方程以及埋地管道的Fick扩散模型，使用专业模拟软件FLUENT选用湍流模型中的Realizablek-ε模型[15-18]以及合适的控制方程并设置合理的边界条件进行数值计算，经迭代收敛得到模拟结果。

1.1 多孔介质

1.1.1 孔隙率

由于土壤这一多孔介质的内孔隙参差错落，故采用孔隙率来度量其空隙性质，即

(1)

式(1)中：Vp、Vt为各孔隙的体积和整体单元体积，并默认在垂直方向孔隙率是均匀的，以便于数值模拟计算。

1.1.2 达西速度

可用达西定律来描述一维流体运动，即

(2)

式(2)中：vx为x轴方向的达西速度；K为流体渗透率；μ为流体的黏度；p为流体压力。

1.1.3 多孔介质的状态方程

只考虑弹性变形范围之内的情况，当压力差不大时，孔隙率变化的状态方程为

φ=φ0(1+CφΔp)

(3)

式(3)中：φ0为参考孔隙率；Cφ为孔隙压缩系数；Δp为相邻孔隙之间的压力差。

1.2 流体力学的基本理论

埋地CO2管道泄漏扩散的影响因素有流速、压力、密度以及温度，所以在直角坐标系下适用于CO2流体的三大守恒定律。

1.2.1 连续方程

可压缩流体在多孔介质中的连续性方程为

(4)

式(4)中：t为时间；ρ为密度；v为速度矢量。

1.2.2 运动方程

将所研究的CO2流体视作可压缩流体，在忽略质量力的情况下，流动气体的N-S(Navier-Stokes)方程为

(5)

式(5)中：ν为速度矢量；μ′与μ分别为体积与动力黏性系数；D为变形张量。

1.2.3 能量方程

要使流体在多孔介质中扩散的模拟过程遵循物理规律，应使用能量方程加以限制，即

(6)

式(6)中：e为内能；f为质量力矢量；P为压力张量；λ为热传导系数；T为绝对温度；q为辐射热量。

1.2.4 状态方程

将CO2流体的压缩性与多孔介质中的空腔结构相结合，得到合适的状态方程。

压缩系数为

(7)

当压力差不大时，式(9)可表示为

ρ=ρ0(1+cfΔp)

(8)

式中：V为流体体积；cf为流体压缩系数；ρ0为参考密度。

1.3 埋地管道的Fick扩散模型

根据Fick扩散模型，气体在土壤中的非稳态泄漏扩散浓度为

(9)

式(9)中：C为扩散物质的浓度，mol/L；M为扩散物质的质量，kg；D为扩散系数，m2/s；x为扩散层的深度，m；t为物质扩散的时间，s。

1.4 k-ε湍流模型

由于CO2黏滞系数较低可以忽略，可以假设流体在多孔介质中为完全湍流。故选用广泛应用于完全湍流过程的k-ε标准模型，其湍动能k和耗散率ε方程形式为

(10)

(11)

式中：速度梯度与浮力对k的影响系数分别为Gk和Gb；YM为由脉动膨胀引起的ε的变化；μt为湍流黏性系数；1/σk、1/σε为k、ε的湍流普朗特数，C1ε、C2ε、C3ε为常参数。

1.5 建立理想化模型并设置边界条件

FLUENT软件在流体力学基础理论的支持下集成了多种解决湍流、传热与相变、多相流等理论模型，可以利用有限体积法等算法计算模拟复杂流体的动态过程，通过其组分输运模型可以得到流体各组分的扩散特性，所以软件非常适用于土壤中泄漏以及扩散的研究。合理建立理想化模型有利于软件的定量计算。

土壤简化为各向同性的多孔介质，孔内为标压空气； 由于小孔泄漏过程所以忽略泄漏过程对土壤结构的改变；由于多孔介质导热系数很低，小孔泄漏规模较小，故忽略气体与土壤之间的热交换；管道输送CO2浓度为99%。

在理想模型的基础上采用疏密结合的方式进行网格划分，其入口边界条件、出口边界条件以及管壁边界条件的具体设置如表1所示。

表1 模型边界条件

将网格参数导入FLUENT，考虑重力影响因而启动非定常求解器，由于处理过程为完全湍流过程故采用k-ε湍流模型，开启组元输运即可追踪CO2的泄漏扩散过程。以土壤颗粒直径Dp=1.6×10-3mm，孔隙率φ=0.45为例，根据Ergun公式计算可得黏性阻力和惯性阻力系数分别为C1=101 035 302.91 m-2，C2=13 203.02 m-1。

2 模拟方案可行性验证

刘正刚[9]将发生小孔泄漏的气体CO2管道埋入沙箱中，经过实验得到了管道发生泄漏时，不同因素的影响下CO2的扩散规律。结合上述理论模型开展与沙箱实验相对应的仿真模拟过程，模拟结果与实验结果相比对，从而验证所建模型的可靠性。

2.1 数值模拟结果

设置好管道泄漏扩散模型的入口边界条件、出口边界条件以及管壁边界条件，选定设置的基本类型，在表2所示工况下，监测泄放口上方地表的固定监测点C1(距泄放口水平距离和竖直距离分别为0.2 m和0.15 m)的CO2浓度(体积分数)随时间的变化情况，并对比在不同泄漏口径下模拟数据以及实验测量数据，对比结果如图1所示。

表2 实验工况参数

图1 不同泄漏口径下C1点二氧化碳浓度随时间的变化Fig.1 Changes of CO2 concentration at C1 point under different leakage caliber with time

2.2 误差计算

将图1所示的实验数据与模拟数据进行比对，分析其中误差，结果如表3所示。

表3 实验、模拟数据对比

2.3 结果分析

由图1和表3可得实验和模拟在不同因素作用下固定点测得的CO2浓度随时间变化曲线是相接近的。实验时泄漏产生的膨胀波会使CO2浓度发生波动，膨胀波会随着泄漏孔径的增大而更加明显、更加剧烈。但本文中模拟的是小孔泄漏的稳态扩散过程，CO2浓度曲线较为平稳，故认为模拟数据和实验结果相比偏差较小，因而反映了本文所建立模型的正确性。

3 孔隙率对泄漏CO2扩散的影响

3.1 建立埋地管道模型

3.1.1 输送CO2管道参数

以中国某油田超临界CO2埋地管道为研究对象，将其输送工况参数转化为模拟实验模型的参数，数据如表4所示。

3.1.2 管道-小孔泄漏模型

研究小孔泄漏初期扩散分析[19]的具体模型如图2所示。如图3所示，结合图2以及表4中的管输参数建立了4 m×3 m×3 m的埋地输气管道小孔泄漏模型，将直角坐标系的XOZ面设置为地表，X轴的正反向为CO2输送方向，则泄漏口的坐标为(0，-1.8，0)m。

表4 CO2输送管道参数汇总表

图2 埋地管道小孔泄漏扩散模型Fig.2 Model of small hole leakage diffusion in buried pipeline

3.1.3 网格划分

采用GAMBIT软件网格划分埋地管道三维模型，计算基于稳态条件地表CO2浓度在不同单元数网格下的扩散情况，验证网格的无关性，结果如图4所示。

图4 网格无关性验证Fig.4 Grid independence verification

由图4可知，四种单元数网格中196万网格数量级就可满足模拟要求，网格划分示意图如图5所示。

图5 网格划分示意图Fig.5 Map of grid

3.2 模拟分析孔隙率对CO2泄漏扩散的影响

3.2.1 扩散云图

埋地管道一旦发生泄漏其输送物质的扩散介质就是土壤，中国土壤资源十分丰富且种类繁多，土壤因疏松程度的差异可分为砂土、黏土和壤土，因渗流特性的不同也可分为均匀土壤和非均匀土壤。在模拟中把土壤视作是各向同性的多孔介质，其疏松程度用孔隙率来表征[20]，具体数据如表5所示。

表5 砂土类别

由中国土壤的特点选取孔隙率φ分别为0.35、0.45、0.55、0.65，导入表4中的数据。各个孔隙率下截面XOY的 CO2浓度云图在t=100、200 s时如图6、图7所示，在时刻t=400、500 s如图8、图9所示。

图7 t=200 s时XOY截面CO2浓度云图Fig.7 The CO2 concentration cloud of XOY section when t=200 s

图8 t=400 s时地表CO2浓度云图Fig.8 The surface CO2 concentration cloud when t=400 s

(1)由图6、图7所示，在相同泄漏扩散时间内孔隙率的增大使得高浓度范围随之增大，CO2会以更短的时间从土壤扩散到地表；随着泄漏扩散时间的增加，土壤中CO2浓度云图以半圆形逐步向外扩散，由于Y轴上CO2扩散速率最大，故监测轴线上CO2的浓度便可测得CO2浓度距地表距离。

(2)从图8、图9所示，扩散到地面的CO2会以平行于地表且泄漏口为圆心按照中间高边缘低的浓度梯度呈圆形分布，在相同泄漏扩散时间内，圆形区域的直径随着孔隙率的增大而增大；随着泄漏扩散时间的增加，圆形区域中的高浓度范围也越来越大，确定浓度在Y=0轴线上的扩散距离即可求出危险区域的范围。

3.2.2 模拟结果及分析

浓度高于5%的CO2足以在短时间内致人窒息，故将5%确定为危险临界浓度[21]。处理云图数据并绘制如图10所示曲线。

绘制各孔隙率Y轴上临界浓度随时间变化曲线如图10所示，可以看到，曲线可以明显分为两段，即地面下0～1.0 m段5%浓度节点位置与孔隙率的相关性很好且与预期相符，即土壤越松散扩散速率越大，CO2浓度达到临界值的时间越短；在地面以下1.0～1.8 m段5%浓度节点位置与孔隙率的相关性并不符合预测规律。

考虑不符合规律的1.0～1.8 m段，绘制该段不同时刻Y轴上压力分布曲线如图11所示。

分析图11可得，地面以下1.0～1.8 m段在泄漏发生的初始时刻高压使CO2在孔隙中以射流的形式扩散，由于多孔介质结构的会降低射流的方向及动量，故在泄漏口处随着压力分层而出现负压区域；压力在泄漏后一段时间内变小并接近大气压，云图呈均匀扩散。

图9 t=500 s时地表CO2浓度云图Fig.9 The surface CO2 concentration cloud when t=500 s

图10 各孔隙率Y轴上5%浓度节点随时间变化曲线Fig.10 Changes of 5% concentration nodes in each porosity Y axis with time

图11 不同时刻Y轴上的压力分布图Fig.11 Distribution of pressure on Y axis at different times

图12 各孔隙率下地表5%等浓度线距Y轴线 距离与时间的关系Fig.12 Correlation between the distance of the surface 5% isoconcentration line to Y axis and diffusion time under each porosity

各孔隙率下地表5%等浓度线距Y轴线距离与时间的关系曲线如图12所示。当土壤孔隙率为0.35时，对应扩散到地表的圆形区域半径最大为1.69 m，当土壤孔隙率为0.65时，对应扩散到地表的圆形区域半径最大为1.975 m，所以以5%为危险浓度可以将已经确定泄漏位置2 m地表范围定为危险区域。

4 结论与建议

基于中国某油田超临界CO2埋地管道的输送工况，利用FLUENT软件模拟预测了超临界CO2管道发生小孔泄漏后在不同孔隙率的土壤中的扩散特性，得出以下结论。

(1)土壤松散程度越高即孔隙率越大，管道泄漏的CO2扩散速率越快，故而越松散的土壤在相同时间内形成的危险区域的面积就越大。

(2)超临界管道为高压输送管道，CO2会在泄漏初期数秒内在土壤孔隙中以射流扩散，多孔介质结构会改变射流的方向并降低射流动量，在泄漏口处由于压力分层而出现负压区域，但压力会在极短的时间内接近大气压，扩散形式变为均匀扩散。故可以将均匀扩散视为为超临界CO2埋地管道发生小孔泄漏后扩散的主要形式。

(3)以计算模拟得到的扩散规律为依据，将5%确定为CO2危险临界浓度，泄漏时间大于800 s时地表CO2浓度分布基本趋于稳定，此时可确定发生泄漏后地表的危险区域。以最大土壤孔隙率0.65为例，模拟可知其对应的危险区域半径为2 m。