有限差分法解热传导问题

偏微分方程工具箱
Step 1 “Draw模式”绘制平面有界区域 ，通过公式把Matlab系统提供的实体模 型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域. Step 2 “Boundary模式”定义边界，声明不同边界段的边界条件. Step 3 “PDE模式”定义偏微分方程，确定方程类型和方程系数c,a,f,d，根据具 体情况，还可以在不同子区域声明不同系数. Step 4 “Mesh模式”网格化区域 ，可以控制自动生成网格的参数，对生成的网 格进行多次细化，使网格分割更细更合理. Step 5 “Solve模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自 适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界 条件后可以求出给定时刻t的解；对于特征值问题，可以求出给定区间上的特征值. 求解完成后，可以返回到Step 4，对网格进一步细化，进行再次求解. Step 6 “View模式”计算结果的可视化，可以通过设置系统提供的对话框，显 示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型 问题的解还可以进行动画演示.
x
x
x
Φ上 Φ下 Φ左＋Φ右 Φv 0
tm1,n tm ,n tm1,n tm ,n t m ,n1 tm ,n tm ,n1 tm ,n y y x x x x y y xy 0 Φ
x y
tm1,n tm1,n tm,n1 tm,n1 4tm,n 4tm,n tm1,n tm1,n tm,n1 tm,n1
Gauss－Seidel迭代
(k 1) 1 (k ) (k ) x ( a x a x b1 ) 12 2 1n n 1 a11 1 (k 1) ( k 1) (k ) (k ) x 2 (a21x 1 a23x 3 a1n x n b2 ) a22 (k 1) 1 (k 1) ( k 1) x n (an1x 1 an n 1x n 1 bn ) ann
NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}]+x,u[x,0]=x*(1x),u[0,t]=0,u[1,t]=0,u,{x,0,1},{t,0,0.3}]
NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=Sin[x]*Sin[x],Derivative[1, 0][u][0,t]=0,Derivative[1,0][u][Pi,t]=0,u,{x,0,Pi},{t,0,1}]
NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=-x*(1x),Derivative[1,0][u][0,t]=0,Derivative[1,0][u][1,t]=3u[1,t],u,{x,0,1},{t,0,0.3}]
建立控制方程及定解条件 设立温度场的迭代初值
xi
( k 1) n 1 i 1 ( k 1) (k ) ( aij x j aij x j bi ) aii j 1 j i 1
Gauss－Seidel迭代
x( k 1) D1 ( Lx( k 1) Ux( k ) b) ( I D1L) x( k 1) D1Ux( k ) D1b x( k 1) ( I D1L)1 D1Ux( k ) ( I D1L)1 D1b x( k 1) ( D L)1Ux( k ) ( D L)1 b
同样可得：
2t tm,n1 2tm,n tm,n1 2 o ( y ) 2 2 y m,n y
对于二维稳态导热问题，在直角坐标中，其导热 微分方程为：
v 2t 2t 2 0 2 x y
其节点方程为：
ti 1, j 2ti , j ti 1, j
qw
y x
x y
4tm,n 2tm1,n 2x m ,n qw tm,n1 tm,n1 Φ x 2


17
(2) 外部角点
qw
y tm1,n tm,n y qw 2 x 2 x x tm,n1 tm,n qw 2 2 y x y Φm ,n 0 2 2
x
2

ti , j 1 2ti , j ti , j 1
y 2
v ,i , j 0
热平衡法
基本思想：对每个有限大小的控制容积应用能量守恒，从 而获得温度场的代数方程组，它从基本物理现象和基本定 律出发，不必事先建立控制方程，依据能量守恒和Fourier
导热定律即可。
能量守恒：流入控制体的总热流量＋控制体内热源生成热 ＝ 流出控制体的总热流量＋控制体内能的增量
内部节点：Φm1,n Φm1,n Φm ,n1 Φm ,n1 0
(m,n+1)
y
Φ上 Φ下 Φ左＋Φ右 0
(m-1,n)
(m, n)
(m+1,n)
y
(m,n-1)
y o
x y
y x
2tm,n
x 2 m ,n tm1,n tm,n1 qw Φ 2
2x
18
(3) 内部角点
qw
tm1,n tm ,n y tm1,n tm ,n y y qw x x 2 2 tm ,n1 tm ,n x tm ,n1 tm ,n x x qw y y 2 2 3xy Φm ,n 0 4
B (D L )1U
, f (D L )1b
Gauss－Seidel迭代
200℃
t1
t2
t3 100℃
100℃
t16
t17
t18
Tf=0 ℃ k=1W/(m*K) h=10W/(m2* ℃) (qw=h*(Tw-Tf))
Ax=b
b=[300,200,200,300,100,0,0,100,100,0,0,100,100,0,0,100,0,0]’ x=[t1,t2,t3,...t18]’
确定节点（区域离散化） 建立节点物理量的代数方程
求解代数方程
改进初场
是否收敛 是 解的分析
来自百度文库
否
N
(m,n)
n
y
y M
x
x
m
泰勒级数展开法：
若取上面式右边的前三项，并将两式相加移项整理即 二阶导数的中心差分：
2t x 2
m ,n
tm1,n 2tm,n tm1,n 2 o ( x ) 2 x
x y
t m ,n
x
y
1 (2tm1,n 2tm ,n1 tm ,n1 tm1,n 6 3x 2 2x 2 qw ) 2
19
写出所有内节点和边界节点的温度差分方程
n个未知节点温度，n个代数方程式：
t1 a11t1 a12t 2 ...... a1n t n b1 t 2 a21t1 a22t 2 ...... a2 n t n b2 .......... .......... .......... .......... .... t n an1t1 an 2t 2 ...... annt n bn
傅里叶定律：在导热现象中，单位时间内通过给定截 面的热量，正比例于垂直于该界面方向上的温度变化 率和截面面积，而热量传递的方向则与温度升高的方 向相反。
导热微分方程：
定解条件：使微分方程获得适合某一问题的的解的附 加条件。
边界条件：
NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=x*(1x),u[0,t]=0,u[1,t]=0,u,{x,0,1},{t,0,0.3}]
x 2

x 2
0 Φ Φ

1.边界节点离散方程的建立：
qw
(1) 平直边界上的节点
tm1,n tm ,n y yqw x x tm ,n1 tm ,n x tm ,n1 tm ,n 2 y 2 y x Φm ,n y 0 2