筒倉(cāng)中小麥密度分布值的有限元分析

王 穎， 程緒鐸， 唐福元

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

我國(guó)是全球第一大小麥生產(chǎn)國(guó)。2019年我國(guó)小麥產(chǎn)量1.33億t，消費(fèi)量為1.28億t[1]。小麥儲(chǔ)藏在糧倉(cāng)中，自身的重力與摩擦力、倉(cāng)壁對(duì)小麥的壓力與摩擦力引起了小麥塊的變形與壓縮[2,3]，導(dǎo)致小麥孔隙率減小、堆密度增大。Thompson等[2]研究發(fā)現(xiàn)小麥堆密度隨著內(nèi)部的壓力變化而變化。小麥儲(chǔ)藏在筒倉(cāng)中，小麥堆的壓應(yīng)力分布是不均勻的，尤其是壓應(yīng)力隨著糧層深度的增加而增加，因此筒倉(cāng)中小麥堆的密度分布也是不均勻的，也是隨著糧層深度的增加而增加。確定筒倉(cāng)中小麥堆的密度分布值是準(zhǔn)確計(jì)算筒倉(cāng)中的糧食質(zhì)量的關(guān)鍵參數(shù)，也是預(yù)測(cè)筒倉(cāng)小麥堆通風(fēng)阻力、小麥堆濕熱傳遞規(guī)律的重要參數(shù)[4,5]。

通常用氣體密度儀來(lái)測(cè)量糧堆的堆密度[6-9]。氣體密度儀測(cè)量的只是小糧堆的堆密度(糧堆體積小于10 L)且測(cè)定的是無(wú)壓縮的糧食密度，它無(wú)法被送達(dá)到糧堆深處，因此不能用來(lái)測(cè)量筒倉(cāng)深處的糧食堆密度。我國(guó)一般使用容重儀測(cè)定糧食的密度，但測(cè)定的是無(wú)壓縮的糧食密度，同樣無(wú)法測(cè)定筒倉(cāng)深處受壓縮的糧食密度。Cheng等[10]建立筒倉(cāng)內(nèi)小麥分層壓縮平衡方程，數(shù)值求解了筒倉(cāng)中小麥各層的堆密度。這個(gè)分層壓縮平衡方程是筒倉(cāng)中各層小麥堆密度的平均值，不是某一糧層上各糧塊的堆密度。

ABAQUS是國(guó)際上功能最強(qiáng)的大型通用有限元軟件之一，擁有能夠真實(shí)反映松軟土體的修正劍橋模型，已應(yīng)用于研究土體中的力學(xué)問(wèn)題。小麥堆的力學(xué)特性與松軟土體相似，本研究使用ABAQUS軟件求解筒倉(cāng)中的小麥堆的修正劍橋模型。通過(guò)有限元方法求解筒倉(cāng)中的小麥堆的修正劍橋模型得到筒倉(cāng)中小麥堆密度的三維空間分布值，即小麥堆密度隨糧塊深度(縱向)、徑向、環(huán)向變化的分布值。

1 材料與方法

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

小麥，煙農(nóng)19號(hào)，產(chǎn)自江蘇淮陰，人工剔除實(shí)驗(yàn)樣品中破損籽粒和雜質(zhì)，測(cè)得原始含水率為12.66%。實(shí)驗(yàn)樣品的含水率分別調(diào)制為10.60%、12.66%、14.22%、16.13%。調(diào)制好含水率的樣品裝在塑料袋中封閉好，然后儲(chǔ)藏在人工氣候箱中，人工氣候箱中的溫度設(shè)置為4 ℃。實(shí)驗(yàn)前從人工氣候箱中取出小麥樣品，采用標(biāo)準(zhǔn)烘箱干燥法在130 ℃下干燥實(shí)驗(yàn)樣品19 h測(cè)定它的含水率(ASAE Standards, 2001)，重復(fù)3次。

1.2 小麥堆應(yīng)力與應(yīng)變關(guān)系模型

儲(chǔ)藏在筒倉(cāng)中的小麥堆在外力作用下產(chǎn)生應(yīng)力，而應(yīng)力的作用又引起應(yīng)變。要求解筒倉(cāng)中小麥堆的應(yīng)力分布與應(yīng)變分布必須建立小麥堆應(yīng)力與應(yīng)變關(guān)系模型。修正劍橋模型將總應(yīng)變分為剪切應(yīng)變與體積應(yīng)變，總應(yīng)力分為平均正應(yīng)力與剪切應(yīng)力，適合求解筒倉(cāng)中小麥堆的體積應(yīng)變(應(yīng)變較大)。故本研究選擇修正劍橋模型來(lái)表征小麥堆的應(yīng)力與應(yīng)變關(guān)系。

修正劍橋模型中的應(yīng)變?cè)隽堪w積應(yīng)變?cè)隽縟εv和剪切應(yīng)變?cè)隽縟εs；應(yīng)力增量為平均正應(yīng)力增量dp和剪切應(yīng)力增量dq，q=σ1-σ3是廣義剪切應(yīng)力，p=(σ1+2σ3)/3是平均主應(yīng)力。

修正劍橋模型的彈塑性矩陣的顯式表達(dá)式為[11]：

(1)

Cpp=

(2)

Cpq=Cqp=

(3)

Cqq=

(4)

由式(2)～式(4)可見(jiàn)，修正劍橋模型有6個(gè)參數(shù)：e為初始孔隙比，ν為泊松比，E為彈性模量/kPa，M為臨界狀態(tài)應(yīng)力比，κ為等向膨脹指數(shù)，λ為對(duì)數(shù)硬化模量；這6個(gè)參數(shù)可以通過(guò)三軸壓縮與剪切實(shí)驗(yàn)而確定[11]。

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

1.3.1 初始孔隙比的確定

1.3.1.1 初始孔隙率ε的測(cè)定

小麥初始孔隙率ε是未壓縮時(shí)小麥堆中孔隙體積與小麥堆體積之比，測(cè)定初始孔隙率是為了得到初始孔隙比。本研究使用LKY-1型糧食孔隙率測(cè)定儀(見(jiàn)圖1)測(cè)定小麥的初始孔隙率。

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

如圖1所示，兩個(gè)容積相等的壓力容器A和B，在容器B中裝滿樣品并將其密封。關(guān)閉閥門2，打開(kāi)閥門1、3，用空氣壓縮機(jī)向容器A中鼓入一定壓力的氣體。當(dāng)壓力表指針達(dá)到一定數(shù)值時(shí)，關(guān)閉閥門1，當(dāng)壓力值穩(wěn)定后，記下壓力表讀數(shù)P1；關(guān)閉閥門3，然后打開(kāi)閥門2，當(dāng)容器A和B中壓力平衡時(shí)，記下此時(shí)壓力表讀數(shù)P2。視空氣為理想氣體(標(biāo)準(zhǔn)大氣壓為Pa)，由理想氣體等溫過(guò)程的特性推出小麥初始孔隙率(未壓縮)為：

(5)

1.3.1.2 初始孔隙比e的確定

初始孔隙比e指未壓縮時(shí)小麥堆中孔隙體積與小麥堆中籽粒體積之比，小麥樣品初始孔隙比(未壓縮)可由初始孔隙率推導(dǎo)得到:

e=ε/(1-ε)

(6)

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

對(duì)樣品進(jìn)行軸向壓縮實(shí)驗(yàn)[12]測(cè)定臨界狀態(tài)應(yīng)力比M。在軸向壓縮實(shí)驗(yàn)中，選定5個(gè)圍壓σ3：30、50、70、90、110 kPa。裝好樣品后，啟動(dòng)儀器對(duì)樣品進(jìn)行剪切(剪切應(yīng)變速率為1.000 mm/min)，位移每增加0.4 mm，記錄1次測(cè)力計(jì)讀數(shù)和樣品體積減少量，直至測(cè)力計(jì)讀數(shù)出現(xiàn)最大值時(shí)，記下峰值時(shí)平均主應(yīng)力p和廣義剪切力q。對(duì)每個(gè)圍壓，實(shí)驗(yàn)重復(fù)3次。將這5組實(shí)驗(yàn)的破壞點(diǎn)(即最大主應(yīng)力差)所對(duì)應(yīng)p和q繪制在p-q平面中，經(jīng)一元線性回歸得直線的斜率M值。

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

對(duì)數(shù)硬化模量λ和等向膨脹指數(shù)κ的測(cè)定實(shí)驗(yàn)為各向等壓實(shí)驗(yàn)[12]，使圍壓σ3從0 kPa增加至200 kPa，圍壓每增加5 kPa，記錄1次樣品體積減少量，待圍壓增至200 kPa，再依次卸載至0 kPa，同樣每減小5 kPa，記錄1次樣品體積增加量。將加載曲線和卸載曲線上的p和所對(duì)應(yīng)的孔隙比e繪制在e-lnp平面中，經(jīng)一元線性回歸得直線的斜率λ和κ值。

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

選定圍壓σ3分別為30、50、70、90、110 kPa，對(duì)樣品進(jìn)行軸向加卸載循環(huán)的常規(guī)三軸壓縮實(shí)驗(yàn)[12]，測(cè)定彈性模量E。在常規(guī)三軸壓縮實(shí)驗(yàn)中，糧堆樣品呈圓柱形，軸向壓力為σ1，圍壓為σ3，主應(yīng)力差為Δσ=σ1-σ3，當(dāng)主應(yīng)力差達(dá)到最大值時(shí)糧堆破壞。對(duì)于每一個(gè)圍壓σ3，可測(cè)定一個(gè)最大主應(yīng)力差。對(duì)選定的圍壓σ3，分級(jí)施加軸向壓力，每級(jí)壓力取試樣最大主應(yīng)力差的10%；施加第1級(jí)壓力，并開(kāi)動(dòng)秒表，記錄加壓1 min后位移計(jì)的讀數(shù)，每1 min施加一級(jí)壓力，記錄位移計(jì)讀數(shù)，直至施加到第4級(jí)壓力為止；在施加第4級(jí)壓力1 min后記錄位移計(jì)的讀數(shù)，并逐級(jí)卸壓，每1 min卸去1級(jí)壓力，并記錄卸壓后1 min時(shí)位移計(jì)的讀數(shù)，直至卸去所有施加的軸向壓力；重復(fù)加卸載荷4次。

彈性模量E按第4次加卸載荷記錄的數(shù)據(jù)計(jì)算公式為：

(7)

式中：Δhe為小麥樣品的軸向彈性變形量/mm；hc為未軸向加載時(shí)小麥樣品高度/mm。

如同測(cè)定等向膨脹指數(shù)κ一樣，對(duì)小麥樣品進(jìn)行各向等壓實(shí)驗(yàn)，體變彈性模量K按卸載記錄的數(shù)據(jù)計(jì)算公式為：

(8)

式中：ΔV為小麥的彈性體積增量/m3；V為未加壓力時(shí)小麥樣品的體積/m3。

泊松比ν按公式計(jì)算：

(9)

1.4 筒倉(cāng)中小麥堆密度分布的有限元求解方法

1.4.1 筒倉(cāng)的幾何尺寸與力學(xué)參數(shù)

筒倉(cāng)的幾何尺寸與力學(xué)參數(shù)見(jiàn)表1。

表1 筒倉(cāng)幾何尺寸及其材料的參數(shù)

1.4.2 筒倉(cāng)與小麥堆單元網(wǎng)格劃分

筒倉(cāng)高31 m, 重力未施加時(shí)筒倉(cāng)內(nèi)小麥堆高度29 m。筒倉(cāng)與小麥堆的形狀和受力都是軸對(duì)稱的，取筒倉(cāng)與小麥堆的一徑向平面表示筒倉(cāng)與小麥堆，選擇四邊形單元?jiǎng)澐衷搹较蚱矫妫倪呅螁卧拿娣e為1×1=1 m2,每一個(gè)四邊形單元代表一個(gè)截面為矩形的環(huán)形體。小麥堆模塊被劃分為29×5個(gè)單元；倉(cāng)壁模塊被劃分成31×1個(gè)單元。

1.4.3 小麥堆與筒倉(cāng)壁的接觸模擬

倉(cāng)壁內(nèi)側(cè)面與小麥堆外側(cè)面有相互力的作用。接觸面之間傳遞切向應(yīng)力與法向應(yīng)力，切向應(yīng)力即摩擦力。在有限元分析中應(yīng)考慮阻止兩表面之間相對(duì)滑動(dòng)而產(chǎn)生的摩擦力。庫(kù)侖模型是常用的摩擦模型，通過(guò)摩擦系數(shù)來(lái)描述兩表面間的摩擦行為。臨界剪應(yīng)力與法向接觸壓力之間的關(guān)系為：

τ=μp

(10)

式中：μ為摩擦系數(shù)；p為法向接觸壓力/kPa。

1.4.4 載荷及邊界約束

筒倉(cāng)中小麥堆受到自重、內(nèi)摩擦、倉(cāng)壁對(duì)小麥的抵抗等多種力的作用，小麥堆的重力由體積力進(jìn)行描述，體積力是主動(dòng)力。小麥堆內(nèi)部作用力通過(guò)小麥堆的應(yīng)力與應(yīng)變的本構(gòu)關(guān)系進(jìn)行描述，小麥堆與倉(cāng)壁之間的相互作用力由摩擦力與正壓力進(jìn)行描述。將小麥堆模型底部邊界條件設(shè)置為U2=0，即對(duì)筒倉(cāng)底板處節(jié)點(diǎn)豎直方向的位移約束。因?yàn)樵趦?chǔ)藏過(guò)程中，筒倉(cāng)是不產(chǎn)生任何運(yùn)動(dòng)的，故將倉(cāng)壁處節(jié)點(diǎn)位移完全約束，以限制其剛體位移。

1.4.5 筒倉(cāng)中小麥堆各處密度的計(jì)算

小麥堆共有n×m個(gè)單元，n(29)為行數(shù)，m(5)為列數(shù)。重力未施加時(shí)各單元密度是相同的；重力加載后，每個(gè)單元在筒倉(cāng)中的位置及其體積都發(fā)生了變化。若i表示行(縱向)編號(hào)，j表示列(橫向)編號(hào)，則第i行j列的單元的密度可由公式計(jì)算：

(11)

2 結(jié)果與分析

2.1 初始堆密度(未壓縮)

根據(jù)國(guó)家標(biāo)準(zhǔn)GB/T 5498—1985測(cè)定得到含水率分別為10.60%、12.66%、14.22%、16.13%的煙農(nóng)19號(hào)表層密度(未壓縮)分別為801.75、787.09、770.47、755.65 kg/m3。

2.2 修正劍橋模型參數(shù)

使用應(yīng)變式三軸儀和糧食孔隙測(cè)定儀測(cè)定煙農(nóng)19號(hào)的修正劍橋模型參數(shù)參數(shù)如表2所示。

表2 修正劍橋模型參數(shù)

2.3 筒倉(cāng)中小麥堆的密度分布值

根據(jù)國(guó)家標(biāo)準(zhǔn)GB/T 5498—1985測(cè)定得到含水率分別為10.60%、12.66%、14.22%、16.13%的小麥堆(煙農(nóng)19號(hào))初始密度(表層密度)分別為801.75、787.09、770.47、755.65 kg/m3。使用修正劍橋模型和有限元軟件ABAQUS計(jì)算得到筒倉(cāng)中小麥堆的密度分布值，見(jiàn)圖2。

圖2 不同含水率小麥堆的密度與半徑的關(guān)系

由圖2可知，筒倉(cāng)中小麥堆的密度隨著糧層深度的增加而逐漸增大，而在筒倉(cāng)拐角處的小麥堆密度隨著糧層深度的增加而逐漸減小。在上部糧層，隨著糧塊與筒倉(cāng)中心軸距離的增加，密度緩慢減小，減小值約為0.07%～0.1%。在中部糧層，隨著糧塊與筒倉(cāng)中心軸距離的增加，密度幾乎不變。在接近倉(cāng)底層，隨著糧塊與筒倉(cāng)中心軸距離的增加，密度逐漸減小，減小值約為0.4%～0.6%。

2.4 筒倉(cāng)中小麥堆平均密度與糧層深度的關(guān)系

重力未施加時(shí)筒倉(cāng)內(nèi)小麥堆高度29 m。筒倉(cāng)與小麥堆的形狀和受力都是軸對(duì)稱的，取筒倉(cāng)與小麥堆的一徑向平面表示筒倉(cāng)與小麥堆，選擇四邊形單元?jiǎng)澐衷搹较蚱矫妫倪呅螁卧拿娣e為1 m2,每一個(gè)四邊形單元代表一個(gè)截面為矩形的環(huán)形體。小麥堆模塊被劃分為29×5個(gè)單元；倉(cāng)壁模塊被劃分成31×1個(gè)單元。這樣小麥堆共有29行5列。筒倉(cāng)中小麥堆第n層(行)平均密度按式(12)計(jì)算：

(12)

依據(jù)修正劍橋模型計(jì)算的數(shù)據(jù)，由式(12)計(jì)算得到小麥堆各糧層的平均密度，如圖3所示。在小麥含水率為10.60%、12.66%、14.22%、16.13%時(shí)，小麥堆的糧層深度的變化范圍分別為0～27.704 8、0～27.661 0、0～27.485 8、0～26.982 1 m；其平均密度的變化范圍分別為800.75～822.08、786.02～807.61、769.32～792.69、754.33～783.23 kg/m3。在同一含水率下，小麥堆的糧層平均密度隨著糧層深度的增大而逐漸增大，逐漸趨于恒定值。

圖3 不同含水率下小麥堆的平均密度與糧層深度的關(guān)系

2.5 筒倉(cāng)中小麥堆平均密度與含水率的關(guān)系

從圖3可以看出，含水率為10.60%、12.66%、14.22%、16.13%的小麥堆在相同的糧層深度下的平均密度是隨著含水率的增大而減小的。

2.6 筒倉(cāng)中小麥堆平均密度與含水率、糧層深度關(guān)系模型

有限元計(jì)算的密度數(shù)據(jù)有這樣的規(guī)律，當(dāng)深度一定時(shí)密度與含水率大致成線性關(guān)系，當(dāng)含水率一定時(shí)密度隨深度增加而增大，增加率隨深度增加而減小，故用1-eah3+bh2+ch表達(dá)密度隨深度增加的變化趨勢(shì)。因此，設(shè)小麥糧層平均密度與糧層深度、含水率的關(guān)系方程為：

ρ=ρ0+(ρmax-ρ0)(1-eah3+bh2+ch)

(13)

式中：ρ為小麥堆糧層平均密度/kg/m3；ρ0為小麥堆初始密度/kg/m3；ρmax為小麥堆壓縮糧層平均密度的最大值/kg/m3；h為小麥堆糧層深度/m。

表3 入倉(cāng)小麥的參數(shù)

小麥堆初始密度與含水率的關(guān)系可擬合為線性方程(14)。

ρ0=892.35-8.56MCρ=0.99

(14)

式中：MC為小麥堆含水率/%。

通過(guò)式(12)計(jì)算得到的小麥堆糧層平均密度的最大值變化范圍為822.08～783.33 kg/m3，其隨含水率的增大而降低。小麥堆糧層平均密度的最大值與含水率的關(guān)系可擬合為線性方程(15)。

ρ0=897.99-7.20 MCR2=0.99

(15)

為了確定模型(13)中a、b、c這3個(gè)模型常數(shù)，將式(13)變換成式(16)：

(16)

以ln(1-(ρ-ρ0)/(ρmax-ρ0))為縱坐標(biāo)，h為橫坐標(biāo)，按式(16)進(jìn)行數(shù)據(jù)擬合從而得到模型常數(shù)a、b、c分別是-0.000 47、0.012、-0.24。

將a、b、c的值及式(14)和(15)代入式(13)，可得小麥堆糧層的平均密度與含水率、糧層深度的關(guān)系模型為：

ρ=892.35-8.56MC+(5.64+1.36MC)(1-e-0.000 47h3+0.012h2-0.24h)

(17)

圖3中的散點(diǎn)是有限元方法計(jì)算的密度，曲線表示模型方程(17)計(jì)算的密度值。由圖3可以看出，經(jīng)模型方程(17)計(jì)算得到的小麥堆糧層密度與限元方法計(jì)算的糧層密度誤差較小，對(duì)于含水率分別為10.60%、12.66%、14.22%、16.13%的小麥，模型方程(17)計(jì)算的小麥堆糧層平均密度與有限元方法計(jì)算的糧層平均層密度的誤差小于0.5%，說(shuō)明模型方程(17)有效。

2.7 實(shí)倉(cāng)驗(yàn)證

實(shí)例一：中央儲(chǔ)備糧溫州直屬庫(kù)11、12、13、14號(hào)筒倉(cāng)內(nèi)徑為10 m，入倉(cāng)小麥的參數(shù)見(jiàn)表3。本研究通過(guò)有限元分析方法求解修正劍橋模型計(jì)算得到小麥的糧層密度分布值，再由積分方程計(jì)算出儲(chǔ)糧總質(zhì)量，并與筒倉(cāng)中實(shí)際儲(chǔ)糧質(zhì)量進(jìn)行比較，結(jié)果見(jiàn)表4。

表4 修正劍橋模型計(jì)算結(jié)果與實(shí)際儲(chǔ)糧數(shù)量的比較

表5 入倉(cāng)小麥的參數(shù)

實(shí)例二：中央儲(chǔ)備糧金華直屬庫(kù)1、2號(hào)淺圓倉(cāng)內(nèi)徑為25 m，入倉(cāng)小麥的參數(shù)見(jiàn)表5。本研究通過(guò)有限元分析方法求解修正劍橋模型計(jì)算得到小麥的糧層密度分布值，再由積分方程計(jì)算出儲(chǔ)糧總質(zhì)量，并與筒倉(cāng)中實(shí)際儲(chǔ)糧質(zhì)量進(jìn)行比較，結(jié)果見(jiàn)表6。

表6 修正劍橋模型計(jì)算結(jié)果與實(shí)際儲(chǔ)糧數(shù)量的比較

由表4與表6可以看出，本研究的修正劍橋模型計(jì)算結(jié)果與實(shí)際賬面數(shù)誤差小于0.3%，說(shuō)明本研究的修正劍橋模型計(jì)算的密度分布值可靠性較高。

入倉(cāng)方式對(duì)糧食分級(jí)有影響，糧食分級(jí)對(duì)密度的分布有影響，糧食入倉(cāng)以后，由于糧食分級(jí)，同一層的糧食密度差別增大，但一層的平均密度變化很小，所以對(duì)計(jì)算糧食總質(zhì)量影響較小，本研究的模型計(jì)算結(jié)果表明影響較小。糧食含雜量對(duì)壓縮密度有影響，本研究選擇的實(shí)倉(cāng)中的小麥含雜量較小(0.2%～1.0%)，所以對(duì)計(jì)算糧食總質(zhì)量影響較小，本研究的模型計(jì)算結(jié)果表明影響較小。

3 結(jié)論

筒倉(cāng)中小麥堆的密度隨著糧層深度的增加而逐漸增大，而在筒倉(cāng)拐角處的小麥堆密度隨著糧層深度的增加而逐漸減小。在上部糧層，糧塊密度隨著糧塊與筒倉(cāng)中心軸距離的增加而緩慢減小，減小值為0.07%～0.10%。在中部糧層，糧塊密度隨著糧塊與筒倉(cāng)中心軸距離的增加幾乎不變。在接近倉(cāng)底糧層，糧塊密度隨著糧塊與筒倉(cāng)中心軸距離的增加而逐漸減小，減小值為0.4%～0.6%。小麥堆的糧層平均密度隨著糧層深度的增大而逐漸增大，逐漸趨于恒定值。相同深度下的糧層平均密度隨著含水率的增大而減小。筒倉(cāng)中小麥堆糧層平均密度與含水率、糧層深度之間關(guān)系方程為：ρ=892.35-8.56MC+(5.64+1.36MC)(1-e-0.000 47h3+0.012h2-0.24h)。實(shí)倉(cāng)驗(yàn)證得到，由修正劍橋模型計(jì)算的筒倉(cāng)中小麥的總質(zhì)量與實(shí)際糧倉(cāng)的賬面小麥總質(zhì)量誤差小，相對(duì)誤差在0.14%～1.25%之間。