用于模拟地质力学储层系统的计算机实现的系统和方法.pdf

6.具体实施方式

6.1地质力学储层模拟方法的例子

图1描绘了用于模拟地质力学储层系统的计算机实现的系统的例子。该模拟系统包括描述地质力学储层系统的各个物理方面的多个模型，诸如储层流体流模型、地质力学模型、和裂缝流模型。储层流体流模型描述例如渗流(porous flow)、生产和注入。地质力学储层模型描述，例如，当将流体注入储层中或从储层中生产流体时以及当将应力施加于储层的边界时引起的应力、应变、和位移。裂缝流模型描述，例如，裂缝中的流和从裂缝到储层的泄漏。非线性偏微分方程组将这些模型的各个方面相互关联。

在接收到代表与要分析的地质力学储层系统相关联的物理性质的数据之后，解算器通过求解偏微分方程组生成预测(例如，压裂预测)。在图1的解算器中，偏微分方程组通过全展开雅可比行列式耦合。偏微分方程组的求解包括在单个时间步长中根据接收的物理性质数据同时求解全展开雅可比行列式。可以将生成的压裂预测输出给各种部件，譬如，输出给用户接口设备、计算机可读存储介质、监视器、用户可访问本地计算机、或作为网络一部分的用户可访问计算机。

非线性偏微分方程组包含与要用在分析地质力学储层系统中的模型相对应的方程。例如，以及如上所述，图1提供了非线性偏微分方程组包括与地质力学储层系统的储层流模型、地质力学模型、和裂缝模型相对应的方程的例子。作为另一种例示，图2提供了非线性偏微分方程组包括与地质力学储层系统的储层流模型、地质力学模型、裂缝模型、和热模型相对应的方程的例子。在图3中，非线性偏微分方程组包括与地质力学储层系统的储层流模型和地质力学模型相对应的方程。在图4中，非线性偏微分方程组包括与地质力学储层系统的储层流模型、地质力学模型、和热模型相对应的方程。在下文的6.3节中将提供与地质力学储层系统的不同模型的每一种相对应的方程的例子。

模拟地质力学储层系统的方法还可以在单个应用中组合水力裂缝生长、多相/多成分达西(Darcy)/非达西渗流、热传导和对流、固体沉积、和多孔弹性/多孔塑性形变。更进一步，可以组合像裂缝宽度变化、裂缝中的层状渠流(laminar channel flow)、储层中的渗流、热传导和对流、和多孔弹性/多孔塑性形变那样的不同机理的微分方程，来生成隐性全耦合公式。该模型可以包括导致弱固结地层中的裂缝生成的多孔塑性，并且能够模拟在裂缝与储层之间流动的复合多相流体注入。一种方法能够模拟牵涉到包括EOR和热处理的复相行为的多孔弹性/多孔塑性物质中的正在生长的裂缝。

模型各个方面的耦合例如可以通过全展开雅可比行列式中的变量来实现。例如，全展开雅可比行列式可以通过一个或多个如下变量将储层中的流体流与地质力学模型耦合：与地质力学模型相关联的有效应力、孔隙率、和一个或多个位移。将地质力学模型与流体流耦合的全展开雅可比行列式中的变量可以是与储层流模型相关联的孔隙率和渗透率。将地质力学模型与裂缝模型耦合的全展开雅可比行列式中的变量可以是与裂缝模型相关联的裂缝宽度。将储层中的流体流与裂缝模型耦合的全展开雅可比行列式中的变量可以是与裂缝模型相关联的流体压力和泄漏率。将热模型与地质力学模型耦合的全展开雅可比行列式中的变量可以是与热模型相关联的热应力。将热模型与储层流模型耦合的全展开雅可比行列式中的变量可以是与热模型相关联的储层中的流体粘度、传导和对流。全展开雅可比行列式可以包括与耦合变量的变化率(即，偏时间导数或时间导数)、空间导数、或偏空间导数有关的项，其中导数可以是任意阶的，例如，一阶导数、二阶导数、三阶导数等。耦合变量的一阶、二阶、三阶和/或更高阶导数(无论是时间导数还是空间导数)可以包括在全展开方程组中。在下文的6.4节中将提供可以耦合不同模型的变量的例子。

一种系统和方法可以通过将不同裂缝延伸模型(描述在下文的6.5节中)应用于地质力学储层系统来预测裂缝生长。例如，可以通过输入到裂缝模型的参数来预测裂缝的轮廓(6.5.1节)。在另一个例子中，根据不同裂缝延伸判据预测裂缝的演变轮廓。也就是说，可以根据临界应力强度因子(6.5.2节)，或根据呈现应变软化行为的内聚单元(6.5.3节)来模拟裂缝生长计算。可以耦合这两种裂缝延伸判据(6.5.4节)。在6.5.6和6.5.7节中将讨论裂缝生长预测的其它特征。

如6.5.5节所讨论，一种系统和方法可以配置成允许干旱带作为求解过程的自然部分出现在裂缝末端处。尤其，该系统和方法能够在裂缝末端处根据两种不同裂缝延伸判据模拟裂缝生长，并且还允许干旱带自然地出现在延伸裂缝的末端处。例如，当沿裂缝存在大压降时，在延伸裂缝的末端处出现干旱带。

全展开雅可比行列式的非线性方程组可以通过像在6.6节中更详细讨论的手段那样的数值手段来求解，其中非线性方程组例如使用所有解变量的完全牛顿-拉夫逊(Newton-Raphson)展开式来求解，这提高了解的稳定性，并且实现非线性叠代的二阶收敛速率。在6.7节中将讨论本文所公开的不同方法的装置和计算机程序实现的例子。

在另一个方面中，一种方法可以包括如下步骤：接收指示与地质力学储层系统相关联的物理性质的数据；定义包含多个三维单元的三维网格，其中每个三维单元包含多个节点；定义包含多个二维单元的二维网格，其中每个二维单元包含多个节点，并且每个二维单元与一个三维单元共享至少四个节点；以及通过求解模拟地质力学储层系统的偏微分方程组生成压裂预测。该偏微分方程组可以包含地质力学储层系统的储层流模型、地质力学模型、和裂缝模型。该偏微分方程组通过全展开雅可比行列式耦合，其中例如使用计算机对该偏微分方程组的求解包括在单个时间步长中根据接收的物理性质数据同时求解全展开雅可比行列式。在三维网格上计算储层流模型和地质力学模型，而在二维网格上计算裂缝模型。将生成的压裂预测输出给用户、用户接口设备、计算机可读存储介质、监视器、本地计算机、或作为网络一部分的计算机。图5例示了将裂缝延伸生长模型用于预测裂缝轮廓，并且还例示了干旱带被模拟。在6.2节中将描述用于本文的仿真方法的网格。

可以将分立网格用于模拟裂缝(以及例如储层等)中的流。例如，使用解变量的完全牛顿-拉夫逊展开式隐性地求解非线性方程组可以提高数值稳定性(例如，当应对裂缝生长、空腔生成、或牵涉到极小网格块的任何仿真时)。模拟与储层流分开的裂缝流使得可以顾及裂缝末端处的干旱带，而使用来自隐性耦合方程组的全展开雅可比行列式使求解过程更加稳定。

其它示例性手段显示在所附的美国临时申请第61/098,604号的附录A中。

6.2仿真方法

将分立的网格用于模拟裂缝中和其它系统中(例如储层中)的流。图6例示了用于计算裂缝流模型的二维(2D)网格602和用于计算其它模型的三维(3D)网格600的例子。例如，可以在3D网格上计算多点通量(例如针对储层和渗流)模型、地质力学模型(例如用于计算应力和位移)、和热模型中的一个或多个。二维(2D)网格用于计算沿裂缝的流(即，水力裂缝仿真)。用于裂缝计算和其它模型的分立网格可以使用并行处理手段来耦合网格。

3D网格可以是包含六面体单元的结构化或非结构化六面体网格。六面体网格单元具有八个角、十二个边缘(或边)、六个面。六面体网格单元每一个都可以包括至少八个节点(每个角上一个)，或更多个并且最多达二十七(27)个节点(即，每个面的中心上，每个边的中心上，每个边缘的中心上，和单元中心处的节点)。不同六面体单元可以包括不同数量的节点。

2D网格可以是包含四边形单元的结构化或非结构化四边形网格。每个四边形网格单元具有四个角和四个边缘。四边形网格单元每一个都包括至少四个节点(每个角上一个)，并且最多达五个节点(即，中心处的节点)。2D网格的节点可以配置成与3D网格的外边界之一上的节点重合，并且在公共节点上耦合诸如裂缝宽度和3D位移的某些计算。此外，2D网格上的计算可以是轴对称二维(2D)计算。

通过用2D网眼(mesh)表示单个2D裂缝表面，可以使用所公开的各种方法仿真裂缝，其中示使用3D坐标指定表面的位置。2D裂缝网眼可以是一般2D表面，无需位于单个平面内。2D裂缝网格的输入数据格式可以非常类似于3D储层网格的格式。通过指定网格单元的维度或坐标，可以将一个网格标识成裂缝网格。2D裂缝网格与3D储层网格之间的连接性通过将定义裂缝单元的四个节点与定义每个储层单元的每个表面的四个节点相比较来确定。图7例示了用在计算中的储层和裂缝网格的视图，其中距离以英尺为单位。

对于某些计算，将2D裂缝网格中的四边形裂缝单元与3D地质力学/渗流网格中的六面体单元的面对齐，以便使诸如流体流、位移、裂缝宽度、和牵引力的参数在分立的网格之间是一致的。作为考虑到裂缝网格与储层网格之间的相互作用的程序(用于实现本文所公开的方法)中的算法的结构的例子，这样的结构可以类似于在已经开发出来运行在具有分布式存储器的并行处理系统上的应用程序中实现的算法的结构。

可以使用例如几个2D网格在单次仿真中引入多个水力裂缝。2D裂缝网格可以被限制成沿着三维网格的边、顶部或底部。

6.3模型

下面提供与地质力学储层系统的每种不同模型相对应的微分方程的例子。组合包括在计算中的模型的微分方程以生成隐性全耦合公式。将一组一致的单位用于包括在计算中的方程中的所有变量。

6.3.1储层模型

用于渗流的方程组包括质量守恒：

∂∂t(ρφ)=-▿(ρU‾)+qw,......(1)]]>

其中，φ是孔隙率，而ρ是可以是压力的函数的流体密度。该模型允许在裂缝单元和储层单元两者中完成井，而上述方程中的q_w考虑到对储层单元的注入。例如，如图8所示，井可以位于储层和裂缝网格(显示在图7中)的边缘上。如果单位不一致，可以提供适当转换因子。

速度是相对于多孔物质的达西速度，可以通过下式来定义：

U→=-Kμ(▿p-ρg▿h),.......(2)]]>

其中，K是渗透率张量，μ是可以是压力的函数的粘度，p是流体压力，是重力项。

包括在流体流方程中的地质力学变量强调流模型与形变模型(在下面的6.3.2节中将描述一些地质力学项的定义)之间的耦合。

对于牵涉到温度变化的计算，可以以不同方式输入与温度相关的水性质。对于全耦合计算，可以作为压力(P)和温度(T)两者的函数输入水性质。对于叠代耦合计算，可以作为压力的函数输入水性质，然后将修正因子用于温度效应。在6.3.4节中将更详细地说明与温度相关的性质的处理。

流体的热行为也可以通过使用修正因子修正流体性质来模拟(描述在下面的6.3.4.3节中)。

6.3.1.1多相渗流

储层模型实现了从单相到基于重油到易挥发油的复合物的若干相行为模型。可以针对水相、非水液相、和非水气相和nc成分模拟达西流。任何相行为模型都可以与本文公开的渗流模型一起使用。流体流方程可以用一般复合公式来表示。代表多相流的成分质量平衡的偏微分方程是：

∂∂t(φNic)=-Σα▿xicαραU→α+qic---(3)]]>

其中，N_ic是由给出的单位孔隙体积的成分ic的浓度，是α相中成分ic的摩尔分数，ρ_α是α相的摩尔密度，而q_ic是单位储层体积的成分ic的摩尔流速(flow rate)。α相的速度通过下式给出：

Uα=-Kkrαμα·(▿Pα-γα▿D)---(4)]]>

相压通过下式定义：

P_α＝P+P_cα    (5)

其中，P_cα是毛细压力，和P是基准压力。基准压力用于PVT计算、井计算、和地质力学计算。基准压力对于二相模型是非水相压力，而对于三相模型是非水液相压力。

对于渗流模型，将孔隙率定义成：

φ＝φ_o[l+c_r(P-P_o)]    (6)

其中，φ₀、c_r、和初始压力P₀是位置的函数。

对于包括地质力学模型的计算，相对于初始未变形容积的孔隙率通过方程15给出，从中看出，c_r可以与孔隙率的毕奥(Biot)常数有关。

6.3.1.2非达西流

在某些例子中，可以将福施海默(Forschheimer)方程用于修正压力梯度与流体速度之间的关系来模拟非达西流。在其它例子中，通过指定流体速度与压力梯度之间的一般关系来模拟非达西流。

对于牵涉到各向异性介质中的3-D多相流的系统，非达西速度的福施海默方程(将取代上面的方程2)通过下式给出：

-krαμαK‾▿→pα=(1+βKραkrαμα||v→α||2)v→α---(7)]]>

其中，μ_α是α相的粘度，K是渗透率张量，k_rα是α相的相对渗透率，v_α是α相的达西速度，表达式||...||₂是通过定义的L₂-范数，而参数β_k是非达西系数。参数β_k可以遍及储层而变化，因此，可以针对每个网格块输入β_k的值。非达西系数β_k与过渡常数的倒数有关，即，β_k＝1/F。有关内容请参见Barree等人，″Beyond Beta Factors，A Complete Model for Darcy，Forschheimer，and Trans-Forschheimer Flow in Porous Media，″SPE 89325，SPE annual Technical Conference and Exhibition，Houston，TX，September 26-29，2004。α相的雷诺(Reynolds)数通过如下方程给出：

Nreα=βkραkrαμα||v→α||2---(8)]]>

方程8中的各项的单位应该选择成使结果是无量纲的。在组合了方程7和8之后，福施海默方程变成：

v→α=-krαμα(1+Nreα)-1K‾▿→pα---(9)]]>

其中，非达西流被表达成随该相的雷诺数而变的与相相关的、渗透率修正函数。这不同于标准渗透率修正函数，因为这个非达西公式对于每个相具有分立值。

在一些情况下，福施海默方程不提供非达西流的适当近似。如下形式的修正函数可以用于近似非达西流：

v→α=-krαμαfα(Nreα)K‾▿→pα---(9A)]]>

可以构建满足和f_α(0)＝1.0约束的修正函数。对于标准福施海默方程，可以指定如下函数：

fα(Nreα)=(1+Nreα)-1---(9B)]]>

对于广义福施海默方程，可以指定如下函数：

fα(Nreα)=Kratio+1-Kratio1+Nreα---(9C)]]>

J.L.Miskimins等人，″Non-Darcy Flow in Hydraulic Fractures：Does It Really Matter？″SPE paper 96389，SPE annual Technical Con-ference and Exhibition，Dallas，TX，October 9-12，2005中的方程3是可应用于本文所公开的方程的非达西公式的另一种形式。

6.3.1.3储层模型中的计算

储层模型的计算是在3D网格(用于地质力学模型的网格)上。3D储层网格可以包括渗流计算。针对储层单元之间的流以及针对储层单元与裂缝单元之间的流计算流体速度项。渗流方程的主要变量可以是可以在每个六面体单元(基于单元)的中心处评估的流体压力或流体组分。在某些计算中，将多点通量算法用于非结构化储层流计算，因此，当八个单元共享一个公共角时，对于3D网格的一般六面体单元，所得计算模板是27个点。

6.3.2地质力学模型

6.3.2.1多孔弹性物质

将线性关系(小应变)用于应变-位移关系。将应力、应变、温度、和孔隙压力相关联的耦合流/位移模型基于毕奥多孔弹性理论。平衡方程基于总应力，并假设准静态平衡。

多孔弹性方程用总应力、体应变、温度、和孔隙压力表述。总应力通过在储层的一个平坦部分上观察到的平均牵引力来定义，其中该平坦部分包括固体承担的负载和来自流体的孔隙压力。体应变是如果将应变计附连在变形多孔物质上，从应变计上观察到的应变。

线性多孔弹性位移的方程组包括应变-位移方程：

ϵij=12(ui,j+uj,i),---(10)]]>

其中，逗号指微分，u_i是沿着i方向的位移，ε_ij是多孔物质的体应变，展开式对应于正应变。总应力满足如下平移方程：

σ_ij，j+f_i＝0(11)

其中，应力张量是对称的，重力项f_i是固体密度、流体密度、和孔隙率的函数。可以在计算模型的三维网格的所有六个边界上沿着所有三个方向指定牵引力或位移边界条件。

当不考虑温差时，将总应力、体应变、和孔隙压力相关联的本构方程是：

σij=σijo+E1+v(ϵij+v1-2vϵkkδij)-α(p-po)δij---(12)]]>

其中，张力是正的，重复索引指求和，是初始原地应力，p₀是初始压力，E是弹性模量，v是泊松比，α是应力/应变方程中的毕奥常数，δ_ij当i＝j时是1，而当i≠j时是0。可以假设当时，应变是零。

在考虑温差的例子中，本构方程是：

σij=σijo+E1+v(ϵij+v1-2vϵkkδij)-α(p-po)δij-αTK(T-To)δij---(13)]]>

其中，α_T是应力/应变方程的热体积膨胀系数，K是体积弹性模量(elastic bulk modulus)。压力p₀是初始孔隙压力，温度T₀是初始温度。

如果将应力和压力项组合成那么，该方程变成应力已经被有效应力取代的标准热线性弹性本构方程。如果初始原地应力和初始孔隙压力是零时，那么，该方程呈现标准形式：

σije=E1+v(ϵij+v1-2vϵkkδij)-αTK(T-To)δij---(14)]]>

孔隙率(相对于未变形容积(undeformed bulk volume))与应变和流体压力之间的关系(当不考虑温差时)通过下式给出：

φ=φo+αϵii+1M(p-po)---(15)]]>

其中，方程15假设初始应变是零，φ₀是初始孔隙率，M^-1是孔隙率方程中的孔隙压力的毕奥常数。

当考虑温差和沉积分数时，相对于未变形容积的孔隙率定义成：

φ=φo+α(ϵv-ϵvo)+1M(P-Po)-αV(T-To)-σ---(16)]]>

其中，α和M^-1是毕奥常数，P是相压力(对于多相流)，α_V是孔隙率的热体积膨胀系数，σ是沉积分数(例如，计算中的网格单元的单位容积沉积的固体废料的体积分数)。固体废料随着废料穿过储层可以沉积在孔隙内，当废料沉积物堆积在孔隙空间中时，可以降低孔隙率和渗透率。计算网格单元中的孔隙率可以是流体压力、温度、和形变的函数，可以将由固体沉积引起的孔隙率减小量设置成等于σ。

对于各向同性物质，在应用地质力学方程模拟地质力学储层系统之前，确定六个多孔弹性物质参数：E，v，α，以及M^-1，α_T，和α_V。

在某些例子中，对于每个点上的渗透率和对于每个时间步长，可以将储层渗透率表达成初始定向渗透率(K_abs)乘以渗透率乘数f：

K＝K_abs×f

其中，f是诸如流体压力、总应力、容积应变(bulk volumetric strain)、孔隙压力、初始基准压力、主应力、有效塑性应变、当前孔隙率、初始孔隙率、和沉积分数的一个或多个其它参数的函数。

6.3.2.2多孔塑性物质

多孔塑性物质呈现非线性行为，因为它可能经历永久(即，塑性)体积应变，并因此经历孔隙率变化。大流体压力可以使多孔塑性物质屈服。其结果是，多孔塑性物质的地质力学计算可以预测到流体孔隙率的突然大变化。孔隙率的这些突然大变化可以引起严重的稳定性问题，并且还产生负流体压力。负压力一般出现在流体压缩率低，渗透率低，和孔隙率突然大增的时候。在上面的6.3.2.1节中讨论的多孔弹性物质的方程也可应用于多孔塑性物质。但是，要修正孔隙率，以便顾及可以针对多孔塑性物质预测的孔隙率的变化。

在某些例子中，使用一个方程来抑制孔隙率突然变化，以便提高计算的数值稳定性和降低遇到负压力的频率。在这些计算中，储层模型中的孔隙率被定义成流体孔隙率(φ_fluid)，并且被视作不同于地质力学模型中的孔隙率(φ_geomech)。在计算期间将通过如下阻尼方程来支配流体孔隙率与地质力学孔隙率之间的关系：

使用方程17计算流体孔隙率(φ_fluid)，并将它用在流体流方程中，而在地质力学模型中计算地质力学孔隙率(φ_geomech)，τ是规定时间常数。在(φ_geomech)阶跃变化之后，φ_fluid与φ_geomech之间的相对差值在5τ的时段之后小于1％。τ的值可以选择成使计算是稳定的，并且可以选择成使时间间隔尽可能短，例如，通过将τ设置成比计算中的时间步长的时帧或计算的总时间短的值。举例来说，如果预计计算持续几天，那么τ可以是分钟的数量级；如果预计计算持续几分钟，那么τ可以是毫秒的数量级。为了对许多多孔塑性物质进行计算，可以将τ的值设置成大约一分钟。在其它例子中，可以将τ的值设置成零。用于给定计算的τ值对于本领域的普通技术人员来说是显而易见的。

6.3.2.3地质力学模型中的计算

地质力学模型的计算是在3D网格(用于储层模型的网格)上。例如，标准有限元方法可以用于在每个单元内部的八个高斯(Gaussian)点上积分应力和在每个单元的中心上积分流体压力的地质力学方程。当八个单元共享一个公共角时，离散化为位移产生27点模板，在每个节点上存在三个位移。在某些计算中，位移是地质力学方程的主要变量，其中在六面体单元的角(基于节点)上评估位移。

6.3.3裂缝模型

裂缝中的流体流可以被模拟成平行板之间的渗流，或被模拟成平行板之间的层流(laminar flow)。对于裂缝中的单相渗流，单位长度的体积流速(volumetric flow rate)是：

q=-kfwμ(▿p-ρg▿h)---(18)]]>

其中，w是感兴趣点上的总裂缝宽度。裂缝宽度w是在裂缝网格中的所有节点位置上确定的。在裂缝模型的计算中，通过2-D裂缝单元的一个边的体积速率是通过将方程18乘以裂缝单元的该边的长度给出的(该边上分布着单位法线)。

在裂缝模型的计算中，裂缝单元是在2-D四边形上的，其中通过四边形的每个边和裂缝张开来计算流，或在每个单元的节点上计算宽度。在计算流体流部分期间，对四边形的每个边缘上的两个裂缝宽度求平均，以获得通过该面的流速。对于井来说也包括宽度计算，以顾及裂缝张开较大的生产率提高。

对于平行板之间的层流(渠流)，为裂缝中的流假设润滑近似，其中单位长度的体积流速是：

Q→=-w312μ(▿p-ρg▿h)---(19)]]>

并且可以使计算局限于牛顿流。裂缝中的质量守恒被写成：

∂∂t(ρw)=-▿(ρQ→)+qw-qL---(20)]]>

其中，q_w代表进入到裂缝单元中的井注入。项q_L代表从裂缝面到储层单元的泄漏，它被包括在用于储层中的质量守恒的流体速度项中。裂缝中的流体压力控制裂缝中的流速，并且用作包括地质力学计算的3D网格上的法向牵引力条件的作用。上面的方程适用于整个裂缝宽度，当只模拟一半裂缝宽度时，在适当项中应该包括1/2的倍数。

当使用全耦合手段时，可以在全展开雅可比行列式中包括方程(19)和(20)中的裂缝宽度的导数。当使用叠代耦合手段时，一般不包括这些项。在方程(19)和(20)中包括这些有关裂缝宽度的项可以提高叠代过程的数值稳定性。

6.3.3.1裂缝模型中的计算

几种裂缝网格可以用于在仿真时计算分立的裂缝。但是，在某些例子中，每个裂缝网格只与井和储层单元联系，并不与其它裂缝网格联系。裂缝中的流使用渠流近似(平行板之间的张开层流)或使用渗流方程(达西或非达西)来模拟。对于裂缝中的渗流，可以与为相邻储层单元指定的值无关地指定裂缝单元的诸如孔隙率、渗透率等的储层性质。也可以为裂缝单元指定岩石表和PVT表，这些表格不同于为相邻储层单元指定的表格。

在某些例子中，裂缝单元是四边形，其中每个四边形与3-D储层模型中的六面体单元的面对齐。如果裂缝单元在3-D多孔固体的外部，每个裂缝单元与零个3-D单元联系，如果裂缝单元在3-D多孔固体的边界上，每个裂缝单元与一个3-D单元联系，或如果裂缝单元在3-D多孔固体的内部，每个裂缝单元与两个3-D单元联系。如果两个3-D单元共享一个公共裂缝单元，则在3-D网格中关闭两个3-D单元之间的联系(流动和传导)。

在某些例子中，通过指定经过裂缝节点的井路径，将裂缝与注入或生产井连接。井可以在包含这些节点的裂缝和3-D单元两者中完成。

在某些例子中，将多点通量算法用于裂缝流计算，因此，当四个单元共享一个公共角时，对于一般四边形单元，所得计算模板是9个点。裂缝宽度从地质力学方程中的位移中确定。裂缝的流速和质量累积项对裂缝宽度w的变化高度敏感。对于渗流和渠流两者，方程19和20中的w的变化包括在全耦合裂缝/流/地质力学解或裂缝/流/地质力学/热解的雅可比行列式中。

多相流可以通过将相对渗透率函数用作上面方程18，19和20中的倍数来模拟。渠流的典型相对渗透率函数在饱和区中是线性的。在上面的体积速率方程18和19中也可以包括与压力相关的渗透率(对于渗流)和非达西流。

6.3.4热模型

热模型计算当将冷或热流体注入储层中或从储层中生产冷或热流体时发生的温度变化，并且还计算来自井筒的传导(可以是循环冷/热流体但不是实际注入或生产流体)。温度计算可以包括渗流和地质力学，但可以配置成不包括蒸汽注入。假设注入水处于液相。

热模型计算由储层中的传导和对流，井上的冷或热流体注入/生产，井上的传导，以及流体与固体之间的热和机械相互作用引起的温度变化。将这些机制组合在与渗流方程和固体形变方程一起求解的单个能量方程中。能量方程用多孔固体的拉格朗日(Lagrangian)手段表述，所有流体运动都是相对于多孔固体的运动。

对于这种类型的公式，当评估储层单元的能量平衡时多孔固体的质量是常数，而流体的质量随着流体流入和流出多孔固体而发生变化。

6.3.4.1组合累积项

流体和固体中的能量变化的组合表达式(当在计算中包括地质力学系统时)通过下式给出：

其中，是相对于未变形容积的孔隙率，α是相，ρ_α是相的密度，S_α是α相的饱和度，h_α是α相的比焓，α_T是应力/应变方程中的体积展开系数，K是渗透率，α_v是孔隙率方程中的体积展开系数，T₀是流体温度。传导和对流项可以影响组合累积项。在方程21中，作出了相压力相同的近似，用近似表示，而C_r项是热容C_v是在恒定孔隙体积和恒定容积下测量的多孔固体的热容。可以证明，如上定义的热容C_r是恒定流体压力和恒定容积下的热容。方程21是一个单元中内部能量的总变化率的一般公式，当在计算中考虑包括流体温度T₀的项并且计算包括热模型和地质力学模型两者时，用于累积项。

如果在计算中考虑流体温度T₀时，方程21呈现较简单形式。

对于包括热模型和地质力学模型，以及在计算中不考虑包括流体温度T₀和的项的计算，简化表达式是：

对于包括热模型，但不包括地质力学模型的计算，累积项是：

而对于包括热模型，但不包括地质力学模型的计算，以及流体内部能启用流体焓近似表示的计算，该表达式是：

6.3.4.2能量方程

上面对累积项的选择决定了包括在最终能量守恒方程中的方程。该能量方程被表达成：

其中，在注入井处应用恒温边界条件。由于方程25不包含除了输运之外的任何机械项，当焓只是温度的函数时，网格块中的温度不应该受流体相或固体的膨胀或收缩影响。但是，方程25的确仍然考虑到从液体到蒸汽的蒸发热，因为在一种成分在液相中的焓与它在汽相中的相应焓之间的差中可以包括潜热。

如果一种仿真包括地质力学计算，以及在计算中考虑包括流体温度T₀的项，那么，一般能量方程通过下式给出：

项出现在能量方程的两侧，即，在累积项中，也作为对整个固体所作的功，因此，在方程26中不包括

6.3.4.3对流体密度和粘度的修正

存在几个计算流体的与温度相关的密度和粘度的选项，可用于计算密度和粘度的方法随相行为模型和正在使用的数值技术而变。在某些例子中，作为温度和压力的函数直接计算流体性质。在某些其它例子中，作为当前单元压力和初始单元温度(等温急流)的函数计算流体性质，并且将修正函数用于校正初始单元温差与当前单元温度之差。

流体性质可以对于隐性单相运行、隐性双相运行、和隐性复合运行作为温度和压力的函数直接计算。这些牵涉到地质力学模型的计算可以叠代耦合或全耦合。

可以修正在重油、K-值、和复合PVT包中计算的粘度和密度，以顾及由温度变化引起的流体性质的变化。在计算期间初始储层温度和当前流体压力可以用于在每个PVT包(等温急流)中计算流体性质，然后，通过乘以用户提供的修正因子，即，来修正这些性质，其中是作为表格输入的，而μ(p，T₀)是在PVT包中计算的。

对于许多热-地质力学研究，兴趣主要集中在温度如何影响流体粘度和热应力上。对于这些类型的应用，可以使用流体性质的修正因子或直接计算，使用两种技术的结果非常类似。但是，对于温度严重影响密度和/或相分离的应用，可以作为温度和压力的函数使用流体性质的直接计算。

6.4各种模型的相互关系

一种系统和方法可以耦合不同地质力学储层系统模型之间(例如，渗流，地质力学位移和应力，和裂缝的流体流和生长之间)的相互作用。例如，地质力学解可以通过孔隙率和渗透率项影响渗流计算，并且可以通过裂缝宽度和裂缝延伸判据影响裂缝计算。举另一个例子来说，储层中的流体解可以通过它在有效应力中的作用影响地质力学计算，并且可以通过调节裂缝与储层之间的流体流来影响裂缝计算。在又一个例子中，裂缝中的流体解可以通过出现在裂缝的面上的法向牵引力影响地质力学计算，并且可以通过泄漏影响储层中的流。

6.5裂缝延伸模型

如下面的6.5.1节所讨论，一种系统和方法可以实现预测裂缝生长，以便调查在网格的边界上水平或垂直裂缝如何随时间演变。在本例中，裂缝的轮廓和裂缝宽度通过输入裂缝模型中的参数来预测。

在裂缝中的流与裂缝宽度之间存在强相互作用，并且可以预测在裂缝前沿裂缝宽度很小。裂缝生长使用地质力学模型来预测，其中使用流体流的隐含时间步进的全耦合地质力学模型。在这些例子中，如下面的6.5.2和6.5.3节所讨论，根据应力强度因子或内聚单元预测裂缝的演变轮廓(即，裂缝延伸)。可以预测沿着3-D储层网格的边缘的裂缝延伸。两个裂缝延伸判据允许确定裂缝在何时以及在哪里生长。裂缝延伸判据可以用作岩石的裂缝韧性的度量。这两个裂缝延伸判据均使用裂缝前沿附近的流体压力、应力、和位移状态来确定裂缝何时延伸。通过释放裂缝末端处以前具有零位移约束的节点来模拟延伸。当以前张开的裂缝节点上的位移变成负值时，假设在该点上裂缝已经闭合，并重新施加位移约束。闭合节点具有零强度，并且如果该节点受到作用在该节点上的张力，则以后可以重新张开。6.5.4节将描述可能出现在这两个判据之间的关系互动的类型。

6.5.5节将描述出现在裂缝末端处的干旱带。在6.5.6和6.5.7节中将讨论裂缝生长预测的其它特征。

6.5.1规定裂缝轮廓

裂缝网格中的裂缝的宽度轮廓可以通过输入裂缝模型中的参数来规定。例如，通过首先指定裂缝网格的布局和性质，然后在裂缝网格中的每个节点上指定裂缝宽度，可以将2-D裂缝输入到3-D储层网格。可以在瞬时输入数据中改变裂缝网格中的裂缝宽度的值，以仿真正在生长的裂缝的影响。

图9示出了已经插入储层网格中的两条裂缝。对于这个网格，为每个裂缝网格中的每个节点输入总裂缝宽度。但是，应该为位于网格边缘上的裂缝输入裂缝半宽度。图中的裂缝实际上位于储层网格的内部，但为了绘制而偏离了储层网格。图9示出了注入井处的裂缝和生产井处的另一个裂缝。对于这种情况，每个裂缝单元与两个储层单元联系。被裂缝单元分开的储层单元之间的联系被设置成零。这牵涉到对于流体流将传输率(transmissibility)设置成零，以及对于热计算将传导率设置成零。

对于规定的裂缝轮廓，裂缝网格中的每个节点的默认裂缝宽度被设置成用户指定值。可以分别为网格中的每个裂缝节点规定裂缝宽度。也就是说，裂缝宽度可以在初始数据中以及在计算期间的各个点上(例如，在输入数据的瞬时数据部分中)指定，或者只包含在初始数据部分中(其中作为时间和位置的函数地指定裂缝宽度)。

将一组流计算用于裂缝网格中裂缝宽度大于预定阈值的单元，而将另一组流计算用于裂缝网格中裂缝宽度小于预定阈值的那些单元。裂缝宽度小于预定阈值的具有所有角节点的单元被认为是闭合裂缝单元。对于这些闭合裂缝单元，单元间的流体流被模拟成渗流，其中从相邻储层单元获取单元中的渗透率并将位移约束应用在适当节点上。

对于裂缝宽度大于预定阈值的单元(即，张开的裂缝单元)之间的流，使用渠流模型或渗流模型。对于渗流模型，将输入裂缝渗透率用于流计算，其中将所有张开裂缝单元的裂缝渗透率设置成用户指定的默认值，或者用户可以为各个单元指定裂缝渗透率的值。另外，对于这种渗流模型，将每个裂缝单元中的孔隙率和孔隙压缩率设置成指定值。对于渠流模型，为裂缝网格中的所有单元将孔隙率设置成100％，并且，对于这种情况，没有孔隙压缩率项。

6.5.2应力强度因子

第一延伸选项使用基于应力强度因子的临界裂缝张开判据。临界应力强度因子基于裂缝末端处的渐近应力/应变状态。在某些例子中，这种延伸判据局限于线性多孔弹性应用。

对于在模式I加载条件下裂缝末端处的线性弹性、平面应变条件，末端附近的裂缝张开可以用模式I应力强度因子K_I表达成：

w(Δa)=4(1-v2)EKI2Δaπ---(27)]]>

其中，Δa是沿着与裂缝的前缘垂直的方向相对于裂缝末端的距离，w是总裂缝宽度，K_I是几何形状和加载条件的函数。每当裂缝末端处的流体压力大到足以使新裂缝区张开，以及相邻节点处的裂缝宽度超过如下宽度时，K_IC裂缝延伸选项就释放裂缝末端处的节点：

wneighbor≥4(1-v2)EKIC2Δaneighborπ---(28)]]>

其中，Δ_aneighbor是从相邻节点到裂缝末端处的节点的距离。图10例示了相对于裂缝末端的临界裂缝张开判据。在本例中，在采用这种裂缝延伸模型的计算中只检查与裂缝末端处的节点共享四边形的一个公共边的相邻节点。每当K_IC被设置成零时，该计算就模拟零韧度的裂缝。

6.5.3内聚单元

第二裂缝延伸判据使用四边形内聚单元来模拟裂缝末端处的内聚带。内聚模型由Hillerborg等人用于仿真裂纹形成(Hillerborg等人，1976，Cement and Concrete Research 6：773-782)。基于混凝土中的生长和有限元的内聚单元已经用在金属的弹性-塑性裂缝的研究中(Ortiz等人，1999，Int.J.Num.Methods in Engr.，44：1267-1282)。尤其，该内聚单元基于能量释放率和拉伸应力。将这种裂缝延伸模型用于多孔弹性和多孔塑性应用。

在一些例子中，假设裂缝流体可以流到内聚带的末端，以及在内聚带中存在试图将裂缝面结合在一起的张力。内聚单元可以用于模拟应变定域，也可以用于多孔塑性计算。内聚单元的另一个优点是它们可以用于裂缝产生和生长这两者，而K_IC手段要求裂缝已存在，或将某种其它逻辑用于裂缝产生。

在延伸水力裂缝末端处的计算中，内聚单元包括拉伸强度和能量释放率。这种裂缝延伸模型使用物质的拉伸强度和能量释放率来确定裂缝延伸的时间。假设裂缝平面上的物质服从如下形式的软化方程：

σ=σc(1-wδc)---(29)]]>

其中，σ是与裂缝平面垂直的应力，σ_c是拉伸强度，w是总裂缝宽度，δ_c是临界间隔。上面的方程假设w正增大。

σ_c、δ_c、和能量释放率之间的关系显示在图11中。每当w减小时，应力/位移曲线沿着不同路径。在某些例子中，当w减小时，该计算使应力/位移曲线向原点返回。图12例示了当节点的裂缝宽度从0增大到δ_a(路径1)并减小回到0(路径2)，然后增大回到δ_a(路径3)，最终生长到δ_c上的总间隔(路径4)时，内聚单元所遵循的加载/卸载路径。

当与裂缝平面垂直的应力等于拉伸强度σ_c时，在裂缝平面上开始分离。然后，界面上的物质随着物质进一步分离和与裂缝平面垂直的应力减小而开始变弱。一旦裂缝间隔超过δ_c的临界值，拉伸应力就变成零。δ_c的临界值被确定并且是临界能量释放率满足如下方程的值：

GC=12σcδc.........(30)]]>

图13例示了相对于裂缝末端的内聚单元裂缝延伸判据。

6.5.4应力强度因子和内聚单元的关系互动

对于线性弹性、平面应变、模式I问题，能量释放率通过下式与应力强度因子相关联：

GC=1-v2EKIC2........(31)]]>

临界应力延伸判据可以通过使用具有有限拉伸应力值和极小能量释放率值的内聚选项来实现，但观察到与网格高度相关的解，除非使干旱带出现在裂缝的末端处，或除非使用包括内聚带的影响的非常精细的网格。在临界应力判据的严格实现中，将在裂缝前面的固定距离处评估应力，以便随着网格的精细化保证算法的收敛。通过对平面应变、模式I裂缝的裂纹宽度使用渐近解，可以获得内聚带尺寸的粗略估计：

w(r)=8(1-v2)KICEr2π,.......(32)]]>

其中，r是相对于末端的距离。如果将方程32与方程30和31结合，那么，可以通过如下方程估计内聚带的长度：

rc=π8(KICσc)2,.......(33)]]>

这是与Dugdale-Barenblatt模型(在裂缝前面具有恒定拉伸应力的带屈服模型)中的塑性带相同的尺寸。方程33是内聚带尺寸的粗略估计，因为它基于不包含内聚带的裂缝的渐近解。另外，如果渐近解在感兴趣范围上无效，则该估计可能不准确。仿真结果可能表明，应该通过一个不同的方程估计内聚带的尺寸。例如，如结合特定例子的6.1.4节所示，内聚带的尺寸可以通过两倍于方程33的方程(π/4)(K_IC/σ_c)²来估计。

6.5.5裂缝末端处的干旱带

裂缝可能比裂缝流体可以充满在延伸裂缝的末端处产生的新裂缝单元更快地延伸。如果裂缝流体必须完全充满末端，这可以引起在延伸裂缝的末端处出现负压力。因为裂缝跑在粘性流体前面，所以在物理应用中，可能在裂缝末端处出现干旱带，或者，由于用于裂缝延伸的数值技术，也可能引起干旱带。干旱带作为本文所公开的方法的人工产物而出现，例如，当在一个时间步长期间释放节点，以及仿真中的时间步长太短使流体难以充满新产生的裂缝单元时。取决于裂缝流体粘度和裂缝宽度，裂缝末端处的干旱带对于一些应用来说可能是必需的，有关内容请参见Yew等人，1993，SPE 22875，SPE Production &Facilities，8：171-177。

为了使干旱带出现，一种系统和方法可以允许裂缝中的单元被流体部分饱和。部分饱和单元的示意图显示在图14中，其中流体从左往右传播而裂缝宽度从左往右变窄。裂缝流体是图14中的阴影区，虚线代表每个裂缝单元之间的边界。单元3被流体完全饱和，而单元2部分饱和，以及单元1没有流体。在一个例子中，允许裂缝流体从图14中的单元3流到单元2，但在单元2充满流体之前，不允许裂缝流体从单元2流到单元1。尽管在本例中不允许裂缝流体直接从单元2流到单元1，但如果储层具有非零渗透率，则储层流体仍然可以流入单元1中。在这个特定例子中的计算期间，对于单元1和2，单元中的流体体积是变量，而对于单元3，流体压力是变量。

部分饱和单元中的流体压力被规定为具有用户指定值，直到单元完全充满流体。规定压力的典型值是初始流体压力和零压力。但是，干旱带中的压力可以选择成小于周围储层单元中的压力的某个值。部分饱和单元被包括在质量守恒方程中。举一个例子来说，这些单元的主要变量是单元中的流体体积，而不是流体压力。通过让附加流体流入单元中，或通过让裂缝宽度减小到裂缝体积等于流体体积的点，可以使单元变成完全饱和。在仿真期间，干旱带可以作为求解过程的自然部分出现，并且，除了用于干旱带中的流体压力的值之外，该方法的实现不需要对干旱带作任何预定的指定。

6.5.6裂缝中的地质力学裂缝生长和渗流

在仍然模拟裂缝中的渗流的同时，对地质力学模型的计算可以被用于确定裂缝轮廓和裂缝宽度。这个选项用于位于网格边缘的裂缝，并且在运算期间输入的裂缝网格的所有裂缝宽度是裂缝半宽度。

初始裂缝半宽度的值可以被输入到输入数据集的初始数据部分中，其中为井附近的单元输入裂缝半宽度的非零值，以便产生以后随时间生长的初始小裂缝。在运算期间可以将应力强度因子或能量释放率用于计算演变裂缝轮廓。

此外，可以根据对地质力学模型的计算来计算裂缝半宽度。如果这些地质力学裂缝半宽度大于为初始裂缝半宽度指定的值，将它们用于流体流和孔隙体积计算。否则，将裂缝半宽度的输入值用于张开裂缝单元中的流和孔隙体积计算。

在本例中，闭合裂缝单元被识别为所有角半宽度小于裂缝宽度的预定阈值的那些单元。对于闭合裂缝单元，在计算时使用渗流方程，但当计算单元的容积和单元之间的传输率时，将位移约束应用在适当节点上。另外，将相邻储层单元中的储层渗透率用于裂缝网格中的闭合裂缝单元。

6.5.7裂缝中的地质力学裂缝生长和渠流

可以利用渠流进行对地质力学模型的计算，以预测裂缝轮廓。对于这个选项，张开裂缝单元中的流体流被模拟成平行板之间的层流，并且作了闭合裂缝单元无流零裂缝宽度的假设。在这个模型中未将渗透率或孔隙率用在裂缝网格中。这个选项还使半干旱带出现在裂缝的前缘附近。这个选项用于位于网格边缘的裂缝，并且在运算期间输入的裂缝网格的所有裂缝宽度是裂缝半宽度。

初始裂缝半宽度的值可以输入到输入数据集的初始数据部分中，其中为井附近的单元输入裂缝半宽度的非零值，以便建立以后随时间生长的初始小裂缝。在运算期间可以将应力强度因子或能量释放率用于计算演变裂缝轮廓。

此外，可以根据对地质力学模型的计算来计算裂缝半宽度。如果这些地质力学裂缝半宽度大于为初始裂缝半宽度指定的值，将它们用于流体流和孔隙体积计算。否则，将裂缝半宽度的输入值用于张开裂缝单元中的流和容积计算。

在本例中，闭合裂缝单元被识别为所有角半宽度小于裂缝宽度的预定阈值的那些单元。对闭合裂缝单元不进行计算。也就是说，假设闭合裂缝单元具有零容积和并且没有流。

6.6压裂预测方法的实现

完全展开包括在计算中的模型的整个方程组的雅可比行列式，并在线性解算器中同时求解所有方程。各种机制使用单向，显性，叠代，或全耦合技术来耦合。更费力的问题(例如，形成砂石，空腔完成，和水力压裂)使用最稳定耦合技术，并且在本文描述的整个水力压裂工作中使用全耦合求解技术。

在前面部分中讨论的方程可以在单个程序中实现，该程序计算所有方程的隐式解，也就是说，可以将反向欧拉(backward-Euler)技术用于近似瞬时(即，时间)导数，并且在时间步长结束时评估方程中的所有主要变量和系数。主要变量的非限制性例子是储层中的流体压力、储层中和储层边界上(如果未受约束)的位移、裂缝中的流体压力(或如果部分饱和，裂缝中的流体体积)、和井压力。方程组可以使用牛顿-拉夫逊技术来求解，其中为整个方程组生成全展开雅可比行列式，并且使用包括所有解变量的线性解算器求出修正增量。

通过计算穿过现有2D裂缝的流，并且与计算在相邻3D储层中出现的有效应力和位移一起计算到相邻3D储层的泄漏，可以仿真裂缝延伸。然后，将裂缝延伸判据与裂缝前沿上的应力和位移状态结合，以确定裂缝前沿延伸的时间。通过释放裂缝末端处以前具有零位移约束的节点来模拟延伸。当以前张开的裂缝节点上的位移变成负值时，程序假设在那个点上裂缝已经闭合，并重新施加位移约束。如果该闭合节点经受到张力，则该节点以后可以重新张开。可以对处在3D网格的边缘上的裂缝进行裂缝延伸预测。

本文讨论过的某些例子使裂缝可以在单个时间步长期间跨过多个网格块延伸。每个时间步长从来自前一时间步长的裂缝配置、流体状态、和地质力学状态开始。然后，对假设裂缝不延伸的质量守恒和应力平衡方程进行叠代计算以获得收敛解。该计算允许裂缝宽度在这组叠代期间发生变化，但不允许裂缝前沿移动。在计算达到收敛之后，将裂缝延伸判据用于裂缝网格中的所有节点，看闭合节点是否应该张开或张开节点是否应该闭合。如果裂缝前沿发生变化，那么重复对质量守恒方程和应力平衡方程的叠代计算，直到获得新收敛解。一旦获得新收敛解，就再次检查延伸判据，看裂缝前沿是否存在进一步的变化。可以对方程继续进行这一系列叠代并检查延伸判据许多次，或直到在移动到下一个时间步长之前对于一个时间步长再也看不到进一步的延伸为止。

在一个时间步长期间裂缝可以生长或收缩。全展开雅可比行列式的解可以使用考虑这种变化的三步过程来计算：

步骤1，第一步骤假设裂缝边界是固定的，并针对给定裂缝配置求解方程组(裂缝可以张开得更宽或变得更窄，但裂缝的边缘不发生变化)。该方程组使用牛顿-拉夫逊技术来求解，其中为整个方程组生成雅可比行列式，并且使用包括所有解变量的线性解算器求出修正增量。该程序生成整个组的雅可比行列式，以便提高健壮性，和保证非线性叠代二阶收敛。

例如，当使用全耦合手段时，可以在全展开方程组中包括方程(19)和(20)中的裂缝宽度的导数。在使用叠代耦合手段的传统方法中一般不包括这些项。包括这些方程(19)和(20)中的裂缝宽度的导数可以提高叠代过程的数值稳定性。

步骤2，在针对给定裂缝几何形状使解收敛之后，该程序接着使用上述的延伸判据之一确定裂缝边缘是否移动，并检查位移解，以确定裂缝表面的一些部分是否闭合。如果在边界上未检测到变化，那么，该程序继续前进到下一个时间步长；否则，该程序转到步骤3，然后对新裂缝边界重新从步骤1开始。这个过程继续直到裂缝边界不发生另外的变化。

步骤3，如果超过了延伸判据以及裂缝在生长，那么，在边界上的适当节点处撤消位移约束，并允许裂缝在那些点处张开，使裂缝中的流体压力变成新产生裂缝表面上的牵引力条件。如果裂缝收缩(即，表面位移预测裂缝宽度是负的)，那么，在适当节点上重新施加位移约束。

图5例示了这种手段的实现，其中裂缝延伸生长模型用于预测裂缝轮廓，步骤1，步骤2和步骤3用在压裂预测方法的实现中，并模拟干旱带。

这些方法的实现可能限制了可以在每组非线性叠代结束时出现的裂纹延伸的数量，因此，如果在单个时间步长期间裂缝跨过大量单元延伸，则该时间步长可以采用几个非线性叠代。许多CPU计算时间可能花费在每次非线性叠代期间生成的线性方程组的求解上。

6.7装置和计算机程序实现的例子

本文所公开的方法可以按照如下程序和方法，使用装置，例如，像本节所描述的计算机系统那样的计算机系统来实现。这样的计算机系统还可以存储和操纵指示与地质力学储层系统相关联的物理性质的数据、包括在计算中的模型的整个方程组的全展开雅可比行列式、全展开雅可比行列式的解、生成的压裂预测、和可以供利用本文所述的分析方法实现的计算机系统使用的测量结果。该系统和方法可以实现在各种类型的计算机架构上，例如，在单台通用计算机，并行处理计算机系统，或工作站上，或在联网系统(例如，像显示在图15中那样的客户机-服务器配置)上。

如图15所示，实现本文所公开的一种或多种方法和系统的模拟计算机系统可以与网络链路链接，网络链路可以是，例如，与其它本地计算机系统连接的局域网(“LAN”)的一部分，和/或与其它远程计算机系统连接的像因特网那样的广域网(“WAN”)的一部分。该系统和方法可以包括用于与一个或多个数据处理设备通信、经由网络(例如，局域网、广域网、因特网、它们的组合等)、光纤媒体、载波、无线网络等传送的数据信号。

模拟系统包含本文所述的任何方法。例如，软件组件可以包括使一个或多个处理器实现如下步骤的程序：接受指示与地质力学储层系统相关联的物理性质的多个参数、和/或生成的压裂预测；以及将指示与地质力学储层系统相关联的物理性质的参数、和/或生成的压裂预测存储在存储器中。例如，该系统可以接受用户人工输入的(例如，通过用户界面)、接收指示与地质力学储层系统相关联的物理性质的参数、和/或所生成压裂预测的参数的命令。该程序可以使该系统从数据存储装置(例如，数据库)中检索指示与地质力学储层系统相关联的物理性质的参数、和/或所生成压裂预测的参数。这样的数据存储装置可以存储在大型存储设备(例如，硬盘驱动器)或其它计算机可读介质上，并且被装载到计算机的存储器中，或者，该数据存储装置可以由计算机系统通过网络访问。

本公开的一个实施例提供了存储执行本文公开的任何方法的步骤的可由计算机执行的计算机程序的计算机可读介质。可以提供与具有一个或多个存储器单元和一个或多个处理器单元的计算机结合在一起使用的计算机程序产品，该计算机程序产品包含具有编码在上面的计算机程序机制的计算机可读存储介质，其中计算机程序机制可以装载到计算机的一个或多个存储单元中，使计算机的一个或多个处理单元执行显示在图5中的方块图的过程的各种步骤。

7.数值结果

在本节中，将讨论处理结果的例子。

7.1裂缝预测结果

这些例子展示了程序的不同方面，并例示了程序对分析解的确认。尤其，这些例子展示了模拟本文所述的方法的地质力学储层系统的预测能力。这些例子还强调了应用两个裂缝延伸判据时的差异，并且例示了预计干旱带可以出现在裂缝末端处的时间。这些例子还允许将使用地质力学储层系统模拟方法而获得的压裂预测与分析解相比较，以确认使用的裂缝延伸算法。

该系统和方法还能够模拟在裂缝与储层之间流动的复合多相流体注入。但是，在本节中给出的数值结果假设单个裂缝中的稍可压缩单相流未泄漏到储层中，以便于与分析解比较，以及将重点集中在程序的裂缝延伸和干旱带方面。在本节中给出的数值结果局限于稍可压缩、等温、单相流体。

用在下面7.1.1，7.1.2和7.1.3节中的流体是稍可压缩的，其密度可以写成ρ_o(1+c_w(P-P_r))，其中压缩率c_w是3.0E-6psi^-1，和基准强度P_r是14.7psi。7.1.1节中的结果是关于恒速注入流体粘度极低的平面应变裂缝中的。对于这个问题，通过应力强度因子，或通过能量释放率和裂缝末端处的拉伸应力来控制裂缝生长。在7.1.1节中的计算中应用两种裂缝延伸模型，以展示两种延伸模型之间的差异。7.1.2节中的结果是关于类似于7.1.1节的计算的，但使用高粘度流体。7.1.2节表明，当在裂缝中存在明显压降时，在裂缝的末端处出现干旱带。7.1.3节中的结果是关于展示程序模拟裂缝平面中沿着两个方向的生长的能力的四分之一币状裂缝的。7.1.3节中的计算还使用了低粘度流体，因此，可以将这些结果与分析结果相比较。在7.1.1，7.1.2和7.1.3节中，将分别等于1.0和0.0的α和M^-1的值用在地质力学模型中。

7.1.1具有低粘度流体的平面应变裂缝

这个例子使用基于应力强度因子的裂缝延伸判据和基于内聚模型的裂缝判据来审视裂缝生长预测。这个例子被设置成具有保证裂缝预测对裂缝延伸判据非常敏感的储层性质。注入井位于沿着x，y和z方向分别是1000ft x 1000ft x 20ft的正方形图案的一角中。该问题在z方向采用平面应变条件，因此，沿着z方向的位移是零。还假设裂缝沿着x方向延伸，以及沿着y方向的最小原地应力是4000psi。储层中的孔隙率是零，所有注入流体都保留在裂缝中。网格具有每一个在裂缝方向长度为1ft的55个网格块，在55ft以外间隔开始变大。仿真只包括裂缝的一翼，并且只模拟总裂缝宽度的一半。裂缝附近的储层网格与裂缝网格一起显示在图7中。裂缝网格被离散成55个一英尺网格块，以便与储层网格相容。以25bbl/min的表面速率(总井速率是100bbl/min)注入低粘度流体(0.1cp)到达10分钟。图16示出了注入压力与裂缝长度之间的关系的曲线图。

储层的弹性模量和泊松比分别是2.0E6psi和0.35，临界应力强度因子是1500psi√in。假设沿裂缝的压降是零以及流体充满整个裂缝，那么，裂缝长度可以表达成：

L=(QEtπ(1-v2)KICΔz)2/3,.......(34)]]>

其中，上面的Q是裂缝一翼的总储层速率(50bbl/d)，而Δz是裂缝高度(20ft)。方程34中的假设是均匀受压裂缝的应力强度因子通过给出，其中Δp是过压力；裂缝是椭圆形的，注入体积与裂缝一翼的裂缝体积之间的关系通过Qt＝πw₀LΔz/4给出，其中w₀是井筒上的总裂缝宽度，而L是半宽度；以及井筒上的裂缝宽度通过4(1-v²)LΔp/E给出。在10分钟注入之后，方程34预测裂缝应该在43.8ft上，即，裂缝应该跨过裂缝和储层网格中的约44个网格块延伸。对裂缝宽度和流体压力的类似分析方程预测井筒上的裂缝宽度和裂缝中的流体压力在10分钟之后，应该分别是约0.002836ft和4037psi。注意，在上面的方程中，由于应力强度因子变为零，所以裂缝长度是无界的。

使用1500psi√in的应力强度因子，以及使用能量释放率为0.987188psi-in的内聚模型来仿真该系统。内聚模型还需要拉伸应力值，并且使用100，300，500，和1000psi的拉伸应力值进行四次运算。由于内聚带通过裂缝末端处的内聚单元来模拟，所以不能期望使用内聚模型获得裂缝延伸的准确预测，除非网格精细到足以求解裂缝末端处的内聚带。关于这里使用的拉伸应力的值，对在方程33中给出的r_c使用粗略近似获得的内聚带的预测值分别是7.4，0.8，0.3，和0.07。图17示出了应力强度因子的预测结果、和具有不同拉伸应力值的内聚模型。

对于内聚运算，图17假设裂缝末端处在内聚带的末端处。对于与分析方程中的假设一致的仿真，沿裂缝的压降可忽略不计；例如，对于使用300psi拉伸应力的内聚运算，在除了裂缝末端处流体压力是4037psi的单元之外的所有单元中，压力都是4038psi。这些压力与4037psi的分析流体压力十分相当。另外，拉伸应力为300psi的内聚情况预测有利地与0.002836ft的分析值相当的0.002870ft的总裂缝宽度。

对于使用应力强度因子的仿真和对于使用300和500psi拉伸应力的内聚运算，裂缝长度的数值预测相当好地与分析结果一致。但是，对于使用100和1000psi拉伸应力的内聚情况，与分析结果存在显著偏差。对于100psi的拉伸情况，模型预测内聚带长度是15ft和裂缝长度为51ft。15ft的内聚带长度是以前获得的内聚带的粗略估计的两倍。对于这种运算，内聚带的长度与裂缝长度相比不算小，因此，不能期望拉伸应力为100psi的内聚模型与分析结果一致。这是因为分析结果假设任何塑性带(或内聚带)相对于裂缝尺度都是小的，拉伸应力为100psi的内聚模型不是这种情况。对于拉伸应力为300psi的内聚模型，仿真结果预测内聚带长度为约2ft。并且，这大约是以前获得的粗略估计的尺寸的两倍。使用0.25ft的网格块尺寸对拉伸应力为300psi的内聚模型进行第二次仿真。对于这种更精细的运算，预测裂缝长度是44.5ft(取代在较粗运算中看到的45.0ft)，而内聚带的尺寸是1.75ft(即，在末端处存在7个内聚单元)。此外，对于更精细的运算，井筒上的预测裂缝宽度是与0.002836ft的分析值十分相当的0.002834ft。

对于使用1000psi拉伸应力的内聚运算，粗网格(1ft单元)仿真预测最终裂缝长度为32ft。这与在10分钟注入之后预测裂缝长度为43.8ft的分析方程不太相当。如果假设内聚带的尺寸的早先估计偏离一半，那么，这种情况的内聚带尺寸应该是近似0.15ft。这么小的内聚带不能在1ft的网格间隔上求解。对于这种运算的1ft的网格间隔，内聚模型仅仅变成拉伸应力延伸选项(由于不求解内聚带)，在网格变得比内聚带的尺寸更精细之前，预计会看到强网格相关性。由于这个原因，利用0.25ft和0.0625ft的网格尺寸重复拉伸应力为1000psi的内聚运算。对于最精细网格，裂缝跨过不止680个网格块延伸。图18示出了三种不同网格尺寸的预测裂缝长度。可以看到，随着网格被细化，结果接近分析方程。对于这种情况，预计极精细网格结果应该与分析预测非常相似，因为内聚带相对于裂缝长度非常小。最精细网格的内聚带跨过与0.125ft的长度相对应的两个网格块延伸。为本例进行的精细网格仿真不是恰好渐近到分析结果，因为该仿真具有流体压缩率，以及储层中的流体的体积比表面条件下的流体的体积小近似1％。对于最精细网格(0.0625ft单元)，最终裂缝宽度，流体压力，和裂缝长度分别是0.002886ft，4039psi，和42.6ft，它们分别有利地与0.002836ft，4037psi，和43.8ft的分析结果相当。

比较两种不同延伸判据的速度、稳定性、和准确性。应力强度因子手段运算得较快，呈现出比内聚模型小的稳定性问题，而对于在本研究中所使用的离散化，内聚模型看起来更加准确。对于1ft网格块的仿真，拉伸应力为300psi的内聚模型的运算比使用应力强度因子的手段多花费大约40％的时间，并且，对于使用大时间步长的运算，内聚模型也更不稳定。当将相同时间步长大小用于两种运算时，内聚运算的附加CPU时间似乎是由线性叠代的次数增加引起的。尽管内聚选项运算比应力强度因子选项更慢，但内聚运算看起来产生更准确的结果。例如，对于应力强度因子运算和内聚运算，粗网格运算(1ft单元)分别预测出0.003014ft和0.002870ft的裂缝宽度。重复更精细网格(0.25ft单元)的运算，结果分别是0.003018ft和0.002824ft。应力强度因子手段随着更精细网格划分看起来并不向0.002836ft的分析值收敛，而内聚手段在更精细网格上呈现出改善的结果。这或许是因为应力强度因子手段对于裂缝末端处的位移存在平方根奇点，以及正好在裂缝末端处的线性近似不随着网格精细化而改善。如果在末端后面的固定距离处而不是在与末端最接近的节点处评估裂缝张开，那么，或许可以看到应力强度因子手段的改善的收敛。在一些程序中，通过使用裂缝末端附近的特殊单元来避免这种收敛问题。内聚单元或许在更精细网格上呈现出改善的结果，因为它们在裂缝末端处不经历相同的奇异行为。内聚单元通过限制裂缝末端附近的应力的幅度来避免奇异行为。

7.1.2具有高粘度流体的平面应变裂缝

这个例子使用大流体粘度来审视裂缝生长预测。具有大流体粘度增加了干旱带出现在裂缝末端处的机会。本例中的网格与前面的精细网格例子(1/4ft网格块)中的相同，以及注入井位于沿着x，y，和z方向分别是1000ft x 1000ft x 20ft的正方形图案的一角中。所有参数都与前例中的参数相同，但弹性刚度是15.0E6psi，以及应力强度因子是500psi√in。

当进行包含干旱带的仿真时，将干旱带中的流体压力的输入参数供应给程序。给定干旱带中的流体压力，程序在仿真期间自动确定干旱带的尺寸。本例的最小初始原地应力是4000psi，因此，运算使用3000，2000，1000，和0psi干旱带压力的情况。对于干旱带压力为3000psi的情况，程序预测裂缝长度为51.25ft，以及裂缝末端处的干旱带长1.9ft。对于0.25ft的网格块尺寸，这对应于有7个网格块完全干燥，以及在流体的前沿上有一个网格块充着41％流体。图19示出了裂缝末端附近区域中的流体压力和裂缝半宽度。程序生成单元中心处的流体压力，但图19中的值被内插成处在每个单元的边缘上，以便在相同位置上画出压力和宽度。裂缝末端处的平滑截止是由干旱带的存在引起的。对于这种情况，井上的流体压力是4443psi，因此，沿着裂缝存在1443psi压降，大部分压降出现在裂缝端部附近。

重复相同运算，但应力强度因子为1500psi√in，这给出50.75ft的裂缝长度和1.45ft的干旱带尺寸。在这个问题中，应力强度因子对结果没有显著影响。对于应力强度因子为500psi√in和干旱带压力为2000，1000，和0psi的运算，裂缝长度都是49.0ft，而干旱带长度对于2000和1000psi情况小于两个网格块(0.5ft)，和对于0psi情况小于一个网格块(0.25ft)。

因此，可以表明，当在裂缝中存在明显压降时，干旱带可以自然地出现在延伸裂缝的末端处。干旱带的尺寸取决于有关干旱带中的压力的假设，但在大多数情况下，干旱带的长度比裂缝长度小得多。另外，本例表明，对于牵涉到沿裂缝的大压降的问题，裂缝延伸判据可能起很小的作用。

7.1.3具有低粘度流体的币状裂缝

本例仿真币状裂缝的生长。它例示了程序模拟平面中的一般平面状裂缝生长的能力。将1cp的粘度用于流体，以便流体压力在整个裂缝中几乎是均匀的。程序能够管理非结构化网格，但将正方形网格块用于这种仿真，以展示即使在使用规则离散化时，程序也可以再现裂缝的圆形轮廓的能力。该仿真只包括四分之一裂缝，并且只模拟总裂缝宽度的一半，因为裂缝处在3D网格的边界上。储层的3D网格在三个方向的每一个上具有40个网格块，其中每个网格块是5ft x 5ft x 5ft。储层是非多孔的，取消(即，不考虑)重力，并且在网格的角上完成井。裂缝网格在x和z方向上也具有40个网格块，并且与储层网格的xz面匹配。弹性模量和泊松比分别是1.0E6psi和0.25，最小原地应力是2000psi，以及应力强度因子是1500psi√in。

以0.25bbl/min的表面速率(总井速率是2bbl/min)注入流体达12分钟。那时，裂缝生长到85ft的半径。图20示出了仿真中途(在6分钟时)的裂缝半宽度。图21示出了仿真结束时(在12分钟时)的裂缝半宽度。即使网格由正方形单元组成，该程序在近似圆形裂缝方面也做得很出色。

当使用币状裂纹的施奈登(Sneddon)解以及币状裂纹的应力强度因子满足方程(其中Δp是过压力)的事实时，以均匀压力延伸币状裂缝的半径由下式给出：

R=(3QEt2π(1-v2)KIC)2/5,.......(35)]]>

其中，R是裂缝半径，和Q是四分之一裂缝的表面上的注入速率(两倍于用在仿真中的速率)。方程35中的假设是均匀受压币状裂缝的应力强度因子通过给出，其中Δp是过压力；裂缝是椭圆形，注入体积与四分之一裂缝的裂缝体积之间的关系通过Qt＝πw₀R²/6给出，其中w₀是井筒上的总裂缝宽度；以及井筒上的裂缝宽度通过8(1-v²)RΔp/πE给出。图22使用方程35来将预测裂缝长度与分析裂缝长度相比较，其中虚线对应于分析结果。图22中的仿真结果是裂缝末端沿着x轴的位置。图22中的点代表仿真时每个时间步长的裂缝长度。该仿真使用5ft网格块，因此，裂缝以5ft增量前进，这是在曲线图上看到离散裂缝阶跃的原因。尽管对于这个应用来说，网格相当粗，但在仿真与分析结果之间匹配得很好。在12分钟之后，计算结果表明在井筒上具有2038.3psi的最大压力，和在没有干旱带迹象的末端处具有2037.4psi的最小压力。在12分钟上数值结果有利地与预测裂缝中2040.6psi的流体压力的分析结果相当。

有关币状裂缝生长的本例表明了程序模拟平面中的一般裂缝生长的能力。即使将长方形网络用于裂缝，该程序也可以产生良好的圆形裂缝表示，并且预测裂缝生长很好地与简单分析结果相当。

7.1.4数值结果的分析

在下面的数值结果中，基于应力强度因子的裂缝延伸判据看起来是预测裂缝生长的最稳定并且最迅速的技术，但当网格划分精细到足以求解裂缝末端处的内聚带时，看起来没有内聚技术那么准确。当网格精细到足以求解裂缝末端处消除了应力奇点的内聚带时，基于内聚带的裂缝延伸判据产生准确结果。但是，对于内聚带小的问题，可能需要极精细网格。仿真结果表明，内聚带的尺寸可以通过两倍于方程33中的粗略估计的方程(π/4)(K_IC/σ_c)²来估计。当网格精细到足以求解内聚带时，以及当内聚带与裂缝长度相比很小时，两种裂缝延伸判据产生相似结果。

8.参考文献引用

本文引用的所有参考文献通过全文引用并入本文中，并且同等程度地用于所有目的，就像每个单独公告、专利或专利申请被具体和分别指示成通过全文引用并入本文中用于所有目的一样。本文对参考文献的讨论或引用不应该理解为承认这样的参考文献是本发明的现有技术。

9.修改

对于本领域的普通技术人员来说，显而易见，可以不偏离本发明的精神和范围对本发明作许多修改和改变。通过例子只提供了本文所述的特定例子，而本发明只受所附权利要求书，以及与这样权利要求书同等权利的等同物的整个范围限制。

作为本文所述的系统和方法的宽范围的例示，本文所述的系统和方法可以通过包含设备处理子系统可执行的程序指令的程序代码，在许多不同类型的处理设备上实现。该软件程序指令可以包括源代码、目标代码、机器码、或可起使处理系统执行本文所述的方法和操作作用的任何其它存储数据。但是，也可以使用其它实现，譬如，固件，或甚至配置成实现本文所述的方法和系统的适当设计硬件。

该系统和方法的数据(例如，关联、映射、输入数据、输出数据、中间数据结果、最终数据结果等)可以存储和实现在一种或多种不同类型的计算机实现数据存储装置中，譬如，不同类型的存储设备和编程结构(例如，RAM、ROM、闪速存储器、平面文件、数据库、编程数据结构、编程变量、IF-THEN(或类似类型)语句结构等)。应该注意到，数据结构描述用在组织数据和将数据存储在数据库、程序、存储器、或其它计算机可读介质中供计算机程序使用中的格式。

该系统和方法可以配备在许多不同类型的计算机可读介质上，该计算机可读介质包括计算机存储机构(例如，CD-ROM、软盘、RAM、闪速存储器、计算机硬盘驱动器等)，该计算机存储机构包含用在处理器的执行中以便执行本文所述的方法的操作和实现本文所述的系统的指令(例如，软件)。

本文所述的计算机部件、软件模块、函数、数据存储装置和数据结构可以相互直接或间接连接，以便使它们的操作所需的数据流动。还应该注意到，模块或处理器不局限于此地包括执行软件操作的代码的单元，并且可以，例如，实现成代码的子例程单元，实现成代码的软件函数单元，实现成对象(像在面向对象的范式中那样)、实现成小应用程序，用计算机脚本语言实现，或实现成另一种类型的计算机代码。取决于未来处境，软件部件和/或函数可以处在单台计算机或分布式跨越多台计算机上。

应该懂得，正如用在本文的描述中和所附的整个权利要求书中那样，除非上下文另有明确指明，单数形式“一个”、“一种”、和“该”的含义也包括复数指代物。此外，正如用在本文的描述中和所附的整个权利要求书中那样，除非上下文另有明确指明，“在...中”的含义包括“在....中”和“在...上”。最后，正如用在本文的描述中和所附的整个权利要求书中那样，除非上下文另有明确指明，“和”和“或”的含义包括连接性和分离性两者，并且可互换使用；短语“不可兼的或”可以用于指示只可应用分离性含义的状况。

附录A

概要

●程序描述和应用

●裂缝和储层的耦合

●裂缝延伸判据和例子

●币状裂缝例子

程序描述

●重油和复合渗流

●多孔弹性/多孔塑性力学

●热传导和对流

●水力裂缝生长

●废料处理

●非结构化网格

应用

●水力和热力压裂

●沉陷

●废料处理

●砂石产生和空腔生成

●完成策略

●CO₂处理

裂缝模型

●裂缝中的平行板之间的层流

●从裂缝面到储层的渗流

●裂缝压力将垂直牵引力施加在储层上

●裂缝前面的裂缝平面中的零垂直位移

●基于应力强度因子或内聚单元的节点释放

确关内聚模型比较的问题

●弹性比较需要相对于裂缝长度的小内聚带

对于σ_c＝100psi：       L_c＝15ft

但                      L_f＝44ft

●准确性需要将内聚带离散化

对于σ_c＝1000psi：L_c＝0.125ft

但                ΔX＝1.0fta

有关裂缝延伸判据的注释

●对于弹性问题，两种技术产生相似结果

●K_IC运算需要比内聚运算更短的CPU时间

并且K_IC运算比内聚运算更稳定

●在当前实现的网格精细化期间内聚模型运算产生更准确结果

●K_IC运算局限于弹性问题，或者塑性带相对于裂缝长度很小的塑性问题

注释

●已经开发出模拟水力裂缝生长的全耦合地质力学模型

●在适当条件下应力强度因子和内聚单元产生相似结果

●内聚单元具有更大适用范围，但需要更多CPU时间并且更不稳定

●可以使用简单长方形网格再现径向形状