CN112307653B - 一种页岩气藏产能数值模拟方法 - Google Patents

Abstract

本发明公开了一种页岩气藏产能数值模拟方法，涉及地质勘探技术领域，包括建立非结构四面体网络模型；建立基质系统流动方程；建立天然裂缝骨架体积应变和大尺度裂缝骨架体积应变的求解方程组；建立页岩气非线性渗流模型，并确定页岩气非线性渗流模型的定解条件；获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式；获得所述目标区域的水平井产量分布。本发明公开的页岩气藏产能数值模拟方法，针对现有技术中储层产量模拟方法中忽略地层骨架体积应变的问题，提供一种页岩气藏产能数值模拟方法，在气藏产能模拟中引入骨架体积应变，提高模拟结果的准确性。

Description

一种页岩气藏产能数值模拟方法

技术领域

本发明涉及气藏开发技术领域，特别涉及一种页岩气藏产能数值模拟方法。

背景技术

随着我国国民经济的快速增长，对于清洁能源的需求不断增大，深层页岩气储层的勘探与开发成为今年来的热点与难点。现场实践表明：对于深层页岩储层，由于基质渗透率低、储层致密，采用多段压裂水平井技术，扩散泄气面积，成为深层页岩储层高效开发的重要手段。同时，由于储层埋深大，储层上覆压力大，在开发过程中随着地层压力的降低，储层所受到的有效应力远远大于浅层页岩储层。应力变化引起的骨架体积应变对生产效果的影响日益突出。而现有技术中却缺少深层页岩储层的产能模拟方法，更少有研究地层骨架体积应变对储层产量的影响。

发明内容

本申请针对现有技术中储层产量模拟方法中忽略地层骨架体积应变的问题，提供一种页岩气藏产能数值模拟方法，在气藏产能模拟中引入骨架体积应变，提高模拟结果的准确性。

为了实现上述发明目的，本申请提供了以下技术方案：一种页岩气藏产能数值模拟方法，包括以下步骤：

根据目标区域内的地质模型数据体生成非结构四面体网络模型；

建立基质系统流动方程；并引入天然裂缝骨架体积应变和大尺度裂缝骨架体积应变，建立天然裂缝流动方程、大尺度裂缝流动方程；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程建立页岩流动方程组；

建立天然裂缝骨架体积应变和大尺度裂缝骨架体积应变的求解方程组；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程和求解方程组建立页岩气非线性渗流模型，并确定页岩气非线性渗流模型的定解条件；

根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式；

根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得所述目标区域的水平井产量分布。

在上述技术方案中，考虑了不同尺度介质流动规律以及储层骨架体积应变对流体流动的影响，在页岩气非线性渗流模型中引入了天然裂缝骨架体积应变和大尺度裂缝骨架体积应变；并采用有限体积方法，基于非结构四面体网格，获取深层多尺度页岩气藏流动过程数值解，进而获得所述目标区域的产能数值模拟结果；解决了传统结构网格模拟器中出现网格取向性强，无法表征复杂裂缝结构，忽略裂缝应力敏感性以及容易产生极小化网格的问题。

进一步地，所述基质系统流动方程为：

其中，

为基质孔隙度；q_m为页岩解吸气量；q_m-f为基质与天然裂缝窜流量；M_g为气体摩尔质量；p_m为基质压力；R为理想气体常数；T为温度；Z为真实气体偏差因子；α为稀薄气体系数；δ^m为页岩气基质孔隙阻塞系数；μ_g为气体黏度；τ^m为页岩气基质孔隙迂曲度；c_g为气体压缩系数；D_g为气体扩散系数；k_m为基质的有效渗透率；p_L为Langmuir压力；ρ_m为吸附相密度；V_m为吸附气量；V_std为标准状况下气体摩尔体积；V_L为Langmuir体积；p_f为天然裂缝压力。

进一步地，所述天然裂缝流动方程为：

其中，k_f为天然裂缝的有效渗透率、

为天然裂缝孔隙度；K_f为天然裂缝骨架弹性系数；β为Biot弹性系数；ε_f为天然裂缝骨架的体积应变。

进一步地，所述大尺度裂缝流动方程为：

其中，ε_hf为大尺度裂缝骨架的体积应变；q_well为水平井产量；k_hf为大尺度裂缝的有效渗透率；

为大尺度裂缝孔隙度；K_hf为大尺度裂缝骨架弹性系数；S_hf为大尺度裂缝综合弹性系数；d_hf为大尺度裂缝无因次开度因子；q_f-h为天然裂缝向大尺度裂缝窜流量；p_hf为大尺度裂缝压力。

进一步地，所述页岩流动方程组通过联立所述基质系统流动方程、天然裂缝流动方程和大尺度裂缝流动方程获得，其具体为：

进一步地，所述求解方程组根据补充广义胡克定律、质量守恒方程和应力平衡方程建立。

进一步地，所述求解方程组为：

其中，δ_ij为克罗内克函数(即，若i＝j，则δ_ij＝1；i≠j若，则δ_ij＝0)；ε_ij岩石骨架垂直于ij平面方向(i,j＝1,2,3)的体积应变；μ为气体黏度；σ_ij为岩石骨架垂直于ij平面方向(i,j＝1,2,3)的应力分量；σ_ij,j为应力张量的散度；F_i为岩石骨架所受体积力；u_i,j和u_j,i为应变引起的位移。

进一步地，根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，具体包括以下步骤：

S31)任意选取一四面体网络作为控制体单元，大尺度裂缝处理相邻两个控制提单元之间的表面；

S32)根据有限体积方法分别对基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程进行线性离散，分别获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式，建立代数方程；所述代数方程为：

其中，a_C、a_F分别为选定的控制体单元和与其相邻的控制体单元中的系数；T_C、T_F分别为选定的控制体单元和与其相邻的控制体单元中的未知数；b_C为余量项；C表示选取的控制体单元；指标F～NB(C)表示所有与控制体单元C相邻的控制体单元；

S33)循环重复步骤S32)，获得每一四面体网络的代数方程，并建立由a_C、a_F组成的全局系数矩阵和由T_C、T_F组成的未知数向量矩阵。

进一步地，根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得所述目标区域的水平井产量分布，具体包括以下步骤：

S41)预设当前时间步为k，获得在时间步为k+1时，所述基质系统流动方程和大尺度裂缝流动方程的半离散格式；

S42)采用Newton-Raphson迭代方法求解，获得

和

的解；并将求解得到的

和

的解带入所述天然裂缝流动方程的半离散格式进行显示求解，获得天然

的解；

S43)将步骤S41)和步骤S42)获得的

的带入所述求解方程组，更新应力状态，获得下一时间步的

的解；

S44)循环步骤S43)即可获得每一时间步的

的解，体积应力、水平井产量值，进而获得所述目标区域的水平井产量分布情况。

与现有技术相比，本发明的具有以下有益效果：

本申请公开了一种页岩气藏产能数值模拟方法，充分考虑了不同尺度介质流动规律以及储层骨架体积应变对流体流动的影响，在页岩气非线性渗流模型中引入了天然裂缝骨架体积应变和大尺度裂缝骨架体积应变；并采用有限体积方法，基于非结构四面体网格，获取深层多尺度页岩气藏流动过程数值解，进而获得所述目标区域的产能数值模拟结果；解决了传统结构网格模拟器中出现网格取向性强，无法表征复杂裂缝结构，忽略裂缝应力敏感性以及容易产生极小化网格的问题。

附图说明

图1为本发明公开的页岩气藏产能数值模拟方法进行求解的流程图示意图；

图2为本发明公开的页岩气藏产能数值模拟方法中建立的网格离散模型；

图3为本发明公开的页岩气藏产能数值模拟方法四面体网格单元示意图；

图4本发明公开的页岩气藏产能数值模拟方法控制体单元及其相邻控制体单元示意图；

图5本发明公开的页岩气藏产能数值模拟方法三维空间中三角形位置向量及表面向量计算方法；

图6本发明公开的页岩气藏产能数值模拟方法中网格编号储存示意图；

图7本发明公开的页岩气藏产能数值模拟方法的模拟结果与商业软件对比图；

图8本发明公开的页岩气藏产能数值模拟方法中骨架体积应变对页岩气产量的影响；

图9本发明公开的页岩气藏产能数值模拟方法中页岩骨架弹性参数对产量的影响。

具体实施方式

下面结合试验例及具体实施方式对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。

需要说明的是，虽然本申请的背景技术中记载了深层页岩储层的气藏产能数值模拟方法的相关现有技术，但本领域技术人员应该知道，本申请公开的深层页岩储层的气藏产能数值模拟方法也可以应用在其他地质构造的气藏，如采用水平井开发的致密砂岩气藏和致密碎屑岩气藏以及埋深较大的深盆气藏。

本申请中提出了一种页岩气藏产能数值模拟方法，参阅图1，该方法具体包括以下步骤：

S1：根据目标区域内的地质模型数据体生成非结构四面体网络模型；

S2：建立基质系统流动方程；并引入天然裂缝骨架体积应变和大尺度裂缝骨架体积应变，建立天然裂缝流动方程、大尺度裂缝流动方程；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程建立页岩流动方程组；

建立天然裂缝骨架体积应变和大尺度裂缝骨架体积应变的求解方程组；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程和求解方程组建立页岩气非线性渗流模型，并确定页岩气非线性渗流模型的定解条件；

S3：根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式；

S4：根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得所述目标区域的水平井产量分布。

需要说明的是，所述地质模型数据体包括所述目标区域的孔隙度、渗透率、储层厚度、天然裂缝分布图、吸附气量、上覆层压力、页岩杨氏模量、泊松比；所述孔隙度、渗透率、吸附气量可通过对现场取岩心样品并进行室内物理实验获得；储层厚度可通过微地震解释资料获得；天然裂缝分布图可通过研究层位地面露头测绘、成像测井解释以及微地震解释资料获得；上覆压力可通过储层埋深关系曲线获得；页岩杨氏模量、泊松比可通过室内岩石力学三轴压缩实验获得。需要说明的是，将上述所获的数据采用多点地质统计学方法进行地质建模生成地质模型数据场从而获得相应的地质模型数据体，其目的区域的边界即为所述地质模型数据体的边界。

需要说明的是，所述非结构四面体网络模型根据采用Matlab软件开源非结构网格剖分工具进行剖分得到。其具体为，根据Matlab软件开源非结构网格剖分工具包中的Distmesh的数据结构要求输入地质模型数据体所包括的数据，生成网络结构模型，并根据自由四面体算法对网络结构模型进行剖分，生成非结构四面体网络模型。

需要说明的是，所述基质系统流动方程为：

其中，

为基质孔隙度，通过孔隙度测试实验获得；q_m为页岩解吸气量，通过解吸附测试实验获得；q_m-f为基质与天然裂缝窜流量；M_g为气体摩尔质量；p_m为基质压力；R为理想气体常数；T为温度；Z为真实气体偏差因子，通过气体膨胀实验获得；α为稀薄气体系数，经验参数可查表获得；δ^m为页岩气基质孔隙阻塞系数；μ_g为气体黏度；τ^m为页岩气基质孔隙迂曲度，根据模型假设条件设定(δ_m/τ_m≤1)；μ_g为气体黏度；c_g为气体压缩系数，通过气体膨胀实验获得；D_g为气体扩散系数，通过扩散能力测试实验获得；k_m为基质的有效渗透率，通过渗透率测试实验获得；p_L为Langmuir压力；ρ_m为吸附相密度；V_m为吸附气量；V_L为Langmuir体积，上述解吸附参数均通过解吸附实验获得；V_std为标准状况下气体摩尔体积； p_f为天然裂缝压力。

在一些实施例中，所述基质系统流动方程通过以下步骤建立：

根据单相不稳定瞬态流动连续性方程，并根据页岩基质系统中吸附解析作用，获得页岩基质中流动过程满足的方程式：

其中，ρ_gm为基质内气体密度；v_m为基质表观流动速度向量；

引入真实气体状态方程、考虑滑脱效应和分子扩散的表观质量流量校正表达式，获得考虑了滑脱效应、分子扩散、吸附解析作用以及真实气体状态方程的基质系统流动方程，即为式(8)。

其中，所述真实气体状态方程为：

其中，ρ_gi(i＝m,f,hf)分别为基质、天然裂缝、大尺度裂缝中的气体密度；p_i(i＝m,f,hf) 分别为基质、天然裂缝和大尺度裂缝压力。

其中，所述表观质量流量校正表达式为：

其中，

为考虑滑脱效应的气体校正速度；

为考虑分子扩散的气体校正速度；其值分别为：

其中，α为稀薄气体系数；K_n为Knudsen数；

其中，页岩解吸气量q_m为：

其中，页岩解吸气量q_m根据Langmuir等温吸附方程获得。

其中，基质与天然裂缝窜流量q_m-f为：

式中符号定义：σ为裂缝形状因子。

需要说明的是，所述天然裂缝流动方程为：

其中，k_f为天然裂缝的有效渗透率，通过渗透率测试获得；

为天然裂缝孔隙度，通过孔隙度测试获得；K_f为天然裂缝骨架弹性系数，通过三轴岩石压缩测试获得；β为Biot弹性系数，用于描述骨架整体弹性变化能力，由经验公式获得；ε_f为天然裂缝骨架的体积应变，由三轴岩石压缩测试获得杨氏模量及泊松比后求解应力平衡方程获得。

在一些实施例中，所述天然裂缝流动方程通过以下步骤建立：

根据守恒原理获得页岩天然裂缝中流动方程的基本式：

其中，S_f为综合弹性系数；v_f为天然裂缝表观流动速度；Q_f为源汇项；

其中，所述综合弹性系数S_f通过以下方程获得；

所述天然裂缝表观流动速度v_f通过以下方程获得：

其中

其中，

为天然裂缝初始参考渗透率；λ为拉梅常数；

分别为天然裂缝的第一、二、三主应力。

所述汇源项Q_f通过以下方程获得：

其中：

其中，q_f-h为天然裂缝向大尺度裂缝窜流量；

为单位时间天然裂缝骨架变形引起的体积应变；β为Biot弹性系数；ε_f为天然裂缝骨架的体积应变。

需要说明的是，所述大尺度裂缝流动方程为：

其中，ε_hf为大尺度裂缝骨架的体积应变；q_well为水平井产量；k_hf为大尺度裂缝的有效渗透率；

为大尺度裂缝孔隙度；K_hf为大尺度裂缝骨架弹性系数；S_hf为大尺度裂缝综合弹性系数；d_hf为大尺度裂缝无因次开度因子；q_f-h为天然裂缝向大尺度裂缝窜流量；p_hf为大尺度裂缝压力。

在一些实施例中，所述大尺度裂缝流动方程通过以下步骤建立：

根据守恒原理获得页岩大尺度裂缝中流动方程的基本式：

S_hf为大尺度裂缝综合弹性系数；d_hf为大尺度裂缝无因次开度因子；v_hf为大尺度裂缝表观流动速度；Q_hf为源汇项；

需要说明的是，d_hf＝d/d_std，其中d为大尺度裂缝开度；d_std为标准单位长度，根据所研究的裂缝开度尺度确定(如所研究的大尺度裂缝开度为微米级，则d_std＝1μm)；d取值满足Weierstrass-Mandelbrot分形分布函数；所述Weierstrass-Mandelbrot分形分布函数为：

其中d(x,y)为粗糙裂缝表面轮廓的开度；D为粗糙度的分形维数，由开度测量结果进行盒维数拟合获得，且2<D<3；γ为与断面轮廓频谱密度有关的变量，一般取1.5,L为测量样品的特征尺度，M为曲面褶皱的重叠数；φ_m,n为随机相位，取值范围[0,2π]；n为频率序数，最低为0，最高为n_max＝int[log(L/L_S)/logγ]，L_S为样品最大长度，则有d_hf＝d(x,y)/d_std。

其中，大尺度裂缝表观流动速度v_hf为：

其中

式中符号定义：

为大尺度裂缝初始参考渗透率；λ为拉梅常数；

分别为大尺度裂缝的第一、二、三主应力。

源汇项

式中符号定义：q_well为水平井产量；

为单位时间大尺度裂缝骨架变形引起的体积应变；ε_f为大尺度裂缝骨架的体积应变。

因此，联立上述式(8)、(15)和(21)，即可获得所述页岩流动方程组；所述页岩流动方程组为：

由于页岩流动方程组中引入了天然裂缝体积应变和大尺度裂缝的体积应变，为了保证页岩流动方程组的封闭性，分别对天然裂缝流动方程和大尺度裂缝流动方程补充广义胡克定律、质量守恒方程以及应力平衡方程，来求解天然裂缝体积应变和大尺度裂缝体积应变。其中补充的广义胡克定律、质量守恒方程以及应力平衡方程即为求解方程组，其具体为：

其中，δ_ij为克罗内克函数(即，若i＝j，则δ_ij＝1；i≠j若，则δ_ij＝0)；ε_ij岩石骨架垂直于ij平面方向(i,j＝1,2,3)的体积应变；μ为气体黏度；σ_ij为岩石骨架垂直于ij平面方向(i,j＝1,2,3)的应力分量；σ_ij,j为应力张量的散度；F_i为岩石骨架所受体积力；u_i,j和u_j,i为应变引起的位移。

需要说明的是，为了保证页岩流动方程组的可解性，需要建立初始条件和边界条件，其初始条件为：当t＝0时，模型内基质、天然裂缝以及大尺度裂缝内各点压力(p_m,p_f,p_hf) 等于原始地层压力，即：

由此，获得页岩流动方程组内各岩石骨架的初始位移(u)及速度

均为0，即：

预设所述目标区域的气藏为定容封闭气藏，则外边界法线方向的压力梯度为0；水平井定压生产，则内边界压力等于井底流压，存在以下压力边界：

于此同时，设定储层底部基岩骨架无法移动，存在固定位移边界：

模型侧面岩石仅能垂直法向方向移动，存在垂直法向位移边界：

模型受上覆岩石压力，在顶面存在体积力载荷边界：

式中符号定义：Γ_i(i＝out，in，底面，侧面，顶面)分别为模型外、内、底面、侧面和顶面边界；p_原始为原始地层压力；F_i为岩石骨架所受体积力；n为法向向量；p_上覆为上覆压力；p_wf为水平井井底流压。

联立上述式(8)、(15)、(21)和(22)，即可获得所述页岩气非线性渗流模型，而上述式(23)～(28)则为页岩气非线性渗流方程的定解条件。

在一些实施例中，根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，具体包括以下步骤：

S31)任意选取一四面体网络作为控制体单元，大尺度裂缝处理相邻两个控制体单元之间的表面；

S32)根据有限体积方法分别对基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程进行线性离散，分别获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式，建立代数方程；所述代数方程为：

其中，a_C、a_F分别为选定的控制体单元和与其相邻的控制体单元中的系数；T_C、T_F分别为选定的控制体单元和与其相邻的控制体单元中的未知数；b_C为余量项；C表示选取的控制体单元；指标F～NB(C)表示所有与控制体单元C相邻的控制体单元；

S33)循环重复步骤S32)，获得每一四面体网络的代数方程，并建立由a_C、a_F组成的全局系数矩阵和由T_C、T_F组成的未知数向量矩阵。

需要说明的是，所述大尺度裂缝处理两个控制体单元之间的表面意为：由于裂缝面长宽比极大的特点(一般裂缝长度可达米级，而开度一般仅为毫米或微米级)，如果将裂缝面处理成三维空间实体，那么需要采用微米或毫米级网格对其进行剖分，将造成极大的计算资源浪费，并且难以保证模型收敛性，因此将所述目标区域内的大尺度裂缝降维处理成二维平面，对相邻的两个控制体单元进行剖分，即可获得包含多尺度特征的非结构四面体网络模型，减少计算量，提高计算效率。

需要说明的是，所述基质系统流动方程、天然裂缝流动方程和大尺度裂缝流动方程的半离散格式为：

所述页岩表观渗透系数κ定义为：

所述基质系统流动方程、天然裂缝流动方程和大尺度裂缝流动方程的半离散格式通过以下步骤获得，以基质系统流动方程为例：

在当前时间步下，将所述目标地层中的渗流视为稳定流动，忽略方程中的瞬时项，即可获得：

其中，V_c为选定的控制体单元体积；指标f～faces(V_c)表示构成控制体单元的所有界面； S_f为界面的表面向量；

任选一控制体单元作为当前控制体单元，参阅图4，即可获得控制体单元C中的渗流方程的稳定项：

其中，下指标1，2，3，4表示构成控制体单元C的4个界面，位置如图4所示；

根据求解方程组的定解条件，即可获得压力梯度算子及边界处压力：

其中，i＝m或f或hf。

将式(33)、(34)代入式(32)得到

参阅图4，A、B、D为控制体单元C的相邻控制体单元，其即为计算节点，由此式(35)等号右侧表达式中，下标为A、B、C、D的变量是可以直接调用的。因此，表面向量S₁、S₂、 S₃、S₄可以由式(36)计算，

其中：r_i(i＝1,2,3)，参阅图5，为三角形位置向量。

需要说明的是，为了提高计算效率和降低所述非结构四面体网络模型的模型容量，在所述非结构四面体网络模型中只储存计算计算节点处的数据，其他位置的数据通过相邻计算节点进行插值计算获得，而计算节点的数据是上个时间步保存的，因此可以直接调用。

将不能用计算节点表示的部分合并为一个参数团，对式(35)等号右侧表达式进行整理，可以得到形如式(37)的关系式：

其中，a_C、a_A、a_B、a_D、a_E分别为选定的控制体单元和与其相邻的控制体单元中的系数；T_C、T_F分别为选定的控制体单元和与其相邻的控制体单元中的未知数；b_C为余量项；C 表示选取的控制体单元；指标F～NB(C)表示所有与控制体单元C相邻的控制体单元；

其选定的控制体单元和与其相邻的控制体单元中的系数a_C、a_A、a_B、a_D、a_E和b_C为余量项分别为：

以图6为例，选定的控制体单元C相邻的三个控制体单元A、B、D，将式(39)展开即可控制体单元C的代数方程为：

当对所述目标区域完成网格剖分后，需要对每一个网格所对应的控制体进行编号以确定控制体之间的连接关系，任选取一网格，以图6为例，则可以将代数方程改写为：

a₃₃T₃+a₃₁T₁+a₃₂T₂+a₃₅T₅＝b₃ (41)

对系统内所有控制体单元重复上述过程，并将得到的代数方程组合成代数方程组，得到图6所示网格对应的系统方程组如下：

设定所述目标区域被分为了N个网格，由于每个代数方程描述当前网格和其他N-1个网格的连接关系，也就是每个方程里都有N项，当前网格与其他没有连接关系的网格所对应的项则为0；对应N个网格就有N个代数方程，就可以构造一个N×N的代数方程组；对于具有N个网格的系统，可以构造具有N×N系数矩阵的线性方程组:

A_N×NT_N×1＝B_N×1 (43)

式中符号定义：A为系数矩阵；T为未知数向量；B为余量向量。

由于四面体网格中每一个控制体单元最多有4个相邻控制体单元，因此A中绝大部分元素为0，一定是稀疏矩阵，可以采用共轭梯度法对该系数矩阵进行求解，获得该时间步内的压力解。

其中，由a_C、a_F组成的全局系数矩阵和由T_C、T_F组成的未知数向量矩阵并将每个代数方程的系数存储在对应的行列位置。单个控制体单元的偏微分方程离散是在控制体单元及其相邻控制体单元的局部范围内完成的，通过系统方程中的系数矩阵反应每个控制体的真实位置与控制体间的拓扑信息，利用额外储存的网格拓扑信息对系数矩阵进行建立。

通过所述目标区域的根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得的线性方程组即可得到所述目标区域的水平井产量分布。其具体求解过程如下：

S40)利用差分原理，用差商代替微分，即将

写作

的形式，表示时间步长内压力变化率，其中△t为时间步长，△p_m为时间步长内压力变化；

S41)添加时间步上标k(k＝0时为初始时刻)，预设当前时间步为k，获得在时间步为 k+1时，所述基质系统流动方程和大尺度裂缝流动方程的半离散格式；

S42)首先采用共轭梯度法尝试求解代数方程组，由于代数方程组存在大量余量项，可能导致线性化程度不足以满足收敛要求，此时采用Newton-Raphson迭代方法求解，获得

和

的解；并将求解得到的

和

的解带入所述天然裂缝流动方程的半离散格式进行显示求解，获得天然裂缝压力在k+1时间步的解

S43)将步骤S41)和步骤S42)获得的

的带入所述求解方程组，更新应力状态，获得下一时间步的

的解；

S44)循环步骤S43)即可获得每一时间步的

的解，体积应力、水平井产量值，进而获得所述目标区域的水平井产量分布情况。

需要说明的是：在时间步为k+1时，所述基质系统流动方程和大尺度裂缝流动方程的半离散格式为：

其中，在时间步为k+1时，所述天然裂缝流动方程的半离散格式为：

参阅图7，将采用上述步骤计算得到的水平井产量值与斯伦贝谢公司的Eclipse油藏数值模拟软件的模拟结果进行对比，其对比结果如图7所示，由此可见，本发明公开的页岩气藏产能模拟方法模拟得到的水平井产量值与Eclipse油藏数值模拟软件的模拟结果吻合度良好，说明本发明公开的页岩气藏产能模拟方法可以较好的模拟气藏产能；参阅图8、图9，考虑骨架应力应变后发现，骨架收缩可以增大页岩气产量(附图8)，修改页岩骨架弹性参数对产气量影响如图9所示，在深层页岩气藏中，天然裂缝和大尺度裂缝对水平井常量影响显著，不能忽略。

以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种页岩气藏产能数值模拟方法，其特征在于，包括以下步骤：

根据目标区域内的地质模型数据体生成非结构四面体网络模型；

建立基质系统流动方程；并引入天然裂缝骨架体积应变和大尺度裂缝骨架体积应变，建立天然裂缝流动方程、大尺度裂缝流动方程；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程建立页岩流动方程组；

根据补充广义胡克定律、质量守恒方程和应力平衡方程建立天然裂缝骨架体积应变和大尺度裂缝骨架体积应变的求解方程组；根据基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程和求解方程组建立页岩气非线性渗流模型，并确定页岩气非线性渗流模型的定解条件；

根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式；

根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得所述目标区域的水平井产量分布；

所述基质系统流动方程为：

其中，

为基质孔隙度；

为页岩解吸气量；

为基质与天然裂缝窜流量；M_g为气体摩尔质量；p_m为基质压力；R为理想气体常数；T为温度；Z为真实气体偏差因子；α为稀薄气体系数；δ^m为页岩气基质孔隙阻塞系数；μ_g为气体黏度；τ^m为页岩气基质孔隙迂曲度；c_g为气体压缩系数；D_g为气体扩散系数；k_m为基质的有效渗透率；p_L为Langmuir压力；ρ_m为吸附相密度；

为吸附气量；V_std为标准状况下气体摩尔体积；V_L为Langmuir体积；p_f为天然裂缝压力；Kn为Knudsen数；σ为裂缝形状因子；

所述天然裂缝流动方程为：

其中，k_f为天然裂缝的有效渗透率、

为天然裂缝孔隙度；K_f为天然裂缝骨架弹性系数；β为Biot弹性系数；ε_f为天然裂缝骨架的体积应变；σ为裂缝形状因子；q_f-h为天然裂缝向大尺度裂缝窜流量；

所述大尺度裂缝流动方程为：

其中，ε_hf为大尺度裂缝骨架的体积应变；q_well为水平井产量；k_hf为大尺度裂缝的有效渗透率；

为大尺度裂缝孔隙度；K_hf为大尺度裂缝骨架弹性系数；S_hf为大尺度裂缝综合弹性系数；d_hf为大尺度裂缝无因次开度因子；q_f-h为天然裂缝向大尺度裂缝窜流量；p_hf为大尺度裂缝压力；

所述页岩流动方程组为：

所述天然裂缝骨架体积应变和大尺度裂缝骨架体积应变的所述求解方程组为：

其中，δ_ij为克罗内克函数，若i＝j，则δ_ij＝1，若i≠j，则δ_ij＝0；ε_ij岩石骨架垂直于ij平面方向的体积应变，i,j＝x,y,z；μ为气体黏度；σ_ij为岩石骨架垂直于ij平面方向的应力分量，i,j＝x,y,z；σ_ij,j为应力张量的散度；F_i为岩石骨架所受体积力；u_i,j和u_j,i为应变引起的位移；

所述基质系统流动方程、天然裂缝流动方程和大尺度裂缝流动方程的半离散格式为：

其中，所述页岩表观渗透系数κ定义为：

V_c为以C点为形心的控制体体积；V_m为吸附气量；f～faces(V_C)为组成控制体V_C的各个界面；f～faces(V_C-F)为控制体V_C与大尺度裂缝共用的界面；S_f为界面f上指向外侧的方向向量；σ为裂缝形状因子；

为裂缝孔隙度；λ_f为裂缝流度；λ_hf大尺度裂缝流度；b～bounds(f)为构成界面f的各个边界；S_b边界f上指向外侧的方向向量；q_m-f为基质与天然裂缝窜流量；q_f-h为天然裂缝与大尺度裂缝窜流量；q_well为气井产量。

2.根据权利要求1所述的页岩气藏产能数值模拟方法，其特征在于，根据所述页岩气非线性渗流模型求解所述非结构四面体网络模型中每一四面体网络，具体包括以下步骤：

S31)任意选取一四面体网络作为控制体单元，大尺度裂缝处理相邻两个控制提单元之间的表面；

S32)根据有限体积方法分别对基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程进行线性离散，分别获得所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式，建立代数方程；所述代数方程为：

其中，a_C、a_F分别为选定的控制体单元和与其相邻的控制体单元中的系数；T_C、T_F分别为选定的控制体单元和与其相邻的控制体单元中的未知数；b_C为余量项；C表示选取的控制体单元；指标F～NB(C)表示所有与控制体单元C相邻的控制体单元；

S33)循环重复步骤S32)，获得每一四面体网络的代数方程，并建立由a_C、a_F组成的全局系数矩阵和由T_C、T_F组成的未知数向量矩阵。

3.根据权利要求2所述的页岩气藏产能数值模拟方法，其特征在于，根据所述基质系统流动方程、天然裂缝流动方程、大尺度裂缝流动方程的半离散格式和求解方程组获得所述目标区域的水平井产量分布，具体包括以下步骤：

S41)预设当前时间步为k，获得在时间步为k+1时，所述基质系统流动方程和大尺度裂缝流动方程的半离散格式；

S42)采用Newton-Raphson迭代方法求解，获得

和

的解；并将求解得到的

和

的解带入所述天然裂缝流动方程的半离散格式进行显示求解，获得天然

的解；

S43)将步骤S41)和步骤S42)获得的

的带入所述求解方程组，更新应力状态，获得下一时间步的

的解；

S44)循环步骤S43)即可获得每一时间步的

的解，体积应力、水平井产量值，进而获得所述目标区域的水平井产量分布情况。

Images