筒仓粮堆内部竖向压力计算方法

李东桥 韩 阳 陈家豪 静 行 段君峰

(河南工业大学土木建筑学院，郑州 450001)

粮食等散体物料具有复杂的力学特性。不同于固体，散粒体不能承受或只能承受较小的拉力，但可以承受较大的压力和剪力；也不同于液体，散粒体虽可以向各方向传递压强，但并不相等。由此造成粮堆内部存在如边坡稳定性、筒仓结拱、筒仓超压等诸多问题[1-3]。在处理这些问题时，国内外学者大多把研究集中在在如何确定粮堆边界上的压力[4-6]，较少粮堆内部的空间压力分布情况。然而，要解决这些问题，仅考虑粮堆边界上的压力是远远不够的，因为这些问题或者其本质就是粮堆的空间压力问题，或者是以粮堆的空间压力为基础的问题。目前，鲜有学者对粮堆内部任意一点的压力进行研究。如果把对粮堆压力的研究从粮堆的边界转向“空间”，给出粮堆的空间压力场模型，大有可能解决很多仓储方向尚未解决的问题。

由实仓测试[7-9]和模型仓实验[10]的结果可知，粮堆内同一水平面上及仓底处的竖向压力都是不均匀分布的，并且其不均匀分布情况随着粮堆深度、高径比的增加逐渐明显。同时，在静态储粮下，粮食与仓壁间的实际摩擦力尚未达到最大静摩擦力，实际摩擦力与侧压力之比随着储粮深度的变化而变化。

Janssen公式[11]及我国规范都是在“仓内粮堆竖向压力在水平面上均匀分布”的假定下得到的，不能计算粮堆内部的空间压力。同时，储藏状态下储料与仓壁之间并未达到极限平衡状态，此时假定“外摩擦系数沿深度方向不变”并且以外摩擦系数计算仓壁摩擦力是不合适的。并且，粮堆内部储料之间的相互作用同样会影响同一深度处的粮堆压力分布情况。

为准确计算粮堆空间竖向压力，便不能沿用Janssen公式“竖向压力在水平面上均匀分布”和“外摩擦系数不变”这两条假定，同时需要确定仓壁摩擦力的分布规律及储料间的相互作用机理。以此，为粮仓结构设计提供准确的储料压力荷载，为仓储研究难题提供参考。

1 筒仓有效摩擦系数的分布规律

筒仓模型实验表明，随着高径比的增加，粮食对仓壁的摩擦逐渐显著，当高径比达到1.0～2.0时，仓壁摩擦力占粮食总质量的30%～48%。同时，粮堆内部还存在储料间的相互摩擦等复杂行为。本研究认为仓壁摩擦力及储料间的相互作用共同造成了粮堆同一水平面上竖向压力的不均匀分布现象，掌握了它们的分布规律，便能解决粮堆竖向压力的计算问题。

1.1 储料与仓壁间的有效摩擦系数

由于静态储粮时储料与仓壁之间并未达到极限平衡状态，导致实际摩擦力并不等于最大外摩擦系数与侧压力的乘积，并且外摩擦系数的值随着粮堆深度的变化而变化。为更准确研究摩擦力沿仓壁的分布情况，本研究引入“仓壁的有效摩擦系数”μwg这一概念，其等同于同一深度处仓壁摩擦力τwg与仓壁侧压力σr的比值，即μwg=τwg/σr。不同于外摩擦系数μ0，μwg并非定值，而是沿粮堆深度z变化的函数，其具体数学表达式μwg(z)可通过有限元算例获得。

文献[10]提出，仓壁的摩擦力分布与仓径大小、装粮高度和外摩擦系数等因素有关。以下三种情况分析了在不同仓径D、装粮高度H及外摩擦系数μ0影响下，仓壁有效摩擦系数μwg沿深度变化的有限元结果。

1.1.1 装粮高度H不变，筒仓不同直径D影响下仓壁有效摩擦系数μwg沿深度z的变化

算例1：筒仓直径分别为4、10、16、22 m，装粮高度均为6.35 m，高径比分别等于1.6、0.6、0.4和0.3。容重为804 kg/m3，内摩擦角为25°(内摩擦系数0.47)，储料与仓壁的摩擦系数为0.40。有限元结果如图1所示。

图1 筒仓不同直径D情况下仓壁有效摩擦系数与深度z的关系

算例2：筒仓直径分别为16、22 m和32 m，装粮高度均为20 m，高径比分别等于1.3、0.9和0.6。容重为804 kg/m3，内摩擦角为25°(内摩擦系数0.47)，储料与仓壁的摩擦系数为0.40。有限元结果如图2所示。

图2 筒仓不同直径D情况下仓壁有效摩擦系数与深度z的关系

1.1.2 筒仓直径D相同，不同装粮高度H影响下仓壁有效摩擦系数μwg沿深度z的变化

算例3：筒仓装粮高度分别为6.35、10、20、30 m和40 m，直径均为16 m，高径比分别等于0.4、0.6、1.3、1.9和2.5。容重为804 kg/m3，内摩擦角为25°(内摩擦系数0.47)，储料与仓壁的摩擦系数为0.40。有限元结果如图3所示。

图3 筒仓不同装粮高度H情况下仓壁有效摩擦系数与深度z的关系

1.1.3 筒仓直径D、装粮高度H相同，不同外摩擦系数μ0影响下仓壁有效摩擦系数μwg沿深度z的变化

算例4：外摩擦系数μ0分别为0.30、0.40和0.45，筒仓装粮高度H均为6.35 m，直径D均为16 m，筒仓高径比均等于0.4。容重为804 kg/m3，内摩擦角为25°(内摩擦系数0.47)。有限元结果如图4所示。

图4 不同外摩擦系数μ0情况下仓壁有效摩擦系数与深度z的关系

可以看到，在粮堆深度为0时，有限元计算得到的有效摩擦系数μwg并不为0，这是由于有限元的特点是先计算变形，即网格单元的变形(压缩、剪切变形)，深度为0处的数值是根据深度为0附近处单元计算得到的。若单元足够小，则深度0处仓壁侧压力、摩擦力接近真实值0，但它们的比值τwg/σr并不如此。为更准确分析仓壁有效摩擦系数沿深度变化的规律，应将摩擦系数初始值考虑在内。本研究参与了文献[10,12]的模型仓实验及有限元计算工作，图5为实验结果与有限元结果的对比，可以看到，除粮堆深度为0处实验结果与有限元结果存在一定偏差，其他深度处两者基本一致。

图5 实验结果与有限元结果对比

从图1至图4可以看出，筒仓直径D、装粮高度H、外摩擦系数μ0及深度比z/H等因素均对仓壁有效摩擦系数μwg有影响。

对于不同的仓型，仓壁有效摩擦系数的变化规律相近。粮堆深度较小时，μwg随深度的增大近似线性增长；粮堆深度到达某一点之后，μwg逐渐达到峰值，并趋于不变；在粮堆底部，略有减少。针对上述规律，可将仓壁有效摩擦系数与深度的关系图形作双折线图形处理，即一条线性斜线和一条铅锤线的组合。

在仓壁有效摩擦系数μwg到达峰值前，其与深度的关系均近似于线性函数μwg=kz+b的形式，并有以下规律：(1)装粮高度H不变，随着仓径D的增加，线性函数的斜率k变化不大，但是其与水平轴的截距b逐渐减小；(2)仓径D不变，随着装粮高度H的增加，线性函数的斜率k略微减小，与水平轴的截距b逐渐增大；(3)随着外摩擦系数μ0的增加，线性函数的斜率k、其与水平轴的截距b均增大。

到达峰值点后，仓壁有效摩擦系数μwg基本保持不变。但是，对于不同高径比的筒仓，峰值点不同：高径比越大，仓壁有效摩擦系数μwg越早到达峰值，即到达峰值点时深度比(深度与装粮高度的比值)z/H越小。

根据以上规律，分析μwg与各变量间的关系，可假定仓壁有效摩擦系数μwg的函数服从式(1)分布：

(1)

式中：zmax指有效摩擦系数达到峰值时的粮堆深度，其选取与筒仓的高径比有关，高径比越大，zmax/H越小；A和B为待拟合常数。

根据有限元及模型仓实验结果，假定：H/D≥1时，zmax=0.13H；1>H/D≥0.6时，zmax=0.25H；H/D<0.6时，zmax=0.3H。

使用Minitab软件基于最小二乘法进行数据拟合，得到A、B两个常数，其结果如表1所示。

表1 常数A、B拟合结果

(2)

将A、B代入μwg的分布函数，并与有限元结果对比，其结果如图6～图8所示。

算例1：

图6 不同仓径影响下本文公式与有限元结果比较

算例3：

图7 不同装粮高度H影响下本文公式与有限元结果比较

算例4：

图8 不同外摩擦系数μ0影响下本文公式与有限元结果比较

可以看到，本文公式与有限元结果比较接近，这样的双折线图形基本可以反映不同仓型在不同装粮高度及外摩擦系数等因素影响下的仓壁有效摩擦系数沿深度的分布规律。

1.2 粮堆内部有效摩擦系数的分布规律

粮堆内部存在储料间的相互摩擦、土拱效应、多场耦合等诸多复杂行为，研究粮堆竖向压力分布时将这些因素均考虑在内会使计算变得十分困难。如果把上述问题看作为粮堆内部储料间摩擦力的另一种表现形式，将同一深度处储料间摩擦力(剪切力)沿径向的不均匀分布当做影响储料竖向压力不均匀分布的主要因素，只考虑粮堆空间任意位置储料微元体所受实际摩擦力与水平压力的比值μgg的分布情况，竖向压力的计算会简单许多。

本研究引入“粮堆内部的有效摩擦系数”μgg这一概念，其等同于粮堆空间任意位置储料微元体所受实际摩擦力τgg与所受水平压力σr的比值，即μgg=τgg/σr。μgg同样不是常数，而是与粮堆深度z和测点到仓中心距离x有关的函数，其数学表达式μgg(x,z)可通过有限元算例获得。本研究将储料与仓壁间的有效摩擦系数μwg与粮堆内部有效摩擦系数μgg统称为筒仓的有效摩擦系数μ(x,z)。μ(x,z)是与测点深度z及测点位置到筒仓中线的水平距离x有关的函数，表示储料任意一点处实际摩擦力与竖向压力的比值。

本研究利用有限元结果分析了筒仓高径比、装粮高度、外摩擦系数、有效摩擦系数及测点至仓心距离等因素对粮堆内部有效摩擦系数μgg的影响，有限元算例参数见表2。

表2 算例参数

对比有限元结果，总结出以下结论：粮堆内部的有效摩擦系数μgg与筒仓直径、装粮高度、仓壁有效摩擦系数、深度以及该点距仓中心水平距离等因素有关，但主要是与装粮高度、仓壁有效摩擦系数和测点距仓中心水平距离均等因素有关。

同一深度位置，有效摩擦系数随着筒仓直径的增大而减小，随着装粮高度的增加而增大，随着外摩擦系数的增大而增大。

同一筒仓的相同深度处的有效摩擦系数随着测点距仓中心水平距离的增加而增大，有效摩擦系数的分布形式近似呈指数函数。

根据以上规律，陈家豪[9-10]提出可将粮堆内部摩擦力与水平压力的比值μgg的分布假定为式(3)函数。

(3)

式中，μwg为仓壁有效摩擦系数，x为测点到同一水平面的筒仓中心的水平距离，R为筒仓半径；H为装粮高度；C为待拟合常数，通过数值拟合，C=3.48。

当测点至筒仓中线的水平距离x=R时，式(3)变为μgg=μwg，即μgg变为仓壁处有效摩擦系数μwg；x=0时，μgg=0，即仓心处有效摩擦系数为0。

综上，粮堆内任意一点的有效摩擦系数可表示为：

(4)

H/2R≥1时，zmax=0.13H；1>H/2R≥0.6时，zmax=0.25H；H/2R<0.6时，zmax=0.3H。

式中：A=0.194 2，B=0.42，C=3.48。

2 粮堆内部空间竖向压力的计算方法

2.1 计算思想及平衡方程的确定

如果将粮堆横截面划分为若干个同心圆环，如图9a所示。沿纵向切割，粮堆被划分为一个实心圆柱和若干个空心圆柱单元，当划分单元足够小时，可以认为竖向压力在同一水平面上的同一圆环内是均匀分布的，第n个空心圆柱竖向受力情况如图9b所示。

图9 计算思想示意图

图9b中，2(x-t)为圆环内径,2x为圆环外径，t为圆环宽度，σz指第n个圆环储料所受竖向压力，G为圆环内储料重力，τ(x,z)指粮堆深度z，距筒仓中线水平距离x处储料所受摩擦力。

由实仓测试和模型仓实验的结果可知，同一水平面上的水平压力趋于均布，且随着粮堆深度的增加其不均匀分布的程度并不明显，由此，可假定粮堆同一深度的水平压力处处均等，且等于仓壁处侧压力。结合本研究所得有效摩擦系数的分布规律，在环宽t足够小时，可以准确计算粮堆空间任意一点出的竖向压力。

综合以上观点，本研究计算方法做出如下假设：

(1)粮堆同一深度的水平压力均匀分布，且等于仓壁处侧压力；

(2)粮堆内部有效摩擦系数服从一定分布规律，且粮堆深度到达某一数值后，仓壁有效摩擦系数接近为定值；

(3)粮食为均匀介质，容重为定值。

则图7(b)中竖向平衡方程为：

γπ[x2-(x-t)2]dz+σzπ[x2-(x-t)2]+τ(x-t,z)2π(x-t)dz=(σz+dσz)π[x2-(x-t)2]+τ(x,z)·2πxdz

(5)

式(5)中，x2-(x-t)2=2xt-t2，若环宽t极小，为简化平衡方程，可设t2=0，则原方程变为：

γxtdz+τ(x-t,z)xdz-τ(x-t,z)tdz=xtdσz+τ(x,z)·xdz

(6)

化简得：

(7)

式(7)中，τ(x,z)可进一步表示为：

τ(x,z)=σr(z)·μ(x,z)

(8)

式中：σr(z)表示粮堆深度为z处的储料所受水平压力，其仅与深度z有关，与测点到筒仓中线的水平距离无关。有效摩擦系数μ(x,z)不同于外摩擦系数μ0，是实际摩擦力与水平压力比值的函数，故式(8)中τ(x,z)为任意一点处实际摩擦力。

将式(8)代入式(7)，得：

(9)

等号两边同时沿z方向积分得:

(10)

对于水平压力σr(z)，我国规范GB J77—85上筒仓的水平压力计算方法为：

(11)

式中：Cr为水平压力修正系数；z为储料顶面或储料锥体重心至所计算截面的距离；K为侧压力系数；μ0为外摩擦系数。结合本研究假设，σr(z)在同一水平面上均布。

将式(4)、式(11)带入式(10)：

z≤zmax时，解之得：

(12)

当t趋于0时，有：

(13)

其中，

(14)

将式(14)带入式(12)，得：

(15)

带入边界条件：σzz=0=0，即：

(16)

将式(16)带入式(15)得：

(17)

z>zmax时，式(10)变为：

(18)

解之得：

(19)

带入边界条件：σzz=zmax=σz(zmax)，得：

C2=σz(zmax)-γzmax

(20)

将式(19)带入式(18)，得：

(21)

其中，Δz=z-zmax。

综上，σz的表达式为：

(22)

H/2R≥1时，zmax=0.13H；1>H/2R≥0.6时，zmax=0.25H；H/2R<0.6时，zmax=0.3H。

式(22)即为粮堆内部任意一点竖向压力的计算方法。

2.2 公式结果与验证

算例10：筒仓直径为16 m，装粮高度为6.35 m，储料容重为804 kg/m3，内摩擦角25°，储料与仓壁间摩擦系数0.40，高径比为0.4，属浅仓。本文公式结果与有限元结果对比如图10所示。

图10 算例10不同深度处竖向压力

算例11：筒仓直径为16 m，装粮高度为10 m，储料容重为804 kg/m3，内摩擦角25°，储料与仓壁间摩擦系数0.40。高径比0.6，属浅仓。本文公式结果与有限元结果对比如图11所示。

图11 算例11不同深度处竖向压力

算例12：筒仓直径为16 m，装粮高度为20 m，粮堆密度为804 kg/m3，内摩擦角25°，储料与仓壁间摩擦系数0.40。高径比1.3，属浅仓。本文公式结果与有限元结果对比如图12所示。

本研究结果与有限元结果比较接近，基本反映了粮堆空间竖向压力的分布规律。

可以看到，粮堆竖向压力的不均匀程度随着深度的增加而增大。当粮堆深度较小时，竖向压力较为均匀；当粮堆深度较大时，同一深度竖向压力呈现了“中间大，两边小”的分布特征，且其不均匀分布情况随着粮堆深度的增加越发明显，边界处与中心处竖向压力差值随粮堆深度的增加而增大。

图12 算例12不同深度处竖向压力

2.3 本文公式的化简

本研究将仓壁有效摩擦系数μwg当做一关于深度z的函数进行推导，使式(22)结果相对复杂，不利于实际计算。但是，若把μwg当做一常数处理，计算过程会简便很多。由此，本研究提出了粮堆空间竖向压力的简化算法。

有限元结果显示，当粮堆深度到达某一数值zmax时，仓壁有效摩擦系数μwg趋于稳定，随深度的增加基本不变。由于在粮堆上部，同一水平面上水平压力σr较小，有效摩擦系数μwg的变化对水平压力σr与μwg的乘积影响不大，这里可将仓壁有效摩擦系数μwg假设为μ′，使计算简化。μ′为定值，表示粮堆深度到达zmax之后的仓壁有效摩擦系数。

依据假设，式(12)变为：

(23)

解之得：

(24)

由于储料与仓壁间并未达到极限平衡状态，仓壁有效摩擦系数应小于外摩擦系数。而其峰值μ′的选取与外摩擦系数、高径比有关。沿用“1.1”中所得仓壁有效摩擦系数的分布规律：在筒仓半径大于1 m时：μ′=Aμ0·zmax+Bμ0·e1/R，H/2R≥1时，zmax=0.13H；1>H/2R≥0.6时，zmax=0.25H；H/2R<0.6时，zmax=0.3H。

2.4 实验及简化公式的验证

算例10见图13。

图13 算例10不同深度处竖向压力

算例11见图14。算例12见图15。

模型仓实验：模型仓仓壁由有机玻璃制成，直径为0.5 m，装粮高度为1 m，高径比2，属深仓。仓底板分为一个中心圆和三个同心圆环，中心圆与同心圆环之间、各同心圆环之间及最外围同心圆环与仓壁之间均设有小于1 mm的缝隙，以防止粮食流出，如图16a所示。中心圆和每个圆环的底部都各有一组力传感器，能够测量粮堆底部不同区域的竖向压力，如图16b所示。由于最大装粮高度仅为1 m，粮堆内各压力(竖向压力、水平压力和摩擦力)都不大，因此在讨论某一深度处的粮食压力时，可以忽略边界条件(该深度处是刚性底板或粮面)对压力的影响，可近似将某一装粮高度下的仓底竖向压力表述为粮面以下某一深度处水平面上的竖向压力。

图14 算例11不同深度处竖向压力

图15 算例12不同深度处竖向压力

图16 模型仓示意图

储料容重为804 kg/m3，内摩擦角25°，储料与仓壁间摩擦系数0.40。由于模型仓半径小于1 m，根据实验数据：μ′=0.36。

实验结果与本文简化公式结果对比如图17所示：

图17 模型仓不同深度处竖向压力

整体上看，简化公式与有限元结果及实验数据都有较好的一致性。

3 结论

实验及有限元结果表明，粮堆同一水平面上竖向压力是不均匀分布的，其沿深度方向变化规律如下：深度较小时，竖向压力分布较为均匀，平均竖向压力随深度的增加近似线性增长；随着粮堆深度的增大，同一水平面上竖向压力逐渐呈现出“中间大，两边小”的不均匀分布特征，其不均匀分布程度随着粮堆深度的增加越发明显，同一深度仓壁出竖向压力与粮堆中心处竖向压力的差值逐渐增大。基于数值模拟和实验数据，本研究进行了筒仓空间竖向压力的计算方法的推导，并提出以下结论：

3.1 本研究提出筒仓内部的摩擦力(仓壁摩擦力和粮堆内部竖向切应力)沿仓径方向和深度方向的变化是造成粮堆同一深度竖向压力不均匀分布的主要因素。为准确分析实际摩擦力的分布情况，本研究引入仓壁有效摩擦系数和粮堆内部的有效摩擦系数这两个概念，分别代表同一深度处仓壁摩擦力和侧压力的比值及粮堆内部竖向切应力与水平压力的比值，将两者统称为筒仓的有效摩擦系数。基于数值模拟结果和实验数据，讨论了在筒仓直径、装粮高度、外摩擦系数、粮堆深度、深度比等因素影响下，有效摩擦系数在仓壁处和粮堆内部的分布规律，并以此提出了筒仓有效摩擦系数的数学表达式。

3.2 根据筒仓有效摩擦系数的分布规律，本研究将粮堆沿纵向分割为n个同心圆环，假定每个同心圆环上竖向压力均匀分布，根据平衡方程推导出粮堆空间任意一点处竖向压力的计算方法。该计算方法与数值模拟结果具有较好的一致性。

3.3 基于仓壁有效摩擦系数分布规律，提出了公式的简化算法，并进行了实验。通过与实验结果和有限元结果对比，验证了本文简化算法的合理性。简化的计算方法较为简便，可为粮仓结构设计中储料压力荷载的确定和仓储方向尚未解决的难题提供理论支持。