筒倉糧堆內部豎向壓力計算方法

李東橋 韓 陽 陳家豪 靜 行 段君峰

(河南工業大學土木建筑學院，鄭州 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 基于倉壁有效摩擦系數分布規律，提出了公式的簡化算法，并進行了實驗。通過與實驗結果和有限元結果對比，驗證了本文簡化算法的合理性。簡化的計算方法較為簡便，可為糧倉結構設計中儲料壓力荷載的確定和倉儲方向尚未解決的難題提供理論支持。