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

热传导与热辐射大作业报告热传导与热辐射大作业报告目录一、作业题目- 1 -二、作业解答- 2 -个人感想- 17 -附件.计算中所用程序- 18 -一、作业题目一矩形平板, ，内有均匀恒定热源，在及处绝热，在及处保持温度，初始时刻温度为，如右图1所示：1、求时，矩形区域内的温度分布的解析表达式；2、若，热传导系数，热扩散系数。请根据1中所求温度分布用MATLAB软件绘出下列结果，加以详细物理比较和分析：(a) 300s内，在同一图中画出点、（单位：m）温度随时间的变化；(b) 200s内，画出点、（单位：m）处，分别沿x、y方向热流密度值随时间的变化；(c) 画出时刻区域内的等温线；(d) 300s内，在同一图中画出点（单位：m）在分别等于，情况下的温度变化；(e) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为、和的温度、热流密度变化；(f) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化；3、运用有限差分法计算2中(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） 在中特征值问题的解可以直接从表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）其中0xa,0yb。当t=0时，上式变为： （1-25）其中0xa,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、若，热传导系数，热扩散系数。请根据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方向热流密度值随时间的变化； 图2.200s内x方向不同点的热流密度曲线（解析解） 图3.200s内y方向不同点的热流密度曲线（解析解）图4.200s内x方向不同点的热流密度曲线（数值解） 图5.200s内y方向不同点的热流密度曲线（数值解） 图7.不同点x方向热流密度数值解与解析解相对误差 图6.不同点x方向热流密度数值解与解析解相对误差分析：右边界（18,4）和（18,8）这两点开始时X方向两侧温差较大，故热流密度也会特别大，开始时内部温度较边界温度低，向内部导热，热流密度为负值。后来内部温度大于边界温度，向外散热，热流密度为正值。而上边界点温度相同，故在X方向不存在热传导，故导热系数为零。而中间点开始时从右向左导热，热流密度为负，随着边界层温度影响的深入，热流密度绝对值越来越大，但由于有内热源，会使此影响逐渐减弱，故在一段时间后待热流密度达到一个顶峰以后会逐渐减小，后来由于内热源的作用，导热由内向外进行，热流密度也由负值变为正值。Y方向分析类似。由于（9,6）离上边界更近，故沿Y方向达到的下边界峰值更大。（c）.画出时刻区域内的等温线； 图8.50s时区域内的等温线 图9.75s时区域内的等温线 图10.100s时区域内的等温线 图11.125s时区域内的等温线 图12.150s时区域内的等温线分析：开始时刻，尽管有内热源的存在，但边界温度比内部温度高，此时边界向内部传热，故开始时靠近边界的温度比内部高，这就是为什么50、75、100s时等温线呈现由坐下到右上温度逐渐升高。过一段时间后，中间部分由于内热源和边界热传导的共同作用，而坐下边界此时收到的内热源和边界热传导的作用小于中间部分，故造成了中间部分温度反而比其他部分高。一段时间后，内热源起主导作用，向外散热，这事等温线上的温度由左下到右上逐渐降低。（d）.300s内，在同一图中画出点（单位：m）在分别等于，情况下的温度变化；图13.不同内热源下温度变化曲线（解析解） 图14.不同内热源下温度变化曲线（数值解）图15.不同内热源下数值解与解析解相对误差分析：内热源越大，单位时间内内部产生的能量越多，节点温度升高的越快。在其它条件相同的情况下，内热源越大，最终内部温度也越高。开始时，由于温度变化剧烈，此时解析解和数值解的误差也相对较大，一段时间以后温度趋于稳定，这个时候相对误差也趋于一个较小的稳定值。（e）.300s内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为、和的温度、热流密度变化； 图16.不同导热数下温度变化曲线（解析解） 图17.不同导热数下温度变化曲线（数值解）图18.不同导热系数下数值解和解析解的相对误差 图20.不同导热系数下X方向热流密度曲线（数值解）图19.不同导热系数下X方向热流密度曲线（解析解） 图20.不同导热系数下X方向热流密数值解与解析解相对误差图22。不同导热系数下Y方向热流密度曲线（数值解）图21.不同导热系数下Y方向热流密度曲线（解析解） 图23.(9,6)点不同导热系数下y方向数值解和解析解的相对误差分析：导热系数K越大，内部温度越能快速的传递给外界，这就是问什么导热系数越大，节点最终温度低。根据热流密度方程，可知。K越大，热流密度越大，这就是为什么K越大，热流密度最低点峰值越大。而最后由于内热源相同，根据能量守恒，最后导热系数也必然趋近于一个定值。开始时由于温度变化剧烈，在不同的导热系数下同一点温度随之间变化值得数值解和解析解的相对误差较大，一段时间后温度趋于稳定，此时数值解和解析解的相对温差是一个较小值。（f）.300s内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化； 图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)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)个人感想经过一个多星期的连续奋战，终于搞定了这“万恶”的热传导与热辐射的大作业。首先真诚的感谢在作业中帮助过我的老师和同学。本来以为求温度场并不会是一件特别难的事情，可是等到实践时却发现里面有很多自己意想不到的困难。自己的MATLAB零基础确实也增添了不少困难。好不容易把程序编出来了，带入运行却是出问题了，总是比时间值少很多，花了一晚上一点一点的查却没有任何结果。知道第二天早上才发现是自己在循环中占用了原先定义的一个量。让人崩溃又让人欣喜：悲的是半天没有结果，喜的是终于找到了问题的根源。这样的事情还有很多很多。有时候为了查一个错误总需要花很长时间，但是经过奋战后终于把问题弄明白的那种欣喜确实很快乐的。在数值解的过程中，出现了一些令人感觉崩溃的问题。比如，步长取大了难以保证精度，取小了计算特别慢，而且出现一个让人再也做不下去的感觉“out of memory”。曾经一次计算了十几个小时最后得出了一个这样的结果，最后只能两者中和取，得出最终结果。从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为边界温度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*btm3*cosh(btm*b);endwt=(a2-x2)*g/(2*k)+T1-wtjfor 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*btn2)+g*gmv2/(k*(gmv2+btn2)*btn2)*cos(btn*x)*cos(gmv*y)*exp(-ar*(btn2+gmv2)*t); endendA(p,q,1)=wt+fwt;end endend附件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为时间节点数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+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 endend for i=1:sjz/xhs TZ(i+50*(xh-1),1)=T(31,21,(i-1)/dt+1);endfor m=1:mx for 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); 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: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*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;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:15 for k0=1:15 btn=(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(gmv*y)*exp(-ar*(btn2+gmv2)*t); endend 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方向总长度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)+(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 endend 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); %（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:mx for 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(:,:,1)=T0; for