在“双碳”目标和能源安全的双重驱动下，地下天然气储气库作为调节季节性和区域性供需差异的重要手段受到广泛关注^[1-2]。传统储气库通常使用天然气作为垫层气，建库成本高昂。而二氧化碳作为垫层气不仅节省大量成本，还能实现碳封存^[3-5]，具有经济和环境双重效益^[6-10]。然而，在注、采气过程中，二氧化碳与天然气的混气现象导致采出天然气物质的量降低，当前对其影响机制的研究仍不足^[11]。二氧化碳作为垫层的应用最早由Oldenburg^[12]提出，认为CO₂在接近临界压力时的密度敏感性与可压缩性有助于提高库容与采出比。近年来，关于CO₂垫层的研究持续推进，重点探讨了其在储层内的混合行为及对储气库性能的影响^[13]。相关研究表明，CO₂与天然气混气量较小，可保持较高的采收率，但不同研究中的混气扩散系数和渗流特性存在差异，难以形成统一的标准^[14-17]。国内外学者从实验和模拟两个方面对混气特性进行了广泛探讨。研究表明，孔隙结构、渗流机制及温度压力等因素共同影响有效扩散系数，这也是混气强度的关键^[18-21]。然而，现有的混气模型在统一性、气体物性差异、计算稳定性等方面仍存在不足，影响模拟结果的准确性和模型的适用性^[22-24]。因此笔者以枯竭气藏型地下储气库为研究对象，提出一种高效的数值模拟方法，研究地下储气库注采过程中的混气规律，并为二氧化碳作为垫层的天然气储气库提供合理的运行边界；考虑储气库多孔介质骨架形变对孔隙率和渗透率的影响，提出分步耦合方法解决工作气和垫层气压缩性差异导致的计算不稳定问题；提出保证全局物质的量守恒的修正方法，确保长期运行的物质的量守恒；探讨不同比例垫层气对建库、注气和采气阶段的影响，并提出初步运行计划。

1 数学模型

储气库运行方案的设计受工作气（天然气）和垫层气（二氧化碳）渗流和混气影响，涉及多孔介质的结构特征、气体在多孔介质中的传输机制、压力等因素，因而具有高度复杂性和强耦合性。

1.1 多组分输运模型

采用多孔介质中多组分可压缩气体输运方程描述组分物质的量浓度随渗流速度和气体扩散系数的变化规律，表示为

式中，φ为孔隙率；n_i为组分i的物质的量浓度，kmol·m^-3；u为渗流速度（采用达西公式计算），m·s^-1；D为多孔介质中的气体扩散系数，m²·s^-1；q_i为组分i的源项，kmol·m^-3·s^-1。

叠加所有组分的输运方程（式1）得到摩尔物质的量守恒方程为

1.2 变形储层模型

储气库多孔介质骨架结构随着气体压力变化而发生形变，从而导致孔隙率变化，因此将储气库多孔介质骨架结构随压力变化纳入模型，能够更准确地反映储气库的储气能力和渗流特性^[25]，气体压力与孔隙率的关系为

式中，φ₀为参考孔隙率；p为气体压力，Pa；p₀为参考压力，Pa；c_φ 为多孔介质的压缩系数，Pa^-1。

1.3 混气模型

在储气库运行过程中工作气与垫层气之间不可避免存在成分交换。多孔介质气体扩散系数是衡量混气在多孔介质混合和传输过程中的重要参数。根据Pisani^[18]的研究，多孔介质中的气体扩散系数与孔隙空间中的气体扩散系数和多孔介质迂曲度有关，其关系为

式中，D_porous和D_void分别为多孔介质和孔隙空间的气体扩散系数，m²·s^-1；τ为多孔介质迂曲度。

1.4 物性模型

在储气库运行压力范围（通常在5~40 MPa）内，天然气保持超临界态，而二氧化碳则在特定压力下发生气态到超临界态的转变。考虑气体物性模型以反映不同运行条件下气体物理特性，采用BWRS方程计算气体物质的量浓度，利用中高压下的非极性气体黏度模型计算气体黏度。

2 数值模拟方法

采用基于三维交错网格的有限容积法隐式离散控制方程，延迟修正格式用于离散对流项，用中心差分格式离散扩散项。由于天然气与二氧化碳的压缩性差异，强注强采工况极易导致计算不稳定，同时迭代误差会破坏离散方程的物质的量守恒性，为解决这些问题，提出分步耦合模拟方法和全局物质的量修正方法。

2.1 分步耦合模拟方法

具体步骤如下：①通过压力与物质的量浓度的关系，将物质的量守恒方程离散为压力扩散方程，隐式求解压力；②通过渗流模型计算渗流速度，显式更新组分输运方程组分的物质的量浓度。分步求解方法的难点在于有效预估压力与物质的量浓度的关系，因此从增量和全量两方面出发，提出两种格式。

2.1.1 格式A

基于SIMPLE算法的增量角度，压力与物质的量浓度之间的增量方程为

式中，p^*为压力增量项，Pa；$n_i^{*}$为组分物质的量浓度增量项，kmol/m³。

将分压定律引入到模型中，得到物质的量浓度的增量格式A为

对于复杂的气体状态方程，$\frac{\partial n_i}{\partial p_i}$的形式非常复杂甚至难以确定。因此，隐式格式无法求解，而改用显式线性插值估计，结合分压定律，利用当前迭代步和上一时步的变量计算，代入式（7）得

其中

$p^{\left(k\right)}=p\left(n_i^{\left(k\right)}\right).$

式中，k 为迭代步。

虽然格式A通过引入分压定律，简化了组分之间的相关性，得到组分物质的量浓度与压力之间的线性关系，但该线性关系在组分变化较小时近似成立，仅使用格式A计算二氧化碳垫层的地下储气库无法完全保证计算稳定性。

2.1.2 格式B

基于SIMPLER算法的全量角度，压力与物质的量浓度之间的全量方程为

显式线性插值仍然可用。与格式A不同，格式B中直接采用当前迭代步变量计算，表示为

由于 n ^（k）/ p^（k）恒为正数，格式 B 可以保证离散方程组系数矩阵严格对角占优，从而保证计算稳定性。但是格式 B 也存在缺点，其对$\frac{\partial n_i}{\partial p_i}$的近似（ 式（12））较格式 A（式（9））精度低。因此，将格式 A 在组分变化较小时的高精度与格式 B 的绝对稳定结合，提出一种新的兼顾稳定性与精度的方法，即组分变化大于 1%的网格及其相邻网格采用格式 B，保证计算稳定性，其他网格采用格式 A，以准确描述压力与组分物质的量浓度之间的关系，表示为

式中，H_A和H_B分别为格式A和B的布尔变量（H_A+H_B=1）。

2.2 全局物质的量修正方法

时步形式的组分总物质的量变化和时步内注入或采出气体总物质的量的全局物质的量守恒方程表示为

其中

$N_i={\int }_{\mathrm{\Omega }}\phi n_i\mathrm{d}V=\sum _{k=1}^{N_{\text{cell}}}\phi n_i\mathrm{\Delta }{V}_k.$

式中，Ω 为储气库计算域; N_cell 为网格单元总数; ΔV_k 为第 k 个网格单元体积; N_i 为组分 i 在储气库的总物质的量，kmol; N_well 为注采井数量; $\mathrm{\Delta }t\sum _{j=1}^{N_{\text{well}}}{\left(q_i^{（n）}\right)}_j$为时步内注入或采出气体的总物质的量，kmol; ΔV 为模拟网格体积，m ³ 。

由于上一时步已经修正了物质的量$N_i^{（n）}$，只需修正当前时步的物质的量$N_i^{（n+1）}$。参考上节方法，在增量角度下引入压力修正倍数${\beta }_{pi}^{*}$，带入全局物质的量守恒方程式（14），得

式中，${\tilde{N}}_i^{（n+1）}$为当前时步计算得到的总物质的量，kmol。

然而，基于增量的全局物质的量修正方法存在一些局限性，不适合在实际计算中使用。有两个问题:

（1）在关井阶段，$\mathrm{\Delta }t\sum _{j=1}^{N_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}}{\left(q_i^{（n）}\right)}_j=0$，导致修正倍数${\beta }_{pi}^{*}$为-1，修正方法失效;

（2）由于储层孔隙内气体物质的量远远大于一个时步内的注、采气量，可能会出现$\left|{\tilde{N}}_i^{（n+1）}-N_i^{（n）}\right|\gg \left|\mathrm{\Delta }t\sum _{j=1}^{N_{\text{well}}}{\left(q_i^{（n）}\right)}_j\right|$，这导致压力修正倍数${\beta }_{pi}^{*}$不接近 0，与预期相悖。

由于基于增量的方法有局限性，通过直接引入压力修正倍数${\beta }_{pi}^{*}$到${\tilde{N}}_i^{（n+1）}$，得到全量的全局物质的量修正式，表示为

其中

${\beta }_{pi}^{*}=\frac{\mathrm{\Delta }t\sum _{j=1}^{N_{\text{well}}}{\left(q_i^{\left(n\right)}\right)}_j+N_i^{\left(n\right)}-{\tilde{N}}_i^{(n+1)}}{{\tilde{N}}_i^{(n+1)}}.$

由于${\tilde{N}}_i^{（n+1）}$和$N_i^{（n）}$的数量级相近，可保证${\beta }_{pi}^{*}$接近 0，推荐使用这种方法修正全局物质的量。全局物质的量修正方法解决了传统方法在计算过程中多组分物质的量不守恒的问题，保证了储气库长期运行的各组分物质的量守恒。

3 数值模拟方法分析

3.1 分步求解方法的计算稳定分析

由于注气阶段计算更不稳定，采用注气算例评估求解方法的计算稳定性。初始储气库压力设置为 10 MPa，垫层气的注气量为 5 kmol·s^-1。在相同迭代次数下以时间步长表征效率，用求解过程中仍可保持收敛的最大松弛系数表征稳定性裕度。迭代次数越少且最大松弛系数越大，表示该方法在给定工况下更高效，稳定性也更好。图1 为不同时步的迭代次数和最大松弛系数，分步方法和全耦合方法的物质的量浓度迭代值。由图1 可知:当时间步长为 7200 s 时，全耦合方法无法收敛，而分步方法仍然可以收敛（图1（a））; 分步方法在迭代次数和松弛系数方面均优于全耦合方法，这表明分步方法在稳定性方面有显著提升。在保证相同迭代次数情况下，分步方法可以将时间步长（效率）扩大到 20 倍。为了比较两种方法在收敛性方面的差异，将时间步长设为 360 s，计算第一个时步的井底处组分物质的量浓度迭代值，可以看到分步方法的最终计算稳定性更好（图1（b））。

3.2 全局物质的量修正方法的物质的量守恒分析

为了验证所提出修正方法的优点，设置算例：①0~30 d，以注气流量5 kmol·s^-1注入垫层气（二氧化碳）；②31~60 d，以注气流量5 kmol·s^-1注入工作气（天然气）；③61~90 d，关井；④91~120 d，以采气流量2 kmol·s^-1采气。

模拟得到传统方法的全局物质的量误差，见图2（a），在整个运行周期内，工作气和垫层气的最大物质的量误差分别为2754 kmol（约50 t）、1449 kmol（约64 t）。采用更严格的全局修正方法重新模拟算例，见图2（b），修正后的物质的量误差数量级为10^-6 kmol，远小于传统方法10³ kmol量级的全局物质的量误差，可以严格保证全局物质的量守恒。

4 混气规律

4.1 算例设置

地下储气库储层模型的尺寸为510 m×510 m×210 m，储层温度保持在350 K不变，孔隙率为10%，渗透率为100×10^-3 μm²。在储层底部中心位置（255 m，255 m，35 m）注入二氧化碳垫层气，顶部中心位置（255 m，255 m，185 m）注采工作气，储气库初始压力设为5 MPa，运行压力控制在15~30 MPa。

在初始条件下，储层气体仅含有甲烷，初始压力为5 MPa，储气库中初始含有天然气总量为9.91×10⁶ kmol。储气库的最低运行压力为15 MPa，若储层全部为天然气，天然气总量为31.08×10⁶ kmol。选取4种不同的二氧化碳垫层气比例（60%、40%、20%、0%）进行模拟，系统分析不同垫层气比例对建库、注气与采气阶段的影响规律。

在建库、注气和采气阶段，分别需要关注垫层气的最高位置、最低位置以及垫层气的最高位置是否高于井底。为了便于分析，引入两个参数：二氧化碳垫层气的最大和最小相对高度，其定义如下：

式中，γ_max和γ_min分别为最大和最小相对高度；H_max和H_min分别为气体在储气库中的最大和最小高度；H为储气库厚度，m。

在开展网格加密与收敛性分析时，设置4档网格：粗（13125单元）、细（28749单元）、较细（54621单元）与极细（111909单元）。在二氧化碳垫层气比例40%的工况下，计算采气阶段采出气体纯度（甲烷体积分数）见图3。以极细网格为参考，评估各档网格的相对最大误差和相对平均误差：粗、细和较细网格相对最大误差分别为0.765%、0.763%、0.464%，相对平均误差分别为0.146%、0.127%、0.087%。这些误差均低于1%，表明该量对网格不敏感。综合精度与计算量，选用较细网格（54621单元）作为后续计算的基准网格。

4.2 建库阶段混气

在建库阶段，需要分别注入二氧化碳和天然气垫层气将储气库压力提升至15 MPa，分析不同二氧化碳垫层气比例对建库阶段的影响，建库阶段持续180 d，之后关井180 d。图4中小图上半部为体积分数50%的二氧化碳等值面，小图下半部为横截面图，红色区域为二氧化碳，蓝色区域为天然气。建库过程可分为4个阶段，阶段1：锥形分布形成。在建库初期，二氧化碳注气流速较大，由于其密度显著高于天然气，二氧化碳迅速下沉并在储气库底部聚集，形成锥形分布。阶段2：凸起分布形成。随着注气流速进一步降低，二氧化碳垫层气的下沉趋势逐渐增大；垫层气开始呈现凸起的分布形态，随着储气库内部压力上升，工作气与垫层气的密度差加大，重力作用增强，垫层气的最大相对高度随之减小。阶段3：分层分布形成。在注气停止后，二氧化碳垫层气仅受重力驱动，逐渐在储气库中呈现出更加明显的分层现象；二氧化碳开始完全占据下层，最大相对高度迅速下降，分层分布结构逐步成型。阶段4：完全分层状态。储气库内气体达到动态平衡后，二氧化碳与天然气之间形成明显的分层结构，储气库内部的气体状态趋于稳定，建库阶段也随之结束。

4.3 注气阶段混气

在注气阶段，通过注入工作气将储气库压力由15 MPa提升至30 MPa，注气持续180 d，注气流量为1.61 kmol·s^-1，之后关井45 d。注气过程可分为3个阶段:阶段1为凹陷分布形成，由于天然气从储气库顶部注入，天然气向下流动并挤压二氧化碳垫层气，形成明显的凹陷分布；阶段2为凹陷分布减弱，随着注气流速减少，天然气对垫层气的挤压作用减弱，凹陷形态逐渐趋于平缓；阶段3为分层分布形成，在注气停止后，垫层气仅受重力影响，逐渐从凹陷形态转变为分层分布，垫层气在储气库底部稳定聚集，注气阶段也随之结束（图5）。

图6（a）中垫层气比例直接影响注气阶段各阶段的持续时间。随着垫层气比例减小，混气界面远离注气井井底，渗流速度对界面的影响减小，导致阶段1和阶段2的持续时间均缩短。如图6（b）所示，比较注气阶段混气量变化，与建库相似，随着垫层气比例增加，混气量显著增加。在注气初期，混气量迅速增加，随后趋于平稳；这是由于初期注气速度较快，混气界面运移速度大，混气加剧，而后期注气速度减缓，混气量趋于稳定。最终混气中的天然气物质的量分别为3.02×10⁶、3.88×10⁶、3.72×10⁶ kmol，远小于建库阶段结束时天然气总量，混气中的天然气分别占总天然气的8.05%、8.89%、7.46%。

4.4 采气阶段混气

在采气阶段通过采出工作气将储气库压力由30 MPa降至15 MPa，由于采气流量过大会直接采出垫层气，因此，采气分为两轮，每轮采气时间持续45 d，第一轮采气流量为4.845 kmol·s^-1（满足冬季高供应需求），第二轮采气流量为1.615 kmol·s^-1（满足春季低供应需求），关井时间为22.5 d，两轮采气共135 d。采气过程可分为两个阶段，阶段1：锥形分布形成。由于从储气库顶部采出气体，二氧化碳垫层气开始向上移动，受到重力和采气流速的共同影响，形成锥形分布；这种分布状态下，垫层气集中在井底附近，容易随着天然气一同被采出，对气体纯度产生不利影响。阶段2：分层分布逐步恢复。随着采气停止，垫层气的锥形分布逐步减弱，二氧化碳重新向储气库底部聚集，逐渐恢复分层状态；然而，由于关井时间较短，采气过程中的一部分锥形分布尚未完全消失，这使得在第二轮采气开始时混气风险较大（图7）。

如图8（a）所示，垫层气比例直接影响采气阶段中各阶段的持续时间，由于垫层气比例减小，采气井井底距离混气界面的距离增大，渗流速度对界面的影响减小，从而垫层气不容易被采出。但值得注意的是，由于二氧化碳压缩性优于天然气，第一轮关井结束后，垫层气比例较高的储层压力大，更容易达到“分层结构”。同时，比例60%（15 MPa）的垫层气在采气过程中最大相对高度高于采气井井底，会导致直接采出垫层气，这是高比例垫层气的局限性，不利于高强度采气。因此，为了提高采气可靠性，应适当降低垫层气比例，采用垫层气比例40%（15 MPa）的垫层气更合适。比较采气阶段的混气量变化，采气阶段有一部分垫层气被采出，见图8（b）。最终混气中的天然气物质的量分别为3.55×10⁶、3.54×10⁶、3.15×10⁶ kmol，远小于建库阶段结束时天然气总量，混气中的天然气分别占总天然气的28.57%、19.00%、12.71%。由于高比例垫层气更易采出混气，导致比例40%和60%垫层气的混气量下降。采气阶段结束时的混气量略高于建库阶段结束的混气量，这说明循环注采过程中二氧化碳垫层气会不断与天然气发生混气。

比较采出气中二氧化碳的体积分数（图9），可以得到不同比例二氧化碳垫层气的高比例垫层气，垫层气比例60%在采气过程中会导致较高的二氧化碳体积分数（15.74%）。这意味着在高强度采气时，二氧化碳容易被直接采出，不利于天然气的纯度。中等比例垫层气，40%的垫层气比例在采气过程中表现出较为平衡的二氧化碳体积分数（5.29%），这种比例在保证采气效率同时，也能较好地控制二氧化碳混入。低比例垫层气，20%的垫层气比例在采气过程中二氧化碳体积分数最低（1.03%），有助于保持采出气的高纯度，但会影响储气库的建库投资。

综上，二氧化碳在储层底部能够形成稳定分层，其抑制天然气扩散的机制主要体现在3方面：①由于二氧化碳密度远高于天然气，在重力作用下快速下沉并聚集在储层底部，形成天然气难以穿透的“高密度屏障层”；②分层结构显著减小了天然气与二氧化碳的接触界面面积，从而降低了分子扩散驱动力；③二氧化碳黏度较大，使得混气界面扰动受到阻碍，减缓了天然气向下部渗透。这些作用共同抑制了天然气扩散，保证了储气库运行的稳定性。

5 结论

（1）模拟模型充分考虑了储层和气体性质随压力变化情况，特别是针对变形储层模型，新提出的更为稳定的分步耦合方法从增量和全量角度得到了兼顾精度和稳定性的计算格式及适用于多组分的全局物质的量修正方法。

（2）模拟建库、注气和采气3个阶段运行周期，比例40%（15 MPa）的二氧化碳是最理想的选择，能够保持储气库性能平衡稳定，二氧化碳作为垫层气体应用于天然气地下储气库是可行的且具有显著优势，二氧化碳与天然气之间的密度差异使二氧化碳能够在储气库底部形成稳定的分层结构，从而有效减少天然气垫层气使用量，提高储气库的经济效益。

参考文献