有限差分和Matlabpde求解一维稳态传热问题.(优选)

有限差分和pde 函数求解一维定态热传导方程

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

22

T T t x α∂∂=∂∂ （1） 求解过程：

1. 用Tylaor 展开法推导出FTCS 格式的差分方程

首先对T 进行泰勒展开得到如下两式子：

2

3

1231

2

3

...

232!

3!

2

3

...

232!3!n

n

n

n n j j j

j j n

n

n

n n j j

j

j

j

t

t

T T t x x T

T

x T T T t t t T T T x x x ++∆∆=+∆+

+

+∆∆=+∆+

+

+⎛⎫⎛⎫∂∂∂⎛⎫

⎪ ⎪

⎪∂∂∂⎝⎭

⎝⎭⎝⎭

⎛⎫⎛⎫∂∂∂⎛⎫

⎪

⎪ ⎪∂∂∂⎝⎭

⎝⎭

⎝⎭

上述两个方程变换得：

()11223

23...23n

n

n n n n n j j j j j j j

T T T T T t T t T o t t t t t t ++--⎛⎫⎛⎫

∂∆∂∆∂⎛⎫=

--=+∆ ⎪ ⎪ ⎪∂∆∂∂∆⎝⎭⎝⎭⎝⎭

（2）

223

123...23n

n

n n n j j j

j j T T T x T x T x x x x --⎛⎫⎛⎫

∂∆∂∆∂⎛⎫=

-- ⎪ ⎪ ⎪∂∆∂∂⎝⎭⎝⎭⎝⎭ ()1232422

342222...3!4!n

n

n

n n n j j j j j j T T T T x T x T x x x x x x +-⎛⎫

⎛⎫⎛⎫∂∂∆∂∆∂⎛⎫=--- ⎪ ⎪ ⎪ ⎪∂∆∆∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭

()()2112

22

22-n n n j j j T T T T

o x x x

+--+⎛⎫∂=+∆ ⎪∆∂⎝⎭

（3）

将上述式子（2）（3）代入（1）得：

12112

22

2()n

n

n n

n n n j j j j j j j

T T T T T T T O t x t x t x αα+-+--+⎛⎫∂∂⎛⎫-=-+∆∆ ⎪ ⎪∂∂∆∆⎝⎭⎝⎭，

（4）

2.

方程的相容性和稳定性讨论：

上述方程截项为：

2233242334()...4...23!3!4!n n

n n j j j j

t T t T x T x T O t x t t x x α⎛⎫⎛⎫⎛⎫

⎛⎫⎛⎫∆∂∆∂∆∂∆∂ ⎪∆∆=--++ ⎪ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭

， 由于(),0

,0lim x t o t x ∆∆→∆∆=所以方程有相容性

其经过傅里叶变换后，只有当,t x ∆∆满足下列条件时，方程具有较好的稳定性：

22

20sin 12

m k x

t x α∆∆≤

≤∆ 其中m N

k L

π=

N 为节点个数，L 为边界长度由于：

2

2sin sin 122m k x N x L π∆∆⎛⎫== ⎪⎝⎭

所以当2

0.5t

s x α∆=

≤∆时方程具有稳定性

3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 由式（4）得：

11

1

11

2

2(12)n n n j j j n n n n n j

j

j j j T T T T

T t

sT s T sT x α-++-+-+≈+∆=+-+∆

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

所以需要知道的是在各个时刻x 轴两侧的温度，以及初始时刻x 轴上各点温度。下述各题求解过程中，均采用0max 0.000010(1,2...1)

j T j j α===-，

0100n T =，作为边界条件。边界条件是由外界给定，本身对过程求解，精

度无影响

4. 编写M 文件求解上述方程 M 文件的编写见附录所示。

在matlab 中调用heat_conduct(0.00001,15000,500,0.1)，得到如下结果： rms error is rms = 0.7959

∆∆分别为0.1和500，时间最终给定为15000 说明：1.上述结果采用的,x t

当()

,x t

∆∆为（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