裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法

摘要

本发明公开了一种裂隙‑孔隙结构双重介质煤储层气水两相流数值模拟方法，研究了煤层气的运移和产出机理，建立了煤储层的地质模型。通过对渗流方程推导，建立气、水两相流数学模型，进一步对数学模型进行差分计算。该方法使用全隐式数值计算方法形成一套完整的煤储层气水两相流数值模拟方法。

说明书

技术领域

本发明属于煤层气的开采技术领域，具体地涉及一种裂隙-孔隙结构双重介质煤层气储层气水两相流数值模拟方法。

背景技术

煤层气是一种自生自储的非常规天然气，它的开发和利用对改善能源结构、环境保护和煤矿安全均具有重要意义。我国煤层气资源遍布全国，埋深小于2000m的煤层气资源总量为22.5×10¹²m³，具有广阔的开发前景。虽然我国煤层气的储藏量较大，但是与美国相比，我国绝大多数的煤层气储藏开采条件困难，普遍性的难点为：

1.低渗透率：煤层裂隙发育不完善，煤层原始渗透率低，与美国圣胡安盆地煤层渗透率1～50mD相比，我国难采煤层气区块渗透率仅为0.1～0.5mD，无强制作用下煤层气溢出速度缓慢。

2.欠饱和吸附煤层：地质结构变化频繁，导致特定阶段煤层气溢出地面，或者煤系发育不完全导致煤层气含量不丰富。目前已发现的煤层气高产富集区块均为高吸附饱和度，如圣胡安盆地为90％～98％，黑勇士盆地为92％～99％，沁水煤层气田为85％～95％，而我国难采煤层气区块煤层气吸附饱和度则低于30％。

3.难压储层，有效地应力高企：煤层有效地应力低的地区，其煤 层渗透率比相同条件下的高应力区的煤层渗透率要高，煤层有效地应力愈大，其压裂难度愈大，煤层压裂后裂隙闭合速度越快。而我国难采煤层气区块有效地应力相对较高，与圣胡安盆地高产区域地应力为3～8MPa，我国难采煤层气区块有效地应力为12～15Mpa以上。

4.煤层结构碎裂，煤质较软：长期地质扰动造成煤层反复挤压造成煤层内颗粒状结构，导致煤层整体抗压性较差，在压裂后形成的煤层裂隙稳定性较差。

以上均是阻碍我国难采煤层气区块煤层气开发较为常见的技术难题。

与常规油气相比，目前对煤层气的数值模拟研究相对较少，理论和方法相对薄弱。

发明内容

针对上述技术问题，本发明旨在提供一种裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法，可以指导煤层气开采。

为达到上述目的，本发明的技术方案是：

一种裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法，其特征在于，包括以下步骤：

S01：建立煤层气储集和运移的地质模型；

S02：根据煤层气吸附-扩散模型，基于渗流力学理论，建立煤储层数值模拟的数学模型，所述数学模型包括裂隙系统的基本微分方程、煤基质微孔隙系统的解吸吸附方程、辅助方程、定解方程；

S03：对上述数学模型进行离散差分求解，获取全隐式格式的数值模型；

S04：对上述数值模型采用预处理正交极小化方法进行求解，得出每个时段的气和水产量值。

优选的，所述步骤S02中裂隙系统的基本微分方程包括气相基本微分方程和水相基本微分方程，气相基本微分方程为：

水相基本微分方程为：

煤基质微孔隙系统的解析吸附方程为：

式中：其中G为几何因子，Vm表示煤基质中吸附气体的平均含量，为气体定产量项，其中h表示产层厚度，P_wfg表示煤层气井井底气体流压，r_e为排泄半径，r_w为井筒半径，S为表皮系数，K_x，K_y，K_z分别表示裂隙系统的绝对渗透率在X，Y，Z三个方向上的分量，Krg表示气相的相对渗透率，ρ_g表示气体密度，μ_g表示气体粘度，分别表示x，y，z三个方向上的气相压力梯度，D表示由某一基准面算起的深度，向下为正，g表示重力加速度，Sg表示气相饱和度，φ表示裂隙系统的储层孔隙度；

式中：为液体产量项，P_wfw表示煤层气井井底液体流压，Krw表示水相的相对渗透率，ρ_w表示地层水密度，μ_w表示地层水粘度，分别表示x，y，z三个方向上的水相压力梯度，S_w表示水相饱和度；

式中：σ表示形状因子，D表示煤基质的气体扩散系数，V_L表示兰缪尔体积，P_L表示兰缪尔压力；P_g表示自由气体压力，V_e与自由气体处于平衡状态的吸附气体含量；

辅助方程：

S_g+S_w＝1；

P_cgw(S_w)＝P_g-P_w，其中Pcgw为毛管压力，P_g和P_w分别为气相和水相压力；

ρ_g＝ρ_g(P_g)，气相密度；

ρ_w＝ρ_w(P_w)，水相密度；

K_rg＝K_rg(S_g)，气相相对渗透率；

K_rw＝K_rw(S_w)，水相相对渗透率；

μ_g＝μ_g(P_g)，气体粘度；

μ_w＝μ_w(P_w)，地层水粘度；

初始条件

P_w(x，y，z，t＝0)＝P_wi(x，y，z)，初始地层压力；

S_w(x，y，z，t＝0)＝S_wi(x，y，z)，初始含水饱和度；

内边界条件

P_rw＝P_wf(x，y，z，t)，内边界水相压力条件；

定压外边界条件

P_g|_E1＝P_g1(x，y，z，t)，外边界气相压力；

P_w|_E1＝P_w1(x，y，z，t)，外边界水相压力；

封闭外边界条件

外边界气相压力梯度；

外边界水相压力梯度。

优选的，所述步骤S03中数值模型为：

优选的，所述气相基本微分方程和水相基本微分方程由气、水相运动方程代入到气、水相连续性方程中得到；所述气相运动方程为：

所述水相运动方程为：

所述气相连续性方程为：

水相连续性方程为：

优选的，所述步骤S04包括，首先将三维中线z方向选定为“线”方向；再固定x坐标，形成一块x-z板；然后在板上按照二维剖面模型的扫描方法进行扫描；完成后进行第二块板的扫描，这样依次将所有的板都扫描完，并反复进行，直到得到满足精度要求的解，从而完成方程的一次求解过程，要在每一条线上的解得出后，施加松弛因子算出该线上各节点未知量的该次迭代值，随后再进行下一条线的计算。

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

1)考虑了裂隙-孔隙双重介质的煤储层中气、水两相流的情况；2)使用气体吸附-扩散模型；3)采用全隐式的数值插分方程；4)使用预处理压 交极小化方法求解方程组。解的稳定性好，结果精度高，可以更好地指导煤层气的开采。

附图说明

图1为本发明裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法的流程图；

图2为煤储层中的微小单元体结构模型图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式以及附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

实施例：

1.1建立地质模型，煤层气储集和运移的地质模型可概化如下：

(1)煤储层是非均质、各向异性的；

(2)煤储层具有典型的裂隙一孔隙双重孔隙结构，是由煤基质块中的微孔隙系统和裂隙系统构成的双孔隙系统；

(3)煤储层中存在气、水两相流体；

(4)煤基质块中的微孔隙是煤层气的主要储集空间，煤层气主要以吸附状态赋存于微孔隙的内表面上。微孔隙很小，水不能进入其中，所以微孔隙系统中只存在单相气体。裂隙系统既是气体和水的储集空间，也是气、水两相流体的渗流通道；

(5)煤层气体从煤储层运移产出经历渗流、解吸和扩散三个阶段，

①渗流：当煤层气井排水降压时，裂隙系统中的煤层气和水在压力梯度作用下以层流形式向井筒流动。这一流动服从Darcy定律；

②解吸：由于气、水从裂隙产出，导致煤储层压力下降，煤基质块表面微孔隙系统中吸附的煤层气发生解吸，使得煤基质块表面微孔隙系统中吸附的煤层气含量即煤层气浓度降低。此煤层气含量或煤层气浓度大小可用兰缪尔方程来描述；

③扩散：随着煤基质块表面微孔隙中煤层气浓度的降低，煤基质块的内部微孔隙系统与表面微孔隙系统之间形成煤层气浓度梯度，即其内部微孔隙系统中煤层气浓度相对较高，而其表面微孔隙中煤层气浓度相对较低。在浓度梯度的作用下，内部微孔隙系统中的煤层气通过微孔隙系统沿着浓度梯度降低方向由内部向表面发生扩散运移，在此视为拟稳态运移，受Fick第一定律控制。

(6)煤储层内气、水的运移过程是等温的。

(7)在裂隙系统中的气体是自由气体，表现为真实气体特性；水是可压缩的流体。在微孔隙系统中，各种气体组分不存在选择性吸附和扩散现象。

(8)基质块表面的解吸作用非常快，足以维持自由气和吸附气之间的平衡。在煤基质块内部的吸附气与自由气处于不平衡状态。

1.2煤层气储层模拟数学模型的建立，实际上就是根据上述概化的地质模型，用数学方法推导出描述煤储层中煤层气和水储集和运移规律的基本偏微分方程组，并给出定解条件。下面分别从裂隙系统、基质的微孔隙系统和煤层气井三个方面进行推导。

1.2.1裂隙系统基本微分方程包括连续性方程和运动性方程。

1.2.1.1裂隙系统连续性方程

在煤储层中取一个六面体的微小单元体，如图2所示，其中心点的坐标为(x，y，z)，其长、宽、高分别为Δx、Δy、Δz，各个侧面分别与x、y、z轴平行在裂隙系统中，流体从前面流入，后面流出；左面流入，右面流出；底面流入，顶面流出。流体的速度为V(x，y，z)，密度为ρ(x，y，z)，饱和度为S(x，y，z)，单元体内的孔隙度为

假设，在x方向，气体的流入速度和流出速度分别为V_gx|_x-Δx/2和V_gx|_x+Δx/2，在y方向，气体的流入速度和流出速度分别为V_gy|_y-Δy/2和V_gy|_y+Δy/2，在z方向，气体的流入速度和流出速度分别为V_gx|_x-Δx/2和V_gz|_z+Δz/2。取微小时间段Δt，考虑气体在微小单元体内的流入流出情况。

在Δt时间内，气体沿x方向流入和流出单元体的质量流量的差值为：

在Δt时间内，气体沿y方向流入和流出单元体的质量流量的差值为：

在Δt时间内，气体沿z方向流入和流出单元体的质量流量的差值为：

在Δt时间内，气体流入和流出单元体引起单元体内气体饱和度发生变化，导致单元体内气体的质量变化为：

根据质量守恒原理，气体流入和流出单元体的质量流量的差值应等于因单元体内气体饱和度变化引起的气体的质量增量：

用ΔxΔyΔzΔt同除方程两边，并令Δx→0，Δy→0，Δz→0，Δl→0，取极限，得到单位时间里裂隙系统的单位体积内气体质量变化的微分方程，即气相的连续性方程：

同理，可得水相的连续性方程为：

1.2.1.2裂隙系统的运动方程

在三维空间、各向异性介质、考虑重力影响情况下，Darcy定律推广形式为：

式中，v为渗流速度，是一个空间向量，k为孔隙介质的绝对渗透率；μ为流体的粘度，γ为流体的重度，γ＝ρg，ρ为流体的密度，g为重力加速度，D为由某一基准面算起的深度，向下为正，▽为Hamilton算子，负号表示沿流动方向的压力是下降的。

这时的渗透率K是一个2阶张量，写成矩阵形式为：

渗流速度v在三个方向上的分量分别为：

通常，渗透率张量的分量K_xy＝K_yz，K_yz-K_zy，K_zx＝K_xz，所以渗透率张量是对称张量。将坐标轴方向取得与介质中某点渗透率张量的主方向一致，则渗透率张量矩阵具有对角线形式：

用这种形式表示的张量称为对角线张量。在这种情况下，Darcy定律的推广形式为：

在气、水两相情况下，气、水的相对渗透率分别表示为K_rg、K_rw，气、水的密度分别表示为ρ_g、ρ_w，气、水的粘度分别表示为μ_g、μ_w，则气相、水相的Darcy定律的推广形式分别为：

1.2.1.3裂隙系统的基本微分方程

将运动方程(1-4)、(1-5)代入到连续性方程(1-1)、(1-2)中，得到裂隙系统中气、水两相渗流的基本微分方程：

1.2.2煤基质微孔隙系统的解析吸附方程推导

裂隙系统中的气体是自由气体，基质微孔隙中的气体则主要是吸附气体。在基质中，只有靠近裂隙面的基质微孔隙中的气体，解吸作用足够快，与自由气体处于平衡状态；而远离裂隙的基质微孔隙中气体与裂隙中的自由气体处于非平衡状态。与自由气体处于平衡状态的吸附气体含量可用兰缪尔模型求得：

式中，V_L为兰缪尔体积；P_L为兰缪尔压力；P_g为自由气体压力；V_e为与自由气体处于平衡状态的吸附气体含量。

在基质块内部和表面之间存在的气体浓度差作用下，基质块内部微孔隙中的气体以扩散方式向外部运移，进入裂隙系统中，可视为点源项来处理。

按拟稳态条件考虑，根据Fick第一定律，煤基质块的平均气含量对时间的变化率与煤基质块平均气含量和其表面吸附气体含量之差成正比，而单位时间内由单位煤基质解吸扩散进入到裂隙系统的气体量与煤基质块平均气含量的变化率成正比，即：

式中：V_m为煤基质中吸附气体的平均含量；

V_e为裂隙面上与自由气体压力处于平衡状态的吸附气体含量；

σ为Arren和Root的形状因子，与基质单元的尺寸大小和形状有关；

D为煤基质的气体扩散系数；

T为吸附时间常数，定义为τ＝1/σD；G为几何因子。

由于煤中割理系统的正交性，柱状基质几何体对于煤层是最合适的。对这种几何体，其形状因子定义为：

σ＝8/a²＝8π/s²

这里，s是割理的平均间距，a是柱状体的等效半径。因此，吸附时间τ可表示为：

τ＝s²/8πD

为了理解吸附时间的物理意义，对上述微分方程进行分离变量，并给出初边条件：

求解得：

V_m＝V_e(P_g)+[V_m-V_e(P_g)]e^-t/τ(1-11)

V_L-V_m＝V_L-V_e(P_g)-[V_L-V_e(P_g)]e^-t/τ

当t＝τ时

由此可见，吸附时间τ是指当解吸气量占总气量的63％时所对应的时间。

1.2.3辅助方程

为了完整地描述和求解气、水在煤储层中的运移过程，除了微分方程组外，还必须提供某些辅助方程来完善数学模型，它们是饱和度方程和毛管压力方程。

S_g+S_w＝1(1-13)

P_ugw(S_w)＝P_g-P_w(1-14)

鉴于求解变量为P_g、P_w、S_g、S_w，共4个，方程组的方程数也是4个，所以这个方程组是封闭的。当然，这个方程组的各项参数还需要由下述密度、相对渗透率、粘度和毛管压力的辅助方程来确定：

ρ_g＝ρ_g(P_g)

ρ_w＝ρ_w(P_w)

K_rg＝K_rg(S_g)

K_rw＝K_rw(S_w)

μ_g＝μ_g(P_g)

μ_w＝μ_w(P_w)

P_cgw＝P_cgw(S_w)

(1-15)

1.2.4边界条件

解上述方程组，还需要根据具体的情况给定边界条件和初始条件。边界条件和初始条件统称为定解条件。

煤层气储层数值模拟中的边界条件分为外边界条件和内边界条件两大类，其中外边界条件是指煤储层外边界所处的状态，内边界条件是指煤层气生产井所处的状态。

外边界条件：外边界条件一般有以下三类：

定压边界条件：外边界E上每一点在每一时刻的压力分布都是已知的，即为一已知函数，在数学上也称第一类边界条件，或称Dirchlet边界条件，表示为：

P_E1＝f₁(x，y，z，t)

定流量边界：外边界E上有流量流过边界，而且每一点在每一时刻的值都是已知的，在数学上也称第二类边界条件，或称Neumann边界条件，表示为：

式中，——表示边界E上压力关于边界外法线方向导数。

实际上，最简单、最常见的定流量边界是封闭边界，也叫不渗透边界，如尖灭或断层遮挡，即在此边界上无流量通过。

第三类边界条件为前两类的混合形式：

内边界条件：当有煤层气生产井时，由于井的半径与井间距离相比很小，所以可把它视作点汇当内边界来处理。在煤层气储层数值模拟中，可考虑两种内边界条件。

定产量条件：当给定井的产量时，可在微分方程中增加一个产量项。根据裘皮产量公式，煤层气井的气、水产量分别为：

式中：h为产层厚度；

P_wfg，P_wfw为分别为煤层气井的井底气、水流压；

r_e为排泄半径；

r_w为井筒半径；

s为表皮系数，已知压裂井的裂缝半长时，可用下式计算：

定井底流压P_wf条件

P_rw＝P_wf(x，y，z，t)

初始条件：给定在煤层气开发的初始时刻t＝0，煤储层内的压力分布和饱和度分布可表示为：

P_w(x，y，z，t＝0)＝P_wi(x，y，z)

S_w(x，y，z，t＝0)＝S_wi(x，y，z)

式中P_wi(x，y，z)和S_wi(x，y，z)是已知函数。

1.2.5煤层气储层模拟的数学模型

综上所述，描述煤储层中煤层气解吸、扩散、运移、产出的完整数学模型如下：

式中：

式中：

S_g+S_w＝1

P_cgw(S_w)＝P_g-P_w

ρ_g＝ρ_g(P_g)

ρ_w＝ρ_w(P_w)

K_rg＝K_rg(S_g)

K_rw＝K_rw(S_w)

μ_g＝μ_g(P_g)

μ_w＝μ_w(P_w)

P_cgw＝P_cgw(S_w)

P_w(x，y，z，t＝0)＝P_wi(x，y，z)

S_w(x，y，z，t＝0)＝S_wi(x，y，z)

P_rw＝P_wf(x，y，z，t)

P_g|_E1＝P_g1(x，y，z，t)

P_w|_E1＝P_w1(x，y，z，t)

1.3裂隙-孔隙双重孔隙结构煤储层气水两相流数值模型

上述建立的描述煤层气在煤储层中运移规律的数学模型是一个复杂的非线性偏微分方程(组)，无法用解析法直接求解。求解这类复杂的偏微分方程的通用方法是将方程及其定解条件离散化，然后采用数值法求解。

目前，在工程中应用的数值方法有限差分法、有限元法、变分法及有限边界元法。由于在油气藏数值模拟中，有限差分法应用最广泛，有关的理论和方法也比其它数值方法更趋成熟，所以本文采用有限差分法来建立描述煤储层内煤层气运移规律的数值模型，即差分方程组。

1.3.1离散差分

用有限差分法求偏微分方程的数值解，就是对连续问题进行空间和时间离散，用有限差商代替微商，得到在一系列离散空间网格或离散时间点上连续解的近似解。

对空间进行离散，网格划分有块中心网格和点中心网格两种格式。本文采用直角坐标下的块中心差分网格，即用网格把求解区域剖分成小块，用块的几何中心作为节点。假设节点的坐标以(i，j，k)来表示，则在x、y、z三个方向上其前、后、左、右、上、下邻块的中心坐标标号分别为(i-1，j，k)、(i+1，j，k)、(i，j-1，k)、(i，j+1，k)、(i，j，k-1)、(i，j，k+1)，这个块的前、后、左、右、上、下边界的坐标标号相应地分别为(i-1/2，j，k)、(i+1/2，j，k)、(i，j-1/2，k)、(i，j+1/2，k)、(i，j，k-1/2)、(i，j，k+1/2)。

对时间进行离散，就是将整个计算时间剖分成多个时间段，通过前后时间段间的数据传递进行计算。

1.3.2非线性差分方程组的建立

在不均匀网格条件下，采用块中心差分格式，对气、水相偏微分方程(1-23)和(1-24)的左端项进行空间差分，右端项进行时间差分。

最终得到描述煤储层中煤层气、水两相流体运移规律的差分方程组，即数值模型为：

上述方程组中包含的四个未知变量P_g、P_w、S_g和S_w，实际上只有两个是独立变量，其余变量可作为这两个独立变量的函数处理。我们选择气相压力P_g和S_w水相饱和度为独立变量进行求解，因为其对毛管压力曲线没有什么特殊的要求，可以完全忽略毛管力的影响。如果选择P_g、P_w为独立变量，要求毛管压力不能为零。

1.4数值解法研究

由于差分方程组(1-27)的左端项系数一传导率T_g和T_w是求解变量气相压力P_g和水相饱和度S_w的函数，因此该差分方程组是非线性的。对这种非线性差分方程组求解时，需要先采用一定的方法进行线性化，使之转化为线性的差分方程组，然后再用某种方法迭代地进行求解。在油气藏数值模拟中，线性化方法有显式方法、IMPES方法、隐式交替方法、半隐式方法和全隐式方法等，它们的区别主要是对传导系数及毛管压力的非线性问题在时间方面的处理方法不同。其中最简单、常用的方法是IMPES方法，即隐式压力显式饱和度方法。IMPES方法对差分方程组左端达西项的变量作隐式处理；对达西项的系数作显式处理，即对所有与时间有关的参数一律取用上一时间步tⁿ的值，对毛管压力P_c也同样取上一时间步的值，其采用的非线性差分方程组形式为：

虽然IMPES方法计算简单、工作量少、易理解，但是计算结果精度低，解的稳定性差、只适用于简单、时间步长较小的模拟计算，因此具有较大的局限性。

目前，在油气藏数值模拟中广泛采用的是全隐式方法，其优点是稳定性好，收敛快，时间步长大。因此，本次对差分方程组左端的达西项、解吸扩散项和井的产量项均作隐式处理，得出的全隐式非线性差分方程组形式如下式所示：

1.4.1全隐式解法的基本原理

全隐式方法求解的基本原理，基于数学上解非线性方程的牛顿法。用此方法求解方程时，每一时间步都要经过一系列迭代过程。其基本做法是：在第一个时间步开始时，先选择一组迭代初值。随后，在每一个n+1时间步的开始，先按第n时间步末所得的求解变量的值，求出方程组内各系数的值，接着解方程组，开始迭代。每一次迭代都求出求解变量的一组新值，以此求出新的系数来更新原来的近似系数值，然后再解方程组，这样逐次迭代下去，直至求出一组满足精度要求的值为止，这最后的一组迭代值便可作为n+1时间步的终值。然后再转入下一个时间步的迭代，这样一步一步地迭代下去。

对任意变量X，设从tⁿ到tⁿ⁺¹的Δt时间步长内，变量X的差值定义为：

δX＝Xⁿ⁺¹-Xⁿ

同时，定义任意两次迭代k和k+1之间，变量X的差值为：

当k＝0，X⁰＝Xⁿ。经过多次迭代，当|X^k+1-X^k|＜ε时

X^k+1＝Xⁿ⁺¹

在迭代过程中

因此，从tⁿ到tⁿ⁺¹的时间步的迭代过程中，方程组变成：

1.4.2全隐式线性化

分别对方程组的左端项和右端项进行全隐式线性化展开。虽然上文已述，选择的求解变量是P_g和S_w，但是为了求解方便，在实际解方程时不是直接求解k+1迭代步的和而是变量增量和

1.5全隐式线性差分方程组的求解

全隐式线性差分方程组的求解包括方程组的空间展开和参数取值和边界条件处理。边界条件处理包括参数取值、导数项处理、井参数处理、外边界条件处理。

1.5.1方程组的求解

按照上游权原则，将上述方程组中所有下标为的求解变量及系数归并到相应的 节点上，并且将每个节点上的两个求解变量按和的次序排列，经过整理，得：

方程组(1-79)中，系数a、b、c、d、e、f、g分别表示左节点、前节点、下节点、本节点、上节点、后节点、右节点的求解变量的系数，上角标1、2分别表示的系数，g、w分别表示气、水方程。上述方程组也可写出如下的矩阵形式：

在要求解的三维区域将方程组(1-79)开展，对M×N×S网格系统，共有M×N×S个节点，每个节点上都有两个方程，展开后共得到2×M×N×S阶方程组。写成矩阵形式有：

Au＝H

其中A是一个七对角分块矩阵，每个块都是二阶的子矩阵。采用线松弛迭代方法对该方程组进行求解。首先将三维中线z方向，即垂直方向选定为“线”方向。再固定x坐标，形成一块x-z板。然后在板上按照二维剖面模型的扫描方法进行扫描。完成后进行第二块板的扫描，这样依次将所有的板都扫描完，并反复进行，直到得到满足精度要求的解，从而完成方程的一次求解过程。为了加快收敛速度，要在每一条线上的解得出后，施加松弛因子算出该线上各节点未知量的该次迭代值。随后再进行下一条线的计算。

在求得k+1次的的值后，由以下公式计算出时的值，即：

由毛管理压力公式，水、气饱和度的关系式求出的值，即：

重复上述计算过程，直到足够小，满足精度要求为止，此时的值即为所求的值。再计算出本时段的气、水产量值。由于方程组的系数也参与迭代，在每次迭代求解后都被更新，因而全隐式解法的系数隐式程度高，解的稳定性好，结果精度高。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。