function u=fai(x) u=4*x*(1-x); function u=g1(x) u=0; function u=g2(x) u=0;
9.3 二维热传导方程的差分解法
一、二维热传导方程 各向同性介质中无热源的二维热传导方程为：
u u u K （ ）， , t T （9.9） t x x c x l , y s,
9.3 二维热传导方程的差分解法
例 9.2 求热传导方程混合问题：
u u u x , y , t y t x u ( x, y,) u ( x,, t ) u ( x,, t ) x , t u (, y, t ) y , t x u (, y, t ) y M , M y , t x u (, y, t ) M y M , t
9.2 一维热传导方程的差分解法
例 9.1 求热传导方程混合问题：
u u t x u ( x,) x( x) u (, t ) , u (, t )
x , t x t
的数值解，取 N=10，h=0.1,计算到 k=36 为止。
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 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
将差分格式（9.10）代入偏微分方程中得：
（9.10）
ui , j ,k ( )u i , j ,k (ui , j ,k u i , j ,k u i , j ,k ui , j ,k )


h

i ,,..., N , j ,,...,M , k ,,,...
k ,,,... i ,,,..., N j ,,,...,M
9.3 二维热传导方程的差分解法
对于节点 (i, j ) ，在 k 时刻有：
u i , j ,k u i , j ,k u i , j ,k t u i , j ,k u i , j ,k u i , j ,k u i , j ,k x h u i , j ,k u i , j ,k u i , j ,k u i , j ,k y h
第九章 热传导方程的差分解法
9.1 热传导方程概述
一、传入热量
t 时间内通过 S 横截面积传导的热量为： u Q K ( x, y, z , t )tS n u 其中，K ( x, y, z, t ) 是介质的热传导系数。 是温度沿 S 面的 n
法向微商，即温度梯度的法向分量。 通过包面 S 传入 V 的热量为：
9.1 热传导方程概述
五、三维齐次热传导方程 当介质均匀（ c 、 和 K 为常数） 内无热源（ F ( x, y, z, t ) ）时： 、V
u c Ku, t
上式可表示为：
其中 x y z
u K u u u ( ) t c x y z
4、 一维热传导方程的差分格式
u i ,k u i ,k

整理得：

u i ,k u i ,k u i ,k h
u i ,k u i ,k ( - )u i ,k u i ,k


h

i ,,...,N , k ,,,...,M
上式为三维齐次热传导方程。
（9.2）
9.2 一维热传导方程的差分解法
一、一维热传导方程 各向同性介质中无热源的一维热传导方程为：
u u t x
二、初始、边界条件 初始条件：

u( x,) ( x)
xl
（9.4）

3、具体步骤

h


（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),u N , j ,k g (k , jh) ； ui ,,k g (k , ih),ui , N ,k g (k , ih) ； （4）用差分格式计算 ui , j ,k 。
(ui 1, j ,k ui 1, j ,k ui , j 1,k ui , j 1,k ) i 0,1,..., N , j 0,1,...,M i 0,1,2,..., N j 1,2,...,M 1 1, M 2 1,...,M 1 j 0,1,...,M j M 1 , M 1 1,...,M 2
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 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 end surf(u(:,:,k+1)); pause(0.05); end
的数值解，取 N=20，M=20，h=0.05,计算到 k=100 为止。
9.3 二维热传导方程的差分解法
解：将初始边界混合问题转化为差分格式
ui , j ,k 1 (1 4 )ui , j ,k ui , j , 0 0 u i , 0 , k ui , M , k 0 u0, j ,k u1, j ,k u N 1, j , k u N , j , k u0, j ,k 1
u i ,k x

和
u i ,k t


u i ,k u i ,k


u i ,k u i ,k h
u i ,k x

3、 二阶中心差商
u i ,k ui ,k x


x


h

ui ,k ui ,k u i ,k h
9.2 一维热传导方程的差分解法
Q dt F ( x, y, z, t )dV
t V
t
其中， F ( x, y, z, t ) 为内部热源密度，表示单位时间单位体积所产生 的热量。 三、消耗热量 V 内消耗热量：
Q
t
t
其中， c 为介质的比热容， 为质量密度。
u dt c dV V t
（9.11）
9.3 二维热传导方程的差分解法
初始、边界条件的差分格式：
u i , j , (ih, jh) u , j ,k g (k , jh) u N , j ,k g (k , jh) u i , ,k g (k , ih) u i , N ,k g (k , ih)
Q dt
t
t
S
u K ( x, y, z , t ) ds n
由矢量分析（高斯散度定理）可得：
Q dt [K ( x, y, z, t )u]dV
t V
t
其中， 是哈密顿算子。
9.1 热传导方程概述
二、产生热量 V 内所有热源产生的热量：
9.1 热传导方程概述
四、三维非齐次热传导方程 由能量守恒定律，即
Q Q Q
可得：
或
u t dt V [c t ( Ku) F ]dV u c ( Ku ) F ( x, y, z, t ) t
t
（9.1）
式（9.1）称为各向同性介质有热源的热传导方程，也叫三维 非齐次热传导方程。
t T
（9.6）
9.2 一维热传导方程的差分解法
边界条件： 3、 第三类边界条件：
u (, t ) x (t )u (, t ) g (t ) u (l , t ) (t )u (l , t ) g (t ) x
其中， (t )
i ,,..., N , j ,,...,M k ,,,..., j ,,...,M k ,,...,M , i ,,..., N