CN109979009B - 一种具有功能梯度的胶结颗粒材料三维细观数字模型重构方法 - Google Patents

Abstract

本发明公开了一种具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，依次通过获取骨料颗粒的STL格式文件、将每个骨料颗粒在三维虚拟空间中重建并将相应的颗粒信息汇总得到骨料文件数据库、调用骨料文件数据库内的骨料颗粒随机投放到指定的样本空间中以得到在骨料颗粒之间填充饱和胶结相的模型、以及利用模拟退火算法对胶结相和随机标记在非骨料颗粒占据点的孔洞相进行重排以得到能量趋于最低的三维细观模型这四个步骤实现；该模型重构方法通过三维数字化复杂颗粒堆积算法，并通过采用两相分布的能量最小原理来生成胶结颗粒材料，其胶凝材料在空间分布上遵守能量最小原理和两相流的物理规律，符合复杂胶结颗粒材料在自然中的实际赋存状态。

Description

一种具有功能梯度的胶结颗粒材料三维细观数字模型重构 方法

技术领域

本发明涉及模型重构技术领域，特别涉及一种具有功能梯度的胶结颗粒材料三维细观数字模型重构方法。

背景技术

随着科技的发展，胶结颗粒材料的相关技术越来越成熟，其已被成功应用于诸如岩土工程、生物工程、海洋工程、冶金工程等重要的工程领域。由于胶结颗粒材料是一个同时涉及到骨料相、胶结相和孔洞三相的复杂系统，所以对于它的物理力学性质的研究是一个巨大的挑战。

近几年，计算机科学技术的迅猛发展，使得数值模拟计算方法成为研究胶结颗粒材料物理力学性质的一个新的强有力工具。目前，胶结颗粒材料三维细观模型的建立主要有以下三种方法：(1)、利用微焦点工业CT扫描技术或者扫描电子显微镜结合数字图像处理技术，在不损伤材料的情况下，将获得的数字图像矢量化，再与相应的数值计算软件结合得到材料的三维细观模型；(2)、蒙特卡罗算法结合MATLAB等可视化软件，将球状颗粒随机分布在一定范围的区域内，形成胶结颗粒材料的三维随机骨料细观模型；(3)、Weibull分布、正态分布等分布函数并结合MATLAB程序，利用统计学原理获得材料的三维随机力学细观模型。以上三种方法都有各自的优势，但也都有很大的局限性。

上述方法(1)构建的模型虽然比较符合实际，但是其制备的设备要求比较高，而且CT图像的识别有一定的精度限制且数字图像处理技术在数据转换的过程中会丢失大量的有效信息，不能完全再现真实模型；而上述方法(2)和方法(3)由于构建过程中采用的蒙特卡罗算法和Weibull分布等分布函数都是以随机数学和统计学的理论来描述胶结颗粒材料的不均匀性，但它们对胶结颗粒材料中胶结相的三维空间分布的预测不能精确反映。因此需开发一种既能够实现骨料相在胶结颗粒材料内部的分布具有一定的功能梯度又可以使胶结相与孔洞的分布呈现能量最小而非随机分布的状态，符合其在自然界中的实际赋存状态。

发明内容

本发明的目的是提供一种基于两相分布最小势能原理实现的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法。

为此，本发明技术方案如下：

一种具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，步骤如下：

S1、采用CT或3D扫描仪对骨料颗粒进行扫描，并利用STL格式将骨料颗粒轮廓离散化为由若干个三角形网片组成的点与面的集合，并采用二进制存储方式或标准ASII码存储方式对STL格式文件进行存储；

S2、在三维虚拟空间内读取骨料颗粒的STL格式文件，使每个骨料颗粒以其三角形面与顶点作为约束条件在其对应的三维虚拟空间内填充，每个三维虚拟空间填充一个骨料颗粒，实现将每个骨料颗粒在三维虚拟空间中重建，汇总每个三维虚拟空间内重建的骨料颗粒的颗粒信息，得到骨料文件数据库；

S3、采用具有功能梯度性的颗粒堆积方法，将已重建过的骨料颗粒随机投放到指定的样本空间中，得到在骨料颗粒之间填充饱和胶结相的模型；

S4、在样本空间中的标记非骨料颗粒占据的点集作为孔洞相，直至满足预定的孔洞饱和度要求，并利用模拟退火算法对孔洞相和胶结相进行重排，得到能量趋于最低的三维细观模型，即骨料模型。

进一步地，步骤S2中，对每个骨料颗粒在三维虚拟空间中重建的具体步骤为：

S201、构建由规则空间点集组成的三维虚拟空间，该虚拟空间的尺寸大于所有待填充的骨料颗粒的尺寸即可；

S202、在三维虚拟空间内，将骨料颗粒的STL格式文件按照二进制形式或者ASII码形式的写入规则读出骨料颗粒的所有三角面及其对应的顶点，并将全部空间点作为待查询点利用自待查询点发射射线的方式判断每个待查询点是否在三角形内；

S203、根据步骤S202将所有的待查询点分为位于骨料颗粒内部和位于骨料颗粒外部两类，进而在三维虚拟空间中标记出全部位于骨料颗粒内部的待查询点；

S204、遍历三维虚拟空间中所有待查询点，依次判断每一个待查询点与其上、下、左、右、前、后六个面接触的点中是否存在有在步骤S203中被标记为位于内部的点，且该点并非六个面全部与被标记为位于内部的点相接触；采用不同的标记方式标记出全部符合上述条件的待查询点，作为外包络面的点。这样，全部的外包络面的点即构成每个骨料颗粒的外包络面，即实现将每个骨料颗粒用一个不同属性的外包络面包裹起来；

S205、将由每个骨料颗粒的坐标阵列及其阵列长度，和外包络面的坐标阵列以及骨料颗粒的中心坐标组成的骨料颗粒信息保存到骨料文件数据库中。

更进一步地，步骤S202中，判断每个待查询点是否在三角形内的具体方法为：

步骤1)将步骤S201中构建形成点阵形式的三维虚拟空间的所有的点作为待查询点，依次地，由待查询点Ο发射的射线R(t)可以表示为：

R(t)＝Ο+tD 0≤t≤∞

其中，t为待查询点Ο距射线与三角面的交点之间的距离，Ο为待查询点，D为射线R(t)的方向向量。

三角面中的任意一点均可表示为：

T(μ,υ)＝(1-μ-υ)V₀+μV₁+υV₂

其中，(μ、υ)被称为点T(μ,υ)的重心坐标参数，且满足0≤μ，0≤υ，υ+μ≤1；V₀、V₁、V₂分别为三角形面的三个坐标。

计算射线R(t)与三角形面的交点，只需联立公式(1)、(2)且满足，R(t)＝T(μ,υ)即：

Ο+tD＝(1-μ-υ)V₀+μV₁+υV₂

变换后，可得：

求解此线性方程组可得t、μ、υ，且当μ、υ满足(2)中条件时，交点在三角形面内，也即射线与三角形面相交；

步骤2)统计每个待查询点发出的射线与三角形面相交的交点个数，判断相应待查询点是否在三角形内；

当从待查询点发射的任意方向的射线与三角面的交点个数为奇数时，则该待查询点位于骨料颗粒内部；

当从待查询点发射的任意方向的射线与三角面的交点个数为偶数时，则该待查询点位于骨料颗粒外部。

进一步地，该方法的步骤S3的具体操作方法为：

S301、构建由规则空间点集组成的样本空间Ω，样本空间的尺寸大于全部骨料颗粒的体积总和；一般情况下，虚拟空间的体积一般设定为所有待填充的骨料颗粒的体积的2倍；

S302、在样本空间Ω的点集中任取一点作为颗粒填充的标志点，其坐标为M(X_i,Y_i,Z_i)，通过位置函数F(X_i,Y_i,Z_i),计算出此标志点对应的颗粒体积参数V_pi；

S303、在步骤S2建立的骨料文件数据库中调用一个与标志点M匹配的颗粒体积参数V_pi相对应的颗粒P_i，通过该标志点将颗粒P_i布置到Ω立方区域内，并进行随机的旋转变换；

S304、颗粒接触判断：以标志点M为中心建立一个略大于被填充颗粒体积的立方体空间，遍历该立方体空间中每一点：

当没有检测到颗粒的外包络面，向样本空间中投放该颗粒；

当检测到外包络面且没有检测到骨料颗粒，向样本空间中投放该颗粒；

当既检测到外包络面又检测到骨料颗粒，则放弃向样本空间中投放该颗粒；

S305、重复上述步骤S302～S304，直至达到所需要的充填密度或者指定空间区域中无法充填入新的颗粒为止，此时以外包络面为基准，将样本空间内除骨料颗粒占据空间外的部分定义为胶结相，计算填充率以判断该颗粒堆积模型是否满足预设的填充要求：

当填充率满足目标填充率要求时，输出颗粒堆积模型；

当填充率低于目标填充率时，继续步骤S306；

S306、采用上下移动方式、左右移动方式、中心移动方式中的一种或多种对骨料颗粒进行至少一次重排，并通过剔除重排后的边界空间不断提高填充率，直至达到目标填充率的要求。

更进一步地，在步骤S303中，旋转变换采用骨料颗粒的各个三角形面作为处理单元，将每个三角形面的顶点坐标扩展为四行一列的矩阵，并通过在其前左乘一个相应的四阶变换矩阵，以表示颗粒变化，具体矩阵公式如下：

(i)当颗粒旋转变换为绕X轴的旋转变换时：

(ii)当颗粒旋转变换为绕Y轴的旋转变换时：

(iii)当颗粒旋转变换为绕Z轴的旋转变换时：

更进一步地，步骤S306的具体步骤为：

步骤1)在模型中插入一个基准面或标记一个中心点，记录模型中每个颗粒重心距离该基准面的垂直间距或距离中心点的间距，并且把所有颗粒按照距离的大小排序并且编号；

步骤2)设置循环移动次数的控制参数i，使一旦循环次数超过i，则自动跳出循环；

步骤3)按照间距由近及远，依次地将骨料颗粒朝向基准面或中心点的方向移动，当全部骨料颗粒依次地完成一次移动即视为完成一次循环移动；其中，骨料颗粒的移动应符合：

I、每个骨料颗粒每一次的移动距离为一个单位，即点集中相邻两个点之间的间距，当判断出骨料颗粒没有移动一个单位的空间时，则该骨料颗粒不发生移动，继续对下一个骨料颗粒进行移动；

II、在每一次移动结束后判断已经发生的移动次数，及步骤2)所述的循环移动次数，当该骨料颗粒的循环移动次数等于i时，该骨料颗粒不再发生移动；

III、对于任意骨料颗粒来说，其移动的终点为模型中插入的水平面，即当骨料颗粒向上移动或向下移动至水平面时，该骨料颗粒不在发生移动；

步骤4)剔除每一次循环移动后由于颗粒移动而产生的顶部空间和底部空间，使全部骨料颗粒的占据空间减小，计算填充率是否达到目标填充率，若仍低于目标填充率，则返回步骤3)进行下一次循环移动，直至填充率满足目标填充率。

进一步地，该方法的步骤S4的具体步骤如下：

S401、利用蒙特卡罗算法随机将样本空间中的除骨料颗粒占据的点集以外的点标记为孔洞相，直至骨料模型中孔洞相的点满足预定的孔洞饱和度；

S402、将整个胶结颗粒料三维模型自下而上划分为若干区域；

S403、对于每个区域，分别用0、1、-1分别代表骨料相、胶结相和孔洞相体素的空间位置，得到数字三维矩阵；随机交换数字三维矩阵中胶结相和孔洞相两个体素的位置，使胶结颗粒材料的状态发生变化，计算每一次胶结颗粒材料状态变化后的界面能E，进而得到能量变化ΔE：

当ΔE≤0时，系统会自动接受这个新的能量变化；

当ΔE＞0时，系统会按照一定的概率来接受新的能量变化；

其中，每个体素迭代的次数由迭代参数N控制；接受概率P由基于固体退火过程中温度的变化并结合Monte Carlo随机方法所形成的Metropolis准则产生；其中，

E_λ是由“冷却进度表”给出的参考能量，E_λ＝μ^mE⁰，μ是冷却参数，m是马尔科夫链的数量，E⁰是初始能量；

当两个马尔科夫链之间的能量变化满足：E^m-E^m-1/E^m-1＜10^-7时，认为胶结颗粒材料内三相分布是处于平衡状态；继续迭代若干次，并且平均每个马尔科夫链的界面能，使其产生光滑的交界面。

与现有技术相比，该具有功能梯度的胶结颗粒材料三维细观数字模型重构方法通过三维数字化复杂颗粒堆积算法，并通过采用两相分布的能量最小原理来生成胶结颗粒材料，其胶凝材料在空间分布上遵守能量最小原理和两相流的物理规律，符合复杂胶结颗粒材料在自然中的实际赋存状态。

附图说明

图1为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中颗粒文件数据库建立的流程图；

图2为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中颗粒排布过程的示意图；

图3为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中颗粒堆积过程的流程图；

图4为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中颗粒移动的流程图；

图5为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中胶结颗粒材料内部细观三相分布图；

图6为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中胶结颗粒材料三维细观结构建立的流程图；

图7(a)为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中引入空隙后的初始模型结构示意图；

图7(b)为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中引入空隙后的过渡模型结构示意图；

图7(c)为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法中引入空隙后的最终模型结构示意图；

图8为本发明的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法的流程图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

如图8所示，该具有功能梯度的胶结颗粒材料三维细观数字模型重构方法的具体操作步骤如下：

S1、采用CT或3D扫描仪对骨料颗粒进行扫描，并利用STL格式将骨料颗粒轮廓离散化为由若干个三角形网片组成的点与面的集合；为了便于后续读取，将上述得到的STL格式文件通过二进制存储方式或标准ASII码存储方式将点与面进行有规则的存储；其中，二进制存储格式具有存储文件小，占用内存少等优点；而ASII码存储格式具有读、写操作简单，信息表达更加直观等特点；

S2、在三维虚拟空间内读取骨料颗粒的STL格式文件，使每个骨料颗粒以其三角形面与顶点作为约束条件在其对应的三维虚拟空间内填充，每个三维虚拟空间填充一个骨料颗粒，实现将每个骨料颗粒在三维虚拟空间中重建，汇总每个三维虚拟空间内重建的骨料颗粒的颗粒信息，得到骨料文件数据库；

该步骤为后续将骨料颗粒投放到样本空间中进行颗粒排布做准备，同时也使得投放的颗粒具备二次操作的特点。

如图1所示，每个骨料颗粒在三维虚拟空间中重建的具体步骤如下：

S201、构建由规则空间点集组成的三维虚拟空间，该虚拟空间的尺寸大于所有待填充的骨料颗粒的尺寸即可；

S202、在三维虚拟空间内，将骨料颗粒的STL格式文件按照二进制形式或者ASII码形式的写入规则读出骨料颗粒的所有三角面及其对应的顶点，并将全部空间点作为待查询点利用自待查询点发射射线的方式判断每个待查询点是否在三角形内；

该步骤由于是将骨料颗粒在三维虚拟空间中重建，因此该过程的关键是准确的判断出位于骨料颗粒内部的空间点，其具体步骤如下：

步骤1)将步骤S201中构建形成点阵形式的三维虚拟空间的所有的点作为待查询点，依次地，由待查询点Ο发射的射线R(t)可以表示为：

R(t)＝Ο+tD 0≤t≤∞ 式(1)

其中，t为待查询点Ο距射线与三角面的交点之间的距离，Ο为待查询点，D为射线R(t)的方向向量。

三角面中的任意一点均可表示为：

T(μ,υ)＝(1-μ-υ)V₀+μV₁+υV₂ 式(2)

其中，(μ、υ)被称为点T(μ,υ)的重心坐标参数，且满足0≤μ，0≤υ，υ+μ≤1；V₀、V₁、V₂分别为三角形面的三个坐标。

计算射线R(t)与三角形面的交点，只需联立公式(1)、(2)且满足，R(t)＝T(μ,υ)即：

Ο+tD＝(1-μ-υ)V₀+μV₁+υV₂ 式(3)

变换后，可得：

求解此线性方程组可得t、μ、υ，且当μ、υ满足(2)中条件时，交点在三角形面内，也即射线与三角形面相交；

步骤2)统计每个待查询点发出的射线与三角形面相交的交点个数，判断相应待查询点是否在三角形内；

当从待查询点发射的任意方向的射线与三角面的交点个数为奇数时，则该待查询点位于骨料颗粒内部；

当从待查询点发射的任意方向的射线与三角面的交点个数为偶数时，则该待查询点位于骨料颗粒外部；

S203、根据步骤S202将所有的待查询点分为位于骨料颗粒内部和位于骨料颗粒外部两类，进而在三维虚拟空间中标记出全部位于骨料颗粒内部的待查询点；

S204、遍历三维虚拟空间中所有待查询点，依次判断每一个待查询点与其上、下、左、右、前、后六个面接触的点中是否存在有在步骤S203中被标记为位于内部的点，且该点并非六个面全部与被标记为位于内部的点相接触；采用不同的标记方式标记出全部符合上述条件的待查询点，作为外包络面的点。这样，全部的外包络面的点即构成每个骨料颗粒的外包络面，即实现将每个骨料颗粒用一个不同属性的外包络面包裹起来；该步骤S204得到的外包络面能够有效避免骨料颗粒在样本空间中进行骨料颗粒排布时发生骨料颗粒相互接触碰撞的问题；

S205、将骨料颗粒信息包括每个骨料颗粒的坐标阵列及其阵列长度(即像素意义下的颗粒体积)，和外包络面的坐标阵列以及骨料颗粒的中心坐标统一保存到骨料文件数据库中，便于后续调用；

S3、采用具有功能梯度性的颗粒堆积方法，将已重建过的骨料颗粒随机投放到指定的样本空间中。颗粒堆积的流程图如图4所示，其具体步骤如下：

S301、构建由规则空间点集组成的样本空间Ω，样本空间Ω为一个立方形区域，其体积为全部骨料颗粒体积的2倍；在本实施例中，由于实际的骨料颗粒的最大直径为10～15cm，因此样本空间Ω的空间点集的单位距离为1cm，即相邻两个空间点之间的间距为1cm；

由于样本空间是目标模型的重建空间，其本质上是由规则的点集构成的空间，因此构建样本空间就是构建一系列规则的空间点阵；此外，该步骤中的样本空间与步骤S201中构建的三维虚拟空间没有关系，步骤S201中的虚拟空间只起到骨料重建的作用，当骨料颗粒数据库建立完成以后就不再需要三维虚拟空间。样本空间是进行骨料颗粒填充的场所。

S302、在样本空间Ω的点集中任取一点作为颗粒填充的标志点，其坐标为M(X_i,Y_i,Z_i)，通过位置函数F(X_i,Y_i,Z_i),计算出此标志点对应的颗粒体积参数V_pi；在该步骤中，位置函数F(Xi,Yi,Zi)为包含颗粒体积参数V_pi的函数，其具体形式与填充时选择的标志点的空间位置有关，是由其三维坐标组成的分布函数；体积参数V_pi是根据位置函数F计算出的与待填充颗粒的体积相对应的系数。由于在颗粒堆积时采用了位置函数F(X_i,Y_i,Z_i)和体积参数V_pi控制颗粒充填时的规律，使得此时生成的模型中颗粒充填的类型与其所在的空间位置具有直接的关系，由此可以使得材料具有一定的功能梯度。由于未在颗粒堆积时引入孔洞(气相)，故而此时的颗粒堆积模型是饱和的具有一定功能梯度的胶结颗粒料堆积模型。

S303、在步骤S2建立的骨料文件数据库中调用一个与标志点M匹配的颗粒体积参数V_pi相对应的颗粒P_i，通过该标志点将颗粒P_i布置到Ω立方区域内，并进行随机的旋转变换；其中，与标志点M匹配的颗粒体积参数V_pi相对应的颗粒P_i是指：根据选定的标志点M计算出体积参数V_pi，然后在颗粒形状数据库中寻找一个具有相同或最相近似的V_pi的颗粒P_i，并将该颗粒P_i投放至标志点M处；

在对任意骨料颗粒进行随机的旋转变换时，将骨料颗粒的各个三角形面作为处理单元，将每个三角形面的顶点坐标扩展为四行一列的矩阵，并通过在其前左乘一个相应的四阶变换矩阵，以表示颗粒变化，具体矩阵公式如下式(5)所示：

当颗粒变化为旋转变换时，颗粒的旋转变换将分解为针对三个坐标轴的变换，包括绕X轴的旋转变换、绕Y轴的旋转变换和绕Z轴的旋转变换：

(i)当颗粒旋转变换为绕X轴的旋转变换时：

(ii)当颗粒旋转变换为绕Y轴的旋转变换时：

(iii)当颗粒旋转变换为绕Z轴的旋转变换时：

S304、颗粒接触判断：以标志点M为中心建立一个略大于被填充颗粒体积的立方体空间，遍历该立方体空间中每一点：

当没有检测到颗粒的外包络面，向样本空间中投放该颗粒；

当检测到外包络面且没有检测到骨料颗粒，向样本空间中投放该颗粒；

当既检测到外包络面又检测到骨料颗粒，则放弃向样本空间中投放该颗粒；

S305、重复上述步骤S302～S304，直至达到所需要的充填密度或者指定空间区域中无法充填入新的颗粒为止，然后此时以外包络面为基准，将样本空间内除骨料颗粒占据空间外的部分定义为胶结相，计算填充率以判断该颗粒堆积模型是否满足预设的填充要求：

当填充率满足目标填充率要求时，步骤S3结束，输出如图2所示的颗粒堆积模型；

当填充率低于目标填充率时，继续步骤S306；

S306、由于在上述不规则颗粒填充过程中没有考虑颗粒形状及随机函数选择的随机充填位置对模型整体填充率的影响，所以填充后的充填率往往达不到要求，因此当经过步骤S305输出的颗粒堆积模型的填充率低于目标填充率时，利用该颗粒移动的方法进一步重排，以提高模型的填充率。

该步骤S306的具体步骤如下：

步骤1)在模型中插入一个水平面，记录模型中每个颗粒重心距离该水平面的竖直距离，并且把所有颗粒按照距离的大小排序并且编号；

步骤2)设置循环移动次数的控制参数i，使一旦循环次数超过i，则自动跳出循环；

步骤3)对于位于水平面以上的颗粒，按照距离该水平面的远近，由近及远依次向下移动骨料颗粒；对于选定水平面以下的颗粒，按照距离该水平面的远近，由近及远依次向上移动颗粒；当全部骨料颗粒依次地完成一次移动即视为完成一次循环移动；

在步骤3)中，骨料颗粒的具体的移动过程应符合以下三个要求：

I、每个骨料颗粒每一次的移动距离为一个单位，即点集中相邻两个点之间的间距，当判断出骨料颗粒没有移动一个单位的空间时，则该骨料颗粒不发生移动，继续对下一个骨料颗粒进行移动；

II、在每一次移动结束后判断已经发生的移动次数，及步骤2)所述的循环移动次数，当该骨料颗粒的循环移动次数等于i时，该骨料颗粒不再发生移动；

III、对于任意骨料颗粒来说，其移动的终点为模型中插入的水平面，即当骨料颗粒向上移动或向下移动至水平面时，该骨料颗粒不在发生移动；

步骤4)剔除每一次循环移动后由于颗粒移动而产生的顶部空间和底部空间，使全部骨料颗粒的占据空间减小，计算填充率是否达到目标填充率，若仍低于目标填充率，则进行下一次循环移动，直至填充率满足目标填充率。

具体的颗粒移动流程图如图5所示。

为了增加模型材料的各向异性，步骤S306的颗粒移动的方式还包括左右移动方式和中心移动方式。

当采用左右移动方式时，步骤S306的具体步骤如下：

步骤1)在模型中插入一个竖直面，记录模型中每个颗粒重心距离该竖直面的水平距离，并且把所有颗粒按照距离的大小排序并且编号；

步骤2)设置循环移动次数的控制参数i，使一旦循环次数超过i，则自动跳出循环；

步骤3)对于位于竖直面左侧的颗粒，按照距离该竖直面的远近，由近及远依次向右移动骨料颗粒；对于选定竖直面右侧的颗粒，按照距离该水平面的远近，由近及远依次向左移动颗粒；当全部骨料颗粒依次地完成一次移动即视为完成一次循环移动；

步骤4)剔除每一次循环移动后由于颗粒移动而产生的左侧空间和右侧空间，使全部骨料颗粒的占据空间减小，计算填充率是否达到目标填充率，若仍低于目标填充率，则进行下一次循环移动，直至填充率满足目标填充率。

当采用中心移动方式时，步骤S306的具体步骤如下：

步骤1)在模型中选择一个点作为颗粒移动的中心点，记录每个颗粒重心的位置距离此中心点的距离，并把所有颗粒按照距离的大小排序并且编号；

步骤2)设置循环移动次数的控制参数i，使一旦循环次数超过i，则自动跳出循环；

步骤3)每个颗粒均按照距离该中心点的远近，由近及远按照朝向中心点的方向依次移动骨料颗粒；当全部骨料颗粒依次地完成一次移动即视为完成一次循环移动；

步骤4)剔除每一次循环移动后由于颗粒移动而产生的顶部、底部、左侧和右侧的空间，使全部骨料颗粒的占据空间减小，计算填充率是否达到目标填充率，若仍低于目标填充率，则进行下一次循环移动，直至填充率满足目标填充率。

以上三种颗粒移动方式既可以单独使用，也可以交叉使用。具体的操作流程可以如图4所示。上述颗粒移动是按照距离远近有秩序的移动，从模型整体的角度看，模型的功能梯度性并没有因为移动而发生实质性的变化。

在完成上述步骤S3后，全部骨料颗粒排布完成，此时骨料模型可以视为一个在骨料颗粒之间填充饱和胶结材料的模型。

S4、采用模拟退火法在模型的胶结相中引入孔洞，即采用能量总是趋于最低的最小能量原理来模拟胶结颗粒材料的细观三维结构；该模拟退火算法相比于其他全局优化方法，最大的优点是它以一定的概率接受在优化过程中出现的不满足其优化条件的解，这样就能充分避免把局部最优解当成全局最优解的问题。

由于在步骤S3中，骨料模型是以规则空间点集构成的样本空间为基础生成的，因此除骨料颗粒占据的点集以外，剩余的点集则应由胶结料和孔洞按照饱和度按界面能最小原理进行分配。因此引入孔洞的具体方法是随机选择未被骨料占据的点，并且赋予相应的属性，按照孔洞饱和度的要求选择随机分配点的数目，被选中的点集即为孔洞，剩余未被选中的点即为胶结相，得到胶结颗粒材料的初始三相分布，如图3所示。

该步骤S4的具体步骤如下：

S401、利用蒙特卡罗算法随机将样本空间中的除骨料颗粒占据的点集以外的点标记为孔洞相，直至骨料模型中孔洞相的点满足预定的孔洞饱和度为止；

S402、将整个胶结颗粒料三维模型自下而上划分为若干区域，通过在每个区域中分别以不同的饱和度进行模拟退火算法优化，使整个模型由若干个饱和度不同的区域组成，模型三相分布的不均匀性得到了很大的加强；

S403、用0、1、-1分别代表骨料相、胶结相和孔洞相体素的空间位置，得到数字三维矩阵；随机交换数字三维矩阵中胶结相和孔洞相两个体素的位置，使胶结颗粒材料的状态发生变化，计算每一次胶结颗粒材料状态变化后的界面能E；界面能E的计算方法如下式(10)所示：

其中，i表示任意体素的空间位置，j表示距离i最近的所有体素的空间位置，如果体素位置在相位K中，则

否则为0；/>

表示距离j体素最近的体素位置在/>

相中/>

取1，否则为0；/>

是不同相之间接触的界面自由能，如下式(10)所示：/>

其中，饱和度Sr为胶结相体积与剔除骨料后模型总体积的比值，其表示不同胶凝材料含量的三维细观结构模型；由于胶结颗粒材料的每个状态都会对应于某个特定的界面能，当材料总体趋于稳定时，系统的总界面能最小；因此在该过程中，通过随机交换数字三维矩阵中胶结相和孔洞相中任意两个体素的位置，进而获得该模型系统的总界面能最小时骨料颗粒、胶结相和孔洞的位置分布方式。

在上述模拟的过程中，两个随机选择的体素相互交换位置，系统的能量会产生一定的变化ΔE：

当ΔE≤0时，系统会自动接受这个新的能量变化；

当ΔE＞0时，系统会按照一定的概率来接受新的能量变化。

每个体素迭代的次数由迭代参数N控制；其中，接受概率P是由基于固体退火过程中温度的变化并结合Monte Carlo随机方法所形成的Metropolis准则产生，如公式(12)所示：

其中，E_λ是由“冷却进度表”给出的参考能量，可表示为：

E_λ＝μ^mE⁰ 式(12)

其中，μ是冷却参数，m是马尔科夫链的数量，E⁰是初始能量；

当两个马尔科夫链之间的能量变化满足：E^m-E^m-1/E^m-1＜10^-7时，认为胶结颗粒材料内三相分布是处于平衡状态。

为了进一步消除由于模拟退火算法中局部搜索优化的随机性导致的胶结相分布呈锯齿状的影响，再继续迭代若干次，并且平均每个马尔科夫链的界面能，使其产生光滑的交界面。

胶结颗粒材料三维细观模型重构流程图如图6所示。

本实施例的胶结颗粒材料经过初步模拟退火优化以后的两相流分布过程示意图如图7(a)～图7(c)所示，其中深灰色代表骨料相，浅灰色代表胶结相，白色代表孔洞，模型的具体参数如表1。

表1胶结颗粒材料模型具体参数：

解析度	饱和度	充填率	颗粒形状
				100*100*100	0.8	0.60	球形

Claims (5)

1.一种具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，其特征在于，步骤如下：

S1、采用CT或3D扫描仪对骨料颗粒进行扫描，并利用STL格式将骨料颗粒轮廓离散化为由若干个三角形网片组成的点与面的集合，并采用二进制存储方式或标准ASII码存储方式对STL格式文件进行存储；

S2、在三维虚拟空间内读取骨料颗粒的STL格式文件，使每个骨料颗粒以其三角形面与顶点作为约束条件在其对应的三维虚拟空间内填充，每个三维虚拟空间填充一个骨料颗粒，实现将每个骨料颗粒在三维虚拟空间中重建，汇总每个三维虚拟空间内重建的骨料颗粒的颗粒信息，得到骨料文件数据库；

在步骤S2中，对每个骨料颗粒在三维虚拟空间中重建的具体步骤为：

S201、构建由规则空间点集组成的三维虚拟空间，该虚拟空间的尺寸大于所有待填充的骨料颗粒的尺寸即可；

S202、在三维虚拟空间内，将骨料颗粒的STL格式文件按照二进制形式或者ASII码形式的写入规则读出骨料颗粒的所有三角面及其对应的顶点，并将全部空间点作为待查询点利用自待查询点发射射线的方式判断每个待查询点是否在三角形内；

S203、根据步骤S202将所有的待查询点分为位于骨料颗粒内部和位于骨料颗粒外部两类，进而在三维虚拟空间中标记出全部位于骨料颗粒内部的待查询点；

S204、遍历三维虚拟空间中所有待查询点，依次判断每一个待查询点与其上、下、左、右、前、后六个面接触的点中是否存在有在步骤S203中被标记为位于内部的点，且该点并非六个面全部与被标记为位于内部的点相接触；采用不同的标记方式标记出全部符合上述条件的待查询点，作为外包络面的点；这样，全部的外包络面的点即构成每个骨料颗粒的外包络面，即实现将每个骨料颗粒用一个不同属性的外包络面包裹起来；

S205、将由每个骨料颗粒的坐标阵列及其阵列长度，和外包络面的坐标阵列以及骨料颗粒的中心坐标组成的骨料颗粒信息保存到骨料文件数据库中；

S3、采用具有功能梯度性的颗粒堆积方法，将已重建过的骨料颗粒随机投放到指定的样本空间中，得到在骨料颗粒之间填充饱和胶结相的模型；

步骤S3的具体操作方法为：

S301、构建由规则空间点集组成的样本空间Ω，样本空间的尺寸大于全部骨料颗粒的体积总和；

S302、在样本空间Ω的点集中任取一点作为颗粒填充的标志点，其坐标为M(X_i,Y_i,Z_i)，通过位置函数F(X_i,Y_i,Z_i),计算出此标志点对应的颗粒体积参数V_pi；

S303、在步骤S2建立的骨料文件数据库中调用一个与标志点M匹配的颗粒体积参数V_pi相对应的颗粒P_i，通过该标志点将颗粒P_i布置到Ω立方区域内，并进行随机的旋转变换；

S304、颗粒接触判断：以标志点M为中心建立一个略大于被填充颗粒体积的立方体空间，遍历该立方体空间中每一点：

当没有检测到颗粒的外包络面，向样本空间中投放该颗粒；

当检测到外包络面且没有检测到骨料颗粒，向样本空间中投放该颗粒；

当既检测到外包络面又检测到骨料颗粒，则放弃向样本空间中投放该颗粒；

S305、重复上述步骤S302～S304，直至达到所需要的充填密度或者指定空间区域中无法充填入新的颗粒为止，此时以外包络面为基准，将样本空间内除骨料颗粒占据空间外的部分定义为胶结相，计算填充率以判断该颗粒堆积模型是否满足预设的填充要求：

当填充率满足目标填充率要求时，输出颗粒堆积模型；

当填充率低于目标填充率时，继续步骤S306；

S306、采用上下移动方式、左右移动方式、中心移动方式中的一种或多种对骨料颗粒进行至少一次重排，并通过剔除重排后的边界空间不断提高填充率，直至达到目标填充率的要求；

S4、在样本空间中的标记非骨料颗粒占据的点集作为孔洞相，直至满足预定的孔洞饱和度要求，并利用模拟退火算法对孔洞相和胶结相进行重排，得到能量趋于最低的三维细观模型，即骨料模型。

2.根据权利要求1所述的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，其特征在于，在步骤S202中，判断每个待查询点是否在三角形内的具体方法为：

步骤1)将步骤S201中构建形成点阵形式的三维虚拟空间的所有的点作为待查询点，依次地，由待查询点Ο发射的射线R(t)可以表示为：

R(t)＝Ο+tD 0≤t≤∞

其中，t为待查询点Ο距射线与三角面的交点之间的距离，Ο为待查询点，D为射线R(t)的方向向量；

三角面中的任意一点均可表示为：

T(μ,υ)＝(1-μ-υ)V₀+μV₁+υV₂

其中，(μ、υ)被称为点T(μ,υ)的重心坐标参数，且满足0≤μ，0≤υ，υ+μ≤1；V₀、V₁、V₂分别为三角形面的三个坐标；

计算射线R(t)与三角形面的交点，只需联立公式(1)、(2)且满足，R(t)＝T(μ,υ)即：

Ο+tD＝(1-μ-υ)V₀+μV₁+υV₂

变换后，可得：

求解此线性方程组可得t、μ、υ，且当μ、υ满足(2)中条件时，交点在三角形面内，也即射线与三角形面相交；

步骤2)统计每个待查询点发出的射线与三角形面相交的交点个数，判断相应待查询点是否在三角形内；

当从待查询点发射的任意方向的射线与三角面的交点个数为奇数时，则该待查询点位于骨料颗粒内部；

当从待查询点发射的任意方向的射线与三角面的交点个数为偶数时，则该待查询点位于骨料颗粒外部。

3.根据权利要求1所述的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，其特征在于，在步骤S303中，旋转变换采用骨料颗粒的各个三角形面作为处理单元，将每个三角形面的顶点坐标扩展为四行一列的矩阵，并通过在其前左乘一个相应的四阶变换矩阵，以表示颗粒变化，具体矩阵公式如下：

(i)当颗粒旋转变换为绕X轴的旋转变换时：

(ii)当颗粒旋转变换为绕Y轴的旋转变换时：

(iii)当颗粒旋转变换为绕Z轴的旋转变换时：

4.根据权利要求1所述的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，其特征在于，步骤S306的具体步骤为：

步骤1)在模型中插入一个基准面或标记一个中心点，记录模型中每个颗粒重心距离该基准面的垂直间距或距离中心点的间距，并且把所有颗粒按照距离的大小排序并且编号；

步骤2)设置循环移动次数的控制参数i，使一旦循环次数超过i，则自动跳出循环；

步骤3)按照间距由近及远，依次地将骨料颗粒朝向基准面或中心点的方向移动，当全部骨料颗粒依次地完成一次移动即视为完成一次循环移动；其中，骨料颗粒的移动应符合：

I、每个骨料颗粒每一次的移动距离为一个单位，即点集中相邻两个点之间的间距，当判断出骨料颗粒没有移动一个单位的空间时，则该骨料颗粒不发生移动，继续对下一个骨料颗粒进行移动；

II、在每一次移动结束后判断已经发生的移动次数，及步骤2)所述的循环移动次数，当该骨料颗粒的循环移动次数等于i时，该骨料颗粒不再发生移动；

III、对于任意骨料颗粒来说，其移动的终点为模型中插入的水平面，即当骨料颗粒向上移动或向下移动至水平面时，该骨料颗粒不在发生移动；

步骤4)剔除每一次循环移动后由于颗粒移动而产生的顶部空间和底部空间，使全部骨料颗粒的占据空间减小，计算填充率是否达到目标填充率，若仍低于目标填充率，则返回步骤3)进行下一次循环移动，直至填充率满足目标填充率。

5.根据权利要求1所述的具有功能梯度的胶结颗粒材料三维细观数字模型重构方法，其特征在于，步骤S4的具体步骤如下：

S401、利用蒙特卡罗算法随机将样本空间中的除骨料颗粒占据的点集以外的点标记为孔洞相，直至骨料模型中孔洞相的点满足预定的孔洞饱和度；

S402、将整个胶结颗粒料三维模型自下而上划分为若干区域；

S403、对于每个区域，分别用0、1、-1分别代表骨料相、胶结相和孔洞相体素的空间位置，得到数字三维矩阵；随机交换数字三维矩阵中胶结相和孔洞相两个体素的位置，使胶结颗粒材料的状态发生变化，计算每一次胶结颗粒材料状态变化后的界面能E，进而得到能量变化ΔE：

当ΔE≤0时，系统会自动接受这个新的能量变化；

当ΔE>0时，系统会按照一定的概率来接受新的能量变化；

其中，每个体素迭代的次数由迭代参数N控制；接受概率P由基于固体退火过程中温度的变化并结合Monte Carlo随机方法所形成的Metropolis准则产生；其中，

E_λ是由“冷却进度表”给出的参考能量，E_λ＝μ^mE⁰，μ是冷却参数，m是马尔科夫链的数量，E⁰是初始能量；

当两个马尔科夫链之间的能量变化满足：E^m-E^m-1/E^m-1＜10^-7时，认为胶结颗粒材料内三相分布是处于平衡状态；继续迭代若干次，并且平均每个马尔科夫链的界面能，使其产生光滑的交界面。

Images