第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
一维非稳态导热问题的数值解.pdf
计算传热学程序报告
题目：一维非稳态导热问题的数值解
姓名：
学号：
学院：能源与动力工程学院
专业：工程热物理
日期： 2014 年 5 月 25 日
一维非稳态导热问题数值解
求解下列热传导问题：
2
T
2
x
T
1
T
t
0)0,(
tLT
),
1
xT
,1),0(
t
L
,1
0
0(
(
x
L
)
0
1.方程离散化
对方程进行控制体积分得到：
t
t
t
e
w
2
T
2
x
d x d t
1
t
t
t
e
w
T
t
d x d t
t
([
t
t
T
x
)
e
(
T
x
)
w
]
dt
1
e
w
t
T
(
t
t
T
)
dx
非稳态项：选取 T 随 x 阶梯式变化，有
e
w
t
T
(
t
t
T
)
dx
t
t
T
p
(
t
T
p
)
x
扩散项：选取一阶导数随时间做显示变化，有
t
t
t
[(
T
x
)
e
(
T
x
)
]
dt
w
[(
t
e
)
(
T
x
T
x
t
])
w
t
进一步取 T 随 x 呈分段线性变化，有
(
T
x
)
e
T
E
(
P
T
)
e
x
，
(
T
x
)
w
T
P
(
T
W
)
w
x
整理可以得到总的离散方程为：
1
t
T
E
t
t
T
P
t
T
E
t
t
T
W
t
T
2
P
2
x
2.计算空间和时间步长
取空间步长为：
网格 Fourier 数为：
h=L/N
F
0
t
2
x
t
2
x
（小于 0.5 时稳定）
时间步长为：
2
h
Fn
0
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
T
t
E
t
TF
(0
E
t
T
2
P
t
T
W
)
t
T
P
转化为相应的温度矩阵形式：
kmT
(
,
)1
F
0
[
mT
(
,1
k
)
kmT
,1
(
2)
kmT
(
,
)]
kmT
(
,
)
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');
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
