计算传热学课程设计摘要：利用数值模拟技术构建应拟实脸是提高传热学教学及实脸效率、降低实脸成本的一条有效途径。
本学期基于Mtalab强大的数位计算功能、圈形处理能力,学习了传热学应用模拟实。
本文主要介绍了在本学期，老师在课堂上所讲的一些基本内容和我在本学期主要学习到了那些东西和课下作业在计算机上进行的数值模拟。
关键字：计算传热学Matlab 数值模拟一、介绍计算传热学又称数值传热学：是指对描写流动与传热问题的控制方程采用数值解法通过计算机予以求解的一门学科。
数值解法是一种离散近似的计算方法。
它所能获得的解不像分析解那样是被研究区域中未知的量的连续函数，而只是某些代表性的点（称为节点）上的近似解【1】。
MATLAB是在上个世纪80年代,美国几个主要从事数学研究的学者,针对矩阵运算的复杂性,开发的专门用于矩阵运算的程序,结果应用效果很好【2】。
二、内容本学期只要学习了五讲内容。
第一讲是微分方程近似解的概述，主要内容是一维翅片传热问题及其解析解、求连续的近似解析函数的方法包括了RITZ法，权余法，配置法，以及迦辽金法。
第二讲是导热的有限体积法，主要内容有一维导热问题的数值解和一维非稳态问题。
要求一维导热问题的数值解首先要将求解区域离散化，离散的方法有外节点法和内节点法；其次在把源项线性化控制方程的离散化；再次进行交界面参数的计算和跃界面的处理，交界面参数计算的方法有线性插值法、调和平均法、待求变量插值和Kirchhoff变换法；最后进行边界条件的处理和差分方程的求解。
一维非稳态问题主要需要进行非稳态项的处理、控制方程的离散化、边界条件的处理和求解。
第三讲是离散方程的求解，主要内容是一维离散方程的求解和对位离散方程的求解。
一维离散方程求解的方法有TDMA、Jacobi迭代、Gaus-Seidel迭代和SOR迭代；多维离散方程的求解方法有点迭代、线迭代和ADI。
第四讲对流-扩散的有限体积法，主要内容是稳定的一维对流-扩散、非稳定的一维对流扩散和多维的对流扩散。
解决稳定的和非稳定的一维对流-扩散问题主要运用的格式有中心差分格式、一阶逆风格式、混合格式、乘方格式和指数格式。
这几种格式各有优缺点，在运用是要根据具体情况选用，以使模拟的精度更高，更接近具体的数据，以便更好的解决在工程运用中的具体问题。
第五讲是对流-扩散的有限体积法,只要内容是在假设流畅已知的条件下，解决稳定的一维对流-扩散、非稳定的一维对流扩散和多维的对流扩散问题。
他们是通过加错网格的压力修正方法和同位网格的压力修正方法。
第六讲是流场计算，这一讲是通过有限体积计算不可压缩流体的流场。
求解不可压缩流动流场数值解的方法有联立求解个变量代数方程组的方法和分离式求解代数方程组的方法。
而流场求解的关键问题是压力梯度的离散和压力的求解。
三、课后作业地源热泵柱热源模型，孔壁半径65mm，岩土导热系数1.8，体积比热3x106，孔壁Q=50W/m无穷远处温度18 给出导热微分方程，用一维数值模型求岩土径向温度数值解与解析解比较这个问题我们可以通过matlab得出各点温度分布图形。
在本题中节点的划分我使用的是外节点法，采用元体平衡法进行方程分析，分别对节点1、2、…、n-1、n进行分析，此处为计算方便直接将节点N的温度等于无穷远处温度。
通过matlab得出的图形如下所示：隐式格式显示格式通过matlab模拟的两个图形的对比可知显式格式和隐式格式的图形对比可知在两个温度模拟的差别并不大，说明两种格式都能反映岩土径向温度的变化情况。
参考文献「1〕陶文栓.传热学[M〕.北京:高等教育出版社,2005.「2〕彭芳辟.数学物理方程的Matlab解法与可视化〔M」.北京:清华大学出版社,2004.附录非稳态隐式clear;qb=50;k=1.8;ri=0.065;T8=18;ro=2.065;rc=3000000;Sc=0;Sp=0;NumNode=100;dx=(ro-ri)/(NumNode-1);Xo=ro;Xi=ri;X=[Xi:dx:Xo];tstep=1000;dtmax=rc*dx*dx/2/k;dt=60;Ar=ri;AWi(1)=0;AEi(1)=(k*Ar*2)/dx+k;APi0(1)=rc*((Ar+dx/2)*(Ar+dx/2)-Ar*Ar)/dt; %°ë¸ö¿ØÖÆÌåAPi(1)=APi0(1)+AWi(1)+AEi(1)-Sp*((Ar+dx/2)^2-Ar^2)*dx; %ÒþʽBi0(1)=((Ar+dx/2)*(Ar+dx/2)-Ar*Ar)*Sc+qb/3.14;%½Úµã2 ---- NumNode-1·½³ÌϵÊý ´Ë²¿·Ö²»Ëæʱ¼ä±ä»¯for i=2:NumNode-1Ar=ri+(i-1)*dx; %iµãÃæ»ýÒò×Ó% Ar+dx/2; %e½çÃæÃæ»ýÒò×Ó% Ar-dx/2; %w½çÃæÃæ»ýÒò×ÓAWi(i)=(2*k*Ar)/dx-k;AEi(i)=(2*k*Ar)/dx+k;APi0(i)=(rc*((Ar+dx/2)*(Ar+dx/2)-(Ar-dx/2)*(Ar-dx/2)))/dt;%AP(i)=AP0(i); %ÏÔʽAPi(i)=APi0(i)+AWi(i)+AEi(i)-Sp*Ar*dx; %ÒþʽBi0(i)=((Ar+dx/2)*(Ar+dx/2)-(Ar-dx/2)*(Ar-dx/2))*Sc;%BP(i)=B0(i)+(AP0(i)-AW(i)-AE(i)+dx*dx*Sp)*T0(i);endAr=ri+(NumNode-1)*dx;AWi(NumNode)=0;AEi(NumNode)=0;APi0(NumNode)=1;%AP(NumNode)=AP0(NumNode);APi(NumNode)=APi0(NumNode);Bi0(NumNode)=0;%B0=B0';TT0(1:NumNode)=18; %ÈÎÒâ³õʼֵfor t=1:tstepBB(1,1)=Bi0(1)+APi0(1)*TT0(1); %BBÿһ²½±ä»¯AA(1,1:2)=[-APi(1) AEi(1)]; %AE*Tw-AP*Tp+AE*TE=-BP %µÚ2~NumNode-1½Úµãfor i=2:NumNode-1AA(i,i-1:i+1)=[AWi(i) -APi(i) AEi(i)]; %AE*Tw-AP*Tp+AE*TE=-BP BB(i,1)=Bi0(i)+APi0(i)*TT0(i);end%AE*Tw-AP*Tp+AE*TE=-BPAA(NumNode,NumNode-1:NumNode)=[AWi(NumNode) -APi(NumNode)]; BB(NumNode,1)=Bi0(NumNode)+APi0(NumNode)*18;BB=-BB; %AE*Tw-AP*Tp+AE*TE=-BPT2=inv(AA)*BB; %Ö±½Ó½â AA*T=BB £¿£¿£¿½á¹û²»¶Ô£¿%A3=[AW' -APi' AE']%T2=tdma1(A3,BB);TT0=T2'; %¸üгõÖµTT(t,:)=T2'; %±£´æµ±Ç°Öµ Òþʽ½á¹ûendfigure;plot(X,TT(500,:),'-x');figure;plot([1:500]*dt/60,TT([1:500],1),'k-');非稳态显式clear;qb=50;k=1.8;ri=0.065;T8=18;ro=1.065;rc=3000000;NumNode=101;dx=(ro-ri)/(NumNode-1);Xo=ro;Xi=ri;X=[Xi:dx:Xo];tstep=1000;dtmax=rc*dx*dx/2/k;dt=60;Ar=ri+dx/2;AW(1)=0;AE(1)=(2*3.14*k*Ar)/dx;AP0(1)=(rc*(Ar*Ar-ri*ri)*3.14)/dt; AP(1)=AP0(1);B0(1)=qb;for i=2:NumNode-1Ar=ri+(i-1)*dx;AW(i)=(2*3.14*k*(Ar-dx/2))/dx;AE(i)=(2*3.14*k*(Ar+dx/2))/dx;AP0(i)=(rc*2*3.14*Ar*dx)/dt;AP(i)=AP0(i);B0(i)=0;endAr=ri+(NumNode-1)*dx;AW(NumNode)=(2*3.14*k*(Ar-dx/2))/dx;AE(NumNode)=(2*3.14*k*(Ar+dx/2))/dx;AP0(NumNode)=(rc*2*3.14*Ar*dx)/dt;AP(NumNode)=AP0(NumNode);B0(NumNode)=0;B0=B0';T0(1:NumNode)=18;for t=1:tstepAr=ri+dx/2;yy=AP0(1)-AW(1)-AE(1);B(1,1)=B0(1)+yy*T0(1);T1(1)=(AE(1)*T0(1+1)+B(1))/AP(1); for i=2:NumNode-1Ar=ri+(i-1)*dx;yy=AP0(i)-AW(i)-AE(i);B(i,1)=B0(i)+yy*T0(i);T1(i)=(AW(i)*T0(i-1)+AE(i)*T0(i+1)+B(i))/AP(i);endAr=ri+(NumNode-1)*dx;yy=AP0(NumNode)-AW(NumNode)-AE(NumNode);B(NumNode,1)=B0(NumNode)+yy*T0(NumNode);T1(NumNode)=(AW(NumNode)*T0(NumNode-1)+AE(NumNode)*T8+B(NumNode))/AP(NumNode); T0=T1;T(t,:)=T1;endfigure;hold onplot(X,T(1,:),'k-');plot(X,T(300,:),'m--*');plot(X,T(600,:),'r:p');plot(X,T(900,:),'g-v');figure;plot([1:500]*dt/60,T([1:500],1),'k-');。