平房倉中小麥堆的孔隙率分布值研究

邵小龍，程緒鐸

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

我國自1998年以來新建的儲備糧庫及地方糧庫的倉型一般為淺圓倉、立筒倉和高大平房倉，除港口和交通樞紐中轉庫以淺圓倉和立筒倉為主外，其余多為高大平房倉，其儲存量約占總儲存量的70%[1]。

小麥儲藏在平房倉中，受到自重、內摩擦力以及倉壁的支持力，小麥堆各處產生壓應力與切應力，從而產生體積應變與形狀應變，體積縮小，孔隙率減小。平房倉中的小麥堆的應力與應變分布是非均勻的，因此孔隙率的分布也是非均勻的?？紫堵实拇笮∨c分布決定了糧食通風的阻力與均勻性，還影響糧堆中的熱量與水分的遷移速率[2-4]。

糧食孔隙率的測定始于20世紀60年代，Thompson等[5]使用氣體比重瓶測定了燕麥、小麥、高粱、大豆、黑麥、大麥、玉米和紫花苜蓿種子的孔隙率。之后，Gustafson等[6]和Chang[7]也利用氣體比重瓶測定了玉米、小麥和高粱的孔隙率。郝倩等[8]開發(fā)了一種糧食孔隙度就倉測量系統(tǒng)。還有學者提出了光學法[9-10]、聲學法[11]等。但光線與聲波無法到達倉內糧堆深處，不適合用來測定糧倉內糧堆孔隙率。

到目前為止，還未見平房倉中糧堆的孔隙率分布的理論與實驗成果報道。本文選定修正劍橋模型為小麥堆的應力與應變關系本構模型，使用有限元方法求解平房倉中小麥堆修正劍橋模型，計算出平房倉內小麥堆的應變分布值，由應變分布值與表層孔隙率（無壓縮）計算出平房倉中小麥堆的孔隙率分布值。

1 材料與方法

1.1 實驗材料

實驗小麥品種名為淮麥44號，產地江蘇南京，手工剔除破碎與不成熟顆粒，初始濕基含水率為12.03%，小麥的最大粒徑、中等粒徑和最小粒徑分別為5.88、3.07、2.79 mm（隨機取樣100顆小麥粒，用游標卡尺測得其3個粒徑，精確到0.01 mm，結果取其平均值）。實驗前，取20 kg小麥樣品，放置在太陽光下晾曬72 h，測得其濕基含水率為9.72%；再取20 kg小麥樣品，分成兩份，各自加適當的蒸餾水，然后將其裝入塑料袋中封閉好，置于恒溫箱中，溫度設定5 ℃，1周后從恒溫箱中取出，測得兩份小麥樣品濕基含水率分別為14.61%和16.55%。用標準烘箱干燥法將10 g樣品在130 ℃下干燥19 h，重復3次，測定小麥的含水率。實驗用的小麥的濕基含水率分別為：9.72%、12.03%、14.61%、16.55%。

1.2 實驗儀器

TSZ-6A型應變控制式三軸儀：南京土壤儀器廠有限公司；LKY-1型糧食孔隙率測定儀：南京土壤儀器廠有限公司。

1.3 小麥堆的應力與應變關系本構模型

平房倉內的小麥堆在重力、倉壁抵抗力的作用下，小麥堆內部各處產生了應力，由于應力的作用，小麥堆各處又產生了應變，小麥堆的應力與應變相互關聯(lián)，它們之間關系的微分方程就稱為應力與應變關系的本構模型。要準確地求解平房倉內小麥堆的應變分布必須選定恰當的應力與應變關系的本構模型。小麥籽粒松軟，小麥堆孔隙大，在應力的作用下產生較大的體積應變與形狀應變。修正劍橋模型中以體積應變和剪切應變作為應變分量，能表征小麥堆大體積壓縮的特性。所以本研究選定修正劍橋模型為小麥堆的應力與應變關系的本構模型。

修正劍橋模型的彈塑性矩陣的增量形式如下[12]：

式中：q為廣義剪切應力，Pa；p為平均主應力，Pa；εv為體積應變；εs為剪切應變；K為體積模量，kPa ；G為剪切模量kPa。

由式（1）～式（5）可知，修正劍橋模型中需要確定的參數為臨界狀態(tài)應力比M、等向膨脹指數κ、對數硬化模量λ、彈性模量E、泊松比υ、初始孔隙比e0。

1.4 修正劍橋模型參數的確定

1.4.1 初始孔隙比e0的測定

小麥堆孔隙率ε為小麥堆中的孔隙體積與整個小麥堆體積之比?？紫侗萫指小麥堆中孔隙體積與小麥堆籽粒體積之比。本實驗通過LKY-1型糧食孔隙率測定儀（見圖1）測定小麥堆的初始孔隙率，初始孔隙比可由初始孔隙率計算得到。

圖1 孔隙率測定儀示意圖

如圖1所示，容器A和B的容積相等，將小麥樣品裝入容器B中，裝滿后將容器B密封。將閥門2關閉，然后將閥門1和3打開。將一定壓力的氣體通過空氣壓縮機充滿容器A。待壓力表達到一定數值后，將閥門1關閉，等待壓力穩(wěn)定后，記錄下此時壓力表讀數（P1）；關閉閥門3，然后打開閥門2，待容器A和B的壓力達到平衡后，記錄壓力表讀數（P2）。設空氣是理想氣體，根據理想氣體等溫特性原理可以推導出小麥樣品的初始孔隙率（一個大氣壓下的孔隙率）為：

式中：Vs為樣品籽粒的體積，m3；V0為小麥樣品初始體積，m3。

1.4.2 臨界狀態(tài)應力比M的測定

臨界狀態(tài)應力比M可通過軸向壓縮實驗測定[12]。使用應變控制式三軸儀對小麥樣品進行不同圍壓（30、50、70、90、110 kPa）的剪切壓縮實驗，每當樣品位移增加0.4 mm，則記錄下樣品體積變化量和測力計讀數；當測力計讀數產生最大值時，即為最大主應力差，記錄下對應的p和q值。將記錄的p和q值擬合成臨界狀態(tài)應力直線，M的值即為直線的斜率，進行3次平行實驗。

1.4.3 對數硬化模量λ和等向膨脹指數κ的測定

通過應變控制式三軸儀對小麥樣品進行各向等壓壓縮，將圍壓σ3由0 kPa增加到200 kPa，每當σ3增加了5 kPa，則記錄一次樣品體積減少量。然后，再由200 kPa卸載到0 kPa，當σ3減小5 kPa，記錄一次體積變化量。最后，將p（圍壓）和孔隙比e繪制成加載曲線和卸載曲線，其直線斜率分別為λ和κ的值[12]，進行3次平行實驗。

1.4.4 彈性模量E和泊松比υ的測定

（1）彈性模量E的測定：參照SL237—1999[13]。選擇實驗測定的最大主應力差的10%作為施壓增量，對糧堆進行每級遞增加壓和遞減卸壓。加壓過程中，在對糧堆施加第一級壓力后，進行計時，對糧堆每隔1 min增加一級壓力，并且記錄下相應的壓縮位移，直至增加到第4級壓力；然后由第4級壓力逐級卸壓，同樣每隔1 min減去一級壓力，并記錄糧堆壓縮位移，直至施加的軸向壓力為0；反復加卸載4次之后，對糧堆進行持續(xù)壓縮直至糧堆被破壞，進行3次平行實驗。

由最后一次加卸載循環(huán)中得到的實驗數據計算出彈性模量。

式中：E為糧堆的彈性模量，kPa；ΔP為軸向荷載，kPa；Δhe為糧堆的彈性變形量，mm；hc為試樣固結后的高度，mm。

（2）泊松比υ的測定：糧堆進行等向壓縮實驗，由式（11）計算出體變彈性模量[14]。

式中：B為小麥堆的體變彈性模量，kPa；ΔV為小麥堆的彈性體積增量，m3；V為小麥堆裝樣后的體積，m3。

則泊松比υ為：

1.5 平房倉中小麥堆應力與應變的有限元計算方法

平房倉內儲藏小麥后，小麥堆各處都產生了應力與應變，小麥堆各處的孔隙率與初始孔隙率及各處的應變密切相關。只要測定出小麥堆的初始孔隙率并計算出小麥堆各處的應變，就能計算出小麥堆各處的孔隙率。小麥堆各處的應變可使用有限元方法求解小麥堆的修正劍橋模型而得出。

1.5.1 平房倉倉壁的幾何與力學參數的選定

混凝土倉壁可視為剛性體。平房倉尺寸參考南京鐵心橋國家糧食儲備庫糧倉的實際大小，具體倉壁材料的幾何與力學參數見表1。

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

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

（1）在 Part 步驟中采用三維可形變模型（3D Deformable）建模，創(chuàng)建平房倉和小麥堆部件，尺寸見表1。

（2）在Property步驟中分別設定倉壁的材料屬性和糧堆的塑彈性屬性（修正劍橋模型選擇塑性模型Clay Plasticity和彈性模型Porous Elastic）。

（3）在Assembly 步驟中將平房倉倉壁和糧堆部件組合起來。

（4）在 Step步驟中設定分析步。

（5）Interaction步驟是定義各個部件之間的相互接觸關系。糧堆與倉壁之間存在相互作用，在ABAQUS中設定一個完整的接觸，包括兩個部分：一是接觸對，選擇點對面離散方法，選擇倉壁內側面為主接觸面，糧堆外側面為從屬接觸面；二是定義相互作用力，兩個接觸面之間會產生摩擦力。其中，庫侖模型是較為常用的摩擦模型，即通過摩擦系數來描述兩表面間的摩擦關系。

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

（7）設小麥堆處于無壓縮狀態(tài)時堆長26 m，寬24 m，高9 m。在 Mesh 模塊中對模型進行網格劃分。小麥堆模型共被劃分為9層，將每層糧堆部件均等分為1 m2的單元。本文采用節(jié)點法，提取各個節(jié)點應力應變值進行計算。

（8）在 Job 步驟中對構建的模型進行運算。

（9）待運算成功后，在Visualization步驟中查看運算結果提取所需數據。

1.6 平房倉中小麥堆孔隙率的計算

1.6.1 平房倉中小麥堆單元孔隙率的計算

小麥堆經修正劍橋模型的有限元計算后提取各個單元節(jié)點的應變值，其單元孔隙率為：

式中：Vv0為小麥堆初始（表層）單元孔隙的體積，m3；V0為小麥堆初始（表層）單元體積，m3；ΔVijk為第ijk個單元的體積增量，m3。

式（13）分子分母同除以V0得到：

式中：ε0為小麥堆的初始（表層）孔隙率。

1.6.2 平房倉中小麥堆的糧層平均孔隙率

有限元軟件計算后提取小麥堆中各個單元的應變值，小麥堆的第i層平均孔隙率可由式（16）計算出：

式中：Vijk為第ijk個單元所在處的單元體積，m3；m為沿寬度方向的單元數；n為沿長度方向的單元數。

2 結果與分析

2.1 平房倉中小麥堆的孔隙率分布特性

由有限元模型計算后可以提取小麥堆中各個單元節(jié)點的應變值，而糧堆各單元的孔隙率可以根據式（15）計算得到。為了研究其同一糧層的小麥堆孔隙率分布特性，本文選取不同含水率小麥堆的第2、5、8層的孔隙率進行具體分析，結果如圖2～圖5所示（設平房倉的拐角為原點）。

由圖2～圖5可知，對于同一含水率，在相同糧層深度下的小麥堆的孔隙率分布是不均勻的。在同一層中，距倉壁距離增加，小麥的孔隙率減小，在該層的中心，在拐角處的孔隙率最大。不同糧層之間的小麥堆的堆積密度的分布也存在明顯差異，在平房倉中小麥堆的糧層深度越深，該層的小麥孔隙率分布越不均勻。

圖2 濕基含水率9.72%的小麥堆不同糧層的孔隙率分布值

圖3 濕基含水率12.03%的小麥堆不同糧層的孔隙率分布值

圖5 濕基含水率16.55%的小麥堆不同糧層的孔隙率分布值

在各個含水率中，同一糧層深度下的小麥堆的孔隙率分布是不均勻的。在同一糧層中，小麥堆的孔隙率由倉壁至糧層中間逐漸減小，在該層中間的孔隙率最小，在平房倉的拐角處孔隙率陡然增大并達到該層的最大值。不同層之間的孔隙率也存在明顯差異。糧層深度增加，同一水平面位置處的孔隙率減小；糧層深度增加，糧層各處的孔隙率的差值增大。同一糧層各處的孔隙率差值很小，濕基含水率9.72% 的小麥堆底層（第8層）的孔隙率最大值為42.69%，最小值為42.58%；濕基含水率12.03%的小麥堆底層（第8層）的孔隙率最大值為43.50%，最小值為43.39%；濕基含水率14.61%的小麥堆底層（第8層）的孔隙率最大值為44.00%，最小值為43.90%；濕基含水率16.55%的小麥堆底層（第8層）的孔隙率最大值為44.46%，最小值為44.36%；孔隙率的差值約為0.1%。

圖4 濕基含水率14.61%的小麥堆不同糧層的孔隙率分布值

2.2 平房倉中小麥堆各糧層的平均孔隙率

由有限元模型計算后，提取小麥堆中各單元節(jié)點的應變值，小麥堆糧層平均孔隙率可以根據式（16）計算得到，結果見圖6。

2.2.1 平房倉中小麥堆糧層平均孔隙率與深度的關系

由圖6可知，在同一含水率下，小麥堆的糧層平均孔隙率隨糧層深度增加而減小。這是由于糧層深度增加，糧層壓力增加，體積縮小，導致孔隙率減小。

2.2.2 平房倉中小麥堆糧層平均孔隙率與含水率的關系

由圖6可知，對于同一糧層，糧層平均孔隙率隨含水率的增加而增大。這是由于含水率增加，糧粒膨脹，引起孔隙率增加。

圖6 小麥堆的孔隙率隨糧層深度變化的實驗值與預測值對比

2.3 平房倉中小麥堆糧層孔隙率與含水率和深度的關系模型

設平房倉中小麥堆（淮麥44號）糧層平均孔隙率與糧層深度和含水率的方程為：

式中：εd為小麥堆糧層平均孔隙率，%；ε0為小麥堆表層孔隙率，%；εmin為平房倉中小麥堆糧層平均孔隙率的最小值，%；h為小麥堆糧層深度，m。

依據表2的數據，小麥堆的表層孔隙率與含水率的關系可擬合為線性方程為：

式中：MC為小麥堆的濕基含水率，%。

依據表2的數據，小麥堆最小糧層平均孔隙與含水率的關系可擬合為線性方程：

將式（17）轉化為式（20）。

將式（18）、式（19）及a、b、c、d的值代入式（17），得平房倉中小麥堆平均孔隙率與含水率和糧層深度的關系模型為：

如圖6所示，由式（21）模擬得到的小麥堆各糧層的平均孔隙率（虛線）與實驗計算的平均孔隙率相比相對誤差小于0.15%，RMSE的值在0.07% ~ 0.15%，誤差較小。這個結果表明式（21）表征平房倉中小麥堆糧層孔隙率與含水率和糧層深度的關系是有效的。

3 結 論

（1）在平房倉內的同一糧層中，小麥堆孔隙率隨距倉壁的距離減小而增大，在糧倉的拐角處小麥堆孔隙率達到該層的最大值，在該層的中間位置出現(xiàn)最小值；糧層深度越深，該層的孔隙率分布越不均勻。

（2）平房倉中小麥堆糧層平均孔隙率隨著糧層深度的增加而減小。

（3）平房倉內小麥堆糧層平均孔隙率隨含水率的增加而增大。

（4）平房倉中小麥堆糧層孔隙率與糧層深度和含水率之間關系方程為：