气体在双重多孔介质中运移的分数阶导数建模研究

胡 韬 吴本霞 余星吉 靳 勇

(南京市江南小化工集中整治工作现场指挥部，江苏 南京 210038)

在我国经济快速发展和工业化进程不断前进中，全球的能源如石油、煤和天然气等占据重要地位，也是国民经济发展的重要基础[1]。随着油气工业的快速发展以及油气资源的不断消耗，非常规油气资源如致密油、致密砂岩气、煤层气和页岩油气等正逐步取代常规油气资源[2]。目前，随着自然资源日益减少，环境污染与能源开采问题是能源可持续发展的两大关键问题[3]，越来越受到整个社会的高度重视。众所周知，双重多孔介质系统如天然土壤、裂缝岩石、油气储藏、煤气层等，其裂缝和基质孔隙自身结构具有不同的几何特性。大多数油气藏结构也是由多孔介质与树状分叉网络镶嵌而成的双重多孔介质结构[4]。双重多孔介质中的流体流动广泛应用于水文、地质、石油工程、农业工程、环境工程等自然界和工业生产中，如地下水源污染预测、二氧化碳地下埋存和油气资源开采输运等，并越来越受到更多领域研究学者的高度关注和重视[5]。就当前研究现状而言，研究较多的是地下水、污染物及溶质的运移过程，可为水污染预测及治理提供理论依据[6]。而在油气开采运输过程中，油气储藏存在分布范围较广、渗透性能较低和开采难度较大等特点，常用的注水驱替技术会损坏套管、发生水窜效应且采收率较低。而气驱法最早源于1958年美国进行的注CO2混相驱项目；1996 年，美国又在圣胡安盆地 Burlington Allison试验区进行了CO2驱替煤层气现场试验，并达到了提高采收率的效果。1963年我国将该技术引入到大庆油田；2003 年我国在山西的沁水盆地南部进行了CO2驱替煤层气的单井试验。此外，高等[7]利用不同高温高压实验系统研究了CO2与原油相互作用的机理；黄等[8]建立了CO2驱替原油的多相渗流模型。总之，将CO2注入到油气藏中，不仅可以实现CO2的地质封存，减少大气中CO2的排放量，也可以降低油的粘度且提高油和气的产量，从而达到提高油气采收率和改善环境的双赢目的[9]。因此本文主要关注气体在复杂油气藏中的运移过程，搞清气体运移规律，对油气勘探、开发都有重要的指导意义。

双重多孔介质中的气体运移过程由通过大尺度裂缝的渗流流动和基质内部的扩散运动两种不同的物理机制构成[10]。Bear[11]第一个将对流—扩散方程(Convection-Diffusion Equation, CDE)模型用于描述宏观多孔介质中溶质的运移过程。而在非均质土壤或储藏结构中，由于流体自身特有的物理属性和赋存压力以及外界环境因素的干扰，使得流体流动出现反常扩散现象，即所测得的穿透曲线(Breakthrough Curves, BTCs)通常表现出早到达和明显的长拖尾现象[12]，CDE模型无法充分描述该现象。本文引入了分数阶导数模型的应用，分数阶微分方程具有记忆性或非局部性，在工程、信息处理、物理、金融、流体等领域的应用越来越广泛，比整数阶微分方程能更好地模拟自然物理过程和动力系统过程[13]。2000年，Benson等[14]最早将分数阶导数模型应用于地下水和溶质运移中，模拟Cape Cod实验场地中溴离子的反常扩散现象；Meerschaert等[15]通过研究发现，初始阶段污染物在土壤中运移的BTCs会出现提前穿透，而在运移后期会出现长拖尾现象。此外，Ali等[16]建立了用于描述致密介质中气体运移规律的时间分数阶对流—扩散模型。El Amin等[17]提出了空间分数阶导数模型用于探讨气体运移过程。总之，目前研究高效开采油气资源及突破技术瓶颈是学者们急需解决的科学和工程问题，而将力学、数学与石油工程、地球科学等学科实现深度交叉融合，可以更加有效地推动非常规油气资源的勘探开发。

1 模型与方法介绍

1.1 分数阶对流—扩散模型

从以往的研究中可知，对流—扩散方程可以准确地揭示流体在均匀多孔介质中的运移机理和时间、空间对运移过程的影响，同时还可以描述流体的反应扩散、质量分布和热量运输等物理现象，且形式简单便于计算。通常用于描述均匀介质中的粒子运移问题的基本方程即经典对流—扩散方程，它的表达式如式1所示。

(1)

式1中C为浓度，v为平均流速，K为扩散系数。

在自然界和实际工程应用中，由于复杂多孔介质(土壤、裂隙岩体)结构和物理性质的非均匀性导致了反常扩散现象的出现，利用CDE模型描述气体运移问题时存在不足。因为流体或粒子的运移过程是随着时间、空间发生变化的，并不是恒定不变的。大量的实验观测和数值模拟结果也显示气体运移的浓度穿透曲线具有时间方向上的拖尾现象，和空间上的非高斯分布特征，而分数阶导数算子具有历史记忆性和空间依赖性，因此引入了分数阶对流—扩散方程模型。本文主要开展了双重多孔介质中气体运移机理、运移模型、简化求解、模型参数、模型应用以及运移规律等研究。

分数阶对流—扩散方程模型的表达式如下：

(ⅰ) 时间分数阶对流—扩散模型

(2)

(ⅱ) 时空分数阶对流—扩散模型

(3)

这里α为时间分数阶阶数，β为空间分数阶阶数。

分数阶导数的定义由时间或空间的非局部卷积算子表示，刻画粒子反常扩散过程。常用的分数阶定义包括Riemanm-Liouville、Caputo和Gründwald-Letnikov。下面简要介绍几种常用的分数阶定义。

(ⅰ) Riemanm-Liouville定义[18]

(4)

和

(5)

(ⅱ) Caputo定义

(6)

(ⅲ) Grünwald-Letnikov定义

——法国农业部称：预计2018年法国葡萄酒总产量为46.1亿升，比2017年增长27%，超出过去5年平均产量达6%。在波尔多，尽管部分酒庄受到冰雹和霜霉病的影响，但预计产量还将达到历史正常水平。

(7)

(8)

1.2 数值模拟方法

通常分数阶微分方程的解析解含有特殊函数或者复杂的级数，难以获得方程的解析解，因而数值求解方法引起了学者们的广泛关注。多种求解对流—扩散方程的数值方法，如有限差分法(Finite Difference Method)、有限体积法(Finite Volume Method)、有限元方法(Finite Element Method)和无网格方法(Meshless Method)等被利用[19-21]。

本文主要采用的是较为灵活、简便的数值计算方法，即有限差分法求解分数阶微分方程，且求解分数阶偏微分方程具有适用范围较广、精度较高和稳定性较好等特点。Liu等[22]建立了一种显式差分格式求解时间分数阶扩散方程，并证明了方法的稳定性和收敛性。Meerschaert和Tadjeran[23]利用有限差分法有效地获得了分数阶微分方程的数值解。此外，Chen等[24]也较好地将有限差分法应用于求解分数阶对流—扩散方程中，且方法的稳定性和收敛性被有效地证明。丁志清[25]利用有限差分法求解了空间分数阶导数方程，并进行了收敛性分析。本文采用隐式有限差分法建立数值求解方案，该方法的稳定性和收敛性已被证明，这里不再赘述，感兴趣的读者请参考文献[26-28]。

2 数值模拟结果分析

近年来的研究发现，整数阶导数的非线性和变系数模型能够描述某些特殊情况下的反常扩散过程，但依然存在适用面窄、参数物理意义不清晰和长时间历程模拟不准确等缺点[18]。本文引入分数阶导数模型，并对其进行数值求解。将时间分数阶导数项进行离散，可以写成如下式9的形式[29]。

(9)

2.1 算例1

假设给定入口边界处一个连续点源Cin且持续向模拟样品入口端通入气体，并在出口边界定义了对流边界条件，气体浓度曲线呈不断上升趋势直至逐渐达到稳定。给出初始浓度值为0，且出口边界条件设置为自由边界条件。

通过数值模拟方法对气体粒子在一定介质区域内的运移过程进行数值模拟，如图1所示。

图1 算例1：已知初始条件和边界条件下，给定不同分数阶阶数α求解时间分数阶对流—扩散方程(FCDE)模型在一定时间内模拟获得的气体粒子的几种无量纲浓度穿透曲线。模型参数如下：v=0.2，K=0.05，T=240和L=120。

在非均匀介质结构中，分数阶导数模型是模拟粒子反常扩散的有效工具。图1显示，时间分数阶对流—扩散模型模拟结果中的分数阶阶数α越小，粒子运移过程越缓慢且浓度达到最大值所需的时间越长。结合模拟结果推断，介质结构是影响浓度穿透曲线变化的主导因素，分数阶模型中的阶数有效反映了介质结构对粒子运移的作用。由于油气储藏基质具有较强的非均质性结构，双重介质结构中存在死孔隙，或部分运移通道被阻塞、填充，这些都阻碍了气体粒子的运移过程。从图1可以看出，随着分数阶模型中的时间阶数α变小，则显示气体运移过程越慢。而α=1时，经典对流—扩散模型在达到最大浓度之前以最快的增长趋势变化，直至浓度达到峰值。

2.2 算例2

假定在入口边界处给定一初始值为C0，而设置入口边界条件为0，同时设置出口边界为自由边界条件。

分数阶导数模型的数值模拟结果展示如下。

图2显示在时间分数阶导数模型中，通过数值方法求得几种不同时间阶数α下的浓度穿透曲线。在分数阶对流—扩散模型中，当0<α<1时，分数阶阶数α对气体运移过程具有强有力的影响。在固定的控制位置1/5L处，随着阶数α的减小，浓度峰值越小，并且反常扩散会出现缓慢的衰减。图2中可以看出浓度穿透曲线的尾部斜率随阶数α的改变而改变，阶数越小，则幂律拖尾现象越严重，拖尾就越长，这是由于粒子运移通道被阻塞且阻碍了气体的运动。图2(b)是半对数坐标图，可以观察到不同时间阶数α下，浓度穿透曲线的尾部斜率不同，说明了时间阶数α是影响拖尾现象的主要参数。

图2 算例2：已知初始条件和边界条件，由时间分数阶对流—扩散模型数值求解获得的不同阶数α下的浓度穿透曲线(a)。模型参数如下：v=0.2，K=0.01，T=200和L=100。(b)是半对数坐标图。

2.3 算例3

在边界处给定一瞬时点源，利用时空分数阶对流—扩散模型(3)模拟粒子运移规律。初始值为C0，入口边界值为0，右边边界条件为自由边界。数值模拟结果如下。

从图3可以得出，空间分数阶对流—扩散模型能够较好地模拟粒子的早到达现象。随着空间阶数的减小，早到达现象越明显。结合实际情况分析，若粒子在双重多孔介质中运移时，介质中存在大尺度裂缝通道，促使粒子加速运动，因此出现了早期的超扩散现象。气体在双重多孔介质中的运移过程与孔隙结构和空间分布密切相关。简而言之，分数阶导数模型是描述具有记忆效应的气体运移过程的有效工具，下面将通过实验数据验证该模型。

图3 展示的是时空分数阶对流—扩散模型模拟的数值结果(a)，(b)图是半对数下的坐标图。模型参数为：C0=20，v=0.1，K=0.05，T=200和L=20。

3 结果讨论

基于文献[30]中的一系列实验研究，探讨不同条件下包括压力、温度、气体流速和介质不均匀性等，CO2驱替页岩上CH4的动力学理论机制和CO2驱替的有利条件。经典对流—扩散模型在模拟CO2在页岩储藏的裂隙和基质系统中的传质过程时，存在明显不足。因此借助分数阶导数模型研究气体运移规律，模拟气体的浓度穿透曲线。

在非均匀页岩储藏中，CO2置换CH4的过程较为复杂，CO2的运动会出现反常扩散过程，其本质是时间上有记忆性和空间上有非局域性的过程，而整数阶导数极限定义具有局域性不能准确地描述这类反常扩散过程。从图4可得出，分数阶导数模型相比于传统模型更具准确性和实用性。分数阶对流—扩散模型能够较好地描述驱替页岩过程中CO2在非均匀储层中的次扩散过程。分数阶导数模型因其具有记忆和遗传、路径依赖性，能够较准确地刻画气体粒子的反常扩散过程。

图4 利用时间分数阶对流—扩散方程(FCDE)模型和经典对流—扩散方程(CDE)模型模拟CO2的浓度穿透曲线。模型参数分别为：L=50 cm，K=10-4 cm2/s，v=0.15 cm/s。

4 总结

分数阶导数是微分—积分形式表示的非局域性算子，能准确地刻画复杂介质结构或流场中粒子运动的全域相关和历史依赖性，已成为一种流行的建模工具。

分数阶对流—扩散方程模型在描述地下含水层、河流、土壤和裂隙岩体等复杂介质中的流体扩散或迁移过程时，具有参数少、物理意义清楚和描述准确等优点。

从数值模拟结果得出，分数阶导数模型能够较准确地刻画粒子的反常扩散现象，通过改变分数阶模型中的时间或空间阶数，则可反映不同程度的次扩散或超扩散现象。

从实验拟合结果得出，整数阶导数的局部极限定义不适合描述地质储存中气体粒子运动过程具有的路径和时间依赖特征。而分数阶导数模型能够较好地拟合实验数据，准确刻画了气体的运移过程。