周 正 孫麗萍 姜 濱

(東北林業(yè)大學(xué)，哈爾濱，150040)

基于控制體積有限元方法的木材干燥過程含水率分布模型1)

周 正 孫麗萍 姜 濱

(東北林業(yè)大學(xué)，哈爾濱，150040)

對木材干燥過程含水率分布進(jìn)行了數(shù)學(xué)建模。首先建立了液體守恒方程和能量守恒方程，然后利用控制體積有限元方法對守恒方程進(jìn)行離散，通過對控制體積層面的通量計(jì)算和離散方程的求解，建立了木材干燥過程含水率分布的仿真圖，驗(yàn)證了數(shù)學(xué)模型的準(zhǔn)確性。

木材干燥；木材含水率模型；控制體積有限元方法

Wood drying; Wood moisture content model; The control volume finite element method (CVFE)

木材含水率是衡量木材干燥質(zhì)量最重要的指標(biāo)之一[1]，如果能準(zhǔn)確地掌握木材干燥過程中含水率的變化情況，及時(shí)調(diào)整干燥窯的性能指標(biāo)，對于木材干燥的效果和減少原材料的損耗具有重大的意義[2]。目前研究木材含水率及其分布規(guī)律普遍利用二維模型[3]，李賢軍[4]等人的研究結(jié)果表明，在常規(guī)干燥和微波干燥時(shí)，木材內(nèi)部存在著整體性的內(nèi)高外低的含水率梯度場；徐兆軍[5]等人利用斷層掃描技術(shù)，通過數(shù)學(xué)模型描述了木材含水率的分布特征；還有一些研究，模擬了木材的木質(zhì)結(jié)構(gòu)和扭曲特性[6]；這些方法和數(shù)學(xué)模型，為研究木材干燥學(xué)，特別是木材含水率提供了基礎(chǔ)，同時(shí)不斷涌現(xiàn)的新方法也增強(qiáng)了該領(lǐng)域的創(chuàng)新性。目前，國內(nèi)外對木材干燥過程建立完整的耦合非線性數(shù)學(xué)模型的研究很少，本文提出的能夠較為精確描述木材干燥過程的三維模型，是在二維模型的基礎(chǔ)上建立的。本文利用控制體積有限元(CVFE)方法，建立了木材干燥過程中的數(shù)學(xué)模型；該三維模型從空間分布的角度，全面地描述了木材干燥過程的水分分布和遷移特性，為控制干燥窯的環(huán)境參數(shù)提供了參考依據(jù)。

1 數(shù)學(xué)模型的建立

1.1 守恒方程

這里建立的數(shù)學(xué)模型由2個(gè)非線性偏微分方程構(gòu)成，即：液體守恒方程、能量守恒方程。在試驗(yàn)中，為了減少計(jì)算量，只研究板材的四分之一，板材內(nèi)部所產(chǎn)生的氣壓忽略不計(jì)。

液體守恒方程：

(1)

能量守恒方程：

(2)

式中：X為木材含水率，其值取決于由自由水含水率、結(jié)合水含水率Xb和溫度T；ρ為密度；ε為體積系數(shù)；v為相速；ω為質(zhì)量分?jǐn)?shù)；h為焓；下標(biāo)a、b、g、s、v、w，分別表示空氣、結(jié)合水、氣體、細(xì)胞壁、水蒸氣、自由水；ρ0為木材密度；Db為自由水?dāng)U散系數(shù)；Dv為有效水蒸氣擴(kuò)散系數(shù)，Keff為有效導(dǎo)熱系數(shù)。

1.2 邊界條件

板材干燥平面的邊界條件為：

(3)

(4)

1.3 方程的離散

利用控制體積有限元方法離散化的數(shù)學(xué)模型由守恒方程式(1)和式(2)推導(dǎo)出來，改寫之后的守恒方程為：

(5)

(6)

式中：Ψw為液體的守恒量；Ψe為能量的守恒量；Jw為液體的通量向量；Je為能量的通量向量。

Ψw和Ψe可表示為：

Ψw=ρ0X+εgρv；

(7)

(8)

Jw和Je可表示為：

Jw=ρwvw-ρ0DbXb-ρgDvωv；

(9)

Je=ρwhwvw-hbρ0DbXb-ρg(hv-ha)Dvωv- Kefft。

(10)

偏微分方程，式(1)和式(2)的離散方程式(5)和式(6)，是通過對控制體積積分和應(yīng)用高斯散度定理求得的，求得的表面上的積分相當(dāng)于表面上的離散和，又有：

(11)

(12)

1.4 離散方程的解法

解離散方程的方法采用牛頓迭代算法。應(yīng)用控制體積有限元方法解在網(wǎng)格每個(gè)節(jié)點(diǎn)上的守恒方程，能夠得到一個(gè)2×N的非線性離散方程組：

F(u)=[Fw1(u)，F(xiàn)e1(u)，F(xiàn)w2(u)，F(xiàn)e2(u)，…，F(xiàn)wN(u)，F(xiàn)eN(u)]T=0。

(13)

式中：Fwi(u)和Fei(u)分別為網(wǎng)格中第i個(gè)控制體積的液體離散守恒方程組和能量離散守恒方程組；解向量u=(X1,T1，X2，T2，…，XN,YN)，u包括網(wǎng)格中每個(gè)節(jié)點(diǎn)上成對的主要變量。

解線性方程組采用穩(wěn)定雙共軛梯度算法，該算法常用于解決大規(guī)模矩陣方程組[7]；這種算法經(jīng)常與預(yù)處理算法ILU(0)一起使用[8]，這樣能大大提高系統(tǒng)的運(yùn)算速度[9]。解線性方程組時(shí)，首先要求出雅各比矩陣，雅各比矩陣是對積分項(xiàng)和通量項(xiàng)進(jìn)行求導(dǎo)得出的，計(jì)算導(dǎo)數(shù)的方法為一階有限差分逼近，可用公式(14)表示。

(14)

式中：ewj為與主要變量X有關(guān)的第j個(gè)單元的向量；τ為位移值，τ必須足夠大以避免產(chǎn)生舍入誤差，并足夠小以避免產(chǎn)生導(dǎo)數(shù)的錯(cuò)誤估計(jì)值，式(14)中τ=10-6。

為了盡可能的減小牛頓算法中的浮點(diǎn)誤差，尤其是構(gòu)造雅各比矩陣時(shí)產(chǎn)生的誤差，非線性方程組可改寫為：

(15)

(16)

式中：D1和D2為對角矩陣。

1.5 控制體積表面的通量計(jì)算

1.5.1 擴(kuò)散項(xiàng)的計(jì)算

本文中對附屬控制體積表面的擴(kuò)散張量和二次變量進(jìn)行計(jì)算的方法為平均法，即取2個(gè)節(jié)點(diǎn)上賦值的平均值。這2個(gè)節(jié)點(diǎn)構(gòu)成的線段，是截?cái)喔綄倏刂企w積的網(wǎng)格邊界(見圖1)，圖1中所標(biāo)注的點(diǎn)ab和cd就是上述2個(gè)節(jié)點(diǎn)。

圖1 擴(kuò)散通量的計(jì)算方法

1.5.2 平流項(xiàng)和對流項(xiàng)的計(jì)算

計(jì)算平流項(xiàng)和對流項(xiàng)的方法為極限通量算法，這種算法能夠有效避免自激振蕩引起的高階空間離散，從而消除了可能在解域中出現(xiàn)的基波或震蕩。極限通量算法計(jì)算附屬控制體積表面液體遷移率的公式為：

λf=λu+σ(r)(λn-λu)/2。

(17)

式中：σ為限制器函數(shù)，0≤σ≤2；當(dāng)σ=0時(shí)，為迎風(fēng)狀態(tài)；當(dāng)σ=1時(shí)，為平衡狀態(tài)；當(dāng)σ=2時(shí)，為順風(fēng)狀態(tài)。r為流向指示器的比率。在P. Arminjon[10]的研究中，σ(r)表示為：

(18)

流向指示器的比率為：

r=I2/I。

(19)

式中：I為處于迎風(fēng)節(jié)點(diǎn)和順風(fēng)節(jié)點(diǎn)中間的控制體積表面；I2為第一個(gè)迎風(fēng)節(jié)點(diǎn)和第二個(gè)迎風(fēng)節(jié)點(diǎn)中間的控制體積表面。

2 仿真結(jié)果與分析

試驗(yàn)中，把板材的立體模擬圖劃分為若干個(gè)三角棱柱體，每個(gè)子單元都能在三維坐標(biāo)系中表達(dá)出來。圖2是四分之一板材的網(wǎng)格模型，橫斷面上共有1022個(gè)子單元，徑向的幾何比率為1.05，整體總共有12270個(gè)網(wǎng)格。

圖2 三維空間的板材有限元網(wǎng)格模型

圖3、圖4分別為板材(四分之一)干燥1 h 30 min時(shí)的仿真圖和3 h 10 min時(shí)的仿真圖。圖3、圖4中的木材含水率范圍不同，是由于隨著木材干燥時(shí)間的增加，板材中的水分不斷地遷移并蒸發(fā)。圖3、圖4充分體現(xiàn)了具有非均勻結(jié)構(gòu)的木材，在干燥時(shí)的特點(diǎn)。從圖3、圖4可見：木材中的水分，優(yōu)先選擇縱向流動(dòng)，這是因?yàn)槟静脑诳v向上的滲透率是最高的。

3 結(jié)束語

本文建立了木材干燥過程含水率分布的數(shù)學(xué)模型。首先，給出了木材干燥過程含水率空間分布的液體守恒方程和能量守恒方程；然后，利用控制體積有限元原理將方程離散化，解離散方程采用了牛頓算法，其中，容限采用10-8；解線性方程時(shí)，采用了穩(wěn)定雙共軛梯度算法和預(yù)處理算法ILU(0)。計(jì)算通量時(shí)，擴(kuò)散張量和二次變量的計(jì)算應(yīng)用了平均值法，平流項(xiàng)和對流項(xiàng)的計(jì)算應(yīng)用了極限通量算法。最后，根據(jù)該數(shù)學(xué)模型，建立了木材干燥過程含水率分布的仿真圖，驗(yàn)證了利用控制體積有限元方法分析木材干燥含水率空間分布規(guī)律的可行性和準(zhǔn)確性。

圖3 1 h 30 min時(shí)板材含水率分布情況

圖4 3 h 10 min時(shí)板材含水率分布情況

[1] 朱政賢.干燥鋸材最終含水率與木材平衡含水率的初步探索[J].東北林業(yè)大學(xué)學(xué)報(bào)，1986，14(2)：1-10.

[2] 孫麗萍.木材含水率在線檢測融合體系及仿真技術(shù)研究[D].哈爾濱：東北林業(yè)大學(xué)，2008.

[3] Dedic A D， Mujumdar A S， Voronject D K. A three dimensional model for heat and mass transfer in convective wood drying[J]. Drying Technology，2003，21(1)：1-15.

[4] 李賢軍，喬建政，蔡智勇，等.微波干燥與常規(guī)干燥中木材內(nèi)含水率動(dòng)態(tài)分布[J].中南林業(yè)科技大學(xué)學(xué)報(bào)，2009，29(6)：98-113.

[5] 徐兆軍，丁建文，丁濤，等.基于斷層掃描圖像技術(shù)的木材纖維飽和點(diǎn)以上水分分布和遷移特性研究[J].木材加工機(jī)械，2010(1)：24-25，10.

[6] Plumb O A， Gong L Gong. Modelling the effect of heterogeneity on wood drying[M]//Turner I W， Mujumdar A S. Mathematical Modeling and Numerical Techniques in Drying Technology. New York： Chandan Kumar Ray，1996：221-258.

[7] 張恩澤，彭樹生，何小祥，等.超松弛迭代-雙共軛梯度在三維電磁問題有限元分析中的應(yīng)用[J].淮陰師范學(xué)院學(xué)報(bào)：自然科學(xué)版，2005，4(4)：292-295.

[8] 李熙銘，歐陽丹彤，白洪濤.基于GPU的混合精度平方根共軛梯度算法[J].儀器儀表學(xué)報(bào)，2012，33(1)：97-104.

[9] 金巍巍，陶文銓，何雅玲.代數(shù)方程求解方法收斂速度比較及對算法健壯性的影響[J].西安交通大學(xué)學(xué)報(bào)，2005，39(9)：966-970.

[10] Arminjon P， Dervieux A. Construction of TVD-like artificial viscosities on two-dimensional arbitrary FEM grids[J].Journal of Computational Physics，1993，106(1)：176-198.

周正，女，1989年6月生，東北林業(yè)大學(xué)機(jī)電工程學(xué)院，碩士研究生。E-mail：1060061063@qq.com。

孫麗萍，東北林業(yè)大學(xué)機(jī)電工程學(xué)院，教授。E-mail：zdhslp@163.com。

2013年10月16日。

S781.3

A Mathematical Model for Moisture Content of Board during Wood Drying Based on CVFE Method/Zhou Zheng, Sun Liping, Jiang Bin(Northeast Forestry University, Harbin 150040, P. R. China)//Journal of Northeast Forestry University.-2014,42(4).-124～126

1) 國家林業(yè)公益性行業(yè)科研專項(xiàng)(201304502)。

責(zé)任編輯：張 玉。

The complicated mathematical model during wood drying was built to describe moisture content, including the conservation equations, the discretisation process, flux approximations, and the solution of a nonlinear system. The simulation results were presented to investigate the performance of the mathematical model.