二维导热物体温度场的数值模拟
一、物理问题
有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸如下图1-1所示，假设在垂直于纸面方向上用冷空气及砖墙的温度变化很小，可以近似地予以忽略。
在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每米长度上通过砖墙的导热量。
第一种情况：内外壁分别均匀维持在0℃及30℃；
第二种情况：内外壁均为第三类边界条件，且已知：
K
m K m W h C t K
m W h C t •=•=︒=•=︒=∞∞/35.0/93.3,10/35.10,302
22211λ砖墙导热系数
二、数学描写
由对称的界面必是绝热面，可取左上方的四分之一墙角为研究对象，该问题为二维、稳态、无内热源的导热问题。
控制方程：
02
222=∂∂+∂∂y t
x t
边界条件： 第一种情况：
由对称性知边界1绝热： 0=w q ； 边界2为等温边界，满足第一类边界条件： C t w ︒=0； 边界3为等温边界，满足第一类边界条件： C t w ︒=30。
第一种情况：
由对称性知边界1绝热： 0=w q ；
边界2为对流边界，满足第三类边界条件： )()(
2f w w w t t h n t
q -=∂∂-=λ； 边界3为对流边界，满足第三类边界条件： )()(2f w w w t t h n
t
q -=∂∂-=λ。
1
-1图2
-1图
三、方程离散
用一系列与坐标轴平行的间隔0.1m 的二维网格线将温度区域划分为若干子区域，如图1-3所示。
采用热平衡法，利用傅里叶导热定律和能量守恒定律，按照以导入元体（m,n ）方向的热流量为正，列写每个节点代表的元体的代数方程，
第一种情况： 边界点：
边界1（绝热边界）：
5~2)2(4
1
1,11,12,1,m =++=
+-m t t t t m m m ， 11~8)2(4
1
1,161,16,15,16=++=+-n t t t t n n n n
， 边界2（等温内边界）： 7,16~7;7~1,6,0,=====n m n m t n
m
边界3（等温外边界）： 12,16~2;12~1,1,30,=====n m n m t n m
内节点：
11
~8,15~6;11~2,5~2)(41
1,1,,1,1,====+++=
-+-+n m n m t t t t t n m n m n m n m n m
第二种情况 边界点：
边界1（绝热边界）： 5~2)2(4
1
1,11,12,1
,m =++=+-m t t t t m m m ， 11~8)2(4
1
1,161,16,15,16=++=+-n t t t t n n n n ，
边界2（内对流边界）：
6~1)
2(2221
11,61,6,5,6=++++=
∆∆-+n Bi t Bi t t t t n n n n ，
3
-1图
16~7)
2(22211
17,17,18,7
,=++++=
∆∆-+m Bi t Bi t t t t m m m m ，
边界3（外对流边界）：
11~1)2(22222
21,11,1,2,1=++++=
∆∆-+n Bi t Bi t t t t n n n n
，
16~2)
2(22222
212,112,111,12,=++++=
∆∆-+m Bi t Bi t t t t m m m m ，
内角点： )
3(22)(21116,67,78,67,57,6+++++=
∆∆Bi t Bi t t t t t
外角点：
)
1(222
211,112,212
,1+++=∆∆Bi t Bi t t t
内节点：
11
~8,15~6;11~2,5~2);
(41
1,1,,1,1,====+++=-+-+n m n m t t t t t n m n m n m n m n m
（10,22121==∆=∞∆t t x
h Bi λ
；30,21212==∆=
∞∆t t x
h Bi λ
）
四、编程思路及流程图
编程思路为设定两个二维数组t(i,j)、ta(i,j)分别表示本次迭代和上次迭代各节点的温度值，iter （实际编程时并未按照此名称来命名迭代步长）表示迭代进行的次数, 1Q 、2Q 分别表示外边界、内边界的散热量。
开始时，给t(i,j ）、ta(i,j) 赋相同的初始值，t(i,j)根据内节点和各边界节点的离散方程进行迭代，迭代后比较t(i,j)、ta(i,j)各个节点之间温度之差，若两个温度之差小给定的精度，则此时迭代完成，t(i,j)就是所求的温度场分布，若两温度之差不满足精度要求，则将t(i,j ）的值赋给ta(i,j)，t(i,j)继续迭代，直到二者各个点的温度之差满足精度要求，记下此时的迭代次数，并根据所得到的温度场分布计算内外边界上散热量以及偏差。
图
4-1
五、结果讨论
本次实验的实际边界条件为等温边界条件，因此以第一种情况为依据进行实验编程。
程序（见附注）运行结果如图1-5所示。
使用MATLAB7.0软件将各节点温度作为原始数据做出四分之一墙角的温度场分布图
5-1图
6-1图
将实验测量结果和数值计算结果进行比较如表1-1
实验测量值 数值计算值 四分之一墙角外侧散热量W / 60.712 60.4286728 四分之一墙角内侧散热量W / 60.341 60.4287262 单位高度墙壁总散热量W /
242.106 241.715 相对平衡偏差%/
0.613
0.0000884
对比两种结果分析数值计算方法得到的散热量比实验测得的散热量相近，而前者的相对
平衡偏差比较小。
前者的误差主要来自于迭代次数和网格数的划分，但也可能存在网格独立解，可以进一步细化网格进行求解观察数值模拟结果；后者的误差主要来自于电阻的阻值误差和读数仪表的基本误差，提高标准电阻的准确性和读数仪表的精度等级可以减小实验误差。
7
-1图。