第5讲 离散方程的求解

n A Bi i 1 Aij k ij j j k Aii j 1 Aii j i 1 Aii
k+1 为新迭代级 对于高斯-赛德尔迭代方法
k 为上一步迭代级
收敛性判据
i
k 1
n A Bi i 1 Aij k 1 ij j j k Aii j 1 Aii j i 1 Aii
• 内部节点的差分方程
aiTi 1 biTi ciTi 1 di
i 2, 3, 4,......,( N - 1)
计算机实现：算例
• 其中，
ai x w ci x e
bi ai ci (S p )i x
为了避免数值振荡解
0 Rcell 2
30
带有中心差分的数值振荡
31
迎风格式III
例子如下：
y 20 2
x 0.1
i 1
然而
i
i 1
如果 u 0
i i 1
x u (
使用向后差分
u 0
aPTP a E TE aW TW b
3
控制方程向代数方程的转化III
对于节点2和节点3，系数为 ： kA kA aE aW aP aE aW x x 对于节点1，系数为：
b qx
kA aE x
aW 0
aP aE
2kA x
2kA b qx TA x b qx 2kA TB x
b1T2 c1T3 d1
计算机实现：算例
• 当i＝3，4，…，N-2时，
w ai 1 x e ci 1 x
di 1 (Sc )i 1 x
2i 1i e i 1 i
bi 1 ai 1 ci 1 (S p )i x
max i( k 1) i( k )
9
迭代法 – 雅可比方法
利用雅可比迭代方法求解方程组的演示
10
迭代方法– 高斯-赛德尔方法
利用高斯-赛德尔迭代方法求解方程组示例
11
二维热传导– 高斯-赛德尔方法
利用高斯-赛德尔迭代方法求解方程组
12
差分方程的求解
aPTP aW TW aETE bP
迎风格式I
包含对流项和扩散项的控制方程：
2 u 0 2 x x
u(
i 1 i 1
2x
) (
i 1 2i i 1
x
2
)0
Rcell
网格雷诺数
ux
ux (0.5i 1 0.5i 1 ) i 1 2i i 1 0
aN 2TN 2 bN 2TN 1 d N 2
计算机实现：算例
• 最后得到由（N－2）个方程构成的方程组为
b1T2 c1T3 d1 ai 1Ti 1 bi 1Ti ci 1Ti 1 di 1
(i 3, 4, 5,......,N - 2)
aN 2TN 2 bN 2TN 1 d N 2
• 将前面得到的差分方程改写为，
aW TW (aP )TP aETE bP

(1)
或者简单的写成矩阵的形式，
[ A][T ] [ D]
其中，
[T ]T [T1, T2 , T3 ,......,Tn1, Tn ]T [D]T [D1, D2 , D3 ,......,Dn1, Dn ]T
差分方程的求解
b1 a 2 [ A]
c1 b2 a3 c2 b3 . c3 . . . . . . . . an an 1 bn 1
cn 1 bn
差分方程的求解
• 与方程（1）对比，知，
1
反演方法效率不高
高斯消元法并非最好的方法
对于方程个数较多的方程组，可以运用迭代法求解
8
矩阵求解 II
每个节点未知变量 的一般形式可以写成：
A
j 1 ij
i 1
j
Aiii
j i 1
A
ij
n
jห้องสมุดไป่ตู้
Bi
重新排列上述方程，用雅可比方法迭代：
i
k 1
2i (1 0.5Rcell )i 1 (1 0.5Rcell )i 1
29
迎风格式 II
2i (1 0.5Rcell )i 1 (1 0.5Rcell )i 1
aPP anbnb
如果 Rcell 3

anb
aP
1 收敛条件
1 0.5Rcell 1 0.5Rcell 3 2
x/2
对于节点4，系数为：
aE 0
x/2 x
aW
kA x
2kA aP aW x 控制体体积
x
1 TA W w
2 P e
3
4 TB
E
4
x
控制方程向代数方程的转化 IV
代入
k 5W/m2
2 q 500kW/m3 x 0.005 m A 1m
控制体每一节点的系数值
di (Sc )i x
2i 1i w i 1 i 2i 1i e i 1 i
计算机实现：算例
• 注意：采用内节点法划分网格时，近边界节点 与其它内部节点不尽相同，所以必须单独考虑。

当i＝2时
(x)w=½ x w= W= 1 所以，当i＝2时
第五讲 离散方程的求解
屠基元 教授 清华大学 墨尔本皇家理工大学
控制方程向代数方程的转化I
均匀发热无限大平板的稳定导热问题。
控制体体积 x/2
1
x
2 w P
x
3 e 4
x/2
TA=100 oC
A
y
TB=400 B
oC
TA W
TB
E
x
q =500 W/m3
大平板区域的有限体积离散
L
x
z
2
控制方程向代数方程的转化 II
5
=
二维离散化I
一维网格 T Ti 1 Ti x x
I=1
I=1 J = NJ J = NJ
i–1
i–1
i
i
i+1
i+1
I=N
I=N
2T Ti 1 2Ti Ti 1 2 x x 2
x
二维网格
2 2 2
j+1 y
j+1 j i 1, j j j–1 j–1
对于恒定的传热系数和热流量，方程变为：
T T k k qx 0 x e x w
引入线性近似梯度：
T T T T kA E P kA P w qx 0 x x
将上述方程进行重新整理，得到代数方程：
计算机实现：算例
• 求解下面的一维稳态导热问题：
d dT S 0 dx dx
T x0 T0 , T x2 1
aT 3 3 bT 0 x 1 1 x 2
S S (T , x)
计算机实现：算例
• 求解区域的离散化：
– 内节点法：先划分控制容积，在确定节点 – 均匀网格：x=x – 将整个求解区域划分为（N-2）个控制容积，N个节 点（包括2个边界节点）
差分方程的求解
• TDMA法Fortran源程序
SUBROUTINE TDMA(A,B,C,D,N) REAL A(N),B(N),C(N),D(N) DO 5 I=1,N-1 F=A(I+1)/B(I) B(I+1)=B(I+1)-F*C(I) 5 D(I+1)=D(I+1)-F*D(I) D(N)=D(N)/B(N) DO 10 I=N-1,1,-1 10 D(I)=(D(I)-C(I)*D(I+1))/B(I) RETURN END
计算机实现：算例
希望大家用计算机完成上面的计算，并与下面 的分析解结果比较：
x 4 T0 b(1 T04 ) 4 ab T 2 x 1 b(1 T04 ) ab
0 x 1 0 x 1
特别提示
计算机实现的基础地位 关键：掌握循环变量的使用 基础：对算法清晰透彻的把握 保障：细心细心再细心
2i 1i w i 1 i
ai 1Ti 1 bi 1Ti ci 1Ti 1 di 1
计算机实现：算例
• 同样，当i＝N-1时
(x)e=½ x e= E= N 所以，当i＝ N-1时
(x)w (x)e (x)e
1
2
3
N-1
N
计算机实现：算例
x
Ti 1, j 2Ti , j T T T T x 2 y 2 x 2
T T Ti 1, j 2Ti , j Ti 1, j Ti , j 1 2Ti , j Ti , j 1 2 2 2 x y x y 2
2 2

Ti , j y j) 2Ti , j Ti , j 1 (i, 1
(i, j)
y 2
y
y
x
x
J=1 I=1 J=1 I=1
i–1
i–1
i
i
i+1
i+1
I = NI
I = NI
6
二维离散化 II
因为
2T 0
代数方程式
x y
Ti 1, j Ti 1, j Ti, j 1 Ti, j 1 4Ti, j 0
• 所以，当i＝ N-1时，
aN 2
w 2 N 2N 1 x x N 2 N 1
2 cN 2 N x x e
bN 2 aN 2 cN 2 (S p ) N 1 x
d N 2 (Sc ) N 1 x cN 2 (TN 1)

(x)w (x)e (x)e
1
2
3
N-1
N
计算机实现：算例
• 所以，当i＝2时，
21 a1 x x w e 2 23 c1 x x 2 3
b1 a1 c1 (S p )2 x
d1 (Sc )2 x a1T0

求解上面的方程，即可得到（N－2）个未知数， 即，T2, T3, T4,……., TN-1。
计算机实现：算例
• 注意：上面的方程组是非线性的，必须用迭代 法求解 • 求解方法：
假定一个温度分布：Ti，i＝1，2，3，。。。，N 计算i ，i＝1，2，3，。。。，N 计算a, b, c, d 用TDMA法求解方程组，得到新的温度分布： Ti’ 计算：Max{abs(Ti -Ti’), i＝1,2,3,……,N} 判断： abs(Ti -Ti’)是否小于（精度要求） 如果不能满足精度要求，令Ti ＝Ti’，重复上面的计 算 满足精度要求：计算结束
以上方程适用于内部点: 对于节点或点
i3 j4
T4,4 T2,4 T3,5 T3,3 4T3,4 0
T3,4 1/ 4 (T4,4 T2,4 T3,5 T3,3 )
7
矩阵求解 I
对于线性代数：
A1AT A1C
T A C
[a] [aW ], [b] [aP ] [c] [aE ], [ D] [bP ]

由方程（1）系数阵[A]的特殊性，通常称之为三 对角方程（Tri-diagonal equation) 三对角方程可以采用非常高效的追赶法（TDMA 法）求解



基于矩阵分解 属于必须掌握的内容
Node 1 2 3 4 aE 1000 1000 1000 0 aW 0 1000 1000 1000 aP 3000 2000 2000 3000 b 2500 + 2000 TA 2500 2500 2500 + 2000 TB
代数方程的矩阵形式：
3000 -1000 0 0 -1000 2000 -1000 0 0 -1000 2000 -1000 0 0 -1000 3000 T1 T2 T3 T4 2000 TA + 2500 2500 2500 2000 TB + 2500