《计算传热学》由会员分享，可在线阅读，更多相关《计算传热学（6页珍藏版）》请在金锄头文库上搜索。
1、计算传热学课程设计摘要：利用数值模拟技术构建应拟实脸是提高传热学教学及实脸效率、降低实脸成本的一条 有效途径。本学期基于Mtalab强大的数位计算功能、圈形处理能力学习了传热学应用模拟 实。本文主要介绍了在本学期，老师在课堂上所讲的一些基本内容和我在本学期主要学习到 了那些东西和课下作业在计算机上进行的数值模拟。关键字：计算传热学 Matlab 数值模拟一、介绍计算传热学又称数值传热学：是指对描写流动与传热问题的控制方程采用数值解法通过 计算机予以求解的一门学科。数值解法是一种离散近似的计算方法。它所能获得的解不像分 析解那样是被研究区域中未知的量的连续函数，而只是某些代表性的点（称为节点）上
2、的近 似解【1】MATLAB是在上个世纪80年代，美国几个主要从事数学研究的学者,针对矩阵运算 的复杂性,开发的专门用于矩阵运算的程序,结果应用效果很好【2】。二、内容本学期只要学习了五讲内容。第一讲是微分方程近似解的概述，主要内容是一维翅片传 热问题及其解析解、求连续的近似解析函数的方法包括了 RITZ 法，权余法，配置法，以及 迦辽金法。第二讲是导热的有限体积法，主要内容有一维导热问题的数值解和一维非稳态问题。要 求一维导热问题的数值解首先要将求解区域离散化，离散的方法有外节点法和内节点法；其 次在把源项线性化控制方程的离散化；再次进行交界面参数的计算和跃界面的处理，交界面 参数计算的方法
3、有线性插值法、调和平均法、待求变量插值和Kirchhoff变换法；最后进行 边界条件的处理和差分方程的求解。一维非稳态问题主要需要进行非稳态项的处理、控制方 程的离散化、边界条件的处理和求解。第三讲是离散方程的求解，主要内容是一维离散方程的求解和对位离散方程的求解。一 维离散方程求解的方法有TDMA、Jacobi迭代、Gaus-Seidel迭代和SOR迭代；多维离散方程 的求解方法有点迭代、线迭代和ADI。第四讲对流-扩散的有限体积法，主要内容是稳定的一维对流-扩散、非稳定的一维对流 扩散和多维的对流扩散。解决稳定的和非稳定的一维对流-扩散问题主要运用的格式有中心 差分格式、一阶逆风格式、混合
4、格式、乘方格式和指数格式。这几种格式各有优缺点，在运 用是要根据具体情况选用，以使模拟的精度更高，更接近具体的数据，以便更好的解决在工 程运用中的具体问题。第五讲是对流-扩散的有限体积法,只要内容是在假设流畅已知的条件下，解决稳定的一 维对流-扩散、非稳定的一维对流扩散和多维的对流扩散问题。他们是通过加错网格的压力 修正方法和同位网格的压力修正方法。第六讲是流场计算，这一讲是通过有限体积计算不可压缩流体的流场。求解不可压缩流 动流场数值解的方法有联立求解个变量代数方程组的方法和分离式求解代数方程组的方法。 而流场求解的关键问题是压力梯度的离散和压力的求解。三、课后作业地源热泵柱热源模型，孔壁半
5、径65mm，岩土导热系数1.8,体积比热3x106,孔壁Q=50W/m 无穷远处温度 18 给出导热微分方程，用一维数值模型求岩土径向温度 数值解与解析解比较这个问题我们可以通过matlab得出各点温度分布图形。在本题中节点的划分我使用的是外 节点法，采用元体平衡法进行方程分析，分别对节点1、2、n-1、n进行分析，此处为 计算方便直接将节点 N 的温度等于无穷远处温度。通过 matlab 得出的图形如下所示：通过 matlab 模拟的两个图形的对比可知显式格式和隐式格式的图形对比可知在两个温度模 拟的差别并不大，说明两种格式都能反映岩土径向温度的变化情况。参考文献1陶文栓传热学M北京:高等教
6、育出版社,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;%。乙00亜i
7、吕APi(1)=APi0(1)+AWi(1)+AEi(1)Sp*(Ar+dx/2)八2-Ar八2)*dx;%Q总坯Bi0(1)=(Ar+dx/2)*(Ar+dx/2)-Ar*Ar)*Sc+qb/3.14;%Upa2 NumNode-13IIpEy E2d2EEaa_for i=2:NumNode-1Ar=ri+(i-1)*dx;%ipaAyOdxO% Ar+dx/2; %eQAAyddxO% Ar-dx/2; %wQAAyddxOAWi(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)
8、*(Ar-dx/2)/dt;%AP(i)=AP0(i); %IOEAPi(i)=APi0(i)+AWi(i)+AEi(i)-Sp*Ar*dx; %0qE坯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);.y (T y oog ： TJ ) lly09/4P oog ： TJ ) 4OTd I ejn5Tj(,X-,y(： yoog)llyX)40ydI ejn5Tj puenT?340 ri0o5+ri,3+%乙 1= ( ： y4) 11riO
9、QEVi3nT %乙 1=011-(aa yV)TUIP4 =/ ,3V ,TdV- ,MV=V%?30EzOTe?3?3?3 aa=lVV eot+O% -99 (W) aut二乙丄aa-=aada-=313V+didV-M1i13V%81* (epoNumN) 0T3V+ (epoNumN) 0T9= (T VpoNumN) 99i (opoNumN) yay- (opoNumN) tmv = (opoNumN : 1-epoNumN opoNumN) W da-=313V+didV-M1i13V% puoda-=313V+didV-M1i13V% (T) T3Vda-=313V+dldV
10、-M1l13V%_e+zo?vaa%(T)Oil*(T)0T3V+(T)0T9=(T T)99(T)TdV- (T)TMV=(T+T:T-TyT)VVT-opoNumN :乙=T JOj erinT-spoNumN 乙 n% (T)T3V (T)TdV-J = (Z：T yT)W (T)Oil*(T)0T3V+(T)0T9=(T yT)99 deq.sq.: i=l- riO3QEeoi3% 8T= (opoNumN：T) Oil-.09=09%0= (epoNumN) 0T9 -(opoNumN) 0TdV= (epoNumN) yav-(epoNumN) 0dV= (epoNumN) dV
11、%1= (epoNumN) OTdV0= (epoNumN) thv0= (epoNumN) tmy (l-epoNwnN)+tj=jvk=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
12、;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(NumN
13、ode)=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-);
