基于有限元分析的筒仓中稻谷重量的测量方法

摘要

本发明提供基于有限元分析的筒仓中稻谷重量的测量方法，涉及粮食测量领域。该测量方法，包括如下步骤：对稻谷进行三轴轴向压缩试验和各向等压压缩试验，计算得到稻谷修正剑桥模型的参数：临界状态应力比、等向膨胀指数、对数硬化模量和泊松比；将筒仓中未压实状态下的稻谷堆划分成N层；每层稻谷堆由中心向外划分成中心轴线与筒仓中心轴线共线的M个单元；各单元的应变与应力关系符合修正剑桥模型；应用ABAQUS求解修正剑桥模型，得到压实状态下稻谷堆各单元的高度，计算压实状态下N×M个单元体积应变,计算N×M个单元密度计算每层稻谷堆的重量，计算出筒仓内稻谷堆的总重量W。本发明测量方法，简单、耗时短、测量结果准确、精度高。

说明书

技术领域

本发明涉及粮食测量领域，具体涉及一种基于有限元分析的筒仓中稻谷重量的测量方法。

背景技术

筒仓是储藏稻谷的主要仓型之一，是一种机械化程度最高的储粮容器，主要用于储藏散粮。稻谷是一种典型的散粒体，当储藏于筒仓中时，其静态特性表现为弹塑性，稻谷会受到自重、内摩擦力以及仓壁、仓底的支持力。稻谷堆的内部由于这些力的作用产生应力，从而产生弹性和塑性形变，稻谷堆的体积缩小，密度增大；随着粮层深度的增加，稻谷堆的压应力与切应力增大，体积应变增大，密度也随之增大；由于仓壁摩擦的影响，同一粮层不同圆环处的应力不同，应变不同，密度也不同。因此，筒仓中稻谷堆密度不是常量，而是空间位置的函数。粮食储藏数量检查是一项重要的粮食库存检查内容，目前的检测方法主要为称重法和体积密度法。称重法是指通过使用符合法定计量标准的衡器来称量粮食重量的检查方法，其优点是其测量结果一般能够准确反映被查粮食的实际数量，但是称重法具有工作量大、效率低等缺点，难以广泛应用于大规模的库存检查中。体积密度法是指通过测量计算粮堆的体积和粮食的表层密度，粮食的表层密度乘以修正系数来获取粮堆的平均密度，体积乘以平均密度得出仓内粮食总重量的检查方法。

稻谷堆的平均密度因筒仓的几何尺寸、粮堆的高度、压缩特性、内摩擦特性及筒仓壁与粮食的摩擦特性的不同而不同，到目前为止，计算平均密度的修正系数是凭经验估计的。测量粮堆的密度分布的新方法正在研究之中，有研究者利用微波检测粮堆的密度，这个方法是通过测量微波在粮堆中通过时的介电常数，由介电常数与密度的关系而获取粮堆的密度，但这种方法不能测量筒仓深处的密度，因微波到达粮堆深处易被干扰，目前只能测量两米深处的粮堆密度。也有研究者在筒壁和底部安装力传感器，由仓壁和底部应力分布推出储粮总重量，但这种方法存在两个问题，一是由仓壁和底部应力分布堆出储粮总重量需要建立精准的数学模型，否则测算的误差大，二是成本大，难以推广。

发明内容

本发明的目的是提供一种基于有限元分析的筒仓中稻谷重量的测量方法，该方法简单、耗时短、测量结果准确、精度高。

本发明的目的采用如下技术方案实现。

基于有限元分析的筒仓中稻谷重量的测量方法，其特征在于包括如下步骤：

(1)对稻谷进行三轴轴向压缩试验和各向等压压缩试验，计算得到稻谷修正剑桥模型的参数：临界状态应力比M、等向膨胀指数κ、对数硬化模量λ和泊松比υ；

(2)测量筒仓的内径d(单位是m)、高度H(单位是m)；测量筒仓中压实状态下稻谷堆的高度h(单位是m)，估算未压实状态下稻谷堆的高度h′(单位是m)；将筒仓中未压实状态下的稻谷堆划分成N层；每层稻谷堆由中心向外划分成中心轴线与筒仓中心轴线共线的M个单元，所述M个单元由1个圆柱和M-1个圆环柱组成，所述圆柱处于每层的中心；其中每层稻谷堆的各单元按照由中心向外依次编号其列数；共得到N×M个单元；各单元的应变与应力关系符合修正剑桥模型；应用大型有限元软件ABAQUS求解修正剑桥模型，得到压实状态下稻谷堆各单元的高度，根据式(A)分别计算筒仓内压实状态下稻谷堆的N×M个单元体积应变

式(A)，式(A)中，n为层数(是自然数)，m为列数(是自然数)，是第n层、m列单元的体积应变，和是第n层、m列单元的三个主应变；

根据式(B)分别计算N×M个单元密度的ρ_nm：(n＝1,2,…,N；m＝1,2,…,M)式(B)，

式(B)中，n为层数(是自然数)，m为列数(是自然数)，ρ_nm是第n层、m列单元的密度，kg/m³；ρ₀是稻谷堆表层的密度，kg/m³；是第n层、m列单元的体积应变；

(3)根据式(C)计算每层稻谷堆的重量，

式(C)中，W_n是第n层稻谷堆的重量，kg；；ρ_nm是第n层、m列单元的密度，kg/m³；h_nm是压实状态下稻谷堆第n层、m列单元的高度；R_m表示第m列单元的外径(当m＝1时，R₁是指圆柱的直径)；n为层数(是自然数)，m为列数(是自然数)，n＝1,2,…,N；m＝1,2,…,M；

采用式(D)计算出筒仓内稻谷堆的总重量W：

式(D)中W_n是第n层稻谷堆的重量，n＝1,2,…,N。

每层稻谷堆的各单元按照由中心向外依次编号其列数，具体编号方法如下：中心处的圆柱为第1列，然后依次向外编号为第2列，直至第M列。将每层稻谷堆划分成5个单元的示意图见图1。

在本发明中，稻谷堆的压实状态是指稻谷堆装入筒仓后的堆积状态。稻谷堆的未压实状态是指重力未施加时稻谷堆的堆积状态。

在本发明中，所述稻谷堆表层是指稻谷堆表面至0.5米深度处的稻谷层。稻谷堆表层所受重力可以忽略。

在本发明中，N取值为：h′≤N≤2h′，h′为未压实状态下稻谷堆的高度，m；M取值为：0.5R≤M≤2R，R为稻谷堆的半径，m。

在本发明中，h′取值如下：当2m<h≤5m时，1.04<h′/h≤1.06；当5m<h≤10m时，1.06<h′/h≤1.10；当10米<h≤15米时，1.10<h′/h≤1.14；当15米<h≤20米时，1.14<h′/h≤1.17；当20米<h≤25米时，1.17<h′/h≤1.19；当25米<h≤30米时，1.19<h′/h≤1.20；其中h是筒仓中压实状态下稻谷堆的高度。上述h′的取值是采用如下方法估算得到的：根据压实状态下稻谷堆的高度h，假设一h′值，由上述基于有限元分析的筒仓中稻谷重量的测量方法，得到压实状态下稻谷堆的N×M个单元的高度(轴向的长度)，从而得到压实状态下稻谷堆的高度h，若则假设的h′的值对了；若则增大h′，若则减小h′；再一次由上述基于有限元分析的筒仓中稻谷重量的测量方法，得到压实状态下粮堆的高度再比较与h，直到为止。

在本发明中，步骤(1)中稻谷的修正剑桥模型参数M、κ、λ和υ通过SLB-6A型应变控制式三轴仪(南京土壤仪器厂有限公司制造)进行三轴轴向压缩和各向等压压缩试验后计算得到。

采用三轴轴向压缩试验确定临界状态应力比M，包括如下步骤：装入稻谷样品，分别设定围压为30、50、70、90、110kPa，启动SLB-6A型应变控制式三轴仪施加轴向力对样品进行轴向压缩，稻谷样品轴向位移每增加0.4mm，记录一次测力计的应力值读数和样品体积减少量，直至测力计读数出现峰值时，记下测力计读数峰值q及对应的平均压应力p，以p为横坐标q为纵坐标作图，经一元线性回归得直线的斜率，该斜率即为临界状态应力比M；

采用各向等压试验确定对数硬化模量λ和等向膨胀指数κ，包括如下步骤：装入稻谷样品，使围压σ₃从0kPa按照每次增加5kPa依次增加至200kPa，记录加载过程中围压σ₃和对应的样品体积减少量；再按照每次减小5kPa将围压由200kPa依次卸载至0kPa，记录卸载过程中围压σ₃和对应的样品体积增加量；绘制加载曲线和卸载曲线；经一元线性回归，将加载曲线的斜率作为λ，将卸载曲线的斜率作为κ；加载曲线和卸载曲线的横坐标为lnp,纵坐标为e，e是稻谷堆孔隙比，按下式计算：e＝V/[V₀(1-ε₀)]-1；式中，V₀是稻谷样品初始体积，m³；V是稻谷样品受压体积，m³；ε₀是稻谷堆表层孔隙率；p是平均压应力,kPa。

采用三轴轴向压缩试验测量弹性模量E和泊松比υ，包括如下步骤：装入稻谷样品，设置围压，轴向压缩样品直至破坏，测定最大主应力差(破坏时轴向压应力与围压差值)；将破坏的稻谷样品取出，再重新装入新的稻谷样品，进行加压实验，分4级施加轴向压力，第1级压力为所述最大主应力差的1/10，以后每级压力比前一级压力增加所述最大主应力差的1/10；每施加1级压力，记录施加的轴向压力和每次加压后1分钟时样品的轴向位移，每隔1分钟施加一级压力，施加到第4级压力为止；逐级卸压至施加的轴向压力全部卸去，分4级卸压，卸压顺序与加压顺序相反；重复加压卸压4遍后，轴向加压直至破坏；制作最后一次加压卸压实验中的加载曲线和卸载曲线，计算弹性模量E，E＝ΔP/(Δh_e/h_c)；其中，ΔP为轴向压力的递增或递减值，kPa；Δh_e为稻谷的轴向弹性变形量，mm；h_c为稻谷样品装样后的高度，mm；采用三轴各向等压试验确定稻谷样品体变弹性模量B，按公式B＝Δσ/(Δν/V)计算，Δσ是稻谷样品围压增量，kPa；Δv是稻谷样品体积压缩量，m³；V为稻谷样品装样后的体积，m³；泊松比ν按公式ν＝(3B-E)/6B计算。

在本发明中，步骤(2)中应用大型有限元软件ABAQUS求解修正剑桥模型，包括如下步骤：在ABAQUS的Part模块中，采用轴对称CAX4单元，创建筒仓和稻谷部件；在Property模块中，定义筒仓的钢材材料属性和稻谷堆的属性，选择修正剑桥模型表征稻谷堆的应变和应力关系，然后将各材料属性指派到各个部件中；将筒仓和稻谷堆两个部件在Assembly模块中装配成整体；在Step模块中创建分析步；在Interaction模块中定义两个部件之间的约束关系，通过摩擦系数来描述两表面间的摩擦行为；在Load模块中定义荷载以及边界条件；在Mesh模块中按照部件对结构(稻谷堆与筒仓)进行网格划分，将筒仓中稻谷堆划分成所述N×M个单元；在Job模块中提交作业进行分析；分析结束后，在Visualization模块中查看压实状态下稻谷堆各单元的三个主应变及各单元的高度h_nm，根据计算第n层、m列单元的体积应变式中，是第n层、m列单元的三个主应变；h_nm是指压实状态下第n层、m列单元的高度。

创建筒仓和稻谷部件具体是指定义筒仓内径、仓高、仓壁厚度、粮堆未压实高度，筒仓部件是矩形平面，该矩形平面是筒仓外径剖面被中心轴线分割成的一半；稻谷堆部件是矩形平面，该矩形平面是圆柱形稻谷堆直径剖面被中心轴线分割成的一半。

定义筒仓的钢材材料属性是指定义钢材材料的弹性模量与泊松比，定义稻谷堆的属性是指定义稻谷堆的修正剑桥模型参数。

在Interaction模块中定义两个部件之间的约束关系为摩擦接触关系，通过稻谷与筒仓仓壁之间的摩擦系数μ来描述两表面间的摩擦行为。

在Load模块中定义荷载是体积力，边界条件是：与筒仓底接触的稻谷无垂直位移，圆柱形稻谷堆中心轴线上的稻谷无水平位移。

在Load模块中荷载通过采用体积力来描述，其他作用力通过稻谷堆与仓壁、仓底之间的相互作用来描述。

在Interaction模块中定义稻谷堆与仓壁之间的相互作用通过设置接触来模拟，接触对选择点对面离散方法，主控面选定仓壁内侧，从属面选定粮堆模型外侧面。接触面相互作用力学模型选用常用的摩擦模型，通过摩擦系数来描述两表面间的摩擦行为。

在确定稻谷修正剑桥模型的参数前，测定稻谷与筒仓仓壁之间的摩擦系数μ、稻谷堆的内摩擦角稻谷堆表层的密度ρ₀以及稻谷堆的含水率MC；测定稻谷堆表层的孔隙率ε₀，并由ε₀计算得到稻谷堆表层的孔隙比

稻谷堆的内摩擦角稻谷与筒仓仓壁之间的摩擦系数μ采用直剪仪进行测定；测定稻谷堆的内摩擦角时，将稻谷置于上下试样盒之间，在稻谷上施加垂直压应力σ，然后施加水平推应力τ，使稻谷在上下盒之间的水平面上发生剪切直至破坏；根据公式计算内摩擦角其中，σ为在稻谷上施加的垂直压应力，kPa；τ为在稻谷上施加的水平推应力，kPa；tg是正切函数；是内摩擦角，°；测定稻谷与筒仓仓壁之间的摩擦系数μ时，上试样盒填充稻谷，下试样盒填充不锈钢板或混凝土板，在稻谷上施加垂直压应力σ，然后施加水平推应力τ，使稻谷在接触面上发生剪切直至破坏，根据μ＝τ/σ计算摩擦系数μ，σ为在稻谷上施加的垂直压应力，kPa；τ为在稻谷上施加的水平推应力，kPa。

使用LKY-1型粮食孔隙率测定仪(南京土壤仪器厂有限公司制造)测定稻谷堆表层的孔隙率ε₀，稻谷堆表层的孔隙比为

在本发明中，所述稻谷堆表层的密度ρ₀采用如下方法测量：在稻谷堆表层内取样，测定样品的质量m₀和体积V₀，根据式(3)ρ₀＝m₀/V₀计算得到ρ₀。

在本发明中，含水率是指粮食中所含水的质量与粮食总质量的比值，用百分数表示。根据国标GB/T5497-85，含水率采用105℃恒重法测定，并按下式计算稻谷的含水率：w＝(m₁-m₂)/(m₁-m₀)×100％，其中：m₀为铝盒质量，g；m₁为烘前稻谷样品和铝盒质量，g；m₂为烘后稻谷样品和铝盒质量，g；w为含水率，w.b.。

有益效果：本发明测量筒仓中稻谷重量方法简单、耗时短、测量结果准确。采用本发明方法测量了2个仓中的稻谷储藏重量，测量值与粮重实际账面数几乎一致，误差为分别为0.32％和2.00％，说明本发明测量方法准确，精度高。

具体实施方式

实施例1采用本发明方法测定北京市小汤山实验基地7号筒仓中稻谷的储藏重量

采用本发明方法测定北京市小汤山实验基地7号筒仓中稻谷的储藏重量。该筒仓内空腔为圆柱体状。

(1)小汤山实验基地7号筒仓的内径d、高度H见表1。筒仓内稻谷堆压实状态下的高度为7.21m(测定值)；筒仓内未压实状态下稻谷堆的高度h′，根据当5m<h≤10m时，1.06<h′/h≤1.10，h′取值为7.7m。

表1筒仓及其材料的参数

(2)测定稻谷与筒仓仓壁之间的摩擦系数μ、稻谷堆的内摩擦角稻谷堆表层的密度ρ₀以及稻谷堆的含水率MC；测定稻谷堆表层的孔隙率ε₀，并由ε₀计算得到稻谷堆表层的孔隙比

稻谷堆的内摩擦角稻谷与筒仓仓壁之间的摩擦系数μ采用直剪仪进行测定；测定稻谷堆的内摩擦角时，将稻谷置于上下试样盒之间，在稻谷上施加垂直压应力σ，然后施加水平推应力τ，使稻谷在上下盒之间的水平面上发生剪切直至破坏；根据公式计算内摩擦角其中，σ为在稻谷上施加的垂直压应力，kPa；τ为在稻谷上施加的水平推应力，kPa；tg是正切函数；是内摩擦角，°；测定稻谷与筒仓仓壁之间的摩擦系数μ时，上试样盒填充稻谷，下试样盒填充不锈钢板或混凝土板，在稻谷上施加垂直压应力σ，然后施加水平推应力τ，使稻谷在接触面上发生剪切直至破坏，根据μ＝τ/σ计算摩擦系数μ，σ为在稻谷上施加的垂直压应力，kPa；τ为在稻谷上施加的水平推应力，kPa。测定结果：稻谷与筒仓仓壁之间的摩擦系数μ为0.35，稻谷堆的内摩擦角为30°。

使用LKY-1型粮食孔隙率测定仪(南京土壤仪器厂有限公司制造)测定稻谷堆表层的孔隙率ε₀，稻谷堆表层的孔隙比为经计算，含水率为11.90％w.b.稻谷堆的孔隙比e₀为1.555。

稻谷堆表层的密度ρ₀采用如下方法测量：在稻谷堆表层内取样，测定样品的质量m₀和体积V₀，根据式(3)ρ₀＝m₀/V₀计算得到ρ₀。稻谷堆表层的密度ρ₀为585.00kg/m³。

含水率是指粮食中所含水的质量与粮食总质量的比值，用百分数表示。根据国标GB/T5497-85，含水率采用105℃恒重法测定，并按下式计算稻谷的含水率：w＝(m₁-m₂)/(m₁-m₀)×100％，其中：m₀为铝盒质量，g；m₁为烘前稻谷样品和铝盒质量，g；m₂为烘后稻谷样品和铝盒质量，g；w为含水率，w.b.。测定结果：稻谷堆的含水率MC为11.90％w.b.。

(3)稻谷的修正剑桥模型参数M、κ、λ和υ通过SLB-6A型应变控制式三轴仪(南京土壤仪器厂有限公司制造)进行三轴轴向压缩试验和各向等压压缩试验后计算得到。

采用三轴轴向压缩试验确定临界状态应力比M，包括如下步骤：装入稻谷样品，分别设定围压为30、50、70、90、110kPa，启动SLB-6A型应变控制式三轴仪施加轴向力对样品进行剪切，稻谷样品轴向位移每增加0.4mm，记录一次测力计的应力值读数和样品体积减少量，直至测力计读数出现峰值时，记下测力计读数峰值q及对应的平均压应力p，以p为横坐标q为纵坐标作图，经一元线性回归得直线的斜率，该斜率即为临界状态应力比M；

采用各向等压试验确定对数硬化模量λ和等向膨胀指数κ，包括如下步骤：装入稻谷样品，使围压σ₃从0kPa按照每次增加5kPa依次增加至200kPa，记录加载过程中围压σ₃和对应的样品体积减少量；再按照每次减小5kPa将围压由200kPa依次卸载至0kPa，记录卸载过程中围压σ₃和对应的样品体积增加量；绘制加载曲线和卸载曲线；经一元线性回归，将加载曲线的斜率作为λ，将卸载曲线的斜率作为κ；加载曲线和卸载曲线的横坐标为lnp,纵坐标为e，e是稻谷堆孔隙比，按下式计算：e＝V/[V₀(1-ε₀)]-1；式中，V₀是稻谷样品初始体积，m³；V是稻谷样品受压体积，m³；ε₀是稻谷堆表层孔隙率；p是平均压应力,kPa。

采用三轴轴向压缩试验测量弹性模量E和泊松比υ，包括如下步骤：装入稻谷样品，设置围压，轴向压缩样品直至破坏，测定最大主应力差(破坏时轴向压应力与围压差值)；将破坏的稻谷样品取出，再重新装入新的稻谷样品，进行加压实验，分4级施加轴向压力，第1级压力为所述最大主应力差的1/10，以后每级压力比前一级压力增加所述最大主应力差的1/10；每施加1级压力，记录施加的轴向压力和每次加压后1分钟时样品的轴向位移，每隔1分钟施加一级压力，施加到第4级压力为止；逐级卸压至施加的轴向压力全部卸去，分4级卸压，卸压顺序与加压顺序相反；重复加压卸压4遍后，轴向加压直至破坏；制作最后一次加压卸压实验中的加载曲线和卸载曲线，计算弹性模量E，E＝ΔP/(Δh_e/h_c)；其中，ΔP为轴向压力的递增或递减值，kPa；Δh_e为稻谷的轴向弹性变形量，mm；h_c为稻谷样品装样后的高度，mm；采用三轴各向等压试验确定稻谷样品体变弹性模量B，按公式B＝Δσ/(Δν/V)计算，Δσ是稻谷样品围压增量，kPa；Δv是稻谷样品体积压缩量，m³；V为稻谷样品装样后的体积，m³；泊松比ν按公式ν＝(3B-E)/6B计算。

结果见表2。

表2稻谷的修正剑桥模型参数

(4)将筒仓中未压实状态下(h′为7.7m)的稻谷堆划分成N层(N＝11)；每层稻谷堆由中心向外划分成中心轴线与筒仓中心轴线共线的M个单元(M＝5)，所述M个单元由1个圆柱和M-1个圆环柱组成，其中圆柱处于每层的中心(见图1)；每层稻谷堆的各单元按照由中心向外依次编号其列数；共得到N×M个单元；各单元的应变与应力关系符合修正剑桥模型；应用大型有限元软件ABAQUS求解修正剑桥模型，计算压实状态下稻谷堆各单元的高度h_nm(第n层m列单元的高度)，根据式(A)分别计算筒仓内压实状态下稻谷堆的N×M个单元体积应变

式(A)，式(A)中，n为层数(是自然数)，m为列数(是自然数)，是第n层、m列单元的体积应变，和是第n层、m列单元的三个主应变；

根据式(B)分别计算N×M个单元密度的ρ_nm：(n＝1,2,…,N；m＝1,2,…,M)式(B)，

式(B)中，n为层数(是自然数)，m为列数(是自然数)，ρ_nm是第n层、m列单元的密度，kg/m³；ρ₀是稻谷堆表层的密度，kg/m³；是第n层、m列单元的体积应变；

每层稻谷堆的各单元按照由中心向外依次编号其列数，具体编号方法如下：中心处的圆柱为第1列，然后依次向外编号为第2列，直至第M列。将每层稻谷堆划分成5个单元，见图1。稻谷堆的压实状态是指稻谷堆装入筒仓后的堆积状态。稻谷堆的未压实状态是指重力未施加时稻谷堆的堆积状态。稻谷堆表层是指稻谷堆表面至0.5米深度处的稻谷层。稻谷堆表层所受重力可以忽略。

应用大型有限元软件ABAQUS求解修正剑桥模型，包括如下步骤：在ABAQUS的Part模块中，采用轴对称CAX4单元，创建筒仓和稻谷部件；在Property模块中，定义筒仓的钢材材料属性和稻谷堆的属性，选择修正剑桥模型表征稻谷堆的应变和应力关系，然后将各材料属性指派到各个部件中；将筒仓和稻谷堆两个部件在Assembly模块中装配成整体；在Step模块中创建分析步；在Interaction模块中定义两个部件之间的约束关系，通过摩擦系数来描述两表面间的摩擦行为；在Load模块中定义荷载以及边界条件；在Mesh模块中按照部件对结构(稻谷堆与筒仓)进行网格划分，将筒仓中稻谷堆划分成所述N×M个单元；在Job模块中提交作业进行分析；分析结束后，在Visualization模块中查看压实状态下稻谷堆各单元的三个主应变及各单元的高度h_nm，根据计算第n层、m列单元的体积应变式中，是第n层、m列单元的三个主应变；h_nm是指压实状态下第n层、m列单元的高度。

创建筒仓和稻谷部件具体是指定义筒仓内径、仓高、仓壁厚度、粮堆未压实高度，筒仓部件是矩形平面，该矩形平面是筒仓外径剖面被中心轴线分割成的一半；稻谷堆部件是矩形平面，该矩形平面是圆柱形稻谷堆直径剖面被中心轴线分割成的一半。定义筒仓的钢材材料属性是指定义钢材材料的弹性模量与泊松比，定义稻谷堆的属性是指定义稻谷堆的修正剑桥模型参数。在Interaction模块中定义两个部件之间的约束关系为摩擦接触关系，通过稻谷与筒仓仓壁之间的摩擦系数μ来描述两表面间的摩擦行为。在Load模块中定义荷载是体积力，边界条件是：与筒仓底接触的稻谷无垂直位移，圆柱形稻谷堆中心轴线上的稻谷无水平位移。在Load模块中荷载通过采用体积力来描述，其他作用力通过稻谷堆与仓壁、仓底之间的相互作用来描述。在Interaction模块中定义稻谷堆与仓壁之间的相互作用通过设置接触来模拟，接触对选择点对面离散方法，主控面选定仓壁内侧，从属面选定粮堆模型外侧面。接触面相互作用力学模型选用常用的摩擦模型，通过摩擦系数来描述两表面间的摩擦行为。

根据式(C)计算每层稻谷堆的重量，

式(C)中，W_n是第n层稻谷堆的重量，kg；；ρ_nm是第n层、m列单元的密度，kg/m³；h_nm是压实状态下稻谷堆第n层、m列单元的高度；R_m表示第m列单元的外径(当m＝1时，R₁是指圆柱的直径)；n为层数(是自然数)，m为列数(是自然数)，n＝1,2,…,N；m＝1,2,…,M。

稻谷堆每层的重量，见表3。

表3筒仓中稻谷堆在每层的重量

(5)采用式(D)计算出筒仓内稻谷堆的总重量W：

式(D)中W_n是第n层稻谷堆的重量，n＝1,2,…,N。

筒仓内稻谷堆的总重量为121786.25kg。

(6)将计算结果同体积密度法计算值以及实际账面数进行比较，其中，体积密度法的修正系数为1.015，对比结果见表4.

表4计算值与实际账面数误差比较

结果显示，本发明方法测量方法得到的稻谷储藏总重量比常用的体积密度法误差小，且更接近实际账面数，体现出了本计算方法的优势。

实施例2采用本发明方法测定2号筒仓中稻谷的储藏重量

(1)中央储备粮温州直属库2号筒仓的内径d、仓高H见表5。该筒仓内空腔为圆柱体状。筒仓内压实状态下稻谷堆的高度为24.61m(测定值)；筒仓内未压实状态下稻谷堆的高度h′，根据当20m<h≤25m时，1.17<h′/h≤1.19，h′取值为29.0m。按照实施例1中方法，得到：稻谷与筒仓仓壁之间的摩擦系数μ为0.3，稻谷堆的内摩擦角为28°，稻谷堆的表层密度ρ₀为629.68kg/m³以及稻谷堆的含水率为13.50％w.b.。

表5 2号筒仓及其材料的参数

(2)按照实施例1中方法测定稻谷堆表层的孔隙率ε₀，并根据实验数据计算得到稻谷堆的孔隙比e₀，经计算，含水率为13.50％w.b.稻谷堆的孔隙比为1.527。

(3)对稻谷堆进行三轴轴向压缩和各向等压压缩试验，并根据实验数据计算得到稻谷堆的修正剑桥模型参数M、κ、λ与υ，具体方法同实施例1。结果见表6；

表6稻谷的修正剑桥模型参数

(4)将未压实状态下的稻谷堆划分成N层(N＝30)，每层划分成M个单元(M＝5)，各单元的应变与应力关系符修正剑桥模型，应用大型有限元软件ABAQUS求解修正剑桥模型，计算压实状态下各单元的高度h_nm(第n层m列单元的高度)及稻谷堆在筒仓中的空间密度分布值(即筒仓中稻谷堆150个单元的密度)，计算稻谷堆每层的重量，具体方法同实施例1，结果见表7。

表7筒仓中稻谷堆在每层的重量

(5)将每层的稻谷堆的重量求和得出筒仓内稻谷堆的总重量，总重量为7813825.08kg。

(6)将计算结果同体积密度法计算值以及实际账面数进行比较，其中，体积密度法的修正系数为1.05，对比结果见表8。

表8计算值与实际账面数误差比较

结果显示，本文计算方法得到的稻谷储藏总重量比常用的体积密度法误差小，且更接近实际账面数，体现出了本发明测量方法的优势。