基于有限元方法的平房仓中小麦堆的密度分布研究

陈 雪 程绪铎 胡美珠

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

当小麦储藏在平房仓中时，粮堆受到自身质量、内摩擦力以及仓壁与粮堆之间的摩擦力作用，产生压应力和切应力[1]。小麦的籽粒松软，粮堆各处的应力作用导致粮堆产生形变和体积压缩，体积压缩导致粮堆密度的增加，平房仓内小麦堆的应力分布是不均匀的，所以小麦堆密度分布也是不均匀的。每年我国均需投入大量人力、物力和资金开展清仓查库工作[2]。粮食储藏数量检查是粮食库存检查的一项重要内容[3]，即通过一定的手段获取粮仓中所储存粮食的数量，主要包括称重法和体积密度法[4]。称重法工作量大、效率低，在大规模的库存检查中很难广泛使用。体积密度法是指通过测量粮堆的体积和平均密度计算粮仓中粮食数量的检查方法。现有的体积密度法是使用粮食的表层密度乘以修正系数来获取粮堆的平均密度，修正系数是凭经验给出的，此方法误差较大。因此，只有准确计算或测量出平房仓中的密度分布值才能准确快捷地计算平房仓中小麦堆的总质量。

前人通过在粮仓底部安装了压力传感器[5,6]，测定粮仓底部的压力分布，从而估算出粮堆的总质量，此方法成本高且只能得出粮堆总质量而无法得出粮仓内粮堆密度分布值。程绪铎等[7]开发了一个分层压缩模型来预测平房仓中的粮堆密度和压力分布，但其假设平房仓中同一粮层的密度是相同的，而实际中平房仓内同一粮层的堆密度分布是不均匀的，该分层压缩模型不能准确给出平房仓中的粮堆的密度分布值。

自20世纪70年代开始，有限元法被用来研究粮仓中物料的应力分布及仓壁的压力问题。Jofriet[8]采用线弹性模型，用有限元方法研究了直筒仓内物料静态时的应力分布问题。Ooi[9]使用有限元方法研究了柔性筒仓内散体的应力分布。国际上一些研究人员构建了许多本构方程，例如：Mohr-Coulomb模型、Drucker-Prager模型和Lade-Duncan模型等，通过有限元方法研究了粮仓中物料的应力分布问题。这些研究很好地计算出了储藏在筒仓中压缩性很小的沙子、碎石、工业物料的应力分布，这些弹塑性本构方程也可以较准确地求解出粮仓内粮堆的应力分布。粮食与土壤、沙子等物料的力学特性有很大的差异，小麦籽粒松软，籽粒间的孔隙大，小麦堆在受压时体积变化大且体积变化同时发生在弹性变形和塑性变形两个阶段。但研究所使用的本构方程均假定物料的体积变化很小且只发生在弹性变形阶段，形状变化发生在塑性变形阶段。因此这些模型难以反映小麦堆变形时的应力与应变关系，用它们求解平房仓内小麦堆的应变分布是不准确的。1968年，Roscoe[10]提出了修正剑桥模型(Modified Cam Cay)，修正剑桥模型既考虑了弹性阶段的体积与形状变化又考虑了塑性阶段体积与形状的变化，修正剑桥模型以体积应变与剪切应变作为应变量，且屈服面为平均主应力屈服面，适合松软的散粒体的大体积变形特征，被广泛应用于土壤力学[11,12]。

目前鲜有将修正剑桥模型应用于平房仓中小麦堆密度分布的研究。本实验将使用有限元软件求解修正剑桥模型计算平房仓内的小麦堆密度分布值。

1 材料与方法

1.1 实验材料

小麦，品种名为淮麦44号，产地江苏南京，原始含水率为12.03%w.b., 小麦籽粒的大直径、中直径和小直径分别为5.88、3.07、2.79 mm。取 10 kg小麦样品，放置在太阳光下晾晒 24 h，测得其含水率为 9.72% w.b.，再取 20 kg小麦样品，分成 2 份，各自加适当的蒸馏水，然后将其装入塑料袋中封闭好，放置恒温箱中，温度选定 5 ℃，1周后从恒温箱中取出，测得2份小麦样品含水率分别为16.55%w.b.和14.61%w.b.。用标准烘箱干燥法将 10 g 样品在 130 ℃下干燥 19 h，重复 3 次，测定小麦的含水率。这样，实验用的小麦的含水率分别为 16.55%、14.61%、12.03%和9.72% w.b.。

1.2 实验仪器

HG202-2(2A/2AD)电热干燥箱，LKY-1型粮食孔隙率测定仪，SLB-6A型应变控制式三轴仪，HGT-1000A 容重仪。

1.3 修正剑桥模型

修正剑桥模型是一个应力-应变关系的弹塑性本构模型[13]。受到外界力作用时，物体产生的应变增量包括体积应变增量dεv和剪切应变增量dεs，而它们分别是由弹性应变增量和塑性应变增量组成：

(1)

(2)

修正剑桥模型为[11]：

(3)

(4)

(5)

(6)

(7)

式中：q为广义剪切力/kPa；p为平均主应力/kPa；K=E/(3(1-2ν)),K为体积模量/kPa ；G=E/(2(1+ν)),G为剪切模量/kPa。

由式(4)～式(7)可知，修正剑桥模型的参数为临界状态应力比M、等向膨胀指数κ、对数硬化模量λ、弹性模量E、泊松比υ、初始孔隙比e0。

1.4 小麦的修正剑桥模型参数测定

使用 HGT-1000A 容重器测定小麦堆无压缩密度(粮仓表层小麦密度)ρ0；使用 LKY-1 型粮食孔隙率测定仪测定小麦堆的修正剑桥模型参数初始孔隙比e0；使用 TSZ-6A 应变控制式三轴仪测定小麦的修正剑桥模型参数E、υ、M、κ、λ。

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

混凝土仓壁可视为刚性体，在有限元分析过程中不产生变形。平房仓尺寸参考南京铁心桥国家粮食储备库中粮仓的实际大小，具体仓壁材料的几何与力学参数见表1。

表1 平房仓及其材料的参数

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

在Part步骤中采用三维可形变模型(3D Deformable)建模，创建平房仓和小麦堆部件，尺寸见表1；在Property步骤中分别设定仓壁的材料属性和粮堆的塑弹性属性(修正剑桥模型选择塑性模型Clay Plasticity和弹性模型Porous Elastic)；在Assembly 步骤中将平房仓仓壁和粮堆部件组合起来；在 Step步骤中设定分析步； Interaction步骤是定义各个部件之间的相互接触关系。粮堆与仓壁之间存在相互作用，在ABAQUS中设定一个完整的接触必须包括两个部分：一是接触对，选择点对面离散方法，选择仓壁内侧面为主接触面，粮堆外侧面为从属接触面；二是定义相互作用力，两个接触面之间会产生摩擦力。其中，库仑模型是较为常用的摩擦模型，即通过摩擦系数来描述两表面间的摩擦关系。

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

在 Mesh 模块中对模型进行网格划分。小麦堆模型共被划分为9层，将每层粮堆部件均等分为1 m2的单元(图1)。本实验采用节点法，提取各个节点应力应变值进行计算；在Job步骤中对构建的模型进行运算；待运算成功后，在Visualization步骤中查看运算结果提取所需数据。

建立有限元模型，平房仓内小麦堆长26 m，宽24 m，高9 m(重力未施加，粮堆无压缩状态)。

1.7 小麦堆节点密度的计算

有限元软件计算后提取小麦堆中各个节点的应变值，小麦堆节点密度可以由式(8)计算得到。

(8)

1.8 平房仓中小麦堆的粮层平均堆密度

限元软件计算后提取小麦堆中各个节点的应变值，小麦堆的第i层平均堆密度可由式(9)计算。

(9)

式中：ρi为第i层小麦层的平均堆密度/kg/m3；Vijk为ijk个节点所在处的单元体积/m3；m是沿宽度方向的单元数，n是沿长度方向的单元数。

2 结果与分析

2.1 小麦的修正剑桥模型参数及无压缩密度

实验测定的修正剑桥模型的参数及无压缩密度见表2。

表2 修正剑桥模型参数表及无压缩密度

2.2 平房仓中小麦堆的密度分布

2.2.1 同一粮层小麦堆密度分布

平房仓中的小麦堆是一个轴对称的长方体，选择其中的四分之一小麦堆来分析其堆密度。由式(8)计算出平房仓内小麦堆的各节点的堆密度，本文给出了平房仓中小麦堆(含水率为12.03% w.b.)的第2、5、8层的堆密度分布，见图1(设平房仓的拐角为原点)。

由图1可知，平房仓内在同一粮层的小麦堆的堆密度分布是不均匀的。在同一粮层中，随着距仓壁越远，小麦的堆密度越大，并且在该层的中间趋于稳定值。在同一粮层中，在平房仓的拐角附近的小麦的堆密度迅速下降，拐角处的小麦堆密度达到该层的最低值。不同粮层之间的小麦堆密度也存在明显差异，小麦粮层深度越深，其拐角的堆密度变化越大，该粮层密度分布越不均匀。

图1 含水率12.03%w.b.的小麦堆不同粮层的堆密度

2.2.2 小麦堆的粮层平均堆密度

根据式(9)计算得到小麦的4个不同含水率各粮层平均密度，见表3。

表3 不同含水率的小麦堆的不同粮层平均密度

由表3可知，在平房仓内，在同一含水率下，小麦堆的粮层平均密度随深度增加而增加，但增加率随深度增加而减小。当含水率分别为9.72%、12.03%、14.61%、16.55% w.b.时，小麦的粮层深度分别从0增至7.62、7.65、7.73、7.81 m时，小麦的堆密度分别从802.3增至841.6 kg/m3、787.3增至829.8 kg/m3、752.5增至794.9 kg/m3和739.8增至782.7 kg/m3。

2.2.3 小麦堆的平均堆密度与含水率的关系

由表3中的数据计算出各个含水率下的平房仓内小麦的平均堆密度。含水率分别为9.72%、12.03%、14.61%、16.55% w.b.时，小麦的平均堆密度分别为827.8、816.1、782.1、770.5 kg/m3，其与无压缩时的密度相比，增加率分别为4.39%、4.75%、5.12%、5.33%。平房仓中小麦堆的平均密度随含水率的增加而减小，平均密度相对于无压缩密度的增加率随含水率增加而增大。

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

设关于平房仓中小麦堆(淮麦44号)粮层平均密度与粮层深度、含水率的关系方程见式(10)。

ρi=ρ0+(ρmax-ρ0)(1-eah3+bh2+ch+d)

(10)

式中：ρi为小麦堆第i层平均密度/kg/m3；ρ0为小麦堆表层平均密度/kg/m3；ρmax为小麦堆最深层平均密度/kg/m3；h为小麦堆粮层深度/m。

首先求出小麦的表层密度和最深粮层平均密度与于含水率的关系式。由表3中的数据，小麦堆表层密度与含水率的关系拟合为线性方程(11)。

ρ0=-9.686 6MC+898.61，R2=0.978 3

(11)

式中：MC为小麦的含水率/%w.b.。

平房仓中小麦堆最深层平均密度与含水率的关系可拟合为线性方程见式(12)。

ρmax=-9.2213MC+934.22，R2=0.968 5

(12)

将式(10)两端取对数式：

(13)

将a、b、c、d的值和式(11)、式(12)代入式(10)，得到平房仓中小麦堆粮层平均密度ρ与含水率MC、粮层深度h的关系模型为：

ρ=-9.688 66MC+898.61+(35.61+0.465 3MC)(1-e-0.016 5h3+0.150 5h2-0.614 3h-0.01)

(14)

图2为由式(14)预测的小麦粮层的平均密度与与含水率、粮层深度的关系曲线。方程(14)的预测值与有限元方法的计算值的相对误差低于0.05%，RMSE的值在0.9～2 kg/m3之间，误差较小。

图2 粮层平均密度与含水率、粮层深度的关系

2.4 实仓验证

以南京铁心桥国家粮食储备库提供的小麦的实仓数据作为验证参数[14]，所得到账面实际小麦的总质量与修正剑桥模型计算得到小麦的总质量进行对比，见表4。

用本修正剑桥模型计算所得的小麦总质量与南京铁心桥国家粮食储备库提供的实仓中小麦总质量比较接近，误差在0.93～1.64%之间，误差较小，说明修正剑桥模型计算出的密度分布值的精度高。

表4 修正剑桥模型模型计算的小麦堆质量与实际储粮质量的对比

3 讨论

3.1 装粮方式对粮堆密度的影响

不同的装粮方式粮食下落的速度不同，下落的粮食对粮堆产生的冲击力不同，粮食落点处的密度与冲击力的大小有关。对于小型粮仓或装粮容器，由于粮堆高度不高，装粮方式对装粮后的粮堆密度有一定的影响。对于高大的平房仓，尽管粮食冲击力改变了粮食落点密度(相对于粮食下落速度为零的粮堆落点密度)，后续下落的粮食压住了先前的粮食落点，随着粮食的不断下落，落点上方的压力越来越大，落点上方的压力改变粮食落点密度，这种改变取代了冲击力对粮食落点密度的改变量，即冲击力对粮食落点密度的改变量有压力引起的密度改变量代替了。对于4 m以上堆高的粮堆，除了粮堆表层，装粮方式对粮堆各处的密度影响很小，可以略去不计。所以本研究可适用于各种装粮方式。

3.2 平方仓中粮堆密度分布值的验证

本研究计算出平房仓内粮堆的密度三维分布值。目前还无法直接测量出粮堆内部各点的密度值，因此，无法使用实验来验证平房仓内粮堆各处的密度值。本研究以南京铁心桥国家粮食储备库的平房仓为例，将计算出的粮堆各点粮块的密度乘以粮块体积，然后求和得到粮堆总质量。将这个总质量与实际账面储粮质量比较，两者几乎一致(误差≤1.64%)，间接说明本研究的计算结果是可靠的。

4 结论

以剑桥修正剑桥模型为小麦堆应力与应变关系本构模型，有限元方法求解平房仓内小麦堆(长26 m，宽24 m，高7.62～7.81 m。)的修正剑桥模型，从而计算出平房仓中的小麦堆密度分布值。结果表明，小麦储藏在平房仓中时，其堆密度分布是不均匀的。小麦平均层密度随着粮层深度的增加而增大；当含水率分别为9.72%、12.03%、14.61%、16.55% w.b.时，小麦的粮层深度和堆密度增加。在同一粮层下，小麦的堆密度随距仓壁的距离减小而减小，在粮仓的拐角处粮堆密度达到该层的最小值，在该层的中间位置密度到达最大值，当粮层深度越深，该层的粮堆密度分布越不均匀。平房仓内小麦的平均堆密度随着含水率的增加而减小，平均堆密度相对于无压缩密度的增加率随含水率增加而增大。平房仓中小麦堆平均层密度与粮层深度、含水率之间关系方程为：ρ=-9.689MC+888.61+(35.61+0.4653MC)(1-e-0.016 5h3+0.150 5h2-0.614 3h-0.01)。

实仓验证得到，由该模型计算的平房仓中小麦的总质量与实际粮仓的账面小麦总质量误差小，相对误差在0.94%～1.64%之间。