有限差分法解热传导问题

傅里叶定律：在导热现象中，单位时间内ቤተ መጻሕፍቲ ባይዱ过给定截
面的热量，正比例于垂直于该界面方向上的温度变化
率和截面面积，而热量传递的方向则与温度升高的方 向相反。
导热微分方程：
定解条件：使微分方程获得适合某一问题的的解的附 加条件。
边界条件：
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}]
x2
y 2

热平衡法
基本思想：对每个有限大小的控制容积应用能量守恒，从 而获得温度场的代数方程组，它从基本物理现象和基本定 律出发，不必事先建立控制方程，依据能量守恒和Fourier 导热定律即可。
能量守恒：流入控制体的总热流量＋控制体内热源生成热 ＝ 流出控制体的总热流量＋控制体内能的增量
内部节点：Φm1,n Φm1,n Φm,n1 Φm,n1 0
x


y 2
tm1,n tm,n x

y 2 qw


x tm,n1 tm,n
y



x 2
tm,n1 tm,n y

x 2
qw


Φ m,n
3xy 4

0
x y
1 tm,n 6 (2tm1,n 2tm,n1 tm,n1 tm1,n
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}]
建立控制方程及定解条件
确定节点（区域离散化）
知识回顾 Knowledge Review
Gauss－Seidel迭代

x1(k

x

(k 2
1) 1)

a111(a12x
(k 2
a212(a21x
(k 1
)
a1n
x
(k n
)
b1
)
1)
a23x
(k 3
)
a1n
x
(k n
)
b2
)


x

(k n
1)

an1n (an1x 1(k 1) an
设立温度场的迭代初值
建立节点物理量的代数方程
求解代数方程
改进初场
是否收敛 否
是 解的分析
(m,n) N
n
y
y
x x
M m
泰勒级数展开法：
若取上面式右边的前三项，并将两式相加移项整理即 二阶导数的中心差分：
2t
tm1,n 2tm,n tm1,n o(x2 )
x2 m,n
偏微分方程工具箱
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
y
y
Φxy 0
x y
tm1,n

tm1,n

tm,n1

tm,n1

4tm,n

x2

Φ

0
4tm,n

tm1,n
tm1,n
tm,n1
tm,n1

x2

Φ
1.边界节点离散方程的建立：
qw
(1) 平直边界上的节点
y tm1,n tm,n
x 2
同样可得：
2t tm,n1 2tm,n tm,n1 o(y2 )
y2 m,n
y 2
对于二维稳态导热问题，在直角坐标中，其导热 微分方程为：
2t x 2

2t y 2
v
0
其节点方程为：
ti1, j 2ti, j ti1, j ti, j1 2ti, j ti, j1 v,i, j 0
x b (k1)
n1 n1
n
)
x (k1) i

1 aii
(
i 1 j 1
aij
x
j
(k
1)

n
aij x j (k )
j i 1
bi )
Gauss－Seidel迭代
x(k1) D1(Lx (k1) Ux (k ) b) (I D1L)x(k1) D1Ux (k ) D1b x(k1) (I D1L)1 D1Ux (k ) (I D1L)1 D1b x(k1) (D L)1Ux (k) (D L)1b
(m,n+1)
y
Φ上 Φ下 Φ左＋Φ右 0
y
y o
(m-1,n)
x
x
(m, n) (m,n-1)
x
(m+1,n)
Φ上 Φ下 Φ左＋Φ右 Φv 0
y tm1,n tm,n y tm1,n tm,n x tm,n1 tm,n x tm,n1 tm,n

3x2
2


2x2

qw
)
19
写出所有内节点和边界节点的温度差分方程 n个未知节点温度，n个代数方程式：
t1 a11t1 a12t2 ...... a1ntn b1 t2 a21t1 a22t2 ...... a2ntn b2 ........................................ tn an1t1 an2t2 ...... anntn bn
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]’
x
yqw
x tm,n1 tm,n x tm,n1 tm,n
2 y
2 y
qw

Φ m,n
x 2
y

0
y x
x y
4tm,n

2tm1,n

2x
qw
tm,n1
tm,n1
Φm,n
x2

17
(2) 外部角点
qw
y x

y tm1,n tm,n 2 x

y 2
qw

x 2
qw


x 2
tm,n1 y
tm,n

Φ m,n
x 2

y 2

0
x y
2tm,n

tm1,n

tm,n1

2x

qw
Φm,n
x2
2
18
(3) 内部角点
y x
qw
y tm1,n tm,n
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}]