#CP090-计算物理热传导方程的差分解法

function chafen2 %二维热传导方程的差分解法 lda=1;l=1;s=2;h=0.05;alpha=0.25;tao=alpha*h^2/lda; T=tao*1000;N=l/h;M=s/h;K=T/tao;M1=0.4*M;M2=0.6*M; for i=1:N+1
for j=1:M+1 u(i,j,1)=0; for k=1:K u(i,1,k)=0; u(i,M,k)=0; if M1<=j & j<=M2 u(1,j,k)=1; end end
Q
t dt
t
S
K(x, y, z,t) u ds n
由矢量分析（高斯散度定理）可得：
Q
t dt
t
[K(x, y, z,t)u]dV
V
其中， 是哈密顿算子。
9.1 热传导方程概述
二、产生热量
V 内所有热源产生的热量：
Q
（4）用差分格式计算 ui, j,k 。
9.3 二维热传导方程的差分解法
例 9.2 求热传导方程混合问题：
u

t

u x

u y
x , y , t
u(x, y,)
u(x,,t) u(x,,t) x , t
end surf(u(:,:,k+1)); pause(0.05); end
练习：
求热传导方程混合问题：
u u u

t

x
y
x , y , t
u(x, y,)
u(x,,t) u(x,,t) x , t
k ,,,..., j ,,...,M
ui,,k g (k , ih) ui,N ,k g (k , ih)
k ,,...,M ,i ,,...,N
9.3 二维热传导方程的差分解法
三、差分计算 1、思路 根据初始条件和边界条件计算 出各层的的温度，如图 9-2 所示。 2、差分格式的稳定条件
9.2 一维热传导方程的差分解法
一、一维热传导方程 各向同性介质中无热源的一维热传导方程为：
u t


u x
K , t T （9.3） c
二、初始、边界条件 初始条件：
u(x,) (x)
xl
（9.4）
9.2 一维热传导方程的差分解法
边界条件： 1、 第一类边界条件：
function main %热传导方程的差分解法 lda=1;l=1;h=0.05;alpha=0.5;tao=alpha*
h^2/lda;T=tao*100;N=l/h;M=T/tao; for i=1:N+1
u(i,1)=fai((i-1)*h); end for k=1:M
u(1,k)=g1(k*tao); u(N,k)=g2(k*tao); end for k=1:M for i=2:N
章热传导方程的差分解法
9.1 热传导方程概述
一、传入热量
t 时间内通过 S 横截面积传导的热量为：
Q K(x, y, z,t)tS u n
其中，K (x, y, z,t) 是介质的热传导系数。u 是温度沿 S 面的 n
法向微商，即温度梯度的法向分量。
通过包面 S 传入 V 的热量为：
t dt
t
F(x, y, z,t)dV
V
其中，F(x, y, z,t) 为内部热源密度，表示单位时间单位体积所产生
的热量。
三、消耗热量
V 内消耗热量：
Q
t dt
t
c udV
V t
其中， c 为介质的比热容， 为质量密度。
9.1 热传导方程概述
四、三维非齐次热传导方程 由能量守恒定律，即
h
3、具体步骤
（1）给定 , h,,T , XN,YM ；
（2）计算 h / , N XN / h, M YM / h, K T / ；
（3）计算初始值： ui, j, (ih, jh) ；
计算边界值： u, j,k g (k , jh), uN , j,k g (k , jh) ； ui,,k g (k , ih), ui,N ,k g (k , ih) ；
x
h
h
9.2 一维热传导方程的差分解法
4、 一维热传导方程的差分格式
ui,k ui,k ui,k ui,k ui,k

h
整理得：
ui,k ui,k ( - )ui,k ui,k



h
i
,,...,
N
, k
9.2 一维热传导方程的差分解法
2、差分格式的稳定条件 3、具体步骤

h
（1）给定 ,l, h,,T ；
（2）计算 h / , N l / h, M T / ；
（3）计算初始值： ui, (ih) ； 计算边界值： u,k g (k ), u N ,k g (k ) ；
u(, y,t)

x


y , t
u(, y,t)

x


y M, M y , t
u(, y,t) M y M , t
的数值解，取 N=20，M=20，h=0.05,计算到 k=100 为止。
9.3 二维热传导方程的差分解法
对于节点 (i, j) ，在 k 时刻有：
ui, j,k

t
ui, j,k ui, j,k

ui, j,k

x
ui, j,k
ui, j,k h
ui, j,k


ui y
,j
,k
1、 一阶向前差商
ui,k ui,k ui,k 和 ui,k ui,k ui,k
x
h
t



2、 一阶向后差商
3、 二阶中心差商
ui,k ui,k ui,k
x
h

ui,k ui,k
ui,k x
x ui,k ui,k ui,k
ui, j,k
ui, j,k h
ui, j,k
（9.10）
将差分格式（9.10）代入偏微分方程中得：
ui, j,k ( )ui, j,k (ui, j,k ui, j,k ui, j,k ui, j,k )

h
Q Q Q
可得：
t dt [c u (Ku) F]dV
t
V t
或
c u (Ku) F(x, y, z,t)
（9.1）
t
式（9.1）称为各向同性介质有热源的热传导方程，也叫三维
非齐次热传导方程。
9.1 热传导方程概述
五、三维齐次热传导方程
u(i,k+1)=alpha*u(i+1,k)+(12*alpha)*u(i,k)+alpha*u(i-1,k); end plot([0:h:l],u(:,k+1)); hold on; pause(0.05); end
function u=fai(x) u=4*x*(1-x);
function u=g1(x) u=0;
end end for k=1:K-1
for j=2:M for i=2:N u(i,j,k+1)=(1-4*alpha)*u(i,j,k)+alpha*(u(i+1,j,k)+u(i1,j,k)+u(i,j+1,k)+u(i,j-1,k)); end u(N+1,j,k+1)=u(N,j,k+1); if (2<=j&j<=M1) || (M2<=j&j<=M) u(1,j,k+1)=u(2,j,k+1); else u(1,j,k+1)=1; end
当介质均匀（ c 、 和 K 为常数）、V 内无热源（ F(x, y, z,t) ）时：
c u Ku,
t
其中

x

y

z
上式可表示为：
u t

K
c
u ( x

u y

u z )
（9.2）
上式为三维齐次热传导方程。
i ,,...,N , j ,,...,M , k ,,,...
（9.11）
9.3 二维热传导方程的差分解法
初始、边界条件的差分格式：
ui, j, (ih, jh) i ,,...,N , j ,,...,M
u, j,k g (k , jh) uN, j,k g (k , jh)

,,,...,
M
5、 初始、边界条件的差分格式
ui, (ih) i ,,...,N
u,k g (k ) uN,k g (k )
k ,,...,M
9.2 一维热传导方程的差分解法
四、差分计算 1、思路 根据初始条件和边界条件计算出各层的的温度，如图 9-1 所示。
（4）用差分格式计算 ui,k 。
9.2 一维热传导方程的差分解法
例 9.1 求热传导方程混合问题：
u(utx,)xux( x) u(,t) ,u(,t)
x , t
x t
的数值解，取 N=10，h=0.1,计算到 k=36 为止。
u(, y,t)

x


y , t
（9.10）
方程（9.9）称为二阶抛物线偏微分方程。
9.3 二维热传导方程的差分解法
二、差商格式
设时间步长为 ，空间步长为 h ，二维平
面 xoy 分为 N M 的网格，并使
则有
t k x ih y jh
Nh l Mh s
k ,,,... i ,,,...,N j ,,,...,M
u(,t) g(t) u(l,t) g (t)
2、 第二类边界条件：
t T
u(, t ) x

g (t)

u(l,
t
)
x

g
(t)
t T
（9.5） （9.6）
9.2 ຫໍສະໝຸດ Baidu维热传导方程的差分解法
边界条件： 3、 第三类边界条件：
u(, t ) x
function u=g2(x) u=0;
9.3 二维热传导方程的差分解法
一、二维热传导方程 各向同性介质中无热源的二维热传导方程为：
u t
（xu

xu ），

K
c
, t
T
（9.9）
x l, y s,
初始条件为
u(x, y,) (x, y)
ui,0,k ui,M ,k 0 u0, j,k u1, j,k
i 0,1,2,...,N j 1,2,...,M1 1, M 2 1,...,M 1
uN 1, j,k uN , j,k
j 0,1,...,M
u0, j,k 1
j M1, M1 1,...,M 2

(t )u(, t )

g (t)
u(l , t ) x


(t )u(l , t )

g
(t)
t T （9.7）
其中， (t) 和 (t) ，且不能同时为零。
9.2 一维热传导方程的差分解法
三、差商格式
设空间和时间步长分别为： h 和 ，空间和时间步序号记为 i 和 k 。
9.3 二维热传导方程的差分解法
解：将初始边界混合问题转化为差分格式
ui, j,k 1 (1 4 )ui, j,k (ui1, j,k ui1, j,k ui, j1,k ui, j1,k )
ui, j,0 0
i 0,1,...,N, j 0,1,...,M