偏微分一维热传导问题

偏微分大作业

一维热传导方程问题

——运用隐式格式求解数值解

目录

问题描述 ..................................................................................................... 3 1 解析解——分离变量法 ......................................................................... 3 2 数值解——隐式格式 ............................................................................. 5 3 证明隐式格式的相容性与稳定性 ......................................................... 5 4 数值解——分析与Matlab实现 ........................................................... 6 5 数值解与解析解的比较 ......................................................................... 9 6 随时间变化的细杆上的温度分布情况 ............................................... 11 7稳定后细杆上的温度分布情况 ............................................................ 12 参考文献 ................................................................................................... 13 附录 ........................................................................................................... 14

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

有限长杆的一维热传导问题

问题描述

一根单位长度的细杆放入100℃的沸水中，当细杆的温度达到100℃时取出。假设细杆四周绝热；在时间t=0时，细杆两端浸入0℃的冰水中。一维热传导方

2uauxx0，现在令a21，从而可知本题：utuxx0。现在程：t要求细杆温度分布：u(x,t)。

1 解析解——分离变量法

热传导偏微分方程：

其中，

utuxx0 (1)

首先令：

u(x,t)X(x)T(t) (2)

将(2)式带入(1)式得： 于是可得：

可以得到两个微分方程：

先求解空间项：

当0时， X(x)AexBex

由于u(0,t)u(1,t)0,t.

可知:由于解的收敛性，B0

则此时是平庸解。

当0时， X(x)ABx

0时， X(x)AcoskxBsinkx，其中k。

n1,2,3...

则此时是平庸解。

当 所以，X(x)Bnsin(nx)，

因为n22

所以，T(t)Cnen22t，

n1,2,3...

则，

文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

0)初始条件：u(x，最终，

(x)

200nn22tu(x,t)(1(1))esin(nx)， n1,2,3...

n1n2 数值解——隐式格式

目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式1格式。

利用u(x,t)，关于t进行向前差商：1k1k1Uk2UUj1jj11kUkUjjt ；关于x进行二阶中心差：

(x)2；

代入偏微分方程可以得到隐式差分格式：

1kUkUjjt1k1k1Uk2UUj1jj1(x)2 （1）

3 证明隐式格式的相容性与稳定性

（1）相容性

代入隐式格式得：

2U12UU22t(t)=+(t)(x) (2) 22t2tx 将（2）与原微分方程相减，得到截断误差：

所以此隐式格式与原微分方程相容。

（2）稳定性

令网格比为rt2x，则可以将（1）式改写得到：

1k1k1krUk(12r)UrUU 3） j1jj1j （

首先令：

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

1k1I(j-1）UkUej11k1IjUkUejUUekjkIj （4）

1k1I(j+1）UkUej+1将（4）代入（3）式，根据欧拉公式化简得： U故得放大因子是：

所以根据Fourier方法，隐式格式恒稳定。

k1 （1+2r2rcos)Uk （5）

4 数值解——分析与Matlab实现

（1） 边值与初值离散化

将边值与初值离散化，与式（3）联立得差分线性方程组：

1k1k1krUk(12r)UrUUj1jj1j ， j(0,1,2,,M-1)

ULk0， k(0,1,2,B的形式：

,N)

再将方程组改写成AU本题的边界条件均为零。所以可以将上式改写。

（2） Matlab的实现

➢ 杆长1米，时间2秒。

设计空间步长h=0.1和时间步长t=0.01，网格比是rth2。

从而得到划分的空间网格点数是M1+1，时间网格点数是M2+1。先设初始的温度矩阵U(M2+1,M1+1)。再将边界条件和初始条件编写到表示温度分布的矩阵中。具体代码可见最后附录。 ➢ 编写矩阵A

核心代码：对角线：A(i,i) = 1+2r

对角线的右方和下方：A(i,i+1) = -r; A(i+1,i) = -r; ➢ 下面就要运用A*U(k1,j)U(k,j)进行迭代。

当k=1时，A*U(2,j)=U(1,j) 当k=2时，A*U(3,j)=U(2,j) 当k=3时，A*U(4,j)=U(3,j)

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

以此迭代下去直到k=M2。就可以得到整个温度随时间和空间的分布矩阵U。

➢ 数值解画图，如图1(a)和图1(b)所示。

图1(a) 数值解的温度分布图

现在将着色平稳过渡。

图1(b) 着色平稳过渡的数值解的温度分布图

5 数值解与解析解的比较

➢

首先，我们需要将解析解离散化，解析解中有一项e里的工作空间运算。

n22t，当n越来越大

时，会快速趋于0，故我们可以取n=8000。现在来证明可行性，在matlab将解析解的温度分布画出来，数值解画图

2，如图2所示。

图2 解析解的温度分布图

将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a)上可以看出空间上的误差，在边界处误差比较大。

图3（a） 数值解与解析解空间误差

我们从图3(a)上可以看出时间的误差，在时间的最开始，处误差最大，然后又有一个小的波动，最后就误差渐渐变小，最后趋于0。

图3（b） 数值解与解析解时间误差

6 随时间变化的细杆上的温度分布情况

从数值解的温度分布三维图，如图4(a)和图4(b)可以看出随着时间的增加，细杆温度下降最后趋于0℃。

从物理角度来说：细杆的温度会不断地向两端扩散，热量会慢慢散失，最终随着时间的增加，细杆的温度会趋于0℃。

图4(a) 细杆温度随时间的变化图

现取细杆中心处一点，观看它随时间的温度变化情况。

图4(b) 细杆中央（x=0.5）温度随时间的变化图

7稳定后细杆上的温度分布情况

从图像上可以看出，最后稳定的情况下，细杆的温度是0℃。

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

参考文献

[1] 冯立伟．热传导方程几种差分格式的MATLAB数值解法的比较[J]．沈阳化工大学，辽宁沈阳．2011(6)．

[2] 一维热传导方程数值解法及Matlab实现[EB/OL]．2014-11-20 -945.html

附录

代码：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %此程序用于解决一维热传导方程:ut-a^2uxx = 0 %

%边界条件：u(0,t) = u(L,t) = 0 % %初始条件：u(x,0) = 100, x!=0和L % % u(0,0) = 0 % % u(L,0) = 0 % %其中,a^2 = 1, L = 1 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc; clear all; %区域及划分网格

L = 1; %单位长度的细杆• T = 2; %时间

h = 0.1; %%%% 空间的划分 %%%% t = 0.01; %%%% 时间的划分%%%% r = t/(h*h); %网格比 %设计步长 M1 = L/h; M2 = T/t;

%构造边界条件 %构造的矩阵：U(时间，空间） U = zeros(M2+1,M1+1); %编程包含边值，如U(k,1)=u(0,t)

for k = 1:M2+1 %时间划分了M2份，有M2+1个节点 U(k,1) = 0; %两个边界处温度恒为零 U(k,M1+1) = 0; end; %构造初始条件

for j = 2:M1 %位置划分了M1份，有M1+1个节点 U(1,j) = 100; end;

U(1,1) = 0; U(1,M1+1) = 0;

%差分格式的矩阵形式 A*U(k+1,j)=U(k,j)

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

%构造矩阵A

A = zeros(M1-1); for i = 1:M1-1 A(i,i) = 1+2*r; end;

for i = 1:M1-2 A(i,i+1) = -r; A(i+1,i) = -r; end;

%构造AU=B中的B %本题边值的特殊，矩阵B大大简化了 B = zeros(M1-1,1); for k = 1:M2 j = 2:M1;

B(j-1,1) = U(k,j); x = A\\B; for j = 2:M1

U(k+1,j) = x(j-1); %k+1时刻的不同位置的温度分布 end; end; %作图

x = 0:h:1; y = 0:t:2;

[xx,yy]=meshgrid(x,y); figure(1); surf(xx,yy,U); shading flat

title('一维热传导方程--数值解--温度分布图'); xlabel('位置x'); ylabel('时间t'); zlabel('温度T'); figure(2) s = 0; for i= 1:8000

s = s+(200*(1-(-1)^i))/(i*pi)*sin(i*pi*xx).*exp(-i^2*pi^2*yy); end;

surf(xx,yy,s);

title('一维热传导方程--解析解--温度分布图'); xlabel('位置x'); ylabel('时间t'); zlabel('温度T'); figure(3) x = 0:h:1; y = 0:t:2;

[xx,yy] = meshgrid(x,y);

1文档来源为:从网络收集整理.word版本可编辑.

文档收集于互联网，已重新整理排版.word版本可编辑.欢迎下载支持.

dd = U-s; surf(xx,yy,dd);

title('一维热传导方程--误差--温度分布图'); xlabel('位置x'); ylabel('时间t');

zlabel('误差（数值解减解析解）'); figure(4)

z = zeros(M2+1,1); if mod((M1+1),2)~=0 i = 1:M2+1;

z(i,1) = U(i,M1/2); else

i = 1:M2+1;

z(i,1) = U(i,(M1+1)/2); end; plot(z);

title('温度随时间增加的趋势图'); xlabel('时间t'); ylabel('温度T');

1文档来源为:从网络收集整理.word版本可编辑.