试系统102通过该过程执行颗粒分配和颗粒接触识别，以识别如由图像体积128表示的岩石样品104的每个颗粒之间的单独颗粒以及接触区域。接触区域相应于与其它颗粒接触的个体颗粒的那些表面部分。在本发明的一些实施例中，颗粒之间的接触区域以及它们的诸如接合度、粗糙度等的特性的分析是有用的，因为这些接触特性能够影响岩石的应力-应变响应。图3C和3D示出在可选过程212的情况下，在图3B的分割2D切片图像300上执行的颗粒分割和颗粒接触识别的示例。如图3C中所示，2D切片图像中的每个不同颗粒都被随机地遮蔽为(shaded)不同的灰度值，以将颗粒彼此区分。每个个体颗粒被遮蔽的特殊灰度值反映了用于在固体基质中识别个体颗粒的独特数值标签。如图3D中所示，作为可选过程212的结果，以与它们相应颗粒的主体不同的灰度值强调每个不同颗粒的颗粒间接触。
[0054]因而，过程210 (包括如果执行的可选过程212)将所分割的数字图像体积内的体素与岩石样品104中的相应位置处的特殊材料(或者在可能情况下的孔隙空间)相关联。在过程210 (和如果执行的可选过程212)中，一些或者全部体素每个都标记以一种或者更多材料属性，该一种或者更多材料属性对应于通过过程210、212被分配给该体素的特殊材料成分(constituent)的，这些成分包括孔隙空间、基质材料、粘土比率、个别颗粒、颗粒接触点、矿物质类型等。那些识别的成分的特殊弹性或者其它材料属性，以对将执行的分析有用的程度，与相应的体素相关联，即体积内的颗粒和矿物质被分配适当的密度和弹性属性。
[0055]例如，当假定个体颗粒、矿物质和接触根据线性弹性表现时，对被标记为个体颗粒、矿物质或者接触的每个体素分配杨氏模量E和泊松比的值是有用的。在本领域中已知，杨氏模量为材料经历线性(即，作为应变的函数的应力关系为线性，斜率等于杨氏模量E的值)的单轴应力变形时的刚度的量度。在本领域中也已知，泊松比为单轴应力行为的条件下的侧向和横向应变的量度。可替选地，可以向材料中的颗粒、矿物质和接触分配体积模量K和剪切模量G的值，以描述那些组分的弹性行为。在本领域中已知，体积模量为材料对静水压力的弹性响应的量度，而剪切模量为材料对剪切应变的弹性响应的量度。在本领域中已知，所有这些弹性系数都通过众所周知的变换而彼此相关。预期的是，对于其中考虑线性弹性材料的那些情况，通常将杨氏模量和泊松比归属于材料的组分(component)，因为能够通过实验直接确定这些参数的值。
[0056]在假定矿物质、颗粒或者接触点展现黏弹性行为使得响应于所施加的应力或者应变的变形是率相关的情况下，必需指定适当的模型参数，像刚度和粘性，例如如果假定是麦克斯韦体。在本领域中已知适合黏弹性和塑性材料并且可以用于描述各种类型的应力/应变行为的多种其它的本构模型。在任何情况下，指定给材料的模型参数应为适合所假定的特定本构模型的那些参数。
[0057]然后由测试系统102执行过程214，以在岩石样品104的所分割的3D图像体积内产生用于固体材料(或者对于来自过程212的所分隔的识别的颗粒和接触区域)的有限元网格。在本发明的实施例中，计算装置120执行测试工具130，以产生这种有限元网格，作为施加给所分割的3D图像体积的非结构网格。这种有限元网格是“非结构的”，因为其由不规则图案(即，具有不规则连接)中的多个多边形单元组成，与处于规则图案(即，具有规则连接)中的多边形单元的“结构化”网格相反。在可选过程212中识别颗粒接触的本发明的实施例中，能够在所识别的接触区域中或者附近改进(即，更精细地图案化)非结构网格。然后，计算装置120也在过程214中将每个体素的每个标记组分的材料属性分配给非结构网格的相应单元。
[0058]图3E示出在过程214中从在过程210、212中产生的3D分割图像体积创建的非结构网格的示例。图3E中所示的视图是3D非结构网格的2D表示，其中表示固体材料302的图像切片的部分由彼此尺寸和连接性不同的有限元来表示。这些有限元中的每个也都被分配有相应于其覆盖的标记组分(例如，通常为固体材料302，或者所表示的特殊材料)的材料属性。虽然图3E示出在该视图中单个2D切片图像300和每个有限元的横截面(示出为三角形)，但是非结构网格的有限元被视为已经应用于由一系列这种2D切片图像组成的3D数字图像体积128的三维(四面体)单元。虽然图3E示出使用四面体单元产生网格，但是预期可能使用任何单元类型或者不同单元类型的组合以产生固体材料302的非结构网格。
[0059]在过程216中，测试系统102向3D图像体积128的非结构网格施加与应力、应变、力、位移等中的一种或者多种对应的模拟变形。在本发明的一些实施例中，测试工具130被配置成执行包括有限元(FE)求解器的一个或者更多软件程序912，以模拟岩石样品104在其处于地层中的地下位置的原位处遇到的变形条件。在本领域中已知，FE分析用于通过将解域分为更小的子域或者有限元而求解复杂问题。在非结构网格的背景下，如上所述，在相同解域中采用多种单元形状和尺寸。每个单元与多个节点相关联，相邻单元在节点处通常利用表示单元上的场变量的变化的内插函数(通常称为形状函数)彼此连接。通常基于在节点处建立平衡以及兼容性的物理自变量，用公式表示用于整个系统的联立代数方程组。通过分配相关变量或者节点负荷/力的特定节点值而在解域的边缘上施加边界条件。然后求解该方程组，获得未知节点值，诸如应力、应变、力和位移。在这种情况下，测试工具130被配置成包括FE求解器，其被实现为能够在过程216中，在过程214中定义的非结构网格上执行这种FE分析的必要逻辑、算法等。具体的FE求解器能够为任何类型的传统已知FE求解器，诸如线性直接求解器、迭代求解器、本征(eigen)求解器、非线性方程求解器或者另一 FE求解器。
[0060]在测试工具130采用有限元技术模拟施加给由数字图像体积128表示的岩石的体积的变形的本发明的实施例中，由计算装置120执行过程216，使具有标记的材料属性的有限元的非结构网格经过FE分析，以根据边界条件求解弹性、粘弹性或者其它适当的本构控制方程组，所述边界条件被以表示待被模拟的所期望的原位地下地层条件的方式分配到网格化体积的表面。例如，这些边界条件可能采取所施加的位移的形式，在该情况下，FE求解器计算网格体积的每个有限元的应力和应变。在其它实施方式中，向非结构网格施加牵引(即，应力)，在该情况下，FE求解器计算网格体积的每个有限元的应力和应变。优选地，所施加的变形的幅值和方向与所期望原位应力-应变条件对应，其示例包括静水力学、单轴和三轴应力-应变。在任一情况下，测试工具130都通过下列方式执行过程216，即对于所施加的边界条件，跨(across)由非结构网格表示的体积来数值求解适当的控制方程(即，诸如那些用于线性弹性的控制方程)。通过这些用于线性弹性的应力-应变计算，FE求解器也能够计算整个图像体积128的有效弹性属性(杨氏模量、泊松比、体积模量以及剪切模量等)。通常通过对涉及应力应变的刚度矩阵或者对涉及应力应变的柔度矩阵求解而找到这些弹性参数。模拟变形的效果影响非结构网格的结构和属性。图3F示出模拟变形的示例，其中已经响应于所施加的位移边界条件而在图3E中所示的网格上计算材料应力。通过比较图3E和3F应明白，在该示例中，这种模拟变形影响图像体积300在X方向中的压缩。
[0061]在图3F中，假定固体基质的弹性属性(E、)在体积各处都均匀，并且在变形阶段期间保持恒定。当存在粘土或者其它显著不同的弹性材料时，以被分配给每种矿物质(石英、粘土等)的弹性属性执行模拟是有用的。此外，当认为颗粒接触对诸如具有弱固结砂的岩石的总体力学表现具有显著影响时，考虑柔度/刚度是有用的，颗粒接触的存在导致柔度/刚度出现。通常，使用多种方法，诸如解析模型、实验数据或者启发函数，向随着所施加的应力/应变而变化的弹性属性分配应力相关接触柔度(法向的和切向的)。用于接触行为的解析模型(Hertz、Mindlin、Walton、Digby等等)通常假定球形颗粒接触，并且接触区域为圆形的。能够在模拟中应用这些模型，以在考虑每个颗粒的配位数的情况下调节每个个体颗粒的接触区域的弹性属性，所述颗粒的配位数表示该颗粒的颗粒接触的数量。此外，由于这些模型通常是所施加的应力的函数，所以能够随着变形推进，取决于在接触区域中计算的增量应力或应变，来调节接触弹性属性。如上所述，另一种向接触区域分配弹性属性的方法为采用实验数据，其中使用作为应力的函数测量的动态弹性属性(压缩和剪切波速)，例如通过假定样品128的静态弹性属性(杨氏模量、泊松比)等于从所测量的波速提取的动态弹性属性(杨氏模量、泊松比)，来校准接触柔度。
[0062]为了在所模拟的变形中考虑到接触刚度/柔度影响，必需执行可选过程212，在可选过程212中，向分割的体积施加颗粒和接触分隔。图4A和4B示出这种颗粒和接触分隔的示例。在图4A中，示出在分隔过程212之前的固体基质材料。图4B示出过程212之后的相同材料，以黑色值示出颗粒分隔，并且以浅灰色值识别颗粒之间的接触。图4C示出在过程214中产生的用于图4B的颗粒分隔的体积的网格。在图4C的示例中，相对于固体颗粒内部的较大三角形尺寸，在图4C中以较小的三角形示出接触区域附近的网格的细化。通过执行一系列的变形，以较小的增量数值实现所期望的应力/应变条件，从而达到所期望的原位应力/应变条件。在每个增量变形之后，通常使用过程212，在已变形体积的体素化表示上创建体积的新颗粒和接触分隔。以这种方式在接触内的网格细化通常有用，因为弹性属性在接触区域和颗粒区域之间差异显著。在图2中以虚线示出其中重复过程212、214、216的这种增量网格细化方法。可替选地，接触区域内的行为能够通过对于固体颗粒的内部和接触区域两者使用适当小的网格单元来表征，代价为由于模型内的单元数目更大而需要更多计算。
[0063]如上所述，能够使用从实验近似或者假定行为遵守启发函数的解析模型对接触区域的弹性属性建模。在图4G中，为根据位移(在该示例中表达为百分比体积变化)的函数而变化的接触区域的弹性属性，显示两个不同函数。使用菱形符号的图示假定接触区域的杨氏模量小于固体颗粒的杨氏模量，并且随着变形增大是恒定的。使用交叉符号的图示假定接触区域的杨氏模量随着变形增大而非线性变化。如果怀疑颗粒包含低于图像分辨率的柔性孔隙率，也可能随着变形而改变个体颗粒的弹性属性。在图4D至4F中，正应变被示出在体积中:在变形之前(图4D);无颗粒接触行为情况下在变形中的一个增量步骤之后(图4E);以及假定颗粒接触行为根据非线性启发函数而变化的变形中的一个增量步骤之后(图4F)。这些图4D至4F示出与不包括接触行为的变形相比，包括接触行为的变形导致的颗粒形状和孔隙空间的显著差异。特别地，当考虑到接触行为时，如降低孔隙率以及改变颗粒形状表现的，出现更大的体积变形。图4F也示出与颗粒区域中的那些相比，在接触区域中产生的非常不同的应变。特别地，图4F示出，在变形中的一个增量之后，在增量变形之前不接触的一些颗粒现在接触，要求在执行后续的变形之前重复颗粒分隔过程212。
[0064]在图4H中，为三种不同的变形场景绘制孔隙率。第一种假定是在于体积各处弹性属性是均匀的，且无接触行为(图4E)。第二种假定是用于接触行为的启发函数随着变形是恒定的，并且第三种假定是，用于接触行为的启发函数随着变形是非线性的(图4F)，在图4G中示出两者。如这些图中所示，从当使用启发函数并考虑到接触行为时出现的孔隙率的显著降低，额外