陳清華蘇國(guó)用孫美華姜闊勝劉 萍

(1. 安徽理工大學(xué)機(jī)械工程學(xué)院，安徽 淮南 232001；2. 智能礦山技術(shù)與裝備安徽省重點(diǎn)實(shí)驗(yàn)室，安徽 淮南 232001)

中國(guó)是農(nóng)業(yè)和糧食消費(fèi)大國(guó)，糧食儲(chǔ)存對(duì)保障國(guó)家安全和社會(huì)穩(wěn)定意義重大，如何實(shí)現(xiàn)糧食的安全儲(chǔ)藏一直是研究熱點(diǎn)之一[1-2]。研究[3]表明冷卻干燥通風(fēng)可以控制糧堆溫濕度，進(jìn)而實(shí)現(xiàn)糧食的長(zhǎng)期安全存放，而糧食作為典型的多孔介質(zhì)，干燥通風(fēng)過(guò)程中其內(nèi)部的傳熱傳質(zhì)規(guī)律較為復(fù)雜。導(dǎo)熱系數(shù)、熱擴(kuò)散率等作為糧食的重要熱物性參數(shù)，由于可直接表征糧食傳遞熱量的能力[4-5]，其準(zhǔn)確獲取對(duì)于指導(dǎo)和優(yōu)化糧食干燥和安全儲(chǔ)藏工藝至關(guān)重要。由于多孔介質(zhì)傳熱機(jī)理較為復(fù)雜，試驗(yàn)測(cè)試仍是獲取其熱物性參數(shù)的主要手段[6]。例如谷和平等[7]和龔紅菊等[8]分別運(yùn)用正規(guī)狀況法及Dickerson圓桶瞬態(tài)熱流法，對(duì)顆粒狀堆積物、稻谷的導(dǎo)熱系數(shù)和熱擴(kuò)散率進(jìn)行了試驗(yàn)測(cè)試。但現(xiàn)有測(cè)試方法中，較少的對(duì)試驗(yàn)過(guò)程中糧堆內(nèi)部熱量的側(cè)向流動(dòng)進(jìn)行考慮，從而使測(cè)試精度受到一定的影響。文獻(xiàn)[9]提出基于非穩(wěn)態(tài)導(dǎo)熱乘積法測(cè)顆粒狀材料熱擴(kuò)散系數(shù)，優(yōu)點(diǎn)是消除了側(cè)向熱流誤差帶來(lái)的影響，因而更適用于松散物料熱物性測(cè)試。在此基礎(chǔ)上，岳高偉等[10]利用高低溫交變濕熱試驗(yàn)箱設(shè)定溫度邊界(-50～100 ℃)，將應(yīng)用范圍擴(kuò)展至低溫環(huán)境。由于非穩(wěn)態(tài)導(dǎo)熱乘積法模型只能直接測(cè)算熱擴(kuò)散系數(shù)a，如需獲取比熱容Cp和導(dǎo)熱系數(shù)λ，仍需借助磁力攪拌水卡計(jì)[11]等測(cè)試手段，大大限制了其推廣應(yīng)用性。通過(guò)易于獲取的數(shù)據(jù)信息結(jié)合明確的數(shù)學(xué)模型，反演或估計(jì)被研究對(duì)象的多個(gè)參數(shù)被稱為多宗量傳熱反問(wèn)題[12]，目前已在多個(gè)領(lǐng)域得到應(yīng)用并被認(rèn)可。本試驗(yàn)擬基于多維非穩(wěn)態(tài)導(dǎo)熱乘積法數(shù)學(xué)模型，并結(jié)合參數(shù)估計(jì)法[13-14]對(duì)糧食顆粒的導(dǎo)熱系數(shù)λ、比熱容Cp等進(jìn)行反演計(jì)算，本測(cè)試方法試驗(yàn)操作難度較低，便于獲取高效準(zhǔn)確的糧粒熱物性參數(shù)，進(jìn)而為制定糧食儲(chǔ)藏技術(shù)方案提供科學(xué)依據(jù)。

1 數(shù)學(xué)模型

在某一刻τ將一直徑為D，高為δ，上端面絕熱，初始溫度為t0的短圓柱體，放入溫度為tw的恒溫環(huán)境中。設(shè)其溫度分布與角度無(wú)關(guān)，則其導(dǎo)熱微分方程見(jiàn)式(1)[15]。

(1)

式中：

T——溫度，℃；

r——徑向變量，m；

x——軸向變量，m；

α——熱擴(kuò)散系數(shù)，m2/s。

顯然該問(wèn)題為二維非穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題，一般采用數(shù)值解法，但其步驟復(fù)雜且不易編程實(shí)現(xiàn)。介于此，試驗(yàn)中引用多維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題的乘積解法[15]。如圖1所示，一短圓柱體由直徑為r的無(wú)限長(zhǎng)圓柱和一塊厚度為x的無(wú)限大平板垂直相貫而成，則有：

T(R,x,τ)=[T(x,τ)]×[T(r,τ)]。

(2)

式(1)的解析解：

(3)

式(3)中包含無(wú)窮級(jí)數(shù)，顯然難以直接求解，根據(jù)試算取級(jí)數(shù)前6項(xiàng)即可滿足精度要求，式(3)可改寫(xiě)為：

(4)

式中：

θ(r,x,τ)——短圓柱過(guò)余溫度，℃；

θ(x,τ)——無(wú)限大平板過(guò)余溫度，℃；

θ(r,τ)——無(wú)限長(zhǎng)圓柱過(guò)余溫度，℃；

θ0——初始過(guò)余溫度，θ=t-tw，其中t為任意時(shí)刻試樣溫度，℃；

tw——恒溫箱內(nèi)部溫度，℃；

L——無(wú)限大平板厚度，m；

R——無(wú)限長(zhǎng)圓柱半徑，m；

μn——特征值；

F0——傅里葉數(shù)，F(xiàn)0=ατ/y2，其中α為熱擴(kuò)散系數(shù)(m2/s)；y為特征長(zhǎng)度；對(duì)于無(wú)限大平板y=L，對(duì)于無(wú)限長(zhǎng)圓柱體y=R。

根據(jù)熱擴(kuò)散率公式a=λ/(ρCp)，顯然如果已知試樣密度ρ、導(dǎo)熱系數(shù)λ、比熱容Cp及厚度δ等物理參數(shù)，可計(jì)算得到短圓柱體內(nèi)任意一點(diǎn)x處，τ時(shí)刻的無(wú)量綱過(guò)余溫度θ(r,x,τ)，此為傳熱正問(wèn)題。

圖1 非穩(wěn)態(tài)導(dǎo)熱乘積法原理圖

2 熱物性參數(shù)估計(jì)靈敏度分析

為研究試樣內(nèi)任意一點(diǎn)處的溫度Y(τ,η)變化規(guī)律，某一時(shí)刻τ將試樣放入溫度為tw的環(huán)境中，溫度擾動(dòng)記為u(τ)。考慮試樣的導(dǎo)熱系數(shù)λ、比熱容Cp等熱物性參數(shù)，以及環(huán)境溫度tw、時(shí)間τ、空間位置等因素都會(huì)影響u(τ)的值，故將這些參數(shù)構(gòu)成一個(gè)向量：

η=(η1,η2,η3,……,ηm)T，

(5)

式中：

m——參數(shù)個(gè)數(shù)。

顯然，在試驗(yàn)過(guò)程中對(duì)于一個(gè)確定的溫度測(cè)點(diǎn)，式(5)中未知參數(shù)η1，η2為試樣的λ和Cp。取試樣中同一高度處，距中心軸均為r的k個(gè)離散測(cè)點(diǎn)處的溫度測(cè)量值記為Yi(τ,η)(i=1,……,k)，然后給出一組λ和Cp的估計(jì)值，根據(jù)式(4)計(jì)算出的過(guò)余溫度記為θ(τ,η)，各離散測(cè)點(diǎn)i的過(guò)余溫度記為θi(τ,η)(i=1，……，k)，并進(jìn)行對(duì)比：

(6)

由于測(cè)量誤差以及參數(shù)估計(jì)具有隨機(jī)性，Yi(τ,η)與θi(τ,η)之間必然存在偏差，即ε(η)>0。為此，參數(shù)估計(jì)的目標(biāo)是基于式(4)迭代，使偏差ε(η)→min。參數(shù)向量η中的λ和Cp值要能同時(shí)以足夠的精度估計(jì)出來(lái)，必須滿足在最小二乘估計(jì)意義下，在參數(shù)估計(jì)的時(shí)間區(qū)間內(nèi)，參數(shù)的靈敏度線性無(wú)關(guān)。本試驗(yàn)中靈敏度是指θ(τ,η)對(duì)參數(shù)η的一階偏導(dǎo)。同時(shí)，考慮到方程式的復(fù)雜性，直接求偏導(dǎo)數(shù)較為困難，采用二階中心差商進(jìn)行計(jì)算：

Xij=[θ(τi,η1,……,ηj+△ηj,……,ηm)-θ(τi,η1,……,ηj-△ηj,……,ηm)]/(2△ηj)，

(7)

式中：

Xij——τi時(shí)刻θ(τi,η)對(duì)參數(shù)ηj的靈敏度系數(shù)。

通常取△ηj=0.000 1ηj。

只有在λ和Cp線性無(wú)關(guān)時(shí)，2個(gè)參數(shù)才能同時(shí)估計(jì)，因此在估計(jì)之前進(jìn)行靈敏度相關(guān)性分析。分別選取玉米和稻谷進(jìn)行研究。根據(jù)參考文獻(xiàn)[16]和[17]分別選定玉米的相關(guān)參數(shù)為：密度750 kg/m3，導(dǎo)熱系0.176 5 W/(m·℃)，比熱容1 500 J/(kg·℃)，含水率17%；稻谷參數(shù)為：密度1 480 kg/m3，導(dǎo)熱系數(shù)0.099 5 W/(m·℃)，比熱容911.9 J/(kg·℃)。

由圖2、3可以看出，玉米和稻谷的導(dǎo)熱系數(shù)和比熱容參數(shù)靈敏度線性均不相關(guān)，但比熱容Cp靈敏度均非常小(僅為10-3數(shù)量級(jí))，若直接進(jìn)行參數(shù)估計(jì)，反演得到的比熱容與實(shí)際相差較大[14]。

圖2 玉米試樣導(dǎo)熱系數(shù)及比熱容靈敏度值

圖3 稻谷試樣導(dǎo)熱系數(shù)及比熱容靈敏度值

3 關(guān)鍵技術(shù)

基于多維非穩(wěn)態(tài)導(dǎo)熱乘積法原理，設(shè)計(jì)測(cè)糧食熱物性參數(shù)的測(cè)試系統(tǒng)，需要解決的關(guān)鍵技術(shù)有測(cè)試系統(tǒng)構(gòu)建、多維非穩(wěn)態(tài)傳熱條件、參數(shù)估計(jì)算法等。

3.1 測(cè)試系統(tǒng)設(shè)計(jì)

如圖4所示，系統(tǒng)主要由恒溫箱、信號(hào)采集與傳輸模塊、信號(hào)發(fā)射端和接收端、試樣筒、溫度傳感器和微型計(jì)算機(jī)等構(gòu)成。為了避免在恒溫箱上開(kāi)孔影響恒溫箱效果，以及增加試驗(yàn)操作難度，系統(tǒng)采用深圳蜂匯公司型號(hào)為Z-0004的ZigBee無(wú)線數(shù)據(jù)采集模塊實(shí)現(xiàn)信號(hào)采集與傳輸。恒溫箱溫度范圍為0～100 ℃，控溫精度±0.1 ℃；測(cè)溫?zé)犭娕季取?.2 ℃；ZigBee數(shù)據(jù)采集模塊的射頻芯片CC2530的接收靈敏度為-97 dBm，有效保證了溫度的實(shí)時(shí)傳輸。系統(tǒng)軟件通過(guò)LabView編程建立基本松散物體參數(shù)輸入?yún)^(qū)、采集溫度變化顯示區(qū)、控制按鈕區(qū)和計(jì)算結(jié)果區(qū)等，便于對(duì)測(cè)試結(jié)果的觀察以及對(duì)不同糧食顆粒測(cè)算參數(shù)設(shè)置。

1. 采集信號(hào)接收端 2. 微型計(jì)算機(jī) 3. 信號(hào)采集與傳輸單元 4. 絕熱端蓋 5. 恒溫箱 6. 試樣 7. 試樣筒體 8. 溫度測(cè)點(diǎn) 9. 采集信號(hào)發(fā)射端

圖4 測(cè)試裝置簡(jiǎn)圖

Figure 4 Diagram of testing system

3.2 非穩(wěn)態(tài)導(dǎo)熱乘積物理邊界的實(shí)現(xiàn)

根據(jù)非穩(wěn)態(tài)導(dǎo)熱乘積法原理及物理模型要求，試樣為圓柱形，同時(shí)上表面絕熱，下表面和圓柱面為恒溫邊界。此處利用恒溫箱實(shí)現(xiàn)恒溫環(huán)境，并設(shè)計(jì)圓柱形試樣筒盛放試樣，為保證試樣放入恒溫箱后較短時(shí)間內(nèi)與環(huán)境溫度保持一致，筒體采用黃銅制作，在內(nèi)壁設(shè)置溫度測(cè)點(diǎn)，待溫度與環(huán)境溫度一樣時(shí)開(kāi)始試驗(yàn)測(cè)試，以此實(shí)現(xiàn)圓柱面和下表面的恒溫邊界條件。同時(shí)，試樣筒上端蓋為絕熱棉制作，形成試樣上絕熱表面。

3.3 比熱容Cp靈敏度較低問(wèn)題的處理

針對(duì)試驗(yàn)選用的材料Cp靈敏度較低，比熱容不可直接估算，采用隨機(jī)共軛梯度法進(jìn)行求解[18]，由于式(4)中的F0=aτ/y2，只有一個(gè)未知數(shù)a，因此，首先反演估算a，然后利用a對(duì)估計(jì)的Cp值進(jìn)行修正，具體流程(見(jiàn)圖5)：

① 反演估計(jì)a；

② 估計(jì)材料λ及Cp的原始估計(jì)值；

③ 根據(jù)式(4)、(5)、(6)求得滿足精度要求的a及λ；

④ 根據(jù)a估計(jì)值對(duì)原始估計(jì)值Cp進(jìn)行修正，得到最終值Cp*。

為避免在測(cè)試過(guò)程中，因外界干擾因素造成的個(gè)別數(shù)據(jù)波動(dòng)對(duì)估算結(jié)果造成影響，本試驗(yàn)首先計(jì)算出每個(gè)特定時(shí)刻τ對(duì)應(yīng)的一組數(shù)值λ和Cp*，然后對(duì)計(jì)算出的一系列物性參數(shù)求加權(quán)平均作為最終參數(shù)的估計(jì)值。

圖5 參數(shù)估計(jì)流程圖

3.4 溫度測(cè)試準(zhǔn)確性及可靠性

當(dāng)材料質(zhì)地均勻，滿足各項(xiàng)同性的性質(zhì)時(shí)，溫度測(cè)點(diǎn)的布置對(duì)物性參數(shù)測(cè)算結(jié)果的影響不大。因此，試驗(yàn)前對(duì)稻谷顆粒進(jìn)行篩選，裝入箱體時(shí)進(jìn)行壓實(shí)處理。同時(shí)，在離圓柱側(cè)面和底面相同距離的不同位置處，安裝多個(gè)溫度傳感器，并取加權(quán)平均作為最終測(cè)試結(jié)果。此外，系統(tǒng)采取無(wú)線信號(hào)發(fā)射與接收的方式，有效保證了試驗(yàn)箱體內(nèi)部為恒溫邊界條件，減少外界環(huán)境溫度對(duì)測(cè)試結(jié)果的影響。

4 熱物性測(cè)試試驗(yàn)分析

4.1 糧食顆粒熱物性參數(shù)測(cè)試

選取淮南地區(qū)含水率為10.5%的皖稻121。首先利用孔徑為2 mm的篩網(wǎng)對(duì)稻谷進(jìn)行篩選，然后裝入直徑為400 mm，高為80 mm的短圓柱筒試樣盒中，填充密度580 kg/m3，環(huán)境溫度20 ℃，恒溫箱溫度保持90 ℃。在離試樣筒壁面和底面均為40 mm的位置處布置2個(gè)溫度測(cè)點(diǎn)，取平均值作為最終測(cè)試值，依次編號(hào)為測(cè)點(diǎn)1和測(cè)點(diǎn)2。某一時(shí)刻τ突然將圓柱體試樣盒放入恒溫箱內(nèi)，同時(shí)觀察銅板內(nèi)壁面溫度變化，當(dāng)銅板內(nèi)壁溫度穩(wěn)定且趨近于恒溫箱體內(nèi)部溫度時(shí)，開(kāi)始采集稻谷內(nèi)部測(cè)點(diǎn)溫度。為消除試驗(yàn)初期，測(cè)試溫度的不穩(wěn)定對(duì)測(cè)算結(jié)果的影響，選取100 s以后的采集數(shù)據(jù)為分析對(duì)象，試驗(yàn)共進(jìn)行3 600 s，溫度采樣間隔10 s，各測(cè)點(diǎn)溫度曲線見(jiàn)圖6。

熱擴(kuò)散率α初始猜測(cè)值為2.02×10-7m2/s，導(dǎo)熱系數(shù)λ和Cp初始猜測(cè)值分別為0.099 0 W/(m·℃)和823 J/(kg·℃)，取參數(shù)估計(jì)結(jié)果見(jiàn)表1。

圖6 測(cè)點(diǎn)1、測(cè)點(diǎn)2的溫升曲線

可以看出熱擴(kuò)散率和導(dǎo)熱系數(shù)估計(jì)結(jié)果與文獻(xiàn)[17]中的測(cè)算結(jié)果較為接近，但比熱容誤差較大(>10%)，而修正后的比熱容則較為吻合。為進(jìn)一步驗(yàn)證測(cè)算結(jié)果準(zhǔn)確性，對(duì)同一種試樣在相同條件下進(jìn)行多次試驗(yàn)，通過(guò)相對(duì)偏差△η衡量測(cè)算值彼此接近的程度。

(8)

式中：

由表2可以看出，導(dǎo)熱系數(shù)與熱擴(kuò)散率參數(shù)反演結(jié)果以及修正后的比熱容Cp*與對(duì)應(yīng)平均值的最大相對(duì)偏差均小于7%，即試驗(yàn)滿足可重復(fù)的要求。

表1 參數(shù)估計(jì)結(jié)果

表2 參數(shù)估計(jì)試驗(yàn)可重復(fù)精度分析

4.2 理論計(jì)算與數(shù)值仿真結(jié)果對(duì)比

將表2中獲取的熱物性參數(shù)的平均值代入式(4)中，計(jì)算τ在100～3 600 s時(shí)的溫度理論值，同時(shí)在Fluent中仿真模擬稻谷內(nèi)部的溫度場(chǎng)變化。數(shù)據(jù)曲線對(duì)比見(jiàn)圖7。

由圖7可以看出，稻谷的仿真曲線、理論計(jì)算曲線與實(shí)測(cè)溫升曲線三者較為吻合，進(jìn)一步驗(yàn)證了參數(shù)估計(jì)的準(zhǔn)確性。值得注意的是，試驗(yàn)至約2 600 s后，實(shí)測(cè)溫升與仿真和理論計(jì)算溫升開(kāi)始出現(xiàn)偏差，且隨時(shí)間的延長(zhǎng)，偏差具有變大趨勢(shì)。原因可能是熱量傳遞至試驗(yàn)箱上端絕熱邊界后因無(wú)法及時(shí)轉(zhuǎn)移，造成熱量積聚效應(yīng)，進(jìn)而對(duì)測(cè)點(diǎn)溫度產(chǎn)生影響，使實(shí)測(cè)溫升偏大，且隨著熱量積聚越來(lái)越多，影響也就越來(lái)越明顯。從而為盡量消除此因素的影響，選取100～2 500 s 時(shí)的溫升為最終有效測(cè)算數(shù)據(jù)。

4.3 3種糧食顆粒熱物性參數(shù)測(cè)試分析

為進(jìn)一步驗(yàn)證測(cè)試系統(tǒng)的穩(wěn)定性與準(zhǔn)確性，分別選取材料Y兩優(yōu)900水稻、聯(lián)創(chuàng)11號(hào)和源育15玉米作為試樣，含水率分別為13.1%，5.0%，12.5%，填充密度分別為600，720， 800 kg/m3。

圖7 實(shí)測(cè)、理論與數(shù)值仿真溫升曲線對(duì)比

由表3數(shù)據(jù)可以看出，測(cè)得的熱導(dǎo)率λ與文獻(xiàn)[16]和[17]數(shù)據(jù)較為接近，其中“源育15”號(hào)玉米的導(dǎo)熱系數(shù)與文獻(xiàn)數(shù)據(jù)相對(duì)誤差最大，為4.64%，分析原因可能是不同品種的玉米本身的熱物性參數(shù)存在差異，且玉米粒與稻谷相比顆粒直徑較大，存入試驗(yàn)箱時(shí)，玉米顆粒之間的孔隙較大，測(cè)溫過(guò)程溫度波動(dòng)較大，最終導(dǎo)致測(cè)算結(jié)果與參考值有偏差，但仍小于5%。符合測(cè)試精度要求。

表3 3種糧食試樣測(cè)試結(jié)果綜合分析

5 結(jié)論

(1) 采用恒溫箱結(jié)合黃銅短圓柱試樣盒構(gòu)建的恒溫邊界模型便于控制實(shí)現(xiàn)，且運(yùn)用非穩(wěn)態(tài)導(dǎo)熱乘積法結(jié)合參數(shù)估計(jì)法對(duì)糧食顆粒的熱物性參數(shù)進(jìn)行測(cè)試在原理與技術(shù)上均是可行的。

(2) 采用ZigBee無(wú)線數(shù)據(jù)采集方式，數(shù)據(jù)采集更方便，且不受試驗(yàn)箱體結(jié)構(gòu)的限制，降低了試驗(yàn)操作的難度。

(3) 糧粒的比熱容Cp的靈敏度系數(shù)較小，需先估計(jì)再修正，才能得到較為準(zhǔn)確的估計(jì)值。利用估算的熱物性參數(shù)值通過(guò)理論計(jì)算與數(shù)值仿真溫升變化曲線與實(shí)測(cè)溫度變化曲

線一致，進(jìn)一步說(shuō)明了參數(shù)估計(jì)的準(zhǔn)確性；對(duì)包括皖稻121在內(nèi)的4種糧食進(jìn)行熱物性測(cè)算，結(jié)果與相關(guān)文獻(xiàn)吻合，能夠滿足實(shí)際應(yīng)用要求。