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

权利要求书

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

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；

页岩固体质点的速度为：

式中，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—页岩固体质点的位移量，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—体积应变，无量纲；

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

式中，Δε—基质收缩变形率，无量纲；Δ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为时间步数。