热传导与热辐射大作业

1、热传导与热辐射大作业解答姓名：学号：1 问题描述一矩形平板, ，内有均匀恒定热源，在及处绝热，在及处保持温度，初始时刻温度为，如右图1所示：1、求时，矩形区域内的温度分布的解析表达式；2、若，热传导系数，热扩散系数。请根据1中所求温度分布用MATLAB软件绘出下列结果，加以详细物理比较和分析：(a) 300s内，在同一图中画出点、（单位：m）温度随时间的变化；(b) 200s内，画出点、（单位：m）处，分别沿x、y方向热流密度值随时间的变化；(c) 画出时刻区域内的等温线；(d) 300s内，在同一图中画出点（单位：m）在分别等于，情况下的温度变化；(e) 300s内，比较点(9,6) （单位

2、：m）在其它参数不变情况下热导率分别为、和的温度、热流密度变化；(f) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化；3、运用有限差分法计算2中(b)、(d)和(e)，并与解析解结果进行比较，且需将数值解与解析解的相对误差减小到1以下；4、附上源程序和个人体会；以报告形式整理上述结果，用A4纸打印上交。2数学描述及解析解该问题的数学描述为： 其中定解条件为： ； ； 该问题可分解为一个的非齐次稳态问题，和一个的齐次非稳态问题。2.1 稳态非齐次问题的求解稳态问题数学描述为： ； ； 求解该问题，定义新变量为 将式(2.3)代入问题(2.

3、2)可得： ； ；选常数,则(2.4)变为： ； ；这是不含内热源的二维稳态问题，用分离变量法求解。可分离为,分离后问题为： ；查M.N.奥齐西克教材中表2-2得到 ， ，其中 是方程 的正根，即 又有 解可取为的完全解由下式组成： 利用的非齐次边界条件可得 利用正交性求系数 ： 积分得到： 代入式(2.9)可得问题(2.5)的解为： 代入式(2.3)可以得到稳态问题的解为： 2.2 齐次非稳态问题的求解齐次问题的数学描述为： ； ；同样使用分离变量法求解。同时假定分离为 的形式。对应与函数 和 的分离方程为： ； ； 的解为： 则完全解可表示为： 当 时，代入初始条件： 根据M.N.奥齐西克

4、教材中表2-2可查得： ， ，其中 是方程 的正根，即 ， ，其中 是方程 的正根 利用这些特征函数的正交性，有 将查表所得数据代入求得 将上式代入(2.20)可得齐次问题的解为： 原问题的解可以由下式取得： 即 至此得到了原问题的解析解。3 数值方法求解3.1划分网格将矩形区域沿方向以间隔等分为个节点，沿方向以间隔等分为个节点，形成节点网络，其中任意一个节点的温度可表示为。时间以步长计算，即 写出节点方程：使用隐式有限差分方程。(1) 对于内部节点 若取，则方程整理为 (2) 对于绝热边界点在的绝热边界 若取，则方程整理为 同理在的绝热边界 (3) 对于定温边界附近的节点在定温边界附近，内部

5、节点方程仍适用，但此时可取 具体在边界内侧附近 并取，得 同样在边界内侧附近 (4) 关于四个角点在原点处 若取，则方程整理为 同样在时 当时 当时 3.2 确定相关数据在本题中取，则,。取。则傅里叶数 这里的傅里叶数取得较大，可以看出采用隐式方法求解的无条件收敛性质的便利性。由于存在定温边界，研究范围为 将二维量排列为一维形式，作如下变换 向量 则线性方程组的系数矩阵和常数向量的元素值分别可取 特别的，当，即时 当，即时 当，即时当，即时当，即时 当，即时当，即时当，即时由此确定系数矩阵和向量的的元素。3.3 进行数值求解根据初始条件可以得到。用MATLAB迭代计算矩阵方程 即可得到温度解。

6、计算所使用的程序附于文后。4 计算结果4.1 解析解计算结果(1) 300s内，在同一图中画出点、（单位：m）温度随时间的变化；温度变化图像绘于下图。所用程序见文后附录。从图中可以看出，各点温度变化呈单调上升趋势，最后趋于平缓。这是由于非稳态过程起始时刻区域内存在较大不平衡势差，因此引起温度剧烈变化，当过程进行一段时间后，逐渐趋于平衡，因此各点温度变化也逐渐平稳，趋于稳态值。(2) 200s内，画出点、（单位：m）处，分别沿x、y方向热流密度值随时间的变化；各点热流密度变化如下图所示。程序见文后附录。这里热流密度取与坐标轴正方向相同的方向为正。由于区域内部初始温度较低，热量从定温边界流入，因此

7、初始热流密度均为负值。随着过程进行，内部温度逐渐升高，由于内热源的作用，内部温度逐渐高于定温边界温度，因此热流密度变化为正值。(3) 画出时刻区域内的等温线； 各所求时间等温线图如上图所示。程序见文后附录。(4) 300s内，在同一图中画出点（单位：m）在分别等于，情况下的温度变化；在不同内热源强度下温度随时间的变化如下图所示。程序见文后附录。可以看到，内热源强度提高，区域的温度变化速率和最终能达到的稳态温度都有提高，这与事实是相符的。(5) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热导率分别为、和的温度、热流密度变化；不同热导率条件下温度和热流密度随时间变化情况如下图所

8、示。程序见文后附录。题干中给出了热导率和热扩散系数两个不相互独立的参数，这一问当热导率增加，其他值不发生变化，意味着热扩散系数不变，从而意味着区域的热容量必然增加。因此区域最终达到的稳态温度会随热导率提高有所下降。同时由于热容量提高，物质所传递的热量会有所减小，因此热流密度的极值会相应减小。(6) 300s内，比较点(9,6) （单位：m）在其它参数不变情况下热扩散系数分别为、和的温度、热流密度变化；不同热扩散系数条件下温度和热流密度随时间变化情况如下图所示。程序见文后附录。当热扩散系数升高时，区域达到稳态的时间缩短，这表现了热扩散系数反应热量传递快慢的物理意义。4.2 数值解计算结果对比在这

9、里比较时直接给出误差大小的图像，用于分析误差，而不再给出变化图像，因为他们相差很小。这里的误差是指相对误差，即 需要说明的是，本文控制误差大小在以下主要针对温度的计算误差，对于热流密度的误差涉及误差的传递，保证其在以下对计算量要求巨大，不适合在家用笔记本电脑上运行。以下计算所使用的程序均附在文后。(1) 不同点的热流密度随时间变化的误差图像可以发现图中存在许多剑锋，这是由于在热流密度曲线在接近和穿越零点时，由于其绝对值太小，以至于与MATLAB数据运算精度相近，因此误差不可避免的很大。在计算中应避免以绝对值较小的数值作为分母。(2) 内热源强度变化时温度随时间的变化可以看到温度的误差达到了精度

10、要求。(3) 热导率变化时，温度和热流密度随时间变化。5 心得体会开学之初老师即布置了本次大作业，并督促及早着手。因此我也较早的自学了分离变量解法，并参照课本给出了一个初步的解。但是在进行过程中屡屡发现各种的问题，有系数不正确，级数在初始时刻不收敛等问题，解析解几经易稿，痛苦不堪。直到课堂讲授分离变量法以后，使我对解法有了宏观的了解。在这个过程中，数值解是在课堂讲授数值解法之后完成的，首先使用稀网格计算，写好各问程序，在编程时又会遇到许多问题，程序频繁报错，修改过程也是极其繁琐，亦是苦不堪言。尤其是网格加密后出现了计算进行缓慢，内存不够，异常死机等现象，给我造成了很大困扰，最后采取分段计算的方

11、法，谨慎进行每一步计算才得以得到最后计算结果。解析解计算中的问题主要在于稳态解代入系数中积分的时候，这里使用了三角函数的正交性，使用后应将用替换以消去对的无穷求和，而反之若保留则会出现初始时刻级数不收敛和计算错误的情况。在解析解编程中，对于一元无穷级数求和，使用while条件循环即可控制精度，难点在于保证二重无穷级数和的收敛精度，我在这里使用了循环增加迭代次数的方式保证了级数和的精度。这在文后的程序中有所体现。本文数值解采用隐式算法，这就导致了在网格加密后矩阵过大，运算缓慢的现象。本文所使用的参数为，由此 。的选取利用了隐式格式的无条件收敛的优点，尽可能地增大了时间步长。在本人32位渣电脑的实

12、际计算中，计算这样一个时间步长所用时间接近，计算200秒需要10000个时间步长，这样一次计算将有约20000秒，大概有5个小时，若再提高精度，计算量将无法接受。本文的计算量对于笔记本电脑来说依然是很大的，经常出现内存不足的现象，因此我采用了分段计算的方法，每完成50秒计算一停，存储数据，这样停3到5次即可完成计算，随后将结果综合即可。最后附上源程序。计算源程序程序1 设置常数%定义各常数值为全局变量以方便调用%稳态参数global a; a=18; %mglobal b; b=12; %mglobal g0; g0=1; %w/m3global T1; T1=600; %kglobal T0

13、; T0=200; %kglobal k; k=1; %w/mkglobal alpha; alpha=0.8; %m2/sglobal eps; eps=1e-4; %控制精度%非稳态参数global dt;dt=0.02; %time stepglobal dx;dx=0.25; %spacingglobal Fo;Fo=alpha*dt/dx2; %0.256global t;%未定义时间，在不同问中另行定义程序2 解析解计算函数function T = Txyt(x,y,t)%解析解计算函数%返回值%使用函数以避免过多占用变量%引入稳态变量global a b g0 T1 T0 k a

14、lpha eps;%Ts function Ts = Tsxy( x,y ) summ=0; e=1; %e使用一次即释放 m=0; while abs(e/summ)>eps %求无穷级数，保证精度 bm=(2*m+1)*pi/2/a; e=(-1)m*cos(bm*x)*cosh(bm*y)/bm3/cosh(bm*b); summ=summ+e; m=m+1; end Ts=T1-g0*(x*x-a*a)/2/k-2*g0*summ/a/k; end%Th function Th = Thxyt( x,y,t ) function S=sumnp(N1,N2) %二重级数求和 fu

15、nction c=u(n,p) %级数一般项 yn=(2*n+1)*pi/(2*a); ep=(2*p+1)*pi/(2*b); c=(-1)(n+p). *(T0-T1)/yn/ep-g0/(k*yn3*ep)+(g0/k)*ep/(yn3*(ep2+yn2). *exp(-alpha*(yn2+ep2)*t). *cos(ep*y)*cos(yn*x); end %求和 S=0; for i=0:N1 for j=0:N2 S=S+u(i,j); end end end %保证级数求和的收敛精度 Th=0; e=1; N=20; while abs(e/Th)>eps Te=Th;

16、Th=sumnp(N,N)*(4/a/b); e=Th-Te; N=N+10; end end%主函数if t=0 T=200; %避免过多计算elseif x=18|y=12 T=600;else T=Tsxy(x,y)+Thxyt(x,y,t);endend程序3热流密度计算函数function qx,qy=flux(x,y,t)%计算解析解热流密度global a b g0 T1 T0 k alpha eps; function qx=fx(x,y,t) %summ summ=0; e=1; m=0; while abs(e/summ)>eps %求无穷级数，保证精度 bm=(2*

17、m+1)*pi/2/a; e=(-1)m*(-sin(bm*x)*cosh(bm*y)/bm2/cosh(bm*b); summ=summ+e; m=m+1; end %sumnp function sumnp=snp(N1,N2) %二重级数求和 function c=u(n,p) %级数一般项 yn=(2*n+1)*pi/(2*a); ep=(2*p+1)*pi/(2*b); c=(-1)(n+p). *(T0-T1)/yn/ep-g0/(k*yn3*ep)+(g0/k)*ep/(yn3*(ep2+yn2). *exp(-alpha*(yn2+ep2)*t). *cos(ep*y)*yn*

18、(-sin(yn*x); end %求和 sumnp=0; for i=0:N1 for j=0:N2 sumnp=sumnp+u(i,j); end end end %精度 sumnp=0; e=1; N=20; while abs(e/sumnp)>eps se=sumnp; sumnp=snp(N,N); e=sumnp-se; N=N+10; end %qx if y=12 qx=0; else qx=g0*x+2*g0*summ/a-4*k*sumnp/a/b; end end function qy=fy(x,y,t) %summ summ=0; e=1; m=0; whil

19、e abs(e/summ)>eps %求无穷级数，保证精度 bm=(2*m+1)*pi/2/a; e=(-1)m*cos(bm*x)*sinh(bm*y)/bm2/cosh(bm*b); summ=summ+e; m=m+1; end %sumnp function sumnp=snp(N1,N2) %二重级数求和 function c=u(n,p) %级数一般项 yn=(2*n+1)*pi/(2*a); ep=(2*p+1)*pi/(2*b); c=(-1)(n+p). *(T0-T1)/yn/ep-g0/(k*yn3*ep)+(g0/k)*ep/(yn3*(ep2+yn2). *ex

20、p(-alpha*(yn2+ep2)*t). *ep*(-sin(ep*y)*cos(yn*x); end %求和 sumnp=0; for i=0:N1 for j=0:N2 sumnp=sumnp+u(i,j); end end end %精度 sumnp=0; e=1; N=20; while abs(e/sumnp)>eps se=sumnp; sumnp=snp(N,N); e=sumnp-se; N=N+10; end %qy if x=18 qy=0; else qy=2*g0*summ/a-4*k*sumnp/a/b; end endqx=fx(x,y,t);qy=fy(

21、x,y,t);end程序4 数值解计算函数function Tnum=Tn()%数值非稳态解%返回三维数组%定义各变量global a b g0 T1 T0 k; %引入稳态变量%引入非稳态变量global Fo dt dx tM=a/dx;N=b/dx;%确定A矩阵及向量bA=zeros(M*N,M*N);B=zeros(M*N,1);for m=0:M-1 for n=0:N-1 i=m+M*n; %(0,M*N-1) flag=1; switch flag case m=0&&n=0&&n=N-1 i=i+1; A(i,i)=1+4*Fo; A(i,i-M

22、)=-Fo; A(i,i+1)=-2*Fo; A(i,i+M)=-Fo; B(i)=Fo*g0*dx2/k; case m=0&&m=M-1&&n=0 i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i+1)=-Fo; A(i,i+M)=-2*Fo; B(i)=Fo*g0*dx2/k; case m=M-1&&n=0&&n=N-1 i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i-M)=-Fo; A(i,i+M)=-Fo; B(i)=Fo*T1+Fo*g0*dx2/

23、k; case m=0&&m=M-1&&n=N-1 i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i-M)=-Fo; A(i,i+1)=-Fo; B(i)=Fo*T1+Fo*g0*dx2/k; case m=0&&n=0 i=i+1; A(i,i)=1+4*Fo; A(i,i+1)=-2*Fo; A(i,i+M)=-2*Fo; B(i)=Fo*g0*dx2/k; case m=0&&n=N-1 i=i+1; A(i,i)=1+4*Fo; A(i,i-M)=-Fo; A(i,i+1)=-2*Fo;

24、B(i)=Fo*T1+Fo*g0*dx2/k; case m=M-1&&n=0 i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i+M)=-2*Fo; B(i)=Fo*T1+Fo*g0*dx2/k; case m=M-1&&n=N-1 i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i-M)=-Fo; B(i)=2*Fo*T1+Fo*g0*dx2/k; otherwise i=i+1; A(i,i)=1+4*Fo; A(i,i-1)=-Fo; A(i,i-M)=-Fo; A(i,i+1)=-Fo;

25、A(i,i+M)=-Fo; B(i)=Fo*g0*dx2/k; end endend%Tp为p时刻解向量Tp=zeros(M*N,1);Tp=Tp+T0;%迭代for p=1:t/dt %不存t=0 Tp=A(Tp+B); %矩阵求解 S(:,p)=Tp; %解的形式为S(i,t) %变换为(x,y,t)形式 for i=0:M*N-1 n=floor(i/M); m=i-M*n; Tnum(m+1,n+1,p)=S(i+1,p); %Tnum即为所求数值解 end Tnum(M+1,:,p)=600; Tnum(:,N+1,p)=600; pendend程序5 第二题各小问%第二题第一问 求

26、温度变化clearclcvart=300; %求解时间for p=1:t Ta(1,p)=Txyt(0,4,p); Ta(2,p)=Txyt(0,8,p); Ta(3,p)=Txyt(6,0,p); Ta(4,p)=Txyt(12,0,p); Ta(5,p)=Txyt(9,6,p); p %观察计算进行情况endfigure(1);plot(1:t,Ta);xlabel('时间t/s');ylabel('温度T/K');legend('(0,4)','(0,8)','(6,0)','(12,0)',

27、'(9,6)')grid on; %第二题第二问 求热流密度vart=200; %求解时间for p=1:t qx(1,p),qy(1,p)=flux(18,4,p); qx(2,p),qy(2,p)=flux(18,8,p); qx(3,p),qy(3,p)=flux(6,12,p); qx(4,p),qy(4,p)=flux(12,12,p); qx(5,p),qy(5,p)=flux(9,6,p); p %控制台监视end%坐标轴方向为正%绘图figure(2);plot(1:t,qx);xlabel('时间 t/s');ylabel('x方向热流

28、密度 qx/ W/m2');legend('(18,4)','(18,8)','(6,12)','(12,12)','(9,6)')grid on;figure(3);plot(1:t,qy);xlabel('时间 t/s');ylabel('y方向热流密度 qy/ W/m2');legend('(18,4)','(18,8)','(6,12)','(12,12)','(9,6)')grid on

29、; %第二题第三问 求等温线vart=50,75,100,125,150; %求解时间%求温度分布for p=1:5 for x=0:0.5:a for y=0:0.5:b Tc(x*2+1,y*2+1,p)=Txyt(x,y,t(p); end end t(p) %绘图 figure(p); cs,h=contourf(0:0.5:a,0:0.5:b,Tc(:,:,p)'); clabel(cs,h); colorbar('EastOutside'); xlabel('x /m'); ylabel('y /m');end %第二题第四问

30、 求g0变化vart=300; %求解时间for i=1:3 g0=i; for p=1:t T90(i,p)=Txyt(9,0,p); p endendfigure(4);plot(1:t,T90);xlabel('时间t/s');ylabel('温度T/K');legend('g0=1','g0=2','g0=3')grid on; %第二题第五问 求k变化vart=300; %求解时间for i=1:3 k=i*0.5; for p=1:t T96(i,p)=Txyt(9,6,p); qx96(i,p),qy

31、96(i,p)=flux(9,6,p); p endend%figure(5);plot(1:t,T96);xlabel('时间t/s');ylabel('温度T/K');legend('k=0.5','k=1','k=1.5')grid on;%figure(6);plot(1:t,qx96);xlabel('时间t/s');ylabel('x方向热流密度 qx/ W/m2');legend('k=0.5','k=1','k=1.5'

32、)grid on;%figure(7);plot(1:t,qy96);xlabel('时间t/s');ylabel('y方向热流密度 qy/ W/m2');legend('k=0.5','k=1','k=1.5')grid on; %第二题第六问 求alpha变化vart=300; %求解时间for i=1:3 alpha=i*0.4; T96(i,1)=200; qx96(i,1)=0; qy96(i,1)=0; for p=1:t T96(i,p+1)=Txyt(9,6,p); qx96(i,p+1),qy96

33、(i,p+1)=flux(9,6,p); p endend%figure(8);plot(0:t,T96);xlabel('时间t/s');ylabel('温度T/K');legend('alpha=0.4','alpha=0.8','alpha=1.2');grid on;%figure(9);plot(0:t,qx96);xlabel('时间t/s');ylabel('x方向热流密度 qx/ W/m2');legend('alpha=0.4','alpha

34、=0.8','alpha=1.2');grid on;%figure(10);plot(0:t,qy96);xlabel('时间t/s');ylabel('y方向热流密度 qy/ W/m2');legend('alpha=0.4','alpha=0.8','alpha=1.2');grid on;程序6 第三题各小问%第二问 求热流密度clearclcvart=200; %求解时间%数值解Tnum=Tn();for p=1:t %先求出温度值 %所求点温度 Tbn(1,p)=Tnum(18/

35、dx+1,4/dx+1,p/dt); Tbn(2,p)=Tnum(18/dx+1,8/dx+1,p/dt); Tbn(3,p)=Tnum(6/dx+1,12/dx+1,p/dt); Tbn(4,p)=Tnum(12/dx+1,12/dx+1,p/dt); Tbn(5,p)=Tnum(9/dx+1,6/dx+1,p/dt); %所求点x方向附近温度 Tbnx(1,p)=Tnum(18/dx,4/dx+1,p/dt); Tbnx(2,p)=Tnum(18/dx,8/dx+1,p/dt); Tbnx(3,p)=Tnum(6/dx,12/dx+1,p/dt); Tbnx(4,p)=Tnum(12/dx

36、,12/dx+1,p/dt); Tbnx(5,p)=Tnum(9/dx,6/dx+1,p/dt); %所求点y方向附近温度 dx=dy Tbny(1,p)=Tnum(18/dx+1,4/dx,p/dt); Tbny(2,p)=Tnum(18/dx+1,8/dx,p/dt); Tbny(3,p)=Tnum(6/dx+1,12/dx,p/dt); Tbny(4,p)=Tnum(12/dx+1,12/dx,p/dt); Tbny(5,p)=Tnum(9/dx+1,6/dx,p/dt);end%解析解for p=1:t qx(1,p),qy(1,p)=flux(18,4,p); qx(2,p),qy(

37、2,p)=flux(18,8,p); qx(3,p),qy(3,p)=flux(6,12,p); qx(4,p),qy(4,p)=flux(12,12,p); qx(5,p),qy(5,p)=flux(9,6,p); p %控制台监视endsave Tnumb;%坐标轴方向为正qxn=-k*(Tbn-Tbnx)/dx;qyn=-k*(Tbn-Tbny)/dx;%ex=(qxn-qx)./qx;ey=(qyn-qy)./qy;%绘图figure(1);plot(1:t,ex);xlabel('时间 t/s');ylabel('x方向热流密度误差 epsx');legend('(18,4)','(18,8)','(6,12)','(12,12)','(9,6)')grid on;figure(2);plot(1:t,ey);xlabel('时间 t/s');ylabel('y方向热流密度误差 epsy');legend(