一维非稳态导热问题的数值解

1、计算传热学程序报告题目：一维非稳态导热问题的数值解姓名：学号：学院：能源与动力工程学院专业：工程热物理日期：2014年5月25日一维非稳态导热问题数值解求解下列热传导问题:62T'xr=O (0x) a tT(x) = 0T(0) = l,T(L) = OL = .a = 11方程离散化对方程进行控制体积分得到:X非稳态项：选取T随X阶梯式变化，有(厂+A _ TJdx = (%7 - TlP )x扩散项：选取一阶导数随时间做显示变化，有(¥)(#) J力=(#);(¥)!J进一步取T随X呈分段线性变化，有竺 _TE-TPx e (x)e(OT _TP-TW整理可以

2、得到总的离散方程为:THT； T：-2T；+T：rx22计算空间和时间步长取空间步长为:h=LN网格FOUriCr数为:Fo=(小于o5时稳定)时间步长为：a3. 建立温度矩阵与边界条件T=OnCS(N+1,M+1)T(J)=Ti(初始条件温度都为0)T(l,：)=TO0力界条件x=()处温度为1)T(N+1,:)=TU(边界条件X=L处温度为0)4. 差分法求解温度由离散方程可得到：T=F.(TlE-TrIPV)-TtP转化为相应的温度矩阵形式：Tm.k + l) = FQ T(m + l,k) + T(m-ik)-2T(m,k) + T(myk)5. 输入界面考虑到方程的变量，采用inpu

3、tdig函数设置5个输入变量，对这5个变量设置了默认值，如图1所示。在计算中可以改变不同的数值，得到不同的结果，特别注 意稳定条件的临界值是().5。根据设置的默认值，得到的计算结果如图2所示。图1 matbb变量输入界面图2默认值的计算结果6结果分析根据上面的分析，给出了程序的输入界面，以及默认值状态下的数值解。可以通过改变不同的输入值，得到需要的分析结杲，总结出了下面4点结论：(I)取FR).48,得到一维非稳态导热结果如下图所示图2 Fll=O.48时一维非稳态导热从图中可以看出，对于长度L=I的细杆，初始时刻尸0时温度为0,边界条 件x=0时，T=I,边界条件x=l时，T=OO随着时间

4、的增加，温度从x=()通过导热 的形式传递到x=l,不同时刻不同位置杆的温度都不同，并且随着时间的增加, 杆的温度也逐渐增加。(2)取Fn=O.48,可以得到不同位置的温度响应曲线，如下图所示图中红色曲线代表x=()l位置的温度瞬态响应，黑色曲线代表x=0,2位置的温度瞬态响应，蓝色曲线代表x=()4位置的温度瞬态响应。从图中可以看出，随着X的增加，曲线与X轴的交点值越大，温度开始传递到该位置的所需的时间越 长。随着X的增加，温度响应曲线的变化速率越慢，最终的达到的温度也越低。(3)取FlI=O.25,得到不同位置的温度响应曲线如下图所示图中三条曲线分别是x=0.1, x=0.2, x=().

5、4位置的温度瞬态响应。与图3的 Fo=O.48迸行对比，两种情况下的珂值不同，值越大表明热扩散系数Q的值越 大。从图中可以看出热扩散系数对于导热的影响，F(I=O.25时，与F=0.48相比较, 各位置开始响应时所需的时间较长，而且各位置响应曲线的变化速率较小，最终 的达到的温度也较低，说明了热扩散系数越小，热传导越慢，传递效率越低。(4)取Fll-0.51 ,得到非稳定的数值解如图所示O 1002003004005006007008009001000时间t图6 Fo=0.51时不同X位直处的混度响应从图中可以看出，对于显示格式的离散方程，并不是所有的F。值都能得到有意义的解，必须要求F<

6、;0.5时才能得到稳定的数值解，当Fo>0.5时，会出现物理上不真实的解。附件:（matlab程序）function hcat_conducri（）n0 %一维齐次热传导方程%设亶输入界面OPtionS彳空间杆长LT空间点数Z ,时间点数MT扩散系数/稳定条件的值Fo(临界值0.5); topics维非稳态导热：标题栏显示IinCS= 1;%输入行为1行dcf=vooooo7i7o.48i% 默认值输入f=inputdlg()ptions,topic,lines,dc0i%y 入框设直L=CVal(fl);%设直输入值N=CVal(F2);MF 邮3);JI=CV 狄 1(F4);Fo=

7、CVaKf5);%FO的值必须小于0.5,小于0.5波动%计算空间步长与时间步长h=L/N;%空间步长xl=O:h:L;x=;n=Fo*h2ai%时间步长tm=n*M;%传导总时间tl=O：n：Im;t=t;%计算初始条件与边界条件Ti=x.*O;%初始条件To=l+t*0j%=0的边界条件Tc=t.*O;%x=L的边界条件%逹立温度矩阵TT=OnCS(N+1,M+1);T(:,I)=Ti;%第一列为初始条件T(1,:)=To;%第一行为X=O边界条件T(N+l,:)=Tc;%>后一行为X=L边界条件%利用差分法求解温度矩阵Tfor k=l:Mm=2;WhiIC m<=N;T(mR

8、+l)=Fo%T(m+ljk)+T(m")2*T(mrk)+T(myk);m=m+l;CndCIId%将时间空间的一维坐标转化为二维坐标PiX=mcshgrid(tl ,);%根据温度矩阵T绘图SUbPIot(2,2,1);mcsh(xy,T)；%三维图绘制(1,-1,1);%调整视图角度 山IUC非稳态导热)；图像若称 MlbUlC长度) 轴合称 ylabclf时间tt)%y轴君称 ZIabCIC温度Tt)；%Z轴各称 SUbPlot(2,2,2);A=T(II取矩阵第11列的值 PlOt(A>*);%二维曲线绘制 kgend(Q=(Ilr);%显示函数名称 UdCCX=0.1瞬态响应*)； XlabCIf时间 t*);ybbul(t温度 T); ais(0 1000 0 1)；%坐标轴敌值范围 SUbPIOt(2,2t3);B=T取矩阵第21列PIot(B,*k); ICgend(IB=O.2t);UtlCCX=0.2瞬态响应)； Xlabel(I时间 tf);Vlabclf 温度 T)； ais(O IoOOO 1);SUbPIot(2,2,4);C=T(41,:);%取矩阵第41列 Plot(A