陳 雪 程緒鐸 胡美珠

(南京財經大學食品科學與工程學院；江蘇省現代糧食流通與安全協(xié)同創(chuàng)新中心，南京 210046)

當小麥儲藏在平房倉中時，糧堆受到自身質量、內摩擦力以及倉壁與糧堆之間的摩擦力作用，產生壓應力和切應力[1]。小麥的籽粒松軟，糧堆各處的應力作用導致糧堆產生形變和體積壓縮，體積壓縮導致糧堆密度的增加，平房倉內小麥堆的應力分布是不均勻的，所以小麥堆密度分布也是不均勻的。每年我國均需投入大量人力、物力和資金開展清倉查庫工作[2]。糧食儲藏數量檢查是糧食庫存檢查的一項重要內容[3]，即通過一定的手段獲取糧倉中所儲存糧食的數量，主要包括稱重法和體積密度法[4]。稱重法工作量大、效率低，在大規(guī)模的庫存檢查中很難廣泛使用。體積密度法是指通過測量糧堆的體積和平均密度計算糧倉中糧食數量的檢查方法?，F有的體積密度法是使用糧食的表層密度乘以修正系數來獲取糧堆的平均密度，修正系數是憑經驗給出的，此方法誤差較大。因此，只有準確計算或測量出平房倉中的密度分布值才能準確快捷地計算平房倉中小麥堆的總質量。

前人通過在糧倉底部安裝了壓力傳感器[5,6]，測定糧倉底部的壓力分布，從而估算出糧堆的總質量，此方法成本高且只能得出糧堆總質量而無法得出糧倉內糧堆密度分布值。程緒鐸等[7]開發(fā)了一個分層壓縮模型來預測平房倉中的糧堆密度和壓力分布，但其假設平房倉中同一糧層的密度是相同的，而實際中平房倉內同一糧層的堆密度分布是不均勻的，該分層壓縮模型不能準確給出平房倉中的糧堆的密度分布值。

自20世紀70年代開始，有限元法被用來研究糧倉中物料的應力分布及倉壁的壓力問題。Jofriet[8]采用線彈性模型，用有限元方法研究了直筒倉內物料靜態(tài)時的應力分布問題。Ooi[9]使用有限元方法研究了柔性筒倉內散體的應力分布。國際上一些研究人員構建了許多本構方程，例如：Mohr-Coulomb模型、Drucker-Prager模型和Lade-Duncan模型等，通過有限元方法研究了糧倉中物料的應力分布問題。這些研究很好地計算出了儲藏在筒倉中壓縮性很小的沙子、碎石、工業(yè)物料的應力分布，這些彈塑性本構方程也可以較準確地求解出糧倉內糧堆的應力分布。糧食與土壤、沙子等物料的力學特性有很大的差異，小麥籽粒松軟，籽粒間的孔隙大，小麥堆在受壓時體積變化大且體積變化同時發(fā)生在彈性變形和塑性變形兩個階段。但研究所使用的本構方程均假定物料的體積變化很小且只發(fā)生在彈性變形階段，形狀變化發(fā)生在塑性變形階段。因此這些模型難以反映小麥堆變形時的應力與應變關系，用它們求解平房倉內小麥堆的應變分布是不準確的。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)可知，修正劍橋模型的參數為臨界狀態(tài)應力比M、等向膨脹指數κ、對數硬化模量λ、彈性模量E、泊松比υ、初始孔隙比e0。

1.4 小麥的修正劍橋模型參數測定

使用 HGT-1000A 容重器測定小麥堆無壓縮密度(糧倉表層小麥密度)ρ0；使用 LKY-1 型糧食孔隙率測定儀測定小麥堆的修正劍橋模型參數初始孔隙比e0；使用 TSZ-6A 應變控制式三軸儀測定小麥的修正劍橋模型參數E、υ、M、κ、λ。

1.5 平房倉倉壁的幾何與力學參數

混凝土倉壁可視為剛性體，在有限元分析過程中不產生變形。平房倉尺寸參考南京鐵心橋國家糧食儲備庫中糧倉的實際大小，具體倉壁材料的幾何與力學參數見表1。

表1 平房倉及其材料的參數

1.6 平房倉中小麥堆有限元模型構建與運行步驟

在Part步驟中采用三維可形變模型(3D Deformable)建模，創(chuàng)建平房倉和小麥堆部件，尺寸見表1；在Property步驟中分別設定倉壁的材料屬性和糧堆的塑彈性屬性(修正劍橋模型選擇塑性模型Clay Plasticity和彈性模型Porous Elastic)；在Assembly 步驟中將平房倉倉壁和糧堆部件組合起來；在 Step步驟中設定分析步； Interaction步驟是定義各個部件之間的相互接觸關系。糧堆與倉壁之間存在相互作用，在ABAQUS中設定一個完整的接觸必須包括兩個部分：一是接觸對，選擇點對面離散方法，選擇倉壁內側面為主接觸面，糧堆外側面為從屬接觸面；二是定義相互作用力，兩個接觸面之間會產生摩擦力。其中，庫侖模型是較為常用的摩擦模型，即通過摩擦系數來描述兩表面間的摩擦關系。

在Load步驟進行載荷以及邊界條件的設定，約束平房倉倉底和倉壁的位移。由于平房倉中的糧堆受到自重，通過體積力來加載。

在 Mesh 模塊中對模型進行網格劃分。小麥堆模型共被劃分為9層，將每層糧堆部件均等分為1 m2的單元(圖1)。本實驗采用節(jié)點法，提取各個節(jié)點應力應變值進行計算；在Job步驟中對構建的模型進行運算；待運算成功后，在Visualization步驟中查看運算結果提取所需數據。

建立有限元模型，平房倉內小麥堆長26 m，寬24 m，高9 m(重力未施加，糧堆無壓縮狀態(tài))。

1.7 小麥堆節(jié)點密度的計算

有限元軟件計算后提取小麥堆中各個節(jié)點的應變值，小麥堆節(jié)點密度可以由式(8)計算得到。

(8)

1.8 平房倉中小麥堆的糧層平均堆密度

限元軟件計算后提取小麥堆中各個節(jié)點的應變值，小麥堆的第i層平均堆密度可由式(9)計算。

(9)

式中：ρi為第i層小麥層的平均堆密度/kg/m3；Vijk為ijk個節(jié)點所在處的單元體積/m3；m是沿寬度方向的單元數，n是沿長度方向的單元數。

2 結果與分析

2.1 小麥的修正劍橋模型參數及無壓縮密度

實驗測定的修正劍橋模型的參數及無壓縮密度見表2。

表2 修正劍橋模型參數表及無壓縮密度

2.2 平房倉中小麥堆的密度分布

2.2.1 同一糧層小麥堆密度分布

平房倉中的小麥堆是一個軸對稱的長方體，選擇其中的四分之一小麥堆來分析其堆密度。由式(8)計算出平房倉內小麥堆的各節(jié)點的堆密度，本文給出了平房倉中小麥堆(含水率為12.03% w.b.)的第2、5、8層的堆密度分布，見圖1(設平房倉的拐角為原點)。

由圖1可知，平房倉內在同一糧層的小麥堆的堆密度分布是不均勻的。在同一糧層中，隨著距倉壁越遠，小麥的堆密度越大，并且在該層的中間趨于穩(wěn)定值。在同一糧層中，在平房倉的拐角附近的小麥的堆密度迅速下降，拐角處的小麥堆密度達到該層的最低值。不同糧層之間的小麥堆密度也存在明顯差異，小麥糧層深度越深，其拐角的堆密度變化越大，該糧層密度分布越不均勻。

圖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 裝糧方式對糧堆密度的影響

不同的裝糧方式糧食下落的速度不同，下落的糧食對糧堆產生的沖擊力不同，糧食落點處的密度與沖擊力的大小有關。對于小型糧倉或裝糧容器，由于糧堆高度不高，裝糧方式對裝糧后的糧堆密度有一定的影響。對于高大的平房倉，盡管糧食沖擊力改變了糧食落點密度(相對于糧食下落速度為零的糧堆落點密度)，后續(xù)下落的糧食壓住了先前的糧食落點，隨著糧食的不斷下落，落點上方的壓力越來越大，落點上方的壓力改變糧食落點密度，這種改變取代了沖擊力對糧食落點密度的改變量，即沖擊力對糧食落點密度的改變量有壓力引起的密度改變量代替了。對于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%之間。