金属杆热传导的数值计算研究

康晓宇 吴惠 刘应开

(云南师范大学物理与电子信息学院 云南 昆明 650092)(昭通师专物理系 云南 昭通 657000)

1 引言

随着计算机技术的发展和数值计算方法的改进，求解偏微分方程组等早期不能完成的任务现在逐步变为现实.在数值计算的各种软件和软件平台中，MATLAB简单易学且具有强大的可视化功能，为科学研究提供了极大的方便.特别是MATLAB中的PDE工具箱，无需编程则可直接求解具有特殊边界条件的导热问题，并将结果可视化，从而对物理过程和有关特征量的演化有一个直观的认识.

在热学中，热传导问题是一个基础性的课题，对其进行研究有助于理解热传导的基本规律，澄清基本物理概念.在数学物理方法中，热传导问题是由热传导方程来描述.热传导方程的求解因边界条件不同而难度各异.有的问题没有严格的解析解，只能求得数值解；有的则存在严格的解析解，其解的形式可能是积分形式.也可能是级数形式，从这些解中很难准确知道温度、热流密度等物理量在空间的分布和随时间的演化规律.本文利用MATLAB中的PDE工具箱求解两类一维热传导问题，它们分别属于第一类和第二类边界条件.通过求解两类问题中的温度场、热流密度分布与演化规律，揭示热传导的本质，从而澄清稳定态与平衡态的区别与联系.

2 一维热传导方程

均匀金属杆可看作一维传热系统.如图1所示，杆的左端作为坐标原点，则右端的坐标为L，杆上各点的温度T(x,t)满足偏微分方程

(1)

其中a=κ/cρ，κ，c和ρ分别为金属杆的热导率、比热容和密度.

图1 一维细杆热传导示意图

3 两类边界条件的一维热传导问题

3.1 第一类边界条件的一维热传导问题

长为L均匀铁制金属杆，杆的初始温度为20℃，其左端与0℃的恒温热源接触，右端与100℃恒温热源接触，杆的侧面与外界环境绝热.现讨论该金属杆的温度场分布、热流密度分布及其演化情况，其方程如下

3.2 第二类边界条件的一维热传导问题

长为L均匀铁制金属杆，初始时刻杆的左端温度为0℃，右端温度为100℃，杆上温度梯度均匀分布，整个金属杆与外界绝热.现讨论该金属杆的温度场分布、热流密度分布及其演化情况，其方程如下

其中a=κ/cρ，计算中铁制杆的密度，比热容和热导分别取ρ=7 800 kg/m3,c=500 J/(kg·℃),κ=48.6 W/(m·℃)[1].

3.3 求解的基本方法

问题的求解是基于有限元法[2]，其基本原理是把计算区域划分成一系列的三角形单元，每个单元上取一个节点，选定一个形状函数(抛物线型或双曲线型)，并通过单元中节点上的被求解变量值表示该函数.有限元法的最大优点是对不规则区域的适应性好，故用MATLAB方法求解的结果在边界上也比较精确[3].

MATLAB中的PDE工具箱在传热学中主要用于解如下两类偏微分方程

式中，u为求解变量，t为时间变量，d，a，f，c为常数或变量.

MATLAB中的PDE工具箱定义了两类边界条件，即

式中n为垂直于边界的单位矢量，h,r,q和t为常量或与u有关的变量[1，4].

利用PDE工具箱求解导热问题的一般步骤如图2所示.

图2 MATLAB的计算流程图

4 结果与讨论

下面我们用MATLAB中的PDE工具箱来求它们的数值解.

4.1 第一类边界条件的一维热传导问题中的温度场、热流密度分布与演化规律

将坐标轴的范围设为0≤x≤5，-1≤y≤1，求解区域为0≤x≤4,-0.2≤y≤0.2.

矩形的上下边界为Neuman边界条件，取q=0，g=0， 左右边界均为Dirichlet边界条件， 分别取

h=1，r=0和h=1，r=100.

方程的类型为抛物线类型，其参数分别取c=48.6，a=0，d=7 800×500，f=0，初始条件为u(t0)=20[2].

图3是经过5小时，20小时，250小时后第一类边界条件的金属杆在空间的温度分布云图、热流密度分布云图，这些图还显示了有限元求解的网格.

图3 经过不同时间后第一类边界条件的金属杆温度空间分布云图及热流密度分布云图

由图3(a1)可知，经过5小时后，金属杆两端温度梯度较大，中间温度梯度较小.图3(a2)为，随着时间的增加，金属杆不同部位的温度梯度差异缩小，但温度分布面仍为曲面；图3(a3)为，当时间进一步延长，温度分布面变得平直，说明温度梯度趋于相等，金属杆的温度在空间的分布形成稳定梯度，不同空间位置的温度分布不再随时间而发生明显的变化.

由图3(b1)可知，经过5小时后，金属杆不同部位热流密度的分布是不同的，低温端热流密度比较平稳，热流密度小，而2 m～4 m区域，由于温度从20℃～100℃变化，热流密度变化较大，大量的热由高温热源输运到金属杆中；图3(b2)所示，金属杆不同部位的热流密度分布存在一定差异，即随着时间的增加，空间各处热流密度差异减小；图3(b3)所示，表明经过足够长的时间金属杆各部位的热流密度趋于相等，对于某一给定区域来说，热流流入多少，就几乎流出多少，也就是说，经过足够长的时间后热流稳定地从高温热源输运到低温热源，这与温度在空间的稳定分布一致.

4.2 第二类边界条件的一维热传导问题的温度场、热流密度分布与演化规律

在使用PDE工具箱求解时，只需将上面求解过程中的左右边界条件改为Neuman边界条件，即取g=0，q=0和初始条件设为u(t0)=25×x即可，其余条件不变[2].

图4 经过不同时间后第二类边界条件的金属杆温度空间分布云图及热流密度分布云图

图4是经过5小时，20小时，250小时后第二类边界条件的金属杆在空间的温度分布云图、热流密度分布云图，这些图还显示了有限元求解的网格.

由图4(a1)可知，经过5小时后，金属杆的温度分布已发生变化，温度梯度不再是常量，空间的温度分布平面变弯曲；图4(a2)为，随着时间的增加，斜面变得较平坦，温度梯度进一步减小；图4(a3)为，随着时间的进一步增加，即经过250小时后，温度分布平面变得几乎是平行xOy平面，温度约为50℃，说明金属杆上各点温度近似相等，不存在明显的温度梯度.即经过足够长的时间后，金属杆上的温度为一个恒定的值50℃.

由图4(b1)可知，经过5小时后，金属杆的热流密度在杆的两端近似为零，而杆的中间部位的热流密度相当大并且近似为1 200 W/m2；由图4(b2)可知，经过20小时后，整个金属杆中的热流密度已经大大减小，即随着时间的增加，金属杆各部位的热流密度在减小；由图4(b3)可知，经过250小时后，金属杆的热流密度逐渐趋于零，也就是说，经过足够长的时间后金属杆内不存在能量的输入与输出，这与温度在空间的均匀平衡分布是一致的.

5 结论

综上所述，本文分别对第一类边界条件下和第二类边界条件下的金属杆的一维热传导问题进行了可视化分析，从它们的温度分布云图和热流密度云图可以得出：第一类边界条件下的金属杆经过足够长的时间后，将处在一个宏观性质不随时间变化的稳定状态，即空间各点的温度是定值，但仍有稳定的热流从高温热源流向低温热源，能量在空间各点没有堆积现象，是一个稳定态.而第二类边界条件下的金属杆经过足够长的时间后，金属杆上各点的温度均为50℃，不会随时间的而变化，而且空间各点的热流密度也趋于零，是一个平衡态.从二者的比较可知：

(1)有外界影响的情况下，经过足够长的时间后，金属杆内空间各点的温度有稳定分布，但热流不为零，属于稳定态.

(2)无外界影响的情况下，经过足够长的时间后，金属杆由非平衡态转化为平衡态，空间各点温度相等，也无热流通过.

二者比较更能清楚地区分稳定态与平衡态的区别与联系，有助于对稳定态和平衡态概念认识的深化.

参考文献

1 彭芳麟.数学物理方程的MATLAB解法与可视化.北京:北京师范大学出版社,2004

2 陶文铨.数值传热学(第2版) .西安:西安交通大学出版社,2001

3 李灿,高彦栋,黄素逸.热传导问题的MATLAB计算.华中科技大学学报(自然科学版),2002,30(9):91-93

4 徐梓斌,闵剑青.基于PDE tool的热传导数值计算.佳木斯大学学报(自然科学版),2006,24(4):270-272