筒倉中稻谷的空隙率分布研究

唐福元 許 倩 程緒鐸

(南京財經大學食品科學與工程學院;江蘇省現代糧食流通與安全協同創新中心, 南京 210046)

筒倉中稻谷的空隙率分布研究

唐福元 許 倩 程緒鐸

(南京財經大學食品科學與工程學院;江蘇省現代糧食流通與安全協同創新中心, 南京 210046)

采用LHT-1糧食回彈模量儀測定稻谷(南粳5055)堆的表層密度及壓縮密度，建立筒倉中稻谷堆的密度、應力與糧層深度關系的微分方程組，用數值方法計算筒倉中稻谷密度與糧層深度關系。采用糧食孔隙率測量儀測定表層稻谷(無壓縮)孔隙率，由表層孔隙率，表層密度及筒倉深處的密度計算出筒倉中稻谷孔隙率與糧層深度關系。計算結果表明：在直徑20米的筒倉中，在30米的筒體部分，南粳5055空隙率變化范圍為61.00%～56.32%，在10 m的錐斗中，空隙率變化范圍為56.32%～59.77%；在帶錐斗筒倉的筒體部分，稻谷堆孔隙率隨著糧層深度的增加而減小；到錐斗部分，稻谷堆孔隙率隨著糧層深度的增加而逐漸增大。在不同直徑的筒倉的筒體部分，在同一深度，稻谷堆孔隙率隨著筒倉直徑的增大而減小。

筒倉 稻谷 孔隙率 密度

稻谷是我國單產最高、總產最多的糧食作物，在國家糧食供應與安全中具有重要的地位[1]。筒倉是儲糧倉型中機械化程度最高的一種倉型，被廣泛應用于稻谷的儲藏[2-3]。稻谷儲藏在筒倉中，受到自重、內摩擦力以及倉壁的支持力，稻谷堆的各處產生壓應力與切應力，從而產生體積應變與形狀應變，體積縮小，密度增大，孔隙率減小。隨著糧層深度的增加，稻谷的密度增大[4-5]，孔隙率減小。準確地計算出筒倉中糧堆孔隙率的分布值能給筒倉中糧堆各處的通風阻力的計算提供技術參數。 糧堆的孔隙率的分布還影響糧堆中的熱量與水分的遷移。因此，計算出筒倉中孔隙率的分布值對糧食安全儲藏具有重要意義。

糧食孔隙率的測定始于二十世紀六十年代[6]。早在1967年Thompson等[6]使用氣體比重瓶測定了燕麥、小麥、高粱、大豆、黑麥、大麥、玉米和紫花苜蓿種子的孔隙率，但是他們測定的樣品容積只有50 mL。Gustafson[7]使用氣體比重瓶測定了玉米的孔隙率。 Chang[8]使用氣體比重瓶測定了不同含水率的玉米、小麥和高粱的孔隙率。目前測定糧食孔隙率的方法是氣體置換法，其測定裝置主要有定容型[9]、變容型[10]和比較型[11]。這些裝置的容積都很小(5 L以下)，只能測量小糧堆(5 L以下)的孔隙率，不能測定糧倉內糧堆孔隙率。郝倩等[12]自制小倉測定了倉內小麥、玉米的孔隙率，該倉的容積只有50 L，測定的只是小麥和玉米的平均孔隙率。測定孔隙率還有光學法[13-18]、聲學法[19]等，這些方法也不適合用來測定糧倉內糧堆孔隙率。

由于筒倉中的糧堆結構的復雜性，揭示糧倉中糧堆的孔隙率的分布是個難題。到目前為止，糧倉中糧堆的孔隙率分布的理論與試驗成果鮮有報道。本研究采用LHT-1糧食回彈模量儀來測定出稻谷堆壓縮密度，建立筒倉中稻谷堆的密度、壓應力與糧層深度關系的微分方程組，用數值方法計算得到筒倉中稻谷密度隨糧層深度變化關系，再由表層孔隙率，表層密度及筒倉深處的密度計算出筒倉中稻谷孔隙率與糧層深度關系。

1 材料與方法

1.1 試驗材料

南粳5055稻谷，原始含水率為12.95% (濕基)，2016年產于南京。2016年11月至12月進行壓縮試驗。

試驗時，將其含水率調節為11.77%、12.95%、14.47%、15.75%、17.58% w.b.。稻谷籽粒的最大直徑、中直徑、最小直徑的平均值分別為6.8、3.4、2.3 mm，對應的標準差分別為0.2、0.2、0.1 mm。(隨機取樣100粒稻谷粒，用游標卡尺測得其3個粒徑，精確到0.1 mm，最后求其平均值。)

1.2 試驗儀器

LHT-1型糧食回彈模量測定儀：南京土壤儀器廠有限公司；LKY-1型糧食孔隙率測定儀：南京土壤儀器廠有限公司；HG202-2(2A/2AD)電熱干燥箱：南京盈鑫實驗儀器有限公司。

1.3 試驗方法

1.3.1 稻谷堆壓縮密度的測定

將稻谷裝入LHT-1糧食回彈模量儀的圓筒中，頂部施加壓力，筒壁對稻谷產生向上的摩擦力，裝樣筒內各糧層所受到的豎直應力與側向應力不同。 本試驗用圓筒中稻谷的平均豎直應力與側向應力表示筒倉中某一糧層的豎直應力與側向應力。

采用微元法對裝樣筒內糧堆樣品進行受力分析，假設裝樣筒內半徑為Rc，壓縮樣品高度為Hc，取距離糧面深度為y，高度為dy的微元體進行受力分析[5]，微元體在豎直方向的受力分析見圖1：

圖1 回彈模量儀裝樣筒中糧堆微元層受力圖

稻谷在裝樣筒內受壓時處于主動應力狀態，根據Rankine理論，側向應力與豎直應力的比為：

(1)

式中：ph1為圓筒中糧堆所受側向應力/kPa；pv1為圓筒中糧堆所受豎直應力/kPa；φ為稻谷內摩擦角/°。

糧堆微元體在豎直方向上受力平衡，可得平衡方程為：

(2)

式中：ρ為稻谷密度/kg/m3；μc為稻谷與裝樣筒筒壁的摩擦系數；Ac為裝樣筒的橫截面積/m2。

(3)

式中：p0為糧堆頂部壓應力/kPa。

在裝入稻谷的圓筒頂部逐級加載，測定每次加載后稻谷樣品的高度，計算出對應的稻谷密度，并根據式(3)計算出對應的豎直應力，即可得到稻谷樣品密度與其所受豎直應力(最大主應力)的關系:

ρ=F(pmax)

(4)

將回彈模量儀放在實驗室中，室內溫度設置為20 ℃。將樣品勻速倒入裝樣筒中，并將其表面鋪平，放上傳壓板，保證傳壓板上表面與裝樣筒上端齊平。保持橫梁杠桿垂直，轉動平衡錘調整杠桿至水平以上，用M16螺母固定平衡錘。旋轉傳壓螺釘與傳壓板接觸，調整0～30 mm位移傳感器的觸頭位置，調零百分表。根據試驗要求，對試樣進行加載。隨著試樣的下沉，杠桿向下傾斜，為防止杠桿傾斜影響加載精度，調節調平手輪，使杠桿處于水平位置。試驗進行3 d后記錄樣品高度，計算樣品壓縮后的體積，結束試驗，倒出裝樣筒內的稻谷樣品，稱重并記錄稻谷質量。稻谷堆的密度即為稻谷的質量除以壓縮后的體積[5]。

1.3.2 稻谷堆表層孔隙率測定

糧食孔隙率ε為糧食中孔隙與糧堆體積之比。本試驗使用LKY-1型糧食孔隙率測定儀(見圖2)測定表層稻谷的孔隙率。

圖2 孔隙率測定儀示意圖

如圖2所示，兩個容積相等的壓力容器A和B，在容器B中裝滿樣品并將其密封。關閉閥門2，打開閥門1、3，用空氣壓縮機向容器A中鼓入一定壓力的氣體。當壓力表指針達到一定數值時，關閉閥門1，當壓力值穩定后，記下壓力表讀數P1；關閉閥門3，然后打開閥門2，當容器A和B中壓力平衡時，記下此時壓力表讀數P2。視空氣為理想氣體，由理想氣體等溫過程的特性推出孔隙率如下：

1.4 筒倉內稻谷密度、壓應力與糧層深度關系的微分方程組的構建

設帶錐斗筒倉的直徑為D，圓筒高為H，錐斗高為h，錐斗壁與水平面夾角為α，稻谷的內摩擦角為φ，倉墻與稻谷間的摩擦系數為μ，設水平壓應力與豎直應力大小不隨半徑的大小和方向而變化。帶錐斗的筒倉和稻谷表面下深度y處稻谷的微元層受力示意圖如圖3所示。

筒體部分，其微元層的受力如圖3b所示。微元層在豎直方向上合力為零，則得到：

(pv+dpv)A+μphCdy-pvA-ρgAdy=0

(5)

式中：pv為稻谷豎直應力/kPa；ph為側應壓力/kPa；ρ為稻谷的密度/kg/m3；g為重力加速度/m/s-2；A為圓筒的橫截面面積/m2；C為圓筒的周長/m；且

C=A/R

(6)

式中：R為水力半徑/m；R=D/4。

圖3 筒倉的筒體中微元體受力圖

稻谷儲藏在筒倉中，稻谷處于主動應力狀態，稻谷最小主應力(即筒壁側壓力ph)和最大主應力(即豎直應力pv)之比k為：

(7)

結合式(5) ～ (7)得到筒體部分最大主應力pv的微分方程：

(8)

錐斗部分，其截面圖和微元層受力如圖4所示。在同一水平面上的錐斗部分微元層，其糧堆的最大主應力的大小與方向都是變化的。在錐斗壁處稻谷處于被動應力狀態，此處最大主應力方向垂直于錐斗壁。筒體底層中心的最大主應力與錐斗頂面中心最大主應力是相等的。錐斗內的任一水平微元層,中心豎直應力(pvc)是該處的最大主應力，邊緣豎直應力并不是該處的最大主應力，也不是該處的最小主應力，在錐斗壁處(微元層邊緣處)取一微元體如圖4b、4c 所示，微元體受到的最大和最小壓應力分別為phmax和phmin，phmax為錐斗壁側壓力，且微元體的最小壓應力和最大壓應力之比為k。

圖4 筒倉錐斗部分微元體受力圖

根據Rankine壓力理論得到水平稻谷薄層邊緣處豎直應力pve：

(9)

令

(10)

同一水平微元層的豎直應力取平均值phv：

phv=(pve+pvc)/2=(λ+1)phmax/2

(11)

作用于微元層上豎直方向的合力等于0：

phvAa+ρgV-(phv+dphv)Ab-phmaxScosα=0

(12)

式中：Aa為上底面的面積/m-2；Ab為下底面的面積/m-2；S為側表面積/m-2；V為微元體的體積/m-2。

設水平面上稻谷薄層的半徑為r，則

Aa=πr2

(13)

Ab=π(r-dycotα)2=π(r2-2rcotαdy)

(14)

(15)

(16)

結合(9) ～ (16)式并化簡得到錐斗部分最大主應力(錐斗壁側壓應力)方程：

(17)

結合式(7)、(8)、(4)得到筒體部分壓應力、密度與糧層深度關系的方程組：

(18)

結合式(10)、(11)、(17)、(4)得到錐斗部分壓應力、密度與糧層深度關系的方程組：

(19)

式中：phv為錐斗中糧堆豎直應力/kPa；phmax為錐斗壁側壓應力/kPa。

方程組(18)、(19)的連續條件為：

pv(H)=phmax(H)

(20)

方程組(18)、(19)的邊界條件為：

(21)

1.5 方程組(18)～(21)數值計算步驟

1.5.1 筒倉中稻谷堆分層

筒倉中稻谷堆分為n層, 筒體中n1層，漏斗中n2層。每層的厚度(高)為dy, 層數為：n1=H1/dy;n2=H2/dy; n=n1+n2。其中，H1為筒體的高，H2為錐斗的高。

1.5.2 第一層密度與豎直壓應力

依據方程21得，ρ1=ρ0(表層密度),pv,1=0 (表層豎直壓應力)。

1.5.3 計算筒體中第2層到第n1的壓應力與密度

a.由方程組18的第3式計算壓應力增量dp；

b.計算下一層壓應力：(pv)next layer= (pv)current layer+dp；

c.由方程組18的第1式計算該層密度。

1.5.4 計算漏斗中各層的壓應力與密度

a.由方程(20)知，漏斗中第1層的最大主應力等于筒體中最下層的豎直壓應力；

b.由方程組19的第2式計算壓應力增量dphmax；

c.計算下一層主應力：(phmax)next layer= (phmax)current layer+dphmax

d.由方程組19的第1式計算該層密度。

1.6 糧食孔隙率與密度的轉換關系

糧堆表層孔隙率為：

(22)

由(22)可得:

(23)

糧堆深層密度為：

(24)

由(24)得：

(25)

糧堆深層孔隙率為：

(26)

式中：e0為表層孔隙率，e為糧堆孔隙率，V孔為孔隙的體積/m3；m為糧堆質量/kg；V總0為糧堆總體積/m3；ρ為糧堆密度/kg/m3；ρ0為糧堆表層密度/kg/m3。

2 結果與分析

2.1 稻谷堆豎直壓應力與壓縮密度

2.1.1 稻谷堆豎直壓應力

由直剪儀測得南京5055稻谷的內摩擦角和摩擦系數的試驗結果見表1。

表1 不同豎直壓應力下稻谷的內摩擦角和摩擦系數

使用回彈模量儀通過逐級加載的方式對稻谷堆進行壓縮試驗，加載的頂部豎直應力分別為50、100、150、200、250、300 kPa。根據表1中的參數與試驗測定的樣品壓縮的頂部壓力與樣品高度，根據式(3)計算得到稻谷堆所受的平均豎直應力，結果見表2。

表2 不同含水率、不同頂部壓力下的稻谷堆所受的平均豎直壓應力

2.1.2 稻谷堆的壓縮密度

不同含水率、不同頂部壓應力下的稻谷堆壓縮后測定的密度見表3。

表3 不同含水率、不同頂部壓力下的稻谷堆的壓縮密度

2.2 表層稻谷孔隙率

糧食孔隙率測定儀測定出南粳5055(含水率12.95% w.b.)表層稻谷孔隙率e0為0.610 0，標準差為0.000 3。

2.3 稻谷堆壓縮密度與平均豎直應力的關系

由表2和表3擬合稻谷堆密度與平均豎直應力的關系方程為ρ=apv2+bpv+c[ρ=F(pmax)]，其中，pv為豎直壓應力，為壓縮密度，pmax為最大主應力，方程系數a、b、c隨含水率的變化而變化，結果見表4。

表4 稻谷堆密度與平均豎直壓應力的擬合方程系數

2.4 筒倉中稻谷堆密度與糧層深度的關系

設筒倉的參數條件如表5所示。

表5 筒倉的參數

儲藏稻谷的含水率為12.9%，內摩擦角為40°，與倉壁的摩擦系數為0.4，表面密度為560 kg/m3。頂部糧面水平，選擇糧堆的每一薄層為0.1 m。數值法求解微分方程組(18) 和(19)，從而計算得到不同糧層深度下稻谷堆的密度,如表6所示。

表6 不同糧層深度下稻谷堆的密度

由表6知，在筒倉的筒體部分，稻谷密度隨著糧層深度的增加而增增加；在錐斗部分，隨著糧層深度的增加，稻谷堆密度逐漸減小。

2.5 筒倉中稻谷堆孔隙率與糧層深度的關系

選定筒倉的參數如表5所示，選定稻谷是南粳5055，2015年產于南京，稻谷的含水率為12.9% w.b.，其內摩擦角為40°，與倉壁的摩擦系數為0.4，稻谷的表面密度為560 kg/m3。

表層孔隙率e0可按1.2試驗方法測得，表層密度ρ0可按1.1試驗方法測得，經過公式(26)的轉換可以由稻谷堆密度(表6)計算出筒倉內稻谷堆孔隙率與糧層深度的關系，如表7所示。表6與表7的結果密切相關，密度大處孔隙率就小。

表7 不同糧層深度下稻谷堆的孔隙率

由表7可知，在直徑10米的筒倉中，在30米的筒體部分，南粳5055空隙率變化范圍為61.00%～57.03%，在5米的錐斗中，空隙率變化范圍為57.03%～58.25%；在直徑20米的筒倉中，在30米的筒體部分，南粳5055空隙率變化范圍為 61.00%～56.32%，在10米的錐斗中，空隙率變化范圍為56.32%～59.77%；在直徑40米的筒倉中，在30米的筒體部分，南粳5055空隙率變化范圍為61.00%～55.95%。

在筒倉的筒體部分，稻谷堆孔隙率隨著糧層深度的增加而增減小，這是由于壓應力隨著糧層深度的增加而增加的緣故。在錐斗部分，隨著糧層深度的增加，稻谷堆孔隙率逐漸增大；這是由于壓應力隨著糧層深度的增加而增加的緣故。在不同直徑的筒倉的筒體部分，在同一深度，稻谷堆孔隙率隨著筒倉直徑的增大而減小。這是由于在同一深度，稻谷堆壓應力隨著筒倉直徑的增大而增大的緣故。

3 結論

3.1 在筒倉的筒體部分，稻谷密度隨著糧層深度的增加而增增加；在錐斗部分，隨著糧層深度的增加，稻谷堆密度逐漸減小。

3.2 在直徑20 m的筒倉中，在30米的筒體部分，南粳5055空隙率變化范圍為61.00%～56.32%，在10米的錐斗中，空隙率變化范圍為56.32%～59.77%。

3.3 在筒倉的筒體部分，稻谷堆孔隙率隨著糧層深度的增加而增減小；在錐斗部分，隨著糧層深度的增加，稻谷堆孔隙率逐漸增大。

3.4 在不同直徑的筒倉的筒體部分，在同一深度，稻谷堆孔隙率隨著筒倉直徑的增大而減小。

[1]王瑞元，朱永義，謝健，等. 我國稻谷業加工現狀與展望[J]. 糧食與飼料工業, 2011 (3): 1-5

Wang Ruiyuan, Zhu Yongyi, Xie Jian, et al. Present situation andprospect of rice processing in China [J]. Food and Feed Industry, 2011 (3): 1-5

[2]張承光. 現代化儲糧倉型的選擇[J]. 糧油倉儲科技通訊, 2002 (2): 42-44

Zhang Chengguang. Selection of modern granary type [J]. Grain and Oil Storage Science and Technology, 2002 (2): 42-44

[3]曾丁,黃文彬,華云龍. 筒倉壁壓的有限元分析[J].農業工程學報, 1998, 14 (2): 44-48

Zeng Ding, Huang Wenbin, Hua Yunlong. Finite element analysis of silo wall pressure [J]. Journal of Agricultural Engineering, 1998, 14 (2): 44-48

[4]石翠霞. 筒倉內小麥、玉米堆的壓縮特性、倉壁壓應力及儲糧總重量的研究[D].南京:南京財經大學, 2011

Shi Cuixia. Compression characteristics of wheat and corn pile，compressive stress of silo wall and total weight of stored grain in silo [D]. Nanjing: Nanjing University of Finances and Economics, 2011

[5]嚴曉婕.稻谷堆與籽粒的壓縮特性及其在筒倉中密度與應力分布的研究[D].南京:南京財經大學, 2014

Yan Xiaojie. Compression characteristics of rice paddy and grain and their density and stress distribution in a silo [D]. Nanjing: Nanjing University of Finances and Economics, 2014

[6] Thompson R A, Isaacs G W. Porosity determination of grains and seeds with an air-comparision pycnometer [J]. Trans. ASAE, 1967, 10(2):693-696

[7] Gustafson R J. Density and porosity changes of shelled corn during drying[J]. Trans. ASAE, 1972, 15(2):523-527

[8] Chang C S. Measuring Density and Porosity of Grain Kernels Using a Gas Pycnometer [J]. Cereal Chemistry.1988,65(1):13-15

[9] 張淑珍，巴衛國，徐寶江，等． 散粒物料孔隙率測定裝置的研制[J]. 農業工程學報，1996,12(1):187-191

Zhang Shuzhen, Ba Weiguo, Xu baojiang, et al. granular material porosity determination device [J]. Journal of Agricultural Engineering, 1996,12 (1): 187-191

[10] Marinder B O. A simple apparatus for determining the density of solids[J]. Measurement Science & Technology,1996,7(11) :1569-1573．

[11] Bielders C L，Debacker L W，Delvaux B. Particle density of volcanic soils as measured with a gas pycnometer[J]. Soil Science,1990,54(3) : 822-826

[12] 郝倩，陳弘毅，平學文，等. 一種新型糧食孔隙度就倉測量系統研究[J]. 農機化研究，2015,9:180-184

Hao Qian, Chen Hongyi,Ping Xuewen, et al. Research on a new grain porosity measurement system [J]. Agricultural Mechanization Research, 2015,9:180-184

[13] Cheng J, Sun W. Automatic measurement of pores and porosity in pork ham and their correlations with processing time,water content and texture[J]. Meat science,2006,72(2) : 294-302

[14] Mang J T, Hjelm R P, Francois E G. Measurement of porosity in a composite high explosive as a function of pressing conditions by ultra-small-angle neutron scattering with contrast variation[J]. Propellants, Explosives,Pyrotechnics,2010,35(1): 7-14

[15] Yao Y, Liu D. Comparison of low- field NMR and mercury intrusion porosimetry in characterizing pore size distributions of coals[J]. Fuel,2012,95: 152-158

[16] Nippolainen E, Tuomas E, Ketolainen J, et al. Measuring tablet porosity using multispectral imaging system[J]. Optical Review,2010,17(3) : 323- 326

[17] Podymova N, Karabutov A, Kobeleva L, et al. Laser optoacoustic method of local porosity measurement of particles reinforced composites[J]. Journal of Physics Conference Series,2011, 278(1) : 12-38

[18] Hogekamp S, Pohl M. Porosity measurement of fragile agglomerates[J]. Powder Technology, 2003,130(1-3):385-392

[19] Umnova O, Chul H, Cummings A. Deduction of tortuosity and porosity from acoustic reflection and transmission measurements on thick samples of rigid-porous materials[J]. Applied Acoustics,2005,66(6):607-624.

Research on Porosity Distribution of Paddy in Silos

Tang Fuyuan Xu Qian Cheng Xuduo
(College of Food Science and Engineering, Nanjing University of Finance and Economics/Collaborative Innovation Center for Modern Grain Circulation and Safety, Nanjing 210046)

The LHT-1 rebound modulus tester was used to determine the compression density of paddy(South japonica5055). The differential equations of density，pressure and depth in a silo was built. Numerical method was used to obtain the relationship of density and depth. The porosity of surface layer of paddy (no compression) was measured by grain porosity meter. According to the porosity of surface layer, the density of surface layer and densities calculated in the silo, the relation between grain porosity and grain depth in silo was calculated. Calculation results show that in the 30 meter cylinder part of a silo with a diameter of 20 meters, the change range of the porosity of South japonica5055 paddy is 61.00%～56.32%, and in the cone hopper of 10 meters of the silo, the change range of porosity is 56.32%～59.77%; in a silo with a cone hopper, porosity of paddy in cylinder of silo decreased with the paddy depth, the porosity of paddy in cone hopper increased with the increase of depth. At the same depth, the porosity of paddy decreases with the increase of silo diameter.

silo, paddy, porosity, density

TS210.4

A

1003-0174(2017)12-0110-07

時間：2017-12-20 17:02:53

網絡出版地址：http://kns.cnki.net/kcms/detail/11.2864.TS.20171220.1702.002.html

國家自然科學基金(31371856)，江蘇高校優勢學科建設工程(PAPD)

2017-06-22

唐福元,男,1962年出生,副教授,糧食儲藏工程

程緒鐸,男,1957年出生,教授,糧食儲藏工程