邵小龍，程緒鐸

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

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

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

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

到目前為止，還未見平房倉中糧堆的孔隙率分布的理論與實(shí)驗(yàn)成果報(bào)道。本文選定修正劍橋模型為小麥堆的應(yīng)力與應(yīng)變關(guān)系本構(gòu)模型，使用有限元方法求解平房倉中小麥堆修正劍橋模型，計(jì)算出平房倉內(nèi)小麥堆的應(yīng)變分布值，由應(yīng)變分布值與表層孔隙率（無壓縮）計(jì)算出平房倉中小麥堆的孔隙率分布值。

1 材料與方法

1.1 實(shí)驗(yàn)材料

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

1.2 實(shí)驗(yàn)儀器

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

1.3 小麥堆的應(yīng)力與應(yīng)變關(guān)系本構(gòu)模型

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

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

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

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

1.4 修正劍橋模型參數(shù)的確定

1.4.1 初始孔隙比e0的測(cè)定

小麥堆孔隙率ε為小麥堆中的孔隙體積與整個(gè)小麥堆體積之比。孔隙比e指小麥堆中孔隙體積與小麥堆籽粒體積之比。本實(shí)驗(yàn)通過LKY-1型糧食孔隙率測(cè)定儀（見圖1）測(cè)定小麥堆的初始孔隙率，初始孔隙比可由初始孔隙率計(jì)算得到。

圖1 孔隙率測(cè)定儀示意圖

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

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

1.4.2 臨界狀態(tài)應(yīng)力比M的測(cè)定

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

1.4.3 對(duì)數(shù)硬化模量λ和等向膨脹指數(shù)κ的測(cè)定

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

1.4.4 彈性模量E和泊松比υ的測(cè)定

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

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

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

（2）泊松比υ的測(cè)定：糧堆進(jìn)行等向壓縮實(shí)驗(yàn)，由式（11）計(jì)算出體變彈性模量[14]。

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

則泊松比υ為：

1.5 平房倉中小麥堆應(yīng)力與應(yīng)變的有限元計(jì)算方法

平房倉內(nèi)儲(chǔ)藏小麥后，小麥堆各處都產(chǎn)生了應(yīng)力與應(yīng)變，小麥堆各處的孔隙率與初始孔隙率及各處的應(yīng)變密切相關(guān)。只要測(cè)定出小麥堆的初始孔隙率并計(jì)算出小麥堆各處的應(yīng)變，就能計(jì)算出小麥堆各處的孔隙率。小麥堆各處的應(yīng)變可使用有限元方法求解小麥堆的修正劍橋模型而得出。

1.5.1 平房倉倉壁的幾何與力學(xué)參數(shù)的選定

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

表1 平房倉倉壁的幾何與力學(xué)參數(shù)

1.5.2 平房倉中小麥堆有限元模型構(gòu)建與運(yùn)行步驟

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

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

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

（4）在 Step步驟中設(shè)定分析步。

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

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

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

（8）在 Job 步驟中對(duì)構(gòu)建的模型進(jìn)行運(yùn)算。

（9）待運(yùn)算成功后，在Visualization步驟中查看運(yùn)算結(jié)果提取所需數(shù)據(jù)。

1.6 平房倉中小麥堆孔隙率的計(jì)算

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

小麥堆經(jīng)修正劍橋模型的有限元計(jì)算后提取各個(gè)單元節(jié)點(diǎn)的應(yīng)變值，其單元孔隙率為：

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

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

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

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

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

式中：Vijk為第ijk個(gè)單元所在處的單元體積，m3；m為沿寬度方向的單元數(shù)；n為沿長(zhǎng)度方向的單元數(shù)。

2 結(jié)果與分析

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

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

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

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

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

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

在各個(gè)含水率中，同一糧層深度下的小麥堆的孔隙率分布是不均勻的。在同一糧層中，小麥堆的孔隙率由倉壁至糧層中間逐漸減小，在該層中間的孔隙率最小，在平房倉的拐角處孔隙率陡然增大并達(dá)到該層的最大值。不同層之間的孔隙率也存在明顯差異。糧層深度增加，同一水平面位置處的孔隙率減小；糧層深度增加，糧層各處的孔隙率的差值增大。同一糧層各處的孔隙率差值很小，濕基含水率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 平房倉中小麥堆各糧層的平均孔隙率

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

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

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

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

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

圖6 小麥堆的孔隙率隨糧層深度變化的實(shí)驗(yàn)值與預(yù)測(cè)值對(duì)比

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

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

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

依據(jù)表2的數(shù)據(jù)，小麥堆的表層孔隙率與含水率的關(guān)系可擬合為線性方程為：

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

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

將式（17）轉(zhuǎn)化為式（20）。

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

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

3 結(jié) 論

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

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

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

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