CN111681722A - 一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法 - Google Patents

Abstract

本发明公开了一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，包括如下步骤：构建两相多孔材料代表性体积单元；在孔相中随机生成一颗虚拟的树种子，进行单次随机树生长模拟；修剪不必要的死枝，并保留随机树贯通模型对立边界的主树干；对主树干进行平滑处理，记录下平滑处理后的主树干长度；统计主树干长度的均值；计算每个两相多孔材料代表性体积单元结构中孔相的几何曲折度；得到同样形状颗粒代表性结构单元的几何曲折度。本发明克服了任意复杂的孔隙网络曲折度难以定量提取的技术瓶颈，摆脱了以往理论方法只针对特定几何形貌的孔隙进行计算的理论约束，使得多孔材料内部的复杂孔隙网络曲折度的测定方法更具有普遍性和代表性。

Description

一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法

技术领域

本发明属于材料和土木工程技术领域，具体涉及一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法。

背景技术

多孔材料中因固相硬化颗粒的存在，导致活性介质在孔相中传输不能沿着直线路径到达目的地，于是便引入了曲折度的概念。曲折度早已成为很多实验、理论研究或数值模拟的一个热门主题。曲折度是对于多孔介质内部结构设计和优化所需要考虑的一种必不可少的结构因子。大多数多孔介质的微观结构和拓扑特征十分复杂多变，包括相互连通的不规则孔洞网络，并且孔洞的形状和大小也在较大范围内波动。此外，由于介质传输的路径不是严格的直线，而是树枝状的、弯曲的和相互连接的，直接导致很难通过常规的试验手段精确地测量多孔材料内部孔隙网络的曲折度。活性介质传输通常必须经过一段比它的起点到终点之间的最小距离长几倍的路径才能贯穿样本材料。现今存在诸多的曲折度模型，针对不同的研究对象主要分为四大类：水力曲折度、电力曲折度、几何曲折度以及扩散曲折度。其中几何曲折度是可以对多孔介质中复杂的微观结构进行详细描述的参数之一，这种复杂的结构特征很大程度上决定了材料的宏观性能包括传输和力学特性。几何曲折度主要解释了多孔材料的内部拓扑结构是如何影响介质传输的效率。因此，几何曲折度是联系多孔介质宏观性能与微观孔隙拓扑结构之间的桥梁。熟悉和了解几何曲折度对于人们设计和优化多孔材料起着相当关键的作用。

几何曲折度的表征主要分为数值方法和理论方法。数值模拟方法主要有LBM方法、随机游走方法以及图像处理法等。这类数值模拟方法计算规模宏大、效率不高。比如LBM方法需要不断求解大量的微分方程，并且对边界条件也非常敏感。对于随机游走方法，需要一次投放大量随机行走者同时进行计算，对计算机的硬件有很高的要求，其收敛性不容乐观。图像处理方法一般用于计算机数值模拟与现场试验相结合的方法，需要先在实验室对实体样本进行人工切片、图像采集，然后结合计算机图形学进行相关模拟计算，并且最后的数值计算结果非常依赖于所采集图像的分辨率。可以看出，现有数值模拟方法具有操作繁复、效率低下的弊端。另一方面，理论方法通常要假定的特定几何形貌的孔隙结构，例如：球形孔、柱状孔、椭球孔等，然后在此基础上采用相关的理论和复杂的数学推导给出特定几何形貌的孔隙曲折度模型，但是针对任意复杂的孔隙结构曲折度的理论计算目前尚未有此报道。

发明内容

发明目的：为了克服现有技术中存在的不足，提供一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其摆脱了以往理论方法只针对特定几何形貌的孔隙进行计算的理论约束，突破了数值模拟方法仅能针对松散堆积的简易形貌颗粒进行研究的短板，能够高效准确地测定两相多孔材料几何曲折度。

技术方案：为实现上述目的，本发明提供一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，包括如下步骤：

S1：构建两相多孔材料代表性体积单元，该代表性体积单元由椭球颗粒组成的固相和任意孔相组成，该代表性体积单元结构满足周期性边界条件；

S2：根据步骤S1构建的代表性体积单元结构，在孔相中随机生成一颗虚拟的树种子，树种子生长成若干树枝，每个生长的树枝不与任何椭球颗粒发生碰撞，形成随机树贯通模型，当随机树贯通对立边界时，随机树停止生长，则单个代表性体积单元的单次随机树生长模拟结束；

S3：根据步骤S2生成的随机树贯通模型，修剪不必要的死枝，并保留随机树贯通模型对立边界的主树干；

S4：对主树干进行平滑处理，记录下平滑处理后的主树干长度；

S5：重复进行N₁次步骤S2-S4，统计主树干长度的均值；

S6：根据主树干长度的均值，计算每个两相多孔材料代表性体积单元结构中孔相的几何曲折度；

S7：随机生成N₂个代表性体积单元结构，重复步骤S2到S6，将每次模拟得到的几何曲折度求均值，得到同样形状颗粒代表性结构单元的几何曲折度。

进一步的，所述步骤S2中随机树的具体生长过程为：

S21：在代表性体积单元结构的z＝L/2的截面上随机撒下种子，若该种子生成在颗粒内部，则该种子是死种子，舍弃；若种子生成在孔相中，则进入S22；

S22：在区间(0,1)内随机产生一个随机数q，若q小于等于prob则进入S23，否则进入S24；

S23：在空间中随机产生一个测试点random_p，然后在已经建立的随机树枝中找到距离random_p最近的随机树点，记为随机树点vertex_nearest，并将点random_p与点vertex_nearest连成线段s_temple，若线段s_temple不与空间中任何颗粒相交，则测试点random_p即为有效的随机树点，且线段s_temple记为随机树中新的树枝；若线段s_temple与空间中颗粒发生相交，则计算出所有交点，取所有交点中距离点vertex_nearest最近的一个交点为intersection_p，然后取二者的中点作为有效的随机树点，记为新的随机树点vertex_new，并且由点vertex_nearest与新的树枝点vertex_new所构成的线段作为随机树中新的树枝，进入S25；

S24：在颗粒堆积模型的对立边界面上随机产生S23中的测试点random_p，然后进行与S23相同的操作，完毕后进入S25；

S25：重复步骤S22-S24，直至随机树贯通对立边界面，则随机树停止生长，并保存所有树枝的数据。

进一步的，所述步骤S3中通过燃烧算法修剪不必要的死枝。

进一步的，所述步骤S4中平滑处理的具体过程为：通过不断迭代，每次迭代取相邻两根树枝的中点并相连，借助线段与椭球的碰撞检测算法判断该新取的线段是否与椭球颗粒碰撞，如果没有碰撞则保留为新树枝，否则不保留，当所有路径总长度变化数值小于设定值时停止迭代，得到贯通路径的平滑结构。

进一步的，所述步骤S6中通过公式(1)计算每个两相多孔材料代表性体积单元结构中孔相的几何曲折度：

τ_i＝L_e/L (1)

其中，L_e表示平滑处理后随机树主干长度的统计平均值，L是代表性体积单元立方型结构的尺寸，τ_i为第i个代表性体积单元几何曲折度预测值。

进一步的，所述步骤S7中同样形状颗粒代表性结构单元的几何曲折度的计算公式为：

其中，τ_i为第i个代表性体积单元几何曲折度预测值，N₂是随机紧密堆积代表性体积单元个数，τ为最终的曲折度统计平均值。

本发明首先构建一个两相的多孔材料细观结构代表性体积单元，该两相代表性体积单元结构由高堆积密度椭球颗粒组成的固相和余下的任意复杂孔相组成。其次，在代表性体积单元结构中随机生成单个虚拟的树种子，让种子在孔相中满足一定规律的生长成为若干的树枝，当生长的树枝贯通单元结构的对立边界(如：立方单元结构的上下边界、左右边界、前后边界)则停止生长。然后，通过“燃烧算法”修剪不必要的“死枝”，只保留连通对立边界的主树干。最后，计算平滑处理后的树干总长度。任意生成N₂个代表性体积单元结构，并且每个代表性体积单元经过N₁次的随机树生长模拟，取统计平均值来计算出每个代表性体积单元结构内孔隙网络的几何曲折度。

本发明方法克服了任意复杂的孔隙网络曲折度难以定量提取的技术瓶颈，摆脱了以往理论方法只针对特定几何形貌的孔隙进行计算的理论约束，突破了数值模拟仅能针对松散堆积的简易形貌颗粒进行研究且效率低下、模拟精度难以得到保证的技术壁垒，使得多孔材料内部的复杂孔隙网络曲折度的测定方法更具有普遍性和代表性。

有益效果：本发明与现有技术相比，具备如下优点：

1、本发明提出的几何曲折度测定方法相比于传统的数值模拟方法，如随机游走、LBM以及图像处理方法等，具有效率高、误差低等特点，极具有推广价值。

2、本发明提出的几何曲折度测定方法对孔隙结构的形貌没有任何限制，适用于任意复杂的孔隙结构，克服了以往的理论方法对孔隙形貌的制约，具有普适性。

3、本发明提出的几何曲折度测定方法优化了机器人路径搜寻的算法(rapidlyexploring random trees,RRTs)，避免了传统RRTs严重的随机性和盲目性，极大提高了几何曲折度测定的效率。

附图说明

图1为两相多孔材料几何曲折度示意图；

图2为由高堆积密度椭球颗粒和孔隙组成两相代表性体积单元的三维细观结构模型以及它的二维截面图；

图3为每个代表性体积单元结构内随机树生长模拟流程图；

图4为随机树生长示意图；

图5为死枝修剪对比示意图；

图6为利用本发明计算得到的不同长径比椭球颗粒堆积体系中的孔隙结构曲折度。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明。

本实施例中将本发明方法应用于高堆积密度椭球颗粒多孔材料的几何曲折度测定，其提供了一种几何曲折度计算方法，包括如下步骤：

步骤1：生成两相多孔材料细观结构的立方型代表性体积单元，该两相代表性体积单元由高堆积密度的椭球颗粒组成的固相和余下的任意复杂孔相组成，代表性体积单元结构满足周期性边界条件。

周期性边界条件是相对于材料结构代表性体积单元而说，一般取一个代表性体积单元，它是能代表整体材料性质的，该性质在单元边界上具有周期可复制性，因此，在材料结构的数值建模中，通常称该代表性体积单元满足周期性边界条件，它的对立面即是刚性边界条件。

本步骤中高堆积密度的颗粒增强复合材料两相细观结构模型的构造方法为：首先在较小的正方体容器中随机地投入一定数量的椭球颗粒，赋予每个颗粒随机的中心坐标以及随机的旋转方向，初始时会有大量的重叠；再以松弛迭代的方法根据重叠量对每个颗粒逐一进行平移和转动，逐渐地消除颗粒相互间的重叠，直到重叠量小于一定值停止，最终形成的高堆积密度的两相复合材料代表性体积单元如图2所示。

步骤2：基于步骤1构建一个代表性体积单元结构，使用“改进随机树算法”进行几何曲折度的测定，其步骤分为以下三个分步：

分步1：由于立方型代表性体积单元有三组对立的边界面，对于曲折度模拟是等价的。因此，在每次随机树模拟前随机选定用于模拟的对立边界，如：上下边界、左右边界抑或前后边界。然后，在代表性体积单元结构的z＝L/2(或y＝L/2抑或x＝L/2)的截面上随机撒下一颗种子seed，使其落在该截面的孔隙中，参照图1，z为结构单元的纵坐标，L为代表性体积单元结构的尺寸，即立方体的边长。以基于椭球颗粒一般方程定义的函数F_A(r-r_A,R_i(Ω_A))作为判别依据，判断随机点是否落入颗粒的内部：

当任意点

时分别对应：

其中椭球A中心点(x_A,y_A,z_A)包含于向量r-r_A，其中r_A表示椭球中心点与坐标系的原点连成的向量，r是空间任意点与坐标原点组成的向量。椭球的三个半轴长a,b,c包含于向量R_i。其中Ω_A包含了椭球A的三个欧拉角α,β,γ。

r-r_A可展开为式(4)：

r-r_A＝(x-x_A,y-y_A,z-z_A)^T (4)

F_A(r-r_A,R_i(Ω_A))展开为：

(r-r_A)^TA^-1(r-r_A)(5)

A矩阵可以展开为：

其中R_i(Ω_A)包含椭球的三个空间欧拉角，分别表示椭球的三个半轴向量R₁(Ω_A)，R₂(Ω_A)，R₃(Ω_A)，每个半轴的向量分别与它们的转置矩阵相乘，求和所有的半轴向量，得到A矩阵。S(Ω_A)是空间旋转矩阵，其展开式为：

通过一系列的整理，得到空间任意椭球的一般表达式：

其中，A_ij ^*代表矩阵A的伴随矩阵对应位置的值，|A|代表矩阵A的行列式值；而(x,y,z)代表任意随机点的坐标，(x_A,y_A,z_A)为椭球的中心坐标。通过将随机点的坐标(x,y,z)代入公式(8)，并结合公式(3)可以准确地得到该随机点与空间任意椭球的位置关系。这样可以确保随机种下的种子seed落于模型的孔隙中。

分步2：开始从种子seed出发，通过不断循环生长出若干树枝，具体的算法流程如图3所示，此过程的具体细节如下：

(a)首先在区间(0,1)内预设一个变量prob，本实施例中设定为0.4(这样大概率决定了随机树朝着对立边界方向生长)。然后，在区间(0,1)内随机产生一个随机数q，如果q大于prob则随机树以各一半的几率随机选择对立边界面(以上下边界面为例)，具体操作是在区间(0,1)内随机产生随机数q₁，若q₁大于等于0.5则在上边界面上随机产生一个测试点random_p，否则在下边界面上随机产生一个测试点random_p；反之，如果q小于等于prob则在空间中随机产生一个测试点random_p。然后在已经建立的随机树枝中找到一个随机树点距离random_p最近，记为随机树点vertex_nearest(如果是第一次生长，此点即为种子seed)，并将点random_p与点vertex_nearest连成线段s_temple。若线段s_temple不与空间中任何颗粒相交，则测试点random_p即为有效的随机树点，且线段s_temple记为随机树中新的树枝(如图4中线段1-seed，2-seed和3-seed所示)，进入步骤(d)；若线段s_temple与空间中颗粒发生相交，则需要求解线段s_temple与颗粒相交的交点，具体步骤如下所述。

(b)凭借“线段与椭球碰撞检测算法”判断线段s_temple是否与模型中的颗粒发生碰撞。对于该算法的设计，具体细节是通过设定一个精度值Accuracy，其大小取椭球三个半长轴最大值的1/10，即：

Accuracy＝max(a,b,c)/10 (9)

其中，a,b,c分别是椭球颗粒的三个半轴的长度。

根据此精度值将指定的空间线段划分为若干份，假定其中分成的若干小线段和节点，并且节点的数量N_node1满足以下公式：

其中，l_segment为线段s_temple的长度，函数ceil()为求不小于括号内值的最小正整数值。

通过循环利用分步1中的函数F_A(r-r_A,R_i(Ω_A))逐一判断每个节点与椭球的位置关系，若所有的节点都在椭球外面则该线段与椭球不发生碰撞。反之，若有大于等于1个节点在椭球内部则该线段与椭球发生碰撞。如果发生了碰撞则进入步骤(c)；

(c)凭借“线段与椭球交点计算算法”计算出线段s_temple与模型中的颗粒的所有交点，对于该算法的设计，具体细节是通过设定一个更小的精度值Accuracy^满足公式(11)来近似求得具体交点：

Accuracy^＝max(a,b,c)/20 (11)

其中，a,b,c分别是椭球颗粒的三个半轴的长度。

根据此精度值将指定的空间线段划分为若干个微线段和节点，并且节点的数量N_node2满足以下公式：

其中，l_segment为线段s_temple的长度，函数ceil()为求不小于括号内值的最小正整数值。

首先判断线段s_temple的两个端点vertex_nearest和random_p与椭球的位置关系，并由此关系求出线段s_temple与每颗发生碰撞的颗粒的具体交点，主要由以下两种情况：

vertex_nearest与random_p均在椭球外部，这说明该线段与椭球有两个交点，则分别从两个端点出发，逐一判断各节点与椭球的位置关系，一旦找到第一个在椭球内部的点，则该点即为线段与椭球的交点；

vertex_nearest在外部而random_p在椭球内部，则仅从vertex_nearest出发逐一判断各节点与椭球的位置关系，一旦找到第一个在椭球内部的点，则该点即为线段与椭球的交点；

这样，找出上述所求的所有交点中距离点vertex_nearest最近的一个交点intersection_p，然后取二者的中点作为有效的随机树点，记为新的随机树点vertex_new，并且由点vertex_nearest与新的树枝点vertex_new所构成的线段作为随机树中新的树枝，如图4所示，若random_p为点5的情况，线段3-5与颗粒发生碰撞，则通过“线段与椭球交点计算算法”计算出所有交点，并取距离点3最近的交点F，然后取点3和点F的中点6作为新的随机树点vertex_new，且线段3-6作为随机树中新的树枝。进入步骤(d)。

(d)判断新添加的树枝点vertex_new是否到达对立的边界上，具体细节：当随机树枝节点vertex_new＝(v_x,v_y,v_z)满足v_z∈[0,Z_b]则视为连通样本的下边界。同样地，当v_z∈[Z_max-Z_b,Z_max]则视为连通样本的上边界(另外两组对立边界是否连通依此类推)。当连通对立边界的树枝数量均大于等于1(至少有1条)的时候，随机树成功贯通模型上下边界，结束搜索进程，进入步骤(e)，否则回到步骤(a)继续生长和探索；

其中，(v_x,v_y,v_z)为新添加的节点vertex_new的空间坐标；Z_b为一个预设值，本实施例中取Z_b＝0.05L；Z_max为代表性体积单元在z方向最大坐标值，此处等于L，坐标系见图1。

(e)利用“死枝修剪算法”提取随机树的主干，其原理是通过不断迭代，用“燃烧算法”从外围逐步删除“孤立”树枝。对于“孤立”树枝的定义：只有一端与其他树枝相连的线段。在修剪前，先对进入对立边界的树枝进行特殊标记，避免在修剪过程中被删除。每次迭代记录下本次迭代中检测到孤立树枝的数量N_isolated，当N_isolated＝0时停止修剪。如图5所示为修剪前后的随机树对比图，图中灰色线段即为被修剪的“死枝”而剩下的平滑、连续的黑色枝丫为随机树的主干，即主要渗流通道。进入步骤(f)。

(f)通过将修剪过的随机树利用“E-P”算法进行平滑处理，细节为通过不断迭代，每次迭代取相邻两根树枝的中点并相连，借助线段与椭球的“碰撞检测算法”判断该新取的线段是否与椭球颗粒碰撞，如果没有碰撞则保留为新树枝，否则不保留。当所有路径总长度变化数值小于1.0×10^-2时停止迭代，这样得到贯通路径的平滑网络结构，计算贯通单元体有效随机树主干的总长度，记为l_k。k为每个代表性体积单元成功完成一次随机树生长模拟的计数器，当k达到每个单元体预设的随机树生长模拟总次数N₁时，则进入分步3，否则返回步骤(a)。

分步3：通过公式(13)计算每个代表性体积单元进行N₁次随机树生长模拟后的孔隙几何曲折度：

其中，τ_i为第i个代表性体积单元几何曲折度预测值。

步骤3：生成N2个同样形状颗粒的代表性体积单元，并重复步骤2，将每次模拟得到的几何曲折度求均值，得到同样形状颗粒代表性结构单元的几何曲折度τ：

其中，τ_i为第i个代表性体积单元几何曲折度预测值，N₂是随机紧密堆积代表性体积单元个数，而τ为最终的曲折度预测值。

对于模拟时间，本发明方法单次随机树生长模拟直到贯通的时间均值仅需30秒。另外，我们利用本发明专利得到不同长径比的椭球紧密堆积结构体系中孔隙曲折度数值结果，如图6所示。

Claims (6)

1.一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：包括如下步骤：

S1：构建两相多孔材料代表性体积单元，该代表性体积单元由椭球颗粒组成的固相和任意孔相组成，该代表性体积单元结构满足周期性边界条件；

S2：根据步骤S1构建的代表性体积单元结构，在孔相中随机生成一颗虚拟的树种子，树种子生长成若干树枝，每个生长的树枝不与任何椭球颗粒发生碰撞，形成随机树贯通模型，当随机树贯通对立边界时，随机树停止生长，则单个代表性体积单元的单次随机树生长模拟结束；

S3：根据步骤S2生成的随机树贯通模型，修剪不必要的死枝，并保留随机树贯通模型对立边界的主树干；

S4：对主树干进行平滑处理，记录下平滑处理后的主树干长度；

S5：重复进行N₁次步骤S2-S4，统计主树干长度的均值；

S6：根据主树干长度的均值，计算每个两相多孔材料代表性体积单元结构中孔相的几何曲折度；

S7：随机生成N₂个代表性体积单元结构，重复步骤S2到S6，将每次模拟得到的几何曲折度求均值，得到同样形状颗粒代表性结构单元的几何曲折度。

2.根据权利要求1所述的一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：所述步骤S2中随机树的具体生长过程为：

S21：在代表性体积单元结构的z＝L/2的截面上随机撒下种子，z为结构单元的纵坐标，L为代表性体积单元结构的尺寸，即立方体的边长，若该种子生成在颗粒内部，则该种子是死种子，舍弃；若种子生成在孔相中，则进入S22；

S22：在区间(0,1)内随机产生一个随机数q，若q小于等于变量prob则进入S23，否则进入S24；

S23：在空间中随机产生一个测试点random_p，然后在已经建立的随机树枝中找到距离random_p最近的随机树点，记为随机树点vertex_nearest，并将点random_p与点vertex_nearest连成线段s_temple，若线段s_temple不与空间中任何颗粒相交，则测试点random_p即为有效的随机树点，且线段s_temple记为随机树中新的树枝；若线段s_temple与空间中颗粒发生相交，则计算出所有交点，取所有交点中距离点vertex_nearest最近的一个交点为intersection_p，然后取二者的中点作为有效的随机树点，记为新的随机树点vertex_new，并且由点vertex_nearest与新的树枝点vertex_new所构成的线段作为随机树中新的树枝，进入S25；

S24：在颗粒堆积模型的对立边界面上随机产生S23中的测试点random_p，然后进行与S23相同的操作，完毕后进入S25；

S25：重复步骤S22-S24，直至随机树贯通对立边界面，则随机树停止生长，并保存所有树枝的数据。

3.根据权利要求1所述的一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：所述步骤S3中通过燃烧算法修剪不必要的死枝。

4.根据权利要求1所述的一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：所述步骤S4中平滑处理的具体过程为：通过不断迭代，每次迭代取相邻两根树枝的中点并相连，借助线段与椭球的碰撞检测算法判断该新取的线段是否与椭球颗粒碰撞，如果没有碰撞则保留为新树枝，否则不保留，当所有路径总长度变化数值小于设定值时停止迭代，得到贯通路径的平滑结构。

5.根据权利要求1所述的一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：所述步骤S6中通过公式(1)计算每个两相多孔材料代表性体积单元结构中孔相的几何曲折度：

τ_i＝L_e/L (1)

其中，L_e表示平滑处理后随机树主干长度的统计平均值，L是代表性体积单元立方型结构的尺寸，τ_i为第i个代表性体积单元几何曲折度预测值。

6.根据权利要求5所述的一种含有高堆积密度椭球颗粒的多孔材料曲折度测定方法，其特征在于：所述步骤S7中同样形状颗粒代表性结构单元的几何曲折度的计算公式为：

其中，τ_i为第i个代表性体积单元几何曲折度预测值，N₂是随机紧密堆积代表性体积单元个数，τ为最终的曲折度统计平均值。

Images