CN109284571A - 一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法 - Google Patents

Abstract

本发明公开了一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法，它包括以下步骤：S1、建立多尺度页岩孔缝中气体流动形态分类模式；S2、建立二氧化碳置换页岩气多尺度多场耦合渗流数学模型，具体包括多尺度页岩多组分气体多场耦合流动方程、页岩骨架变形方程、页岩储层能量守恒方程、应力场和孔缝介质变形场作用下的物性参数方程、定解条件；S3、二氧化碳置换页岩气多尺度多场耦合渗流数学模型的求解应用。本发明的有益效果是：本发明在建立的多尺度页岩孔缝中气体流动形态分类模式的基础上，建立了页岩基质纳米孔、基质纳微米孔及裂缝中二氧化碳置换页岩气多组分气体体系多场耦合渗流数学模型，并建立块中心网格块内及网格块间的差分数值模型。

Description

一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法

技术领域

本发明涉及一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法。

背景技术

页岩气在全球范围内具有非常丰富的储量，据EIA(美国能源信息情报署)2015年统计：世界上共有42个国家发现了95个含页岩气盆地、137套页岩地层，页岩气资源量总计约1013×10¹²m³，页岩气的技术可采储量达到220.69×10¹²m³。据2015年国土资源部发布的《中国页岩气资源调查报告(2014)》评价结果显示，我国页岩气资源储量为134×10¹²m³，可采储量为25×10¹²m³。页岩气同时又是一种具有典型的多尺度孔缝，致密的储渗条件，不经压裂改造就不具有工业产能的非常规天然气。

二氧化碳是一种破坏环境的温室气体，减排二氧化碳是当今国际上面临的共同难题。研究表明，注入页岩层的气态(超临界态)二氧化碳具有强吸附、易扩散、增渗增能等优势，利用注入二氧化碳的方式能达到强化页岩气开采同时满足二氧化碳地质埋存的双重效果。因此，二氧化碳置换页岩气技术正在成为当前页岩气开采暨二氧化碳埋存领域的热点之一，二氧化碳置换页岩气渗流模型的建立是其中必须攻克的难点及关键技术之一。

页岩气与煤层气具有相似的渗流过程，有关煤层气地下渗流的数学模型已经发展得较为完善，现有的页岩气地下渗流的数学模型也大多取自对煤层气渗流模型的修改完善。人们也提出过二氧化碳强化煤层气开采技术，同时建立了二氧化碳强化煤层气开采渗流数学模型。

如前所述，二氧化碳置换页岩气技术是当前页岩气开采暨二氧化碳埋存领域的热点之一，并没有发展成成熟技术，也没有完整的表征储层条件下二氧化碳置换页岩气渗流的数学模型。

限于煤岩与页岩的巨大差别，煤岩与页岩的孔缝系统也差异巨大，因此，有关二氧化碳强化煤层气开采渗流数学模型并不能直接用于二氧化碳置换页岩气渗流模拟。

储层条件下二氧化碳置换页岩气，二氧化碳与页岩混合多组分气体的渗流过程必将受到温度场、压力场、应力场、孔缝介质变形场、气体解吸扩散渗流场五场耦合作用的影响，同时，页岩特有的纳米孔隙系统更增加了二氧化碳置换页岩气地下渗流的复杂性。因此，必须解析页岩多尺度孔缝系统，建立每个尺度孔缝系统的流态模式，然后据此建立相应的数学模型，在数学建模过程中，需要考虑温度场、压力场、应力场、孔缝介质变形场、气体解吸扩散渗流场五场耦合作用的影响。

发明内容

本发明的目的在于克服现有技术的缺点，提供一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法。

本发明的目的通过以下技术方案来实现：一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法，它包括以下步骤：

S1、建立多尺度页岩孔缝中气体流动形态分类模式，具体包括以下步骤：

S11、建立克努森数与页岩孔径尺寸间的一一对应关系，解析得到页岩中分别发育的基质纳米孔、基质纳微米孔及人工或天然裂缝三种不同尺度的孔缝中气体的各种流动形态；

S12、将纳米孔、纳微米孔和裂缝中发生的步骤S11的流动形态及其相互作用关系用框图的形式表达出来，构建形成页岩多尺度孔缝中气体流动形态的分类模式；

S2、建立二氧化碳置换页岩气多尺度多场耦合渗流数学模型，具体包括多尺度页岩多组分气体多场耦合流动方程、页岩骨架变形方程、页岩储层能量守恒方程、应力场和孔缝介质变形场作用下的物性参数方程、定解条件。

S21、建立多尺度页岩多组分气体多场耦合流动方程，依据S1获得的多尺度页岩孔缝中气体流动形态分类模式具体建立包括纳米孔中的流动方程、纳微米孔中的流动方程和裂缝中的流动方程；

S211、建立纳米孔中的流动方程：

二氧化碳置换页岩气将产生大量的气体吸附于页岩基质孔隙壁面，气体组分i在页岩基质表面的吸附量与压力、温度之间的关系式为：

Knudsen扩散的质量流量可以表示为：

基质纳米孔中多组分气体在温度场、压力场、气体解吸扩散场耦合作用下的流动方程为：

式中，q_ads,i—单位体积的i组分气体的吸附量，kg/m³；ρ_s—标准状态下混合气体的密度，g/cm³；M_i—气体组分i的摩尔质量，kg/mol；V_std-标准状态下的摩尔体积，m³/mol；q_std,i-标准状态下页岩单位质量的i组分吸附体积，m³/kg；D-吸附特征常数；m-吸附特征常数；R-普适气体常数，8.314J/(mol*K)；p-体系的压力，MPa；T-体系的温度，K；p_ci-i组分的临界压力，MPa；T_ci-i组分的临界温度，K；C-气体的摩尔浓度，mol/m³；J_i-气体组分i的质量流量，kg/m²/s；M_i-气体组分i的摩尔质量，kg/mol；D_Ki-气体组分i的Knudsen扩散系数，无量纲；x_i-气体组分i在气相中的摩尔分数；Z_i-气体组分i的偏差系数，无量纲；Z-混合气体系偏差系数，无量纲；ρ_i为气体组分i的密度，kg/m³；

S212、建立纳微米孔中的流动方程：

基质纳微米孔中多组分气体在温度场、压力场、气体解吸扩散渗流场耦合作用下的流动方程为：

式中，Q_i为气体组分i的基岩-裂缝窜流量，kg/s；k—表观渗透率，μm²；k_n—克努森数，无量纲；k_∞—绝对渗透率，μm²；α—稀疏系数，无量纲；Z_i—气体组分i的偏差系数，无量纲；σ—形状因子，1/m²；k_m—基质渗透率，μm²；μ_i—气体组分i的粘度，Pa·s；p_m—基质压力，MPa；p_f—裂缝压力，MPa；φ_m2—纳微米孔孔隙度，无量纲；b为滑移系数，无量纲；

S213、建立裂缝中的流动方程，包括多场耦合运动方程的建立和多场耦合连续性方程的建立：所述多场耦合运动方程的建立：

流体渗流时的真实速度方程为：

式中，V_a—流体相a的真实运动速度，m/s；V_ra—流体相a相对于页岩固体质点的流动速度，m/s；V_s—页岩固体质点的真实运动速度，m/s；v_a—流体相a渗流的达西速度，m/s；Q—流体通过渗流断面的体积流量，m³/s；f—流体渗流断面上，各个孔隙通道截面积之和，m²；S_a—流体相a的饱和度；φ—岩块的孔隙度，无量纲；

气体的真实速度、达西速度以及页岩固体质点速度三者之间的关系式为：

气体的达西速度为：

式中，k—页岩气层的绝对渗透率，μm²；k_ra—流体相a的相对渗透率，无量纲；μ_a—流体相a的粘度，mPa·s；p_a—流体相a的压力，MPa；ρ_a—地层条件下流体相a的密度，kg/m³；g—重力加速度，m/s²；

页岩固体质点的速度为：

式中，U—页岩固体质点的位移矢量；

其位移分量为：

式中，V_sx,V_sy,V_sz—页岩固体质点在x，y，z方向上的速度，m/s；U_x,U_y,U_z—页岩固体质点在x，y，z方向上的位移量，m；

多场耦合连续性方程的建立：

温度场、压力场、应力场、孔缝介质变形场、气体渗流场耦合作用下二氧化碳置换页岩气体系的连续性方程为：

式中，φ_f-裂缝的孔隙度，无量纲；k_f-裂缝的渗透率，μm²；N_i-气体组分i从裂缝采出的质量流量，kg/s；Q_i为气体组分i的基岩-裂缝窜流量，kg/s；

S22、建立页岩骨架变形方程，具体包括建立页岩骨架的连续性方程、各向同性弹塑性岩石骨架的平衡方程、岩石骨架的几何方程、岩石骨架本构关系；

所述页岩骨架的连续性方程为：

所述各向同性弹塑性岩石骨架的平衡方程为：

σ′_ij＝λεδ_ij+2Gε_ij (2.2-2)

所述岩石骨架的几何方程为：

式中，ε_v-体积应变量；ε—各向同性弹塑性岩石体积应变量；σ_i′_j—有效应力；δ_ij—Kroneker张量；ε_ij—柯西应变张量或小应变张量；u_i,j,u_j,i—位移分量；λ—拉梅常数；G—岩石的剪切模量；

所述岩石骨架本构关系为：

式中，U—页岩固体质点的位移量，m；α—Biot参数；p-孔隙流体压力，MPa；F_i-体力分量，MPa；

S23、建立页岩储层能量守恒方程：

式中：ρ_s为页岩骨架密度；α_s为热膨胀系数；c_s为页岩骨架的热容；k_s为页岩骨架的热传导系数；k_g为气体的热传导系数；k_t＝(φ_m+φ_f)k_g+k_s；ρ_g为气体密度；c_g为气体的热容；

径向气体流量；Q_T为总的热源强度，Q_T＝(φ_m+φ_f)Q_Tg+Q_Ts；Q_Ts为页岩骨架热源强度；

Q_Tg为气体热源强度；q_g为气体的流量；c_p表示定压比热；ΔT为温度变化量；公式(2.3-1)左边第一项为系统内能的变化率；第二项为孔隙度变化使得气体质量变化从而使气体内能发生改变；第三项和第四项分别是由骨架变形产生的变形能以及热对流引起的附加项；第五项是导热项；

S24、建立应力场和孔缝介质变形场作用下的物性参数方程，具体包括建立的孔隙度方程、渗透率方程、孔隙压缩系数方程；

所述孔隙度方程为：

所述渗透率方程为：

式中，φ—变形后的孔隙度(对应于前述的基质孔隙度φ_m)，无量纲；φ_o—原始条件下的孔隙度，无量纲；K—变形后的渗透率(对应于前述的基质及裂缝渗透率k_m，k_f)，μm²；K_o—原始条件下的渗透率，μm²；ε_v—体积应变，无量纲。

所述孔隙压缩系数方程为：

式中，C_o—孔隙压缩系数，无量纲；Δε-基质收缩变形率，无量纲；ΔP—地层压力变化量，MPa；Ln取自然对数；

S25、定解条件的设定，具体包括设定页岩骨架变形初始及边界条件、多组分渗流初始及边界条件、储层温度初始条件；

所述页岩骨架变形初始及边界条件为：

第一类页岩骨架的表面力已知：

σ_ijL_j＝S_i，j(x,y,z) (2.5-1)

式中，L_j-边界的方向导数；S_i,j(x,y,z)-表面力分别在x,y,z三轴方向上的分布函数；第二类页岩骨架在x,y,z三轴方向上的表面位移已知：

u_i,j＝g_i,j(x,y,z) (2.5-2)

所述多组分渗流初始及边界条件为：

初始条件(t₀时间的压力P)：

式中，p_i-原始地层压力，MPa；

第一类定压边界条件(P_c)：

P_c＝f_p(x,y,t) (2.5-4)

式中，f_p(x,y,t)-边界上在t时间x,y方向上的已知函数；

第二类定流量边界条件

式中，f_q(x,y,t)-边界上在t时间x,y方向上的已知函数；

所述储层温度初始条件(t₀时间的温度T)为：

式中，f_k(z)-原始地层深度z处的已知温度函数；

S3、二氧化碳置换页岩气多尺度多场耦合渗流数学模型的求解应用，具体包括网格块内流动方程的差分离散、网格块间流动方程的差分离散；

S31、网格块内流动方程的差分离散

为了简化模型，网格块内温度视为恒定；在网格块内，考虑解吸、扩散和粘性流，基质纳米孔中的气体解吸后直接扩散到基质微纳米孔中，基质微纳米孔中的气体直接流到与之临近的裂缝中，且为拟稳态流动，同一基质块内压力处处相等，而网格块内裂缝的压力也是稳定的。纳微米孔中连续性方程中的气源为纳米孔中气体的解吸扩散，该阶段的扩散可采用拟稳态扩散模型进行描述(2.1-3)。拟稳态扩散忽略空间上气体的浓度扩散，而是以时间变化为核心，每个时间点的气体都能达到平均浓度，而不同时间段内的平均浓度各不相同。由此可以得到储层压力从降低到时，在纳米孔内气体的扩散量可以离散为：

在网格中心点(x_i,y_j,z_k)上对多组分气体在纳微米孔中流动的数学模型(2.1-4)不需要在空间上进行差分离散，仅对公式(2.1-4)右端项进行时间差分离散。

S32、网格块间流动方程的差分离散

流体在裂缝中渗流时，考虑温度和岩石变形对渗流的影响，在网格点(x_i,y_j,z_k)上对多组分气体在页岩裂缝多场耦合连续性流动方程(2.1-12)分别对压力p_f和温度T进行时间和空间上的差分离散，整理后最终得到：

其中：

式中：本发明采用直角坐标下的块中心差分网格，对于结点坐标(i，j，k)，则在x、y、z三个方向上，其前、后、上、下、左、右、相邻块的中心坐标分别为(i-1，j，k)、(i+1，j，k)、(i，j，k-1)、(i，j，k+1)、(i，j-1，k)、(i，j+1，k)，这个块的前后上下左右边界的坐标相应的分别为(i-1/2，j，k)、(i+1/2，j，k)、(i，j，k-1/2)、(i，j，k+1/2)、(i，j-1/2，k)、(i，j+1/2，k)；网格块的中心满足：x_i＝(x_i+1/2+x_i-1/2)/2，Δx_i＝x_i+1/2-x_i-1/2；Δx_i,Δy_j,Δz_k分别为x、y、z三个方向上的网格步长；Δt为时间步长；上标n,n+1为时间步数；别的参数符号与前面相同。

本发明具有以下优点：本发明在建立的多尺度页岩孔缝中气体流动形态分类模式的基础上，建立页岩基质纳米孔、基质纳微米孔及裂缝中二氧化碳置换页岩气多组分气体体系，在温度场、压力场、应力场、孔缝介质变形场、气体解吸扩散渗流场这五场耦合作用下的包括多组分气体流动方程、页岩骨架变形方程、页岩储层能量守恒方程、辅助方程和物性参数耦合方程，以及相应的定解条件的完整的一套数学模型，并建立了块中心网格块内及网格块间的差分数值模型。

附图说明

图1为本发明的建模流程图；

图2为多尺度页岩孔缝中气体流动形态分类模式图；

图3为渗流场、应力场以及温度场之间的耦合关系。

具体实施方式

下面结合附图对本发明做进一步的描述，本发明的保护范围不局限于以下所述：

如图1所示，一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法，它包括以下步骤：

S1、建立多尺度页岩孔缝中气体流动形态分类模式，具体包括以下步骤：

S11、建立克努森数(Kn)与页岩孔径尺寸间的一一对应关系如表1所示，得到页岩中主要发育的孔径＜10nm的基质纳米孔、10nm＜孔径＜1μm的基质纳微米孔(包括部分微裂缝)，及缝宽＞1μm的人工或天然裂缝，刚好分别对应于Kn>0.1、0.001<Kn<0.1和Kn<0.001。

表1页岩三尺度孔缝、克努森数(Kn)及气体流态划分方案

由此将基质纳米孔(简称“纳米孔”，下同)、基质纳微米孔(简称“纳微米孔”，下同)及人工或天然裂缝(简称“裂缝”，下同)作为页岩发育的三种不同尺度的孔缝，利用建立的孔缝尺度与Kn数间的关系，即可解析获得每个尺度的页岩孔缝中气体的各种流动形态，从而形成了页岩储层多孔介质中气体流态划分方案(表1)。

不同的Kn数对应着气体的不同流动形态，按照解析得到的多尺度页岩孔缝与Kn数间的关系，不难发现，气体在页岩的纳米孔、纳微米孔和裂缝三种不同尺度的孔缝空间中的流动具有阶段性，从而将气体在多尺度页岩孔缝中的流动划分为纳米孔、纳微米孔和裂缝三个阶段。这三个阶段气体的运移为串联运移，具体的流动形态阐述如下：

(1)气体在纳米孔中主要发生吸附解吸及少量Knudsen扩散；(2)从纳米孔中扩散到纳微米孔中的气体将与纳微米孔中原有的气体一起在纳微米孔中发生Knudsen扩散与粘性流；(3)进入裂缝中的气体将主要发生达西渗流；

S12、将纳米孔、纳微米孔和裂缝中发生的步骤S11的流动形态及其相互作用关系用框图的形式表达出来，构建形成页岩多尺度孔缝中气体流动形态的分类模式，如图2所示，该模式解析了页岩多尺度孔缝——基质纳米孔、基质纳微米孔及裂缝中的渗流流态，并建立了各渗流流态间的传质、传能关系特征模型；

S2、建立二氧化碳置换页岩气多尺度多场耦合渗流数学模型，具体包括多尺度页岩多组分气体多场耦合流动方程、页岩骨架变形方程、页岩储层能量守恒方程、应力场和孔缝介质变形场作用下的物性参数方程、定解条件。

S21、建立多尺度页岩多组分气体多场耦合流动方程，依据S1获得的多尺度页岩孔缝中气体流动形态分类模式具体建立包括纳米孔中的流动方程、纳微米孔中的流动方程和裂缝中的流动方程；

S211、建立纳米孔中的流动方程：

二氧化碳置换页岩气将产生大量的气体吸附于页岩基质孔隙壁面，气体组分i在页岩基质表面的吸附量与压力、温度之间的关系式为：

Knudsen扩散的质量流量可以表示为：

基质纳米孔中多组分气体在温度场、压力场、气体解吸扩散场耦合作用下的流动方程为：

其中：ρ_i为气体组分i的密度，kg/m³，可表示为式中，q_ads,i-单位体积的i组分气体的吸附量，kg/m³；ρ_s-标准状态下混合气体的密度，g/cm³；M_i—气体组分i的摩尔质量，kg/mol；V_std—标准状态下的摩尔体积，m³/mol；q_std,i—标准状态下页岩单位质量的i组分吸附体积，m³/kg；D—吸附特征常数；m—吸附特征常数；R—普适气体常数，8.314J/(mol*K)；p—体系的压力，MPa；T—体系的温度，K；p_ci—i组分的临界压力，MPa；T_ci—i组分的临界温度，K；C—气体的摩尔浓度，mol/m³；J_i—气体组分i的质量流量，kg/m²/s；M_i—气体组分i的摩尔质量，kg/mol；D_Ki—气体组分i的Knudsen扩散系数，无量纲；x_i—气体组分i在气相中的摩尔分数；Z_i—气体组分i的偏差系数，无量纲；Z—混合气体系偏差系数，无量纲；

S212、建立纳微米孔中的流动方程：

基质纳微米孔中多组分气体在温度场、压力场、气体解吸扩散渗流场耦合作用下的流动方程为：

其中Q_i为气体组分i的基岩-裂缝窜流量，kg/s，表示为

式中，k—表观渗透率，μm²；k_n—克努森数，无量纲；k_∞—绝对渗透率，μm²；α—稀疏系数，无量纲；Z_i—气体组分i的偏差系数，无量纲；σ—形状因子，1/m²；k_m—基质渗透率，μm²；μ_i—气体组分i的粘度，Pa·s；p_m—基质压力，MPa；p_f—裂缝压力，MPa；φ_m2—纳微米孔孔隙度，无量纲；b为滑移系数，无量纲；

S213、建立裂缝中的流动方程，包括多场耦合运动方程的建立和多场耦合连续性方程的建立：

所述多场耦合运动方程的建立：

采用达西定律来描述气体在页岩裂缝中的渗流，建立考虑岩石固体质点的位移，温度场、压力场、应力场、孔缝介质变形场、气体渗流场耦合作用下气体渗流的真实速度。伴随着气体的开采，压力不断改变将会使页岩有效应力发生改变，从而使页岩固体质点发生位移，页岩固体质点的位移将会对流体的流速产生影响，页岩固体质点速度V_s反映了有效应力的改变对页岩储层骨架变形的影响，进一步推导得流体渗流时的真实速度方程为：

式中，V_a—流体相a的真实运动速度，m/s；V_ra—流体相a相对于页岩固体质点的流动速度，m/s；V_s—页岩固体质点的真实运动速度，m/s；v_a—流体相a渗流的达西速度，m/s；Q—流体通过渗流断面的体积流量，m³/s；f—流体渗流断面上，各个孔隙通道截面积之和，m²；S_a—流体相a的饱和度；φ—岩块的孔隙度，无量纲；

气体的真实速度、达西速度以及页岩固体质点速度三者之间的关系式为：

气体的达西速度为：

式中，k—页岩气层的绝对渗透率，μm²；k_ra—流体相a的相对渗透率，无量纲；μ_a—流体相a的粘度，mPa·s；p_a—流体相a的压力，MPa；ρ_a—地层条件下流体相a的密度，kg/m³；g—重力加速度，m/s²；

页岩固体质点的速度为：

式中，U—页岩固体质点的位移矢量；

其位移分量为：

式中，V_sx,V_sy,V_sz—页岩固体质点在x，y，z方向上的速度，m/s；U_x,U_y,U_z—页岩固体质点在x，y，z方向上的位移量，m；

所述多场耦合连续性方程的建立：

气体以游离态的形式存在于页岩裂缝中，且不考虑气体在裂缝中的吸附与解吸，该阶段气体的流动符合达西定律，不考虑毛管压力和重力的影响，且只存在二氧化碳和页岩气气体体系。根据质量守恒原理，即可得到组分i的连续性方程

结合公式(2.1-4)、(2.1-8)、(2.1-9)和(2.1-14)即可推导整理得到温度场、压力场、应力场、孔缝介质变形场、气体渗流场耦合作用下二氧化碳置换页岩气体系的连续性方程为：

可以看出，对于所建立的多组分气体多场耦合渗流数学模型，比常规的渗流数学模型多了一项耦合项该项表示固体质点位移对气体渗流的影响，这也体现了多场作用耦合关系；

式中，φ_f—裂缝的孔隙度，无量纲；k_f—裂缝的渗透率，μm²；N_i—气体组分i从裂缝采出的质量流量，kg/s；

S22、建立页岩骨架变形方程，具体包括建立页岩骨架的连续性方程、各向同性弹塑性岩石骨架的平衡方程、岩石骨架的几何方程、岩石骨架本构关系；

所述页岩骨架的连续性方程为：

所述各向同性弹塑性岩石骨架的平衡方程为：

σ′_ij＝λεδ_ij+2Gε_ij (2.2-2)

所述岩石骨架的几何方程为：

式中，ε_v—体积应变量；ε—各向同性弹塑性岩石体积应变量；σ′_ij—有效应力；δ_ij—Kroneker张量；ε_ij—柯西应变张量或小应变张量；u_i,j,u_j,i—位移分量；λ—拉梅常数；G—岩石的剪切模量；

所述岩石骨架本构关系为：

式中，U—页岩固体质点的位移量，m；α—Biot参数；p—孔隙流体压力，MPa；F_i—体力分量，MPa；

S23、建立页岩储层能量守恒方程：

式中：ρ_s为页岩骨架密度；α_s为热膨胀系数；c_s为页岩骨架的热容；k_s为页岩骨架的热传导系数；k_g为气体的热传导系数；k_t＝(φ_m+φ_f)k_g+k_s；ρ_g为气体密度；c_g为气体的热容；

径向气体流量；Q_T为总的热源强度，Q_T＝(φ_m+φ_f)Q_Tg+Q_Ts；Q_Ts为页岩骨架热源强度；

Q_Tg为气体热源强度；q_g为气体的流量；c_p表示定压比热；ΔT为温度变化量；公式(2.3-1)左边第一项为系统内能的变化率；第二项为孔隙度变化使得气体质量变化从而使气体内能发生改变；第三项和第四项分别是由骨架变形产生的变形能以及热对流引起的附加项；第五项是导热项；

S24、建立应力场和孔缝介质变形场作用下的物性参数方程，页岩储层渗流场、应力场以及温度场之间的耦合关系如图3所示，据此建立了在应力场和孔缝介质变形场作用下的孔隙度、渗透率和岩石压缩系数方程，具体包括建立的孔隙度方程、渗透率方程、孔隙压缩系数方程；

所述孔隙度方程为：

所述渗透率方程为：

式中，φ—变形后的孔隙度(对应于前述的基质孔隙度φ_m)，无量纲；φ_o—原始条件下的孔隙度，无量纲；K—变形后的渗透率(对应于前述的基质及裂缝渗透率k_m，k_f)，μm²；K_o—原始条件下的渗透率，μm²；ε_v—体积应变，无量纲。

所述孔隙压缩系数方程为：

式中，C_o—孔隙压缩系数，无量纲；Δε—基质收缩变形率，无量纲；ΔP—地层压力变化量，MPa；

S25、定解条件的设定，具体包括设定页岩骨架变形初始及边界条件、多组分渗流初始及边界条件、储层温度初始条件；

所述页岩骨架变形初始及边界条件为：

第一类页岩骨架的表面力已知：

σ_ijL_j＝S_i，j(x,y,z) (2.5-1)

式中，L_j—边界的方向导数；S_i,j(x,y,z)—表面力分别在x,y,z三轴方向上的分布函数；第二类页岩骨架在x,y,z三轴方向上的表面位移已知：

u_i,j＝g_i,j(x,y,z) (2.5-2)

所述多组分渗流初始及边界条件为：

初始条件(t₀时间的压力P)：

式中，p_i—原始地层压力，MPa；

第一类定压(P_c)边界条件：

P_c＝f_p(x,y,t) (2.5-4)

式中，f_p(x,y,t)—边界上在t时间x,y方向上的已知函数；

第二类定流量边界条件

式中，f_q(x,y,t)—边界上在t时间x,y方向上的已知函数；

所述储层温度初始条件(t₀时间的温度T)为：

式中，f_k(z)—原始地层深度z处的已知温度函数；

S3、二氧化碳置换页岩气多尺度多场耦合渗流数学模型的求解应用，具体包括网格块内流动方程的差分离散、网格块间流动方程的差分离散；

S31、网格块内流动方程的差分离散

为了简化模型，网格块内温度视为恒定；在网格块内，考虑解吸、扩散和粘性流，基质纳米孔中的气体解吸后直接扩散到基质微纳米孔中，基质微纳米孔中的气体直接流到与之临近的裂缝中，且为拟稳态流动，同一基质块内压力处处相等，而网格块内裂缝的压力也是稳定的。纳微米孔中连续性方程中的气源为纳米孔中气体的解吸扩散，该阶段的扩散可采用拟稳态扩散模型进行描述(2.1-3)。拟稳态扩散忽略空间上气体的浓度扩散，而是以时间变化为核心，每个时间点的气体都能达到平均浓度，而不同时间段内的平均浓度各不相同。由此可以得到储层压力从降低到时，在纳米孔内气体的扩散量可以离散为：

在网格中心点(x_i,y_j,z_k)上对多组分气体在纳微米孔中流动的数学模型(2.1-5)不需要在空间上进行差分离散，仅对公式(2.1-5)右端项进行时间差分离散。

S32、网格块间流动方程的差分离散

流体在裂缝中渗流时，考虑温度和岩石变形对渗流的影响，在网格点(x_i,y_j,z_k)上对多组分气体在页岩裂缝多场耦合连续性流动方程(2.1-15)分别对压力p_f和温度T进行时间和空间上的差分离散，整理后最终得到：

其中：

式中：本发明采用直角坐标下的块中心差分网格，对于结点坐标(i，j，k)，则在x、y、z三个方向上，其前、后、上、下、左、右、相邻块的中心坐标分别为(i-1，j，k)、(i+1，j，k)、(i，j，k-1)、(i，j，k+1)、(i，j-1，k)、(i，j+1，k)，这个块的前后上下左右边界的坐标相应的分别为(i-1/2，j，k)、(i+1/2，j，k)、(i，j，k-1/2)、(i，j，k+1/2)、(i，j-1/2，k)、(i，j+1/2，k)；网格块的中心满足：x_i＝(x_i+1/2+x_i-1/2)/2，Δx_i＝x_i+1/2-x_i-1/2；Δx_i,Δy_j,Δz_k分别为x、y、z三个方向上的网格步长；Δt为时间步长；上标n,n+1为时间步数；别的参数符号与前面相同。

因此本发明在建立的多尺度页岩孔缝中气体流动形态分类模式的基础上，建立页岩基质纳米孔、基质纳微米孔及裂缝中二氧化碳置换页岩气多组分气体体系，在温度场、压力场、应力场、孔缝介质变形场、气体解吸扩散渗流场这五场耦合作用下的包括多组分气体流动方程、页岩骨架变形方程、页岩储层能量守恒方程、辅助方程和物性参数耦合方程，以及相应的定解条件的完整的一套数学模型，并建立了块中心网格块内及网格块间的差分数值模型。

Claims (1)

1.一种二氧化碳置换页岩气多尺度多场耦合渗流数学建模方法，其特征在于：它包括以下步骤：

S1、建立多尺度页岩孔缝中气体流动形态分类模式，具体包括以下步骤：

S11、建立克努森数与页岩孔径尺寸间的一一对应关系，解析得到页岩中分别发育的纳米孔、纳微米孔及裂缝三种不同尺度孔缝中气体的各种流动形态；

S12、将纳米孔、纳微米孔和裂缝中发生的步骤S11的流动形态及其相互作用关系用框图的形式表达出来，构建形成页岩多尺度孔缝中气体流动形态的分类模式；

S2、建立二氧化碳置换页岩气多尺度多场耦合渗流数学模型，具体包括多尺度页岩多组分气体多场耦合流动方程、页岩储层能量守恒方程、页岩骨架变形方程、应力场和孔缝介质变形场作用下的物性参数方程、定解条件。

S21、建立多尺度页岩多组分气体多场耦合流动方程，依据S1获得的多尺度页岩孔缝中气体流动形态分类模式具体建立包括纳米孔中的流动方程、纳微米孔中的流动方程和裂缝中的流动方程；

S211、纳米孔中的流动方程的建立：

二氧化碳置换页岩气将产生大量的气体吸附于页岩基质孔隙壁面，气体组分i在页岩基质表面的吸附量与压力、温度之间的关系式为：

Knudsen扩散的质量流量可以表示为：

J_i＝-M_iD_Ki▽x_iC (2.1-2)

建立得到基质纳米孔中多组分气体在温度场、压力场、气体解吸扩散场耦合作用下的流动方程为：

式中，q_ads,i—单位体积的i组分气体的吸附量，kg/m³；ρ_s—标准状态下混合气体的密度，g/cm³；M_i—气体组分i的摩尔质量，kg/mol；V_std—标准状态下的摩尔体积，m³/mol；q_std,i—标准状态下页岩单位质量的i组分吸附体积，m³/kg；D—吸附特征常数；m—吸附特征常数；R—普适气体常数，8.314J/(mol*K)；p-体系的压力，MPa；T-体系的温度，K；p_ci-i组分的临界压力，MPa；T_ci-i组分的临界温度，K；C-气体的摩尔浓度，mol/m³；J_i-气体组分i的质量流量，kg/m²/s；M_i-气体组分i的摩尔质量，kg/mol；D_Ki-气体组分i的Knudsen扩散系数，无量纲；x_i-气体组分i在气相中的摩尔分数；Z_i-气体组分i的偏差系数，无量纲；Z-混合气体系偏差系数，无量纲；ρ_i为气体组分i的密度，kg/m³；

S212、纳微米孔中的流动方程的建立：

基质纳微米孔中多组分气体在温度场、压力场、气体解吸扩散渗流场耦合作用下的流动方程为：

式中，Q_i为气体组分i的基岩-裂缝窜流量，kg/s；k—表观渗透率，μm²；k_n—克努森数，无量纲；k_∞—绝对渗透率，μm²；α—稀疏系数，无量纲；Z_i—气体组分i的偏差系数，无量纲；σ—形状因子，1/m²；k_m—基质渗透率，μm²；μ_i—气体组分i的粘度，Pa·s；p_m—基质压力，MPa；p_f—裂缝压力，MPa；φ_m2—纳微米孔孔隙度，无量纲；b为滑移系数，无量纲；

S213、裂缝中的流动方程的建立，包括多场耦合运动方程的建立和多场耦合连续性方程的建立：

所述多场耦合运动方程的建立：

流体渗流时的真实速度方程为：

式中，V_a—流体相a的真实运动速度，m/s；V_ra—流体相a相对于页岩固体质点的流动速度，m/s；V_s—页岩固体质点的真实运动速度，m/s；v_a—流体相a渗流的达西速度，m/s；Q—流体通过渗流断面的体积流量，m³/s；f—流体渗流断面上，各个孔隙通道截面积之和，m²；S_a—流体相a的饱和度；φ—岩块的孔隙度，无量纲；

气体的真实速度、达西速度以及页岩固体质点速度三者之间的关系式为：

气体的达西速度为：

式中，k-页岩气层的绝对渗透率，μm²；k_ra—流体相a的相对渗透率，无量纲；μ_a—流体相a的粘度，mPa·s；p_a—流体相a的压力，MPa；ρ_a—地层条件下流体相a的密度，kg/m³；g—重力加速度，m/s²；

页岩固体质点的速度为：

式中，U-页岩固体质点的位移矢量；

其位移分量为：

式中，V_sx,V_sy,V_sz-页岩固体质点在x，y，z方向上的速度，m/s；U_x,U_y,U_z-页岩固体质点在x，y，z方向上的位移量，m；

多场耦合连续性方程的建立：

温度场、压力场、应力场、孔缝介质变形场、气体渗流场耦合作用下二氧化碳置换页岩气体系的连续性方程为：

式中，φ_f—裂缝的孔隙度，无量纲；k_f—裂缝的渗透率，μm²；N_i—气体组分i从裂缝采出的质量流量，kg/s；Q_i为气体组分i的基岩-裂缝窜流量，kg/s；

S22、建立页岩骨架变形方程，具体包括建立页岩骨架的连续性方程、各向同性弹塑性岩石骨架的平衡方程、岩石骨架的几何方程、岩石骨架本构关系；

所述页岩骨架的连续性方程为：

所述各向同性弹塑性岩石骨架的平衡方程为：

σ′_ij＝λεδ_ij+2Gε_ij (2.2-2)

所述岩石骨架的几何方程为：

式中，ε_v-体积应变量；ε-各向同性弹塑性岩石体积应变量；σ′_ij—有效应力；δ_ij—Kroneker张量；ε_ij—柯西应变张量或小应变张量；u_i,j,u_j,i—位移分量；λ—拉梅常数；G—岩石的剪切模量；

所述岩石骨架本构关系为：

ε＝▽U (2.2-5)

式中，U—页岩固体质点的位移量，m；α-Biot参数；p-孔隙流体压力，MPa；F_i—体力分量，MPa；

S23、建立页岩储层能量守恒方程：

式中：ρ_s为页岩骨架密度；α_s为热膨胀系数；c_s为页岩骨架的热容；k_s为页岩骨架的热传导系数；k_g为气体的热传导系数；k_t＝(φ_m+φ_f)k_g+k_s；ρ_g为气体密度；c_g为气体的热容；径向气体流量；Q_T为总的热源强度，Q_T＝(φ_m+φ_f)Q_Tg+Q_Ts；Q_Ts为页岩骨架热源强度；Q_Tg为气体热源强度；q_g为气体的流量；c_p表示定压比热；ΔT为温度变化量；公式(2.3-1)左边第一项为系统内能的变化率；第二项为孔隙度变化使得气体质量变化从而使气体内能发生改变；第三项和第四项分别是由骨架变形产生的变形能以及热对流引起的附加项；第五项是导热项；

S24、建立应力场和孔缝介质变形场作用下的物性参数方程，具体包括建立的孔隙度方程、渗透率方程、孔隙压缩系数方程；

所述孔隙度方程为：

所述渗透率方程为：

式中，φ—变形后的孔隙度(对应于前述的基质孔隙度φ_m)，无量纲；φ_o-原始条件下的孔隙度，无量纲；K-变形后的渗透率(对应于前述的基质及裂缝渗透率k_m，k_f)，μm²；K_o—原始条件下的渗透率，μm²；ε_v—体积应变，无量纲。

所述孔隙压缩系数方程为：

式中，C_o—孔隙压缩系数，无量纲；Δε—基质收缩变形率，无量纲；ΔP—地层压力变化量，MPa；Ln取自然对数；

S25、定解条件的设定，具体包括设定页岩骨架变形初始及边界条件、多组分渗流初始及边界条件、储层温度初始条件；

所述页岩骨架变形初始及边界条件为：

第一类页岩骨架的表面力已知：

σ_ijL_j＝S_i，j(x,y,z) (2.5-1)

式中，L_j—边界的方向导数；S_i,j(x,y,z)—表面力分别在x,y,z三轴方向上的分布函数；

第二类页岩骨架在x,y,z三轴方向上的表面位移已知：

u_i,j＝g_i,j(x,y,z) (2.5-2)

所述多组分渗流初始及边界条件为：

初始条件(t₀时间的压力P)：

式中，p_i—原始地层压力，MPa；

第一类定压边界条件(P_c)：

P_c＝f_p(x,y,t) (2.5-4)

式中，f_p(x,y,t)—边界上在t时间x,y方向上的已知函数；

第二类定流量边界条件

式中，f_q(x,y,t)—边界上在t时间x,y方向上的已知函数；

所述储层温度初始条件(t₀时间的温度T)为：

式中，f_k(z)—原始地层深度z处的已知温度函数；

S3、二氧化碳置换页岩气多尺度多场耦合渗流数学模型的求解应用，具体包括网格块内流动方程的差分离散、网格块间流动方程的差分离散，以及二氧化碳置换页岩气多尺度多场耦合渗流数值模拟器的研制；

S31、网格块内流动方程的差分离散

为了简化模型，网格块内温度视为恒定；在网格块内，考虑解吸、扩散和粘性流，基质纳米孔中的气体解吸后直接扩散到基质微纳米孔中，基质微纳米孔中的气体直接流到与之临近的裂缝中，且为拟稳态流动，同一基质块内压力处处相等，而网格块内裂缝的压力也是稳定的。纳微米孔中连续性方程中的气源为纳米孔中气体的解吸扩散，该阶段的扩散可采用拟稳态扩散模型进行描述(2.1-3)。拟稳态扩散忽略空间上气体的浓度扩散，而是以时间变化为核心，每个时间点的气体都能达到平均浓度，而不同时间段内的平均浓度各不相同。由此可以得到储层压力从降低到时，在纳米孔内气体的扩散量可以离散为：

在网格中心点(x_i,y_j,z_k)上对多组分气体在纳微米孔中流动的数学模型(2.1-4)不需要在空间上进行差分离散，仅对公式(2.1-4)右端项进行时间差分离散。

S32、网格块间流动方程的差分离散

流体在裂缝中渗流时，考虑温度和岩石变形对渗流的影响，在网格点(x_i,y_j,z_k)上对多组分气体在页岩裂缝多场耦合连续性流动方程(2.1-12)分别对压力p_f和温度T进行时间和空间上的差分离散，整理后最终得到：

其中：

式中：本发明采用直角坐标下的块中心差分网格，对于结点坐标(i，j，k)，则在x、y、z三个方向上，其前、后、上、下、左、右、相邻块的中心坐标分别为(i-1，j，k)、(i+1，j，k)、(i，j，k-1)、(i，j，k+1)、(i，j-1，k)、(i，j+1，k)，这个块的前后上下左右边界的坐标相应的分别为(i-1/2，j，k)、(i+1/2，j，k)、(i，j，k-1/2)、(i，j，k+1/2)、(i，j-1/2，k)、(i，j+1/2，k)；网格块的中心满足：x_i＝(x_i+1/2+x_i-1/2)/2，Δx_i＝x_i+1/2-x_i-1/2；Δx_i,Δy_j,Δz_k分别为x、y、z三个方向上的网格步长；Δt为时间步长；上标n,n+1为时间步数；别的参数符号与前面相同。