平房仓中小麦堆的孔隙率分布值研究

邵小龙，程绪铎

（南京财经大学 食品科学与工程学院/江苏省现代粮食流通与安全协同创新中心，江苏 南京 210046）

我国自1998年以来新建的储备粮库及地方粮库的仓型一般为浅圆仓、立筒仓和高大平房仓，除港口和交通枢纽中转库以浅圆仓和立筒仓为主外，其余多为高大平房仓，其储存量约占总储存量的70%[1]。

小麦储藏在平房仓中，受到自重、内摩擦力以及仓壁的支持力，小麦堆各处产生压应力与切应力，从而产生体积应变与形状应变，体积缩小，孔隙率减小。平房仓中的小麦堆的应力与应变分布是非均匀的，因此孔隙率的分布也是非均匀的。孔隙率的大小与分布决定了粮食通风的阻力与均匀性，还影响粮堆中的热量与水分的迁移速率[2-4]。

粮食孔隙率的测定始于20世纪60年代，Thompson等[5]使用气体比重瓶测定了燕麦、小麦、高粱、大豆、黑麦、大麦、玉米和紫花苜蓿种子的孔隙率。之后，Gustafson等[6]和Chang[7]也利用气体比重瓶测定了玉米、小麦和高粱的孔隙率。郝倩等[8]开发了一种粮食孔隙度就仓测量系统。还有学者提出了光学法[9-10]、声学法[11]等。但光线与声波无法到达仓内粮堆深处，不适合用来测定粮仓内粮堆孔隙率。

到目前为止，还未见平房仓中粮堆的孔隙率分布的理论与实验成果报道。本文选定修正剑桥模型为小麦堆的应力与应变关系本构模型，使用有限元方法求解平房仓中小麦堆修正剑桥模型，计算出平房仓内小麦堆的应变分布值，由应变分布值与表层孔隙率（无压缩）计算出平房仓中小麦堆的孔隙率分布值。

1 材料与方法

1.1 实验材料

实验小麦品种名为淮麦44号，产地江苏南京，手工剔除破碎与不成熟颗粒，初始湿基含水率为12.03%，小麦的最大粒径、中等粒径和最小粒径分别为5.88、3.07、2.79 mm（随机取样100颗小麦粒，用游标卡尺测得其3个粒径，精确到0.01 mm，结果取其平均值）。实验前，取20 kg小麦样品，放置在太阳光下晾晒72 h，测得其湿基含水率为9.72%；再取20 kg小麦样品，分成两份，各自加适当的蒸馏水，然后将其装入塑料袋中封闭好，置于恒温箱中，温度设定5 ℃，1周后从恒温箱中取出，测得两份小麦样品湿基含水率分别为14.61%和16.55%。用标准烘箱干燥法将10 g样品在130 ℃下干燥19 h，重复3次，测定小麦的含水率。实验用的小麦的湿基含水率分别为：9.72%、12.03%、14.61%、16.55%。

1.2 实验仪器

TSZ-6A型应变控制式三轴仪：南京土壤仪器厂有限公司；LKY-1型粮食孔隙率测定仪：南京土壤仪器厂有限公司。

1.3 小麦堆的应力与应变关系本构模型

平房仓内的小麦堆在重力、仓壁抵抗力的作用下，小麦堆内部各处产生了应力，由于应力的作用，小麦堆各处又产生了应变，小麦堆的应力与应变相互关联，它们之间关系的微分方程就称为应力与应变关系的本构模型。要准确地求解平房仓内小麦堆的应变分布必须选定恰当的应力与应变关系的本构模型。小麦籽粒松软，小麦堆孔隙大，在应力的作用下产生较大的体积应变与形状应变。修正剑桥模型中以体积应变和剪切应变作为应变分量，能表征小麦堆大体积压缩的特性。所以本研究选定修正剑桥模型为小麦堆的应力与应变关系的本构模型。

修正剑桥模型的弹塑性矩阵的增量形式如下[12]：

式中：q为广义剪切应力，Pa；p为平均主应力，Pa；εv为体积应变；εs为剪切应变；K为体积模量，kPa ；G为剪切模量kPa。

由式（1）～式（5）可知，修正剑桥模型中需要确定的参数为临界状态应力比M、等向膨胀指数κ、对数硬化模量λ、弹性模量E、泊松比υ、初始孔隙比e0。

1.4 修正剑桥模型参数的确定

1.4.1 初始孔隙比e0的测定

小麦堆孔隙率ε为小麦堆中的孔隙体积与整个小麦堆体积之比。孔隙比e指小麦堆中孔隙体积与小麦堆籽粒体积之比。本实验通过LKY-1型粮食孔隙率测定仪（见图1）测定小麦堆的初始孔隙率，初始孔隙比可由初始孔隙率计算得到。

图1 孔隙率测定仪示意图

如图1所示，容器A和B的容积相等，将小麦样品装入容器B中，装满后将容器B密封。将阀门2关闭，然后将阀门1和3打开。将一定压力的气体通过空气压缩机充满容器A。待压力表达到一定数值后，将阀门1关闭，等待压力稳定后，记录下此时压力表读数（P1）；关闭阀门3，然后打开阀门2，待容器A和B的压力达到平衡后，记录压力表读数（P2）。设空气是理想气体，根据理想气体等温特性原理可以推导出小麦样品的初始孔隙率（一个大气压下的孔隙率）为：

式中：Vs为样品籽粒的体积，m3；V0为小麦样品初始体积，m3。

1.4.2 临界状态应力比M的测定

临界状态应力比M可通过轴向压缩实验测定[12]。使用应变控制式三轴仪对小麦样品进行不同围压（30、50、70、90、110 kPa）的剪切压缩实验，每当样品位移增加0.4 mm，则记录下样品体积变化量和测力计读数；当测力计读数产生最大值时，即为最大主应力差，记录下对应的p和q值。将记录的p和q值拟合成临界状态应力直线，M的值即为直线的斜率，进行3次平行实验。

1.4.3 对数硬化模量λ和等向膨胀指数κ的测定

通过应变控制式三轴仪对小麦样品进行各向等压压缩，将围压σ3由0 kPa增加到200 kPa，每当σ3增加了5 kPa，则记录一次样品体积减少量。然后，再由200 kPa卸载到0 kPa，当σ3减小5 kPa，记录一次体积变化量。最后，将p（围压）和孔隙比e绘制成加载曲线和卸载曲线，其直线斜率分别为λ和κ的值[12]，进行3次平行实验。

1.4.4 弹性模量E和泊松比υ的测定

（1）弹性模量E的测定：参照SL237—1999[13]。选择实验测定的最大主应力差的10%作为施压增量，对粮堆进行每级递增加压和递减卸压。加压过程中，在对粮堆施加第一级压力后，进行计时，对粮堆每隔1 min增加一级压力，并且记录下相应的压缩位移，直至增加到第4级压力；然后由第4级压力逐级卸压，同样每隔1 min减去一级压力，并记录粮堆压缩位移，直至施加的轴向压力为0；反复加卸载4次之后，对粮堆进行持续压缩直至粮堆被破坏，进行3次平行实验。

由最后一次加卸载循环中得到的实验数据计算出弹性模量。

式中：E为粮堆的弹性模量，kPa；ΔP为轴向荷载，kPa；Δhe为粮堆的弹性变形量，mm；hc为试样固结后的高度，mm。

（2）泊松比υ的测定：粮堆进行等向压缩实验，由式（11）计算出体变弹性模量[14]。

式中：B为小麦堆的体变弹性模量，kPa；ΔV为小麦堆的弹性体积增量，m3；V为小麦堆装样后的体积，m3。

则泊松比υ为：

1.5 平房仓中小麦堆应力与应变的有限元计算方法

平房仓内储藏小麦后，小麦堆各处都产生了应力与应变，小麦堆各处的孔隙率与初始孔隙率及各处的应变密切相关。只要测定出小麦堆的初始孔隙率并计算出小麦堆各处的应变，就能计算出小麦堆各处的孔隙率。小麦堆各处的应变可使用有限元方法求解小麦堆的修正剑桥模型而得出。

1.5.1 平房仓仓壁的几何与力学参数的选定

混凝土仓壁可视为刚性体。平房仓尺寸参考南京铁心桥国家粮食储备库粮仓的实际大小，具体仓壁材料的几何与力学参数见表1。

表1 平房仓仓壁的几何与力学参数

1.5.2 平房仓中小麦堆有限元模型构建与运行步骤

（1）在 Part 步骤中采用三维可形变模型（3D Deformable）建模，创建平房仓和小麦堆部件，尺寸见表1。

（2）在Property步骤中分别设定仓壁的材料属性和粮堆的塑弹性属性（修正剑桥模型选择塑性模型Clay Plasticity和弹性模型Porous Elastic）。

（3）在Assembly 步骤中将平房仓仓壁和粮堆部件组合起来。

（4）在 Step步骤中设定分析步。

（5）Interaction步骤是定义各个部件之间的相互接触关系。粮堆与仓壁之间存在相互作用，在ABAQUS中设定一个完整的接触，包括两个部分：一是接触对，选择点对面离散方法，选择仓壁内侧面为主接触面，粮堆外侧面为从属接触面；二是定义相互作用力，两个接触面之间会产生摩擦力。其中，库仑模型是较为常用的摩擦模型，即通过摩擦系数来描述两表面间的摩擦关系。

（6）在 Load 步骤进行载荷以及边界条件的设定，约束平房仓仓底和仓壁的位移。由于平房仓中的粮堆受到自重，通过体积力来加载。

（7）设小麦堆处于无压缩状态时堆长26 m，宽24 m，高9 m。在 Mesh 模块中对模型进行网格划分。小麦堆模型共被划分为9层，将每层粮堆部件均等分为1 m2的单元。本文采用节点法，提取各个节点应力应变值进行计算。

（8）在 Job 步骤中对构建的模型进行运算。

（9）待运算成功后，在Visualization步骤中查看运算结果提取所需数据。

1.6 平房仓中小麦堆孔隙率的计算

1.6.1 平房仓中小麦堆单元孔隙率的计算

小麦堆经修正剑桥模型的有限元计算后提取各个单元节点的应变值，其单元孔隙率为：

式中：Vv0为小麦堆初始（表层）单元孔隙的体积，m3；V0为小麦堆初始（表层）单元体积，m3；ΔVijk为第ijk个单元的体积增量，m3。

式（13）分子分母同除以V0得到：

式中：ε0为小麦堆的初始（表层）孔隙率。

1.6.2 平房仓中小麦堆的粮层平均孔隙率

有限元软件计算后提取小麦堆中各个单元的应变值，小麦堆的第i层平均孔隙率可由式（16）计算出：

式中：Vijk为第ijk个单元所在处的单元体积，m3；m为沿宽度方向的单元数；n为沿长度方向的单元数。

2 结果与分析

2.1 平房仓中小麦堆的孔隙率分布特性

由有限元模型计算后可以提取小麦堆中各个单元节点的应变值，而粮堆各单元的孔隙率可以根据式（15）计算得到。为了研究其同一粮层的小麦堆孔隙率分布特性，本文选取不同含水率小麦堆的第2、5、8层的孔隙率进行具体分析，结果如图2～图5所示（设平房仓的拐角为原点）。

由图2～图5可知，对于同一含水率，在相同粮层深度下的小麦堆的孔隙率分布是不均匀的。在同一层中，距仓壁距离增加，小麦的孔隙率减小，在该层的中心，在拐角处的孔隙率最大。不同粮层之间的小麦堆的堆积密度的分布也存在明显差异，在平房仓中小麦堆的粮层深度越深，该层的小麦孔隙率分布越不均匀。

图2 湿基含水率9.72%的小麦堆不同粮层的孔隙率分布值

图3 湿基含水率12.03%的小麦堆不同粮层的孔隙率分布值

图5 湿基含水率16.55%的小麦堆不同粮层的孔隙率分布值

在各个含水率中，同一粮层深度下的小麦堆的孔隙率分布是不均匀的。在同一粮层中，小麦堆的孔隙率由仓壁至粮层中间逐渐减小，在该层中间的孔隙率最小，在平房仓的拐角处孔隙率陡然增大并达到该层的最大值。不同层之间的孔隙率也存在明显差异。粮层深度增加，同一水平面位置处的孔隙率减小；粮层深度增加，粮层各处的孔隙率的差值增大。同一粮层各处的孔隙率差值很小，湿基含水率9.72% 的小麦堆底层（第8层）的孔隙率最大值为42.69%，最小值为42.58%；湿基含水率12.03%的小麦堆底层（第8层）的孔隙率最大值为43.50%，最小值为43.39%；湿基含水率14.61%的小麦堆底层（第8层）的孔隙率最大值为44.00%，最小值为43.90%；湿基含水率16.55%的小麦堆底层（第8层）的孔隙率最大值为44.46%，最小值为44.36%；孔隙率的差值约为0.1%。

图4 湿基含水率14.61%的小麦堆不同粮层的孔隙率分布值

2.2 平房仓中小麦堆各粮层的平均孔隙率

由有限元模型计算后，提取小麦堆中各单元节点的应变值，小麦堆粮层平均孔隙率可以根据式（16）计算得到，结果见图6。

2.2.1 平房仓中小麦堆粮层平均孔隙率与深度的关系

由图6可知，在同一含水率下，小麦堆的粮层平均孔隙率随粮层深度增加而减小。这是由于粮层深度增加，粮层压力增加，体积缩小，导致孔隙率减小。

2.2.2 平房仓中小麦堆粮层平均孔隙率与含水率的关系

由图6可知，对于同一粮层，粮层平均孔隙率随含水率的增加而增大。这是由于含水率增加，粮粒膨胀，引起孔隙率增加。

图6 小麦堆的孔隙率随粮层深度变化的实验值与预测值对比

2.3 平房仓中小麦堆粮层孔隙率与含水率和深度的关系模型

设平房仓中小麦堆（淮麦44号）粮层平均孔隙率与粮层深度和含水率的方程为：

式中：εd为小麦堆粮层平均孔隙率，%；ε0为小麦堆表层孔隙率，%；εmin为平房仓中小麦堆粮层平均孔隙率的最小值，%；h为小麦堆粮层深度，m。

依据表2的数据，小麦堆的表层孔隙率与含水率的关系可拟合为线性方程为：

式中：MC为小麦堆的湿基含水率，%。

依据表2的数据，小麦堆最小粮层平均孔隙与含水率的关系可拟合为线性方程：

将式（17）转化为式（20）。

将式（18）、式（19）及a、b、c、d的值代入式（17），得平房仓中小麦堆平均孔隙率与含水率和粮层深度的关系模型为：

如图6所示，由式（21）模拟得到的小麦堆各粮层的平均孔隙率（虚线）与实验计算的平均孔隙率相比相对误差小于0.15%，RMSE的值在0.07% ~ 0.15%，误差较小。这个结果表明式（21）表征平房仓中小麦堆粮层孔隙率与含水率和粮层深度的关系是有效的。

3 结 论

（1）在平房仓内的同一粮层中，小麦堆孔隙率随距仓壁的距离减小而增大，在粮仓的拐角处小麦堆孔隙率达到该层的最大值，在该层的中间位置出现最小值；粮层深度越深，该层的孔隙率分布越不均匀。

（2）平房仓中小麦堆粮层平均孔隙率随着粮层深度的增加而减小。

（3）平房仓内小麦堆粮层平均孔隙率随含水率的增加而增大。

（4）平房仓中小麦堆粮层孔隙率与粮层深度和含水率之间关系方程为：