1、热传导与热辐射大作业报告热传导与热辐射大作业报告目录一、作业题目 .- 1 -二、作业解答 .- 2 -个人感想 .- 17 -附件.计算中所用程序 .- 18 -热传导与热辐射大作业报告- 1 -一、作业题目一矩形平板 , ，内有均匀恒定热ax0by源 ，在 及 处绝热，在 及 处0gaxby保持温度 ，初始时刻温度为 ，如右图 1 所示：1T0T1、求 时，矩形区域内的温度分布0t的解析表达式；yx,2、若 ， ， ， ， ，热传导系ma8b1230mWgK6T0120数，热扩散系数 。请根据 1 中所求温度分布用KWk0.12.8sMATLAB 软件绘出下列结果，加以详细物理比较和分析：
2、(a) 300s 内，在同一图中画出点 、)4,0(、 、 、 （单位：m）温度随时间的变化；)8,0(,6)0,12(6,9(b) 200s 内，画出点 、 、 、 、 （单位：),18,(12,6),(6,9(m）处，分别沿 x、y 方向热流密度值随时间的变化；(c) 画出 时刻区域内的等温线；sst 502750、(d) 300s 内，在同一图中画出点 （单位：m）在 分别等于,90g， ， 情况下的温度变化；31W323m(e) 300s 内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为 、 和 的温度、热流密度变化；K5.0W0.1K5.1(f) 300s 内，比较
3、点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为 、 和 的温度、热流密度变化；s24.s28.s2.3、运用有限差分法计算 2 中(b)、(d) 和(e) ，并与解析解结果进行比较，且需将数值解与解析解的相对误差减小到 1以下；4、附上源程序和个人体会；以报告形式整理上述结果，用 A4 纸打印上交。热传导与热辐射大作业报告- 2 -二、作业解答1、求 时，矩形区域内的温度分布 的解析表达式；0t tyxT,解答：我们令 ，则可以得到一个方程和边界条件：1T（1-1）tkgyx02T.,da0x，,ydb0，t1，T将上式分解为一个 的稳态问题：)yxs，（1-2）0kgx2ss
4、2.,dsa0sx，,ydsb0s，和一个 的其次问题：),htyx（1-3）t122Th热传导与热辐射大作业报告- 3 -.0,dxhah，0,ydbh，其中 ),(),),(),h yxfyxFs（ 则原问题的解根据下式求得：（1-4）),(),(), tyxtyxhs（发热强度为常数的特解可从表 2-4 中查的，则新变量可定义为：)yxs，（1-5）Axkg20s )y)x，（，（ 将（1-5 ）带入（ 1-2）整理得到：（1-6） ，（ 0y),(x),22x byax0,.,0daa2xAkg，0,ydbxkgy2，若令常数 ，则上式可以变为：0aA（1-7） ，（ 0y),(x),
5、22xbyax0,.,0dax，热传导与热辐射大作业报告- 4 -0,ydbxf)(，其中 )220akg假定 可以分离出如下形式：),yx（1-8）)(),yYxX（对应于 的分离方程为：()yxX和（1-9）0)(d22x,0xaX，（1-10）0)()(d22yY,0ybxfY，)(在 中特征值问题的解可以直接从表 2-2 第 6 条中得到，)(xX只需要用 a 代替 L，（1-11）xXmmcos),（1-12）a2)1N（是下面方程的正根：m（1-13）0acosm方程（1-10）的解可以取为热传导与热辐射大作业报告- 5 -（1-14）yhYmmcos,）（的完全解由下式组成：),
6、yx（1-15）1cosh),(mmxyCyx此式满足热传导问题（1-7）及三个齐次边界条件，其中，系数 可以根据方程的解还应满足非齐次的边界条件来决定。mC利用 的边界条件可得：by（1-16）1cosh)(mmxbCxf a0利用函数 的正交性可以求得系数 ，cos mC（1-17）amdxfbN0m)(cossh）（式中： 3m02a01800 sin)(cos)(cos kagdxakgxdxfmm 将这个表达式带入式（1-15） ，其中范数 在前面已经）（N给出，解得结果为（1-18）301msincoshbco2), mm kagxyayx （则： Axkg20s )y)x，（，（
7、 )(2singcoshbco12 2030m xakxyamm 热传导与热辐射大作业报告- 6 -（1-19）假定 分离成如下表达式),(htyx（1-20）)()(,(hyYxXt对应于函数 和 的分离方程为)xXy： （1-21）)(x 0)(d22xax0,xaX，： （1-22）)(yY0)()(d22yYby,0yYb，的解为：)(t（1-23）tvt)(2ne)上述问题的完全解为：（1-24）1)(h ),(,),(2nmn vnt yYxXeCtyxv其中 0xa,0yb。当 t=0 时，上式变为：（1-25）1h ),(,),(mnvnvyYxXCyx其中 0xa,0yb。确
8、定未知系数 的方法是，在上式两边逐项用如下算子作mn热传导与热辐射大作业报告- 7 -运算：及a0n),(dxXb0v),(dyY并利用这些函数的正交性，得到：（1-26）a0b vnnv ),(,(),1 dyxyxNChvn）（）（最终得到问题的解为：1 vn)(h ),(,1),(2mn vnt yYxXNetyxvn ）（）（1-27）a0b ),(,(dyxyYXhv），（式中出现的特征函数，特征值及范数可以从表 2-2 中直接查得：（1-28）xXnncos)(，（1-29）a2)(1nN且 为如下方程的正根：m（1-30）0cosna满足特征值问题的函数 对应于表 2-2 中的第
9、 6 条，得),(yYn到：（1-31）yYvcos),(v（1-32）b2)(1vN且 是如下方程的正根：n0cosv最后得到：热传导与热辐射大作业报告- 8 -1)(h cos14),(2nmvnt yxabetyxv （1-33）a0b n),(cosdyxyhv令 vna0 ),(cosxIhb 1 20300 )(sincoscok2m mmv dyxakgxybagTy a0 10cosbvndyxyx）（ 3 sincoshco2m mmdyxabakga0 20 )(csbvn yxyx其中令 vnbvn baTdTI sis10a0 101 ）（）（ 令 a01 302 si
10、coshco2cosbm mmvn dyxybakgyxI 10b 3 ssshim vnxkg10b030 cohcocoin2a mvmnm dyydba 根据余弦函数的正交性，只有当 m=n 时积分才不为 0，故上式可以化为： a nvnnnm dyydxbakg0b030 coshcoscoshi2 再令 anaxI021 2sbhdyy vnvnv sicocohnb02 所以 )(iscs)1( 23n213n2 vmVaIbakgI 热传导与热辐射大作业报告- 9 -令 30a0 203 sin)(cos vbvn kabgdyxakgyxI 所以 )()(1220010n321
11、 vnnvTII由（1-4 ） 、 （1-19 ）及（ 1-33）可知 )(2cossh2),( 201coshin30 xakgxytyxT mmbakgm 1v 220010n)( cos)()(si42nm vnvnnvmt yxkgTbaeabv 1T。以上是解析解的全过程，具体值的计算采用 MATLAB 编程计算求取。2、若 ， ， ， ， ，热ma18b2301mWgK6T0120传导系数 ，热扩散系数 。请根据 1 中所求KWk0. 2.8s温度分布用 MATLAB 软件绘出下列结果，加以详细物理比较和分析：300s 内，在同一图中画出点 、)4,0(、 、 、 （单位：m）温度
12、随时间的变化；)8,0(,6)0,12(6,9(热传导与热辐射大作业报告- 10 -图 1.不同点温度随时间变化曲线图分析：开始时刻通过右、上边界向内部导热，这时候尽管有内热源，但谁相对离右、上边界越近，温度曲线越陡。即开始时刻（0,8）点比（0,4）点温度曲线陡， （12,0）点比（6,0）点温度曲线陡，一定时间后由于有内热源，内部温度逐渐高于边界温度，这时内部开始向边界导热。这时谁离两个绝热边交点越近，谁的温度会越高，这就是为什么最后（0,4）点比（0,8）点温度高，（6,0）点温度比（12,0）点温度高。（b）.200s 内，画出点 、)4,18(、 、 、 （单位：m）处，分别沿 x、
13、y 方向热)8,1(2,6)1,(6,9流密度值随时间的变化；图 2.200s 内 x 方向不同点的热流密度曲线（解析解） 图 3.200s 内 y 方向不同点的热流密度曲线（解析解）图 4.200s 内 x 方向不同点的热流密度曲线（数值解） 图 5.200s 内 y 方向不同点的热流密度曲线（数值解） 热传导与热辐射大作业报告- 11 -图 6.不同点 x 方向热流密度数值解与解析解相对误差分析：右边界（18,4）和（18,8）这两点开始时 X 方向两侧温差较大，故热流密度也会特别大，开始时内部温度较边界温度低，向内部导热，热流密度为负值。后来内部温度大于边界温度，向外散热，热流密度为正值
14、。而上边界点温度相同，故在 X 方向不存在热传导，故导热系数为零。而中间点开始时从右向左导热，热流密度为负，随着边界层温度影响的深入，热流密度绝对值越来越大，但由于有内热源，会使此影响逐渐减弱，故在一段时间后待热流密度达到一个顶峰以后会逐渐减小，后来由于内热源的作用，导热由内向外进行，热流密度也由负值变为正值。Y 方向分析类似。由于（9,6）离上边界更近，故沿 Y 方向达到的下边界峰值更大。（c）.画出 时刻区域内的等温线；sst 1502750、图 8.50s 时区域内的等温线 图 9.75s 时区域内的等温线图 10.100s 时区域内的等温线 图 11.125s 时区域内的等温线图 7.
15、不同点 x 方向热流密度数值解与解析解相对误差热传导与热辐射大作业报告- 12 -图 12.150s 时区域内的等温线分析：开始时刻，尽管有内热源的存在，但边界温度比内部温度高，此时边界向内部传热，故开始时靠近边界的温度比内部高，这就是为什么50、75、100s 时等温线呈现由坐下到右上温度逐渐升高。过一段时间后，中间部分由于内热源和边界热传导的共同作用，而坐下边界此时收到的内热源和边界热传导的作用小于中间部分，故造成了中间部分温度反而比其他部分高。一段时间后，内热源起主导作用，向外散热，这事等温线上的温度由左下到右上逐渐降低。（d）.300s 内，在同一图中画出点 （单位：m）在 分0,90
16、g别等于 ， ， 情况下的温度变化；31mW323mW图 13.不同内热源下温度变化曲线（解析解） 图 14.不同内热源下温度变化曲线（数值解）图 15.不同内热源下数值解与解析解相对误差分析：内热源越大，单位时间内内部产生的能量越多，节点温度升高的越快。在其它条件相同的情况下，内热源越大，最终内部温度也越高。开始时，由于温度变化剧烈，此时解析解和数值解的误差也相对较大，一段时间以后温度趋于稳定，这个时候相对误差也趋于一个较小的稳定值。热传导与热辐射大作业报告- 13 -（e）.300s 内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为 、 和 的温度、KW5.00.1KW5
17、.1热流密度变化；图 16.不同导热数下温度变化曲线（解析解） 图 17.不同导热数下温度变化曲线（数值解）图 18.不同导热系数下数值解和解析解的相对误差热传导与热辐射大作业报告- 14 -图 20.不同导热系数下 X 方向热流密数值解与解析解相对误差图 23.(9,6)点不同导热系数下 y 方向数值解和解析解的相对误差分析：导热系数 K 越大，内部温度越能快速的传递给外界，这就是问什么导热系数越大，节点最终温度低。根据热流密度方程 ，可知。K 越大，xtkq热流密度越大，这就是为什么 K 越大，热流密度最低点峰值越大。而最后由于内热源相同，根据能量守恒，最后导热系数也必然趋近于一个定值。开
18、始时由于温度变化剧烈，在不同的导热系数下同一点温度随之间变化值得数值解和解图 19.不同导热系数下 X 方向热流密度曲线（解析解） 图 20.不同导热系数下 X 方向热流密度曲线（数值解）图 21.不同导热系数下 Y 方向热流密度曲线（解析解） 图 22。不同导热系数下 Y 方向热流密度曲线（数值解）热传导与热辐射大作业报告- 15 -析解的相对误差较大，一段时间后温度趋于稳定，此时数值解和解析解的相对温差是一个较小值。（f）.300s 内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为 、 和 的温度、s24.0s28.s2.1热流密度变化；图 24（9,6）点在不同热扩
19、散系数下的温度曲线 图 25.（9,6）点在不同热扩散系数下 X 方向的热流密度图 26.（9,6）点在不同热扩散系数下 Y 方向的热流密度分析：热扩散系数越大，边界对温度越能快速的影响到内部，这就是为什么同一点热扩散系数越大，温度升高的越快。热扩散系数越大，边界对温度越能快速的影响到内部，这导致最低点峰值向左移动。热扩散系数表示“温度扯平能力” ，热扩散系数越高表示其温度扯平能力越大。如果时间趋于无穷大，最终即使热扩散系数不同，最终温度也会趋于同一个值。300s 对于热扩散系数为0.8 和 1.2 值来说已经时间足够趋于同一个稳定值，但对于 0.4 的值来说时间却不是够大，这就是为什么 30
20、0s 时，热扩散系数为 0.8 和 1.2 的趋于同一值，而 0.4 的却比它们的小。三、关于绘图命令的说明热传导与热辐射大作业报告- 16 -绘图命令大致类似，故我们这里只以 X 方向热流密度为例来说明，其它的绘图命令不再赘述。plot(KX1)hold onplot(KX2,r)plot(KX3,k)plot(KX4,y)plot(KX5,g)xlabel(时间 t)ylabel(x方向热流密度)title(不同点 x方向热流密度曲 线（数值解）)legend(（18,4）, （18,8） ,（6,12 ）,(12,12),(9,6) )热传导与热辐射大作业报告- 17 -个人感想经过一个
21、多星期的连续奋战，终于搞定了这“万恶”的热传导与热辐射的大作业。首先真诚的感谢在作业中帮助过我的老师和同学。本来以为求温度场并不会是一件特别难的事情，可是等到实践时却发现里面有很多自己意想不到的困难。自己的 MATLAB 零基础确实也增添了不少困难。好不容易把程序编出来了，带入运行却是出问题了，总是比时间值少很多，花了一晚上一点一点的查却没有任何结果。知道第二天早上才发现是自己在循环中占用了原先定义的一个量。让人崩溃又让人欣喜：悲的是半天没有结果，喜的是终于找到了问题的根源。这样的事情还有很多很多。有时候为了查一个错误总需要花很长时间，但是经过奋战后终于把问题弄明白的那种欣喜确实很快乐的。在数
22、值解的过程中，出现了一些令人感觉崩溃的问题。比如，步长取大了难以保证精度，取小了计算特别慢，而且出现一个让人再也做不下去的感觉“out of memory”。曾经一次计算了十几个小时最后得出了一个这样的结果，最后只能两者中和取，得出最终结果。从 MATLAB 的零基础、从对温度场求解的模糊认识。这种现象伴随着作业的深入，使自己对这些问题有了一个更加清晰地认识。同时也对 MATLAB这个软件有了一定的了解。最后再次感谢在这次作业中帮助过我的各位同学和老师！热传导与热辐射大作业报告- 18 -附件.计算中所用程序附件 1.解析解完整程序clear all; %清除系统中原有的变量clc; %清除屏
23、幕a=18; %x方向 长度b=12; %y方向 长度g=1; %g为内热源k=1; %k为导热系数ar=0.8; %ar为热扩散系数T0=200; % T0为初始温度T1=600; %T1为边界温度for p=1:19x=p-1;for q=1:13y=q-1;wtj=0;for i=1:15btm=(2*i-1)*pi/(2*a);wtj=wtj+2*g*sin(btm*a)*cosh(btm*y)*cos(btm*x)/(a*k*btm3*cosh(btm*b);endwt=(a2-x2)*g/(2*k)+T1-wtjfor i=1:300fwt=0;for j=1:15for k0=1
24、:15btn=(2*j-1)*pi/(2*a);gmv=(2*k0-1)*pi/(2*b);fwt=fwt+4/(a*b)*sin(btn*a)*sin(gmv*b)/(btn*gmv)*(T0-T1-g/(k*btn2)+g*gmv2/(k*(gmv2+btn2)*btn2)*cos(btn*x)*cos(gmv*y)*exp(-ar*(btn2+gmv2)*t);endendA(p,q,1)=wt+fwt;endend热传导与热辐射大作业报告- 19 -end附件 2.数值解完整程序（显式法）clear alldx=0.3; %dx为x方向步长dy=0.3; %dy为y方向步长dt=0.0
25、1; %dt为时间步长k=1; %k为导热系数sjz=300; %sjs为时间值x=18; %x为x 方向总长度y=12; %y为y 方向总长度ar=0.8; %ar为热扩散率q=1; %q为热流密度sjs=sjz/dt; %sjs为时间节 点数mx=x/dx+1;%mx为x 方向节 点数ny=y/dy+1;%ny为n 方向节 点数目T0=200; %T0为初始温度T1=600; %T1为边界温度Fo=ar*dt/(dx2);T=zeros(mx,ny,sjs/xhs+1);%定义初始温度和边界温度T(mx,:,:)=T1;T(:,ny,:)=T1;T(:,:,1)=T0;xhs=6;for
26、xh=1:xhsfor t=1:sjs/xhsT(1,1,t+1)=2*Fo*(T(2,1,t)+T(1,2,t)+(1-4*Fo)*T(1,1,t)+ar*q*dt/k; %(0,0)处温度计算式for m=2:mx-1T(m,1,t+1)=Fo*(2*T(m,2,t)+T(m-1,1,t)+T(m+1,1,t)+(1-4*Fo)*T(m,1,t)+ar*q*dt/k;%下边界处温度计算式endfor n=2:ny-1T(1,n,t+1)=Fo*(2*T(2,n,t)+T(1,n-1,t)+T(1,n+1,t)+(1-4*Fo)*T(1,n,t)+ar*q*dt/k;%左边界处温度计算式en
27、dfor m=2:mx-1for n=2:ny-1T(m,n,t+1)=Fo*(T(m+1,n,t)+T(m,n+1,t)+T(m-1,n,t)+T(m,n-热传导与热辐射大作业报告- 20 -1,t)+(1-4*Fo)*T(m,n,t)+ar*q*dt/k;%内部温度计算式 endendendfor i=1:sjz/xhsTZ(i+50*(xh-1),1)=T(31,21,(i-1)/dt+1);endfor m=1:mxfor n=1:nyT(m,n,1)=T(m,n,t+1);endendend热传导与热辐射大作业报告- 21 -附件 3.X 和 Y 方向的热流密度计算程序（解析解）X
28、方向function qx=qx(x,y,t)a=18; %x方向 长度b=12; %y方向 长度g=1; %g为内热源k=1; %k为导热系数ar=0.8; %ar为热扩散系数T0=200; % T0为初始温度T1=600; %T1为边界温度wtj=0;for i=1:15btm=(2*i-1)*pi/(2*a);wtj=wtj+2*g*sin(btm*a)*cosh(btm*y)*sin(btm*x)/(a*btm2*cosh(btm*b);endwt=g*x-wtj;fwt=0;for j=1:15for k0=1:15btn=(2*j-1)*pi/(2*a);gmv=(2*k0-1)*
29、pi/(2*b);fwt=fwt+4*k/(a*b)*sin(btn*a)*sin(gmv*b)/gmv*(T0-T1-g/(k*btn2)+g*gmv2/(k*(gmv2+btn2)*btn2)*sin(btn*x)*cos(gmv*y)*exp(-ar*(btn2+gmv2)*t);endendqx=wt+fwt;endY 方向function qy=qy(x,y,t)a=18; %x方向 长度b=12; %y方向 长度g=1; %g为内热源k=1; %k为导热系数ar=0.8; %ar为热扩散系数T0=200; % T0为初始温度T1=600; %T1为边界温度wtj=0;热传导与热辐射
30、大作业报告- 22 -for i=1:15btm=(2*i-1)*pi/(2*a);wtj=wtj+2*g*k*sin(btm*a)*sinh(btm*y)*cos(btm*x)/(a*k*btm2*cosh(btm*b);endfwt=0;for j=1:15for k0=1:15btn=(2*j-1)*pi/(2*a);gmv=(2*k0-1)*pi/(2*b);fwt=fwt+4*k*gmv/(a*b)*sin(btn*a)*sin(gmv*b)/(btn*gmv)*(T0-T1-g/(k*btn2)+g*gmv2/(k*(gmv2+btn2)*btn2)*cos(btn*x)*sin(
31、gmv*y)*exp(-ar*(btn2+gmv2)*t);endendqy=wtj+fwt;end热传导与热辐射大作业报告- 23 -附件 4.X 和 Y 方向的热流密度计算程序（数值解）X 方向dx=0.3; %dx为x方向步长dy=0.3; %dy为y方向步长dt=0.01; %dt为时间步长k=1; %k为导热系数sjz=300; %sjs为时间值x=18; %x为x 方向总长度y=12; %y为y 方向总长度ar=0.8; %ar为热扩散率q=1; %q为热流密度sjs=sjz/dt; %sjs为时间节 点数mx=x/dx+1;%mx为x 方向节 点数ny=y/dy+1;%ny为n
32、方向节 点数目T0=200; %T0为初始温度T1=600; %T1为边界温度Fo=ar*dt/(dx2);xhs=6;T=zeros(mx,ny,sjs/xhs+1);%定义初始温度和边界温度T(mx,:,:)=T1;T(:,ny,:)=T1;T(:,:,1)=T0;for xh=1:xhsfor t=1:sjs/xhsT(1,1,t+1)=2*Fo*(T(2,1,t)+T(1,2,t)+(1-4*Fo)*T(1,1,t)+ar*q*dt/k; %(0,0)处温度计算式for m=2:mx-1T(m,1,t+1)=Fo*(2*T(m,2,t)+T(m-1,1,t)+T(m+1,1,t)+(1
33、-4*Fo)*T(m,1,t)+ar*q*dt/k;%下边界处温度计算式endfor n=2:ny-1T(1,n,t+1)=Fo*(2*T(2,n,t)+T(1,n-1,t)+T(1,n+1,t)+(1-4*Fo)*T(1,n,t)+ar*q*dt/k;%左边界处温度计算式endfor m=2:mx-1for n=2:ny-1T(m,n,t+1)=Fo*(T(m+1,n,t)+T(m,n+1,t)+T(m-1,n,t)+T(m,n-1,t)+(1-4*Fo)*T(m,n,t)+ar*q*dt/k;%内部温度计算式 endend热传导与热辐射大作业报告- 24 -endfor i=1:sjz/x
34、hsKX1(i+50*(xh-1),1)=k*(T(60,14,(i-1)/dt+1)-T(61,14,(i-1)/dt+1)/dx; %（18,4）点x 方向热流密度KX2(i+50*(xh-1),1)=k*(T(60,28,(i-1)/dt+1)-T(61,28,(i-1)/dt+1)/dx; %（18,8）点x方向热流密度KX3(i+50*(xh-1),1)=k*(T(20,41,(i-1)/dt+1)-T(22,41,(i-1)/dt+1)/(2*dx); %（6,12）点x 方向热流密度KX4(i+50*(xh-1),1)=k*(T(40,41,(i-1)/dt+1)-T(42,41
35、,(i-1)/dt+1)/(2*dx); %（12,12）点x方向热流密度KX5(i+50*(xh-1),1)=k*(T(30,21,(i-1)/dt+1)-T(32,21,(i-1)/dt+1)/(2*dx); %（9,6）点x方向热流密度endfor m=1:mxfor n=1:nyT(m,n,1)=T(m,n,t+1);endendendY 方向dx=0.3; %dx为x方向步长dy=0.3; %dy为y方向步长dt=0.01; %dt为时间步长k=1; %k为导热系数sjz=300; %sjs为时间值x=18; %x为x 方向总长度y=12; %y为y 方向总长度ar=0.8; %ar
36、为热扩散率q=1; %q为热流密度sjs=sjz/dt; %sjs为时间节 点数mx=x/dx+1;%mx为x 方向节 点数ny=y/dy+1;%ny为n 方向节 点数目T0=200; %T0为初始温度T1=600; %T1为边界温度Fo=ar*dt/(dx2);xhs=6;T=zeros(mx,ny,sjs/xhs+1);%定义初始温度和边界温度T(mx,:,:)=T1;T(:,ny,:)=T1;T(:,:,1)=T0;热传导与热辐射大作业报告- 25 -for xh=1:xhsfor t=1:sjs/xhsT(1,1,t+1)=2*Fo*(T(2,1,t)+T(1,2,t)+(1-4*Fo
37、)*T(1,1,t)+ar*q*dt/k; %(0,0)处温度计算式for m=2:mx-1T(m,1,t+1)=Fo*(2*T(m,2,t)+T(m-1,1,t)+T(m+1,1,t)+(1-4*Fo)*T(m,1,t)+ar*q*dt/k;%下边界处温度计算式endfor n=2:ny-1T(1,n,t+1)=Fo*(2*T(2,n,t)+T(1,n-1,t)+T(1,n+1,t)+(1-4*Fo)*T(1,n,t)+ar*q*dt/k;%左边界处温度计算式endfor m=2:mx-1for n=2:ny-1T(m,n,t+1)=Fo*(T(m+1,n,t)+T(m,n+1,t)+T(m
38、-1,n,t)+T(m,n-1,t)+(1-4*Fo)*T(m,n,t)+ar*q*dt/k;%内部温度计算式 endendendfor i=1:sjz/xhsKY1(i+50*(xh-1),1)=k*(T(61,13,(i-1)/dt+1)-T(61,15,(i-1)/dt+1)/(2*dy); %（18,4)点y方向热流密度KY2(i+50*(xh-1),1)=k*(T(61,27,(i-1)/dt+1)-T(61,29,(i-1)/dt+1)/(2*dy); %（18,8)点y方向热流密度KY3(i+50*(xh-1),1)=k*(T(21,40,(i-1)/dt+1)-T(21,41,(i-1)/dt+1)/dy; %（6,12)点y 方向热流密度KY4(i+50*(xh-1),1)=k*(T(41,40,(i-1)/dt+1)-T(41,41,(i-1)/dt+1)/dy; %（12,12) 点y方向热流密度KY5(i+50*(xh-1),1)=k*(T(31,20,(i-1)/dt+1)-T(31,22,(i-1)/dt+1)/(2*dy); %（9,6)点y 方向热流密度endfor m=1:mxfor n=1:nyT(m,n,1)=T(m,n,t+1);endendend
