仓储粮堆三维非规则孔道网络模型的优化研究

陈鹏枭，朱文学，吴建章，杨 帆，杨德勇，黄 凯

1.河南工业大学 粮油食品学院，河南 郑州 450001

2.中国农业大学 工学院， 北京 100083

3.南京农业大学 工学院， 江苏 南京 210095

仓储粮堆干燥过程涉及粮粒内湿热传递、孔隙之间的气体流动以及孔隙间气体与粮粒之间的湿热传递，其研究涉及生物、化学、热力学以及流体力学等多种学科领域，加之干燥过程中粮堆内各物理参数较难测定，直接对湿热传递机制进行研究十分困难。随着计算机技术的发展，数值模拟方法成为研究粮食干燥过程湿热传递机制的重要手段[1]。

基于不同的模型假设，目前关于仓储粮堆湿热传递的模拟模型可分为两种。

一种是基于连续介质假设的仓储粮堆湿热传递模型[2-5]。Hammami等[6]对圆筒仓粮堆进行了湿热传递模拟，研究发现通风温度对谷物温度变化具有重要作用，并总结了气流特性对仓内粮食温湿度的影响。Wang等[7]建立了描述深层粮堆内部湿热传递的数学模型，研究发现在冷却干燥通风过程中，降温和降水同时存在，冷却前沿的移动速度大大快于干燥前沿，且粮堆中温度和水分分布具有分层的现象。Chen等[8-9]以小麦为研究对象，采用局部非平衡法，建立了仓储小麦三维通风干燥湿热传递模型，通过模拟研究发现固定床干燥有明显的干燥前沿，干燥不均匀度较大，且干燥前沿附近的温度梯度最大。

另一种是基于离散介质假设的孔道网络湿热传递模型[10-13]，其物理模型如图1所示。Kharaghani等[14-16]通过X射线获取球形堆积体三维图片，提取该堆积结构的孔道网络信息，并基于该孔道网络信息模拟了物体液相传递以及干燥过程，研究发现，在第一个干燥阶段，模拟所得干燥曲线与实验结果基本吻合，之后实验曲线发生明显干燥梯度，而模拟结果则表现为无梯度。Yuan等[17-18]以玉米为研究对象，构建了单层籽粒堆积孔道网络物理模型，并结合多尺度湿热传递数学模型对仓储粮堆干燥过程进行了模拟，研究发现物料孔隙率、配位数越大，物料干燥越快，升温也越快，达到相同湿含量所需干燥时间越短。黄凯[19]提出了构建三维非规则孔道网络物理模型的方法，并利用三维规则孔道网络模拟了粮堆通风干燥过程，研究发现任意时刻粮堆孔隙内的气相分压随粮层深度的降低而降低，随着通风的进行，玉米颗粒平均含水率不断降低，温度先迅速降低后缓慢升高，内部不同层含水率由外而内逐步降低。Chen等[20-21]构建了仓储粮堆二维非规则以及三维非规则孔道网络模型，并结合孔道网络湿热传递数学模型，对仓储粮堆通风干燥过程进行了模拟，经干燥实验验证，两种孔道网络模型皆能够描述粮堆机械通风干燥过程，且研究结果表明孔道网络模型中的颗粒配孔数与拓扑结构对模型模拟结果影响较大。

从前述分析可以得出，基于离散介质假设的非规则孔道网络物理模型能够清晰描述仓储粮堆内部各颗粒之间、孔喉之间以及颗粒与孔喉之间的空间位置关系，这在探究粮食湿热传递机理与三传过程模拟领域具有较大优势。作者拟以仓储粮堆为研究对象，利用三维CT扫描设备与Avizo软件对颗粒堆积样件的孔喉进行提取，通过EDEM软件模拟粮食颗粒堆积状态，基于泰森多边形原理，通过voro++软件处理堆积体中的颗粒坐标及半径等数据，制定孔隙划分规则，对不合理孔隙进行处理，结合Avizo软件获取的孔喉信息，获得能够描述仓储粮堆内部空间结构的三维非规则孔道网络物理模型。

图1 粮堆孔道网络模型构建过程图

1 实际颗粒堆积结构孔道网络的提取

1.1 试验材料与方法

试验材料：颗粒半径为4 mm的玻璃珠若干，内边长为3 cm的正方体容器。

试验设备：英华检测(上海)有限公司的phoenix v/ tome/ x m全能型X射线微焦点CT系统。

试验方法: 将备好的玻璃珠装满内边长为3 cm的容器，所有颗粒完全在容器内部，对样品进行扫描，采集原始图像，得到反映样件内部结构的图片和数据，使用Avizo软件对所得样件内部孔喉结构参数进行分析。

1.2 试验结果与分析

图2为玻璃球形颗粒堆积结构的三维CT扫描图片。从图2可以看出，该CT从3个不同的维度对物料进行步进扫描，得到关于物料结构的三维立体图像。

经过无损扫描以及Avizo软件对样件孔道网络的孔喉提取发现，该试验样品中共含有51个颗粒、66个孔隙、205个喉道，其平均配位数为6.246，孔隙率为49.36%，孔隙密度为2.4×106个/m3。图3为球形颗粒堆积结构内部喉道长度分布，从图3可得，该样品的喉道长度主要分布在6.0～14.0 mm之间，大约占总喉道数的96%，喉道长度分布大致为Poisson分布。图4为球形颗粒堆积结构的孔隙半径分布，从图4可得，该结构孔隙半径多在2.0～10.0 mm之间。

图2 球形颗粒堆积结构局部三维扫描图

图3 球形颗粒堆积结构喉道长度分布

图4 球形颗粒堆积结构孔隙半径分布

2 仓储粮堆三维非规则孔道网络物理模型的构建

2.1 基于泰森多边形的孔道网络物理模型的构建

利用EDEM软件模拟颗粒堆积状态，假设容器形状为方仓，铝合金材质，玉米颗粒为球体。颗粒、粮仓仓壁以及颗粒与颗粒之间的物性参数，如表1所示。通过EDEM软件模拟得到堆积体中玉米颗粒的信息，图5所示为经过voro++软件处理过的颗粒堆积结构对应的孔道网络划分。通过上述方法构建的物理模型大小为-0.020 12≤x≤0.020 12、-0.020 12≤y≤0.020 12、0≤z≤ 0.063 36 (单位: m)，总共有179个颗粒，孔隙数为991，喉道数为1 978，颗粒半径为4.0 mm，孔隙配位数主要集中在4。经计算，该模型的孔隙率为52.39%，孔喉连通熵为6.775 4，孔隙密度为1.03×107个/m3。与样件提取所得的参数对比，从规模数上讲，提取所得的孔道网络物理模型的颗粒数、孔隙数与喉道数之比为1∶1.275∶3.980，而构建所得物理模型的三参数之比为1∶5.536∶11.050，相关参数比值相差较大；提取所得的三维非规则孔道网络物理模型的平均配位数为6.246，构建所得三维非规则孔道网络物理模型平均配位数为3.992，相差较大；此处构建获得的孔隙密度为提取孔隙密度的4.3倍左右，相差较大。

表1 模拟所用材料物性参数

图5 颗粒堆积结构孔道网络划分

图6为构建的三维非规则孔道网络模型中孔间距的长度分布规律，从图6可得，孔间距中有较大一部分长度小于5 mm，孔间距等于喉道长度加上喉道两端的两孔隙半径，依此定义可知，小于5 mm的喉道个数会更多，这与提取结果相差较大。

图6 三维非规则孔道网络模型孔间距长度分布

在等径球状颗粒最紧密堆积的多孔介质的结构中，等径球状颗粒(半径为R)多孔介质模型最小单元体由 4 个球组成，呈正四面体结构。其最小单元体的孔隙半径(r)，如式(1)所示。

(1)

如有喉道长度小于2r，则该喉道不符合实际情况。经统计，构建的三维非规则孔道网络模型中孔间距小于2r的有191个。因此,需对该三维非规则孔道网络物理模型进行改进。

2.2 仓储粮堆多尺度物理模型的改进

将孔道网络中不符合实际情况的孔隙按照一定规则进行合并，使其更符合真实孔道网络结构。

2.2.1 孔隙合并规则

(1)合并规则一。在孔道网络模型中，孔心距小于2r的两孔隙需要合并，形成的新孔隙的孔心坐标为此两孔中心连线的中点。

(2)合并规则二。如图7所示，图7a为合并孔隙前的孔喉结构示意图，孔1的相邻颗粒为C1和C3，孔2的相邻颗粒为C1和C2，孔3的相邻颗粒为C1、C2、C3和C4，孔1、2、3两两通过不同的喉道相连接，此时可以认为孔1、2、3都在颗粒C1、C2、C3、C44个颗粒围成的大孔隙内部，可以将其合并成一个大孔隙，为了简化算法，直接将孔1和孔2合并到孔3，将该孔隙拓扑结构改变成图7b，即完成该类孔隙的合并。

2.2.2 孔隙半径提取算法

平均半径算法是指将孔隙周围颗粒中心与孔心的距离减去对应颗粒半径的值的平均值视为孔隙半径的一种算法。以一个配位数为3的孔隙为例，孔隙周围有3个颗粒，3个颗粒的半径分别为r1、r2、r3，从孔隙中心到3个颗粒中心的长度分别为l1、l2、l3，则平均距离法得到孔隙半径rp。

(2)

以平均半径算法获得各个孔隙的半径之后，合并孔间距小于两孔半径之和的孔隙，新得孔隙的中心为两合并孔隙中心连线的中点，再应用平均半径算法，获得新的孔隙半径。

2.3 改进后孔道网络模型结构参数分析

算法改进后，获得的三维非规则孔道网络物理模型颗粒数为179，孔隙数为237，喉道数为702，颗粒半径为4.0 mm，平均配位数为5.924。将上述参数与改进前的模型参数对比发现，提取所得的孔道网络物理模型的颗粒数、孔隙数与喉道数之比为1∶1.275∶3.980，而改进后物理模型的3个参数之比为1∶1.324∶3.922，单从模型规模数上讲，两物理模型相差较小，可以互相表征。图8为改进后三维非规则孔道网络模型配位数分布，从图8可得，配位数主要集中在3、4、5和6，大致呈Poisson分布。图9为改进后三维非规则孔道网络模型喉道长度分布，可知该孔道网络模型的喉道长度主要分布在3.0～13.0 mm之间，大致呈Poisson分布。

图8 改进后三维非规则孔道网络模型配位数分布

图9 改进后三维非规则孔道网络模型喉道长度分布

图10为改进后三维非规则孔道网络模型孔隙半径分布，从图10可知，在该模型中孔隙分布大致在2.0～10.0 mm之间，且存在直径大于10 mm的孔隙，这是由于在孔道网络模型的边界处，孔隙只有一侧有颗粒，限制孔隙半径的条件比仓储粮堆内部的孔隙少，所以此处孔隙半径较大；根据提取信息可得该孔道网络模型的孔喉连通熵为5.712 4。

图11、图12为改进后三维非规则孔道网络模型颗粒配孔数、配喉数分布，颗粒配孔数在7～18之间，颗粒配喉数在4～20之间，其中大部分配孔数集中在8～16之间。

图10 改进后三维非规则孔道网络模型孔隙半径分布

图11 改进后三维非规则孔道网络模型颗粒配孔数分布

图12 改进后三维非规则孔道网络模型颗粒配喉数分布

3 实际结构与构建模型孔喉数据对比分析与验证

3.1 实际结构与构建模型孔喉数据对比

将样件提取数据与三维非规则孔道网络模型构建数据进行对比分析，如表2所示。从表2可以看出，构建的三维非规则孔道网络模型结构参数与提取的结果相符合，可以说明该模型在统计意义上与实际颗粒堆积体的各项孔隙参数都较为吻合，且构建的孔道网络模型结构参数中有颗粒配孔数与颗粒配喉数，能够反映堆积结构中颗粒与孔、喉之间的拓扑关系，说明构建的三维非规则孔道网络模型要比Avizo软件提取所得的三维非规则孔道网络模型更能反映颗粒堆积结构内部的孔隙结构真实情况。

表 2 三维提取信息与三维构建信息对比

3.2 构建所得孔道网络物理模型流通性验证

图13 三维非规则孔道网络模型初始风压与出口流速关系

在粮食通风干燥过程中，空气由进风口进入，穿过粮食区域，从出风口排出，这个过程可以看作是一个流体穿过非固结多孔介质内部的过程，采用达西定律来描述这个过程[20-21]，并利用程序计算初始风压与出口流速之间的关系。图13为三维非规则孔道网络模型初始风压与出口流速关系，从图13可以看出，底部风压与出口风速大致呈线性变化。该结果可以证明两个结论：第一，证明了该物理模型能够有效表征非固结多孔介质的孔隙连通性；第二，证明了该物理模型可以结合湿热传递数学模型用来描述非固结多孔介质内部流体流速分布情况。

4 结论与展望

采用孔道网络方法构建了能够描述仓储粮堆内部结构的孔道网络模型，并运用达西定律验证了物理模型内部孔喉结构的流通性，通过分析研究，可得结论如下。

利用三维CT扫描设备与Avizo软件提取所得颗粒堆积结构中喉道长度主要分布在6.0～14.0 mm之间，喉道长度分布大致为Poisson分布；孔隙半径多在2.0～10.0 mm之间。

通过EDEM软件模拟粮食颗粒堆积状态，基于泰森多边形原理，通过voro++软件处理堆积体中的颗粒坐标及半径等数据，结合孔隙划分规则对不合理孔隙进行处理，所得的颗粒堆积结构参数中配位数主要集中在3、4、5和6，喉道长度主要分布在3.0～13.0 mm之间，孔隙半径分布大致在2.0～10.0 mm之间，3个参数大致均呈Poisson分布；颗粒配孔数在7～18之间，颗粒配喉数在4～20之间，其中大部分配孔数集中在8～16之间。

通过构建所得三维非规则孔道网络模型较Avizo软件所获得的孔隙结构参数更全面，更能准确地描述仓储粮堆内部颗粒、孔隙以及喉道之间的关系。