
文章信息
- 赵景尧, 付宗营, 宦思琪, 蔡英春
- Zhao Jingyao, Fu Zongying, Huan Siqi, Cai Yingchun
- 有限差分法逆求木材导热系数
- Inverse Determination of Thermal Conductivity of Wood Using Finite Difference Method
- 林业科学, 2015, 51(4): 134-140
- Scientia Silvae Sinicae, 2015, 51(4): 134-140.
- DOI: 10.11707/j.1001-7488.20150417
-
文章历史
- 收稿日期:2014-05-21
- 修回日期:2014-06-23
-
作者相关文章
导热系数是木材热学性质(热容量、热传导、热扩散与热膨胀)中最重要的物理参数之一(高瑞堂等,1985;Yang,2001),也是木材热处理(如木材干燥、人造板热压等)过程中的重要工艺参数,尤其在木材干燥方面,直接影响干燥过程中木材升温速率、内部温度分布以及干燥热质迁移模型的准确性等(俞昌铭,2011;Walter et al.,2011)。因此,准确测算木材导热系数,对了解干燥过程中木材内部温度变化规律及科学合理地实施干燥工艺意义重大。
木材是一种天然高分子多孔材料,其内部的复杂化学成分和结构给学者研究木材导热性质带来极大困难。侯祝强(1992)、陈瑞英等(2005)、Hunt等(2006)、俞自涛等(2007)、蔡从中等(2009)、林铭等(2013)采用不同理论与方法对木材内部结构模型化,根据类比原理与数学分析法建立了导热系数的理论表达式并探讨了导热系数的影响因素。实际上导热系数受含水率、温度影响较大,仅由理论导出的平均导热系数来研究木材升温速率与内部温度变化必将与实际产生较大误差。因此,建立导热系数依据含水率、温度变化的回归方程更具实际意义(MacLean,1941;Kollmann et al.,1968;Siau,1984),而导热系数的精确测算则是建立回归方程的重要前提。
导热系数的试验测量方法主要有稳态与瞬态法(俞自涛,2005;Dupleix et al.,2013),近年来,有限差分逆求法也被广泛应用于木材导热系数测算的研究中(Yeung et al.,1996;Chen et al., 1996;Cai et al.,2006)。与传统商用平板导热仪器相比,有限差分法快速、简单、经济,且不受试样尺寸限制;同时,可以测算试样内部任意层位置的导热系数。其主要思想是:在特定含水率下,测量升温过程中木材内部不同位置温度随时间的变化,据各层温度变化运用有限差分法逆求导热系数。但上述相关研究仅探讨了逆求某一含水率时单一温度段的各层间导热系数,没有涉及温度、含水率对导热系数的影响;且边界条件不合理,如Cai等(2006)运用此方法测算冰冻材的导热系数,计算过程中边界节点方程仅1阶展开,导致边界计算精度与内部不匹配,影响导热系数的求解精度。
鉴于此,本文以兴安落叶松(Larix gmelinii)木材为对象,采用有限差分逆求法,合理化边界条件后编程求解其在不同含水率、不同温度段变化时的径向导热系数,并分析导热系数随含水率、温度变化的规律,以期为后续建立导热系数随含水率、温度变化的二元回归方程提供基础。
1 材料与方法试材为采自黑龙江省东方红林业局的人工林兴安落叶松。初选3根不同初含水率(约60%,40%,30%),无开裂、腐朽、变形等明显缺陷,规格为4 200 mm(长)× 225 mm(宽)×42 mm(厚)的锯材,4面刨光并锯解成规格为600 mm(长)×220 mm(宽)×40 mm(厚)的3组弦切板试样(图 1a,每组6块,共计18块),各组试样平均含水率分别为 MC#1=28.3%,MC#2=41.2%,MC#3=62.3%。用环氧树脂和铝箔隔热卷材对各试样端面及侧面封闭后再用保鲜膜全面密封,以保证试验时热量只沿厚度方向传递并防止水分流失(试验前后水分流失不足1.2%)。由于温度巡检仪测点数量限制,每次仅对2块试样进行试验(每组在相同条件下进行3次),其余试样放入冷库中备用。3组不同平均含水率、初始温度为室温的试样分别在50,70,90 ℃的环境温度(湿度均为100%,风速2 m·s-1)下试验并测算。为减小材性对试验结果的影响,每次试验结束待试样温度降到室温后,再置入下一水平温度环境中试验。试验前,在试样侧面用直径为1 mm的钻头钻孔(具体位置见图 1b)并称重,然后将T形热电偶埋入试样中,再将试样放入温度已经稳定为设定值的DS-408型恒温恒湿箱中加热升温。使用NEC Remote Scanner Jr. DC3100多点信号巡检仪通过埋入木材的T形热电偶检测材温(每分钟记录1次)。自冷库中取出的试样,先室温解冻并稳定至室温后再进行上述操作。试验结束后,对试样按GB/T 1933—2009《木材密度测定方法》测量基本密度,采用绝干称重法检测含水率。试样初始温度(室温)为16~18 ℃。
![]() |
图 1 锯解示意与热电偶测温部位分布 Fig. 1 Sawing diagram and thermocouple arrangement |
模型建立前,提出合理的假设条件是必要的。由于试样厚度小于宽度、远小于长度,且其端口和侧面都进行了隔热处理,因此,假定温度仅沿厚度方向传递;同时,假定木材结构、物理化学特性在厚度方向上是对称的,厚度方向上的中心层面即为绝热面;假定试样初始温度分布均匀;假定升温阶段试样含水率保持恒定,即试样表面只存在热交换;假定试样外部环境不随其热量变化而变化,即环境所有状态参数是恒定的。
试样内部温度控制方程(逆求导热系数方程)如下:
${\rho _{\text{d}}}{c_{\text{p}}}\frac{{\partial T}}{{\partial \tau }} = \frac{\partial }{{\partial x}}\left({\lambda \frac{{\partial T}}{{\partial x}}} \right)$ | (1) |
显然,在式(1)中,当基本密度ρd、比热cp、各层温度$\partial T/\partial x$、不同时刻温度$\partial T/\partial x$变化已知的情况下,导热系数可解。
初始条件:
$T = T\left({0,x} \right)$ | (2) |
边界条件:
木材表面处(x=L,L为木材厚度的一半)的边界条件:
$-\lambda \frac{\partial T}{\partial x}=h\left( {{T}_{\operatorname{s}}}-{{T}_{\text{e}}} \right)$ | (3) |
木材厚度中心层(x=0)的边界条件:
$\frac{{\partial T}}{{\partial x}} = 0$ | (4) |
运用有限差分法对式(1)(逆求导热系数方程)进行求解。首先,对空间区域进行离散化,将区域[0,L]分成J个相同的单元,节点依次为j=1,…,J,共J个,空间步长记作Δx,xj=jΔx(图 2);将时间区域分成I个单元,节点依次为i=1,…,I,共I个,时间步长记作Δτ,τi=iΔτ。其次,写出所有方程的有限差分格式。
![]() |
图 2 试样与离散示意 Fig. 2 Diagram of specimen and discretization a.试样Specimen; b.沿径向空间离散 Space discretization in radial direction. |
1)内部节点差分方程(0 < j < J):
$\begin{gathered} 4\rho {c_{\text{p}}}\frac{{{{\left({\Delta x} \right)}^2}}}{{\Delta \tau }}\left({T_j^{i + 1} - T_j^i} \right)= \lambda _{j - 1}^i\left({T_{j - 1}^i - T_{j + 1}^i} \right)+ \hfill \\ 4\lambda _j^i\left({T_{j + 1}^i - 2T_j^i + T_{j - 1}^i} \right)+ \hfill \\ \lambda _{j + 1}^i\left({T_{j + 1}^i - T_{j - 1}^i} \right)\hfill \\ \end{gathered}$ | (5) |
2)为保证模型精度,表面节点采用2阶差分方程(j=J):
$\begin{gathered} \rho {c_{\text{p}}}\frac{{{{\left({\Delta x} \right)}^2}}}{{\Delta \tau }}\left({T_j^{i + 1} - T_j^i} \right)= \lambda _{j - 1}^i\left({T_{j - 1}^i - T_j^i} \right)+ \hfill \\ h\left({{T_{\text{e}}} - T_j^i} \right)\Delta x \hfill \\ \end{gathered}$ | (6) |
3)心层节点差分方程(j=1):
$\rho {c_{\text{p}}}\frac{{{{\left({\Delta x} \right)}^2}}}{{\Delta \tau }}\left({T_j^{i + 1} - T_j^i} \right)= 2\lambda _j^i\left({T_{j + 1}^i - T_j^i} \right)$ | (7) |
令矩阵d=(λj)J×1,矩阵b=(bj)J×1,矩阵a=(aj,k)J×J,以上有限差分方程即可表示为ad=b,如下式所示:
$\begin{gathered} {\left(\begin{gathered} {a_{1,1}}{\text{ }}0{\text{ }}0{\text{ }} \cdots {\text{ }}0 \hfill \\ {a_{2,1}}{\text{ }}{a_{2,2}}{\text{ }}{a_{2,3}}{\text{ }} \cdots {\text{ }}0 \hfill \\ \vdots {\text{ }}{a_{3,2}}{\text{ }} \ddots {\text{ }}{a_{3,4}}{\text{ }} \vdots \hfill \\ 0{\text{ }} \vdots {\text{ }}{a_{j,j - 1}}{\text{ }}{a_{j,j}}{\text{ }}{a_{j,j + 1}} \hfill \\ 0{\text{ }}0{\text{ }}0{\text{ }} \cdots {\text{ }}{a_{J,J}} \hfill \\ \end{gathered} \right)_{J \times J}}{\left(\begin{gathered} \lambda _1^i \hfill \\ \lambda _2^i \hfill \\ \vdots \hfill \\ \lambda _{J - 1}^i \hfill \\ \lambda _J^i \hfill \\ \end{gathered} \right)_{J \times 1}} = \hfill \\ {\left(\begin{gathered} b_1^i \hfill \\ b_2^i \hfill \\ \vdots \hfill \\ b_{J - 1}^i \hfill \\ b_J^i - h\left({{T_{\text{e}}} - T_J^i} \right)\Delta x \hfill \\ \end{gathered} \right)_{J \times 1}} \hfill \\ \end{gathered}$ | (8) |
${b_j} = \rho {c_{\text{p}}}\frac{{{{\left({\Delta x} \right)}^2}}}{{\Delta \tau }}\left({T_j^{i + 1} - T_j^i} \right),j = 1,j;$ | (9) |
${b_j} = 4\rho {c_{\text{p}}}\frac{{{{\left({\Delta x} \right)}^2}}}{{\Delta \tau }}\left({T_j^{i + 1} - T_j^i} \right),1{\text{ < }}j{\text{ < }}J$ | (10) |
${a_{j,j}} = 2\left({T_{j + 1}^i - T_j^i} \right),j = 1$ | (11) |
${a_{j,j - 1}} = T_{j - 1}^i - T_{j + 1}^i,1{\text{ < }}j{\text{ < }}J$ | (12) |
${a_{j,j}} = 4\left({T_{j + 1}^i - 2T_j^i + T_{j - 1}^i} \right),1{\text{ < }}j{\text{ < }}J$ | (13) |
${a_{j,j + 1}} = T_{j + 1}^i - T_{j - 1}^i,1{\text{ < }}j{\text{ < }}J$ | (14) |
${a_{j,j}} = T_{j + 1}^i - T_j^i,j{\text{ = }}J$ | (15) |
显然,矩阵a为“三对角矩阵”,因此可通过追赶法(TDMA,tridiagonal matrix algorithm)求出(Kreith et al.,2010;王正林等,2009)。通过Matlab2010b软件对上述矩阵式编程并运行解出导热系数(矩阵d)。
3 结果与讨论图 3为表 1所示组别#1(平均含水率MC为28.3%的试样分别在50,70,90 ℃的环境温度下)升温过程中不同深度位置(图 1)实测温度随时间的变化曲线,其中,各时刻温度为6块试样同部位温度均值。从图中可清晰看出:升温过程中试样各层温度随时间变化相似,即初期升温较快,随时间的延续变化渐缓,到后期逐渐趋于一致达到目标温度。
![]() |
图 3 试样升温过程各层温度随时间变化曲线(组别#1) Fig. 3 Variation of temperature as a function of time during heating(group#1) |
![]() |
以T1(距表层20 mm的心层温度)与T5(距表层4 mm)为例,初始阶段温差逐渐增大,3个温度水平下均在τ=30 h左右温差达到最大,分别为8.1,13.1,17.9 ℃,随后此温差逐渐变小。
图 4显示由矩阵关系式(8)编程计算得出的各组别不同厚度位置(x=4,8,12,16,20 mm)与时刻(τ=10,25,50,70,100 min)的导热系数(每组6块试样同位置同时刻的导热系数均值)分布,显示了导热系数的较大变异性。以图 4组别#1-50(MCave=28.3%,T=50 ℃)的导热系数分布为例,同时刻不同位置导热系数、同位置不同时刻导热系数都存在不同程度差异,且变化无明显规律。导热系数这2项差异中的前者,主要由不同部位的密度、构造、树脂含量等差异引起(Gu et al.,2005;2006;Vay et al.,2013;成俊卿,1985);而后者,可能是由于试样在加热升温过程中其内部水分产生某种程度迁移引起,而式(1)中并未考虑试样含水率分布的变化,因而产生了误差。总之,前者反映了客观存在即材性差异的影响,后者则显示了试验及求解误差即模型本身的影响,也表明了后续改进研究的必要,如考虑木材内部结构、含水率分布的影响等。
![]() |
图 4 各组别下不同位置与时刻的导热系数分布 Fig. 4 Variation of thermal conductivities with the changes of depths and moments for different groups |
表 1为试样不同平均含水率时在不同温度环境下的平均导热系数及其标准差。从表 1可以看出,总平均导热系数为0.106 1 W·m-1K-1,标准差为0.010 8,与按Gu(2005)几何模型公式计算得出的平均导热系数0.110 9 W·m-1K-1、林铭(2013)的理论计算导热系数0.125 2 W·m-1K-1接近,说明由本模型计算得出的平均导热系数有效。
表 2为含水率与温度对导热系数影响的方差分析。从表 2可以看出,含水率、温度对导热系数均存在显著性影响(P < 0.01),其中,含水率的影响高(FMC=126.942 1>FT=99.008 3);二者交互作用亦对导热系数存在显著影响。将图 4所显示的导热系数按含水率和温度分组求均值(同一温度同一含水率的所有不同时刻、不同层的导热系数均值)并绘制图 5。
![]() |
图 5和表 1表明:1)在试样温度一定时,导热系数随含水率的升高而增大(例如在90 ℃环境温度下,λ#1-90=0.099 9 W·m-1K-1,λ#2-90=0.112 9 W·m-1K-1,λ#3-90=0.139 8 W·m-1K-1),原因为,水分子的导热能力(λ≈0.6 W·m-1K-1)约是木材分子的3倍、空气分子的23倍,随着木材内水分的增加,孔隙内空气被水分替代比例增大,因而产生上述结果;2)在试样含水率恒定时,导热系数亦随温度升高而增大(例如组别#3,即MCave=62.3%状态下,导热系数随温度升高而增大,λ#3-50=0.113 1 W·m-1K-1,λ#3-70=0.121 0 W·m-1K-1,λ#3-90=0.139 8 W·m-1K-1),这是因为,随着温度的升高,构成木材的分子及其内部水分子的热运动加剧,且孔壁间的辐射能增大,导致导热系数增大;3)随着含水率的增高,温度对导热系数的影响更显著,MCave=62.3%状态下,温度每升高1 ℃,导热系数随之增大约0.067%,而MCave=41.2%,28.3%状态下则分别增大约0.044%,0.047%,可能为温度对水分子的影响大于对空气分子的影响所致;4)含水率在28%(约为纤维饱和点)以上时,随着温度的升高,含水率对导热系数的影响程度增大: T=90 ℃环境下,含水率每增大1%,导热系数随之增大约0.117%,而T=70,50 ℃下则分别增大约0.091%,0.094%;低含水率范围内,温度对含水率的影响程度则无明显作用。原因可能为,纤维饱和点之上,温度对细胞腔内自由水分子热运动程度的影响比对细胞壁内结合水的影响大,致使产生本项结果。
![]() |
图 5 不同组别平均导热系数分布 Fig. 5 Average of thermal conductivities at various groups |
1)兴安落叶松径向导热系数,尽管计算值沿厚度方向存在一定波动性,但其平均值0.106 1 W·m-1K-1(标准差0.010 8)符合实际要求,且与理论计算值(0.110 9,0.125 2 W·m-1K-1)较为接近。
2)含水率、温度对导热系数影响显著,前者影响高于后者;同时,含水率、温度共同作用对导热系数亦存在显著影响;导热系数随含水率的升高而增大,随温度的升高亦增大。
3)材性、含水率分布对导热系数的影响较大。
[1] |
成俊卿. 1985.木材学. 北京: 中国林业出版社. (Cheng J Q. 1985. Wood Science. Beijing. China Forestry Publishing House.[in Chinese])( ![]() |
[2] |
蔡从中,温玉锋,朱星键,等. 2009.木材导热系数的支持向量回归预测.重庆大学学报,32(8): 960-964. (Cai C Z, Wen Y F, Zhu X J, et al. 2009. Wood thermal conductivity prediction by using support vector regression. Journal of Chongqing University, 32(8):960-964[in Chinese]).( ![]() |
[3] |
陈瑞英,谢拥群,杨庆贤,等. 2005.木材横纹导热系数的类比法研究.林业科学,41(1): 123-126. (Chen R Y, Xie Y Q, Yang Q X, et al. 2005. Study on wood thermal conductivity in transverse direction by analogism. Scientia Silvae Sinicae, 41(1):123-126[in Chinese]).( ![]() |
[4] |
高瑞堂,刘一星,李文深,等. 1985. 木材热学性质与温度关系的研究.东北林业大学学报,13(4): 22-27. (Gao R T, Liu Y X, Li W S, et al. 1985. Study on the relationship between the thermal properties of wood and the temperature. Journal of Northeast Forestry University. 13(4):22-27[in Chinese]).( ![]() |
[5] |
侯祝强. 1992.木材导热系数的研究.林业科学,28(2): 153-160.[in Chinese])(Hou Z Q. 1992. Studies on thermal conductivity of air-dried wood. Scientia Silvae Sinicae, 28(2):153-160.[in Chinese])(![]() |
[6] |
林铭,谢拥群,饶久平,等.2013.木材横纹导热系数一种新的表达式的推导及与实测值比较.林业科学,49(2): 108-112. (Lin M, Xie Y Q, Rao J P, et al. 2013. Derivation of a new expression for wood transverse thermal conductivity and comparison to experimental value. Scientia Silvae Sinicae, 49(2):108-112.[in Chinese])( ![]() |
[7] |
王正林,龚纯,何倩. 2012.精通MATLAB科学计算.北京: 电子工业出版社. (Wang Z L, Gong C, He Q. 2012. Master MATLAB scientific computation. Beijing. Publishing House of Electronics Industry.[in Chinese]) |
[8] |
俞昌铭. 2011.多孔材料传热传质数值分析. 北京: 高等教育出版社(Yu C M. 2011. Numerical analysis of heat and mass transfer for porous materials. Beijing:Higher Education Press.[in Chinese])(![]() |
[9] |
俞自涛.2005.着火前木材传热传质过程的试验和理论研究.杭州:浙江大学博士学位论文. (Yu Z T. 2005. Experimental and theoretical investigation on heat and mass transfer in wood prior to fire. Hangzhou:PhD thesis of Zhejiang University[in Chinese]).( ![]() |
[10] |
俞自涛,胡亚才,田甜,等.2007.木材横纹有效导热系数的分形模型. 浙江大学学报:工学版,41(2): 351-355(Yu Z T, Hu Y C, Tian T, et al. 2007. Fractal model for predicting effective thermal conductivity perpendicular to fibres of wood. Journal of Zhejiang University:Engineering Science, 41(2): 351-355.[in Chinese])(![]() |
[11] |
Cai L P, Garrahan P. 2006. Inverse determination of thermal conductivity in frozen lumber. Wood Sci Technol,40(8): 665-672.(![]() |
[12] |
Chen H T, Lin J Y, Wu C H, et al. 1996. Numerical algorithm for estimating temperature-dependent thermal conductivity. Numerical Heat Transfer, 29(4):509-522.(![]() |
[13] |
Dupleix A, Kusiak A, Hughes M, et al. 2013. Measuring the thermal properties of green wood by the transient plane source(TPS)technique. Holzforschung, 67(4):437-445.(![]() |
[14] |
Gu H M, Zink-Sharp A.2005. Geometric model for softwood transverse thermal conductivity. Part I. Wood and Fiber Science, 37(4):699-711.(![]() |
[15] |
Gu H M, Hunt J F. 2006. Two-dimensional finite element heat transfer model of softwood. Part II. Macro structural effects. Wood and Fiber Science, 38(4): 599 -608.(![]() |
[16] |
Hunt J F, Gu H M. 2006. Two-dimensional finite element heat transfer model of softwood. Part I. Effective thermal conductivity. Wood and Fiber Science, 38(4): 592-598.(![]() |
[17] |
Kreith F, Manglik R, Bohn M. 2010. Principles of heat transfer. Sixth edition. Thomson Learning. Australia, Canada, A51-A53.(![]() |
[18] |
Kollmann F F P, Cote Jr W A. 1968. Principles of Wood Science and Technology. Vol. I. Solid Wood. Springer-Verlag, Berlin.(![]() |
[19] |
MacLean J D. 1941.Thermal conductivity of wood. Heating, Piping and Air Conditioning,13:380-391.(![]() |
[20] |
Pordage L J, Langrish T A G. 1999. Simulation of the effect of air velocity in the drying of hardwood timber. Drying Tech, 17(1/2):237-255.(![]() |
[21] |
Steinhagen H P, Lee H W. 1988. Enthalpy method to compute radial heating and thawing of logs. Wood Fiber Sci, 20(4):415-421.(![]() |
[22] |
Siau J F. 1984. Transport processes in wood. Berlin, Heideberg, New York:Spring-Verlag, 245.(![]() |
[23] | Sonderegger W, Hering S, Niemz P. 2011. Thermal behaviour of Norway spruce and European beech in and between the principal anatomical directions. Holzforschung, 65(3):369-375. |
[24] |
Vay O, Obersriebnig M, Müller U, et al. 2013. Studying thermal conductivity of wood at cell wall level by scanning thermal microscopy(SThM). Holzforschung, 67(2):155-159.(![]() |
[25] |
Yang Q X.2001.Theoretical expression of thermal conductivity of wood. Journal of Forestry Research, 12(1):43-46.(![]() |
[26] |
Yeung W K, Lam T T. 1996. Second-order finite difference approximation for inverse determination of thermal conductivity. Int J Heat Mass Transfer, 39(17):3685-3693.(![]() |