CN116882255B - 一种基于傅里叶级数随机生成多孔介质模型的方法及系统 - Google Patents

Abstract

本发明公开了一种基于傅里叶级数随机生成多孔介质模型的方法及系统，方法包括以下步骤：对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置；初始化多孔介质模型，并生成颗粒的位置信息；提取颗粒轮廓边缘；获得填充颗粒；使填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；预设循环弹出条件，若判断结果满足循环弹出条件，则继续；否则，更新傅里叶参数；将满足循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；判断生成的模型是否满足预设生成要求，若满足，则输出多孔介质模型。该算法的加入大大提高了该方法的计算速度，最多可以在两个循环中完成单个粒子的重建。

Description

一种基于傅里叶级数随机生成多孔介质模型的方法及系统

技术领域

本发明属于多孔介质模型生成技术领域，具体涉及一种基于傅里叶级数随机生成多孔介质模型的方法及系统。

背景技术

石油、煤层气等化石能源所处于的地质储层本质是上是由颗粒堆积的复杂多孔介质结构。随着成像和数值技术的发展，使得孔隙尺度计算成为了研究储层内流体的时空特性的主要技术。获取孔隙尺度计算中的计算模型，目前主流的技术主要是CT扫描和数字生成。相比较而言，CT扫描的成本高、效率低，而且伴随着很大样本局限性，无法为普适性规律的研究提供大量的计算模型。而对于目前的数字生成的多孔介质的途径，主要随机场模型和颗粒堆积模型。对于随机场模型，可以实现满足Minkowski函数的多孔介质生成，但是在低孔隙度的模型上，很难保证模型的连通性，这对于孔隙尺度模拟时没有意义的。同时随机场模型是很难真实还原地质储层内颗粒的复杂构型的。对于颗粒堆积模型，从原始的球形颗粒演化为非规则性颗粒，颗粒堆积模型逐渐向真实孔隙介质结构逼近。针对非规则性颗粒的生成，主要的方法有复合颗粒和单颗粒的方法。在复合粒子方法中，使用圆盘/球体的团簇来近似不规则形状的粒子。尽管团簇方法在建模各种形状的粒子时提供了灵活性，但该方法固有的不切实际的多个子接触会产生生成误差。在单粒子方法中，可以用闭合曲线或曲面复制粒子形状，包括椭圆/椭球、超二次曲面、多边形/球体/多面体、非均匀有理基样条、水平集函数和傅里叶级数。上述方法中水平集函数和傅立叶级数在理论上可以再现任意粒子形状。然而，水平集函数在计算上是昂贵的，并且需要相当大的存储空间。与水平集函数相比，傅立叶级数的存储成本相对较低，因为它只记录粒子的形状描述符(即傅立叶系数)和位置描述符。因此，作为一种精确的颗粒形态重建方法，傅立叶级数可以为精细的孔隙尺度建模提供良好的基础。虽然傅里叶级数在颗粒重构上展现出良好的数值表现，但是在颗粒堆积时进行的碰撞检测又是一个难以解决的问题。

本发明提出的一种基于傅里叶级数随机生成多孔介质模型的方法及系统，用于对多孔介质模型的重构和生成。

发明内容

本发明旨在解决现有技术的不足，提出一种基于傅里叶级数随机生成多孔介质模型的方法及系统，此方法良好的结合了傅里叶级数在生成不规则颗粒形状上的优势，在此基础上开发了颗粒填充和碰撞检测算法，实现了理论上生成孔隙度接近0％且相互连通的多孔介质。

为实现上述目的，本发明提供了如下方案：

一种基于傅里叶级数随机生成多孔介质模型的方法，包括以下步骤：

S1、对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置。

S2、初始化所述多孔介质模型，并生成颗粒的位置信息；

S3、基于S2的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得所述多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；

S4、对所述离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

S5、对所述颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；

S6、使所述填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

S7、预设循环弹出条件，若S6的判断结果满足所述循环弹出条件，则进行S8；否则，更新所述傅里叶参数，返回S3；

S8、将满足所述循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；

S9、判断S8生成的所述多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回S2。

优选的，所述初始化的参数包括模型孔隙度以及模型尺寸；

所述位置信息包括颗粒中心的位置坐标以及旋转角度。

优选的，所述参数方程的获得方法为：

建立表示孔隙空间的全局坐标系以及表示粒子空间的局部坐标系；

基于所述局部坐标系，展开傅里叶级数，获得单值函数；

基于所述单值函数，获得所述局部坐标系的颗粒轮廓参数方程；

将所述局部坐标系的颗粒轮廓参数方程，转换为所述全局坐标系的颗粒轮廓参数方程；

基于所述全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程。

优选的，所述单值函数的公式为：

其中，r表示粒子中心点P到轮廓上的一个点P'的半径；θ'表示局部坐标系x'-y'的x'轴的半径的极角，0≤θ'<2π；a0、an、bn和N是傅里叶系数。

优选的，所述局部坐标系的颗粒轮廓参数方程如下：

x′(θ′)＝r(θ′)cosθ′

y′(θ′)＝r(θ′)sinθ′

所述全局坐标系的颗粒轮廓参数方程如下：

x＝x_p+x′(θ′)cos(θ₀)＝x_p+r(θ′)cos(θ)

y＝y_p+y′(θ′)sin(θ₀)＝y_p+r(θ′)sin(θ)

其中，x_p和y_p是粒子中心在全局坐标系中的平移距离，变量θ₀是从全局坐标系x轴到局部坐标系x'轴的角度，θ＝θ₀+θ′表示点在全局坐标系中的极角。

优选的，所述碰撞检测的方法为：

将多孔介质对应区域与网格区域进行矩阵相加，判断是否存在异常值；

若存在所述异常值，则存在局部被标记点交叉，实现所述碰撞检测。

优选的，步骤S7中，采用Floyd-Warshall算法，更新所述傅里叶参数；

所述loyd-Warshall算法的迭代方程为：

式中，w_i,j表示从节点i到节点j的路径长度值。

优选的，所述迭代方程的求解方法为：

S71：初始化方程参数，k＝0，对所有节点i和j，令其中，i≠j；若节点i与节点j之间没有弧连接，则w_i，j＝∞；

S72：令k＝k+1，对所有与k节点邻接的流入节点i和流出节点j，若令/>否则令/>

S73：若k＝n，则完成求解；否则，返回S72。

本发明还提供一种基于傅里叶级数随机生成多孔介质模型的系统，包括：模型获取模块、初始化模块、离散模块、边缘提取模块、填充模块、碰撞检测模块、更新模块、模型生成模块以及输出模块；

所述模型获取模块，用于对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置；

所述初始化模块，用于初始化所述多孔介质模型，并生成颗粒的位置信息；所述初始化的参数包括模型孔隙度以及模型尺寸；所述位置信息包括颗粒中心的位置坐标以及旋转角度；

所述离散模块，用于基于所述初始化模块的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得所述多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；所述参数方程的获得过程为：

建立表示孔隙空间的全局坐标系以及表示粒子空间的局部坐标系；

基于所述局部坐标系，展开傅里叶级数，获得单值函数；

基于所述单值函数，获得所述局部坐标系的颗粒轮廓参数方程；

将所述局部坐标系的颗粒轮廓参数方程，转换为所述全局坐标系的颗粒轮廓参数方程；

基于所述全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程；

所述边缘提取模块，用于对所述离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

所述填充模块，用于对所述颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；

所述碰撞检测模块，用于使所述填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

所述更新模块，用于预设循环弹出条件，若所述碰撞检测模块的判断结果满足所述循环弹出条件，则执行所述模型生成模块；否则，更新所述傅里叶参数，返回所述更新模块；

所述模型生成模块，用于将满足所述循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；

所述输出模块，用于判断所述模型生成模块中生成的所述多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回所述初始化模块。

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

本发明提出了一种多孔介质生成方法，该方法利用傅立叶级数在以较低的存储成本描述颗粒形状参数方面的优势，并通过调整傅立叶参数N进行精细建模。例如，可以在高频下描述颗粒表面粗糙度，在低频率下描述颗粒形状。同时，该方法可以使用网格映射将问题从连续转换为离散，这有助于粒子轮廓填充和接触检测。在接触检测方法中，本发明引入了Floyd-Warshall算法来更新可行的参数，而不改变粒子形状描述符。该算法的加入大大提高了该方法的计算速度，最多可以在两个循环中完成单个粒子的重建。

附图说明

为了更清楚地说明本发明的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例一基于傅里叶级数随机生成多孔介质模型的方法流程示意图；

图2是本发明实施例一局部和全局坐标系中不规则粒子轮廓的几何表示；

图3是本发明实施例二参数方程生成的可视化；

图4是本发明实施例二离散的参数方程可视化；

图5是本发明实施例二颗粒边缘提取过程示意图；

图6是本发明实施例二颗粒填充示意图；

图7是本发明实施例二通过基于傅立叶级数的方法显示多孔介质的几何模型；

图8是本发明实施例二从236个生成颗粒中随机选择的50个颗粒与100个天然砂粒的球度和凸度的比较。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例一

如图1所示，一种基于傅里叶级数随机生成多孔介质模型的方法，包括以下步骤：

S1、对预生成的多孔介质模型的孔隙度、分辨率以及模型尺寸进行设置；

S2、生成位置信息包括颗粒中心的位置坐标以及旋转角度。

S3、基于S2的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；

参数方程的获得方法为：

如图2所示，建立表示孔隙空间的全局坐标系x-y以及表示粒子空间的局部坐标系x'-y'；

基于局部坐标系，展开傅里叶级数，获得单值函数；

单值函数的公式为：

其中，r表示粒子中心点P到轮廓上的一个点P'的半径；θ'表示局部坐标系x'-y'的x'轴的半径的极角，0≤θ'<2π；a₀、a_n、b_n和N是傅里叶参数。a₀、a_n和b_n共同决定颗粒尺寸，N是傅立叶级数展开的阶数，其决定颗粒形态。

基于实验统计分析，多孔介质中颗粒尺寸的分布范围和分布模式是已知的(例如，正态、对数正态和双峰正态)。在这些结果的帮助下，决定颗粒尺寸的傅立叶参数(a₀、a_n和b_n)是随机生成的，并遵循分布模式。傅立叶参数N的随机生成启用了粒子形态的设置。对于粒子形态，如果研究对象是粒子形状，则必须在低频区间中随机生成傅立叶参数N。如果研究对象是粒子的表面粗糙度，那么傅立叶参数N必须在高频区间内随机生成。

基于单值函数，获得局部坐标系的颗粒轮廓参数方程；具体的，局部坐标系x'-y'中颗粒轮廓的参数方程也通过在水平和垂直坐标上映射坐标，使用傅立叶级数作为θ′的函数来表达：

局部坐标系的颗粒轮廓参数方程如下：

x′(θ′)＝r(θ′)cosθ′

y′(θ′)＝r(θ′)sinθ′

将局部坐标系的颗粒轮廓参数方程，转换为全局坐标系的颗粒轮廓参数方程；具体的，通过将粒子中心和粒子绕其中心的旋转平移为，将粒子的位置从局部坐标系转换为全局坐标系。

全局坐标系的颗粒轮廓参数方程如下：

x＝x_p+x′(θ′)cos(θ₀)＝x_p+r(θ′)cos(θ)

y＝y_p+y′(θ′)sin(θ₀)＝y_p+r(θ′)sin(θ)

其中，x_p和y_p是粒子中心在全局坐标系中的平移距离，变量θ₀是从全局坐标系x轴到局部坐标系x'轴的角度，θ＝θ₀+θ′表示点在全局坐标系中的极角。通过上述方程，任意的颗粒形状和空间位置信息都可以通过傅立叶参数、平移(x_p和y_p)和旋转(θ₀)来充分描述。

基于全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程。

参数方程按照一定的步长进行分散，步长通常比较小，否则无法确保精度。步长值与网格分辨率相关，以确保每个网格内至少有一个点通过参数方程离散化。

S4、对离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

通过将离散点投影到网格上来提取粒子轮廓，在提取粒子轮廓时，必须仅将预设矩形区域中的网格点与离散化的点进行比较，以提高程序的运行速度。矩形区域的大小可以由x和y方向上离散点的最小值和最大值来确定，包括|x_min|、和/>其中对应的离散点是P₁(x_min，y_P1)、P₂(x_max，y_P2)，P₃(x_p3，y_min)和P₄(x_p4，y_max)。记录这些离散点，将矩形区域划分为四个区域。这些区域的大小如下：

然后遍历四个区域中的每个区域中的网格点，以将这些值与相应的离散点进行比较，如果由网格点a_n(x_n，y_n)和a″_n(x_n+γ_x，y_n-γ_y)界定的区域具有离散点a′_n(x′_n，y′_n)，则该区域被标记为粒子轮廓；否则，则不是。这其中的γ_x和γ_y分别表示在每个网格点的x方向和y方向的尺寸。确定四个区域中是否存在离散点的方法如下：

Region1:x_n≤x′_n&y_n≥y′_n&x_n+γ_x≥x′_n&y_n-γ_y≤y′_n,

Region2:x_n≥x′_n&y_n≥y′_n&x_n-γ_x≤x′_n&y_n-γ_y≤y′_n,

Region3:x_n≤x′_n&y_n≤y′_n&x_n+γ_x≥x′_n&y_n+γ_y≥y′_n,

Region4:x_n≥x′_n&y_n≤y′_n&x_n-γ_x≤x′_n&y_n+γ_y≥y′_n.

S5、对颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；具体的，填充颗粒轮廓以形成单个颗粒。如果标记网格点的x坐标值不连续，则选择步骤S4中沿y方向的标记网格点并单独标记。在这一点上，将基于傅立叶级数方法生成粒子的问题从连续转换为离散，以便于接触检测。

S6、使填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

S7、预设循环弹出条件，若S6的判断结果满足循环弹出条件，则进行S8；否则，更新傅里叶参数，返回S3；

进行碰撞检测以确定生成的颗粒之间是否发生重叠。对于碰撞检测，如果判断结果为真，执行Floyd-Warshall算法，并进入参数方程以更新傅立叶参数，并重复步骤S3-S5。如果结果为假，则执行S8。使用Floyd-Warshall算法，可以及时校正傅立叶参数。最多可以执行两个循环来生成满足接触检测的粒子，从而大大降低了程序的复杂性。

对于Floyd-Warshall算法的嵌入优化了整体模型的计算时间和复杂度。其核心思想是，最短路路径的本质就是比较在两个顶点之间中转点，比较经过与不经过中转点的距离哪个更短。其核心思想是，最短路路径的本质就是比较在两个顶点之间中转点，比较经过与不经过中转点的距离哪个更短。

具体的，碰撞检测的方法为：

将多孔介质对应区域与网格区域进行矩阵相加，判断是否存在异常值；

若存在异常值，则存在局部被标记点交叉，实现碰撞检测。

步骤S7中，采用Floyd-Warshall算法，更新傅里叶参数；

Floyd-Warshall算法的迭代方程为：

式中，w_i,j表示从节点i到节点j的路径长度值。

迭代方程的求解方法为：

S71：初始化方程参数，k＝0，对所有节点i和j，令其中，i≠j；若节点i与节点j之间没有弧连接，则w_i，j＝∞；

S72：令k＝k+1，对所有与k节点邻接的流入节点i和流出节点j，若令/>否则令/>

S73：若k＝n，则完成求解；否则，返回S72。

S8、将满足循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；具体的，存储的参数包括：包括粒子位置坐标、旋转角度和傅立叶参数；并计算孔隙率。

S9、判断S8生成的多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回S2。如果孔隙率值小于设定的目标值，则重复循环。当孔隙率满足设计要求时，输出多孔介质模型的相关信息。

实施例二

为了使本发明所述的基于傅里叶级数随机生成多孔介质模型的方法更加便于理解和接近真实应用，下面根据前人对多孔介质颗粒形状的具体统计，参照前人的结果，生成一个满足测量参数的多孔介质模型。通过测量的颗粒大小的范围设定约束范围，根据测量的颗粒形状参数范围，设置傅里叶参数范围，随后对多孔介质模型进行生成，并将生成数据保存出来。

对于颗粒的尺寸通常通过主轴(d₁)和次轴(d₂)来量化。为了计算d₁和d₂，首先确定粒子的主轴，然后旋转粒子，使主轴遵循全局坐标系，然后再进行计算。d₁和d₂的值由决定颗粒尺寸的傅立叶系数控制。主要尺寸d₁和d₂的比值α被定义为颗粒的纵横比。比率α在0到1之间变化，其中1表示圆形粒子，而接近0的值表示非常细长的粒子。

生成粒子的形状描述符——球度(S)和凸度(C_x)。球度通常用于表示粒子几何体与完美圆的近似程度。凸性可以定义为粒子的面积除以包围粒子的凸壳的面积。

其中A和P分别表示颗粒的面积和周长。A_c是凸壳的面积。

其具体操作步骤如下：

根据对颗粒形状的统计情况，颗粒的长轴的取值范围100-300μm，且呈正态分布，短轴与长轴的比值分布在0.6-1之间，同时根据对傅里叶对颗粒的描述关系，对于颗粒形状的描述一般取傅里叶级数的阶数在4-7之间。因此生成一个分布在5000×3000μm²以及孔隙度为0.24的多孔介质模型作为本发明的实例，并将其的球度和凸度和真实颗粒的进行比较，用来验证程序生成的有效性。

第一步：首先初始化程序，将程序的中的关键参数，孔隙度和模型的物理尺寸分别设定为0.24和5000×3000μm²，以及程序需要的网格分辨率设定为50000×30000。这是为了使的模型的计算足够精细。

第二步：将预计生成的颗粒中心位置设定为0<x<5000和0<y<3000，旋转角度0≤θ₀≤360。该步骤中的颗粒中心位置是随机产生的，当程序进行迭代过程中，该范围变化会及时更新，对被颗粒覆盖的网格单元进行标记，在进行选择生成位置的过程中只会从未被覆盖的颗粒位置进行生成。

第三步，根据傅里叶级数对颗粒表征的特性，设置颗粒的长轴范围分别为100-300，短轴和长轴的比值为0.6-1，上述数据类型为浮点型。设置傅里叶参数的取值范围为4-7，该处的傅里叶参数的数据类型为整型。同时需要引入步骤二中的位置坐标和旋转角度，到这定义好生成的颗粒构型的形状，如图3所示。

第四步，对颗粒构型的参数方程进行离散化，根据离散化的要求，离散的分辨率需比网格的分辨率小，即网格分辨率为50000*30000，表示在5000*3000，最小的单元为0.1，因此在进行离散的时候需要保证在0.1*0.1的网格单元内至少存在1个离散的点，不然会导致该网格内无法搜索到有效的边界单元，导致持续识别边缘错误。在此处设定为5000个点进行离散。由于点过于密集，图4只展示被简化的模型。

第五步，将离散化的参数方程投影到网格中，并且提取颗粒轮廓边缘。调取第二步中生成的位置坐标和旋转角度，将离散的参数方程布置到50000*30000的网格中，因为生成的位置坐标和旋转角度为浮点型，并未对位置信息与网格精度进行转换，在本实施例中，需要将位置进行放大10倍，并进行向下取整，获取在网格中的坐标数值。这样做可以直接获得网格位置的索引，确定位置信息方便快捷。到此，将离散的参数方程置于网格化区域中。在进行边缘检测之前，需要对颗粒的覆盖区域进行标记，标记方法是寻找到离散参数方程中对应的x和y坐标的最大值和最小值，最小值向下取整，最大值向上取整，根据这个操作可以获得一个完全覆盖的离散化参数方程的网格区域。将在此区域的范围内进行边缘检测，提取颗粒轮廓边缘。这样做可以大大提高检测边缘的速度。接着就是进行边缘检测，首先先对标记的网格区域进行划分，划分为四个区域，四个区域的确定方式是根据离散参数方程在X和Y的最大值和最小值，第一个区域为x的最小值和y的最大值确定的，在该区域范围内对网格区域进行由左向右的遍历，判定所在的区域内是否存在参数方程的离散点。该部分的运行嵌入了程序并行的思想，使得计算速度得到很大的提升。在进行边缘检测的时候，对于单个的网格单元，判断在横坐标为x和x+1、纵坐标y和y+1围成的区域内是否存在参数方程的离散点，如果存在，则标记该区域为边界，若不存在，则不是边界。其他的区域分别按照X的最大值和Y的最大值，X的最小值和Y的最小值以及X的最大值和Y的最小值。判定方式相似，只是判定方向有所转变。具体示意图为图5。

第六步：根据第五步提取的边界，进行颗粒填充。分析边界信息，找到对应颗粒边缘y坐标的最大值和最小值。遍历该区域的所有的x坐标，同时进行大小排序，可以获取到在该y坐标下所有被标记的点。随后通过x坐标数组的连续性，获取未被标记点的坐标，进行标记。例如该x坐标的排列顺序为[1,3,4,5,8,9]，就说明在y坐标下存在x坐标为[2,6,7]的网格位置没有被标记，就通过搜索的方式将其标记。特别说明，如果x坐标的数组内只有一个值，则无需标记，直接进入循环。直到可以颗粒边缘所围成的颗粒区域全部被标记。

循环完成。如图6所示

第七步：主要完成对颗粒的碰撞检测。首先利用第六步建立的网格区域，去检查多孔介质区域内的被标记点是否存在交叉。对于这个交叉判定的方法，通过索引到多孔介质对应的区域和第六步建立的区域进行矩阵相加，通过判断是否存在异常值，若存在异常值就说明存在交叉。则将执行第八步Floyd-Warshall算法，对颗粒生成中心即第二步生成的中心坐标点进行距离判断，将中心点与最近的标记点的距离计算。将该数值回传到第三步，重新定义傅里叶参数的随机范围，然后再重新执行步骤4-6。通过这样的操作，该中心坐标的颗粒生成最多两次循环就可以满足生成条件。

第八步：将符合碰撞检测的颗粒位置和标记特征更新到多孔介质模型中，并将结果进行存储。

随后就是重复执行上述步骤，直到多孔介质模型满足预设的孔隙度要求，则结束计算。通过上述计算，完成了在矩形物理域的5000×3000μm²中布置236个不规则颗粒位于，并通过布尔运算提取了流体域，为了验证多孔介质模型的准确性，本发明从236个生成的颗粒中随机选择了50个颗粒，计算了它们的球度和凸度，然后将它们与100个天然砂粒进行比较，获得了良好的比较结果。如图7、图8所示。

与传统的随机叠加算法相比，本发明提出的方法降低了算法的复杂度和迭代序数。基于具有3.4GHz处理器的计算机，使用所提出的方法生成孔隙率＝0.4的多孔介质模型需要大约30秒，而使用传统的随机堆叠算法需要大约3小时。特别是，当生成的模型的孔隙率较低(例如，孔隙率＝0.2)时，传统的算法运行时间变得非常长(超过24小时)，而所提出的方法仅用了2小时就完成了。计算时间差距如此之大的原因是传统算法的大部分计算时间都花在了寻找合适的生成位置上。对于孔隙率小的模型，所提出的方法的执行效率仍然很高，因为所选择的生成位置对于每次执行都是有效的，并且粒子参数是基于接触检测的结果来调整的。此外，当网格分辨率足够大时，生成的模型的理论最小孔隙率趋于零。使用所提出的方法生成具有低孔隙率的模型可以通过提高网格分辨率来实现。增加网格分辨率会导致大量计算消耗。然而，通过分割对模型的生成过程进行了合理的优化。在较小的网格分辨率的情况下，可以排列大粒子的位置，计算结果可以用作下一代的输入，并且可以增加网格分辨率以排列较小的粒子。在理论上，算法可以生成无限接近0％且连通的多孔介质结构。但会造成运算时间的损失。基于本发明的测试，该方法对于孔隙率大于8％的多孔介质模型是高效的，并且适用于生成中渗透和高渗透储层的多孔模型。

实施例三

本发明还提供一种基于傅里叶级数随机生成多孔介质模型的系统，包括：模型获取模块、初始化模块、离散模块、边缘提取模块、填充模块、碰撞检测模块、更新模块、模型生成模块以及输出模块；

模型获取模块，用于对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置；

初始化模块，用于初始化多孔介质模型，并生成颗粒的位置信息；初始化的参数包括模型孔隙度以及模型尺寸；位置信息包括颗粒中心的位置坐标以及旋转角度；

离散模块，用于基于初始化模块的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；参数方程的获得过程为：

建立表示孔隙空间的全局坐标系以及表示粒子空间的局部坐标系；

基于局部坐标系，展开傅里叶级数，获得单值函数；

基于单值函数，获得局部坐标系的颗粒轮廓参数方程；

将局部坐标系的颗粒轮廓参数方程，转换为全局坐标系的颗粒轮廓参数方程；

基于全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程。

边缘提取模块，用于对离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

填充模块，用于对颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；

碰撞检测模块，用于使填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

更新模块，用于预设循环弹出条件，若碰撞检测模块的判断结果满足循环弹出条件，则执行模型生成模块；否则，更新傅里叶参数，返回更新模块；

模型生成模块，用于将满足循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；

输出模块，用于判断模型生成模块中生成的多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回初始化模块。

以上所述的实施例仅是对本发明优选方式进行的描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (9)

1.一种基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，包括以下步骤：

S1、对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置；

S2、初始化所述多孔介质模型，并生成颗粒的位置信息；

S3、基于S2的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得所述多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；

S4、对所述离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

S5、对所述颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；

S6、使所述填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

S7、预设循环弹出条件，若S6的判断结果满足所述循环弹出条件，则进行S8；否则，更新所述傅里叶参数，返回S3；

S8、将满足所述循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；

S9、判断S8生成的所述多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回S2。

2.根据权利要求1所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，

所述初始化的参数包括模型孔隙度以及模型尺寸；

所述位置信息包括颗粒中心的位置坐标以及旋转角度。

3.根据权利要求1所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，所述参数方程的获得方法为：

建立表示孔隙空间的全局坐标系以及表示粒子空间的局部坐标系；

基于所述局部坐标系，展开傅里叶级数，获得单值函数；

基于所述单值函数，获得所述局部坐标系的颗粒轮廓参数方程；

将所述局部坐标系的颗粒轮廓参数方程，转换为所述全局坐标系的颗粒轮廓参数方程；

基于所述全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程。

4.根据权利要求3所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，所述单值函数的公式为：

其中，r表示粒子中心点P到轮廓上的一个点P'的半径；θ'表示局部坐标系x'-y'的x'轴的半径的极角，0≤θ'<2π；a₀、a_n、b_n和N是傅里叶系数。

5.根据权利要求3所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，所述局部坐标系的颗粒轮廓参数方程如下：

x′(θ′)＝r(θ′)cosθ′

y′(θ′)＝r(θ′)sinθ′

所述全局坐标系的颗粒轮廓参数方程如下：

x＝x_p+x′(θ′)cos(θ₀)＝x_p+r(θ′)cos(θ)

y＝y_p+y′(θ′)sin(θ₀)＝y_p+r(θ′)sin(θ)

其中，x_p和y_p是粒子中心在全局坐标系中的平移距离，变量θ₀是从全局坐标系x轴到局部坐标系x'轴的角度，θ＝θ₀+θ′表示点在全局坐标系中的极角；θ′表示局部坐标系x'-y'的x'轴的半径的极角，0≤θ'<2π。

6.根据权利要求1所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，所述碰撞检测的方法为：

将多孔介质对应区域与网格区域进行矩阵相加，判断是否存在异常值；

若存在所述异常值，则存在局部被标记点交叉，实现所述碰撞检测。

7.根据权利要求1所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，步骤S7中，采用Floyd-Warshall算法，更新所述傅里叶参数；

所述Floyd-Warshall算法的迭代方程为：

式中，w_i,j表示从节点i到节点j的路径长度值。

8.根据权利要求7所述的基于傅里叶级数随机生成多孔介质模型的方法，其特征在于，所述迭代方程的求解方法为：

S71：初始化方程参数，k＝0，对所有节点i和j，令其中，i≠j；若节点i与节点j之间没有弧连接，则w_i，j＝∞；

S72：令k＝k+1，对所有与k节点邻接的流入节点i和流出节点j，若令否则令/>

S73：若k＝n，则完成求解；否则，返回S72。

9.一种基于傅里叶级数随机生成多孔介质模型的系统，其特征在于，包括：模型获取模块、初始化模块、离散模块、边缘提取模块、填充模块、碰撞检测模块、更新模块、模型生成模块以及输出模块；

所述模型获取模块，用于对预生成的多孔介质模型的孔隙度、分辨率和尺寸进行设置；

所述初始化模块，用于初始化所述多孔介质模型，并生成颗粒的位置信息；所述初始化的参数包括模型孔隙度以及模型尺寸；所述位置信息包括颗粒中心的位置坐标以及旋转角度；

所述离散模块，用于基于所述初始化模块的设置以及傅里叶级数展开，生成并组合傅里叶参数，获得所述多孔介质颗粒轮廓的参数方程，并进行离散，获得离散的颗粒轮廓数据；所述参数方程的获得过程为：

建立表示孔隙空间的全局坐标系以及表示粒子空间的局部坐标系；

基于所述局部坐标系，展开傅里叶级数，获得单值函数；

基于所述单值函数，获得所述局部坐标系的颗粒轮廓参数方程；

将所述局部坐标系的颗粒轮廓参数方程，转换为所述全局坐标系的颗粒轮廓参数方程；

基于所述全局坐标系的颗粒轮廓参数方程，获得多孔介质颗粒轮廓的参数方程；

所述边缘提取模块，用于对所述离散的颗粒轮廓数据，进行网格映射，提取颗粒轮廓边缘；

所述填充模块，用于对所述颗粒轮廓边缘的局部标记和区域搜索，获得填充颗粒；

所述碰撞检测模块，用于使所述填充颗粒与预设颗粒进行碰撞检测，判断颗粒生成位置的有效性；

所述更新模块，用于预设循环弹出条件，若所述碰撞检测模块的判断结果满足所述循环弹出条件，则执行所述模型生成模块；否则，更新所述傅里叶参数，返回所述更新模块；

所述模型生成模块，用于将满足所述循环弹出条件的颗粒构型添加至模型生成区域，并进行参数存储；

所述输出模块，用于判断所述模型生成模块中生成的所述多孔介质模型是否满足预设生成要求，若满足，则输出多孔介质模型；否则，返回所述初始化模块。