王 颖， 程绪铎， 唐福元

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

我国是全球第一大小麦生产国。2019年我国小麦产量1.33亿t，消费量为1.28亿t[1]。小麦储藏在粮仓中，自身的重力与摩擦力、仓壁对小麦的压力与摩擦力引起了小麦块的变形与压缩[2,3]，导致小麦孔隙率减小、堆密度增大。Thompson等[2]研究发现小麦堆密度随着内部的压力变化而变化。小麦储藏在筒仓中，小麦堆的压应力分布是不均匀的，尤其是压应力随着粮层深度的增加而增加，因此筒仓中小麦堆的密度分布也是不均匀的，也是随着粮层深度的增加而增加。确定筒仓中小麦堆的密度分布值是准确计算筒仓中的粮食质量的关键参数，也是预测筒仓小麦堆通风阻力、小麦堆湿热传递规律的重要参数[4,5]。

通常用气体密度仪来测量粮堆的堆密度[6-9]。气体密度仪测量的只是小粮堆的堆密度(粮堆体积小于10 L)且测定的是无压缩的粮食密度，它无法被送达到粮堆深处，因此不能用来测量筒仓深处的粮食堆密度。我国一般使用容重仪测定粮食的密度，但测定的是无压缩的粮食密度，同样无法测定筒仓深处受压缩的粮食密度。Cheng等[10]建立筒仓内小麦分层压缩平衡方程，数值求解了筒仓中小麦各层的堆密度。这个分层压缩平衡方程是筒仓中各层小麦堆密度的平均值，不是某一粮层上各粮块的堆密度。

ABAQUS是国际上功能最强的大型通用有限元软件之一，拥有能够真实反映松软土体的修正剑桥模型，已应用于研究土体中的力学问题。小麦堆的力学特性与松软土体相似，本研究使用ABAQUS软件求解筒仓中的小麦堆的修正剑桥模型。通过有限元方法求解筒仓中的小麦堆的修正剑桥模型得到筒仓中小麦堆密度的三维空间分布值，即小麦堆密度随粮块深度(纵向)、径向、环向变化的分布值。

1 材料与方法

1.1 实验材料

小麦，烟农19号，产自江苏淮阴，人工剔除实验样品中破损籽粒和杂质，测得原始含水率为12.66%。实验样品的含水率分别调制为10.60%、12.66%、14.22%、16.13%。调制好含水率的样品装在塑料袋中封闭好，然后储藏在人工气候箱中，人工气候箱中的温度设置为4 ℃。实验前从人工气候箱中取出小麦样品，采用标准烘箱干燥法在130 ℃下干燥实验样品19 h测定它的含水率(ASAE Standards, 2001)，重复3次。

1.2 小麦堆应力与应变关系模型

储藏在筒仓中的小麦堆在外力作用下产生应力，而应力的作用又引起应变。要求解筒仓中小麦堆的应力分布与应变分布必须建立小麦堆应力与应变关系模型。修正剑桥模型将总应变分为剪切应变与体积应变，总应力分为平均正应力与剪切应力，适合求解筒仓中小麦堆的体积应变(应变较大)。故本研究选择修正剑桥模型来表征小麦堆的应力与应变关系。

修正剑桥模型中的应变增量包括体积应变增量dεv和剪切应变增量dεs；应力增量为平均正应力增量dp和剪切应力增量dq，q=σ1-σ3是广义剪切应力，p=(σ1+2σ3)/3是平均主应力。

修正剑桥模型的弹塑性矩阵的显式表达式为[11]：

(1)

Cpp=

(2)

Cpq=Cqp=

(3)

Cqq=

(4)

由式(2)～式(4)可见，修正剑桥模型有6个参数：e为初始孔隙比，ν为泊松比，E为弹性模量/kPa，M为临界状态应力比，κ为等向膨胀指数，λ为对数硬化模量；这6个参数可以通过三轴压缩与剪切实验而确定[11]。

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

1.3.1 初始孔隙比的确定

1.3.1.1 初始孔隙率ε的测定

小麦初始孔隙率ε是未压缩时小麦堆中孔隙体积与小麦堆体积之比，测定初始孔隙率是为了得到初始孔隙比。本研究使用LKY-1型粮食孔隙率测定仪(见图1)测定小麦的初始孔隙率。

图1 孔隙率测定仪示意图

如图1所示，两个容积相等的压力容器A和B，在容器B中装满样品并将其密封。关闭阀门2，打开阀门1、3，用空气压缩机向容器A中鼓入一定压力的气体。当压力表指针达到一定数值时，关闭阀门1，当压力值稳定后，记下压力表读数P1；关闭阀门3，然后打开阀门2，当容器A和B中压力平衡时，记下此时压力表读数P2。视空气为理想气体(标准大气压为Pa)，由理想气体等温过程的特性推出小麦初始孔隙率(未压缩)为：

(5)

1.3.1.2 初始孔隙比e的确定

初始孔隙比e指未压缩时小麦堆中孔隙体积与小麦堆中籽粒体积之比，小麦样品初始孔隙比(未压缩)可由初始孔隙率推导得到:

e=ε/(1-ε)

(6)

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

对样品进行轴向压缩实验[12]测定临界状态应力比M。在轴向压缩实验中，选定5个围压σ3：30、50、70、90、110 kPa。装好样品后，启动仪器对样品进行剪切(剪切应变速率为1.000 mm/min)，位移每增加0.4 mm，记录1次测力计读数和样品体积减少量，直至测力计读数出现最大值时，记下峰值时平均主应力p和广义剪切力q。对每个围压，实验重复3次。将这5组实验的破坏点(即最大主应力差)所对应p和q绘制在p-q平面中，经一元线性回归得直线的斜率M值。

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

对数硬化模量λ和等向膨胀指数κ的测定实验为各向等压实验[12]，使围压σ3从0 kPa增加至200 kPa，围压每增加5 kPa，记录1次样品体积减少量，待围压增至200 kPa，再依次卸载至0 kPa，同样每减小5 kPa，记录1次样品体积增加量。将加载曲线和卸载曲线上的p和所对应的孔隙比e绘制在e-lnp平面中，经一元线性回归得直线的斜率λ和κ值。

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

选定围压σ3分别为30、50、70、90、110 kPa，对样品进行轴向加卸载循环的常规三轴压缩实验[12]，测定弹性模量E。在常规三轴压缩实验中，粮堆样品呈圆柱形，轴向压力为σ1，围压为σ3，主应力差为Δσ=σ1-σ3，当主应力差达到最大值时粮堆破坏。对于每一个围压σ3，可测定一个最大主应力差。对选定的围压σ3，分级施加轴向压力，每级压力取试样最大主应力差的10%；施加第1级压力，并开动秒表，记录加压1 min后位移计的读数，每1 min施加一级压力，记录位移计读数，直至施加到第4级压力为止；在施加第4级压力1 min后记录位移计的读数，并逐级卸压，每1 min卸去1级压力，并记录卸压后1 min时位移计的读数，直至卸去所有施加的轴向压力；重复加卸载荷4次。

弹性模量E按第4次加卸载荷记录的数据计算公式为：

(7)

式中：Δhe为小麦样品的轴向弹性变形量/mm；hc为未轴向加载时小麦样品高度/mm。

如同测定等向膨胀指数κ一样，对小麦样品进行各向等压实验，体变弹性模量K按卸载记录的数据计算公式为：

(8)

式中：ΔV为小麦的弹性体积增量/m3；V为未加压力时小麦样品的体积/m3。

泊松比ν按公式计算：

(9)

1.4 筒仓中小麦堆密度分布的有限元求解方法

1.4.1 筒仓的几何尺寸与力学参数

筒仓的几何尺寸与力学参数见表1。

表1 筒仓几何尺寸及其材料的参数

1.4.2 筒仓与小麦堆单元网格划分

筒仓高31 m, 重力未施加时筒仓内小麦堆高度29 m。筒仓与小麦堆的形状和受力都是轴对称的，取筒仓与小麦堆的一径向平面表示筒仓与小麦堆，选择四边形单元划分该径向平面，四边形单元的面积为1×1=1 m2,每一个四边形单元代表一个截面为矩形的环形体。小麦堆模块被划分为29×5个单元；仓壁模块被划分成31×1个单元。

1.4.3 小麦堆与筒仓壁的接触模拟

仓壁内侧面与小麦堆外侧面有相互力的作用。接触面之间传递切向应力与法向应力，切向应力即摩擦力。在有限元分析中应考虑阻止两表面之间相对滑动而产生的摩擦力。库仑模型是常用的摩擦模型，通过摩擦系数来描述两表面间的摩擦行为。临界剪应力与法向接触压力之间的关系为：

τ=μp

(10)

式中：μ为摩擦系数；p为法向接触压力/kPa。

1.4.4 载荷及边界约束

筒仓中小麦堆受到自重、内摩擦、仓壁对小麦的抵抗等多种力的作用，小麦堆的重力由体积力进行描述，体积力是主动力。小麦堆内部作用力通过小麦堆的应力与应变的本构关系进行描述，小麦堆与仓壁之间的相互作用力由摩擦力与正压力进行描述。将小麦堆模型底部边界条件设置为U2=0，即对筒仓底板处节点竖直方向的位移约束。因为在储藏过程中，筒仓是不产生任何运动的，故将仓壁处节点位移完全约束，以限制其刚体位移。

1.4.5 筒仓中小麦堆各处密度的计算

小麦堆共有n×m个单元，n(29)为行数，m(5)为列数。重力未施加时各单元密度是相同的；重力加载后，每个单元在筒仓中的位置及其体积都发生了变化。若i表示行(纵向)编号，j表示列(横向)编号，则第i行j列的单元的密度可由公式计算：

(11)

2 结果与分析

2.1 初始堆密度(未压缩)

根据国家标准GB/T 5498—1985测定得到含水率分别为10.60%、12.66%、14.22%、16.13%的烟农19号表层密度(未压缩)分别为801.75、787.09、770.47、755.65 kg/m3。

2.2 修正剑桥模型参数

使用应变式三轴仪和粮食孔隙测定仪测定烟农19号的修正剑桥模型参数参数如表2所示。

表2 修正剑桥模型参数

2.3 筒仓中小麦堆的密度分布值

根据国家标准GB/T 5498—1985测定得到含水率分别为10.60%、12.66%、14.22%、16.13%的小麦堆(烟农19号)初始密度(表层密度)分别为801.75、787.09、770.47、755.65 kg/m3。使用修正剑桥模型和有限元软件ABAQUS计算得到筒仓中小麦堆的密度分布值，见图2。

图2 不同含水率小麦堆的密度与半径的关系

由图2可知，筒仓中小麦堆的密度随着粮层深度的增加而逐渐增大，而在筒仓拐角处的小麦堆密度随着粮层深度的增加而逐渐减小。在上部粮层，随着粮块与筒仓中心轴距离的增加，密度缓慢减小，减小值约为0.07%～0.1%。在中部粮层，随着粮块与筒仓中心轴距离的增加，密度几乎不变。在接近仓底层，随着粮块与筒仓中心轴距离的增加，密度逐渐减小，减小值约为0.4%～0.6%。

2.4 筒仓中小麦堆平均密度与粮层深度的关系

重力未施加时筒仓内小麦堆高度29 m。筒仓与小麦堆的形状和受力都是轴对称的，取筒仓与小麦堆的一径向平面表示筒仓与小麦堆，选择四边形单元划分该径向平面，四边形单元的面积为1 m2,每一个四边形单元代表一个截面为矩形的环形体。小麦堆模块被划分为29×5个单元；仓壁模块被划分成31×1个单元。这样小麦堆共有29行5列。筒仓中小麦堆第n层(行)平均密度按式(12)计算：

(12)

依据修正剑桥模型计算的数据，由式(12)计算得到小麦堆各粮层的平均密度，如图3所示。在小麦含水率为10.60%、12.66%、14.22%、16.13%时，小麦堆的粮层深度的变化范围分别为0～27.704 8、0～27.661 0、0～27.485 8、0～26.982 1 m；其平均密度的变化范围分别为800.75～822.08、786.02～807.61、769.32～792.69、754.33～783.23 kg/m3。在同一含水率下，小麦堆的粮层平均密度随着粮层深度的增大而逐渐增大，逐渐趋于恒定值。

图3 不同含水率下小麦堆的平均密度与粮层深度的关系

2.5 筒仓中小麦堆平均密度与含水率的关系

从图3可以看出，含水率为10.60%、12.66%、14.22%、16.13%的小麦堆在相同的粮层深度下的平均密度是随着含水率的增大而减小的。

2.6 筒仓中小麦堆平均密度与含水率、粮层深度关系模型

有限元计算的密度数据有这样的规律，当深度一定时密度与含水率大致成线性关系，当含水率一定时密度随深度增加而增大，增加率随深度增加而减小，故用1-eah3+bh2+ch表达密度随深度增加的变化趋势。因此，设小麦粮层平均密度与粮层深度、含水率的关系方程为：

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

(13)

式中：ρ为小麦堆粮层平均密度/kg/m3；ρ0为小麦堆初始密度/kg/m3；ρmax为小麦堆压缩粮层平均密度的最大值/kg/m3；h为小麦堆粮层深度/m。

表3 入仓小麦的参数

小麦堆初始密度与含水率的关系可拟合为线性方程(14)。

ρ0=892.35-8.56MCρ=0.99

(14)

式中：MC为小麦堆含水率/%。

通过式(12)计算得到的小麦堆粮层平均密度的最大值变化范围为822.08～783.33 kg/m3，其随含水率的增大而降低。小麦堆粮层平均密度的最大值与含水率的关系可拟合为线性方程(15)。

ρ0=897.99-7.20 MCR2=0.99

(15)

为了确定模型(13)中a、b、c这3个模型常数，将式(13)变换成式(16)：

(16)

以ln(1-(ρ-ρ0)/(ρmax-ρ0))为纵坐标，h为横坐标，按式(16)进行数据拟合从而得到模型常数a、b、c分别是-0.000 47、0.012、-0.24。

将a、b、c的值及式(14)和(15)代入式(13)，可得小麦堆粮层的平均密度与含水率、粮层深度的关系模型为：

ρ=892.35-8.56MC+(5.64+1.36MC)(1-e-0.000 47h3+0.012h2-0.24h)

(17)

图3中的散点是有限元方法计算的密度，曲线表示模型方程(17)计算的密度值。由图3可以看出，经模型方程(17)计算得到的小麦堆粮层密度与限元方法计算的粮层密度误差较小，对于含水率分别为10.60%、12.66%、14.22%、16.13%的小麦，模型方程(17)计算的小麦堆粮层平均密度与有限元方法计算的粮层平均层密度的误差小于0.5%，说明模型方程(17)有效。

2.7 实仓验证

实例一：中央储备粮温州直属库11、12、13、14号筒仓内径为10 m，入仓小麦的参数见表3。本研究通过有限元分析方法求解修正剑桥模型计算得到小麦的粮层密度分布值，再由积分方程计算出储粮总质量，并与筒仓中实际储粮质量进行比较，结果见表4。

表4 修正剑桥模型计算结果与实际储粮数量的比较

表5 入仓小麦的参数

实例二：中央储备粮金华直属库1、2号浅圆仓内径为25 m，入仓小麦的参数见表5。本研究通过有限元分析方法求解修正剑桥模型计算得到小麦的粮层密度分布值，再由积分方程计算出储粮总质量，并与筒仓中实际储粮质量进行比较，结果见表6。

表6 修正剑桥模型计算结果与实际储粮数量的比较

由表4与表6可以看出，本研究的修正剑桥模型计算结果与实际账面数误差小于0.3%，说明本研究的修正剑桥模型计算的密度分布值可靠性较高。

入仓方式对粮食分级有影响，粮食分级对密度的分布有影响，粮食入仓以后，由于粮食分级，同一层的粮食密度差别增大，但一层的平均密度变化很小，所以对计算粮食总质量影响较小，本研究的模型计算结果表明影响较小。粮食含杂量对压缩密度有影响，本研究选择的实仓中的小麦含杂量较小(0.2%～1.0%)，所以对计算粮食总质量影响较小，本研究的模型计算结果表明影响较小。

3 结论

筒仓中小麦堆的密度随着粮层深度的增加而逐渐增大，而在筒仓拐角处的小麦堆密度随着粮层深度的增加而逐渐减小。在上部粮层，粮块密度随着粮块与筒仓中心轴距离的增加而缓慢减小，减小值为0.07%～0.10%。在中部粮层，粮块密度随着粮块与筒仓中心轴距离的增加几乎不变。在接近仓底粮层，粮块密度随着粮块与筒仓中心轴距离的增加而逐渐减小，减小值为0.4%～0.6%。小麦堆的粮层平均密度随着粮层深度的增大而逐渐增大，逐渐趋于恒定值。相同深度下的粮层平均密度随着含水率的增大而减小。筒仓中小麦堆粮层平均密度与含水率、粮层深度之间关系方程为：ρ=892.35-8.56MC+(5.64+1.36MC)(1-e-0.000 47h3+0.012h2-0.24h)。实仓验证得到，由修正剑桥模型计算的筒仓中小麦的总质量与实际粮仓的账面小麦总质量误差小，相对误差在0.14%～1.25%之间。