賈瀟然，趙景堯，蔡英春(東北林業(yè)大學材料科學與工程學院，哈爾濱150040)

鋸材高頻-對流聯(lián)合加熱干燥傳熱傳質(zhì)數(shù)學模型

賈瀟然，趙景堯，蔡英春
(東北林業(yè)大學材料科學與工程學院，哈爾濱150040)

本文以多孔介質(zhì)傳熱傳質(zhì)理論為基礎(chǔ)，結(jié)合高頻-對流聯(lián)合加熱干燥過程中熱質(zhì)傳遞的機理及特性，通過質(zhì)量、動量、能量守恒方程等建立了高頻-對流聯(lián)合加熱干燥過程中鋸材熱質(zhì)傳遞數(shù)學模型并給出相應定解條件。模型中每個自變量及因變量都有獨立的控制方程，利用控制容積法建立相應控制方程的差分格式，將耦合偏微分方程按步驟解耦，以便于利用Matlab或FORTRAN等編程求解時，達到每個自變量均可以獨立求解的目的。

鋸材；高頻-對流干燥；數(shù)學模型；差分方程

1 前言

高頻-對流聯(lián)合加熱木材干燥，可實現(xiàn)大斷面鋸材高效加熱，易于控制其溫度梯度、周圍干燥介質(zhì)的相對濕度，因而是適用于大斷面木材快速、低成本、高品質(zhì)干燥的先進技術(shù)，具有廣闊應用前景。木材干燥過程是一個典型的非穩(wěn)態(tài)不可逆熱質(zhì)耦合過程，包括多相多組分的混合熱質(zhì)轉(zhuǎn)移效應和水分的相變過程，與木材材性、狀態(tài)、高頻加熱熱流密度、以及空氣流對材料表面?zhèn)鳠醾髻|(zhì)系數(shù)等密切相關(guān)，影響因素眾多、關(guān)系復雜。目前，木材對流加熱干燥方法在我國甚至全世界范圍內(nèi)仍占主導地位，高頻真空干燥技術(shù)應用也日趨廣泛，研究者們已經(jīng)對采用該兩種干燥技術(shù)干燥過程中木材熱質(zhì)遷移規(guī)律進行了大量研究，并能較好地描述干燥過程木材熱質(zhì)遷移規(guī)律[1～6]，但對高頻-對流聯(lián)合加熱干燥傳熱傳質(zhì)數(shù)學模型的建立及求解尚未見報道。本研究根據(jù)木材中水分和熱量的遷移機理，以非平衡熱力學理論和相平衡理論為基礎(chǔ)，根據(jù)質(zhì)量守恒定律、動量守恒定律及能量守恒定律確定描述熱質(zhì)遷移過程的控制方程，從而建立描述木材高頻-對流聯(lián)合加熱干燥過程中木材內(nèi)部熱質(zhì)遷移的數(shù)學模型，應用計算機數(shù)值模擬技術(shù)，可以精確有效地對其求解，以期能及時精確把握該雙熱源干燥過程中木材溫度及水分分布，為干燥過程的精準自動控制奠定基礎(chǔ)，為干燥工藝的改進及可靠實施提供必要的信息。

2 模型建立的基礎(chǔ)

根據(jù)高頻-對流聯(lián)合加熱干燥過程中木材內(nèi)部水分在干燥不同階段呈現(xiàn)不同存在形式，將整個干燥過程分為纖維飽和點以上和纖維飽和點以下兩個階段。根據(jù)木材高頻以及對流加熱干燥熱質(zhì)傳遞機理可知，在纖維飽和點以上，木材內(nèi)部自由水、空氣和水蒸氣在總壓力（毛細管張力、空氣壓力和水蒸氣壓力等）梯度作用下通過由連續(xù)的細胞腔-細胞壁上紋孔所構(gòu)成的大毛細管路徑向外滲流。盡管沿纖維方向滲流阻力小（滲透率大），但由于鋸材長度遠大于厚度和寬度，所以該方向滲流占比小，甚至可以忽略。干燥室內(nèi)由木隔條隔開的鋸材上下表面水分向干燥介質(zhì)中蒸發(fā)，由于鋸材側(cè)面不留間隙，所以無水分蒸發(fā)，常用鋸材可忽略板寬方向水分遷移，但大斷面鋸材（如柱材）則應給予重視。在纖維飽和點以下，吸著水擴散或蒸發(fā)到細胞腔后在混合壓力梯度作用下向外遷移，沿各方向的擴散占比與自由水遷移相近，液態(tài)吸著水在細胞壁內(nèi)的遷移占比很小[7]。物理模型通常包括過程定性分析及若干簡化假定兩部分。對實際過程及材料結(jié)構(gòu)進行簡化假定時，既要注意假定盡可能接近實際，又要便于模型建立，因此引入[6，8，9]并分析建立假定如下。

1）將實際復雜的木材結(jié)構(gòu)視為由諸多不同直徑的毛細管組成的一個網(wǎng)絡(luò)系統(tǒng)，存在于大毛細管內(nèi)的液相水為自由水，其中液相水與氣相水蒸氣始終保持著相平衡，大毛細管內(nèi)水蒸氣處于飽和狀態(tài)。存在于微毛細管內(nèi)的液相水為吸著水，其中液相水與氣相水蒸氣始終保持相平衡。

2）自由水遷移的驅(qū)動力是內(nèi)外混合氣體壓力差以及毛細管力的合力，遷移形式為滲流。吸著水既有滲流，同時也存在擴散，但由于滲流量較小，可忽略不計。

3）細胞溫度對應下的飽和水蒸氣壓力與空氣壓力之和等于細胞腔內(nèi)的總壓力，空氣被視為理想氣體?？諝馀c水蒸氣滲流速度相同，即為混合氣體滲流速度。

4）在流體和木材實質(zhì)之間的熱轉(zhuǎn)移相當快，認為木材中所有相態(tài)處于局部熱力學平衡狀態(tài)，即各相的溫度相同。

5）盡管實際木材干燥過程中吸著水的遷出使木材產(chǎn)生干縮，但便于分析，假定固相骨架是剛性的，即木材的尺寸和內(nèi)部結(jié)構(gòu)不發(fā)生變化。

6）建立和求解模型時，認為熱量和水分僅沿一維方向(材厚方向)進行遷移。

7）能量的轉(zhuǎn)移發(fā)生在各相之間的熱傳導、對流換熱和相變傳熱當中，忽略各種能量損失。

8）木材內(nèi)部的初含水率和初始溫度分布均勻。

3 熱質(zhì)遷移控制方程

3.1 質(zhì)量守恒方程

纖維飽和點以上：木材可以認為是由木材實質(zhì)、自由水、吸著水、水蒸氣和空氣構(gòu)成的多相混合體系，而在干燥過程中參與流動的液相是自由水，氣相是水蒸氣、空氣。

自由水質(zhì)量守恒：

式（1）為含水率M表征的自由水控制方程，含義是木材內(nèi)任意一處單位時間、單位體積內(nèi)自由水的質(zhì)量變化量等于從該處鄰近以滲流形式在單位時間內(nèi)自由水通過該處邊界滲入的質(zhì)量減去其內(nèi)部自由水因邊界內(nèi)蒸發(fā)而減少的質(zhì)量。式（1）中，ρd為木材絕干密度，kg/m3；vf為自由水滲流速度，m/s；mv為體積蒸發(fā)率（相變率），kg/m3?s；M為含水率；τ為時間，s；z為距離，m。

水蒸氣質(zhì)量守恒：

式（2）為水蒸氣密度的控制方程，由于飽和水蒸氣密度由飽和溫度決定，因此，將式（2）變形為式（3）的體積蒸發(fā)率mv的控制方程，含義是木材內(nèi)任意一處單位時間單位體積內(nèi)自由水蒸發(fā)為水蒸汽的質(zhì)量mv等于單位時間單位體積內(nèi)水蒸汽質(zhì)量的增加量加上單位時間單位體積內(nèi)水蒸汽通過該體積邊界滲入的質(zhì)量。式（2）、（3）中，ρv為水蒸氣密度，kg/m3；vg為混合氣體（水蒸氣和空氣）滲流速度，m/s；?為孔隙率。

式（4）為空氣密度控制方程，含義是木材內(nèi)任意一處單位時間、單位體積內(nèi)空氣的質(zhì)量變化量等于從該處鄰近以滲流形式在單位時間內(nèi)空氣通過該處邊界滲入的質(zhì)量。式（4）中，ρa為空氣密度，kg/m3。

纖維飽和點以下：木材內(nèi)部的自由水全部蒸發(fā)完畢，木材中的吸著水向臨近的細胞腔蒸發(fā)，再以水蒸氣的形式向木材表面遷移。

式（5）為含水率表征的吸著水控制方程，含義是木材內(nèi)任意一處單位體積、單位時間吸著水質(zhì)量減少的量等于單位體積單位時間內(nèi)其內(nèi)部吸著水蒸發(fā)為水蒸氣的質(zhì)量；纖維飽和點以下體積蒸發(fā)率控制方程同式（2），含義是木材內(nèi)任意一處單位時間單位體積內(nèi)自由水蒸發(fā)為水蒸氣的質(zhì)量mv等于單位時間、單位體積內(nèi)水蒸氣質(zhì)量的增加量加上單位時間、單位體積內(nèi)水蒸氣通過該體積邊界滲入的質(zhì)量?？諝饷芏瓤刂品匠掏剑?）。

3.2 動量守恒方程

纖維飽和點以上：根據(jù)達西定律[10]，自由水、水蒸氣和空氣混合氣體在一維方向上滲流速度控制方程如下：

纖維飽和點以下：混合氣體滲流速度同式（7），式（6）、式（7）中，μ為動力粘度，Pa?s；Kf為自由水有效滲透率（滲透性），m2/Pa?s；Kg為混合氣體有效滲透率，m2/Pa?s；Pf為自由水壓力；Pg為水蒸氣和空氣混合氣體總壓力。

3.3 熱力學關(guān)系式

纖維飽和點以上：液相水為自由水，與它處于相平衡的水蒸氣為飽和水蒸氣，由溫度T即可得相應的飽和水蒸氣密度ρv與飽和水蒸汽壓力Psv（即水蒸氣壓力Pv），具體運算時通過“飽和水與飽和水蒸氣熱力性質(zhì)表”查得，木材內(nèi)部空氣可視為理想氣體，空氣壓力可通過理想氣體狀態(tài)方程求得，空氣壓力Pa、水蒸氣壓力Pv、水蒸氣密度ρv控制方程如下[11]：

式（8）中，R為通用氣體常數(shù)；Ma為空氣摩爾質(zhì)量，kg/mol。

纖維飽和點以下：液相水為吸著水，與它處于相平衡的水蒸氣為未飽和水蒸氣（過熱水蒸氣）。未飽和水蒸氣的壓力對應著當?shù)禺敃r濕空氣相對濕度φ[12]，水蒸氣壓力Pv可通過當?shù)禺敃r相對濕度及溫度對應的飽和水蒸氣壓力Psv得到，在已知Pv條件下，水蒸氣密度由理想氣體狀態(tài)方程求得，因此控制方程如下：

式（12）中，Mv為水蒸氣摩爾質(zhì)量，kg/mol。

3.4 自由水、混合氣體壓力

細胞腔內(nèi)總壓力見式（14），加熱初期自由水未汽化，自由水的驅(qū)動力等于木材毛細張力Pc[13]，見式（15）；當自由水溫度達到汽化溫度時，毛細張力影響將大大降低，自由水驅(qū)動力為細胞腔內(nèi)總壓力，根據(jù)假設(shè)中假定纖維飽和點以上水蒸氣處于飽和狀態(tài)，因此自由水壓力如式（16）所示：

式（15）中，A、n為常數(shù)；Sf為液相水飽和度。

3.5 能量守恒方程

纖維飽和點以上：

式（17）為纖維飽和點以上溫度T的控制方程，含義是木材單位體積內(nèi)自由水的能量在單位時間內(nèi)的增加量等于以熱傳導形式在單位時間內(nèi)從單位體積邊界傳入的能量加上自由水、水蒸氣以滲流的方式在單位時間內(nèi)從單位體積邊界傳入的能量減去該單位體積內(nèi)液相水單位時間蒸發(fā)而減少的能量加上單位體積吸收的高頻功。式（17）中，ρ為木材實際密度，kg/m3；cp為木材比熱，J/kg?K；λeff為有效導熱系數(shù)，W/m?K；γ為水汽化潛熱，J/kg；q為高頻電功率密度，W?m3；cf為自由水比熱，J/kg?K；cv為水蒸氣比熱，J/kg?K。

纖維飽和點以下：

式（18）為纖維飽和點以下溫度T的控制方程，含義是木材單位體積內(nèi)吸著水的能量在單位時間內(nèi)的增加量等于以熱傳導形式在單位時間內(nèi)從單位體積邊界傳入的能量加上水蒸氣以滲流的方式在單位時間內(nèi)從單位體積邊界傳入的能量減去該單位體積內(nèi)液相水單位時間蒸發(fā)而減少的能量加上單位體積吸收的高頻功。

其他參數(shù)如下：

λeff為有效導熱系數(shù)[10]：

當M≥0.4時：

當M＜0.4時：

cp為木材比熱[14]：

ρ為木材在某一含水率時刻的實際密度[15]：

γ為水汽化潛熱[16]：

纖維飽和點以上：

纖維飽和點以下：由于部分吸著水被吸附在木材細胞壁上的初級吸著點和次級吸著點上，部分凝結(jié)于微毛細管內(nèi)，使其蒸發(fā)，必須克服微毛細管凝結(jié)力、破壞水分子與各級吸著點之間的氫鍵結(jié)合，即與自由水相比，需要額外增加一個潤濕熱F[17]。

式(25)中，q為高頻電功率密度[7]；f為電場頻率，Hz；E為電場強度，V/m；ε為介電常數(shù)；tanδ為介電損耗因子。

3.6 表面蒸發(fā)率及界面蒸發(fā)率

對流干燥過程中，材料表面的液相水將蒸發(fā)，因此表面存在表面蒸發(fā)率ms，kg/m2?s。

式（26）為表面蒸發(fā)率的控制方程，該式應在式（1）的邊界條件中出現(xiàn)，含義是單位時間從材料內(nèi)部以滲流形式到達單位表面積水的質(zhì)量等于單位時間單位表面積上由自由水蒸發(fā)變成水蒸氣的質(zhì)量。當表面處含水率在纖維飽和點以上，自由水浸潤材料整個表面，在表面就地蒸發(fā)的水蒸氣被環(huán)境空氣流以對流形式帶走，則ms的表達式為式（27）。式中hm為對流傳質(zhì)系數(shù)，m/s；ρve為環(huán)境水蒸氣密度，kg/ m3。

當材料表面含水率降至纖維飽和點以下時，該處水是吸著水，其水蒸氣密度不僅與溫度有關(guān)，更有賴于當?shù)氐暮?。因此假定材料表面對環(huán)境的傳質(zhì)驅(qū)動力是表面處含水率與對應于環(huán)境條件下的木材平衡含水率兩者之間的差值[6]。因此ms*的表達式為式（28）。式中hm*為對流傳質(zhì)系數(shù)，m/s；Me為環(huán)境平衡含水率。

由于高頻加熱導致木材內(nèi)部液相水處于汽化狀態(tài)，并且混合氣體流速度較常規(guī)干燥中混合氣體擴散速度要大得多，因此不存在明顯的蒸發(fā)界面向木材內(nèi)部移動的情況，故不考慮移動蒸發(fā)界面（界面蒸發(fā)率）。

4 定解條件

4.1 初始條件

圖1、圖2分別表述了纖維飽和點以上、纖維飽和點以下干燥模型因變量間耦合關(guān)系。加熱過程開始之前（含開始瞬時），根據(jù)圖1應給出含水率M、溫度T、木材內(nèi)氣體總壓力Pg（Pg＝Pe），根據(jù)溫度可算出水蒸氣壓力Pv、水蒸氣密度ρv、空氣壓力Pa、空氣密度ρa、自由水壓力Pf等參數(shù)的初始值，自由水滲流速度vf、混合氣體滲流速度vg、體積蒸發(fā)率mv、高頻功率密度q為0。

圖1 纖維飽和點以上因變量耦合關(guān)系Fig.1 The coupling relationship between the dependent variables above FSP

4.2 邊界條件

木材厚為H，木材中心為原點即z=0，木材表面為外邊界即z=H/2。

z=0處：對于含水率M、溫度T的邊界條件如下：

圖2 纖維飽和點以下因變量耦合關(guān)系Fig.2 The coupling relationship between the dependent variables under FSP

z=H/2處：含水率M的邊界條件為式（26）、式（27）、式（28）。

含水率T的邊界條件為式（31），式中h為換熱系數(shù)，W/m2?K；Te為環(huán)境溫度，℃。

體積蒸發(fā)率mv，其本身無需邊界條件，但在尋求邊界處mv時，需要水蒸氣密度ρv的邊界條件，為此水蒸氣密度ρv的邊界條件為式（32）。

5 差分方程

熱質(zhì)遷移控方程組是非線性的，方程系數(shù)都是熱質(zhì)變量（T，M）的函數(shù)，且溫度場、含水率場、壓力場互相耦合，形成了一組復雜的非線性偏微分方程組，因此若要求出該偏微分方程組的分析解極為困難。然而，借助于計算機的數(shù)值解法，問題可以得到解決。其基本思想是把原來在空間與時間坐標中連續(xù)的物理量的場（如速度場、溫度場、濃度場），用一系列有限個離散點（節(jié)點）上的值的集合來代替，通過一定的原則建立起這些離散點上變量值之間關(guān)系的代數(shù)方程(離散方程)，求解所建立起來的代數(shù)方程以獲得所求解變量的近似值[8，18，19]。本文采用有限差分法及控制容積法對熱質(zhì)遷移控制方程和相應的邊界條件在計算區(qū)域內(nèi)進行離散化。將木材沿厚度方向從中心到木材表面按間距（空間步長）Δx分成JZ層，則節(jié)點j＝1，2，3，…，JZ；干燥總時間=時間步長Δτ*n，n=1，2，3，…。由于方程中存在速度v，它與壓力P的一階導數(shù)成正比，當一階導數(shù)寫成差分格式后速度應在兩個壓力值之間，因此本文使用的是交錯網(wǎng)格，即將壓力等對應的單元所建立的網(wǎng)格系統(tǒng)成為主網(wǎng)格，速度所對應的網(wǎng)格與主網(wǎng)格相差半個空間步長，引入交錯網(wǎng)格時保證計算木材內(nèi)部流體速度的重要條件[6]，如圖3所示。速度被定義在單元的邊界，節(jié)點用小方塊表示，壓力、密度被定義在單元的中心，節(jié)點用圓圈表示。由于邊界條件的存在會使在不同單元差分方程的形式不同，因此選取存在具有代表性的變量差分方程如下。

圖3 交叉網(wǎng)格Fig.3 Cross grid

5.1 含水率

式（1）差分方程

1＜j＜JZ：內(nèi)部單元

j=1：左邊界

j=JZ：右邊界

5.2 體積蒸發(fā)率

式（3）差分方程

1＜j＜JZ：內(nèi)部單元

j=1：左邊界

j=JZ：右邊界

5.3 滲流速度

式（6）、（7）差分方程

j=1：左邊界

1＜j＜JZ+1：內(nèi)部單元

j=JZ+1：右邊界

5.4 空氣密度

式（4）差分方程

j=1：左邊界

1＜j≤JZ：內(nèi)部單元、右邊界

5.5 溫度

式（17）差分方程

1＜j＜JZ：內(nèi)部單元

j=1：左邊界

j=JZ：右邊界

5.6 其他參數(shù)

除上述參數(shù)控制方程外，其他參數(shù)控制方程的差分格式同理可推。

6 差分方程求解

差分方程分顯式與隱式兩種，為了便于求解，上文中采用的是顯式差分格式，每個變量都有獨立的控制方程，這樣將原本的耦合方程解耦，變成單獨的線性方程，每個因變量都可以單獨求解。利用Matlab或者FORTRAN軟件編程，通過n時層的已知值求解n+1時層值，如此反復，直到干燥結(jié)束。

7 結(jié)語

本文結(jié)合高頻-對流聯(lián)合加熱干燥過程中熱質(zhì)傳遞的機理及特性，通過質(zhì)量、動量、能量守恒方程等建立了高頻-對流聯(lián)合加熱干燥過程中鋸材熱質(zhì)傳遞數(shù)學模型并給出相應定解條件，利用控制容積法建立了相應控制方程的差分格式，將耦合偏微分方程按步驟解耦。

作為本研究的繼續(xù)，筆者將通過反復計算模擬，逐步完善控制方程以及變量參數(shù)，在此基礎(chǔ)上進行二維或者三維的數(shù)學模型建立，以期能到達模擬真實干燥過程。最后通過模擬計算出的數(shù)據(jù)與實驗數(shù)據(jù)進行對比驗證，進一步修正模型各項參數(shù)，給出精確的模型及解法。

[1]AnastasiosKounmoutsakos，StavrosAvramidis，SavvasGHatzikiriakos.Radio frequency vacuum drying of wood：I.Mathematical model[J].Drying Technology，2001，19(1)：65-84.

[2]AnastasiosKounmoutsakos，StavrosAvramidis，SavvasGHatzikiriakos.Radio frequency vacuum drying of wood：II.Experimental model evaluation[J].Drying Technology，2001，19(1)：85-98.

[3]盧 濤.毛細多孔介質(zhì)干燥過程中傳熱傳質(zhì)模型研究及應用[D].大連：大連理工大學，2003.

[4]劉 偉.多孔介質(zhì)傳熱傳質(zhì)理論與應用[M].北京：科學出版社，2006.

[5]Jia Donghua，Muhammad T Afzal.Modeling of moisture diffusion in microwave drying of hardwood[J].Drying Technology，2007，25(3)：449-454.

[6]俞昌銘.多孔材料傳熱傳質(zhì)數(shù)值分析[M].北京：高等教育出版社，2011.

[7]蔡英春.木材高頻真空干燥機理[M].哈爾濱：東北林業(yè)大學出版社，2006.

[8]陶文銓.數(shù)值傳熱學[M].二版.西安：西安交通大學出版社，2001.

[9]于建芳.木材微波干燥熱質(zhì)轉(zhuǎn)移及其數(shù)值模擬[D].呼和浩特：內(nèi)蒙古農(nóng)業(yè)大學，2010.

[10]Siau J F.Transport Processes in wood[M].New York：Springer-Verlga，1984.

[11]嚴家綠，王永青.工程熱力學[M].北京：中國電力出版社，2007.

[12]Avramidis S.Evaluation of“three-Variable”models for the pre-diction of equilibrium moisture content in wood[J].Wood Science and Technology，1989，23(3)：251-257.

[13]成俊卿.木材學[M].北京：中國林業(yè)出版社，1985.

[14]Zhangjing C.Primary drying forces in wood vacuum drying[D]. Virginia：Virginia Tech，1997.

[15]楊世銘，陶文銓.傳熱學[M].北京：高等教育出版社，1998.

[16]伊松林，張璧光，常建民.木材浮壓干燥過程的傳熱傳質(zhì)[D].北京：北京林業(yè)大學，2002.

[17]陸金甫，關(guān) 治.偏微分方程數(shù)值解法[M].北京：清華大學出版社，2004.

[18] SPOLEK GA Plumb OA.Capillary pressure in softwoods[J]. Wood Science and Technology，1981，15：189-199.

[19]李賢軍.木材微波-真空干燥特性的研究[D].北京：北京林業(yè)大學，2005.

Modeling of heat and mass transfer for timber during radio-frequency and convective drying

Jia Xiaoran，Zhao Jingyao，Cai Yingchun
(College of Material Science and Engineering，Northeast Forestry University，Harbin 150040，China)

The research based on the theory of porous medium heat and mass transfer，combined with the mechanism and characteristics of heat mass transfer in radio-frequency and convective drying，establish the mathematical model and give the definite condition through mass momentum and energy conservation equations.Each independent and dependent variable has independent control equation.Use the finite volume method to establish difference equation for to decouple the coupled partial differential equations，so that each independent variable can be solved independently when programming in Matlab or FORTRAN.

timber；radio-frequency-convective drying；mathematical model；difference equation

TS652

A

1009-1742-（2014）04-0106-07

2014-01-27

林業(yè)公益性行業(yè)科研專項經(jīng)費項目（201304502）

賈瀟然，1984年出生，男，黑龍江哈爾濱市人，博士研究生，主要研究方向為木材高頻真空干燥；E-mail：80449337@qq.com