含杂质CO2管道输送泄漏扩散的数值模拟

陈 兵，康庆华，肖红亮

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

随着社会科学技术的进步和人们生活水平的提高，环保问题受到广泛的关注，这不仅仅是我国面临的一个严峻问题，也是全球亟待解决的问题。在这样的国际大环境下，碳捕获、利用与封存(Carbon Capture,Utilization and Storage,CCUS)这一项新兴的、具有大规模CO2减排潜力和环境效益的技术[1-7]应运而生。CO2的管道输送是CCUS系统的中间环节，承担着将CO2运输到工作点的任务，是实现CCUS大规模工程化应用的关键技术环节[8]。根据国家职业安全卫生研究所(NIOSH)提供的数据显示：CO2浓度达40 000 ppm(1 ppm=10-6)时，将直接危及生命[8]。尽管CO2管道输送泄漏时无燃烧风险，但仍需将泄漏点周围的超高压风险当作泄漏风险来考虑，尤其是在发生预料之外的管道破裂时[9-10]。

管道破裂时会造成重大压力损失和管道流体的损失[11]，在压力释放期间，减压会使CO2气体较之原先的体积以很大的膨胀系数扩散至大气中，而且减压会导致其发生相变[12-14]，随着部分高蒸汽压CO2汽化为气体，剩下的CO2则变为液体和固体颗粒。通常在减压结束后这种膨胀力会驱使这些CO2固体颗粒以很高的速度弹出，这会给人员和设备造成严重的危害[15-16]。此外，减压致使CO2及其杂质[17](如H2S)从管道中释放出去，在释放位置处，如果不与空气迅速混合，那么在管道破裂位置相对较近范围内的环境中则具有毒性，当管道内所有毒性物质在失压后全部脱离出来时，若不采取适当的措施，这些物质可能会污染并沉淀在所在的区域，给暴露于该环境中的人造成危害，且在合适的温度和压力下，也会形成水合物[18-19]。CO2泄漏后直接排放在基质上会出现干冰“喷砂效应”，在失压后，沙子、土壤、和其他固体碎片与生成的干冰相结合形成研磨能力更大的颗粒，增强了对泄漏位置和管道工程的腐蚀效应[20]。CO2输送管道破裂时，CO2流进入大气后会扩散，在气流消散前有不同的泄漏阶段和主导作用，所以研究泄漏口附近的压力随时间的变化规律、压力分布以及泄漏阶段、持续时间和浓度等级是非常重要的，可以为管道放空CO2流的相态变化、干冰或者水合物的形成、CO2云扩散和管道完整性研究提供一定的参考价值。为此，本文选取某油田CO2管道输送工况作为研究对象进行了动态数值模拟，得出的结论可为相关研究提供一定的技术支撑。

1 某油田管道输送气源指标及成分

经调研，某油田含杂质CO2管道输送气源指标见表1，CO2气流成分见表2。

表1 某油田含杂质CO2管道输送气源指标

表2 某油田管道输送CO2气流成分

2 超临界CO2输送管道泄漏点计算流体动力学(CFD)模拟研究

在设计和运行CO2管道输送时，CO2泄漏模拟是定量风险分析的一部分，而分析随后的CO2扩散规律将有助于对这类风险的影响因素进行评估。为了更深一步地了解CO2的扩散机理，本文重点研究CO2泄漏和扩散过程的特性，包括近场的压力场和浓度场等，即选用Fluent软件中的组元输运模型，将网格“外疏里密”，建立结构化网格，并利用某油田CO2管道输送工况作为研究管道泄漏时CO2气体扩散特性的入口边界条件，选用湍流模型中的k-ε模型和真实气体状态P-R方程进行求解和计算，得到相应的CO2浓度扩散云图，最后经数据处理后得到模拟结果。

2.1 Fluent软件简介

Fluent是世界领先的计算流体动力学(Computational Fluid Dynamics,CFD)软件，在流体建模中被广泛应用。该软件针对各种复杂流动的物理现象，采用不同的离散格式和数值方法，以期在特定的流域内使计算速度、稳定性和精度等方面达到最佳组合，从而高效地解决各个领域的复杂流体流动计算问题[21]。Fluent软件主要基于N-S方程和有限体积法，适用于复杂的流体流动计算，其组元输运模型可用于计算气体的扩散过程，且可操作性较好，因此本文选用该软件对超临界CO2输送管道泄漏过程进行数值模拟。

2.2 数值建模

超临界CO2输送管道泄漏过程主要分为射流阶段和扩散阶段两个阶段，建立的数值模型见图1。

2.3 基本控制方程和状态方程

CO2在泄漏或放空过程中，其流体流动仍遵循动量、能量、质量三大守恒定律，并由此推导出在直角坐标系下CO2流体适用的最基本的连续方程、运动方程和能量方程如下。

(1) 连续方程：

(1)

式中：ρ为气体密度；u、v、ω分别为x、y、z轴的分量，该方程本身守恒。

(2) 运动方程：

(2)

式中：pij为定点压强的张量表示法，xi、xj为坐标矢量的张量表示形式，i、j表示自由指标，分别取1、2。

一般对于气体流动而言，质量力fi可以忽略，即认为fi=0，再将牛顿黏性应力公式代入，可得到N-S方程(纳维斯-托克斯方程)：

(3)

(4)

式中：ui、uj为速度矢量的张量表示形式。

(3) 能量方程：

(5)

该式表明，单位流体微团能量的增量等于单位时间流体面力所做功加上它所获取的能量。

(4) 状态方程：

本文在计算CO2流体物性参数和CO2流体扩散过程时均选用P-R方程：

(6)

式中：ω为偏心因子,ω=-1;μ为动力黏性系数,μ=2；a、b为物质特有常数，与临界参数有关，其关系式如下：

式中：Tc、Pc为临界温度和临界压力，R为通用气体常数。

2.4 网格划分

建立好满足要求的几何模型后，需要进行网格划分。网格划分是CFD 流体力学计算中最为关键的一个环节，网格质量的优劣对最终求解结果的准确度有直接的影响。若网格划分过疏，会使计算结果误差过大或者所得结果根本无法收敛；若网格划分过密，尽管会提高计算结果的准确度，但计算精度增幅不大，反而会极大地延长计算时间，造成计算成本的增加。

某油田CO2管道输送工况中CO2浓度为5 000 ppm的扩散区域很大，考虑到气体在大气中的扩散呈三维对称，可将三维模型简化成在主风向上垂直于地面的一个剖面，并使用GAMBIT软件进行网格划分(见图2)。在进行网格划分时，为了保证计算速度和计算精度，将建立结构化网格，即选用四边形网格划分区域，由于泄漏口射流速度较大，因此整个区域网格由外层向内层逐步细化，由外到内共分为5层。左边设置为风速进口(VELOCITY_INLET)，上边和右边设置为压力出口(PRESSURE_OUTLET)，底边除进口外设置为墙(WALL)。

图2 网格划分示意图Fig.2 Grid generation diagram

2.5 基本设定

启动非定常求解器，重力加速度为9.81 m/s2，方向为y轴方向，压力单位为kPa，环境温度为280.95℃，大气压为101.325 kPa，开启能量方程，选用湍流模型中的k-ε模型，利用水力直径和湍流强度，开启组元输运模型，组分选用某油田CO2管道输送组分(但这里仅以CO2为研究对象，故只追踪CO2)。

2.6 入口边界条件

考虑到输送管道中CO2流体为超临界/密相状态，泄漏后必然会发生相变，因此计算得到孔口处流量、温度、压力随时间的变化规律，并用C语言编写UDF导入Fluent软件作为入口边界条件。

本文使用ASPEN HYSYS V8.4创建Depressuring-Dynamics单元，以储罐表示管道，阀门表示孔口，对某油田CO2管道输送工况进行动态数值模拟。该工况为CO2流体超临界/密相输送，绝热过程，输送压力为10 MPa，初始温度为60℃，最优管径为φ273 mm×10.0 mm。模拟得到孔口处流量、压力和温度随时间的变化曲线，详见图3至图5。

图3 孔口处流量随时间的变化曲线Fig.3 Curve of flow at the pipe orifice with time

图4 孔口处压力随时间的变化曲线Fig.4 Curve of pressure at the pipe orifice with time

图5 孔口处温度随时间的变化曲线Fig.5 Curve of temperature at the pipe orifice with time

经拟合，可得到孔口处流量(Q)、压力(P)和温度(T)与时间(t)的关系式如下：

Q=14.2+1 263.4e(-t/92.3);

P=156.4+9 895.0e(-t/119.7);

T=-93.3-144.6e(-t/307.3).

2.7 无关性检测

由于模型所涉及的扩散区域过大，因此需要对模型无关性区域进行检测，即对CO2浓度波动进行评价。初始场域定为200 m×1 800 m，场域过大会使网格划分极为复杂，为了避免网格奇异、过大纵横比，网格划分基本遵循“外疏里密”的原则，泄漏口位置根据左侧场域变化进行调整。经模拟验证，左侧CO2扩散变化范围为120 m左右，在扩散时间为1 000 s时，上方130 m处、右侧750 m处CO2扩散浓度不再发生变化，因此模型扩散区域定为150 m×1 000 m，泄漏口离左侧边界为200 m。

2.8 模拟结果与分析

利用Fluent软件的模拟计算中，由于未考虑干冰的形成，在离开泄漏口，温度会很快接近大气温度，高压CO2射流决定了其扩散形貌，所以本次模拟不考虑速度场、温度场，主要分析浓度场和近场的压力场[22]。

2.8.1 CO2浓度的变化

本文利用Fluent软件模拟得到的某油田管道泄漏CO2浓度云图，见图6。

图6 某油田管道泄漏CO2浓度云图Fig.6 Cloud map of leaking CO2 concentration of an oil field

由图6可以看出：

(1) 在泄漏最初阶段，CO2浓度云图并未发生倾斜，可见在此阶段初始泄漏高压起主导作用，对应模型中的射流阶段。

(2) 随着时间的推移，CO2扩散高度的增加，动能不断减小，CO2浓度云图开始向右倾斜，即风速开始起主导作用，大气湍流作用增强，云团变化变得复杂化，对应模型中的扩散阶段。

(3) CO2扩散以顺风风向为主，在泄漏源停止泄漏之前(800 s之前)，沿顺风方向越远，CO2扩散浓度越低，泄漏源停止泄漏之后，CO2会继续沿顺风方向运移，其浓度基本保持在10-2数量级(比实际略大)。

2.8.2 压力的变化

本文利用Tecplot软件对Case和Data文件进行后处理，将泄漏口附近3个时刻的压力云图与Tecplot导出的压力等值线图、流线图集于一图，得到3个时刻泄漏口附近的压力云图，见图7至图9。

图7 10 s时泄漏口附近的压力云图局部放大图Fig.7 Partial enlargement of pressure cloud near the leak hole at 10 s

图8 50 s时泄漏口附近的压力云图局部放大图Fig.8 Partial enlargement of pressure cloud near the leak hole at 50 s

图9 200 s时泄漏口附近的压力云图局部放大图Fig.9 Partial enlargement of pressure cloud near the leak hole at 200 s

根据上述给出的3个时刻泄漏口附近的压力云图，现将泄漏口上方0.1 m处的压力导出，得到局部压力变化曲线，见图10。

图10 局部压力变化曲线Fig.10 Curve of local pressure

由图7至图10可以看出：

(1) 高压CO2流体释放到大气中会形成自由射流，在泄漏口处压力很高，有明显的压力分层，并形成负压区域，50 s时压力达到最高，300 s后趋于稳定，即泄漏口压力趋于大气压。

(2) 泄漏口附近压力场几乎呈对称分布，流线呈辐射状分布，如同CO2浓度云图在泄漏口附近并未发生倾斜，此阶段初始泄漏高压起主导作用，对应模型中的射流阶段。

(3) 随着时间的推移，泄漏口上方0.1 m处的压力分布变得复杂化，负压区区域由大变小而后又变大，这是由于高压CO2流体会在泄漏初期快速膨胀又迅速减小；当压力快速上升进入减压阶段时，减压会使CO2气体较之原先的体积以很大的膨胀系数扩散至大气中，具体膨胀程度由管道工况压力与温度协同作用决定，此阶段对应模型中的扩散阶段。

(4) 减压时使负压区区域变大，时间持续较长，使得扩散出来的高压CO2流体与空气更易混合且复杂化。结合图6可知，在300 s之内，CO2浓度等级特别高(从释放位置起，会持续一段距离)，减压致使CO2及其毒性物质(如H2S)从管道中全部脱离出来，使管道破裂位置一定范围内的环境具有较大的毒性，会对意外进入该范围内的人类或动物造成严重的危害。为了安全，应在危险浓度区域边界放置警告牌，即：在橙色区域附近应放置职业浓度警告牌，当人或动物进入该范围时可能会出现危险；在红色区域附近应放置危险浓度警告牌[23]，当人或动物进入该范围时会立即出现危险。

3 结 论

(1) 本文利用Fluent软件，选用组分输运模型对含杂质CO2管道输送泄漏进行了数值模拟计算，得到了CO2浓度和压力云图。

(2) 对应模型中的泄漏阶段划分，由CO2浓度云图可以看出CO2管道输送泄漏主要分为两个阶段：①射流阶段,该阶段输送压力起主导作用，范围小，持续时间较短；②扩散阶段,该阶段风速、大气湍流起主导作用，范围大，持续时间较长。

(3) 在泄漏口附近会形成压力膨胀层，压力值迅速上升又下降，并由高速气流造成负压区，该区域的流线和压力场呈对称分布。

(4) 通过数值模拟得到的含杂质CO2管道输送泄漏扩散规律基本符合实际，但由于Fluent软件的可操作性较大，模拟得到的具体数值与实际略有差别，仍需参照实验值并设定具体函数进一步研究CO2管道输送泄漏过程特性。

(5) 为了安全，根据CO2扩散特性和压力响应应在泄漏点一定范围内(警示区域、危险区域)放置安全警示牌，即：当人或动物进入该范围时可能会出现危险时放置职业浓度警告牌；当人或动物进入该范围时会立即出现危险时放置危险浓度警告牌。在控制好危险范围的同时，更要做好气源地的除杂工作，严格限制气源中H2S、CO、SO2等杂质的含量。