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

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

2.中國農(nóng)業(yè)大學(xué) 工學(xué)院， 北京 100083

3.南京農(nóng)業(yè)大學(xué) 工學(xué)院， 江蘇 南京 210095

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

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

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

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

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

圖1 糧堆孔道網(wǎng)絡(luò)模型構(gòu)建過程圖

1 實際顆粒堆積結(jié)構(gòu)孔道網(wǎng)絡(luò)的提取

1.1 試驗材料與方法

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

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

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

1.2 試驗結(jié)果與分析

圖2為玻璃球形顆粒堆積結(jié)構(gòu)的三維CT掃描圖片。從圖2可以看出，該CT從3個不同的維度對物料進(jìn)行步進(jìn)掃描，得到關(guān)于物料結(jié)構(gòu)的三維立體圖像。

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

圖2 球形顆粒堆積結(jié)構(gòu)局部三維掃描圖

圖3 球形顆粒堆積結(jié)構(gòu)喉道長度分布

圖4 球形顆粒堆積結(jié)構(gòu)孔隙半徑分布

2 倉儲糧堆三維非規(guī)則孔道網(wǎng)絡(luò)物理模型的構(gòu)建

2.1 基于泰森多邊形的孔道網(wǎng)絡(luò)物理模型的構(gòu)建

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

表1 模擬所用材料物性參數(shù)

圖5 顆粒堆積結(jié)構(gòu)孔道網(wǎng)絡(luò)劃分

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

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

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

(1)

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

2.2 倉儲糧堆多尺度物理模型的改進(jìn)

將孔道網(wǎng)絡(luò)中不符合實際情況的孔隙按照一定規(guī)則進(jìn)行合并，使其更符合真實孔道網(wǎng)絡(luò)結(jié)構(gòu)。

2.2.1 孔隙合并規(guī)則

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

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

2.2.2 孔隙半徑提取算法

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

(2)

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

2.3 改進(jìn)后孔道網(wǎng)絡(luò)模型結(jié)構(gòu)參數(shù)分析

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

圖8 改進(jìn)后三維非規(guī)則孔道網(wǎng)絡(luò)模型配位數(shù)分布

圖9 改進(jìn)后三維非規(guī)則孔道網(wǎng)絡(luò)模型喉道長度分布

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

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

圖10 改進(jìn)后三維非規(guī)則孔道網(wǎng)絡(luò)模型孔隙半徑分布

圖11 改進(jìn)后三維非規(guī)則孔道網(wǎng)絡(luò)模型顆粒配孔數(shù)分布

圖12 改進(jìn)后三維非規(guī)則孔道網(wǎng)絡(luò)模型顆粒配喉數(shù)分布

3 實際結(jié)構(gòu)與構(gòu)建模型孔喉數(shù)據(jù)對比分析與驗證

3.1 實際結(jié)構(gòu)與構(gòu)建模型孔喉數(shù)據(jù)對比

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

表 2 三維提取信息與三維構(gòu)建信息對比

3.2 構(gòu)建所得孔道網(wǎng)絡(luò)物理模型流通性驗證

圖13 三維非規(guī)則孔道網(wǎng)絡(luò)模型初始風(fēng)壓與出口流速關(guān)系

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

4 結(jié)論與展望

采用孔道網(wǎng)絡(luò)方法構(gòu)建了能夠描述倉儲糧堆內(nèi)部結(jié)構(gòu)的孔道網(wǎng)絡(luò)模型，并運用達(dá)西定律驗證了物理模型內(nèi)部孔喉結(jié)構(gòu)的流通性，通過分析研究，可得結(jié)論如下。

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

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

通過構(gòu)建所得三維非規(guī)則孔道網(wǎng)絡(luò)模型較Avizo軟件所獲得的孔隙結(jié)構(gòu)參數(shù)更全面，更能準(zhǔn)確地描述倉儲糧堆內(nèi)部顆粒、孔隙以及喉道之間的關(guān)系。