一种考虑颗粒大小分布和孔隙度的致密数字岩心生成方法

技术领域

本发明属于工程热物理、化学和石油工程等领域，尤其是涉及一种以多孔介质孔隙结构为基础的实现微观多相多组分流动、流-固耦合等孔隙尺度数值模拟的数字岩心生成的方法。

背景技术

在工程热物理、化学和石油工程等领域中，随机的球形颗粒填砂岩芯或者由圆形颗粒组成的微观模型作为一个研究工具，在微观多相多组分流动、流-固耦合等研究中具有非常重要的作用。其中，填充岩芯的孔隙度，填充配位数和填充颗粒尺寸分布是颗粒随机填充的三个主要表征参数。孔隙度是填砂岩芯中被孔隙空间占据总体积的比例，配位数是某一颗粒周围相邻的颗粒数，颗粒尺寸分布是填砂岩芯中各种存在的颗粒大小(通常是颗粒直径)的相对数量。许多颗粒状的物质，例如几种类型的土壤和真实的粉末，可以用具有对数正态尺寸分布的随机填充球形颗粒来表示。

随机填充岩心可以通过实验、理论和计算机模拟进行建模。使用计算机模拟方法的优势在于，已知具有不同尺寸大小的球体的配位数，可以使用配位数来计算所生成多孔介质的其他几何特性。已经开发了许多计算机模拟算法来生成球形颗粒的虚拟随机填砂岩芯。传统的计算机球形颗粒填充方法包括顺序增加(或者球体沉降)算法、随机放置算法、生长算法和动态仿真模型等。在顺序增加方法中每个新增加颗粒的向下运动导致垂直方向的各向异性较弱，因为垂直方向的作用不同于水平方向，所以通常情况下，如果没有额外的压制或致密化程序，这种填充不是很致密。在随机放置算法中，如果新引入的颗粒与之前放置的颗粒发生重叠，通过丢弃该新引入的颗粒来消除重叠，由于可能必须丢弃一些颗粒，所以填充球的尺寸分布将与预期的初始分布有所不同；此外，该方法填充致密程度取决于颗粒添加的顺序，例如先添加较大的颗粒球体，与先放置较小的颗粒球体相比，最终填充密度往往会更高。在生长算法中，如果没有进一步的颗粒位置重新排列，则该算法通常不会导致致密的填充。如果采用颗粒位置重新排列，例如通过允许粒子绕其初始位置移动，和/或通过使其收缩或扩展，可以生成更致密的填充岩芯，但也无法精确控制最终填充中的尺寸分布，同时如果进行重排，该方法也可能非常耗时。另外还有动态仿真模型(例如离散元方法)，除了可能面临与其他传统填充算法相同的数学和实现困难之外，与专用填充算法相比，它们往往对计算的要求更高，因为它们还必须正确模拟填充过程的动力学，而不仅仅是最终结果，所以需要的计算时间长，同时随着颗粒复杂度的增加，这也会导致仿真速度下降。

原理上，以上方法不能同时满足：具有大颗粒优先填充，能移动、收缩和扩展颗粒尺寸的特点，使其能生成更致密的多孔介质；在满足致密的前提下能控制填充颗粒的尺寸大小分布；较低的计算时间。在实现方法上，前人采用的溶胀球算法只考虑了考虑颗粒大小分布作为输入参数，未考虑孔隙度的限制，但对于较致密的多孔介质，限制填充岩芯的孔隙度是不可忽略的控制致密的因素。

发明内容

本发明的目的是为了同时满足具有大颗粒优先填充、能移动、收缩和扩展颗粒尺寸大小的特点，能生成更致密的多孔介质；在满足致密的前提下能控制填充颗粒的尺寸分布；在计算过程有较低的计算时间。本发明是一种同时考虑孔隙度和颗粒大小分布这两种孔隙结构参数的致密数字岩心生成方法。这对于研究孔隙尺度下致密多孔介质内多相多组分流动更加方便，能加深孔隙结构和流体流动耦合的理解。具有生成数字岩心的孔隙结构更简洁、高效，孔隙结构建模更便利等优点。

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

一种考虑颗粒大小分布和孔隙度的生成致密数字岩心的方法，根据实际岩石的孔隙结构物性参数，通过随机生长膨胀收缩策略得到球体颗粒堆，针对球体的颗粒实现球体颗粒生长和移动过程的模拟，模拟过程中颗粒的半径不断膨胀，如果颗粒之间出现覆盖的情况，则移动颗粒消除覆盖现象。当通过移动颗粒位置无法消除颗粒间的覆盖时，收缩颗粒半径。

一种考虑颗粒大小分布和孔隙度的致密数字岩心生成方法，包括步骤如下：

(1)输入参数确定：根据目标岩芯或者多孔介质的孔隙结构物性参数，指定孔隙度、颗粒尺寸分布、最大半径为和最小半径，确定颗粒数量、数字岩心的区域尺寸大小和颗粒膨胀次数；

(2)初始设置：根据颗粒数量，在数字岩心的区域内随机分布地设置颗粒点(图2a)；每个点的位置都是构成最终填充球体的初始中心，该点的半径设置为零作为最终填充球体的初始半径；根据步骤(1)的颗粒尺寸分布将最终颗粒半径分配给每个点；根据每个点最终颗粒半径和膨胀次数，将半径增长率的比值分配给这些点；

(3)颗粒位置排序：将步骤(2)生成的随机分布的颗粒，按照球体颗粒中心坐标的位置大小重新排列，生成一个有序链表，然后给每一个颗粒中心代表的球体颗粒指定一个ID；

(4)颗粒半径膨胀：球体开始生长膨胀，球体的半径按初始指定的半径增长率变大(图2b)；在颗粒膨胀步骤中，球体之间可能发生重叠(图2c)；但是，两个刚性颗粒不能占有相同的孔隙空间；所以在每个球体经过膨胀步骤后，必须检查填充盒中的所有球体颗粒之间是否发生重叠，球体中心通过轻微重新排列以消除重叠。

(5)颗粒覆盖判断和移动：当某一个颗粒ID＝n(0<＝n<＝N-1)，与其之后的颗粒ID＝n+i(1<＝in的颗粒不再重叠时开始检查ID＝n+1与ID>n+1间的重叠情况。去除一个重叠可能会导致新的重叠，也有可能还存在其他重叠。因此去除重叠的过程将持续进行，直到在下一个球体生长步骤之前在数字岩心的区域内没有发现重叠情况为止，即从ID＝0与ID＝1开始检查当前结果是否存在球体颗粒间覆盖的情况，一直持续到ID＝N-2和ID＝N-1之间判断完结束。

(6)颗粒半径收缩：随着球体颗粒半径的增加，球体间重叠无法通过移动球体中心解决(图2e)。在达到移动次数Nmax步之后，无法移动的颗粒和ID＝n的颗粒两个球体的半径相对于其各自半径收缩率收缩(图2e)，使其可以适合填充空间，使两个颗粒不再重叠(图2f)，然后对该ID后面重叠的颗粒再次执行步骤(5)。

(7)目标颗粒半径判断：颗粒位置排序、颗粒半径膨胀、颗粒覆盖判断和移动、颗粒半径收缩步骤一直持续到球体半径达到其最终分配值(图2g)。

(8)根据步骤(1)中区域尺寸生成矩形或者长方体，用矩形或者长方体对步骤(7)生成的颗粒做差集，得到目标随机多孔介质(图2h)。

进一步地，所述颗粒尺寸分布包括对数正态分布、均一分布、韦布尔分布、单一分布中的一种。

在步骤(1)中，根据颗粒的颗粒大小分布、填充区域的几何结构，指定颗粒生长区域以及颗粒的尺度分布特征，在数字岩心的区域内进行颗粒生长。

生成的多孔介质图案在各种空间尺度上都可以使用，包括但不限于纳米、微米、厘米和米等尺度。

发明的有益效果是：本发明根据目标油藏或者其他地质结构、填砂岩心等储层岩心的孔隙度和颗粒尺寸分布，通过一种模拟球型颗粒膨胀和收缩的过程的计算机模拟方法，较为精确的生成与地下储层中的孔隙结构异质性和复杂性相符的虚拟致密数字岩心。该生成方法结合孔隙尺度模拟方法，为生成具有表征致密岩心孔隙结构提供了更简洁的方法，为模拟多孔介质内多相多组分流动运移孔隙结构建模提供便利，为加深孔隙结构和流体流动耦合的理解提供基础，推动揭示微观孔隙流动规律研究的发展。

附图说明

图1为考虑颗粒大小分布和孔隙度的致密数字岩心生成方法的算法流程图。

图2实施例不同步骤示意图。

图3颗粒半径分别为1μm，1.25μm，1.5μm，1.75μm，2μm的正态分布随机多孔介质。

图4颗粒半径分别为60μm，65μm，75μm，94μm，100μm的正态分布随机多孔介质。

具体实施方式

以下结合技术方案和附图详细叙述本发明的具体实施例。

实施例1

(1)输入参数确定：根据目标岩芯或者多孔介质的孔隙结构物性参数，指定孔隙度为0.15、颗粒尺寸分布为正态分布，最大半径为2μm，最小半径为1μm，确定颗粒数量为866，和数字岩心的矩形区域尺寸大小为宽120μm高60μm，膨胀次数为30；

(2)初始设置：根据颗粒数量，在数字岩心的区域内随机分布地设置866个点的位置；每个点的位置都是构成最终填充球体的初始中心，该点的半径设置为零作为最终填充球体的初始半径；根据步骤(1)的颗粒尺寸分布将最终颗粒半径分配给每个点，颗粒半径分别为1μm，1.25μm，1.5μm，1.75μm，2μm；根据每个点最终颗粒半径和膨胀次数，将半径增长率的比值分配给这些点；

(3)颗粒位置排序：将步骤(2)生成的随机分布的颗粒，按照球体颗粒中心坐标的位置大小重新排列，生成一个有序链表，然后给每一个颗粒中心代表的球体颗粒指定一个ID(范围为0～865)；

(4)颗粒半径膨胀：球体开始生长膨胀，球体的半径按初始指定的半径增长率变大；在颗粒膨胀步骤中，球体之间可能发生重叠；但是，两个刚性颗粒不能占有相同的孔隙空间；所以在每个球体经过膨胀步骤后，必须检查填充盒中的所有球体颗粒之间是否发生重叠，球体中心通过轻微重新排列以消除重叠。

(5)颗粒覆盖判断和移动：当某一个颗粒ID＝n(0<＝n<＝866)，与其之后的颗粒ID＝n+i(1<＝i<866)，两个颗粒之间的球体中心间距离小于两个球体颗粒的半径之和时，判断两个颗粒存在重叠的情况，如果不是则表明两个颗粒未发生重叠的情况；当确定ID＝n与ID＝n+i发生重叠时，针对ID＝n+i颗粒，向ID＝n+i的中心相对ID＝n的中心相反的方向移动两个球颗粒中心距离的0.12倍，不断移动直到这两个颗粒不再重叠，之后对ID＝n之后的颗粒重新按照空间位置排列指定ID，在ID＝n上一直检查到与ID>n的颗粒不再重叠时开始检查ID＝n+1与ID>n+1间的重叠情况。去除一个重叠可能会导致新的重叠，也有可能还存在其他重叠。因此去除重叠的过程将持续进行，直到在下一个球体膨胀步骤之前在特定区域内没有发现重叠情况为止，即从ID＝0与ID＝1开始检查当前结果是否存在球体颗粒间覆盖的情况，一直持续到ID＝864和ID＝865之间判断完结束。

(6)颗粒半径收缩：随着球体颗粒半径的增加，球体间重叠无法通过移动球体中心解决。在达到移动次数44步之后，无法移动的颗粒和ID＝n的颗粒两个球体的半径相对于其半径收缩0.10，使其适合填充空间，使两个颗粒不再重叠，然后对该ID后面重叠的颗粒再次执行步骤(5)。

(7)目标颗粒半径判断：颗粒位置排序、颗粒半径膨胀、颗粒覆盖判断和移动、颗粒半径收缩步骤一直持续到球体半径达到其最终分配值，运行时间75.638s。

(8)根据步骤(1)中区域尺寸生成尺寸大小为宽120μm高60μm的矩形，用矩形对步骤(7)生成的颗粒做差集，最终的孔隙度为0.32，得到目标随机多孔介质如图3。

实施例2

(1)输入参数确定：根据目标岩芯或者多孔介质的孔隙结构物性参数，指定孔隙度为0.15、颗粒尺寸分布为正态分布，最大半径为100μm，最小半径为60μm，确定颗粒数量为761，和数字岩心的矩形区域尺寸大小为宽6000μm高3000μm，膨胀次数为30；

(2)初始设置：根据颗粒数量，在数字岩心的区域内随机分布地设置761个点的位置；每个点的位置都是构成最终填充球体的初始中心，该点的半径设置为零作为最终填充球体的初始半径；根据步骤(1)的颗粒尺寸分布将最终颗粒半径分配给每个点，颗粒半径分别为60μm，65μm，75μm，94μm，100μm；根据每个点最终颗粒半径和膨胀次数，将半径增长率的比值分配给这些点；

(3)颗粒位置排序：将步骤(2)生成的随机分布的颗粒，按照球体颗粒中心坐标的位置大小重新排列，生成一个有序链表，然后给每一个颗粒中心代表的球体颗粒指定一个ID(范围为0～760)；

(4)颗粒半径膨胀：球体开始生长膨胀，球体的半径按初始指定的半径增长率变大；在颗粒膨胀步骤中，球体之间可能发生重叠；但是，两个刚性颗粒不能占有相同的孔隙空间；所以在每个球体经过膨胀步骤后，必须检查填充盒中的所有球体颗粒之间是否发生重叠，球体中心通过轻微重新排列以消除重叠。

(5)颗粒覆盖判断和移动：当某一个颗粒ID＝n(0<＝n<＝760)，与其之后的颗粒ID＝n+i(1<＝i<760)，两个颗粒之间的球体中心间距离小于两个球体颗粒的半径之和时，判断两个颗粒存在重叠的情况，如果不是则表明两个颗粒未发生重叠的情况；当确定ID＝n与ID＝n+i发生重叠时，针对ID＝n+i颗粒，向ID＝n+i的中心相对ID＝n的中心相反的方向移动两个球颗粒中心距离的0.15倍，不断移动直到这两个颗粒不再重叠，之后对ID＝n之后的颗粒重新按照空间位置排列指定ID，在ID＝n上一直检查到与ID>n的颗粒不再重叠时开始检查ID＝n+1与ID>n+1间的重叠情况。去除一个重叠可能会导致新的重叠，也有可能还存在其他重叠。因此去除重叠的过程将持续进行，直到在下一个球体膨胀步骤之前在特定区域内没有发现重叠情况为止，即从ID＝0与ID＝1开始检查当前结果是否存在球体颗粒间覆盖的情况，一直持续到ID＝759和ID＝760之间判断完结束。

(6)颗粒半径收缩：随着球体颗粒半径的增加，球体间重叠无法通过移动球体中心解决。在达到移动次数1344步之后，无法移动的颗粒和ID＝n的颗粒两个球体的半径相对于其半径收缩0.10，使其适合填充空间，使两个颗粒不再重叠，然后对该ID后面重叠的颗粒再次执行步骤(5)。

(7)目标颗粒半径判断：颗粒位置排序、颗粒半径膨胀、颗粒覆盖判断和移动、颗粒半径收缩步骤一直持续到球体半径达到其最终分配值，运行时间146.374s。

(8)根据步骤(1)中区域尺寸生成尺寸大小为宽6000μm高3000μm的矩形，用矩形对步骤(7)生成的颗粒做差集，最终的孔隙度为0.34，得到目标随机多孔介质如图4。