CN102575510B - 控制地质力学油藏系统中的出砂的计算机实施的系统和方法 - Google Patents

Abstract

提供系统和方法用于在预测地质力学油藏系统中的出砂中使用。出砂预测的计算可包括对模拟地质力学油藏系统的偏微分方程组求解。还提供系统和方法用于在基于出砂预测来操作地质力学油藏系统中使用，以便控制地质力学油藏系统中的出砂。

Description

控制地质力学油藏系统中的出砂的计算机实施的系统和方法

技术领域

本文涉及在预测地质力学油藏系统中的出砂中使用的计算机实施的系统和方法。本文还涉及在基于所述预测来控制地质力学油藏系统中的出砂中使用的系统和方法。

背景技术

油藏系统的生产通常是在油藏开发之后进行的阶段，在该阶段中，诸如碳氢化合物(油或气)之类的油藏流体被从油藏引出。出砂是指伴随油藏流体产生地层固体颗粒的事件。用于描述可能与油藏流体一起产生的小地层(井眼周围的岩石)颗粒的一般术语是“砂”。术语“微粒”已用于一些文献中。油藏地层材料通常包括孔隙度和渗透率足以储存和传输诸如油或水之类的油藏流体的岩石类型。由于沉积岩是多孔的并且在可以保留化石残迹(从其得到碳氢化合物)的温度条件下形成，所以它们是最常见的油藏岩石类型(而不是火成岩和变质岩)。沉积岩的示例包括但不限于砾岩、砂岩、粉砂岩、页岩、石灰岩、白云岩、岩盐(盐)、盐、石膏岩以及硫酸钙硬石膏岩。沉积岩可以包括多种矿物质，包括但不限于石英、长石、方解石、白云石以及黏土组矿物质。

出砂是生产期间油藏流体(例如油)的流动所导致的地层砂的迁移。出砂可能由油藏地层材料内部的剪切破坏或拉伸破坏造成。剪切破坏可能在井眼压力显著减小时发生(例如在井的寿命晚期)，压力减小增大了井眼附近的应力，导致地层破坏。拉伸破坏可能在井眼附近的油藏地层材料的孔隙度和渗透率显著受损或者流速极高时发生。在任一拉伸破坏条件下，流动流体可能对地层中的单个颗粒施加很大的拖曳力，如果拖曳力过大的话，则能导致单个颗粒之间的粘结受到破坏，造成拉伸破坏和出砂。

在很多情况中，由于出砂会限制产能、侵蚀完井组件、阻塞井眼通道、干扰钻井设备操作以及给出严重的处理难题，所以出砂是不期望的。出砂增大了射孔腔或井眼的直径，减小了套管周围的支撑。结果，射孔坍塌，腔变大，最终从井眼的生产可能停止。如果出砂严重，则可能需要补救措施来完全控制或防止出砂，诸如砾石充填或固砂。在极端情况下，可能发生大规模出砂，其中出砂不受控制地增大，最终完全侵蚀形成井基础的油藏材料。

通常，在油藏开发项目早期就会采取措施来尽可能地完全防止出砂。例如，具有多个产油层的疏松地层可以是套管完井。套管完井过程包括构造围绕油藏井眼的屏障系统，以便防止或推迟出砂的开始或程度，这包括将钢管封固到井眼中以及使用网，例如可扩展砂筛技术。然而，这些措施可能对井的产能产生负面影响。也就是说，套管井的产能可能比裸眼井低很多。此外，由于需要用于实施的资源、时间和人工，所以这些措施增大了油藏开发项目的成本。

允许一定数量的出砂会有助于提高井的产能。因此，在预测地质力学油藏系统的出砂中使用的方法将会是有用的，因为它将会允许在井开发项目早期确定是否及何时将发生出砂，以及出砂的程度。通过使用这些预测，可以构造和操作油藏从而产生有限数量的砂，实质上维持油藏中产生的腔且提高井的产能。

发明内容

如这里公开的那样，提供计算机实施的系统和方法用于在预测地质力学油藏系统的出砂中使用。所述方法和系统包括：接收指示地质力学油藏系统内材料的逐渐出砂和大规模出砂的数据；在计算机系统上，基于包括硬化模型的地质力学模型对所接收的大规模出砂数据的拟合来计算临界塑性应变的值；在计算机系统上，基于包括出砂函数的地质力学模型对所接收的逐渐出砂数据的拟合并使用临界塑性应变值以及一个或多个硬化参数中的至少一个值来计算出砂函数的至少一个参数值；其中一个或多个硬化参数中的所述至少一个值基于硬化模型对指示地质力学油藏系统内材料的塑性变形的数据的拟合来计算；其中所述硬化模型地质力学油藏系统内的材料的塑性变形行为建模；且其中所述出砂函数预测地质力学油藏系统的出砂。一个或多个硬化参数中的至少一个值、临界塑性应变值或者出砂函数的至少一个参数值可输出到用户接口设备、监视器、计算机可读存储介质、本地计算机或者是作为网络的一部分的计算机。

还提供计算机实施的系统和方法用于在预测地质力学油藏系统的出砂中使用，包括：接收指示地质力学油藏系统内材料的塑性变形、逐渐出砂和大规模出砂的数据；在计算机系统上，基于硬化模型对所接收的塑性变形数据的拟合来计算一个或多个硬化参数的值；其中所述硬化模型对地质力学油藏系统内材料的塑性变形行为建模；在计算机系统上，基于包括硬化模型的地质力学模型对所接收的大规模出砂数据的拟合来计算临界塑性应变值；在计算机系统上，基于包括出砂函数的地质力学模型对所接收的逐渐出砂数据的拟合且使用硬化参数的至少一个值的值以及临界塑性应变值来计算出砂函数的至少一个参数；其中所述出砂函数预测地质力学油藏系统的出砂。一个或多个硬化参数中的至少一个值、临界塑性应变值或者出砂函数的至少一个参数的值可输出到用户接口设备、监视器、计算机可读存储介质、本地计算机或是作为网络的一部分的计算机。

在前述方法和系统的一个方面中，出砂函数是函数f(x)，其中当x＝0时f(x)＝0，当x＝1时f(x)＝1；且其中x是临界塑性应变的函数。在前述方法和系统的另一个方面中，出砂函数的至少一个参数的值是出砂函数的指数值。所述出砂函数可以由下面的表达式给出：

$f\left(x\right)={\left(\frac{{\mathrm{\ϵ}}^p}{{\mathrm{\ϵ}}_{\lim }^p}\right)}^m$

其中其中ε^p是塑性应变不变量，其中是临界塑性应变，以及其中m是指数值。

硬化模型可以是修正的Bigoni-Piccolroaz模型。塑性变形数据可以从三轴试验获得。逐渐出砂和大规模出砂可以从空心圆柱试验获得。所述硬化模型对所接收的塑性变形数据的拟合可以通过回归来获得。

在前述方法和系统的一个方面中，包括硬化模型的地质力学模型对所接收的大规模出砂数据的拟合是通过求解对地质力学油藏系统建模的偏微分方程组来获得的；其中所述偏微分方程组包括油藏流动模型以及包含硬化模型的地质力学模型；其中所述偏微分方程组通过完全展开的雅可比行列式耦合；且其中对偏微分方程组的求解包括基于所接收的大规模出砂数据来在单个时间步长中同时求解完全展开的雅可比行列式。

在前述方法和系统的另一方面中，包含出砂函数的地质力学模型对所接收的逐渐出砂数据的拟合是通过求解对地质力学油藏系统建模的偏微分方程组来获得的；其中所述偏微分方程组包括油藏流动模型以及包含硬化模型的地质力学模型；其中所述偏微分方程组通过完全展开的雅可比行列式来耦合；且其中对偏微分方程组的求解包括基于所接收的逐渐出砂数据来在单个时间步长中同时求解完全展开的雅可比行列式。

还提供计算机实施的系统和方法用于在预测地质力学油藏系统的出砂中使用，包括：接收指示与地质力学油藏系统内的材料相关的物理或化学属性的数据；通过求解对地质力学油藏系统建模的偏微分方程组来生成出砂预测；其中所述偏微分方程组包括地质力学油藏系统的油藏流动模型和地质力学模型；其中所述地质力学模型包括硬化模型；其中出砂判据应用到所述地质力学模型；其中所述偏微分系统通过完全展开的雅可比行列式耦合；其中求解偏微分方程组包括基于所接收的物理属性数据来在单个时间步长中同时求解完全展开的雅可比行列式；且其中所述生成在计算机系统上实施。所生成的出砂预测可输出到用户接口设备、监视器、计算机可读存储介质、本地计算机或是作为网络的一部分的计算机。出砂判据可以是总应变的临界值、塑性应变不变量的临界值或是最大有效应力。

此外，还提供计算机实施的系统和方法用于在预测地质力学油藏系统中的出砂中使用，包括：接收指示与地质力学油藏系统内的材料相关的物理或化学属性的数据；定义包含多个网格单元的网格；通过求解对地质力学油藏系统建模的偏微分方程组来生成出砂预测；其中所述偏微分方程组包括地质力学油藏系统的油藏流动模型和地质力学模型；其中所述地质力学模型包括硬化模型；其中出砂判据应用到所述地质力学模型；其中所述偏微分方程组通过完全展开的雅可比行列式来耦合；其中求解所述偏微分方程组包括基于所接收的物理属性数据来在单个时间步长中同时求解完全展开的雅可比行列式；其中所述油藏模型和地质力学模型在所述网格上计算；且其中所述生成在计算机系统上实施。所生成的出砂预测可输出到用户接口设备、监视器、计算机可读存储介质、本地计算机或是作为网络的一部分的计算机。出砂判据可以是总应变的临界值、塑性应变不变量的临界值或是最大有效应力。

本公开的一个方面提供一种用于执行包括前述系统和方法在内的在这里公开的任一方法和系统的步骤的计算机系统。所述计算机系统包括一个或多个处理器单元；以及与所述一个或多个处理器单元连接的一个或多个记忆单元，所述一个或多个记忆单元包括一个或多个模块，所述一个或多个模块包含使所述一个或多个处理器单元执行步骤的一个或多个程序，所述步骤包括：执行包括前述系统和方法在内的在这里公开的系统和方法中的任一个的步骤。在前述实施方式中，所述一个或多个记忆单元可以包括一个或多个模块，所述一个或多个模块包括使所述一个或多个处理器单元执行步骤的一个或多个程序，所述步骤包括将所述系统和方法的结果输出至显示器、用户接口设备、有形的计算机可读数据存储产品或是有形的随机存取存储器。例如，如可应用到正在执行的方法那样，所输出的系统和方法的结果可以是一个或多个硬化参数中的至少一个值、临界塑性应变的值、出砂函数的至少一个参数的值或者所产生的出砂预测。

本公开的另一方面提供一种存储计算机程序的计算机可读介质，该计算机程序可被计算机运行以用于执行包括前述系统和方法在内的在这里公开的系统和方法中的任一个的步骤。该计算机程序产品被提供以用于与具有一个或多个记忆单元以及一个或多个处理器单元的计算机结合使用，所述计算机程序产品包括编码有计算机程序机制的计算机可读存储介质，其中所述计算机程序机制可被载入计算机的一个或多个记忆单元中并使计算机的一个或多个处理器单元执行步骤，所述步骤包括：执行包括前述系统和方法在内的在这里公开的任一系统和方法的步骤。在前述实施方式中，计算机程序机制可被载入所述计算机的一个或多个记忆单元中并促使计算机的一个或多个处理器单元执行步骤，所述步骤包括将系统或方法的结果输出至显示器、用户接口设备、有形的计算机可读数据存储产品或是有形的随机存取存储器。例如，如可应用到所执行的方法那样，所输出的系统或方法的结果可以是一个或多个硬化参数中的至少一个值、临界塑性应变的值、出砂函数的至少一个参数的值或是所生成的出砂预测。

还提供用于地质力学油藏系统的出砂的系统和方法，包括：根据包括前述系统和方法在内的在这里公开的任一方法和系统的结果来操作地质力学油藏系统。

还提供系统和方法用于操作地质力学油藏系统以控制地质力学油藏系统的出砂。所述系统和方法包括：基于包含出砂函数的地质力学模型对指示与油藏系统内的材料相关的物理或化学属性的数据的拟合来计算至少一个操作参数的值；其中所述出砂函数预测地质力学油藏系统的出砂；且其中操作参数的所述至少一个值指示地质力学油藏系统的稳定出砂条件；以及根据所述至少一个操作参数的值来操作地质力学油藏系统。

附图说明

图1是在预测地质力学油藏系统的出砂中使用的方法的流程图。

图2A示出在三轴试验中用于对样品200进行偏差加载的应力分量(σ_x、σ_y、σ_z)(在该示例中，σ_x＝σ_z)。

图2B示出用于固结试验的应力分量(σ_y)，其中样品202径向约束在限制环204中，且仅在y方向上被加载(σ_y)。

图3示出空心圆柱试验的示例性示意图。

图4示出空心圆柱试验的结果，绘制为出砂(以立方厘米(cc)为单位)对比时间(秒)的图以及使用临界应力值＝0.014和指数(m)＝1的材料常数对结果的拟合。对曲线的大规模出砂部分和逐渐出砂部分二者都进行了空心圆柱出砂试验数据的拟合。

图5A示出射孔试验的示例性示意图。

图5B示出来自射孔试验的岩芯样品。

图6示出射孔试验的结果，绘制为出砂(以立方厘米(cc)为单位)对比时间(秒)的图以及使用不同的指数(m)值对结果的拟合。圆柱岩芯样品用半球状末端来射孔。

图7是在基于出砂预测结果来操作地质力学油藏系统中使用的方法的流程图。

图8是在对地质力学油藏系统的出砂进行建模中使用的示例方案的框图，包括油藏模型和地质力学模型。

图9是在对地质力学油藏系统的出砂进行建模中使用的示例方案的框图，包括油藏模型、地质力学模型和热模型。

图10A示出在模型计算中使用的三维网格的示例。

图10B示出在模型计算中使用的二维网格的示例。

图11A示出使用剪切硬化的Drucker-Prager硬化模型中屈服面的移动曲线。

图11B示出使用帽硬化的Drucker-Prager硬化模型中屈服面的移动曲线。

图12A示出对于四个不同的偏差(K)值，修正Drucker-Prager模型的八面体(偏差的)平面中屈服面的曲线。

图12B示出对于贝塔值(β)与偏差值(γ)的四种不同组合，修正Bigoni-Piccolroaz模型的八面体(偏)平面中的屈服面的曲线。

图13示出使用模型对出砂的示例计算。

图14示出用于实施这里公开的出砂预测方法的示例计算机系统。

图15示出用于盐水试验(使用例如图3所示的空心圆柱试验设备)的岩芯样品，包括样品尺寸。

图16A示出对于盐水饱和的岩芯样品，三轴试验数据(黑)和来自仿真(灰色)的数值结果(数据拟合)的曲线。

图16B示出图16A所示曲线的一部分的特写，其中岩芯样品材料成为几乎完美塑性的。

图17示出对于煤油饱和的样品，三轴试验结果(黑)和来自仿真(灰色)的数值结果的曲线。

图18示出应用到仿真空心圆柱试验中的加载的岩芯样品仿真的边界条件，其中(A)示出力学边界条件，(B)示出流动边界条件，这些条件应用到仿真中的样品。

图19示出对于盐水饱和样品，与时间对照的所测量的离散约束应力、流速以及流动压力。

图20示出作为流速的函数的所计算的渗透率的曲线。

图21示出作为约束应力的函数的所计算的渗透率的曲线。

图22示出显示了作为时间函数的离散的模拟流速的匹配的曲线图。

图23A示出对于Drucker-Prager型材料模型(贝塔＝0且偏差＝0)，偏差平面中的屈服面的形状。

图23B示出对于Lade型材料模型(贝塔＝0且偏差＝0.95)，偏差平面中的屈服面的形状。

图24示出从岩芯样品的测量得到的出砂体积(cc)与时间(sec)的关系曲线(黑线)以及一模型使用0.05mm、0.01mm和0.005mm的网格尺寸对结果的不同拟合。0.01mm的网格尺寸足够小以捕获模型细节且消除网格依赖性。所使用的材料常数是：偏差＝0.5，贝塔＝0，临界应变值＝0.017以及指数＝2。

图25示出从岩芯样品的测量得到的出砂体积(cc)与时间(sec)的关系曲线(黑线)以及使用偏差＝0.5、贝塔＝0、临界应变值＝0.017以及指数＝2的材料常数对结果的拟合。对曲线的大规模出砂部分和逐渐出砂部分二者都进行了数据的拟合。

图26示出从岩芯样品的测量得到的离散的测量水量(公升(l)/min)与时间(sec)的关系曲线(黑线)以及使用偏差＝0.5、贝塔＝0、临界应变值＝0.017以及指数＝2的材料常数对结果的拟合。

图27示出对于盐水试验，有效塑性应变与硬化(psi)的关系曲线。

图28示出从岩芯样品的测量得到的出砂体积(cc)与时间(sec)的关系曲线(黑线)以及对结果的不同拟合，其中偏差和临界应变限制如所示地改变，贝塔和指数分别保持恒定在0和1。

图29示出从岩芯样品的测量得到的出砂体积(cc)与时间(sec)的关系曲线(黑线)以及对结果的不同拟合，其中偏差＝0.5，贝塔＝0，临界应变值＝0.017，且指数(m)的值如所示地改变。

图30示出从岩芯样品的测量得到的出砂体积(cc)与时间(sec)的关系曲线(黑线)以及对结果的不同拟合，其中偏差＝0.5，贝塔＝0，指数＝1，且临界应变值改变。

图31示出从岩芯样品的顶部到底部取得的水平切片，显示了腔的进展。该切片在出砂试验结束时取得。

图32示出随时间变化的腔进展的数值仿真。每个图是轴向对称的垂直切片，其中灰色区域显示腔的位置和范围，黑色区域显示完好岩体(或间隔物)。仿真中使用的材料常数是：偏差＝0.5，贝塔＝0，临界塑性应变值＝0.017以及指数＝2。

图33A示出利用偏差＝0.5、贝塔＝0、临界应变值＝0.017以及指数＝2的材料常数计算的出砂体积(cc)与时间(sec)的关系曲线。

图33B示出腔尺寸的数值仿真，其是从116675秒开始对于46MPa的恒定约束应力和0.1643398MPa的恒定压差，且利用偏差＝0.5、贝塔＝0、临界应变值＝0.017以及指数＝2的材料常数计算的。灰色区域表示腔的位置，黑色区域表示完好岩体(或间隔物)。

图34示出出砂(cc)与时间(sec)的关系的数值仿真，其是从113230秒开始对于44MPa的恒定约束应力和1.255317MPa的恒定压差，且利用偏差＝0.5、贝塔＝0、临界应变值＝0.017以及指数＝2的材料常数计算的。

具体实施方式

提供系统和方法以用于在预测地质力学油藏系统的出砂中使用。所述系统和方法使用所测量的油藏属性来产生这种出砂预测，由于所述预测允许在油藏开发项目早期确定是否以及何时会发生出砂且以何种速率发生，因此是有用的。

这种出砂预测可以用于在油藏开发项目早期确定完井技术类型，该完井技术可被实施以在井的整个寿命期间彻底消除出砂或是允许一定数量的逐渐出砂(出砂的逐渐时间演变)，但该完井技术不会导致大规模出砂。例如，可以确定在油藏的井眼周围使用的屏障系统的设计，例如但不限于砂筛技术或砾石充填。可以确定具有多个产油层的疏松地层是否应是套管完井。由于与没有出砂预测时使用的出砂缓解设备相比，可以根据出砂允许程度来确定安装更少的出砂缓解设备，所要这种出砂预测能有助于减少与完井相关的成本。出砂预测还可以用于确定如何操作油藏，例如但不限于生产压差、产量、最小井底压、产层温度以及井眼中的流体流动压，以便实现期望的出砂体积。此外，出砂预测可以用于确定井的寿命中使用缓解出砂的技术以及安装缓解出砂的设备的时间点。

这里公开的出砂预测系统和方法是以物理为基础且可以用于预测开始出砂、逐渐出砂以及大规模出砂的现象。出砂预测系统和方法是确定性的，也就是说，它们是以物理原理而不仅是以相关性为基础。例如，基于物理的出砂判据应用到确定是否发生出砂、出砂体积以及出砂速率的计算。

根据出砂预测结果来操作油藏可以通过降低完井成本而节约可观的费用。出砂预测结果可以用于出砂管理。出砂管理可以用于例如决定所要安装的防砂完井的类型和成本，例如套管井和射孔，或是常规的支撑裂缝完井，它们产生有限的砂量。

5.1出砂预测系统和方法的示例

图1的流程图示出在预测地质力学油藏系统的出砂中使用的示例系统和方法中的步骤。

步骤100。在步骤100中，接收指示地质力学油藏系统内的材料的物理或化学属性的数据。例如，可以接收指示在一个或多个材料出砂阶段期间地质力学油藏系统中的材料的物理或化学属性的数据。在图1的示例中，指示地质力学油藏系统中材料的逐渐出砂的数据在步骤100A中被接收。逐渐出砂指的是油藏中材料出砂的逐渐时间演变。在步骤100B中，接收指示地质力学油藏系统中材料的大规模出砂的数据。大规模出砂指的是出砂的急剧增大。此类数据包括但不限于出砂速率(以随时间变化的出砂体积为单位)。此类数据的其他示例包括但不限于地层材料的类型、地层材料的孔隙率、地层材料的渗透性、油藏中流体的类型、孔隙压力、温度、井眼中流体的粘度、产层的温度、井眼中的流体流动压力、井眼中流体流动的拖曳力以及井眼中流体流动的类型。此类数据的其他示例包括但不限于井深、油比重、油粘度、岩石类型的有效厚度、当前油藏压力、最小油含量、油饱和率、岩石的渗透性和孔隙率、系统温度、岩层传输率、水的盐度、已有裂缝系统、气帽、倾角、井距、接受度、碳氢化合物(HC)成分、最低互溶压力、压力比、气饱和度、起泡点压力、临界气饱和度、气体比以及垂向波及因子。

所接收的数据可以指示地质力学油藏系统中材料的塑性变形。也就是说，数据可以是提供地层材料的塑性变形的量度的一个或多个参数值。还可以接收指示油藏系统中材料的弹性行为的数据。此类数据的示例是但不限于杨氏模量、屈服强度、材料的应力-应变曲线、极限强度、应变硬化行为、颈缩行为、断裂点。在一示例中，可以接收指示塑性变形起始(也就是从弹性到塑性行为的转变点)的数据。例如，屈服强度可用于准确确定材料的弹性极限，如果超出该极限，那么材料上的额外应力会导致产生永久(塑性)变形。

材料的屈服判据的确定还可以用于确定材料的塑性变形的起始。屈服判据(其可显示为屈服面)可用于指示在不同的应力组合下的材料弹性极限(以及塑性变形起始)。能应用到各向同性材料(即在所有方向上都具有一致属性的材料)的屈服判据的示例可以是以最大主应力、最大主应变、最大剪应力、总应变能以及畸变能为基础的判据。对于以最大主应力为基础的屈服判据来说，屈服可视为在施加到材料的最大主应力超出了单轴拉伸屈服强度时发生。对于最大主应变判据来说，屈服可视为在材料的最大主应变达到与简单的拉伸试验期间的屈服点对应的应变时发生。对于最大剪切应力屈服判据(Tresca屈服判据)来说，屈服可视为在应用到材料的切应力超出了切变屈服强度时发生。在总应变能屈服判据中，可假设在屈服点与弹性变形相关的储能与具体的应力张量无关，从而屈服在每单位体积的应变能大于简单拉伸中的弹性极限处的应变能时发生。对于畸变能屈服判据(VonMises屈服判据)来说，屈服可视为在材料的形状畸变超出拉伸试验的屈服点时发生。应用于各向同性材料的屈服判据的其他示例是Mohr-Coulomb屈服判据、Drucker-Prager屈服判据以及Bresler-Pister屈服判据。可以如应用到各向异性材料(即塑性屈服行为表现出方向相关性的材料)那样的屈服判据的示例包括但不限于Hill二次屈服判据、广义Hill屈服判据以及Hosford屈服判据。

指示油藏系统中的材料的弹性行为、塑性变形行为、逐渐出砂或大规模出砂的数据可从试验获取，例如但不局限于三轴压缩试验、三轴拉伸试验、单轴应变试验、固结试验以及静水压缩试验。图2A示出三轴试验中可施加到材料200的圆柱的应力分量(σ_x、σ_y、σ_z＝σ_x)的方向和关系。图2B示出在固结试验中能沿y方向施加到材料202的圆柱的应力分量(σ_y)；在压缩加载期间材料可限制在限制环204中，从而没有应力沿x、z方向施加。这两个测试中的任一个或是二者都可用于提供指示油藏系统中的材料的弹性行为、塑性变形行为、逐渐出砂和/或大规模出砂的数据。数据可以从对材料执行的试验获得，试验涉及用于材料的不同加载路径。数据还可以从空心圆柱试验获取，在该试验中，流动试验以及轴向、球面或扭转应力的不同组合被应用到材料的空心圆柱样品。图3示出空心圆柱试验设备的示意图。图4示出从对材料执行的空心圆柱试验获取的数据，且还示出数据中的塑性变形起始点、逐渐出砂区域以及大规模出砂区域。数据还可以从射孔试验获取，在该试验中，在材料圆柱样品的一端已被射孔(即经历聚能射孔弹)之后对材料的圆柱样品进行流动试验。流动试验可用油、气、盐水或这些流体的任何组合来进行。图5A示出射孔试验设置的示意图，图5B示出已经经历射孔试验的材料的样品。图6示出从对材料执行的空心圆柱试验获取的数据，还示出数据中的塑性变形起始点、逐渐出砂区域以及大规模出砂区域。

数据可以从对从实际油井位置或建议油井位置取得的一个或多个岩芯样品执行的试验获取。例如，数据可以从对材料的岩芯样品执行的试验获取，岩芯样品可以从井眼获取。这种岩芯样品可以从代表可能发生出砂的地层的不同油藏部分获取。在另一示例中，数据可以从对合成样品执行的试验获取，合成样品创建为与来自油井(例如来自可能发生出砂的油井地层区域)的实际地层材料具有相似的物理和化学属性。

步骤102。基于至少一个硬化模型对所接收数据的拟合来得到一个或多个硬化参数的值。硬化模型可以用于对地质力学油藏系统内的材料的塑性变形行为进行建模。硬化模型论述于下面的第5.5.2.4节中。硬化模型的示例包括具有切变硬化或帽盖硬化的Drucker-Prager模型，具有板状硬化的修正Drucker-Prager模型，修正Bigoni-Piccolroaz模型以及Matsuoka-Nakai模型。硬化模型的其他示例包括但不限于修正Cam-Clay模型，DiMaggio和Sandler广义帽盖模型，Lade模型，Iwan/Mroz多表面模型以及Fossum和Fredrich连续表面帽盖塑性模型。

可关于修正Drucker-Prager模型得到的硬化参数的示例是但不限于α(总应力的第一不变量的乘数)、Drucker-Prager指数、屈服常数(Г)、有效塑性应变ε^P以及偏差(K)。可关于修正Bigoni-Piccolroaz模型得到的硬化参数的示例是但不限于偏差(γ)、贝塔(β)、α以及屈服常数(Г)。

如图1所示，步骤102可以包括若干步骤。在步骤102A中，硬化模型被选择。在下一个步骤中，所选择的硬化模型拟合到所接收的数据。例如，如步骤102B所示，所选择的硬化模型可以拟合到指示油藏中的材料的塑性变形行为的接收数据(在上文中结合步骤100进行了论述)。在步骤102C中，获得对数据的硬化模型拟合的一个或多个参数的值。

硬化模型对所接收的塑性变形数据的拟合可以使用任何可应用的数据拟合方法来进行。例如，拟合可以使用诸如线性回归和非线性回归之类的回归法来执行。对数据进行回归拟合的回归程序包为本领域所知。回归可以用有限因变量来执行，可以是Bayesian线性回归、分量回归或非参数回归。拟合可以用统计方法执行。可用于执行硬化模型对所接收的塑性变形数据的拟合的程序包的示例包括但不限于(SAS Institute Inc.，Cary，NC)、(The MathWorks，Inc.，Natick，MA)、R(可以经由万维网在用于统计计算的R项目的网站上获取)以及Dap(可以经由万维网在GNU操作系统的网站上获取)。

步骤104。在步骤104中，基于包含硬化模型的地质力学模型对所接收的出砂数据的拟合来计算临界塑性应变的值。在一个示例中，出砂数据可以是指示逐渐出砂数据的接收数据或是指示大规模出砂数据的接收数据，或者是二者的组合。临界塑性应变可以是材料的塑性应变的临界值，它指示材料发生故障的点，也就是材料破碎化(rubblize)而形成砂且产生腔的点。在一示例中，临界塑性应变是材料的有效塑性应变的临界值。

一种可应用的地质力学模型可以对各向同性材料、横向各向同性材料、线性弹性材料、多孔材料、实心材料或者它们的任意组合的应力、应变和/或位移建模。在一示例中，地质力学模型可以对材料的塑性变形行为建模。以下在第5.5.2节将会论述地质力学模型的示例(包括对硬化模型的论述)。

包含硬化模型的地质力学模型对出砂数据的拟合可以使用任何可应用的数据拟合方法来执行。例如，拟合可以用诸如线性回归和非线性回归之类的回归法来执行。在另一示例中，拟合可以通过求解偏微分方程组来执行，其中偏微分方程组包括包含了硬化模型的地质力学模型。偏微分方程组还可以包括油藏流动模型(以下在第5.5.1节进行论述)。如以下在第5.7节所述，偏微分方程组可以通过完全展开的雅可比行列式来耦合。求解偏微分方程组可以包括基于诸如出砂数据之类的接收数据(参见以下的第5.7节)在单个时间步长中同时求解完全展开的雅可比行列式。偏微分方程组还可以包括热模型(以下在第5.5.3节论述)。以下在第5.7节论述可用于求解偏微分方程组的过程。

包含硬化模型的地质力学模型对出砂数据的拟合可以作为对地质力学油藏系统建模的方程组的更广泛的计算的一部分来执行，所述方程组可以用于计算在将流体注入到油藏中或是从油藏产生流体时以及将应力施加到油藏边界时出现的应力、应变和/或位移。在一些计算中，包含地质力学模型(包括硬化模型)、热模型和油藏流动模型的油藏系统模型能求解包含渗流、热流和地质力学的系统。

在一示例中，包含硬化模型的地质力学模型对所接收的大规模出砂数据的拟合可以通过求解对地质力学油藏系统建模的偏微分方程组来获取，其中偏微分方程组包括油藏流动模型以及包含硬化模型的地质力学模型，其中偏微分方程组通过完全展开的雅可比行列式来耦合，且其中求解偏微分方程组包括基于所接收的大规模出砂数据在单个时间步长中同时求解完全展开的雅可比行列式。偏微分方程组还可以包括热模型。

虽然步骤102在步骤104之前论述，但是应注意，二者中的任一个可以首先执行。如图1所示，在流程图的该阶段可以执行步骤102和104中的任一个。步骤102和104可被迭代地执行，其中执行一个步骤得到的结果可以在另一步骤中使用。举例来说，如果首先执行步骤102，那么在步骤104中使用源自步骤102中的拟合的一个或多个参数值来计算临界塑性应变。如果首先执行步骤104，那么在步骤102的拟合中使用在步骤104中计算的临界塑性应变来得出一个或多个硬化参数值。此外，步骤102和104可以重复迭代执行，以实现对数据的更佳拟合。作为非限制性示例，如果硬化模型是修正Bigoni-Piccolroaz模型，那么偏差(γ)、贝塔(β)、α或屈服常数(Г)中的一个或多个的初始值可以在在步骤100中得到并用在步骤104中以计算例如拟合到指示大规模出砂的数据的临界塑性应变值。临界塑性应变的计算值可以在步骤100中的硬化模型参数计算的二次迭代中使用。步骤100的二次迭代可以提供对所接收的数据的更好拟合。步骤100的二次迭代产生的结果可提供给步骤104，以便用于临界塑性应变计算的二次迭代，从而改善对大规模出砂数据的拟合。一旦步骤100和104中的拟合收敛到对数据的最佳匹配，则可以停止迭代。

步骤106。基于包含出砂函数的地质力学模型对所接收的逐渐出砂数据的拟合，且使用临界塑性应变值(在步骤104中论述)和至少一个硬化参数值(在步骤102中论述)，计算出砂函数的至少一个参数值。

出砂函数对在达到临界塑性应变之前材料的出砂体积建模。以下在第5.5.2.5节论述出砂函数。出砂函数可以是其值从x＝0时的f(x)＝0变化至x＝1时的f(x)＝1的任何函数f(x)。出砂函数也可以是其值能缩放以落在从x＝0时的f(x)＝0到x＝1时的f(x)＝1的范围内的任何函数f(x)。在另一示例中，出砂函数也可以是任何函数f(x)，其值可通过应用合适的变换而被变换以落在从x＝0时的f(x)＝0到x＝1时的f(x)＝1的范围内，所述变换可以是而不局限于诸如小波变换和拉格朗日变换。项x可以是临界塑性应变的任何函数g()(即)。函数f(x)可以是x的任何单调函数，包括但不局限于分数函数、幂函数、正弦函数、余弦函数、对数函数、指数函数、S函数或是它们的任何组合。在一些示例中，x可以是材料的临界塑性应变和塑性应变不变量(ε^P)二者的函数，例如但不局限于比值在一个示例中，出砂函数可以由给出，其中m是指数。

值在步骤106中计算的参数可以是能够用于表征出砂函数的任何参数。举个例子，如果出砂函数是幂函数，那么参数可以是幂函数的至少一个指数。在由给出出砂函数的某个示例中，参数可以是指数m。在上述示例中，值在步骤106中计算的出砂函数的至少一个参数是指数。在另一些示例中，参数可以是出砂函数的项的乘数。

出砂函数可以用于预测在发生故障之前材料的出砂体积，也就是在达到临界塑性应变之前在逐渐出砂期间材料的出砂体积。

例如，图4显示了出砂函数到从空心圆柱试验接收的数据的拟合结果。使用临界应变值以及指数(m)＝1的材料常数，对曲线的大规模出砂和逐渐出砂部分进行数据拟合。在另一个示例中，图6显示的是利用不同的指数(m)值，出砂函数对从射孔试验接收的数据的不同拟合的结果。

包含出砂函数的地质力学模型可以使用任何可应用的数据拟合方法拟合到出砂数据。例如，拟合可以用诸如线性回归和非线性回归之类的回归法执行。在另一个示例中，拟合可以通过求解偏微分方程组来执行，其中偏微分方程组包括包含了出砂函数的地质力学模型。偏微分方程组也可以包括油藏流动模型。偏微分方程组可以通过完全展开的雅可比行列式来耦合，求解偏微分方程组可以包括基于诸如出砂数据之类的接收数据在单个时间步长中同时求解完全展开的雅可比行列式。偏微分方程组也可以包括热模型。

在一个示例中，通过求解对地质力学油藏系统建模的偏微分方程组，可以获取包含出砂函数的地质力学模型对所接收的逐渐出砂数据的拟合，其中偏微分方程组包括油藏流动模型以及包含了出砂函数的地质力学模型，其中偏微分方程组通过完全展开的雅可比行列式耦合，且其中求解偏微分方程组包括基于所接收的逐渐出砂数据在单个时间步长中同时求解完全展开的雅可比行列式。

包含出砂函数的地质力学模型对出砂数据的拟合可以作为对地质力学油藏系统建模的方程组的更广泛计算的一部分来执行，所述方程组可用于计算在将流体注入油藏或是从油藏产生流体以及将应力施加到油藏的边界时出现的应力、应变和/或位移。在一些计算中，包含地质力学模型(包括出砂函数)、热模型和油藏流动模型的油藏系统模型能求解包含渗流、热流和地质力学的系统。

在步骤108，输出可以用于预测地质力学油藏系统的出砂的信息。此类信息可以是但不局限于一个或多个硬化参数中的至少一个值、临界塑性应变值和/或出砂函数的至少一个参数值。信息可输出至用户或各种组件，例如用户接口设备、计算机可读存储介质、监视器、用户可访问的本地计算机或是作为网络一部分的用户可访问计算机。举例来说，该输出可以用监视器或用户接口设备(例如包括个人数字助理(PDA)在内的手持图形用户界面(GUI))视觉显示给用户。

5.1.1另一些出砂预测系统和方法

在预测地质力学油藏系统的出砂中使用的另一些示例系统和方法包括将出砂判据应用于包含一个或多个硬化模型的地质力学模型的步骤。所述系统和方法包括以下步骤：接收指示与地质力学油藏系统内的材料相关的物理或化学属性的数据，以及通过求解对地质力学油藏系统建模的偏微分方程组来生成出砂预测。除了地质力学模型之外，偏微分方程组还可以包括地质力学油藏系统的油藏流动模型和/或热模型。

一个或多个出砂判据可被应用于地质力学模型。出砂判据可被确定成是何时达到(1)总应变不变量、(2)塑性应变不变量或(3)最大有效应力的临界值。当达到出砂判据的临界值时，材料将会发生故障，也就是破碎以形成砂且产生腔。以下在第5.5.2.5节论述出砂判据。

在一个示例中，偏微分方程组可以通过完全展开的雅可比行列式来耦合，其中求解偏微分方程组(例如通过使用计算机系统)包括基于接收数据在单个时间步长中同时求解完全展开的雅可比行列式。所产生的出砂预测可被输出到用户、用户接口设备、监视器、计算机可读存储介质、本地计算机或是作为网络一部分的计算机。例如，所产生的出砂预测可以用监视器或用户接口设备(例如包括个人数字助理(PDA)在内的手持图形用户界面(GUI))视觉显示给用户。

5.2用于操作地质力学油藏系统的系统和方法

还公开了用于在操作期间控制地质力学油藏系统的出砂中使用的系统和方法。操作地质力学油藏系统的方法根据的是这里公开的任一出砂预测系统和方法的实施结果。图7的流程图示出在基于出砂预测结果操作地质力学油藏系统中使用的示例系统和方法中的步骤。

在步骤700，接收指示与油藏系统内的材料相关的物理或化学属性的数据。步骤700中接收的数据可以包括以上在步骤100中描述的任何数据。

在步骤702中，基于包含出砂函数的地质力学模型对所接收的数据的拟合，计算至少一个操作参数的值，其中出砂函数预测地质力学油藏系统的出砂，且其中所述至少一个操作参数值指示地质力学油藏系统的出砂的条件。在一个示例中，操作参数的计算值指示地质力学油藏系统的稳定出砂的条件。

从这里公开的任一系统和方法得到的油藏出砂预测均可用于计算或得到用于操作油藏以实现期望的出砂量或期望的出砂行为的至少一个操作参数的值。操作参数的示例包括但不局限于生产压降、生产率、最小井底压力、产层温度、井眼中的流体流动压力、约束应力以及压差。出砂预测的实施结果可用来计算一个或多个操作参数的值，该值导致在油藏中产生和维持基本稳定的腔。如果油藏中的腔不随时间显著长大或者长大程度可以忽略，那么腔可以基本稳定。例如，一个或多个操作参数的值可以基于出砂预测结果得到，该值指示油藏的操作条件，在该条件下油藏产生有限的出砂量，且有限的出砂量产生的油藏中的腔实质上是稳定的。油藏出砂可以通过控制生产变量而被保持在最小，生产变量可以是生产压差、最小井底压力和生产率，但不局限于此。

在步骤704，根据至少一个操作参数的值来操作地质力学油藏系统。也就是说，地质力学油藏系统可以以至少一个操作参数的计算值操作，例如生产压差、生产率、最小井底压力、产层温度、井眼中的流体流动压力、约束应力、压差或是这些参数的任何组合的值。

任一出砂预测系统和方法的实施结果可用于确定能安装在油藏中以在井的整个寿命中实现期望量的出砂的完井技术类型。举个例子，油藏的井眼周围使用的屏障系统的设计(例如但不局限于砂筛技术或砾石充填)可以基于这些结果来选择。出砂预测还可以影响完井策略，例如裸眼井、套管井以及射孔套管井、使用带眼衬管或压裂充填(在水力压裂之后向裂缝中注入支撑剂)。此外，出砂预测可以用于确定在井的寿命中使用缓解出砂技术和安装缓解出砂设备的时点。

5.3关于建模方法的示例

图8和9示出在产油期间预测地质力学油藏系统的出砂中使用的系统的示例。可应用的建模系统包括描述地质力学油藏系统的各种物理方面的一个或多个模型。图8和图9示出包含油藏流体流动模型和地质力学模型的建模系统。油藏流体流动模型描述例如渗流、生产和注入。地质力学油藏模型描述例如在将流体注入到油藏中或者从油藏产生流体时以及将应力施加于油藏的边界时出现的应力、应变和位移。图9的系统还包括热模型。热模型描述热流。非线性偏微分方程组可用于关联这些模型的各个方面。

在接收指示地质力学油藏系统内的材料的物理或化学属性的数据之后(例如但不局限于塑性变形、逐渐出砂或大规模出砂)，解算器通过应用以上第5.1节中的方法的步骤且在相关点通过求解偏微分方程组来产生预测(例如出砂预测)。在图8和9的解算器中，偏微分方程组可以通过完全展开的雅可比行列式来耦合。求解偏微分方程组包括基于接收数据在单个时间步长中同时求解完全展开的雅可比行列式。

如图8和9所示，出砂预测步骤(上述步骤100到106)可以用非线性的偏微分方程组的解迭代执行。也就是说，在包含执行拟合的出砂预测步骤例如步骤102、104和106中的每个中，计算可以通过在单个时间步长中同时求解与特定步骤相关的代表地质力学油藏系统的方程的完全展开雅可比行列式来执行。所产生的出砂预测可输出到各种组件，例如输出到用户接口设备、计算机可读存储介质、监视器、用户可访问的本地计算机或是作为网络一部分的用户可访问计算机。

非线性的偏微分方程组包括与在出砂预测的给定步骤中分析地质力学油藏系统所使用的模型相对应的方程式。举例来说，图8提供一示例，其中非线性的偏微分方程组包括与地质力学油藏系统的油藏流动模型以及地质力学模型相对应的方程式。根据所执行的出砂预测步骤，地质力学模型可以包括一个或多个硬化模型、出砂模型(包含出砂函数)或是二者。作为另一示范，图9提供一示例，其中非线性的偏微分方程组包括与地质力学油藏系统的油藏流动模型、地质力学模型以及热模型相对应的方程式。在第5.5节提供与地质力学油藏系统的每一个不同模型相对应的方程式的示例。

模型的各个方面的耦合可以例如通过完全展开的雅可比行列式中的变量来实施。例如，通过以下变量中的一个或多个，完全展开的雅可比行列式可用来将油藏中的流体流动耦合到地质力学模型：有效应力、与地质力学模型相关联的孔隙度以及一个或多个位移。将地质力学模型耦合至流体流动的完全展开的雅可比行列式中的变量可以是与油藏流动模型相关联的孔隙度和渗透率。将热模型耦合到地质力学模型的完全展开的雅可比行列式中的变量可以是与热模型相关联的热应力。将热模型耦合到油藏流动模型的完全展开的雅可比行列式中的变量可以是与热模型相关联的油藏中的流体粘度、传导和对流。完全展开的雅可比行列式可以包括与耦合变量的变化率(即时间导数)、空间导数或部分空间导数相关的项，其中导数可以具有任何阶数，例如一阶导数、二阶导数、三阶导数等。耦合变量的一阶、二阶、三阶和/或更高阶的导数(无论是时间导数还是空间导数)可包含在完全展开的方程组中。在第5.6节提供可以耦合不同模型的变量的示例。

完全展开的雅可比行列式的非线性方程组可以通过数值方法来求解，例如在第5.7节中更详细论述的方法，其中非线性方程组例如利用所有解变量的完全Newton-Raphson展开来求解，这增强了解的稳定性且允许非线性迭代的二阶收敛速度。在第5.8节将会论述这里公开的不同方法的装置和计算机程序实施的示例。

在另一方面，系统和方法可以包括以下步骤：接收指示与地质力学油藏系统内的材料相关联的物理或化学属性的数据，定义包含多个网格单元的网格，以及通过求解对地质力学油藏系统建模的偏微分方程组来产生出砂预测。除了地质力学模型之外，偏微分方程组还可以包括地质力学油藏系统的油藏流动模型和/或热模型。

一个或多个出砂判据可被应用于地质力学模型。出砂判据可被确定为是何时达到(1)总应变不变量、(2)塑性应变不变量或(3)最大有效应力的临界值。在第5.5.2.5节论述出砂判据。

在使用判据(2)的计算，也就是包括对塑性应变不变量的临界值(临界塑性应变)进行计算的计算中，步骤100-106可被实施，出砂函数的至少一个参数可被计算。出砂函数可用于预测在达到临界塑性应变值之前网格单元的出砂量。

在一个示例中，偏微分方程组可以通过完全展开的雅可比行列式来耦合，其中求解偏微分方程组(例如通过使用计算机系统)包括基于接收数据在单个时间步长中同时求解完全展开的雅可比行列式。油藏模型、热模型和地质力学模型可以在三维网格单元或二维网格单元上计算。在第5.4节论述可用于这里的仿真方法中的三维和二维网格单元。

所产生的出砂预测可被输出到用户、用户接口设备、监视器、计算机可读存储介质、本地计算机或是作为网络一部分的计算机。

隐含地求解非线性方程组(例如使用解变量的完全Newton-Raphson展开)可以增强数值稳定性(例如当处理腔生成或者涉及非常小的网格块的任何仿真时)。使用隐含耦合的方程组的完全展开雅可比行列式为求解过程提供更高稳定性。

5.4仿真方法

图10A示出可用于地质力学模型以及油藏模型和/或热模型的计算的三维(3D)网格的示例。例如，在3D网格上可以计算多点通量(例如用于油藏和渗流)模型、地质力学模型(例如用于计算应力、位移、和/或腔生成(破碎))以及热模型中的一个或多个。地质力学模型以及油藏模型和/或热模型的计算可以使用并行处理方法来耦合网格。网格单元的尺寸可以是大约数英尺、数英寸或零点几英寸。在另一个示例中，网格单元的尺寸可以是大约数米、数厘米、数毫米、数微米或是零点几微米。

3D网格可以是包含六面体元素的结构化或非结构化的六面体网格。六面体网格单元具有八个角、十二个边缘(或边)以及六个面。六面体网格单元中的每一个可以包括至少八个节点(每个拐角处一个)或更多且多达二十七(27)个节点(也就是每个面的中心、每条边的中心、每个边缘的中心以及单元中心都有一个节点)。不同的六面体单元可以包括不同数量的节点。在一个示例中，3D网格可以包括结构化或非结构化的四面体网格元素。在另一个示例中，3D网格可以包括跨越了结构化或非结构化的四面体与六面体网格元素之间范围的其他形态的元素。3D网格可以包括上述网格元素的任何组合。

二维(2D)网格也可用于地质力学模型以及油藏模型和/或热模型的计算。例如，2D网格可用于轴对称计算。图10B示出可以用于地质力学模型、油藏模型和热模型的计算的二维(2D)网格的示例。

2D网格可以是包含四边形元素的结构化或非结构化四边形网格。每个四边形网格单元都具有四个角和四个边缘。四边形网格单元每个都包括至少四个节点(每个角处一个)以及上至五个节点(也就是一个节点处于中心)。

在某些示例中，计算可以在2D网格和3D网格二者上执行。举例来说，取决于油藏构形，某些计算可以在2D网格上执行，而另一些计算则可以在3D网格上执行。在这些计算中，2D网格的节点可被配置成与3D网格的外边界之一上的节点重合，诸如裂缝宽度和3D位移之类的某些计算可以在公共节点处耦合。2D网格的输入数据格式可以类似于3D网格的输入格式。

对于某些计算，参数诸如流体流动、位移、腔生成和牵引力可以跨网格单元的元素进行监视。

5.5模型

下面提供与地质力学油藏系统的每个不同模型相对应的微分方程的示例。用于计算中包含的模型的微分方程可被组合以便产生隐含的完全耦合公式。一致的单位组可用于计算中包含的方程式中的所有变量。

5.5.1油藏模型

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

$\frac{\∂}{\∂}\left(\mathrm{\ρ\φ}\right)=-\▿\left(\mathrm{\ρ}\stackrel{\→}{U}\right)+q_w,......\left(1\right)$

其中φ是孔隙度且ρ是可以作为压力的函数的流体密度。该模型允许在油藏元素中完井，且以上方程中的q_w解释了向油藏元素中的注入。

速度是与多孔材料相关的Darcy速度，并且可以定义为：

$\stackrel{\→}{U}=-\frac{K}{\mathrm{\μ}}(\▿p-\mathrm{\ρg}\▿h),......\left(2\right)$

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

流体流动方程中包含的地质力学变量强调的是流动与变形模型之间的耦合(以下在第5.5.2节描述一些地质力学项的定义)。

对于涉及温度变化的计算，温度相关的水属性可以以不同的方式进入。水属性可以作为压力(P)和温度(T)二者的函数进入以用于完全耦合计算。对于迭代耦合计算来说，水属性可以作为压力的函数进入，然后修正因子可以用于温度效应。在第5.5.3节更详细地说明温度相关属性的处理。

流体的热行为还可以通过使用修正因子修改流体属性来建模(以下在第5.5.3.3节中描述)。

5.5.1.1多相渗流

油藏模型考虑了从单相到复合物(从黑油到基于逸度)范围中的若干相行为模型。可对于含水相、非水液体相以及非水蒸汽相和nc成分来建模达西流动。任何相行为模型都可以与这里公开的渗流模型一起使用。流体流动方程可以按照通用的复合公式给出。用于多相流的代表成分质量平衡的偏微分方程是：

$\frac{\∂}{\∂t}\left(\mathrm{\φ}{N}_{\mathrm{ic}}\right)=-\underset{\mathrm{\α}}{\mathrm{\Σ}}\▿x_{\mathrm{ic}}^{\mathrm{\α}}{\mathrm{\ρ}}_{\mathrm{\α}}{\stackrel{\→}{U}}_{\mathrm{\α}}+q_{\mathrm{ic}},......\left(3\right)$

其中N_ic是由给出的每单位孔隙体积的成分ic的浓度，是处于相α的成分ic的克分子分数，ρ_α是相的摩尔密度，q_ic是每单位油藏体积的成分ic的摩尔流率。相α的速度由下式给出：

$U_{\mathrm{\α}}=-\frac{{\mathrm{KK}}_{\mathrm{r\α}}}{{\mathrm{\μ}}_{\mathrm{\α}}}\·(\▿P_{\mathrm{\α}}-{\mathrm{\γ}}_{\mathrm{\α}}\▿D),......\left(4\right)$

相压力可由下式定义：

P_α＝P+P_cα，......(5)

其中P_ca是毛细管压力，P是基准压力。基准压力用于PVT计算、井计算以及地质力学计算。基准压力可以是用于二相模型的非水相压力以及用于三相模型的非水液相压力。

渗流模型的孔隙度可被定义成：

φ＝φ_O[1+c_r(P-P_O)]......(6)

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

对于包含地质力学模型的计算来说，与初始未变形块体体积相关的孔隙度由方程15给出，其中可以看出，c_r可与孔隙度的Biot常数相关。

5.5.1.2非达西流动

在某些示例中，非达西流动可使用Forschheimer方程来建模以修改压力梯度与流体速度之间的关系。在另一些示例中，非达西流动可以通过规定流体速度与压力梯度之间的一般关系来建模。

对于涉及各向异性介质中的3-D多相流的系统来说，用于非达西速度的Forschheimer方程(其将替换以上的方程2)可以由下式给出：

$-\frac{k_{\mathrm{r\α}}}{{\mathrm{\μ}}_{\mathrm{\α}}}\underset{\‾}{K}\stackrel{\→}{\▿}{p}_{\mathrm{\α}}=(1+{\mathrm{\β}}_K{\mathrm{\ρ}}_{\mathrm{\α}}\frac{k_{\mathrm{ra}}}{{\mathrm{\μ}}_{\mathrm{\α}}}{\left|\right|{\stackrel{\→}{v}}_{\mathrm{\α}}\left|\right|}_2){\stackrel{\→}{v}}_{\mathrm{\α}},......\left(7\right)$

其中μ_a是相α的粘度，K是渗透张量，k_ra是相α的相对渗透度，v_α是相α的达西速度，表达式||...||₂是由定义的L₂范数，参数β_k是非达西系数。参数β_k可以在油藏中变化，由此，可以对每个网格块输入β_k的值。非达西系数β_k与变换常数的倒数相关，即β_k＝1/Г。参见Barree等人的“Beyond Beta Factors，A Complete Model forDarcy，Forchheimer，and Trans-Forschheimer Flow in PorousMedia，”SPE 89325，SPE annual Technical Conference andExhibition，Houston，TX，September 26-29，2004。用于相的雷诺数由以下方程给出：

$N_{\mathrm{re}}^{\mathrm{\α}}={\mathrm{\β}}_k{\mathrm{\ρ}}_{\mathrm{\α}}\frac{k_{\mathrm{r\α}}}{{\mathrm{\μ}}_{\mathrm{\α}}}{\left|\right|{\stackrel{\→}{v}}_{\mathrm{\α}}\left|\right|}_2......\left(8\right)$

方程式8中的项的单位应选择为使结果是无量纲的。在组合方程式7和8之后，Forschheimer方程式变成：

${\stackrel{\→}{v}}_{\mathrm{\α}}=-\frac{k_{\mathrm{r\α}}}{{\mathrm{\μ}}_{\mathrm{\α}}}{(1+N_{\mathrm{re}}^{\mathrm{\α}})}^{-1}\underset{\‾}{K}\stackrel{\→}{\▿}{p}_{\mathrm{\α}},......\left(9\right)$

其中非达西流动表达成相相关的渗透性修正函数，该函数随着用于该相的雷诺数而改变。这不同于标准的渗透性修正函数，因为该非达西公式对于每个相具有单独的值。

在一些情况中，Forschheimer方程式不提供对非达西流动的合适近似。如下形式的修正函数可以用于近似非达西流动：

${\stackrel{\→}{v}}_{\mathrm{\α}}=-\frac{k_{\mathrm{r\α}}}{{\mathrm{\μ}}_{\mathrm{\α}}}{f}_{\mathrm{\α}}\left(N_{\mathrm{re}}^{\mathrm{\α}}\right)\underset{\‾}{K}\stackrel{\→}{\▿}{p}_{\mathrm{\α}},......\left(9A\right)$

可以构造满足约束条件和f_α(0)＝1.0的修正函数。对于标准Forschheimer方程式来说，可以指定下面的函数：

$f_{\mathrm{\α}}\left(N_{\mathrm{re}}^{\mathrm{\α}}\right)={(1+N_{\mathrm{re}}^{\mathrm{\α}})}^{-1},......\left(9B\right)$

对于扩展Forschheimer方程式来说，可以指定下面的函数：

$f_{\mathrm{\α}}\left(N_{\mathrm{re}}^{\mathrm{\α}}\right)=K_{\mathrm{ratio}}+\frac{1+K_{\mathrm{ratio}}}{1+N_{\mathrm{re}}^{\mathrm{\α}}},......\left(9C\right)$

J.L.Miskimins等人的“Non-Darcy Flow in Hydraulic Fractures：Does It Really Matter？”SPE paper 96389，SPE annual TechnicalConference and Exhibition，Dallas，TX，October 9-12，2005中的方程式3是适用于本申请中公开的方法的非达西公式的另一个形式。

5.5.1.3油藏模型中的计算

油藏模型的计算可以在3D网格(该网格可以是用于地质力学模型的网格)上进行。3D油藏网格可以包括渗流计算。流体速度项可以对于油藏单元之间的流动以及对于油藏单元与腔单元(已发生出砂的单元)之间的流动来计算。渗流方程式的主变量可以是流体压力或流体成分，流体压力或流体成分可以在每个六面体元素的中心评估(基于单元)。在某些计算中，可以使用多点通量算法以用于非结构化油藏流动计算，因此当八个元素共享一个公共角时，所得到的计算模板可以是用于3D网格的通用六面体元素的27个点。

5.5.2地质力学模型

5.5.2.1多孔弹性材料

线性关系(小应变)可用于应变-位移相关性。涉及应力、应变、温度和孔隙压力的耦合流动/位移模型可以基于Biot多孔弹性理论。平衡方程可以基于总应力且假设准静态平衡。

多孔弹性方程式可以依据总应力、体应变、温度以及孔隙压力来公式化。总应力可以用平均牵引力定义，牵引力是在油藏的平面区段被观察到，其中所述平面区段包括由固体和来自流体的孔隙压力携带的负载。体应变可以是在将应变仪附着于变形的多孔材料时从应变仪观察到的应变。

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

${\mathrm{\ϵ}}_{\mathrm{ij}}=\frac{1}{2}{(u_{i,j}+u_{j,i})}_{\′}......\left(10\right)$

其中逗号表示微分，u_i是i方向的位移，ε_ij是多孔材料的体应变，扩展部分与正应变相对应。总应力满足平衡方程

σ_ij，j+f_i＝0.......(11)

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

在没有考虑温差时，关联总应力、体应变和孔隙压力的本构方程是：

${\mathrm{\σ}}_{\mathrm{ij}}={\mathrm{\σ}}_{\mathrm{ij}}^o+\frac{E}{1+v}({\mathrm{\ϵ}}_{\mathrm{ij}}+\frac{v}{1-2v}{\mathrm{\ϵ}}_{\mathrm{kk}}{\mathrm{\δ}}_{\mathrm{ij}})-\mathrm{\α}(p-p_o){\mathrm{\δ}}_{\mathrm{ij}},......\left(12\right)$

其中张力是正的，重复索引kk意味着总和，是初始原地应力，p_o是初始压力，E是弹性模量，v是泊松比，α是应力/应变方程式中的Biot常数，当i＝j时δ_ij是1，当i≠j时δ_ij是0。当时，可以假设应变为零。

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

${\mathrm{\σ}}_{\mathrm{ij}}={\mathrm{\σ}}_{\mathrm{ij}}^o+\frac{E}{1+v}({\mathrm{\ϵ}}_{\mathrm{ij}}+\frac{v}{1-2v}{\mathrm{\ϵ}}_{\mathrm{kk}}{\mathrm{\δ}}_{\mathrm{ij}})-\mathrm{\α}(p-p_o){\mathrm{\δ}}_{\mathrm{ij}}-{\mathrm{\α}}_{T}K(T-T_o){\mathrm{\δ}}_{\mathrm{ij}}$

，......(13)

其中α_T是应力/应变方程式的热体积膨胀系数，K是弹性体模量。压力p_o是初始孔隙压力，T_o是初始温度。

如果组合应力和压力项以形成那么该方程式变成标准热线性弹性本构方程，其中应力被有效应力取代。如果初始原地应力和初始孔隙压力为零，那么方程式采取标准形式：

${\mathrm{\σ}}_{\mathrm{ij}}^e=\frac{E}{1+v}({\mathrm{\ϵ}}_{\mathrm{ij}}+\frac{v}{1-2v}{\mathrm{\ϵ}}_{\mathrm{kk}}{\mathrm{\δ}}_{\mathrm{ij}})-{\mathrm{\α}}_{T}K(T-T_o){\mathrm{\δ}}_{\mathrm{ij}},......\left(14\right)$

孔隙度(与未变形的块体积相关)与应变以及流体压力(在不考虑温差时)之间的关系由下式给出：

$\mathrm{\φ}={\mathrm{\φ}}_o+{\mathrm{\α\ϵ}}_{\mathrm{ii}}+\frac{1}{M}(p-p_o),......\left(15\right)$

其中方程式15假设初始应变为零，φ₀是初始孔隙度，M^-1是孔隙度方程式中用于孔隙压力的Biot常数。

当考虑温差和沉淀分数时，与未变形块体积相关的孔隙度定义成：

$\mathrm{\φ}={\mathrm{\φ}}_o+\mathrm{\α}({\mathrm{\ϵ}}_v-{\mathrm{\ϵ}}_v^o)+\frac{1}{M}(P-P_o)-{\mathrm{\α}}_v(T-T_o)-\mathrm{\σ},......\left(16\right)$

其中α和M^-1是Biot常数，P是相压力(多相流动)，α_V是用于孔隙度的热体积膨胀系数，σ是沉淀分数(例如计算中网格元素的每块体体积沉积的固体废弃物的体积分数)。在废弃物移动穿过油藏时固定废弃物可沉积在孔内，随着废弃物沉积在孔隙空间中积累，孔隙度和渗透率可能下降。计算网格元素中的孔隙度可以是流体压力、温度及变形的函数，而由于固体沉积所导致的孔隙度减小量可设置成等于σ。

对于各向同性材料来说，在将地质力学方程式应用到地质力学油藏系统的建模之前，确定六个多孔弹性材料参数：E、v、α以及M^-1、α_T和α_V。

在某些示例中，油藏渗透度可表述成初始的方向渗透度(K_abs)乘以用于每个点的渗透率且用于每次步骤的渗透度乘数：K＝K_abs×f，其中f是一个或多个其他参数的函数，例如流体压力、总应力、块体积应变、孔隙压力、初始基准压力、主应力、有效塑性应变、当前孔隙度、起始孔隙度以及沉积分数。

在用于横向各向同性材料的本构方程中，应变可以用应力表达。在某些计算中，有效应力可以用于计算，初始原地应力可以是非零的。在另一些计算中，方程式可以用总应力且使用初始原地应力可以为零的假设来简化。关联横向各向同性材料的应力和应变的本构方程可表达成：

${\mathrm{\ϵ}}_{\mathrm{xx}}=\frac{1}{E_h}{\mathrm{\σ}}_{\mathrm{xx}}-\frac{v_h}{E_h}{\mathrm{\σ}}_{\mathrm{yy}}-\frac{v_v}{E_v}{\mathrm{\σ}}_{\mathrm{zz}}......\left(17\right)$

${\mathrm{\ϵ}}_{\mathrm{yy}}=\frac{1}{E_h}{\mathrm{\σ}}_{\mathrm{yy}}-\frac{v_h}{E_h}{\mathrm{\σ}}_{\mathrm{yy}}-\frac{v_v}{E_v}{\mathrm{\σ}}_{\mathrm{zz}}......\left(18\right)$

${\mathrm{\ϵ}}_{\mathrm{zz}}=\frac{1}{E_v}{\mathrm{\σ}}_{\mathrm{zz}}-\frac{v_v}{E_v}{\mathrm{\σ}}_{\mathrm{xx}}-\frac{v_v}{E_v}{\mathrm{\σ}}_{\mathrm{yy}}......\left(19\right)$

${\mathrm{\γ}}_{\mathrm{xy}}=\frac{2(1+v_h)}{E_h}{\mathrm{\σ}}_{\mathrm{xy}}......\left(20\right)$

${\mathrm{\γ}}_{\mathrm{xz}}=\frac{1}{G_v}{\mathrm{\σ}}_{\mathrm{xz}}......\left(21\right)$

${\mathrm{\γ}}_{\mathrm{yz}}=\frac{1}{G_v}{\mathrm{\σ}}_{\mathrm{yz}}......\left(22\right)$

方程式17-22假设对称轴是z方向且垂直方向是z方向。在执行涉及横向各向同性固体的计算之前，可以将方程式17-22中的五个弹性常数提供给模型。在执行多孔弹性计算时，可以假设Biot常数是各向同性且热膨胀系数是各向同性的，由此，在分析横向各向同性的多孔材料时，除了五个横向各向同性常数之外，还可以提供两个Biot常数和两个热膨胀系数。对于包含热模型(该热模型具有横向各向同性的弹性常数)的多孔弹性计算而言，方程式17-22中的应力可以是有效应力，方程式17-22中的应变可以是如下定义的有效应变：

${\mathrm{\ϵ}}_{\mathrm{ij}}^e={\mathrm{\ϵ}}_{\mathrm{ij}}-\frac{1}{3}{\mathrm{\α}}_T(T-T_o){\mathrm{\δ}}_{\mathrm{ij}}......\left(23\right)$

其中ε_ij是真实应变。

5.5.2.2多孔塑性材料

多孔塑性材料展现非线性行为，这是因为它可经历永久(即塑性)体积应变及由此的孔隙度变化。大的流体压力可导致多孔塑性材料屈服。结果，用于多孔塑性材料的地质力学计算可以预测流体孔隙度的大的突变。这些大的孔隙度突变能导致严重的稳定性问题，还产生负的流体压力。负压通常在流体压缩性低、渗透率低以及孔隙度扩张突然且大的时候出现。以上在第5.5.2.1节中论述的多孔弹性材料的方程式也可应用于多孔塑性材料。然而，孔隙度可被修改以顾及对多孔塑性材料可以预测的孔隙度变化。

在某些示例中，方程式可以用于抑制突然的孔隙度变化，以便改善计算的数值稳定性以及降低遭遇到负压的频率。在这些计算中，油藏模型中的孔隙度可定义成流体孔隙度(φ_fluid)，且可以被视为不同于地质力学模型中的孔隙度(φ_geomech)。在计算过程中，流体孔隙度与地质力学孔隙度之间的关系由下面的阻尼方程式支配：

流体孔隙度(φ_fluid)可用等式24计算且用于流体流动方程中，而地质力学孔隙度(φ_geomech)可以在地质力学模型中计算，τ是规定的时间常数。在(φ_geomech)的台阶变化之后，φ_fluid与φ_geomech之间的相对差可以在5τ的时段之后小于1％。可以选择τ值使得计算稳定，且可以将τ值选择为尽可能短的时间间隔，例如通过将τ设置成比计算中的时间步长或是计算的总时间的时间帧更短的值。举个例子，如果预期计算持续数天，那么τ可以是大约数分钟；如果预期计算持续几分钟，那么τ可以是大约数毫秒。对于许多多孔塑性材料的计算而言，τ值可以设置成大约一分钟。在另一些示例中，τ值可以设置为零(这去除了所有阻尼)。对本领域技术人员来说，用于给定计算的τ值是显而易见的。

5.5.2.3塑性流动

塑性流动也称为屈服，它表示材料的永久变形。一旦材料开始屈服，则塑性流动既可以持续，也可以不持续。屈服条件是标记从弹性变形到塑性变形的转变的数学表示。塑性流动方程中的假设是单个屈服条件和单个塑性势能用来描述油藏材料，且假设可以使用单个硬化参数来表示屈服面的移动。还可以假设应力和应变速率可以用下式表述：

${\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{ij}}=E_{\mathrm{ijkl}}({\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{kl}}-{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{kl}}^p)......\left(25\right)$

其中张量E_ijkl表示油藏材料的线性弹性属性：

E_ijkl＝λδ_ijδ_kl+μ(δ_ikδ_jl+δ_ilδ_jk).......(26)

其中λ和μ是Lame常数。方程式26可以用于各向同性弹性材料。在某些计算中，塑性模型有可以限制到各向同性弹性属性。在另一些计算中，方程式25和26可以包括对材料的塑性属性建模的其他项。

在某些情况中可以假设存在与屈服判据可以相等的塑性势G，但是塑性势G可以与非关联流动的屈服条件不同。此外，可以假设塑性应变速率由下式给出：

${\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p=\stackrel{\·}{\mathrm{\λ}}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ij}}}......\left(27\right)$

其中与Lame常数不相关的标量的值与可对材料属性设置的某些约束条件(即屈服条件)相关。塑性势G可以用于确定塑性应变速率的方向，而标量可用于确定其大小。

有效塑性应变(ε^p)可定义成：

${\mathrm{\ϵ}}^p={\∫\stackrel{\·}{\mathrm{\ϵ}}}^p\mathrm{dt}=\∫\sqrt{{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p}\mathrm{dt}......\left(28\right)$

有效塑性应变的变化速率与参数之间的关系是：

${\stackrel{\·}{\mathrm{\ϵ}}}^p=\stackrel{\·}{\mathrm{\λ}}\sqrt{\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ij}}}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ij}}}}......\left(29\right)$

如果将单位张量n_ij定义成：

$n_{\mathrm{ij}}=\frac{\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ij}}}}{\sqrt{\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{kl}}}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{kl}}}}}......\left(30\right)$

那么可以将塑性应变速率张量写为：

${\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p={\stackrel{\·}{\mathrm{\ϵ}}}^p{n}_{\mathrm{ij}}......\left(31\right)$

可以假设屈服条件的简化形式可以用应力(σ_ij)和硬化参数(κ)写成：

F(σ_ij，κ)≤0.......(32)

其中F在弹性状态中为负，在塑性状态中为零，并且κ是塑性应变(应变-硬化)的函数。在塑性变形过程中，屈服条件满足以下关系：

$\frac{\∂F}{{\∂\mathrm{\σ}}_{\mathrm{ij}}}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{ij}}+\frac{\∂F}{\∂\mathrm{\κ}}\frac{\∂\mathrm{\κ}}{{\∂\mathrm{\ϵ}}_{\mathrm{ij}}^p}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p=0......\left(33\right)$

方程式33可与方程式25、26和27组合以得到以下表示式：

${\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{ij}}=\{E_{\mathrm{ijkl}}-\frac{E_{\mathrm{ijpq}}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{pq}}}\frac{\∂F}{{\∂\mathrm{\σ}}_{\mathrm{mn}}}{E}_{\mathrm{mnkl}}}{\frac{\∂F}{{\∂\mathrm{\σ}}_{\mathrm{rs}}}{E}_{\mathrm{rsab}}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ab}}}-\frac{\∂F}{\∂\mathrm{\κ}}\frac{\∂\mathrm{\κ}}{{\∂\mathrm{\ϵ}}_{\mathrm{ab}}^p}\frac{\∂G}{{\∂\mathrm{\σ}}_{\mathrm{ab}}}}\}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{kl}}......\left(34\right)$

本构方程将应力速率关联到应变速率，其中系数取决于弹性属性、当前应力、当前塑性应变以及硬化参数。对于关联类型的流动计算来说，可以假设F＝G，这使本构方程是对称的。如果非线性流动方程式用隐含耦合计算来求解，那么即便本构方程是对称的，用于方程组的雅可比行列式也可以略微非对称。

5.5.2.4硬化模型

对材料的后续塑性变形要求以及应力-应变关系的检查可以提供关于塑性流动可以还是不可以持续的指示。材料可以是理想塑性的或经历应变-硬化。理想塑性材料(例如但不限于结构钢)展现出不被塑性变形改变的屈服条件。然而，很多材料被非弹性形变(称为应变-硬化)改变，屈服条件可随着材料对屈服的抵抗性增大而改变。

地质力学包括硬化模型，硬化模型可用于对地质力学油藏系统内的材料的塑性变形建模。硬化模型的示例包括：具有切变硬化或帽盖硬化的Drucker-Prager模型，具有板状(tabular)硬化的修正Drucker-Prager模型，修正Bigoni-Piccolroaz模型以及Matsuoka-Nakai模型。

具有切变硬化的Drucker-Prager

具有切变硬化和正拉伸应力的Drucker-Prager方程式的屈服条件可表述成：

$F=\sqrt{J_2}+\mathrm{\α}{I}_1-k{\left({\mathrm{\ϵ}}^p\right)}^m-\mathrm{\Γ}=0......\left(35\right)$

其中所有标量参数是非负常数，α是常数，m是指数，Г是屈服常数。参数α、k、m以及Г的值对应于可被输入以执行计算的输入参数。常数α可以采用0.0与1√3之间的值。如果α值为零，那么Drucker-Pragermodel变成Von Mises塑性模型。参数α与用于Mohr-Coulomb模型的摩擦角相关。例如，有效塑性应变ε^P是可被计算并报告给例如用户的硬化参数。

在硬化也就是ε^P增大之后，如图11A所示，屈服面可视为从初始表面位置移动到最终表面位置。在该计算中，第一应力不变量I₁(总应力的第一不变量)关于压缩(I₁＝σ_kk，关于k求和)是负的。

用于Drucker-Prager模型的塑性流动方程可表示成：

${\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p=\stackrel{\·}{\mathrm{\λ}}\{\mathrm{\α}{\mathrm{\δ}}_{\mathrm{ij}}+\frac{S_{\mathrm{ij}}}{2\sqrt{J_2}}\}......\left(36\right)$

其中并且如果假设的是相关联的流动，那么方程式35中的α与屈服方程式中使用的α相同。

具有帽盖硬化的Drucker-Prager

具有帽盖硬化的Drucker-Prager模型可以具有两个屈服面。一个屈服面是由下式给出的非硬化Drucker-Prager破坏面(完全塑性)：

$F_s=\sqrt{J_2}+\mathrm{\α}{I}_1-\mathrm{\Γ}=0......\left(37\right)$

其中所有标量参数都是非负常数且拉伸应力是正的，α是参数，Г是屈服常数。第二屈服面是如下形式的椭圆硬化帽盖：

$F_c=J_2-\frac{1}{R^2}[{(X-L)}^2-{(I_1-L)}^2]=0......\left(38\right)$

变量X≤RГ是塑性体积应变的函数(压缩为负)，并且是位于椭圆帽盖与I₁轴交叉处的第一应力不变量I₁的值。变量L≤0也是塑性体积应变的函数，并且是椭圆帽盖与Drucker-Prager表面F_s交叉处的l₁值。椭圆帽盖在与l₁轴交叉处是垂直的，并且在与Drucker-Prager破坏面交叉处是水平的。以下方程式可以用于将X关联到有效塑性体积应变以及强制盖帽在Drucker-Prager屈服面处有零斜率的约束条件：

$X=X_o+\frac{1}{D}\ln (1+\frac{{\stackrel{\‾}{\mathrm{\ϵ}}}_v^p}{W})$ 或 $X=X_o-H(-{\stackrel{\‾}{\mathrm{\ϵ}}}_v^p)......\left(39\right)$

以及

$\frac{X-L}{R}=\mathrm{\αL}-\mathrm{\Γ}......\left(40\right)$

从方程式39和40可以得到两个帽盖硬化模型。第一帽盖硬化模型使用方程式39中的对数用于帽盖处的硬化。第二帽盖硬化模型使用方程式39的列表函数H()，其中H()是严格单调递增的，并且H(0)＝0。在方程式39和40中，有效塑性体积应变以及变量L和X全都为零或负值。对于X的初始大小超出X₀的大小时的计算，的初始值可以为非零。

在图11B中显示了屈服面的移动，其中实线是初始表面，虚线是硬化之后的表面位置，L₁和L₂是初始表面和后面的屈服面的L值。

使用塑性方程式可以显示，塑性体积应变速率可以用来自方程式的塑性参数表达为下面的形式：

${\stackrel{\·}{\mathrm{\ϵ}}}_v^p=\frac{6}{R^2}(I_1-L)\stackrel{\·}{\mathrm{\λ}}......\left(41\right)$

其中不论何时在帽盖上发生塑性变形，都有I₁≤L≤0，并且以及可计算并报告给例如用户的硬化参数是L、以及剪切面上的塑性流动方程式可以与先前对于标准Drucker-Prager模型所描述的方程式相同，而帽盖上的流动可以用下式给出：

${\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p=\stackrel{\·}{\mathrm{\λ}}\{\frac{2}{R^2}(I_1-L){\mathrm{\δ}}_{\mathrm{ij}}+S_{\mathrm{ij}}\}......\left(42\right)$

具有板状硬化的修正Drucker-Prager

用于具有剪切硬化和正拉伸应力的修正Drucker-Prager方程式的屈服条件是

$F=\frac{1}{2}\sqrt{J_2}[1+\frac{1}{K}-(1-\frac{1}{K})\cos \left(3\mathrm{\θ}\right)]+\mathrm{\α}{I}_1-H\left({\mathrm{\ϵ}}^p\right)-\mathrm{\Γ}=0$

.......(43)

其中α是常数，K是偏差，Г是屈服常数，θ对于单轴压缩试验来说是π/3并且可定义成：

$\mathrm{\θ}=\frac{1}{3}{\cos }^{-1}\left(\frac{3\sqrt{3}}{2}\frac{J_3}{J_2^{3/2}}\right)......\left(44\right)$

H()是有效塑性应变的列表函数，l₁是第一应力不变量，J₂和J₃是总应力偏量的第二和第三不变量：

$J_2=\frac{1}{2}{S}_{\mathrm{ij}}{S}_{\mathrm{ij}}......\left(45\right)$

$J_3=\frac{1}{3}{S}_{\mathrm{ij}}{S}_{\mathrm{ij}}{S}_{\mathrm{ki}}......\left(46\right)$

$s_{\mathrm{ij}}={\mathrm{\σ}}_{\mathrm{ij}}-\frac{1}{3}{\mathrm{\σ}}_{\mathrm{kk}}{\mathrm{\δ}}_{\mathrm{ij}}......\left(47\right)$

其中S_ij是偏应力(减去平均应力)，δ_ij是克罗内克函数(kroneckerdelta)(在i＝j时为1，否则为0)，J₂在i和j上的求，且J₃在i、j和k上求和。如果将材料参数K设置成1.0并且如果H()的表值满足H(ε^p)＝k(ε^p)m，那么修正Drucker-Prager模型还原到Drucker-Prager模型。参数K的值可以在从0.78到1.0的范围。有效塑性应变ε^P的硬化参数可被报告给例如用户。常数α，K和Г可以被输入。

在图12A中示出了参数K的影响，图12A示出对于八面体平面中的四个K值(且对于恒定的I₁)，八面体(偏的)平面中修正Drucker-Prager屈服面的曲线图。在图12A的四个图中，可以假设压缩是正的，并且都被归一化以与正垂直轴上的1.0相交。标准Drucker-Prager具有K＝1.0且在平面中是圆形。对小于0.78的K值来说，修正Drucker-Prager表面不再凸起。某些计算可被限制从而K值大于或等于0.78。

用于修正Drucker-Prager方程式的塑性流动方程是：

${\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij}}^p\stackrel{\·}{\mathrm{\λ}}\{{\mathrm{\α}}_f{\mathrm{\δ}}_{\mathrm{ij}}+\frac{S_{\mathrm{ij}}}{4\sqrt{J_2}}[1+\frac{1}{K}-(1-\frac{1}{K})\cos \left(3\mathrm{\θ}\right)]+\frac{3\sqrt{3}}{4}(1-\frac{1}{K})T_{\mathrm{ij}}\}$

.......(48)

其中T_ij由下式定义：

$T_{\mathrm{ij}}=\frac{2}{3}{\mathrm{\δ}}_{\mathrm{ij}}-\frac{S_{\mathrm{ik}}{S}_{\mathrm{kj}}}{J_2}+\frac{\cos \left(3\mathrm{\θ}\right)}{\sqrt{3J_2}}{S}_{\mathrm{ij}}......\left(49\right)$

方程中的α_f与在假设相关流动的屈服方程中使用的相同。如果假设T_ii＝0以及S_ii＝0，那么在方程式49中将α_f设置成等于零将可以导致用于计算的塑性体积应变为零。

具有板状硬化(tabular hardening)的修正Bigoni-Piccolroaz

修正Bigoni-Piccolroaz模型(MBP)是对D.Bigoni等人的“Yieldcriteria for quasi-brittle frictional materials，”Intl.J.of Structures 41：2855-2878中的屈服判据的修正。用于具有正拉伸应力的MBP模型的屈服方程可表达成：

$A^{-1}=\cos \{\frac{1}{3}{\cos }^{-1}(-\mathrm{\γ})-\mathrm{\β}\frac{\mathrm{\π}}{6}\}......\left(51\right)$

其中α是常数，Г是屈服常数，γ是偏差，且其中0≤β≤1，Г≥0，0≤γ≤1是材料常数，H()是有效塑性应变的列表函数。参数J₂和J₃是应力偏量的不变量，I₁是第一应力不变量。可选择常数A使得对于单轴压缩测试是1.0。MBP模型使用应变-硬化来对塑性硬化建模，并且可以假设屈服面在硬化时均匀扩张(没有旋转)。可以是输入参数的常数α、β、Г和γ控制主应力空间中屈服面的形状。主应力空间中的屈服锥体角度由α控制，而I₁轴上的锥体顶端位置在硬化之前由Г/α给出。八面体平面中的屈服面形状(对于常数I₁)可以由参数β和γ控制。就MBP参数的具体选择而言，若干个公共屈服面可被再现。公共屈服面(其是MBP屈服面的特例)是von Mises、Drucker-Prager、Mohr-Coulomb、Lade、Tresca以及Matsuoka-Nakai。

参数β和γ可被修改以确定八面体平面中的屈服面形状。图12B中显示了不同的β和γ值得到的屈服面。在图12B的曲线图中，可假设压缩是正的，所有曲线都被归一化成交叉正垂直轴上的1.0。MBP模型屈服方程未显现出修正Drucker-Prager模型屈服方程可能出现的凸性损失。举个例子，对于β＝0.0且γ＝1.0的情况，MBP模型屈服方程在八面体平面中产生三角形(相比于K＝0.60时的修正Drucker-Prager模型屈服方程的非凸表面)。

用于MBP模型的塑性流动方程是

该流动方程中的α_f与假设了关联流动的屈服方程中使用的α相同。如果在方程52中将α_f设置成等于零，那么对于计算而言塑性体积应变是零。

5.5.2.5出砂模型

来自用于塑性变形计算的硬化模型的一个或多个参数可被用于出砂模型计算。在使用出砂模型的计算中，至少一个出砂判据可被应用于一个或多个网格单元。可以假设，当在网格单元中心处的以下各项达到临界值时，计算网格单元失效：

(1)总应变不变量，或

(2)塑性应变不变量，或

(3)最大有效应力。

对于这些出砂判据中的每个来说，网格单元失效表示油藏材料的出砂(导致腔)。

对判据(1)来说，总应变不变量可表达成：

$\mathrm{\ϵ}=\sqrt{{\mathrm{\ϵ}}_{\mathrm{ij}}{\mathrm{\ϵ}}_{\mathrm{ij}}}.....\left(53\right)$

对判据(2)来说，塑性应变不变量可表达成：

${\mathrm{\ϵ}}^p=\sqrt{{\mathrm{\ϵ}}_{\mathrm{ij}}^p{\mathrm{\ϵ}}_{\mathrm{ij}}^p}......\left(54\right)$

在某些计算中，如果硬化模型是Drucker-Prager帽盖硬化模型，则塑性应变不变量判据可以不被应用。

对判据(3)最大有效应力而言，最大主有效应力可被计算且其数值可以与输入值(假设拉伸应力是正的)相比较。

判据(2)塑性应变不变量可以用于说明在单元完全失效之前来自单元的暂时出砂。如果用于完全单元失效的临界塑性应变是那么在单元完全坍塌之前单元的出砂分数可以由出砂函数表示。出砂函数可以是函数f(x)，其中当x＝0时，f(x)＝0，当x＝1时，f(x)＝1，x是临界塑性应变的函数函数f(x)可以是x的任何单调函数，包括但不限于分数函数、幂函数、正弦函数、余弦函数、对数函数、指数函数、S函数或它们的任何组合。在一些示例中，x可以是临界塑性应变和塑性应变不变量(ε^p)二者的函数，例如但不局限于比值在一个示例中，出砂函数可以由下式给出：

${\left(\frac{{\mathrm{\ϵ}}^p}{{\mathrm{\ϵ}}_{\lim }^p}\right)}^m.....\left(55\right)$

其中m是指数。

5.5.2.6地质力学模型中的计算

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

5.5.3热模型

热模型计算当热或冷流体注入到油藏中或者从油藏中产生时发生的温度变化，还计算来自井眼的传导，所述传导可以循环热/冷流体，而实际上不注入或产生流体。温度计算可包括渗流和地质力学，但可配置成不包括蒸汽注入。注入的水可假设为液相。

热模型计算由于油藏中的传导和对流、井上热或冷流体的注入/产生、井的传导以及流体与固体之间的热和机械相互作用而出现的温度变化。这些机制可以在与渗流方程和固体形变方程一起求解的单个能量方程中组合。能量方程可以依照用于多孔固体的Lagrangian方法来公式化，全部流体移动可以与多孔固体的移动相关。

对于此类公式来说，在评估油藏元素的能量平衡时，多孔固体的质量可以是恒定的，而流体质量随着流体流入和流出多孔固体而改变。

5.5.3.1组合累积项

用于流体和固体中的能量变化的组合表达(当计算中包含地质力学系统时)可以用下式给出：

.......(56)

其中是关于未变形块体容积的孔隙度，α是相，ρ_α是相密度，S_α是相α的饱和度，h_α是相α的比焓，α_T是应力/应变方程中的体积膨胀系数，K是渗透率，α_v是孔隙度方程中的体积膨胀系数，T₀是流体温度。传导和对流项可影响组合累积项。在方程56中，可以近似为相压力相同，项可以用来近似，C_r项可以是热容C_v是以恒定孔隙体积和恒定块体体积测量的多孔固体的热容。可以看出，如上定义的热容C_r是在恒定流体压力和恒定块体体积下的热容。方程56是用于元素中内能变化总速率的通式，并且当在计算中考虑了包含流体温度T₀的项且计算包含热和地质力学模型二者时，方程56用于累积项。

如果计算中考虑流体温度T₀，那么方程56采取更简单的形式。

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

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

而对于包含热模型但不包含地质力学模型的计算，且流体内能通过流体焓来近似的情况而言，表达式是：

5.5.3.2能量方程

以上对累积项的选择确定了最终能量守恒方程中包含的方程。能量方程可表达成：

.......(60)

其中在注入井可以应用恒定温度边界条件。由于方程60不包含输送之外的任何力学项，因此当焓仅是温度的函数时，网格块中的温度不应受到流体相或固体的膨胀或压缩的影响。然而，由于潜热可包含在液相组分的焓与其在汽相中的对应焓之间的差别中，因此，方程60仍考虑了从液体到蒸汽的汽化化热。

如果仿真包括地质力学计算且在计算中考虑了包含流体温度T₀的项，那么通用能量方程可以用下式给出：

项存在于能量方程两端，也就是存在于累积项中且还作为对块固体(bulk solid)做的功，由此方程61中并未包含

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

有若干选项用于计算流体的温度相关密度和粘度，可用于计算密度和粘度的方法根据所使用的相行为模型及数值技术而变化。在某些示例中，流体属性可以直接作为温度和压力的函数来计算。在某些其他示例中，流体属性可以作为当前单元压力和初始单元温度(等温闪蒸(isothermal flash))的函数来计算，并且修正函数可用于校正初始单元温度与当前单元温度之间的差。

对于隐含单相运行(run)、隐含双相运行以及隐含复合运行，流体属性可以直接作为温度和压力的函数来计算。包含地质力学模型的这些计算可被迭代耦合或完全耦合。

在黑油、K值和复合PVT程序包中计算的粘度和密度可被修正以说明由于温度改变引起的流体属性变化。计算期间的当前流体压力和初始油藏温度可以在每个PVT程序包(等温闪蒸)中用于计算流体属性，然后，这些属性通过乘以用户提供的修正因子而被修正，即其中可以作为表格被输入，且μ(p，T_o)可在PVT程序包中被计算。

对于许多热-地质力学研究，关注主要集中在温度如何影响流体粘度和热应力。对此类应用来说，修正因子或直接计算流体属性都可使用，使用这两种技术的结果可非常相似。然而，对于温度强烈影响密度和/或相分离的应用而言，则流体属性的直接计算可以用作温度和压力的函数。

5.6各种模型的相互关系

一种系统和方法可以耦合不同地质力学油藏系统模型之间的相互作用，例如渗流、地质力学位移和应力、出砂以及腔形成之间。例如，地质力学解可通过孔隙度项和渗透率项来影响渗流计算。作为另一个示例，油藏中的流体解可以通过其在有效应力中的作用来影响地质力学计算。在又一个示例中，流体解可以通过在腔的面处发生的法向牵引力来影响包含出砂模型的地质力学计算，且可以影响油藏中的流动。

5.7出砂预测的实施

5.7.1包含出砂模型的计算

图13示出在计算的一时间步长中在出砂预测的示例计算中的若干考虑，在该时间步长中腔已经产生在油藏履盖层(覆盖井中的关注区域的岩层)与下伏岩层(处于井中的关注区域下方的岩层)之间。在每个时间步长末端，可以检查失效网格单元，一旦单元被识别成故障单元，则它可以立即从变形计算(例如包含地质力学模型的计算)移除。失效单元是达到出砂判据的单元(参见以上的第5.5.2.5节)，其可用于指示材料发生故障的点，也就是破碎以形成砂并产生腔的点。一旦单元失效且成为腔的一部分，则它可被从用于计算位移的刚度矩阵排除。如图13所示，在每个时间步长中，可对已有腔与未失效单元之间界面前头的未失效油藏单元进行失效检查。亦如图13所示，已有腔与未失效单元之间界面处的牵引力条件可以取为等于界面处腔单元中的流体压力。标准渗流模型和地质力学模型可以应用于未失效的完好网格单元。

然后失效单元中的砂量可被添加到最近采油井的报告出砂量。在一些计算中，单元的砂量可以是单元的初始块体体积。失效单元可以保留在流体流动计算中，以便允许所产生的腔的壁与井之间的流动，并且单元的孔隙度可以保持恒定。在一些计算中，腔中渗透率可增大以最小化井与腔壁之间的压降(用于流体流动)。一旦单元从变形计算移除，那么它的属性不再影响变形，除了与腔直接相邻的失效单元可对未失效的多孔固体(即未失效单元)施加压负载之外，并且腔壁处的压负载可以直接关联到与腔壁相邻的失效单元中的流体压力。

腔生成计算可以在物理上和/或数值上是不稳定的。为了避免物理不稳定性，在几乎完美弹性的塑性曲线的末端可以添加一些硬化，并且可以为不稳定腔的生长设置时间帧(也就是说，腔的生长速率可以直接关联到时间标尺设置)，作为非限制性示例，可以是一分钟。为了将数值不稳定性减至最小，模型的若干参数可设置成确定值。例如，方程24的时间常数τ(其将流体孔隙度变化关联到地质力学孔隙度变化)可设置成非零值以最小化不稳定性，因为零值去除了所有阻尼。较高的时间常数τ值可以增强数值稳定性。作为另一个示例，为失效的网格单元分配的渗透率(指示出砂)可设置成最小化数值不稳定性的值。作为非限制性示例，渗透率可设置成大约1000达西的值。可为渗透率设置更低或更高的值以增强数值稳定性。

5.7.2偏微分方程组的解

用于计算中包含的模型的整个方程组的雅可比行列式可被完全展开，并且所有方程可以在线性解算器中同时求解。各种机制可以用单向、显性、迭代和完全耦合技术中的任一种来耦合。最稳定的耦合技术可用于出砂预测。

以上章节中论述的方程可以在单个程序中实施，程序计算全部方程的隐含解，也就是说，反向尤拉技术可用于近似时间导数，方程中的全部主变量和系数可以在时间步长末端被评估。主变量的非限制性示例是油藏中的流体压力、油藏中和油藏边界处的位移(如果无约束)、裂缝中的流体压力(或者裂缝中的流体量，如果部分饱和的话)、井压及至少一个生产函数参数。方程组可以使用Newton-Raphson技术求解，其中为整个方程组产生完全展开的雅可比行列式，并且增量校正利用包含所有解变量的线性解算器来发现。

可以通过计算到失效单元的流动以及计算邻接的油藏中出现的有效应力和位移，可以对腔的生成进行仿真。然后，出砂判据与腔界面处的应力和位移状态以及牵引力条件相结合，以便确定腔是否发展。在单个时间步长期间，腔可以跨多个网格单元发展。每个时间步长可以用前一时间步长得到的腔配置、地质力学状态、流体状态(如果适用)以及热状态(如果适用)开始。然后，假设腔不扩展，计算将会迭代至质量守恒和应力平衡方程的收敛解。一旦获取了新的收敛解，则可以检查出砂判据来了解在任何单元中是否达到出砂判据的临界值。这些单元可以如第5.7.1节中描述的那样被处理。该迭代序列可以在方程以及出砂判据检查上继续进行多次，或者一直到在移到下一时间步长之前对于一个时间步长没有看到腔的进一步发展。

5.8关于装置和计算机程序实施的示例

这里公开的方法可以利用依照以下程序和方法的装置例如计算机系统来实施，诸如在本章节中描述的计算机系统。这样的计算机系统还可以存储和操作数据，该数据指示：与地质力学油藏系统相关联的物理属性、用于计算中所包含的模型的全部方程组的完全展开的雅可比行列式、完全展开的雅可比行列式的解、所产生的裂缝预测、或可供使用这里描述的分析方法实施的计算机系统使用的量度。所述系统和方法可以在各种类型的计算机架构上实施，例如在单个通用计算机、并行处理计算机系统、工作站、或联网系统(例如图14所示的客户机-服务器构造)上实施。

如图14所示，用于实施这里公开的一个或多个方法和系统的建模计算机系统可以链接到网络链，网络链可以是例如与其他本地计算机系统相连的局域网(“LAN”)的一部分，和/或连接到其他远程计算机系统的诸如因特网之类的广域网(“WAN”)的一部分。

建模系统包括这里描述的任何方法。例如，软件组件可以包括促使一个或多个处理器实施以下步骤的程序：接受指示与地质力学油藏系统相关联的物理属性的多个参数，和/或所产生的裂缝预测；以及在存储器中存储指示与地质力学油藏系统相关联的物理属性的参数和/或所产生的裂缝预测。举例来说，系统可以接受命令，以便接收指示与地质力学油藏系统相关联的物理属性的参数和/或用户手动输入(例如经由用户接口)的所产生的裂缝预测的参数。程序可以促使系统从数据存储(例如数据库)检索指示与地质力学油藏系统相关联的物理属性的参数和/或所产生的裂缝预测的参数。这种数据存储可以保存在大容量存储器(例如硬盘驱动器)或其他计算机可读介质上，并且可被载入计算机记忆体，或者计算机系统也可以经由网络来访问所述数据存储。

6.结果

6.1出砂预测试验

出砂试验在盐水和煤油饱和的岩芯样品上进行。出砂预测应用于出砂试验结果。未考虑温度影响。假设单相仿真可用于仿真油藏系统。

6.2数据采集

下面论述网格的几何形状、材料特性信息、边界条件以及渗透率函数。

6.2.1几何形状/网格

试验样品长约100毫米，且具有约100毫米的圆直径。中心孔的直径约为20毫米。在表1中给出了盐水和煤油试验样品的确切尺寸。图15示出用于盐水试验的岩芯样品，包括样品尺寸。在试验过程中使用了末端压板来施加轴向负载。还应用了空心圆柱试验(如图3所示)。

基于表1给出的尺寸建立单元的轴对称网格以用于计算。橡胶套未表示在计算网格单元中。但是，仿真模型设置成包含两个压板，压板的半径尺寸与样品的半径尺寸匹配，压板的轴向尺寸是20毫米。

表1：样品尺寸

试验	长度	直径	内部孔直径
					(mm)	(mm)	(mm)
盐水	102.72	99.77	19.56
				煤油	105.25	99.95	19.47

6.2.2材料特性

材料的物理常数用来自三轴试验的数据确定。所测量的材料响应和利用模型计算的材料响应对于盐水饱和芯(在图16中显示)和煤油饱和芯(在图17中显示)二者来说都接近一致。表2中给出了用于对盐水饱和样品建模的物理常数。

表2：盐水饱和样品的物理属性

变量	值
		孔隙率	0.296
泊松比	0.03
		杨氏模量	218823psi
岩石密度	2.65g/cc
		屈服常数Г	43.2
阿尔法	0.36581616

6.2.3边界条件

仿真通过求解完全耦合的地质力学模型和流体流动模型来进行，由此需要设置流体流动和地质力学边界条件(图18所示)。对地质力学边界条件来说，恒定应力被指定给空心圆柱的顶部、外表明和内表面(边界)(用图18的面板(A)中的箭头显示)。地质力学模型是在底部位移受限以防止数值模型移动和旋转(也就是限制整个表面以不在z方向上移动)。为了仿真试验流体流动条件(见图18的面板(B))，在空心圆柱上应用压降，圆柱外侧的压力大于圆柱内侧的压力。在圆柱顶部和底部没有应用流体流动边界条件。间隔物的孔隙度设置成零，这实际上使间隔物的渗透率也等于零。

空心圆柱内侧上的压力和空心圆柱内表面上的应力保持等于0psi。仿真使用绝对压力值，因此零压力值被设置成14.7psi以考虑大气压力。圆柱外侧上的应力和外侧上的压力二者都根据在对盐水饱和样品执行的试验中应用的流动和应力体系而改变。图19示出对于盐水饱和样品，测量且离散化的(仿真)约束应力、流速和流动压力与时间的关系。

6.2.4渗透率函数

使用试验测量，计算径向各向同性均质构型的渗透率。用于芯样品构型的径向流动的Darcy定律可表达成：

$k=\frac{\mathrm{Q\μ}\ln \left(\frac{r_e}{r_i}\right)}{2\mathrm{\πh\ΔP}}......\left(62\right)$

其中k是渗透率，Q是流速，μ是粘度，r_i是内半径，r_e是外半径，h是高度，ΔP是压降。渗透率用芯样品的尺寸以及方程62确定。该计算产生的渗透率分布显示，作为流速的函数，渗透率不是恒定的(如图20所示)。特别地，当流速增大时，渗透率减小(观察到受压力测量的精度影响)。平均渗透率通过将幂律曲线拟合到作为约束应力的函数的计算渗透率的曲线来确定(图21)。该幂律函数用于计算实现流速的压力。换句话说，给定幂律拟合推导的渗透率函数，改变压力从而获取期望的流速。

通常，渗透率是单元属性的函数而不是诸如约束应力之类的边界条件的函数。假设平均应力可用作约束应力的替代值，则平均应力用下式计算：

${\mathrm{\σ}}_{\mathrm{mean}}=\frac{1}{3}({\mathrm{\σ}}_1+{\mathrm{\σ}}_2+{\mathrm{\σ}}_3)......\left(62\right)$

其中σ_mean是平均应力，σ₁到σ₃是主应力。该假设使得计算流速与测量流速不相同。试验中测得的流速与通过数值仿真预测的流速比较。校正因子应用到外边界压力以使流速匹配。所计算的压力乘以1.035以获得流速匹配(如图22所示)，这意味着仅使用3.5％的校正因子。对于没有出砂的样品获得该流速匹配。如以下的第6.3节所述，流速依赖于出砂量而变化。

6.3分析

下面论述所计算的砂匹配及其对某些参数的敏感度。还论述是否可创建可维持的稳定腔的问题。

6.3.1出砂模型计算

用于出砂模型的以下输入将被计算：

1.临界塑性应变值——导致岩层破坏的有效塑性应变的临界值。

2.指数——用于计算达到临界塑性应变值之前单元的出砂量的出砂函数的指数。

3.硬化模型参数——三轴试验数据用修正Bigoni-Piccolroaz(MBP)硬化模型的一个版本进行匹配(拟合)。

为了将MBP硬化模型拟合到三轴试验数据，除了使用表2给出的参数之外，可以使MBP模型允许用户定义八面体面中的屈服面形状(常数I₁)(见图23A和23B)。八面体面中的屈服面形状由MBP硬化模型的两个参数定义的：偏差(γ)和贝塔(β)。偏差和贝塔二者都是值介于0与1.0之间的常数。在MBP模型中将贝塔设置成0及将偏差设置成0产生Drucker-Prager材料模型的八面体面中的屈服面(如图23A所示)。在MBP模型中将贝塔设置成0及将偏差设置成0.95产生Lade型材料模型的八面体面中的屈服面(如图23B所示)。

当样品达到屈服面时，发生屈服。因此，在图22A和22B所示形状的内部区域，线弹性行为占据支配地位。从图22A和22B看出，与Drucker-Prager型模型相比，应力路径命中Lade型模型的偏平面中的屈服面的可能性更高。预期的是，在将两个材料模型应用于一般应力状态时，随着偏差值的增大将会发生更多屈服。注意，在应用于三轴压缩试验时，图23A和23B中的屈服面产生相同的结果。

6.3.1.1网格单元的定义

计算在2D轴对称网格上执行。用于表示样品的网格单元大小被改变以研究模型结果的网格大小相关性。在垂直方向上，样品的单元大小是1.284毫米，在水平方向上，样品单元大小从0.05毫米变成0.01毫米乃至0.005毫米。测量精确至大约0.1cc，由此在模型中使用小的单元大小来获取这个精度水平。在垂直方向上，间隔物中的单元大小是2毫米。这导致分别具有80300、401100和802100个单元的模型。发现运行时间从大约2小时渐增至大约14个小时乃至大约30个小时。图24显示了作为单元大小函数的模型结果的发展。从水平单元大小为0.01毫米到0.05毫米的网格没有观察到显著的预测出砂量变化。减小运行时间到14小时以下，创建额外的网格单元模型，该模型具有与先前给出的垂直单元尺寸相同的垂直单元尺寸，但是只有前7毫米的样品用0.01毫米的单元大小网格化。超过初始的7毫米之外，单元大小增大到0.11035毫米。这个渐变网格单元模型具有100,000个单元，运行大约3-4个小时。用该网格单元模型的计算与来自其他网格单元模型的结果相比较，可以看出，渐变网格单元模型的结果与更精细的网格单元模型重叠，表明结果独立于单元大小(图24)。这种渐变网格构造用于所有后续运行。

6.3.1.2到出砂函数的拟合

用于拟合的出砂函数是仿真通过求解耦合的地质力学模型(包括出砂函数)和流体流动模型的偏微分方程组来进行。对于0.5的偏差、0的β值、2的指数(m)值以及0.017的临界应变限制获得最佳拟合(参见图25)。数据被拟合以用于曲线的大规模出砂部分和逐渐出砂部分。在图25中指示了塑性变形起始。在图26中显示了水速率(water rate)匹配。观察到出砂量和水速率的接近匹配。归因于指定给出砂单元的大渗透率值，在稍后时间观察到水速率的轻微失配。试验结果表明，在数值模型中，渗透率增大稍微高估了水速率。此外还观察到，0.017的临界应变值表明出砂发生在远超样品的初始屈服时(图27)，由此表明在将初始屈服点作为出砂起始时，可能略微高估样品的出砂趋势。

为了获取对出砂量的拟合，首先，偏差将被改变，其中指数(m)保持恒定在大小为1的值。对于每个偏差值，选择临界应变值以与大规模出砂时的屈服值相似。为了获取该值，使用包含硬化模型而不是出砂模型的地质力学模型来执行拟合。使用多个偏差值、临界应变值和指数值以研究解的唯一性。

首先，偏差值从0.92变到0.5。如前所述，通过使用包含硬化模型而不是出砂模型的地质力学模型来执行拟合，为每个偏差确定临界塑性应变值。大规模出砂起始处的屈服值用作临界塑性应变的初始值。有时进行对临界塑性应变值的调整以获得给定偏差值下的最佳可能匹配。使用的最大偏差值是0.92。这个值从Matsuoka-Nakai材料模型确定(MBP模型的一个版本)。如图28所示，0.5的偏差给出最接近的匹配，因此，在所有后续仿真中都使用了这个值。

指数(m)的变化表明，指数值只影响出砂量，但是它不改变发生大规模出砂的点。较大的指数值减小出砂总量(图29)。改变临界塑性应变的值表明，该常数支配大规模出砂的起始，较大的指数值导致较晚时间的出砂(即处于较大的约束应力下)(参见图30)。减小临界塑性应变值导致以更小的约束应力发生出砂。

图31显示进行空心圆柱试验后源自出砂的材料故障的形状演变。图31中的每个面板是在岩芯样品中不同深度处获取的水平切片。从图31观察到，出砂在样品的中心岩芯切片处规模最大(也就是说，腔在岩芯中心处最深)，并且朝着岩芯样品末端减小。在数值仿真运行中也观察到这种形状。图32显示了为数值仿真计算的材料故障形状随时间的仿真演变(图31中的每个切片被垂直取得)。虽然初始破坏处于岩石-间隔物界面处，但是腔发展得使最深腔发生在岩芯的中心(见图32)，也就是说，腔朝向仿真岩芯的中心变深，这与空心圆柱试验结果中观察到的相匹配。

6.3.2腔生长

从图32可以看出，在116675秒处，砂仍被产生。为了确定出砂是否稳定化，在116675秒处的负载(即46MPa的约束应力和0.1643398MPa的压差)将对额外的116675秒保持恒定。然而，仿真设置成如果超出50％的样品作为砂产生则终止。图33A显示出砂(cc)与时间(s)的关系曲线，图33B显示数值仿真的腔大小，其中从第116675秒开始约束应力恒定在46MPa，压差恒定在0.1643398MPa。使用连续的0.01毫米大小的网格。在图B中，灰色区域表示腔的位置，暗区域表示完好岩石(或间隔物)。在第126618.836995秒，50％的样品作为砂产生，仿真终止(见图33A和33B)。在体积计算中，还计算了间隔物的体积。如图33B所示，腔在中心处很深，几乎贯穿整个仿真岩石样品。由此断定，46MPa的约束应力和0.1643398MPa的压差导致不稳定的出砂。

为了确定是否可以发展和维持稳定的腔，负载被改变。所使用的材料常数与图25中相同，即偏差＝0.5、贝塔＝0、临界塑性应变值＝0.017以及指数(m)＝2。结合图32概述的类似加载过程被使用。特别地，在第113230秒负载达到44MPa的约束应力和1.255317MPa的压差，并且负载对额外的120120秒保持恒定(最大仿真时间233350秒)。该仿真显示，产生了最大10.2966cc的砂，腔是稳定的(图34)。该观察结果表明存在这样的条件，在该条件下，可以预期一定程度的出砂，但是该条件不会导致灾祸性故障(由于大规模出砂引起)。

6.4总结

结果表明，使用包含出砂函数的地质力学模型，可以唯一地拟合空心圆柱出砂试验的出砂起始和大规模出砂二者。用于拟合数据的参数是偏差、贝塔、临界应变和指数。对于Landana盐水试验，发现以下参数值导致最佳匹配：

偏差＝0.5

贝塔＝0

临界应变值＝0.017

指数＝2

结果还表明存在有条件，在该条件下可以预期一定数量的出砂，但是该条件不会导致灾难性故障(因为大规模出砂而引起)。用于第6.3.1节的拟合过程的工作流程如下：

1)基于总砂体积以及测量砂体积的精度，选择恰当的网格大小。对于该研究来说，如果产生总量约6cc的砂，那么用于获取预期精度的单元大小是0.01毫米。

2)使用从到三轴试验数据的拟合获取的材料常数，创建具有变化的偏差值的多个模型。在拟合试验数据的计算中，地质力学模型不包括出砂模型。

3)使用步骤2的结果，记录大规模出砂时的屈服值。该屈服值设置成与临界塑性应变限制相等，指数值保持等于1。在该步骤中，可以用所获取的值附近范围的临界塑性应变值来进行多次运行。

4)使用步骤3的结果，用于拟合大规模出砂起始的仿真被用于改变指数值，其中观察到指数值的增大减少了产生的砂量。与来自试验的数据最紧密拟合的仿真被选定成为代表性的。

7.引用的参考文献

这里提及的所有参考文献都通过引用整体合并于此且用于所有目的，程度上如图每个单独公开或专利或专利申请被具体或单独地描述为通过引用整体合并于此用于所有目的。这里对参考文献的论述或提及不应被解释成是对所述参考文献是本发明的现有技术的认可。

8.修改

如对本领域技术人员来说显然的那样，在不脱离本发明的实质和范围的情况下，可以对本发明进行许多修改和变化。这里描述的具体示例仅作为示例提供，本发明只受附加权利要求项以及这些权利要求享有的全部范围的等价物的限制。

作为这里描述的系统和方法的宽范围示例，这里描述的系统和方法可以通过包含了可供设备处理子系统运行的程序指令的程序代码而在众多不同类型的处理设备上实施。软件程序指令可以包括源编码、目标代码、机器码或是可以操作来促使处理系统执行这里描述的方法和操作的任何其他存储数据。然而，其他实施方式也可被使用，例如配置来执行这里描述的方法和系统的固件或甚至恰当设计的硬件。

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

所述系统和方法可以提供在众多不同类型的计算机可读介质(例如CD-ROM、磁盘、RAM、闪存、计算机硬盘驱动器等)上，该介质包括可供处理器运行以便执行方法操作以及实施这里描述的系统的指令(例如软件)。

这里描述的计算机组件、软件模块、函数、数据存储和数据结构可以直接或间接彼此连接，以便允许它们的操作所需的数据流动。还应指出，模块或处理器包括但不局限于执行软件操作的代码单元，可实施为例如代码的子程序单元、代码的软件功能单元、作为对象(如在面向对象范例中一样)、作为小应用程序、在计算机脚本语言中或是作为别的类型的计算机代码。根据当前状况，软件组件和/或功能可以位于单个计算机上，或分布在多个计算机上。

Claims (10)

1.一种在预测地质力学油藏系统的出砂中使用的方法，包括：

接收指示地质力学油藏系统内的材料的塑性变形和出砂数据；

在计算机系统上，基于硬化模型对所接收的塑性变形数据的第一拟合来计算一个或多个硬化参数的值；其中所述硬化模型对地质力学油藏系统内的材料的塑性变形行为进行建模；

在计算机系统上，基于包含所述硬化模型的地质力学模型对所接收的出砂数据的第二拟合来计算临界塑性应变的值；

在计算机系统上，基于包含出砂函数的地质力学模型对所接收的出砂数据的第三拟合并使用所述临界塑性应变的所述值和所述硬化参数的至少一个值的值来计算所述出砂函数的指数值；

其中所述出砂函数预测所述地质力学油藏系统的出砂，所述出砂函数由下面的表达式给出：

$f\left(x\right)={\left(\frac{{\mathrm{\ϵ}}^p}{{\mathrm{\ϵ}}_{\lim }^p}\right)}^m$

其中其中ε^p是塑性应变不变量，其中是所述临界塑性应变；以及其中m是所述指数值；以及

将所述一个或多个硬化参数的至少一个值、所述临界塑性应变的所述值、或所述出砂函数的所述指数值输出到用户接口设备、监视器、计算机可读存储介质、本地计算机或作为网络的一部分的计算机。

2.根据权利要求1所述的方法，其中所述出砂函数是函数f(x)，其中当x＝0时，f(x)＝0，并且当x＝1时，f(x)＝1；且其中x是临界塑性应变的函数。

3.根据权利要求1所述的方法，其中所述硬化模型是修正的Bigoni-Piccolroaz模型。

4.根据权利要求1所述的方法，其中所述塑性变形数据从三轴试验获得。

5.根据权利要求1所述的方法，其中所述出砂数据从空心圆柱试验获得。

6.根据权利要求1所述的方法，其中所述硬化模型对所接收的塑性变形数据的第一拟合通过回归获得。

7.根据权利要求1所述的方法，其中：

包括所述硬化模型的所述地质力学模型对所接收的出砂数据的所述第二拟合通过求解对地质力学油藏系统建模的偏微分方程组获得；

所述偏微分方程组包括油藏流动模型以及包含了所述硬化模型的所述地质力学模型；其中所述偏微分方程组通过完全展开的雅可比行列式耦合；以及

所述偏微分方程组的求解包括基于所接收的出砂数据在单个时间步长中同时求解完全展开的雅可比行列式。

8.根据权利要求1所述的方法，其中：

包含所述出砂函数的所述地质力学模型对所接收的出砂数据的所述第三拟合通过求解对地质力学油藏系统建模的偏微分方程组获得；

所述偏微分方程组包括油藏流动模型以及包含了所述硬化模型的所述地质力学模型；其中所述偏微分方程组通过完全展开的雅可比行列式耦合；以及

所述偏微分方程组的求解包括基于所接收的出砂数据在单个时间步长中同时求解完全展开的雅可比行列式。

9.根据权利要求1所述的方法，其中所输出的所述一个或多个硬化参数的至少一个值、所述临界塑性应变的所述值、或所述出砂函数的所述指数值视觉显示在用户接口设备或监视器上。

10.一种用于地质力学油藏系统的出砂的方法，包括：根据权利要求1的方法的结果来操作所述地质力学油藏系统。