陳鵬梟，朱文學，吳建章，楊 帆，楊德勇，黃 凱

1.河南工業(yè)大學 糧油食品學院，河南 鄭州 450001

2.中國農業(yè)大學 工學院， 北京 100083

3.南京農業(yè)大學 工學院， 江蘇 南京 210095

倉儲糧堆干燥過程涉及糧粒內濕熱傳遞、孔隙之間的氣體流動以及孔隙間氣體與糧粒之間的濕熱傳遞，其研究涉及生物、化學、熱力學以及流體力學等多種學科領域，加之干燥過程中糧堆內各物理參數較難測定，直接對濕熱傳遞機制進行研究十分困難。隨著計算機技術的發(fā)展，數值模擬方法成為研究糧食干燥過程濕熱傳遞機制的重要手段[1]。

基于不同的模型假設，目前關于倉儲糧堆濕熱傳遞的模擬模型可分為兩種。

一種是基于連續(xù)介質假設的倉儲糧堆濕熱傳遞模型[2-5]。Hammami等[6]對圓筒倉糧堆進行了濕熱傳遞模擬，研究發(fā)現通風溫度對谷物溫度變化具有重要作用，并總結了氣流特性對倉內糧食溫濕度的影響。Wang等[7]建立了描述深層糧堆內部濕熱傳遞的數學模型，研究發(fā)現在冷卻干燥通風過程中，降溫和降水同時存在，冷卻前沿的移動速度大大快于干燥前沿，且糧堆中溫度和水分分布具有分層的現象。Chen等[8-9]以小麥為研究對象，采用局部非平衡法，建立了倉儲小麥三維通風干燥濕熱傳遞模型，通過模擬研究發(fā)現固定床干燥有明顯的干燥前沿，干燥不均勻度較大，且干燥前沿附近的溫度梯度最大。

另一種是基于離散介質假設的孔道網絡濕熱傳遞模型[10-13]，其物理模型如圖1所示。Kharaghani等[14-16]通過X射線獲取球形堆積體三維圖片，提取該堆積結構的孔道網絡信息，并基于該孔道網絡信息模擬了物體液相傳遞以及干燥過程，研究發(fā)現，在第一個干燥階段，模擬所得干燥曲線與實驗結果基本吻合，之后實驗曲線發(fā)生明顯干燥梯度，而模擬結果則表現為無梯度。Yuan等[17-18]以玉米為研究對象，構建了單層籽粒堆積孔道網絡物理模型，并結合多尺度濕熱傳遞數學模型對倉儲糧堆干燥過程進行了模擬，研究發(fā)現物料孔隙率、配位數越大，物料干燥越快，升溫也越快，達到相同濕含量所需干燥時間越短。黃凱[19]提出了構建三維非規(guī)則孔道網絡物理模型的方法，并利用三維規(guī)則孔道網絡模擬了糧堆通風干燥過程，研究發(fā)現任意時刻糧堆孔隙內的氣相分壓隨糧層深度的降低而降低，隨著通風的進行，玉米顆粒平均含水率不斷降低，溫度先迅速降低后緩慢升高，內部不同層含水率由外而內逐步降低。Chen等[20-21]構建了倉儲糧堆二維非規(guī)則以及三維非規(guī)則孔道網絡模型，并結合孔道網絡濕熱傳遞數學模型，對倉儲糧堆通風干燥過程進行了模擬，經干燥實驗驗證，兩種孔道網絡模型皆能夠描述糧堆機械通風干燥過程，且研究結果表明孔道網絡模型中的顆粒配孔數與拓撲結構對模型模擬結果影響較大。

從前述分析可以得出，基于離散介質假設的非規(guī)則孔道網絡物理模型能夠清晰描述倉儲糧堆內部各顆粒之間、孔喉之間以及顆粒與孔喉之間的空間位置關系，這在探究糧食濕熱傳遞機理與三傳過程模擬領域具有較大優(yōu)勢。作者擬以倉儲糧堆為研究對象，利用三維CT掃描設備與Avizo軟件對顆粒堆積樣件的孔喉進行提取，通過EDEM軟件模擬糧食顆粒堆積狀態(tài)，基于泰森多邊形原理，通過voro++軟件處理堆積體中的顆粒坐標及半徑等數據，制定孔隙劃分規(guī)則，對不合理孔隙進行處理，結合Avizo軟件獲取的孔喉信息，獲得能夠描述倉儲糧堆內部空間結構的三維非規(guī)則孔道網絡物理模型。

圖1 糧堆孔道網絡模型構建過程圖

1 實際顆粒堆積結構孔道網絡的提取

1.1 試驗材料與方法

試驗材料：顆粒半徑為4 mm的玻璃珠若干，內邊長為3 cm的正方體容器。

試驗設備：英華檢測(上海)有限公司的phoenix v/ tome/ x m全能型X射線微焦點CT系統(tǒng)。

試驗方法: 將備好的玻璃珠裝滿內邊長為3 cm的容器，所有顆粒完全在容器內部，對樣品進行掃描，采集原始圖像，得到反映樣件內部結構的圖片和數據，使用Avizo軟件對所得樣件內部孔喉結構參數進行分析。

1.2 試驗結果與分析

圖2為玻璃球形顆粒堆積結構的三維CT掃描圖片。從圖2可以看出，該CT從3個不同的維度對物料進行步進掃描，得到關于物料結構的三維立體圖像。

經過無損掃描以及Avizo軟件對樣件孔道網絡的孔喉提取發(fā)現，該試驗樣品中共含有51個顆粒、66個孔隙、205個喉道，其平均配位數為6.246，孔隙率為49.36%，孔隙密度為2.4×106個/m3。圖3為球形顆粒堆積結構內部喉道長度分布，從圖3可得，該樣品的喉道長度主要分布在6.0～14.0 mm之間，大約占總喉道數的96%，喉道長度分布大致為Poisson分布。圖4為球形顆粒堆積結構的孔隙半徑分布，從圖4可得，該結構孔隙半徑多在2.0～10.0 mm之間。

圖2 球形顆粒堆積結構局部三維掃描圖

圖3 球形顆粒堆積結構喉道長度分布

圖4 球形顆粒堆積結構孔隙半徑分布

2 倉儲糧堆三維非規(guī)則孔道網絡物理模型的構建

2.1 基于泰森多邊形的孔道網絡物理模型的構建

利用EDEM軟件模擬顆粒堆積狀態(tài)，假設容器形狀為方倉，鋁合金材質，玉米顆粒為球體。顆粒、糧倉倉壁以及顆粒與顆粒之間的物性參數，如表1所示。通過EDEM軟件模擬得到堆積體中玉米顆粒的信息，圖5所示為經過voro++軟件處理過的顆粒堆積結構對應的孔道網絡劃分。通過上述方法構建的物理模型大小為-0.020 12≤x≤0.020 12、-0.020 12≤y≤0.020 12、0≤z≤ 0.063 36 (單位: m)，總共有179個顆粒，孔隙數為991，喉道數為1 978，顆粒半徑為4.0 mm，孔隙配位數主要集中在4。經計算，該模型的孔隙率為52.39%，孔喉連通熵為6.775 4，孔隙密度為1.03×107個/m3。與樣件提取所得的參數對比，從規(guī)模數上講，提取所得的孔道網絡物理模型的顆粒數、孔隙數與喉道數之比為1∶1.275∶3.980，而構建所得物理模型的三參數之比為1∶5.536∶11.050，相關參數比值相差較大；提取所得的三維非規(guī)則孔道網絡物理模型的平均配位數為6.246，構建所得三維非規(guī)則孔道網絡物理模型平均配位數為3.992，相差較大；此處構建獲得的孔隙密度為提取孔隙密度的4.3倍左右，相差較大。

表1 模擬所用材料物性參數

圖5 顆粒堆積結構孔道網絡劃分

圖6為構建的三維非規(guī)則孔道網絡模型中孔間距的長度分布規(guī)律，從圖6可得，孔間距中有較大一部分長度小于5 mm，孔間距等于喉道長度加上喉道兩端的兩孔隙半徑，依此定義可知，小于5 mm的喉道個數會更多，這與提取結果相差較大。

圖6 三維非規(guī)則孔道網絡模型孔間距長度分布

在等徑球狀顆粒最緊密堆積的多孔介質的結構中，等徑球狀顆粒(半徑為R)多孔介質模型最小單元體由 4 個球組成，呈正四面體結構。其最小單元體的孔隙半徑(r)，如式(1)所示。

(1)

如有喉道長度小于2r，則該喉道不符合實際情況。經統(tǒng)計，構建的三維非規(guī)則孔道網絡模型中孔間距小于2r的有191個。因此,需對該三維非規(guī)則孔道網絡物理模型進行改進。

2.2 倉儲糧堆多尺度物理模型的改進

將孔道網絡中不符合實際情況的孔隙按照一定規(guī)則進行合并，使其更符合真實孔道網絡結構。

2.2.1 孔隙合并規(guī)則

(1)合并規(guī)則一。在孔道網絡模型中，孔心距小于2r的兩孔隙需要合并，形成的新孔隙的孔心坐標為此兩孔中心連線的中點。

(2)合并規(guī)則二。如圖7所示，圖7a為合并孔隙前的孔喉結構示意圖，孔1的相鄰顆粒為C1和C3，孔2的相鄰顆粒為C1和C2，孔3的相鄰顆粒為C1、C2、C3和C4，孔1、2、3兩兩通過不同的喉道相連接，此時可以認為孔1、2、3都在顆粒C1、C2、C3、C44個顆粒圍成的大孔隙內部，可以將其合并成一個大孔隙，為了簡化算法，直接將孔1和孔2合并到孔3，將該孔隙拓撲結構改變成圖7b，即完成該類孔隙的合并。

2.2.2 孔隙半徑提取算法

平均半徑算法是指將孔隙周圍顆粒中心與孔心的距離減去對應顆粒半徑的值的平均值視為孔隙半徑的一種算法。以一個配位數為3的孔隙為例，孔隙周圍有3個顆粒，3個顆粒的半徑分別為r1、r2、r3，從孔隙中心到3個顆粒中心的長度分別為l1、l2、l3，則平均距離法得到孔隙半徑rp。

(2)

以平均半徑算法獲得各個孔隙的半徑之后，合并孔間距小于兩孔半徑之和的孔隙，新得孔隙的中心為兩合并孔隙中心連線的中點，再應用平均半徑算法，獲得新的孔隙半徑。

2.3 改進后孔道網絡模型結構參數分析

算法改進后，獲得的三維非規(guī)則孔道網絡物理模型顆粒數為179，孔隙數為237，喉道數為702，顆粒半徑為4.0 mm，平均配位數為5.924。將上述參數與改進前的模型參數對比發(fā)現，提取所得的孔道網絡物理模型的顆粒數、孔隙數與喉道數之比為1∶1.275∶3.980，而改進后物理模型的3個參數之比為1∶1.324∶3.922，單從模型規(guī)模數上講，兩物理模型相差較小，可以互相表征。圖8為改進后三維非規(guī)則孔道網絡模型配位數分布，從圖8可得，配位數主要集中在3、4、5和6，大致呈Poisson分布。圖9為改進后三維非規(guī)則孔道網絡模型喉道長度分布，可知該孔道網絡模型的喉道長度主要分布在3.0～13.0 mm之間，大致呈Poisson分布。

圖8 改進后三維非規(guī)則孔道網絡模型配位數分布

圖9 改進后三維非規(guī)則孔道網絡模型喉道長度分布

圖10為改進后三維非規(guī)則孔道網絡模型孔隙半徑分布，從圖10可知，在該模型中孔隙分布大致在2.0～10.0 mm之間，且存在直徑大于10 mm的孔隙，這是由于在孔道網絡模型的邊界處，孔隙只有一側有顆粒，限制孔隙半徑的條件比倉儲糧堆內部的孔隙少，所以此處孔隙半徑較大；根據提取信息可得該孔道網絡模型的孔喉連通熵為5.712 4。

圖11、圖12為改進后三維非規(guī)則孔道網絡模型顆粒配孔數、配喉數分布，顆粒配孔數在7～18之間，顆粒配喉數在4～20之間，其中大部分配孔數集中在8～16之間。

圖10 改進后三維非規(guī)則孔道網絡模型孔隙半徑分布

圖11 改進后三維非規(guī)則孔道網絡模型顆粒配孔數分布

圖12 改進后三維非規(guī)則孔道網絡模型顆粒配喉數分布

3 實際結構與構建模型孔喉數據對比分析與驗證

3.1 實際結構與構建模型孔喉數據對比

將樣件提取數據與三維非規(guī)則孔道網絡模型構建數據進行對比分析，如表2所示。從表2可以看出，構建的三維非規(guī)則孔道網絡模型結構參數與提取的結果相符合，可以說明該模型在統(tǒng)計意義上與實際顆粒堆積體的各項孔隙參數都較為吻合，且構建的孔道網絡模型結構參數中有顆粒配孔數與顆粒配喉數，能夠反映堆積結構中顆粒與孔、喉之間的拓撲關系，說明構建的三維非規(guī)則孔道網絡模型要比Avizo軟件提取所得的三維非規(guī)則孔道網絡模型更能反映顆粒堆積結構內部的孔隙結構真實情況。

表 2 三維提取信息與三維構建信息對比

3.2 構建所得孔道網絡物理模型流通性驗證

圖13 三維非規(guī)則孔道網絡模型初始風壓與出口流速關系

在糧食通風干燥過程中，空氣由進風口進入，穿過糧食區(qū)域，從出風口排出，這個過程可以看作是一個流體穿過非固結多孔介質內部的過程，采用達西定律來描述這個過程[20-21]，并利用程序計算初始風壓與出口流速之間的關系。圖13為三維非規(guī)則孔道網絡模型初始風壓與出口流速關系，從圖13可以看出，底部風壓與出口風速大致呈線性變化。該結果可以證明兩個結論：第一，證明了該物理模型能夠有效表征非固結多孔介質的孔隙連通性；第二，證明了該物理模型可以結合濕熱傳遞數學模型用來描述非固結多孔介質內部流體流速分布情況。

4 結論與展望

采用孔道網絡方法構建了能夠描述倉儲糧堆內部結構的孔道網絡模型，并運用達西定律驗證了物理模型內部孔喉結構的流通性，通過分析研究，可得結論如下。

利用三維CT掃描設備與Avizo軟件提取所得顆粒堆積結構中喉道長度主要分布在6.0～14.0 mm之間，喉道長度分布大致為Poisson分布；孔隙半徑多在2.0～10.0 mm之間。

通過EDEM軟件模擬糧食顆粒堆積狀態(tài)，基于泰森多邊形原理，通過voro++軟件處理堆積體中的顆粒坐標及半徑等數據，結合孔隙劃分規(guī)則對不合理孔隙進行處理，所得的顆粒堆積結構參數中配位數主要集中在3、4、5和6，喉道長度主要分布在3.0～13.0 mm之間，孔隙半徑分布大致在2.0～10.0 mm之間，3個參數大致均呈Poisson分布；顆粒配孔數在7～18之間，顆粒配喉數在4～20之間，其中大部分配孔數集中在8～16之間。

通過構建所得三維非規(guī)則孔道網絡模型較Avizo軟件所獲得的孔隙結構參數更全面，更能準確地描述倉儲糧堆內部顆粒、孔隙以及喉道之間的關系。