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

计算传热学程序报告
题目：一维非稳态导热问题的数值解姓名：
学号：
学院：能源与动力工程学院
专业：工程热物理
日期：2014年5月25日
一维非稳态导热问题数值解
求解下列热传导问题：
⎪⎪⎩
⎪⎪⎨⎧=====≤≤=∂∂-∂∂1,10),(,1),0(0)0,()0(01T 22ααL t L T t T x T L x t T x 1.方程离散化
对方程进行控制体积分得到：
dxdt t T dxdt x T t
t t e
w t
t t e w ⎰⎰⎰⎰∆+∆+∂∂=∂∂α122 ⎰⎰-=∂∂-∂∂∆+∆+e
w t t t w e t
t t dx T T dt x T x T )(1])()([α
非稳态项：选取T 随x 阶梯式变化，有
x T T dx T T t p t t p e w t t t ∆-=-∆+∆+⎰
)()( 扩散项：选取一阶导数随时间做显示变化，有 t x
T x T dt x T x T t
w t
e w e t
t t ∆∂∂-∂∂=∂∂-∂∂⎰∆+])()[(])()[( 进一步取T 随x 呈分段线性变化，有
e P E e x T T x T )()(δ-=∂∂ ， w
W P w x T T x T )()(δ-=∂∂ 整理可以得到总的离散方程为：
2
21x T T T t T T t W t P t E t P t t E ∆+-=∆-∆+α 2.计算空间和时间步长
取空间步长为：
h=L/N
网格Fourier 数为：
220x
t x t
F ∆∆=∆∆=α（小于0.5时稳定）
时间步长为： α20h F n =
3.建立温度矩阵与边界条件
T=ones(N+1,M+1)
T(:,1)=Ti (初始条件温度都为0）
T(1,:)=To (边界条件x=0处温度为1）
T(N+1,:)=Te (边界条件x=L 处温度为0）
4.差分法求解温度
由离散方程可得到：
t P t W t P t E t t E T T T T F T -+-=∆+)2(0
转化为相应的温度矩阵形式：
),()],(2),1(),1([)1,(0k m T k m T k m T k m T F k m T +*--++*=+
5.输入界面
考虑到方程的变量，采用inputdlg 函数设置5个输入变量，对这5个变量设置了默认值，如图1所示。

在计算中可以改变不同的数值，得到不同的结果，特别注意稳定条件的临界值是0.5。

根据设置的默认值，得到的计算结果如图2所示。

图1 matlab 变量输入界面
图2 默认值的计算结果
6.结果分析
根据上面的分析，给出了程序的输入界面，以及默认值状态下的数值解。

可以通过改变不同的输入值，得到需要的分析结果，总结出了下面4点结论：
（1）取F0=0.48，得到一维非稳态导热结果如下图所示
图2 F0=0.48时一维非稳态导热
从图中可以看出，对于长度L=1的细杆，初始时刻t=0时温度为0，边界条件x=0时，T=1,边界条件x=1时，T=0。

随着时间的增加，温度从x=0通过导热的形式传递到x=1，不同时刻不同位置杆的温度都不同，并且随着时间的增加，杆的温度也逐渐增加。

（2）取F0=0.48，可以得到不同位置的温度响应曲线，如下图所示
图3 F0=0.48时不同x位置处的温度响应
图中红色曲线代表x=0.1位置的温度瞬态响应，黑色曲线代表x=0.2位置的温度瞬态响应，蓝色曲线代表x=0.4位置的温度瞬态响应。

从图中可以看出，随着x的增加，曲线与x轴的交点值越大，温度开始传递到该位置的所需的时间越长。

随着x的增加，温度响应曲线的变化速率越慢，最终的达到的温度也越低。

（3）取F0=0.25，得到不同位置的温度响应曲线如下图所示
图4 F0=0.25时不同x位置处的温度响应
图中三条曲线分别是x=0.1，x=0.2，x=0.4位置的温度瞬态响应。

与图3的F0=0.48进行对比，两种情况下的F0值不同，F0值越大表明热扩散系数 的值越大。

从图中可以看出热扩散系数对于导热的影响，F0=0.25时，与F0=0.48相比
较，各位置开始响应时所需的时间较长，而且各位置响应曲线的变化速率较小，最终的达到的温度也较低，说明了热扩散系数越小，热传导越慢，传递效率越低。

（4）取F0=0.51，得到非稳定的数值解如图所示
图5 F0=0.51时一维非稳态导热
图6 F0=0.51时不同x位置处的温度响应
从图中可以看出，对于显示格式的离散方程，并不是所有的F0值都能得到有意义的解，必须要求F0<0.5时才能得到稳定的数值解，当F0>0.5时，会出现物理上不真实的解。

附件：（matlab程序）
function heat_conduction() %一维齐次热传导方程
%设置输入界面
options={'空间杆长L','空间点数N' ,'时间点数M','扩散系数a','稳定条件的值Fo(临界值0.5)',}; topic='一维非稳态导热';%标题栏显示
lines=1;%输入行为1行
def={'1','100','1000','1','0.48'};%默认值输入
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的值必须小于0.5，小于0.5波动
%计算空间步长与时间步长
h=L/N;%空间步长
x1=0:h:L;
x=x1';
n=Fo*h^2/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
for k=1:M
m=2;
while m<=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=0.1');%显示函数名称title('x=0.1瞬态响应');
xlabel('时间t');
ylabel('温度T');
axis([0 1000 0 1]);%坐标轴数值范围subplot(2,2,3);
B=T(21,:);%取矩阵第21列
plot(B,'k');
legend('B=0.2');
title('x=0.2瞬态响应');
xlabel('时间t');
ylabel('温度T');
axis([0 1000 0 1]);
subplot(2,2,4);
C=T(41,:);%取矩阵第41列
plot(A,'r');
hold on;%多条曲线绘制
plot(B,'k');
plot(C);
hold off;
title('瞬态响应');
xlabel('时间t');
ylabel('温度T');
axis([0 1000 0 1]);
legend('A=0.1','B=0.2','C=0.4');。