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

周 正 孙丽萍 姜 滨

(东北林业大学，哈尔滨，150040)

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

周 正 孙丽萍 姜 滨

(东北林业大学，哈尔滨，150040)

对木材干燥过程含水率分布进行了数学建模。首先建立了液体守恒方程和能量守恒方程，然后利用控制体积有限元方法对守恒方程进行离散，通过对控制体积层面的通量计算和离散方程的求解，建立了木材干燥过程含水率分布的仿真图，验证了数学模型的准确性。

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

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

木材含水率是衡量木材干燥质量最重要的指标之一[1]，如果能准确地掌握木材干燥过程中含水率的变化情况，及时调整干燥窑的性能指标，对于木材干燥的效果和减少原材料的损耗具有重大的意义[2]。目前研究木材含水率及其分布规律普遍利用二维模型[3]，李贤军[4]等人的研究结果表明，在常规干燥和微波干燥时，木材内部存在着整体性的内高外低的含水率梯度场；徐兆军[5]等人利用断层扫描技术，通过数学模型描述了木材含水率的分布特征；还有一些研究，模拟了木材的木质结构和扭曲特性[6]；这些方法和数学模型，为研究木材干燥学，特别是木材含水率提供了基础，同时不断涌现的新方法也增强了该领域的创新性。目前，国内外对木材干燥过程建立完整的耦合非线性数学模型的研究很少，本文提出的能够较为精确描述木材干燥过程的三维模型，是在二维模型的基础上建立的。本文利用控制体积有限元(CVFE)方法，建立了木材干燥过程中的数学模型；该三维模型从空间分布的角度，全面地描述了木材干燥过程的水分分布和迁移特性，为控制干燥窑的环境参数提供了参考依据。

1 数学模型的建立

1.1 守恒方程

这里建立的数学模型由2个非线性偏微分方程构成，即：液体守恒方程、能量守恒方程。在试验中，为了减少计算量，只研究板材的四分之一，板材内部所产生的气压忽略不计。

液体守恒方程：

(1)

能量守恒方程：

(2)

式中：X为木材含水率，其值取决于由自由水含水率、结合水含水率Xb和温度T；ρ为密度；ε为体积系数；v为相速；ω为质量分数；h为焓；下标a、b、g、s、v、w，分别表示空气、结合水、气体、细胞壁、水蒸气、自由水；ρ0为木材密度；Db为自由水扩散系数；Dv为有效水蒸气扩散系数，Keff为有效导热系数。

1.2 边界条件

板材干燥平面的边界条件为：

(3)

(4)

1.3 方程的离散

利用控制体积有限元方法离散化的数学模型由守恒方程式(1)和式(2)推导出来，改写之后的守恒方程为：

(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)，是通过对控制体积积分和应用高斯散度定理求得的，求得的表面上的积分相当于表面上的离散和，又有：

(11)

(12)

1.4 离散方程的解法

解离散方程的方法采用牛顿迭代算法。应用控制体积有限元方法解在网格每个节点上的守恒方程，能够得到一个2×N的非线性离散方程组：

F(u)=[Fw1(u)，Fe1(u)，Fw2(u)，Fe2(u)，…，FwN(u)，FeN(u)]T=0。

(13)

式中：Fwi(u)和Fei(u)分别为网格中第i个控制体积的液体离散守恒方程组和能量离散守恒方程组；解向量u=(X1,T1，X2，T2，…，XN,YN)，u包括网格中每个节点上成对的主要变量。

解线性方程组采用稳定双共轭梯度算法，该算法常用于解决大规模矩阵方程组[7]；这种算法经常与预处理算法ILU(0)一起使用[8]，这样能大大提高系统的运算速度[9]。解线性方程组时，首先要求出雅各比矩阵，雅各比矩阵是对积分项和通量项进行求导得出的，计算导数的方法为一阶有限差分逼近，可用公式(14)表示。

(14)

式中：ewj为与主要变量X有关的第j个单元的向量；τ为位移值，τ必须足够大以避免产生舍入误差，并足够小以避免产生导数的错误估计值，式(14)中τ=10-6。

为了尽可能的减小牛顿算法中的浮点误差，尤其是构造雅各比矩阵时产生的误差，非线性方程组可改写为：

(15)

(16)

式中：D1和D2为对角矩阵。

1.5 控制体积表面的通量计算

1.5.1 扩散项的计算

本文中对附属控制体积表面的扩散张量和二次变量进行计算的方法为平均法，即取2个节点上赋值的平均值。这2个节点构成的线段，是截断附属控制体积的网格边界(见图1)，图1中所标注的点ab和cd就是上述2个节点。

图1 扩散通量的计算方法

1.5.2 平流项和对流项的计算

计算平流项和对流项的方法为极限通量算法，这种算法能够有效避免自激振荡引起的高阶空间离散，从而消除了可能在解域中出现的基波或震荡。极限通量算法计算附属控制体积表面液体迁移率的公式为：

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

(17)

式中：σ为限制器函数，0≤σ≤2；当σ=0时，为迎风状态；当σ=1时，为平衡状态；当σ=2时，为顺风状态。r为流向指示器的比率。在P. Arminjon[10]的研究中，σ(r)表示为：

(18)

流向指示器的比率为：

r=I2/I。

(19)

式中：I为处于迎风节点和顺风节点中间的控制体积表面；I2为第一个迎风节点和第二个迎风节点中间的控制体积表面。

2 仿真结果与分析

试验中，把板材的立体模拟图划分为若干个三角棱柱体，每个子单元都能在三维坐标系中表达出来。图2是四分之一板材的网格模型，横断面上共有1022个子单元，径向的几何比率为1.05，整体总共有12270个网格。

图2 三维空间的板材有限元网格模型

图3、图4分别为板材(四分之一)干燥1 h 30 min时的仿真图和3 h 10 min时的仿真图。图3、图4中的木材含水率范围不同，是由于随着木材干燥时间的增加，板材中的水分不断地迁移并蒸发。图3、图4充分体现了具有非均匀结构的木材，在干燥时的特点。从图3、图4可见：木材中的水分，优先选择纵向流动，这是因为木材在纵向上的渗透率是最高的。

3 结束语

本文建立了木材干燥过程含水率分布的数学模型。首先，给出了木材干燥过程含水率空间分布的液体守恒方程和能量守恒方程；然后，利用控制体积有限元原理将方程离散化，解离散方程采用了牛顿算法，其中，容限采用10-8；解线性方程时，采用了稳定双共轭梯度算法和预处理算法ILU(0)。计算通量时，扩散张量和二次变量的计算应用了平均值法，平流项和对流项的计算应用了极限通量算法。最后，根据该数学模型，建立了木材干燥过程含水率分布的仿真图，验证了利用控制体积有限元方法分析木材干燥含水率空间分布规律的可行性和准确性。

图3 1 h 30 min时板材含水率分布情况

图4 3 h 10 min时板材含水率分布情况

[1] 朱政贤.干燥锯材最终含水率与木材平衡含水率的初步探索[J].东北林业大学学报，1986，14(2)：1-10.

[2] 孙丽萍.木材含水率在线检测融合体系及仿真技术研究[D].哈尔滨：东北林业大学，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] 李贤军，乔建政，蔡智勇，等.微波干燥与常规干燥中木材内含水率动态分布[J].中南林业科技大学学报，2009，29(6)：98-113.

[5] 徐兆军，丁建文，丁涛，等.基于断层扫描图像技术的木材纤维饱和点以上水分分布和迁移特性研究[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] 张恩泽，彭树生，何小祥，等.超松弛迭代-双共轭梯度在三维电磁问题有限元分析中的应用[J].淮阴师范学院学报：自然科学版，2005，4(4)：292-295.

[8] 李熙铭，欧阳丹彤，白洪涛.基于GPU的混合精度平方根共轭梯度算法[J].仪器仪表学报，2012，33(1)：97-104.

[9] 金巍巍，陶文铨，何雅玲.代数方程求解方法收敛速度比较及对算法健壮性的影响[J].西安交通大学学报，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月生，东北林业大学机电工程学院，硕士研究生。E-mail：1060061063@qq.com。

孙丽萍，东北林业大学机电工程学院，教授。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) 国家林业公益性行业科研专项(201304502)。

责任编辑：张 玉。

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.