热传导及热辐射大作业报告

热传导与热辐射大作业报告

专业整理分享

完美WORD格式

目录

一、作业题目 ............................................................. - 1 - 二、作业解答 ............................................................. - 2 - 个人感想................................................................ - 17 - 附件.计算中所用程序 ..................................................... - 18 -

专业整理分享

完美WORD格式

一、作业题目

一矩形平板0xa, 0yb，内有均匀恒定热源g0，在x0及y0处绝热，在xa及yb处保持温度T1，初始时刻温度为T0，如右图1所示：

1、求t0时，矩形区域内的温度分布Tx,y,t的解析表达式；

2、若a18m，b12m，g01Wm3，T1600K，T0200K，热传导系数 k1.0WmK，热扩散系数0.8m2s。请根据1中所求温度分布用

MATLAB软件绘出下列结果，加以详细物理比较和分析：

(a) 300s内，在同一图中画出点(0,4)、(0,8)、6,0、(12,0)、(9,6)（单位：m）温度随时间的变化；

(b) 200s内，画出点(18,4)、(18,8)、6,12、(12,12)、(9,6)（单位：m）处，分别沿x、y方向热流密度值随时间的变化；

(c) 画出t50s、75s、100s、125s、150s时刻区域内的等温线；

(d) 300s内，在同一图中画出点9,0（单位：m）在g0分别等于1Wm3，

2Wm3，3Wm3情况下的温度变化；

(e) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为0.5WmK、1.0WmK和1.5WmK的温度、热流密度变化；

(f) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为0.4m2s、0.8m2s和1.2m2s的温度、热流密度变化； 3、运用有限差分法计算2中(b)、(d)和(e)，并与解析解结果进行比较，且

需将数值解与解析解的相对误差减小到1‰以下；

4、附上源程序和个人体会；

以报告形式整理上述结果，用A4纸打印上交。

专业整理分享

完美WORD格式

二、作业解答

1、求t0时，矩形区域内的温度分布Tx,y,t的解析表达式；

解答：我们令θTT1，则可以得到一个方程和边界条件：

22g01T 22xykt （1-1）

ddx 0，xa

0,x0.

d0,y0 dy 0，yb T0T1，t0

（ 将上式分解为一个sx，y)的稳态问题：

2s2sg0 220 （1-2）

xyk dsdx0,x0.

s0，xa

ds0,y0 dy s0，yb 和一个（hx,y,t)的其次问题：

2h2h1T 22xyt （1-3）

dh0,x0. dx h0，xa

dh0,y0 dy 专业整理分享

完美WORD格式

h0，yb

（ 其中（hx,y)F(x,y)sx,y)f(x,y)

则原问题的解根据下式求得：

（x,y,t)s(x,y)h(x,y,t) （1-4）

（发热强度为常数的特解可从表2-4中查的，则新变量sx，y)可定义为：

g02 （x，y)（x，y)xA （1-5） s2k将（1-5）带入（1-2）整理得到：

2（x,y)2(x,y)0， 0xa,0yb （1-6） 22xy ddx0,x0. g02aA，xa 2k 

d0,y0 dyg02xA，yb 2kg若令常数A0a2，则上式可以变为：

2k 2（x,y)2(x,y)0， 0xa,0yb （1-7） x2y2 d0,x0. dx 0，xa

d0,y0 dyf(x)，yb

其中f(x)g02(xa2) 2k假定（x,y)可以分离出如下形式：

专业整理分享

完美WORD格式

（x,y)X(x)Y(y) （1-8） 对应于X(x)和Y(y)的分离方程为：

d2X(x)2X(x)0 （1-9） 2dx dX

0,x0 dx X0，xa

d2Y(y) 2Y(y)0 （1-10） 2dy

dY0,y0 dy Yf(x)，yb

在X(x)中特征值问题的解可以直接从表2-2第6条中得到，只需要用a代替L，

X（m,x)cosmx （1-11）

12 （1-12） N（m)am是下面方程的正根：

cosma0 （1-13） 方程（1-10）的解可以取为

coshmy （1-14） Y（m,y）（x,y)的完全解由下式组成：

(x,y)Cmcoshmycosmx （1-15）

m1

专业整理分享

完美WORD格式

此式满足热传导问题（1-7）及三个齐次边界条件，其中，系数Cm可以根据方程的解还应满足非齐次的边界条件来决定。利用yb的边界条件可得：

f(x)Cmcoshmbcosmx 0xa （1-16）

m1利用函数cosmx的正交性可以求得系数Cm， Cm式中：

a1'''cosxf(x)dxmN（m）coshmb0 （1-17）

a0cosmxf(x)dxcosmx'0'''18g0'22'g0sinma(xa)dx32kkm

将这个表达式带入式（1-15），其中范数N（m）在前面已经给出，解得结果为

（x,y)2则：

（sx，y)（x，y)g02xA 2k1gsinacoshmycosmx03m （1-18）

kmm1acoshmb 2

1gsinagcoshmycosmx03m0(a2x2)

km2km1acoshmb （1-19） 假定h(x,y,t)分离成如下表达式

h(x,y,t)(t)X(x)Y(y) （1-20） 对应于函数X(x)和Y(y)的分离方程为

专业整理分享

完美WORD格式

d2X(x)2X(x)0 0xa （1-21） X(x)： 2dx dXdx0,x0

X0，xa

d2Y(y)2Y(y)： Y(y)0 0yb （1-22） 2dy

dY0,y0 dy Y0，yb (t)的解为： (t)e22(nv)t （1-23）

上述问题的完全解为： h(x,y,t)Cmne(m1n12n2v)tX(n,x)Y(v,y) （1-24）

其中0h(x,y)CnvX(n,x)Y(v,y) （1-25）

m1n1其中0确定未知系数Cmn的方法是，在上式两边逐项用如下算子作运算：

0X(n,x)dx及0Y(v,y)dy 并利用这些函数的正交性，得到：

ab1CnvX(n,x')Y(v,y')h(x',y')dx'dy' （1-26） N（n）N（v）00ab 专业整理分享

完美WORD格式

最终得到问题的解为：

h(x,y,t)e(m1n122nv)t1X(n,x)Y(v,y)

N（n）N（v）a00bX（n，x'）Y(v,y')h(x',y')dx'dy' （1-27）

式中出现的特征函数，特征值及范数可以从表2-2中直接查得：

X(n，x)cosnx （1-28）

12 （1-29） N(n)a且m为如下方程的正根：

cosna0 （1-30） 满足特征值问题的函数Y(n,y)对应于表2-2中的第6条，得到：

Y(v,y)cosvy （1-31）

12 （1-32） N(v)b且n是如下方程的正根： cosvb0 最后得到：

h(x,y,t)4e(m1n12n2v)t1cosnxcosvy ab 00cosnx'cosvy'h(x',y')dx'dy' （1-33） 令I00cosnx'cosvy'h(x',y')dx'dy'

a00ababbcosnxcosvy[T0T12''g0sinmag02coshmy'cosmx'(ax2)]dx'dy'3m2km1akcoshmb 专业整理分享

完美WORD格式

a00b'cosnx'cosvy（T0T1）dx'dy'

0a0bcosnx'cosvy'22kg0sinma'' coshmy'cosmx'dxdy3mm1akcoshmb 0a0bcosnx'cosvy'g0(a2x2)dx'dy'

'（T0T1）dx'dy'（T0T1）其中令I100cosnx'cosvyabsinnasinvbnv

令I20a0bcosnx'cosvy'2 g0sinma'' coshmy'cosmx'dxdy3akcoshbm1mm2g0sinma3m1akcoshmbm0ab0cosnx'cosvy'coshmy'cosmx'dx'dy'

b2g0sinma 3m1akcoshmbma0cosnxcosmxdxcosvy'coshmy'dy'

'''0根据余弦函数的正交性，只有当m=n时积分才不为0，故上式可以化为：

b2g0sinmaa'''cosxcosxdxcosvy'coshny'dy' nn30akcoshnbn0再令I210cosnx'cosnx'dx'a

a2 I220cosvy'coshny'dy'bvn2v2coshnbsinvb

2g0(1)nVsinmasinvb所以I2 II21223akcoshnbnn3(n2v2)na 令I30a0bcosnx'cosvy'g0(a2x2)dx'dy'g0sinvbsin32kkvng0g0v2所以II1I2I3 [(T0T12)222nvknkn(nv)1由（1-4）、（1-19）及（1-33）可知

g0sinmaT(x,y,t)2akcoshmycosmx3coshbm1mmg02(ax2) 2k 专业整理分享

完美WORD格式

4(n2v2)tsinmasinvbg0g0v2e[(T0T12)22]cosnxcosvyT12abm1v1nvknkn(nv)。

以上是解析解的全过程，具体值的计算采用MATLAB编程计算求取。

2、若a18m，b12m，g01Wm3，T1600K，T0200K，热传导系数k1.0WmK，热扩散系数0.8m2s。请根据

1中所求温

度分布用MATLAB软件绘出下列结果，加以详细物理比较和分析：

300s内，在同一图中画出点(0,4)、(0,8)、6,0、(12,0)、(9,6)（单位：m）温度随时间的变化；

图1.不同点温度随时间变化曲线图

分析：开始时刻通过右、上边界向内部导热，这时候尽管有内热源，但谁相对离右、上边界越近，温度曲线越陡。即开始时刻（0,8）点比（0,4）点温度曲线陡，（12,0）点比（6,0）点温度曲线陡，一定时间后由于有内热源，内部温度逐渐高于边界温度，这时内部开始向边界导热。这时谁离两个绝热边交点越近，谁的温度会越高，这就是为什么最后（0,4）点比（0,8）点温度高，（6,0）点温度比（12,0）点温度高。

（b）.200s内，画出点(18,4)、(18,8)、6,12、(12,12)、(9,6)（单位：m）处，分别沿x、y方向热流密度值随时间的变化；

专业整理分享

完美WORD格式

图2.200s内x方向不同点的热流密度曲线（解析解） 图3.200s内y方向不同点的热流密度曲线（解析

解）

图4.200s内x方向不同点的热流密度曲线（数值解） 图5.200s内y方向不同点的热流密度曲线（数值解）

图6.不同点x方向热流密度数值解与解析解相对误差

分析：右边界（18,4）和（18,8）这两点开始时X方向两侧温差较大，故热流密度也会特别大，开始时内部温度较边界温度低，向内部导热，热流密度为负值。后来内部温度大于边界温度，向外散热，热流密度为正值。而上边界点温度相同，故在X方向不存在热传导，故导热系数为零。而中间点开始时从右向左导热，热流密度为负，随着边界层温度影响的深入，热流密度绝对值越来越大，但由于有内热源，会使此影响逐渐减弱，故在一段时间后待热流密度达到一个顶峰

专业整理分享

7.不同点x方向热

图

完美WORD格式

以后会逐渐减小，后来由于内热源的作用，导热由内向外进行，热流密度也由负值变为正值。Y方向分析类似。由于（9,6）离上边界更近，故沿Y方向达到的下边界峰值更大。

75s、100s、125s、150s时刻区域内的等温线； （c）.画出t50s、

图8.50s时区域内的等温线 图9.75s时区域内的等温线

图10.100s时区域内的等温线 图11.125s时区域内的等温线

图12.150s时区域内的等温线

分析：开始时刻，尽管有内热源的存在，但边界温度比内部温度高，此时边

界向内部传热，故开始时靠近边界的温度比内部高，这就是为什么50、75、100s时等温线呈现由坐下到右上温度逐渐升高。过一段时间后，中间部分由于内热源和边界热传导的共同作用，而坐下边界此时收到的内热源和边界热传导的作用小于中间部分，故造成了中间部分温度反而比其他部分高。一段时间后，内热源起主导作用，向外散热，这事等温线上的温度由左下到右上逐渐降低。

（d）.300s内，在同一图中画出点9,0（单位：m）在g0分

专业整理分享

完美WORD格式

3331Wm2Wm3Wm别等于，，情

况下的温度变化；

图13.不同内热源下温度变化曲线（解析解） 图14.不同内热源下温度变化曲线（数值解）

图15.不同内热源下数值解与解析解相对误差

分析：内热源越大，单位时间内内部产生的能量越多，节点温度升高的越快。在其它条件相同的情况下，内热源越大，最终内部温度也越高。开始时，由于温度变化剧烈，此时解析解和数值解的误差也相对较大，一段时间以后温度趋于稳定，这个时候相对误差也趋于一个较小的稳定值。

（e）.300s内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为0.5W流密度变化；

mK、1.0WmK和1.5WmK的温度、热

专业整理分享

完美WORD格式

图16.不同导热数下温度变化曲线（解析解） 图17.不同导热数下温度变化曲线（数值解）

图18.不同导热系数下数值解和解析解的相对误差

（数值解） 图19.不同导热系数下X方向热流密度曲线（解析解） 图20.不同导热系数下X方向热流密度曲线

专业整理分享

完美WORD格式

图20.不同导热系数下X方向热流密数值解与解析解相对误差

图21. 不同导热系数下Y方向热流密度曲线（解析解） 图22。不同导热系数下Y方向热流密度曲线（数值解）

图23.(9,6)点不同导热系数下y方向数值解和解析解的相对误差

分析：导热系数K越大，内部温度越能快速的传递给外界，这就是问什么导热系数越大，节点最终温度低。根据热流密度方程qkt，可知。K越大，热

x流密度越大，这就是为什么K越大，热流密度最低点峰值越大。而最后由于内热源相同，根据能量守恒，最后导热系数也必然趋近于一个定值。开始时由于温度变化剧烈，在不同的导热系数下同一点温度随之间变化值得数值解和解析解的相对误差较大，一段时间后温度趋于稳定，此时数值解和解析解的相对温差是一个较小值。

专业整理分享

完美WORD格式

（f）.300s内，比较点(9,6) （单位：m）在其它参数不变

2220.4ms0.8ms1.2ms的温度、热流情况下热扩散系数分别为、和

密度变化；

图24（9,6）点在不同热扩散系数下的温度曲线 图25.（9,6）点在不同热扩散系数下X方向的热流密度

图26.（9,6）点在不同热扩散系数下Y方向的热流密度

分析：热扩散系数越大，边界对温度越能快速的影响到内部，这就是为什么同一点热扩散系数越大，温度升高的越快。热扩散系数越大，边界对温度越能快速的影响到内部，这导致最低点峰值向左移动。热扩散系数表示“温度扯平能力”，热扩散系数越高表示其温度扯平能力越大。如果时间趋于无穷大，最终即使热扩散系数不同，最终温度也会趋于同一个值。300s对于热扩散系数为0.8和1.2值来说已经时间足够趋于同一个稳定值，但对于0.4的值来说时间却不是够大，这就是为什么300s时，热扩散系数为0.8和1.2的趋于同一值，而0.4的却比它们的小。

三、关于绘图命令的说明

绘图命令大致类似，故我们这里只以X方向热流密度为例来说明，其它的绘图命令不再赘述。 plot(KX1)

专业整理分享

完美WORD格式

hold on

plot(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)')

专业整理分享

完美WORD格式

个人感想

经过一个多星期的连续奋战，终于搞定了这“万恶”的热传导与热辐射的大作业。首先真诚的感谢在作业中帮助过我的老师和同学。

本来以为求温度场并不会是一件特别难的事情，可是等到实践时却发现里面有很多自己意想不到的困难。自己的MATLAB零基础确实也增添了不少困难。好不容易把程序编出来了，带入运行却是出问题了，总是比时间值少很多，花了一晚上一点一点的查却没有任何结果。知道第二天早上才发现是自己在循环中占用了原先定义的一个量。让人崩溃又让人欣喜：悲的是半天没有结果，喜的是终于找到了问题的根源。这样的事情还有很多很多。有时候为了查一个错误总需要花很长时间，但是经过奋战后终于把问题弄明白的那种欣喜确实很快乐的。

在数值解的过程中，出现了一些令人感觉崩溃的问题。比如，步长取大了难以保证精度，取小了计算特别慢，而且出现一个让人再也做不下去的感觉“out of memory”。曾经一次计算了十几个小时最后得出了一个这样的结果，最后只能两者中和取，得出最终结果。

从MATLAB的零基础、从对温度场求解的模糊认识。这种现象伴随着作业的深入，使自己对这些问题有了一个更加清晰地认识。同时也对MATLAB这个软件有了一定的了解。

最后再次感谢在这次作业中帮助过我的各位同学和老师！

专业整理分享

完美WORD格式

附件.计算中所用程序

附件1.解析解完整程序

clear all; %清除系统中原有的变量 clc; %清除屏幕 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:19 x=p-1; for q=1:13 y=q-1; wtj=0; for i=1:15

btm=(2*i-1)*pi/(2*a);

wtj=wtj+2*g*sin(btm*a)*cosh(btm*y)*cos(btm*x)/(a*k*btm^3*cosh(btm*b)); end

wt=(a^2-x^2)*g/(2*k)+T1-wtj for i=1:300 fwt=0; for j=1:15

for k0=1:15

btn=(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*btn^2))+g*gmv^2/(k*(gmv^2+btn^2)*btn^2))*cos(btn*x)*cos(gmv*y)*exp(-ar*(btn^2+gmv^2)*t);

end end

A(p,q,1)=wt+fwt; end end end

附件2.数值解完整程序（显式法）

专业整理分享

完美WORD格式

clear all

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方向节点数目 T0=200; %T0为初始温度 T1=600; %T1为边界温度 Fo=ar*dt/(dx^2);

T=zeros(mx,ny,sjs/xhs+1); %定义初始温度和边界温度 T(mx,:,:)=T1; T(:,ny,:)=T1; T(:,:,1)=T0; xhs=6;

for xh=1:xhs for t=1:sjs/xhs

T(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-1

T(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;%下边界处温度计算式 end

for n=2:ny-1

T(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;%左边界处温度计算式 end

for m=2:mx-1 for n=2:ny-1

T(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;%内部温度计算式 end end end

专业整理分享

完美WORD格式

for i=1:sjz/xhs

TZ(i+50*(xh-1),1)=T(31,21,(i-1)/dt+1); end

for m=1:mx for n=1:ny

T(m,n,1)=T(m,n,t+1); end end end

附件3.X和Y方向的热流密度计算程序（解析解）

①X方向

function qx=qx(x,y,t)

专业整理分享

完美WORD格式

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:15

btm=(2*i-1)*pi/(2*a);

wtj=wtj+2*g*sin(btm*a)*cosh(btm*y)*sin(btm*x)/(a*btm^2*cosh(btm*b)); end

wt=g*x-wtj; fwt=0; for j=1:15

for k0=1:15

btn=(2*j-1)*pi/(2*a); gmv=(2*k0-1)*pi/(2*b);

fwt=fwt+4*k/(a*b)*sin(btn*a)*sin(gmv*b)/gmv*((T0-T1-g/(k*btn^2))+g*gmv^2/(k*(gmv^2+btn^2)*btn^2))*sin(btn*x)*cos(gmv*y)*exp(-ar*(btn^2+gmv^2)*t); end end

qx=wt+fwt; end

②Y方向

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; for i=1:15

btm=(2*i-1)*pi/(2*a);

wtj=wtj+2*g*k*sin(btm*a)*sinh(btm*y)*cos(btm*x)/(a*k*btm^2*cosh(btm*b)); end

fwt=0; for j=1:15

for k0=1:15

btn=(2*j-1)*pi/(2*a);

专业整理分享

完美WORD格式

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*btn^2))+g*gmv^2/(k*(gmv^2+btn^2)*btn^2))*cos(btn*x)*sin(gmv*y)*exp(-ar*(btn^2+gmv^2)*t); end end

qy=wtj+fwt; end

附件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方向总长度

专业整理分享

完美WORD格式

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/(dx^2); xhs=6;

T=zeros(mx,ny,sjs/xhs+1); %定义初始温度和边界温度 T(mx,:,:)=T1; T(:,ny,:)=T1; T(:,:,1)=T0;

for xh=1:xhs for t=1:sjs/xhs

T(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-1

T(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;%下边界处温度计算式 end

for n=2:ny-1

T(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;%左边界处温度计算式 end

for m=2:mx-1 for n=2:ny-1

T(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;%内部温度计算式 end end end

for i=1:sjz/xhs

KX1(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); %

专业整理分享

完美WORD格式

（6,12）点x方向热流密度

KX4(i+50*(xh-1),1)=k*(T(40,41,(i-1)/dt+1)-T(42,41,(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方向热流密度 end

for m=1:mx for n=1:ny

T(m,n,1)=T(m,n,t+1); end end end

②Y方向

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方向节点数目 T0=200; %T0为初始温度 T1=600; %T1为边界温度 Fo=ar*dt/(dx^2); xhs=6;

T=zeros(mx,ny,sjs/xhs+1); %定义初始温度和边界温度 T(mx,:,:)=T1; T(:,ny,:)=T1; T(:,:,1)=T0;

for xh=1:xhs for t=1:sjs/xhs

T(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-1

专业整理分享

完美WORD格式

T(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;%下边界处温度计算式 end

for n=2:ny-1

T(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;%左边界处温度计算式 end

for m=2:mx-1 for n=2:ny-1

T(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;%内部温度计算式 end end end

for i=1:sjz/xhs

KY1(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方向热流密度 end

for m=1:mx for n=1:ny

T(m,n,1)=T(m,n,t+1); end end end

专业整理分享