在木材热加工过程中，只有准确预测温度随时间的变化规律以及温度在木材中分布梯度，才能有效地控制热应力与结构单元的变形。目前，计算温度场的方法主要有数值法和有限元法。由于木材是正交各向异性材料(李坚，2002)，数值法无法计算木材传热过程中瞬态温度场的变化，而有限元法以其几何边界适应性好、剖分单元灵活、计算精度高等优点，已成功应用于各种温度场的计算(冯武等，2003；张国新，2004；陈庆军等，2006)。木材在热处理时纵向温度场变化相同，而径向和弦向温度场呈梯度变化(俞自涛等，2006)，针对木材的热各向异性，本文按照二维温度场进行了瞬态计算。本文选取东北林业大学帽儿山实验林场落叶松(Larix gmelinii)为试材，应用二维双线性矩形单元形函数，通过Galerkin方法求解热扩散方程，导出木材热分析的有限元算法(Lewis et al., 1996)。根据上述理论推导，利用ANSYS软件仿真，求解木材中的温度分布，再加以验证研究。

1 材料与方法 1.1 材料

落叶松试材含水率为(12±0.5)%，密度为650 kg·m^-3，径木的小头直径为15~20 cm，选取规则排列、无暇疵试材采用刨切机进行加工，获取10 cm×10 cm×15 cm(径向×切向×纵向)长方体试样。

1.2 导热系数和换热系数的测定

采用NETZSCH公司推出的新型HFM436系列热流导热测定仪，采用热流法，将厚度一定的长方形样品插入2个平板间，设置一定的温度梯度，使用校正过的热流传感器测量通过样品的热流，传感器在平板与样品之间和样品接触，通过测量样品厚度、温度梯度与通过样品的热流来计算导热系数。经过多次有效测量，可得落叶松木材试样的弦向导热系数为k_x＝0.128 9 W·m^-1K^-1，径向导热系数为k_y＝0.139 9 W·m^-1K^-1。内置热源100 ℃，外表面裸露在空气中，设空气恒温为30 ℃，相应的自然热传导系数h_x ＝14.1 W·m^-2K^-1，h_y＝16.7 W·m^-2K^-1。

1.3 试样内部温度测定

在试样的芯材施加线热源，由于正方形截面的木材温度场具有对称性，只分析试材1/4就可获知整个截面温度场的分布情况。选择分析的面被分成8个单元，15个节点，在每个节点处打上直径与铂电阻温度传感器相当，深度位于木材纵向的中部且大于传感器的感温传感头的长度，埋设位置如图 1所示。通过微处理器控制的自动巡测电路，实时依次测量各个节点的温度，并与上位机串行通信，保存所测得的各点温度数据。

2 木材瞬态热传导分析的基本理论 2.1 木材热扩散方程的建立

温度场随时间发生变化的传热过程称为非稳态传热，这类传热按其过程特点分为周期性和非周期性2种。非周期性传热过程中，物体内的温度随时间升高或降低，经历一段时间后趋于周围介质的温度达到平衡状态，这类传热过程又称为瞬态传热。工程中绝大部分都是瞬态传热，稳态是瞬态的最后工作步，瞬态与稳态分析的主要区别是瞬态分析中的载荷是随时间变化的。木材在应用时通常被取为长方形板，可以认为木材沿长度方向每个端面的温度场的变化相同(Mackerle，2005；Younsi et al., 2006)，因此，可近似为二维温度场处理，在笛卡儿坐标系下对系统应用能量守恒原理，得到热扩散方程

(1)

式中：表示单位体产生的热量，k_x为弦向导热系数，k_y为径向导热系数，T为木材温度，ρ为材料密度，C_p为比热容。

2.2 矩形单元的公式推导

瞬态分析中使用的单元与稳态相同，二维双线性矩形单元(图 2)可用节点的温度和形函数表达(Moaveni，2003；张汝清，2004)，即：

(2)

式中：[T]为节点温度，[S]为形函数。

将热扩散方程应用Galerkin方法得：

(3)

式中:dA为矩形单元积分变量，该方程由4个主要的积分组成：

(4)

公式(4)里的第一项和第二项形式相同，可以应用链式法则将二阶化为一阶，以(4)式的第一项为例：

(5)

对该(5)式整理后得到：

(6)

将(6)代入(4)式，简化结果为：

(7)

同理得

(8)

其中

(9)

(10)

(7) 和(8)式右边第一项使用格林理论，积分项可以由绕单元边界的积分表示：

(11)

(12)

式中：τ表示单元的边界，θ表示单位法线的角度。

考虑传导边界条件及木材的各向异性：

(13)

(14)

式中：h_x, h_y分别为x，y方向上的换热系数，T_f为流动流体的温度为温度梯度。

将(13)、(14)式代入(11)、(12)式得：

(15)

(16)

由(9)和(10)式代入(7)和(8)式右边第二项整理得到：

(17)

同理可得：

(18)

接着计算热负荷项，即(4)式的第三项：

(19)

对(15)和(16)两式应用矩形单元不同边的边界条件积分，经整理得到双线性单元的热传导矩阵：

(20)

对(4)式的最后项积分可以得到总的热容矩阵[N]，最后将单元矩阵组合成总体矩阵，控制方程可以表示为，具体的流程将在数值计算中给出。

3 有限元数值模拟计算

由于选取试材为正方形截面，所以可以利用问题的对称性，对1/4面进行网格划分(图 1)，选择分析的面被分成8个单元、15个节点，每一个单元内部节点都使用i、j、m、n进行编号，见表 1。刚度矩阵和热容矩阵的求解方法相同，本文以刚度矩阵为例给出具体的计算方法。

将单元信息和测得的导热系数带入矩阵单元传导矩阵得：

任意边界条件对传导矩阵和负荷都会产生影响，受h_x的影响对单元(2)、(3)、(8)矩阵记录为

对单元(6)、(7)、(8)的受h_y的影响记录为

传导边界条件对单元(2)、(3)、(8)的热负荷矩阵产生的影响为：

传导边界条件对单元(6)、(7)、(8)影响为

将所有的单元矩阵组合起来，应用每个单元相邻的所给的节点信息，整体刚度矩阵为：

用同样的方法可以得到总的热容矩阵[N]，合成后的刚度矩阵和热容矩阵均是15×15维矩阵。显然，单元划分得越密，节点数目越大，刚度矩阵和热容矩阵也越大，具体可以从选定的温度场初值开始，选择一定的时间步长编程计算。利用MATLAB 7.0，编写相应的计算程序，进行有限元模拟计算，流程图如图 3所示。

4 试验验证及分析 4.1 试验验证

为了检验温度数值仿真结果的有效性，同时进行了对木材试样验证试验，通过微处理器依次测量各个节点的温度，分别获得了温度场在30、40和50 s后的实际测量值，试验测量值和理论计算值的比较结果见表 2、3和4。经过理论值与测量值的对比，发现数值模拟计算值在开始时会出现震荡，减少时间步长能使稳定性和精度有一定的提高，同时也要加密单元的划分，否则效果不明显，但是随着加热时间的加长，精度会逐步提高。木材热加工的过程中，纤维饱和点时的温度控制是关键，所以本研究有工程应用价值，误差可以控制在-2.0%~2.6%之间，可满足工程计算的精度要求。

4.2 有限元分析

采用ANSYS大型有限元仿真软件，选用四节点的四边形单元PLANN55用于二维木材温度场的建模。单元的每一个节点具有一维的自由度，传导和热流可以看作是施加在单元表面的表面负荷，单元的输出包括数据节点温度和单元数据，例如热梯度和热流量。通过设定单元参数(密度、径向切向导热、换热系数大小如前所述)、建模、网格划分、加载、求解，得到t＝50s时的温度场的云图如图 4所示，从图中可清楚看出，其温度分布与数值求解是相符的，而且x(切向)方向的温度梯度要稍微小于y(径向)的温度梯度，这是由于弦向传导系数k_x低于径向传导系数k_y。切向与径向导热系数的大小与木材射线的容积和早材、晚材密度的差异有关，图 5的温度梯度变化曲线验证了这一结论。

5 结论

本文对木材的相关材料热参数进行了测量，利用传热学基本原理和有限元基本方法推导出了基于木材瞬态温度场的有限元算法，采用MATLAB编写计算程序进行数值计算，并设计了木材温度场的测量试验，试验表明2种方法结果基本相符，可满足实际木材工业的需要。同时，又采用ANSYS有限元仿真软件对木材瞬态温度场进行了建模，进一步验证了有限元方法对于木材热传导问题的分析的可行性，仿真结果同时也说明木材各向异性对温度场分布有一定的影响。在今后的研究中，对于木材温度场的传递过程中对木材含水率、导热系数的影响以及木材的非均质性、非线性的特点有待于进一步深入讨论。