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

1、.wd.wd.wd.热传导与热辐射大作业报告目录TOC o 1-1 h u HYPERLINK l _Toc5749 一、作业题目 PAGEREF _Toc5749 - 1 - HYPERLINK l _Toc18253 二、作业解答 PAGEREF _Toc18253 - 2 - HYPERLINK l _Toc13344 个人感想 PAGEREF _Toc13344 - 17 - HYPERLINK l _Toc11721 附件.计算中所用程序 PAGEREF _Toc11721 - 18 -一、作业题目一矩形平板, ，内有均匀恒定热源，在及处绝热，在及处保持温度，初始时刻温度为，如右图1

2、所示：1、求时，矩形区域内的温度分布的解析表达式；2、假设，热传导系数，热扩散系数。请根据1中所求温度分布用MATLAB软件绘出以下结果，加以详细物理对比和分析：300s内，在同一图中画出点、单位：m温度随时间的变化；200s内，画出点、单位：m处，分别沿x、y方向热流密度值随时间的变化；画出时刻区域内的等温线；300s内，在同一图中画出点单位：m在分别等于，情况下的温度变化；300s内，对比点(9,6) 单位：m在其它参数不变情况下热导率分别为、和的温度、热流密度变化；300s内，对比点(9,6) 单位：m在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化；3、运用有限差分法计算2

3、中(b)、(d)和(e)，并与解析解结果进展对比，且需将数值解与解析解的相对误差减小到1以下；4、附上源程序和个人体会；以报告形式整理上述结果，用A4纸打印上交。二、作业解答1、求时，矩形区域内的温度分布的解析表达式；解答：我们令，那么可以得到一个方程和边界条件： 1-1 将上式分解为一个的稳态问题： 1-2 和一个的其次问题： 1-3 其中那么原问题的解根据下式求得： 1-4发热强度为常数的特解可从表2-4中查的，那么新变量可定义为： 1-5将1-5带入1-2整理得到： 1-6 假设令常数，那么上式可以变为： 1-7 其中假定可以别离出如下形式： 1-8对应于的别离方程为： 1-9 1-10

4、在中特征值问题的解可以直接从表2-2第6条中得到，只需要用a代替L， 1-11 1-12是下面方程的正根： 1-13方程1-10的解可以取为 1-14的完全解由下式组成： 1-15此式满足热传导问题1-7及三个齐次边界条件，其中，系数可以根据方程的解还应满足非齐次的边界条件来决定。利用的边界条件可得： 1-16利用函数的正交性可以求得系数， 1-17式中：将这个表达式带入式1-15，其中范数在前面已经给出，解得结果为 1-18那么： 1-19假定别离成如下表达式 1-20对应于函数和的别离方程为： 1-21： 1-22的解为： 1-23上述问题的完全解为： 1-24其中0 xa,0yb。当t=

5、0时，上式变为： 1-25其中0 xa,0yb。确定未知系数的方法是，在上式两边逐项用如下算子作运算：及并利用这些函数的正交性，得到： 1-26最终得到问题的解为： 1-27式中出现的特征函数，特征值及范数可以从表2-2中直接查得： 1-28 1-29且为如下方程的正根： 1-30满足特征值问题的函数对应于表2-2中的第6条，得到： 1-31 1-32且是如下方程的正根：最后得到： 1-33令其中令令根据余弦函数的正交性，只有当m=n时积分才不为0，故上式可以化为：再令所以令所以由1-4、1-19及1-33可知。以上是解析解的全过程，具体值的计算采用MATLAB编程计算求取。2、假设，热传导系

6、数，热扩散系数。请根据1中所求温度分布用MATLAB软件绘出以下结果，加以详细物理对比和分析：300s内，在同一图中画出点、单位：m温度随时间的变化；图1.不同点温度随时间变化曲线图分析：开场时刻通过右、上边界向内部导热，这时候尽管有内热源，但谁相对离右、上边界越近，温度曲线越陡。即开场时刻0,8点比0,4点温度曲线陡，12,0点比6,0点温度曲线陡，一定时间后由于有内热源，内部温度逐渐高于边界温度，这时内部开场向边界导热。这时谁离两个绝热边交点越近，谁的温度会越高，这就是为什么最后0,4点比0,8点温度高，6,0点温度比12,0点温度高。b.200s内，画出点、单位：m处，分别沿x、y方向热

7、流密度值随时间的变化；图2.200s内x方向不同点的热流密度曲线解析解 图3.200s内y方向不同点的热流密度曲线解析解图4.200s内x方向不同点的热流密度曲线数值解 图5.200s内y方向不同点的热流密度曲线数值解 图7.不同点x方向热流密度数值解与解析解相对误差 图6.不同点x方向热流密度数值解与解析解相对误差分析：右边界18,4和18,8这两点开场时X方向两侧温差较大，故热流密度也会特别大，开场时内部温度较边界温度低，向内部导热，热流密度为负值。后来内部温度大于边界温度，向外散热，热流密度为正值。而上边界点温度一样，故在X方向不存在热传导，故导热系数为零。而中间点开场时从右向左导热，热

8、流密度为负，随着边界层温度影响的深入，热流密度绝对值越来越大，但由于有内热源，会使此影响逐渐减弱，故在一段时间后待热流密度到达一个顶峰以后会逐渐减小，后来由于内热源的作用，导热由内向外进展，热流密度也由负值变为正值。Y方向分析类似。由于9,6离上边界更近，故沿Y方向到达的下边界峰值更大。c.画出时刻区域内的等温线； 图8.50s时区域内的等温线 图9.75s时区域内的等温线 图10.100s时区域内的等温线 图11.125s时区域内的等温线图12.150s时区域内的等温线分析：开场时刻，尽管有内热源的存在，但边界温度比内部温度高，此时边界向内部传热，故开场时靠近边界的温度比内部高，这就是为什么

9、50、75、100s时等温线呈现由坐下到右上温度逐渐升高。过一段时间后，中间局部由于内热源和边界热传导的共同作用，而坐下边界此时收到的内热源和边界热传导的作用小于中间局部，故造成了中间局部温度反而比其他局部高。一段时间后，内热源起主导作用，向外散热，这事等温线上的温度由左下到右上逐渐降低。d.300s内，在同一图中画出点单位：m在分别等于，情况下的温度变化；图13.不同内热源下温度变化曲线解析解 图14.不同内热源下温度变化曲线数值解图15.不同内热源下数值解与解析解相对误差分析：内热源越大，单位时间内内部产生的能量越多，节点温度升高的越快。在其它条件一样的情况下，内热源越大，最终内部温度也越

10、高。开场时，由于温度变化剧烈，此时解析解和数值解的误差也相对较大，一段时间以后温度趋于稳定，这个时候相对误差也趋于一个较小的稳定值。e.300s内，对比点(9,6) 单位：m在其它参数不变情况下热导率分别为、和的温度、热流密度变化；图16.不同导热数下温度变化曲线解析解 图17.不同导热数下温度变化曲线数值解图18.不同导热系数下数值解和解析解的相对误差图20.不同导热系数下X方向热流密度曲线数值解图19.不同导热系数下X方向热流密度曲线解析解 图20.不同导热系数下X方向热流密数值解与解析解相对误差图22。不同导热系数下Y方向热流密度曲线数值解图21.不同导热系数下Y方向热流密度曲线解析解

11、图23.(9,6)点不同导热系数下y方向数值解和解析解的相对误差分析：导热系数K越大，内部温度越能快速的传递给外界，这就是问什么导热系数越大，节点最终温度低。根据热流密度方程，可知。K越大，热流密度越大，这就是为什么K越大，热流密度最低点峰值越大。而最后由于内热源一样，根据能量守恒，最后导热系数也必然趋近于一个定值。开场时由于温度变化剧烈，在不同的导热系数下同一点温度随之间变化值得数值解和解析解的相对误差较大，一段时间后温度趋于稳定，此时数值解和解析解的相对温差是一个较小值。f.300s内，对比点(9,6) 单位：m在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化； 图249,6点

12、在不同热扩散系数下的温度曲线 图25.9,6点在不同热扩散系数下X方向的热流密度图26.9,6点在不同热扩散系数下Y方向的热流密度分析：热扩散系数越大，边界对温度越能快速的影响到内部，这就是为什么同一点热扩散系数越大，温度升高的越快。热扩散系数越大，边界对温度越能快速的影响到内部，这导致最低点峰值向左移动。热扩散系数表示“温度扯平能力，热扩散系数越高表示其温度扯平能力越大。如果时间趋于无穷大，最终即使热扩散系数不同，最终温度也会趋于同一个值。300s对于热扩散系数为0.8和1.2值来说已经时间足够趋于同一个稳定值，但对于0.4的值来说时间却不是够大，这就是为什么300s时，热扩散系数为0.8和

13、1.2的趋于同一值，而0.4的却比它们的小。三、关于绘图命令的说明 绘图命令大致类似，故我们这里只以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)个人感想经过一个多星期的连续奋战，终于搞定了这“万恶的热传导与热辐射的大作业。首先真诚的感谢在作业中帮助过我的教师和同学。本来以为求温度场并不会是一件特别难的事情，

14、可是等到实践时却发现里面有很多自己意想不到的困难。自己的MATLAB零根基确实也增添了不少困难。好不容易把程序编出来了，带入运行却是出问题了，总是比时间值少很多，花了一晚上一点一点的查却没有任何结果。知道第二天早上才发现是自己在循环中占用了原先定义的一个量。让人崩溃又让人欣喜：悲的是半天没有结果，喜的是终于找到了问题的根源。这样的事情还有很多很多。有时候为了查一个错误总需要花很长时间，但是经过奋战后终于把问题弄明白的那种欣喜确实很快乐的。在数值解的过程中，出现了一些令人感觉崩溃的问题。比方，步长取大了难以保证精度，取小了计算特别慢，而且出现一个让人再也做不下去的感觉“out of memory

15、。曾经一次计算了十几个小时最后得出了一个这样的结果，最后只能两者中和取，得出最终结果。从MATLAB的零根基、从对温度场求解的模糊认识。这种现象伴随着作业的深入，使自己对这些问题有了一个更加清晰地认识。同时也对MATLAB这个软件有了一定的了解。最后再次感谢在这次作业中帮助过我的各位同学和教师！附件.计算中所用程序附件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为边界温度fo

16、r 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*btm3*cosh(btm*b);endwt=(a2-x2)*g/(2*k)+T1-wtjfor i=1:300 fwt=0;for j=1:15for 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)

17、*(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;endendend附件2.数值解完整程序显式法clear alldx=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

18、为时间节点数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 xh=1:xhsfor 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+

19、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-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;%左边界处温度计算式endfor m=2:mx-1for 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;%内部温度

20、计算式 endendendfor i=1:sjz/xhs TZ(i+50*(xh-1),1)=T(31,21,(i-1)/dt+1);endfor m=1:mxfor n=1:ny T(m,n,1)=T(m,n,t+1);endendend附件3.X和Y方向的热流密度计算程序解析解X方向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:15 btm=(2*i-1)*pi/(2*a)

21、; 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: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*btn2)+g*gmv2/(k*(gmv2+btn2)*btn2)*sin(btn*x)*cos(gmv*y)*exp(-ar*(btn2+gmv2)*t);endendqx=

22、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;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*btm2*cosh(btm*b);end fwt=0;for j=1:15for k0=1:15 btn=(2*j-1)*pi/(2*a); gmv=(

23、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(gmv*y)*exp(-ar*(btn2+gmv2)*t);endendqy=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方向总

24、长度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);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/xhs T(1,1,t+1)=2*Fo*(T(2,1,t)+T(1,2,t)+(

25、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;%下边界处温度计算式endfor 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;%左边界处温度计算式endfor m=2:mx-1for n=2:ny-1 T(m,n,t+1)=Fo*(T(m+1,n,t)+T(m,n

26、+1,t)+T(m-1,n,t)+T(m,n-1,t)+(1-4*Fo)*T(m,n,t)+ar*q*dt/k;%内部温度计算式 endendendfor 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,

27、(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,(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:ny T(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为热扩散率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(:,: