技术领域
[0001] 本
发明涉及油藏模拟技术领域，具体是一种结合连接对理论的嵌入式离散裂缝模型的数学推导和数值计算方法。
背景技术
[0002] 在对非常规油气资源的开发时，通常需要结合了多级压裂
水平井技术。由于裂缝性储层渗流规律极其复杂，因此研究这类多级压裂水平井的储层渗流规律通常采用数值模拟方法。多重介质数值模拟方法分为连续介质模型和离散裂缝模型，但是由于连续介质模型只适合大量微裂缝联通的储层，并不适合具有强
各向异性的水
力压裂储层，因此多用离散裂缝模型对其进行研究。嵌入式离散裂缝模型属于离散裂缝模型的分支，其裂缝网格是由基质网格边界剖分得到的，因此裂缝附近并不需要如同常规离散裂缝模型那样的网格加密，所以具有较快的运算速度。然而，在运用嵌入式离散裂缝模型进行数值模拟研究时，由于具有基质网格、裂缝网格、井网格多种流动系统的网格，且各类网格内部和各类网格之间都会发生
质量交换，因此每个网格的渗流方程都不相同，这也就造成了在运用嵌入式离散裂缝模型进行数学推导和数值计算过程中将面临较大困难。
发明内容
[0003] 本发明目的在于结合网格连接对理论，提出一种嵌入式离散裂缝模型的数学推导和数值计算方法，为嵌入式离散裂缝模型的应用以及复杂裂缝性油气藏的开发提供坚实的理论
基础。
[0004] 针对上述问题，本发明提出的嵌入式离散裂缝模型的数学推导和数值计算方法的技术方案如下：
[0005] 首先，为了方便
数据处理及方法阐述，本发明对嵌入式离散裂缝模型进行假设：
[0006] 1)储层中为单相渗流；
[0007] 2)考虑重力的影响；
[0008] 3)
生产井为多级压裂水平井模型，忽略基质流入井筒的
流体；
[0009] 4)生产制度为定压生产，且忽略生产井中的压降，即所有井网格压力恒定且相同；
[0010] 5)基质系统和裂缝系统渗透率恒定。
[0011] 同时对数值模拟方法也进行限定：
[0012] 1)采用
块中心网格，即基质和裂缝网格几何中心处的参数代表整个网格的属性；
[0013] 2)采用全隐式方法离散控制方程；
[0014] 3)采用中心加权处理传导率中与压力有关的弱非线性项。
[0015] 为了对嵌入式离散裂缝模型进行数学推导和数值计算，本发明的主要步骤包括：
[0016] 步骤一：结合连接对理论与全隐式有限差分方法，对渗流微分方程进行推导和离散，本步骤的具体操作为：将基质网格和裂缝网格的渗流微分方程写成如下形式，其中，方程左侧为质量交换项，方程右侧为质量累积项(以下流量交换都以流入该网格为正，流出为负)：
[0017]
[0018]
[0019] 其中：
[0020] Vb表示基质网格块的体积，
[0021] αc表示体积转换因子，在采用公制单位时，取1，
[0022] Δt表示两个时间点的时间差，
[0023] φm、φf表示基质网格和裂缝网格内的孔隙体积与基质网格块的体积之比，[0024] Bm、Bf表示基质网格和裂缝网格内流体的体积系数，
[0025] ∑qmm表示Δt时间内从与该基质网格相邻的所有基质网格流入该基质网格的总流量，
[0026] ∑qmf表示Δt时间内从嵌入该基质网格的所有裂缝网格流入该基质网格的总流量，
[0027] ∑qf表示Δt时间内从该裂缝网格所在裂缝面内的所有相邻裂缝网格流入该裂缝网格的总流量，相邻的裂缝网格即在同一条裂缝内且具有相同的
顶点的裂缝网格，[0028] qfm表示Δt时间内从该裂缝网格嵌入的基质网格流入该裂缝网格的流量，[0029] ∑qff表示Δt时间内从该裂缝网格相交的所有裂缝网格流入该裂缝网格的总流量，其中，相交的裂缝网格需要在该裂缝网格嵌入的基质网格内，
[0030] qfw表示Δt时间内从该裂缝网格穿过的井网格流入该裂缝网格的流量，当裂缝与生产井相连时，流体从裂缝网格流入井网格，本项为负数，
[0031] 表示Δt时间内该基质网格的质量净增加量，
[0032] 表示Δt时间内该裂缝网格的质量净增加量。
[0033] 将质量交换项写成传导率与势差之积的统一格式，其中势差由压力造成的势差和重力造成的势差两个部分组成，将两部分引起的质量交换项分开，从而可以将质量交换项写成以下统一格式：
[0034] q＝T(Φ2-Φ1)
[0035] ＝Tp(p2-p1)+TZ(Z2-Z1)
[0036] ＝qp+qZ
[0037] 其中：
[0038] q表示两网格间的质量交换，
[0039] Φ2-Φ1表示两个网格之间的势差，
[0040] T表示势差的传导率，
[0041] p2-p1表示两个网格之间的压差，
[0042] Tp表示压差的传导率，
[0043] Z2-Z1表示两个网格之间的高度差，
[0044] TZ表示高度差的传导率，
[0045] qp表示两网格间压差而产生的质量交换，
[0046] qZ表示两网格间高度差而产生的质量交换；
[0047] 渗流方程中的质量交换项q表示发生在n时步(已知时步)与n+1时步(未知时步)之间的Δt时间内的质量交换，全隐式方法是将求取质量交换q的参数全设置为n+1时步的未知量，而n+1时步的参数值由
迭代得到，当满足迭代收敛条件时，则可认为第v+1次的迭代值为n+1时步的值。基于全隐式方法，将上述qp和qZ2个部分分别展开，质量交换项展开为如下形式：
[0048]
[0049] 其中：
[0050] qn+1表示n+1时步的质量交换值，
[0051] q(v+1)、q(v)表示全隐式法中第v+1次和第v次迭代步的质量交换值，
[0052] 表示全隐式法中第v次迭代步与第v+1次迭代步间，由压差而产生的质量交换qp在两次迭代步间的差值，
[0053] 表示全隐式法中第v次迭代步与第v+1次迭代步间，由高度差而产生的质量交换qZ在两次迭代步间的差值；
[0054] 按照连接对理论，对上述公式中的 和 继续展开，如下所示：
[0055]
[0056]
[0057] 并将渗流方程右侧的积累项按照全隐式格式进行展开，可得：
[0058]
[0059] 其中
[0060]
[0061] 步骤二：对网格属性与基本参数进行赋值与计算，其中需要赋值和计算的参数包括：
[0062] 1)储层信息
[0063] 储层的长度、宽度、高度；
[0064] 2)网格信息
[0065] 基质网格的大小、基质网格中心点的坐标，裂缝网格的多边形端点坐标、水平井网格的两个线段端点坐标、基质网格编号、裂缝网格编号、井网格编号；
[0066] 3)网格及流体属性
[0067] 基质网格流体压力、基质网格流体
粘度、基质网格流体压缩系数、基质网格流体
密度、基质网格孔隙度、裂缝网格流体压力、裂缝网格流体粘度、裂缝网格流体压缩系数、裂缝网格流体密度、裂缝网格孔隙度。
[0068] 步骤三：基于网格连通性理论，根据网格
位置建立网格连通性表，包括相邻基质网格的连接对、裂缝网格与基质网格之间的连接对、一条裂缝内部的两个相邻裂缝网格之间的连接对、不同裂缝中的裂缝网格相交的连接对、裂缝网格和井网格之间的连接对，并记录计算网格间传导率需要的参数。
[0069] 步骤四：结合连通性表，根据嵌入式离散裂缝模型理论推导并计算网格间传导率的几何参数部分，根据中心加权方法推导并计算网格间传导率的弱非线性部分，其具体方法如下：
[0070] 单相流的总传导率T、压差的传导率Tp和高度差的传导率TZ的计算通式都可展开为：
[0071] T＝G·fp(p)＝G·fp
[0072] 式中G为几何因子，fp为与压力有关的弱非线性项；
[0073] G和fp的具体定义如下：
[0074]
[0075]
[0076] 其中：βc表示传导率转换因子，当采用公制单位时βc＝84.4×10-6；
[0077] k表示网格的渗透率；
[0078] A表示网格交面的面积；
[0079] Δl表示网格间的距离；
[0080] ρ表示流体密度；
[0082] μ表示流体粘度；
[0083] B表示流体体积系数。
[0084] 对于不同类型连接对的几何参数G计算方法各有不同，计算方法如下：
[0085] (1)基质网格与其相邻基质网格连接对的几何参数Gmm：
[0086]
[0087] 其中：km1、km2——基质网格1和基质网格2的渗透率；
[0088] Amm——两基质公共面面积；
[0089] dm1、dm2——基质网格1和基质网格2的形心到公共面的的平均距离；
[0090] (2)基质网格和其网格内部的裂缝网格连接对的几何参数Gmf：
[0091]
[0092] 其中：Amf——裂缝网格与基质网格交面的面积，由于裂缝面两侧均与基质相连，因此交面面积等于裂缝段面单侧面积的两倍；
[0093] kmf——裂缝和基质连接对的，由于裂缝渗透率远大于基质渗透率，因此连接对的渗流能力主要由渗透率较小的基质所决定，即kmf≈km；
[0094] dmf——平均
正交距离；表示基质网格与裂缝网格间的等效距离，即基质网格内所有点到裂缝段垂直距离的平均值；
[0095] (3)一条裂缝内部的两个相邻裂缝网格之间连接对的几何参数Gf：
[0096]
[0097] 其中： ——裂缝网格1和裂缝网格2的渗透率；
[0098] Af——两裂缝网格公共面面积；
[0099] ——裂缝网格1和裂缝网格2的形心到公共面的的平均距离；
[0100] (4)不同裂缝中的裂缝网格相交连接对的几何参数Gff：
[0101]
[0102] 其中： ——裂缝网格1和裂缝网格2的渗透率；
[0103] ——裂缝网格1和裂缝网格2的开度；
[0104] Lint——裂缝交线长度；
[0105] ——裂缝网格1和裂缝网格2的裂缝中心到两裂缝交线的距离；
[0106] (6)裂缝网格与井网格的连接对的几何参数Gwf：
[0107]
[0108] 其中：kf——裂缝渗透率；
[0109] ωf——裂缝开度；
[0110] rw——井半径；
[0111] Lf、hf——裂缝段的长度、高度；
[0112] Δθ——径向井包含在裂缝内的圆心
角角度，一般情况下，井眼在某一条裂缝中只穿过一个裂缝网格时，Δθ＝2π。
[0113] fp采用中心加权的方法，即对于网格1和网格2形成的连接对，两网格界面处的压力参数由连接对两个网格中心处的压力得到，因此界面上的弱非线性项fp为两侧网格压力的函数，即：
[0114]
[0115] 其中：
[0116] 对于 的计算采用割线代替切线(梯度)，即：
[0117]
[0118] 其中Δp表示压力p一个极小的变化量。
[0119] 步骤五：结合连通性表，组装差分方程的系数矩阵。对单相流数值模拟来说，待求未知量为所有网格在两个迭代步之间的压力差。因此需要对未知量进行排序，确定未知量的书写顺序，即对所有网格排序。本方法按照基质网格编号从小到大，再裂缝网格编号从小到大的顺序，对所有网格进行排序，从而完成了对所有网格未知压差的排序。由于本模型假设是定压生产，所以井网格的压力恒定，即所有时步的压力差变化都已知，且变化量为零，井网格压差为常数，不需要对井网格进行排序。再改写为矩阵方程的形式：
[0120]
[0121] 其中A表示系数矩阵， 表示未知压力差所组成的向量，b表示常数所组成的向量。
[0122] 若第i和第j个未知数所代表的两个网格有连接对关系，则此连接对之间的质量交换在全隐式方法中将展开为以下的形式：
[0123]
[0124] 上式的展开结果显示出第i和第j个网格间的质量交换可以展开由第i个未知数及其系数、第j个未知数及其系数、常数组成的线性方程。其中 和 分别位于系数矩阵A的第i行的第i列和第i行的第j列，bij位于常数向量b的第i行。
[0125] 按照以上方法对所有连接对进行展开，并将展开后的系数和常数放入
指定位置，即完成了对质量交换项的处理。
[0126] 接下来再对每个网格的质量累积项进行处理，例如第i个网格的展开形式如下[0127]
[0128] 上式表示第i个网格的质量累积项可以展开为第i个未知数及其系数、常数组成的线性方程。其中ai位于系数矩阵A第i行的第i列，bi位于常数向量b的第i行。
[0129] 由以上步骤对所有网格的质量累积项进行展开，并将展开后的系数和常数放入指
定位置，即完成了对质量累积项的处理；
[0130] 以上即完成了矩阵方程的组装，之后就可以对嵌入式离散裂缝模型进行计算，从而通过迭代得到下一个时步各个网格的压力。
[0131] 本发明的有益效果是：在嵌入式离散裂缝模型中，每个裂缝网格、基质网格连接的网格种类和数量都未知，因此每个网格所对应的渗流方程也各不相同，每个渗流方程的有限差分形式也各不相同，因而在组装矩阵时无法形成常规的矩阵形式。本发明从连接对理论入手，不需要获得每个网格的渗流方程的有限差分形式，只需要推导出几种网格连接对的质量交换项和网格积累项的有限差分形式，并通过连接对信息，直接组装矩阵方程，大大简化了难度和工作量，并形成了一套对于嵌入式离散裂缝模型数学推导和数值计算的具有普适性的方法，解决了嵌入式离散裂缝模型计算难度大的问题。
具体实施方式
[0133] 实施例：
[0134] 在全隐式有限差分方法中，渗流方程都将会被离散化、线性化处理为多项式方程组，方程组的未知数为各个网格相邻时步的压力差，方程组中的每一个多项式都由一个网格的与其他网格的质量交换以及本网格最终的质量积累所组成。通常会将多项式方程组改写成矩阵方程的形式。
[0135] 对于将控制方程转化为矩阵方程的整个过程来说，常规的是对所有网格依次写出其离散后的差分方程，再列出其多项式形式，最后组装为矩阵。然而，对嵌入式离散裂缝模型来说，其中含有基质和裂缝两套网格，且每个网格所连接的网格类型、数量都未知，因此对每个网格直接写出或通过网格连接规律间接写出渗流方程和多项式都是无法完成。
[0136] 基于以上难点，本发明结合连接对理论，分别对网格间的质量交换项和网格积累项进行推导与展开，并直接利用连接对的网格联通信息组装矩阵方程的系数矩阵，而不在意每个网格的渗流差分方程的具体形式，由此进行模型推导以及数值计算。基于以上思路方法，本方法分为以下五个步骤。
[0137] 步骤一
[0138] 本步骤将结合连接对理论与全隐式有限差分方法，对渗流微分方程进行推导和离散。由于嵌入式离散裂缝模型中包括多种流动系统，且流动系统之间都存在窜流，因此将基质网格和裂缝网格的渗流微分方程写成如下形式，方程左侧为质量交换项，方程右侧为质量累积项，以下流量交换都以流入该网格为正，流出为负：
[0139]
[0140]
[0141] 其中：
[0142] Vb表示基质网格块的体积，
[0143] αc表示体积转换因子，在采用公制单位时，取1，
[0144] Δt表示两个时间点的时间差，
[0145] φm、φf表示基质网格和裂缝网格内的孔隙体积与基质网格块的体积之比，[0146] Bm、Bf表示基质网格和裂缝网格内流体的体积系数，
[0147] ∑qmm表示Δt时间内从与该基质网格相邻的所有基质网格流入该基质网格的总流量，
[0148] ∑qmf表示Δt时间内从嵌入该基质网格的所有裂缝网格流入该基质网格的总流量，
[0149] ∑qf表示Δt时间内从与该裂缝网格所在裂缝的所有相邻裂缝网格流入该裂缝网格的总流量，
[0150] qfm表示Δt时间内从该裂缝网格嵌入的基质网格流入该裂缝网格的流量，[0151] ∑qff表示Δt时间内从与该裂缝网格相交的所有裂缝网格流入该裂缝网格的总流量，其中，相交的裂缝网格需要在该裂缝网格嵌入的基质网格内，
[0152] qfw表示Δt时间内从该裂缝网格穿过的井网格流入该裂缝网格的流量，[0153] 表示Δt时间内该基质网格的质量净增加量，
[0154] 表示Δt时间内该裂缝网格的质量净累加量；
[0155] 将质量交换项写成传导率与势差之积的统一格式，其中势差由压力造成的势差和重力造成的势差两个部分组成，将两部分引起的质量交换项分开，从而可以将质量交换项写成以下统一格式：
[0156] q＝T(Φ2-Φ1)
[0157] ＝Tp(p2-p1)+TZ(Z2-Z1)
[0158] ＝qp+qZ
[0159] 其中：
[0160] q表示两网格间的质量交换，
[0161] Φ2-Φ1表示两个网格之间的势差，
[0162] T表示势差的传导率，
[0163] p2-p1表示两个网格之间的压差，
[0164] Tp表示压差的传导率，
[0165] Z2-Z1表示两个网格之间的高度差，
[0166] TZ表示高度差的传导率，
[0167] qp表示两网格间压差而产生的质量交换，
[0168] qZ表示两网格间高度差而产生的质量交换；
[0169] 渗流方程中的质量交换项q表示发生在n时步(已知时步)与n+1时步(未知时步)之间的Δt时间内的质量交换，全隐式方法是将求取质量交换q的参数全设置为n+1时步的未知量，而n+1时步的参数值由迭代得到，当满足迭代收敛条件时，则可认为第v+1次的迭代值为n+1时步的值。基于全隐式方法，将上述qp和qZ2个部分分别展开，质量交换项展开为如下形式：
[0170]
[0171] 其中：
[0172] qn+1表示n+1时步的质量交换值；
[0173] q(v+1)、q(v)表示全隐式法中第v+1次和第v次迭代步的质量交换值；
[0174] 表示全隐式法中第v次迭代步与第v+1次迭代步间，由压差而产生的质量交换qp的差值；
[0175] 表示全隐式法中第v次迭代步与第v+1次迭代步间，由重力差而产生的质量交换qZ的差值；
[0176] 按照连接对理论，对上述公式中的 和 继续展开，如下所示：
[0177]
[0178]
[0179] 并将渗流方程右侧的积累项按照全隐式格式进行展开，可得：
[0180]
[0181] 其中
[0182]
[0183] 步骤二
[0184] 嵌入式离散裂缝模型的网格中包括基质网格系统、裂缝网格系统和井网格系统三个流体流动系统，各个系统之间会发生流体交换，为了计算得到系统间的质量交换，必须先计算其传导率。而计算传导率之前又必须先构建连通性表。因此在本步骤中应赋值或计算的参数包括：
[0185] 1)储层信息
[0186] 储层的长度、宽度、高度；
[0187] 2)网格信息
[0188] 基质网格的大小、基质网格中心点的坐标，裂缝网格的多边形端点坐标、水平井网格的两个线段端点坐标、基质网格编号、裂缝网格编号、井网格编号；
[0189] 3)网格及流体属性
[0190] 基质网格流体压力、基质网格流体粘度、基质网格流体压缩系数、基质网格流体密度、基质网格孔隙度、裂缝网格流体压力、裂缝网格流体粘度、裂缝网格流体压缩系数、裂缝网格流体密度、裂缝网格孔隙度。
[0191] 步骤三
[0192] 本步骤需要建立连通性表，连通性表由所有连接对及其相关信息组成。因此在建立连通性表前需要找到所有网格连接对。对于各类连接对的确定，分别采用了不同的方法。
[0193] (1)基质网格与其相邻的基质网格连接对
[0194] 基质网格都为结构网格，因此当两个基质网格相连时，则只会出现为上-下、左-右、前-后三种联通方式。若i、j两个基质网格的位置坐标为(xi,yi,zi)、(xj,yj,zj)，即i、j网格在x、y、z方向上分别是第xi、yi、zi和xj、yj、zj个网格，则当两个网格相邻时，在数学上有以下关系：
[0195] |(xi,yi,zi)-(xj,yj,zj)|＝1
[0196] (2)基质网格和其网格内部的裂缝网格连接对
[0197] 由于裂缝网格是由基质网格边界剖分出来的，因此若裂缝网格的中点在某个基质网格内部，则该裂缝网格就和这个基质网格形成连接对。
[0198] (3)一条裂缝内部的两个相邻裂缝网格之间的连接对
[0199] 对此类裂缝网格之间的连接对，首先需要判断两个裂缝网格是否在同一条裂缝内部。对于三维模型，因为裂缝都被
降维成二维，因此只需判断同一条裂缝内部的两个网格是否有具有两个相同的顶点。若有两个相同顶点，则这两个裂缝网格形成连接对。
[0200] (4)不同裂缝中的裂缝网格相交的连接对
[0201] 此连接对的两裂缝网格处于不同的裂缝中，因此首先判断两裂缝面是否相交，若相交，再判断两裂缝面的交线是否处于裂缝边界内，若在裂缝边界内，则两条裂缝相交。若两个裂缝网格所在的裂缝相交，且交线穿过这两个裂缝网格，则这两个裂缝网格形成连接对。
[0202] (5)裂缝网格与井网格的连接对
[0203] 首先判断某裂缝面与井所在直线是否相交，若相交，再判断该交点是否在井始末两点构成的线段内，若在线段内，则裂缝与井相交。然后再确定交点属于哪一个裂缝网格和井网格，如此就得到了该裂缝网格和井网格形成连接对。
[0204] 步骤四
[0205] 本步骤为计算传导率，因此需要计算传导率中G和fp两个参数。
[0206] 得到步骤三中的连接对之后，记录连接对中基质网格、裂缝网格、井网格的编号，并记录传导率计算过程中需要的参数。单相流的总传导率T、压差的传导率Tp和高度差的传导率TZ的计算通式都可展开为：
[0207] T＝G·fp(p)＝G·fp
[0208] 式中G为几何因子，fp为与压力有关的弱非线性项；
[0209] G和fp的具体定义如下：
[0210]
[0211]
[0212] 其中：βc表示传导率转换因子，当采用公制单位时βc＝84.4×10-6；
[0213] k表示网格的渗透率；
[0214] A表示网格交面的面积；
[0215] Δl表示网格间的距离；
[0216] ρ表示流体密度；
[0218] μ表示流体粘度；
[0219] B表示流体体积系数。
[0220] 对于不同类型连接对的几何参数G计算方法各有不同，计算方法如下：
[0221] (1)基质网格与其相邻基质网格连接对的几何参数Gmm：
[0222]
[0223] 其中：km1、km2——基质网格1和基质网格2的渗透率；
[0224] Amm——两基质公共面面积；
[0225] dm1、dm2——基质网格1和基质网格2的形心到公共面的的平均距离；
[0226] (2)基质网格和其网格内部的裂缝网格连接对的几何参数Gmf：
[0227]
[0228] 其中：Amf——裂缝网格与基质网格交面的面积，由于裂缝面两侧均与基质相连，因此交面面积等于裂缝段面单侧面积的两倍；
[0229] kmf——裂缝和基质连接对的，由于裂缝渗透率远大于基质渗透率，因此连接对的渗流能力主要由渗透率较小的基质所决定，即kmf≈km；
[0230] dmf——平均正交距离；表示基质网格与裂缝网格间的等效距离，即基质网格内所有点到裂缝段垂直距离的平均值；
[0231] (3)一条裂缝内部的两个相邻裂缝网格之间连接对的几何参数Gf：
[0232]
[0233] 其中： ——裂缝网格1和裂缝网格2的渗透率；
[0234] Af——两裂缝网格公共面面积；
[0235] ——裂缝网格1和裂缝网格2的形心到公共面的的平均距离；
[0236] (4)不同裂缝中的裂缝网格相交连接对的几何参数Gff：
[0237]
[0238] 其中： ——裂缝网格1和裂缝网格2的渗透率；
[0239] ——裂缝网格1和裂缝网格2的开度；
[0240] Lint——裂缝交线长度；
[0241] ——裂缝网格1和裂缝网格2的裂缝中心到两裂缝交线的距离；
[0242] (7)裂缝网格与井网格的连接对的几何参数Gwf：
[0243]
[0244] 其中：kf——裂缝渗透率；
[0245] ωf——裂缝开度；
[0246] rw——井半径；
[0247] Lf、hf——裂缝段的长度、高度；；
[0248] Δθ——径向井包含在裂缝内的圆心角角度，一般情况下，井眼在某一条裂缝中只穿过一个裂缝网格时，Δθ＝2π。
[0249] fp采用中心加权的方法，即对于网格1和网格2形成的连接对，两网格界面处的压力参数由连接对两个网格中心处的压力得到，因此界面上的弱非线性项fp为两侧网格压力的函数，即：
[0250]
[0251] 其中：
[0252] 对于 的计算采用割线代替切线(梯度)，即：
[0253]
[0254] 其中Δp表示压力p一个极小的变化量。
[0255] 步骤五
[0256] 本步骤将结合连通性表，组装差分方程的系数矩阵。对单相流数值模拟来说，待求未知量为所有网格在两个迭代步之间的压力差。因此需要对未知量进行排序，确定未知量的书写顺序，即对所有网格排序。本方法按照基质网格编号从小到大，再裂缝网格编号从小到大的顺序，对所有网格进行排序，从而完成了对所有网格未知压差的排序。由于本模型假设是定压生产，所以井网格的压力恒定，即所有时步的压力差变化都已知，且变化量为零，井网格压差为常数，不需要对井网格进行排序。再改写为矩阵方程的形式：
[0257]
[0258] 其中A表示系数矩阵， 表示未知压力差所组成的向量，b表示常数所组成的向量。
[0259] 若第i和第j个未知数所代表的两个网格有连接对关系，则此连接对之间的质量交换在全隐式方法中将展开为以下的形式：
[0260]
[0261] 上式的展开结果显示出第i和第j个网格间的质量交换可以展开由第i个未知数及其系数、第j个未知数及其系数、常数组成的线性方程。其中 和 分别位于系数矩阵A的第i行的第i列和第i行的第j列，bij位于常数向量b的第i行。
[0262] 由以上方法对所有连接对进行展开，并将展开后的系数和常数放入指定位置，即完成了对质量交换项的处理。
[0263] 接下来再对每个网格的质量累积项进行处理，例如第i个网格的展开形式如下[0264]
[0265] 上式表示第i个网格的质量累积项可以展开为第i个未知数及其系数、常数组成的线性方程。其中ai位于系数矩阵A第i行的第i列，bi位于常数向量b的第i行。
[0266] 由以上方法对所有网格的质量累积项进行展开，并将展开后的系数和常数放入指定位置，即完成了对质量累积项的处理。
[0267] 以上即完成了矩阵方程的组装，之后就可以对嵌入式离散裂缝模型进行计算，从而得到下一个时步各个网格的压力。
[0268] 以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单
修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
