计算传热学程序报告
题目：一维非稳态导热问题的数值解
姓名：
学号：学院：能源与动力工程学院专业：工程热物理日期：2014年5月25 日
一维非稳态导热问题数值解
求解下列热传导问题:
1. 方程离散化
对方程进行控制体积分得到: 非稳态项：选取T随x阶梯式变化，有扩散项：选取一阶导数随时间做显示变化，有进一步取T随x呈分段线性变化，有
T ()e
x T
E
T
P (T \ T P T W
，( )w (x)x ( X)
整理可以得到总的离散方程为:
2. 计算空间和时间步长
取空间步长为:
h=L/N
网格Fourier数为：
F。
一^ —V (小于时稳定)
x x
时间步长为：
3. 建立温度矩阵与边界条件
T=o nes(N+1,M+1)
T(:,1)=Ti(初始条件温度都为0)
T(1,:)=To(边界条件x=0处温度为1)
T(N+1,:)=Te(边界条件x=L处温度为0)
4. 差分法求解温度
由离散方程可得到：
转化为相应的温度矩阵形式：
5. 输入界面
考虑到方程的变量，采用inputdlg函数设置5个输入变量，对这5个变量设置了默认值，如图1所示。
在计算中可以改变不同的数值，得到不同的结果，特别注意稳定条件的临界值是。
根据设置的默认值，得到的计算结果如图2所示。
图1matlab变量输入界面
图2 默认值的计算结果
6. 结果分析根据上面的分析，给出了程序的输入界面，以及默认值状态下的数值解。
可以通过改变不同的输入值，得到需要的分析结果，总结出了下面4 点结论：
（1）取F o=，得到一维非稳态导热结果如下图所示
图2F。
二时一维非稳态导热
从图中可以看出，对于长度L=1 的细杆，初始时刻t=0 时温度为0，边界条件
x=0时，T=1,边界条件x=1时，T=0。
随着时间的增加，温度从x=0通过导热的形式传递到x=1，不同时刻不同位置杆的温度都不同，并且随着时间的增加，杆的温度也逐渐增加。
（2）取F o=，可以得到不同位置的温度响应曲线，如下图所示
图3F o=时不同X位置处的温度响应
图中红色曲线代表x=位置的温度瞬态响应，黑色曲线代表x=位置的温度瞬态响应，蓝色曲线代表X=位置的温度瞬态响应。
从图中可以看出，随着X的增加，曲线与X 轴的交点值越大，温度开始传递到该位置的所需的时间越长。
随着x 的增加，温度响应曲线的变化速率越慢，最终的达到的温度也越低。
（3）取F o=，得到不同位置的温度响应曲线如下图所示
图4F o=时不同X位置处的温度响应
图中三条曲线分别是X=，X=，x=位置的温度瞬态响应。
与图3的F o=进行对比，两种情况下的F o值不同，F o值越大表明热扩散系数的值越大。
从图中可以
看出热扩散系数对于导热的影响，尸0=时，与F o=相比较，各位置开始响应时所需的时间较长，而且各位置响应曲线的变化速率较小，最终的达到的温度也较低，说明了热扩散系数越小，热传导越慢，传递效率越低。
（4）取F o=，得到非稳定的数值解如图所示
图5F o二时一维非稳态导热
图6F o=时不同X位置处的温度响应
从图中可以看出，对于显示格式的离散方程，并不是所有的F o值都能得到有意义的解，必须要求F o＜时才能得到稳定的数值解，当F o＞时，会出现物理上不真实的解。
附件：( matlab 程序) functionheat_conduction()% 一维齐次热传导方程%设置输入界面
options={' 空间杆长L',' 空间点数N',' 时间点数M',' 扩散系数a','稳定条件的值Fo(临界值',};
topic=' 一维非稳态导热';% 标题栏显示
lines=1;% 输入行为1 行def={'1','100','1000','1',''};% 默认值输入
f=inputdlg(options,topic,lines,def);% 输入框设置
L=eval(f{1});% 设置输入值
N=eval(f{2});
M=eval(f{3});
a=eval(f{4});
Fo=eval(f{5});%Fo 的值必须小于，小于波动%计算空间步长与时间步长
h=L/N;%空间步长
x1=0:h:L;
x=x1';
n=Fo*h A2/a;%时间步长
tm=n*M;%传导总时间t1=0:n:tm; t=t1'; %计算初始条件与边界条件
Ti=x.*0;% 初始条件To=1+t.*0;%x=0 的边界条件Te=t.*0;%x=L 的边界条件%建立温度矩阵T T=ones(N+1,M+1);
T(:,1)=Ti;% 第一列为初始条件T(1,:)=To;% 第一行为x=0 边界条件T(N+1,:)=Te;% 最后一行为x=L 边界条件%利用差分法求解温度矩阵T fork=1:M m=2;
whilem<=N;
T(m,k+1)=Fo*(T(m+1,k)+T(m-1,k)-2*T(m,k))+T(m,k); m=m+1;
end end %将时间空间的一维坐标转化为二维坐标
[Y,X]=meshgrid(t1,x);
%根据温度矩阵T 绘图subplot(2,2,1);
mesh(X,Y,T);% 三维图绘制view([1,-1,1]);% 调整视图角度title(' 非稳态导热');% 图像名称xlabel(' 长度x');%x 轴名称ylabel(' 时间
t');%y 轴名称zlabel(' 温度T');%z 轴名称subplot(2,2,2);
A=T(11,:);% 取矩阵第11 列的值plot(A,'r');% 二维曲线绘制
legend('A=');% 显示函数名称title('x= 瞬态响应');
xlabel(' 时间t');
ylabel(' 温度T'); axis([0100001]);% 坐标轴数值范围subplot(2,2,3); B=T(21,:);% 取矩阵第21 列plot(B,'k');
legend('B=');
title('x= 瞬态响应'); xlabel(' 时间t');
ylabel(' 温度T');
axis([0100001]); subplot(2,2,4);
C=T(41,:);% 取矩阵第41 列plot(A,'r'); holdon;% 多条曲线绘制plot(B,'k');
plot(C);
holdoff;
title(' 瞬态响应'); xlabel(' 时间t'); ylabel(' 温度T'); axis([0100001]); legend('A=','B=','C=');。