有限差分和Matlabpde求解一维稳态传热问题.docx

有限差分和pde函数求解一维定态热传导方程
分别用有限差分方法和pde函数求解一维定态热传导方程，初始条件和边界
条件，热扩散系数α=0.00001，
2
2
TTtx




（1）

求解过程：
1. 用Tylaor展开法推导出FTCS格式的差分方程
首先对T进行泰勒展开得到如下两式子：

23
123123...232!3!23...232!3!nnnnnjjjjjnnnnnjjjjjttTTtxxTTxTTTtttTTTxxx




































上述两个方程变换得：

1122323...23nnnnnnnjjjjjjjTTTTTtTtTotttttt














（2）

223
123...23nnnnnjjjjjTTTxTxTxxxx
















1232422342222...3!4!nnnnnnjjjjjjTTTTxTxTxxxxxx

















2
1122222-nnnjjjTTTToxxx












（3）

将上述式子（2）（3）代入（1）得：
12112222()nnnnnnnjjjjjjjTTTTTTTOtxtxtx







，

（4）
2. 方程的相容性和稳定性讨论：
上述方程截项为：
223324
2334
()...4...23!3!4!nnnnjjjjtTtTxTxTOtxttxx，

由于,0,0limxtotx所以方程有相容性
其经过傅里叶变换后，只有当,tx满足下列条件时，方程具有较好的稳
定性：
2
2

20sin12mkxtx




其中mNkLN为节点个数，L为边界长度由于：

22
sinsin122mkxNxL





所以当20.5tsx时方程具有稳定性
3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。
由式（4）得：

1111122(12)nnnjjjnnnnnjjjjjTTTTTtsTsTsTx








要想知道下一个时刻1nt下的温度，则必须知道在nt时刻下该点以及相邻
两点的温度，如下图所示：

所以需要知道的是在各个时刻x轴两侧的温度，以及初始时刻x轴上各点
温度。下述各题求解过程中，均采用
0
max0.000010(1,2...1)j
Tjj，
0
100nT
，作为边界条件。边界条件是由外界给定，本身对过程求解，精
度无影响
4. 编写M文件求解上述方程
M文件的编写见附录所示。
在matlab中调用heat_conduct(0.00001,15000,500,0.1)，得到如下结果：
rms error is
rms =
0.7959

说明：1.上述结果采用的,xt分别为0.1和500，时间最终给定为15000
当,xt为（0.01,5）时的结果
附录1，传热方程有限差分的 M文件编写
function heat_conduct(alpha,Timemax,dt,dx)
%方程需要输入的参数有，最大时刻数，alpha，dt，dx,实际
Tmax=Timemax/dt+1;%时刻数
JMAX=1/dx+1;%x轴上的差分点个数
J=JMAX-1;MAXEX=J;
T=zeros(Tmax,JMAX);%构造温度矩阵
TE=zeros(Tmax,JMAX);%实际温度矩阵
%设定边界条件
T(:,1)=100;T(:,JMAX)=100;
s=alpha*dt/dx/dx;
%计算各个点的温度
for i=1:(Tmax-1)
for j=2:J
T(i+1,j)=s*T(i,j+1)+(1-2*s)*T(i,j)+s*T(i,j-1);
end
end
%计算实际值
%画图
z=0:dx:1;
t=0:(Tmax-1);
a=3000/dt+1;b=9000/dt+1;
subplot(1,2,1)
plot(z,T(a,:),'g*',z,T(a,:),'g',z,T(b,:),'o',z,T(b,:),'r',z,T(Tmax,:),'b')
xlabel('x方向')
ylabel('温度T')
title('3000,9000,15000时刻的温度分布')
%画整体分布图
for i=1:JMAX
X(:,i)=t;
end;
for j=1:Tmax
Y(j,:)=z;
end
subplot(1,2,2)
mesh(X,Y,T);
view([-1 1 1]);
xlabel('时刻t');
ylabel('方向x');
zlabel('温度T');
title('整体分布')
二、调用MATLAB函数完成上述计算
1. 编写M文件求解上述方程，并用适当的文字对程序做出说明。
Matlab中自带的pdepe函数可直接求解一般偏微分方程。其调用格式为：
T=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);其中T为输出结果，在本题中，其输出格
式为一个二维数组，分别表示一维空间和时间上的温度；x，t，分别为所需输入
的时间和空间数组；
为达到该pdepe函数求解要求，将上述方程转化为如下函数：即@pdefun
所要求的函数的描写格式：

,,,,,,,,mmTTTcxtxxfxtusxtuxtxxx







本题中，在描述@pdefun时，

,,1Tcxtx；10mxm；,,,TTfxtuxx
；,,,0sxtux




@pdeic为初值条件，其描述的标准格式为：

0
0,0Tx

@pdebc为边界条件，本题要求边界温度为100℃，其描述的标注格式为：

,,,,*,,,0TpxtTqxtTfxtux






所以左右边界的,,pxtT与,,qxtT分别表示为：,,100pxtTT边界边界；

,,=0qxtT

边界

上述程序编写见附录2
heat_conduct1(0.1,500,15000)得到如下结果所示：

附录2:Matlab PDE求解上述一维传热方程
%带求解函数
function [c,f,s]=pdefun(x,t,T,dT)