第四章扩散方程的数值求解

中国科学院研究生院2010年春季
5
按这种方式编制程序时，只要设置一个变量MODE， 它按以下方式取值，程序即可自动处理三种坐标系。
MODE R sx
1(x-y) 1 1
2(x-r) r 1
3(theta-r) r r
各种商业软件大致按照类似方式处理不同坐标系中的离散方 程系数的计算问题。
中国科学院研究生院2010年春季
迭代法
块迭代法
迭
迭交
代
代替
逐 线 迭 代 法
法
法方
向
逐
线
逐
线
中国科学院研究生院2010年春季
多维导热差分方程求解 差分方程
非线性差分方程
假定代求量分布，并计算非 线性方程的系数，使其线化
公式
TP − TP0 ∆t
=
a (TE 2
− 2TP + TW ∆x2
+ TE0
− 2TP0 ∆x2
+ TW0 ) 中国科学院研究生院2010年春季
T
T
T
T1
T1
T1
相邻时层 T0
T0
T0
温度分布
t0
t1
显式格式
t
t
t0 t1 C-N格式
t t0 t1 全隐格式
正系数准则 有条件满足 有条件满足 无条件满足
中国科学院研究生院2010年春季
2. 圆柱轴对称
3. 极坐标
aPTP = aETE + aWTW + aNTN + aSTS + b
aE
=
rP∆r (δ x)e
λe
aE
=
∆r rP (δθ )e
λe
中国科学院研究生院2010年春季
三种二维正交坐标系中离散方程的统一表达式
中国科学院研究生院2010年春季
求解线性代数方程的迭代法，不以有限步运算获得线性代数 方程组的真解为目的，而是以经过有限步运算获得线性代数方 程一定程度上的逼近解为目的
中国科学院研究生院2010年春季
6
Gauss-Seidel Jacobi
Gauss-Seidel Jacobi
2. 迭代法 迭代法类型
点迭代法
迭
代
法
逐
点
迭
代
逐
法
点
中国科学院研究生院2010年春季
非线性代数方程组的求解流程
给定节点上的温度值T* aPTP = aETE + aW TW + b
计算差分方程的系数和源项 a*PTP = a*ETE + a*WTW + b*
求解线性代数方程组，得到新的温度分布T
T → T*
Max T − T * < ε? Y 结束
N
中国科学院研究生院2010年春季
Pi
=
ai
bi − ciPi−1
Qi
=
di ai
+ ciQi−1 − ciPi−1
**
显然，回代或代入方程系数L公式是一个递推关系式， 为此，还需给出Pi和Qi的起始值P1和Q1。
由节点1的差分方程可得到递推公式的起始值如下：
P1
=
b1 a1
Q1
=
d1 a1
*
中国科学院研究生院2010年春季
当i=N时，利用i=N-1得到的TN-1和TN的关系代入i=N的差分方 程，则可把TN表示成TN+1的已知函数
中国科学院研究生院2010年春季
由于TN+1并不存在，这就意味着到i=N时，已经得到待求 TN的值了。
将解出的TN代入i=N-1时得到的TN-1和TN的关系中，即可 得到TN-1的值了。
将解出的TN-1代入i=N-2时得到的L TN-2和TN-1的关系中，即 可得到TN-2的值了。
例子
d2T d2x
+
f
(
x)
=
0
x
中国科学院研究生院2010年春季
3
T(K) T(K)
T(K) T(K)
x= 0 ∆x
1
i= 1 2 … i-1 i i+1 … ∆x=1/(n-1)
…
n
中国科学院研究生院2010年春季
500
prediction accuracy solution D-D,T(1)=300,T(100)=500,f(x)=600x
*
B B 中国科学院研究生院2010年春季
线对矩一 性角阵维 方特具嗥 程征有系 组的三数 TDMA
线性方程组 划组特具二 分嗣征有维 有带的带嗥 关宽线状系 与性稀数 网方疏矩 格程阵阵
划组特具三 分嗣征有维 有带的带嗥 关宽线状系
与性稀数 网方疏矩 格程阵阵
直
迭
直
迭
接
代
接
代
法
法
法
法
中国科学院研究生院2010年春季
aPTP = aETE + aW TW + aN TN + aSTS + aTTT + aBTB + b
1. 直接法
aPTP = aETE + a W TW + aN TN + aSTS + b
AT = b
A是一个主对角线及 其邻近若干层元素不 为零，而其他元素都 为零的带状稀疏矩阵
Cramer法则
外节点法
非紧邻边界节 点或紧邻边界 节点的差分方 程
边界节点差
分方程引入 边界条件
aB
TB
=
a E TE
+b
(a B TB
=
a W TW
+ b)
aPTP = aETE + b (aPTP = a W TW + b)
内节点法
紧邻边界节点差分方程，引入边界条件 中国科学院研究生院2010年春季
所有系数 与待求温 度无关
( Γϕ
∂ϕ ∂x j
)
+ Sϕ
不考虑对流
∂ρϕ ∂t
=
∂ ∂x j
(Γϕ
∂ϕ ∂x j
) + Sϕ
典型导热问题 的控制方程
基于导热问题的数值计算，来认识 如何处理边界问题及实现源项负线 化和差分方程求解
中国科学院研究生院2010年春季
源项负线化 源项负线化的目的是尽可能的把源项的影响通过 中心节点温度的系数反映出来
三种格式的比较
加权 系数f
相关 节点
显式格式 0
仅与0时层相应 邻点相关
C-N格式 0.5
既与0时层相应 邻点相关，又 与当前时层邻
点相关
全隐格式 1
只与当前时层邻点 相关
TP − TP0
计算 ∆t
= a TE0
− 2TP0 ∆x2
+ TW0
TP − TP0 ∆t
= a TE
− 2TP + TW ∆x2
aPTP = aETE + a W TW + b
系数是待 求温度的 函数
线性代数方程组
线性代数 方程组的 求解方法
非线性代数方程组 线化
温度场
假定温度场
中国科学院研究生院2010年春季
非线性代数方程组的求解步骤
1、在所有各个网格节点上，猜测或估计或假定一个T值 2、用这些估计的T值去计算差分方程中的所有系数，从而差 分方程中的所有系数变成了已知量，而使差分方程变成了线性 方程。 3、求解上边的线性方程组，得到各离散点新T值。 4、用新得到的T值去计算差分方程中的所有系数，并返回第 3步求解系数发生了变化的线性差分方程。重复3、4这个过程， 直到重复计算不在引起T值任何有意义的变化为止。
求解方法
直接求解
TDMA
TDMA
中国科学院研究生院2010年春季
4.4 多维导热基本方程与差分方程
1. 二维直角坐标系
ρc ∂T = ∂ (λ ∂T ) + ∂ (λ ∂T ) + S
∂t
∂x
∂xL
∂y
(δx)_e
∂y
(δx)+e
NnBiblioteka (δy)+s (δy)-s
∆y
Ww P e E
s
S ∆x (δx)e
将解出的T3代入i=2时得到的T2和T3的关系中，即可得到 T2值了。
将解出的T2代入i=1时得到的T1和T2的关系中，即可得到
T1的值了。
中国科学院研究生院2010年春季
TDMA的原理
所有待求变量T1,T2,…,TN就全部求出了。从而完成这个三 对角阵的线性代数方程组的求解
这种求解三对角矩阵的方法称之为三对角矩阵算法，常常 又称之为TDMA或追赶法
0.2
0.4
0.6
0.8
1.0
中国科学院研究生院2010年春季
x
4.3 一维非稳态导热基本方程与差 分方程
ρc∂T = 1 d [λA(x)dT]+S
∂t A(x) dx
dx
L (δx)w
w W 一维非稳态导热方程离散
(δx)e
(δx)e− (δx)e+
e
P ∆x
E
中国科学院研究生院2010年春季
4
计算传热学 第四章 扩散方程的数值求解
中科院研究生院2010年春季
4.1 导热问题
建立描述导热过程的基本方程 计算区域和控制方程的离散化
得到离散点满足的方程，即差分方程
存在一些问 题有待解决
目标：利用计算的方法获得导热过程的数值解
中国科学院研究生院2010年春季
尚待解决的问题 边界问题 源项负线化问题 差分方程如何求解问题
600 580 560 540 520 500
0.0
accuracy solution prediction N-D, T(100)=500, f(x)=600x
0.2
0.4
0.6
0.8
1.0
x
中国科学院研究生院2010年春季
500
450
400
350
300 0.0
accuracy solution prediction R-D,T(100)=500,alf=10.0,tf=350,f(x)=600x
多维导热差分方程的求解
aPTP = aETE + aW TW + aN TN + aSTS + aTTT + aBTB + b
多维导热差分方程 中的系数和源项本 身是待求变量函数
非线性方程组
线化、迭代
线性方程组
a*PTP
=
a*ETE
+
a*WTW
+
a*NTN
+
a*STS
+
a*TTT
+ a T + b *
中国科学院研究生院2010年春季
2
当i=1时，T1是T2的已知函数
当i=2时，利用i=1得到的T1和T2的关系代入i=2的差分方程，则 可把T2表示成T3的已知函数
当i=3时，利用i=2得到的T2和T3的关系代入i=3的差分方程，则 可把T3表示成T4的已知函数
L
当i=N-1时，利用i=N-2得到的TN-2和TN-1的关系代入i=N-1的差 分方程，则可把TN-1表示成TN的已知函数
Gauss消元法 计算量
计算次数正比于(N+1)!
所用乘法的次数是N**3次
当代数方程组的未知量数目很大时直接法 的计算量极大，甚至是不可接受的
中国科学院研究生院2010年春季
对有26个未知数的代数方程组，在每秒109次的亿次计算机 上用Cramer法则进行求解，大约需要1013年才能求完
当方程个数很大时，存储系数矩阵所需要的内存也是十分可 观的
450
400
350
300
0.0
0.2
0.4
0.6
0.8
1.0
x
中国科学院研究生院2010年春季
500
accurary solution prediction 450 D-D,T(1)=300,T(100)=500,f(x)=-600x
400
350
300
0.0
0.2
0.4
0.6
x
0.8 1.0 中国科学院研究生院2010年春季
正系数问题
不论什么过程的数值 计算中，均存在这些 问题
有对流过程的数值计算中存在 正系数问题，无对流过程则没 有正系数问题
鉴于一些问题与对流有关，因此本着循序渐进的原 则，先不考虑对流，在简单过程中先来认识前三类问 题如何解决
中国科学院研究生院2010年春季
∂ρϕ ∂t
+
∂ρu jϕ ∂x j
=
∂ ∂x j
由上边的分析可以看出，TDMA由代入和回代两个过程组 成，其核心是代入或回代方程，L 即反映Ti-1和Ti的函数关系式
TDMA的原理代入或回代方程 TDMA方法程序的实现
中国科学院研究生院2010年春季
假设Ti-1和Ti的函数关系如下所示： Ti−1 = Pi−1Ti + Qi−1
则将其代入i节点差分方程中 ai Ti = biTi+1 + ci (PL i−1Ti + Qi−1 ) + di
源项线化
S = Sc + S p TP
aPTP = aETE + a W TW + b
aP = aE + aW − Sp∆x
b = Sc∆x
线化系数为负
Sp < 0
中国科学院研究生院2010年春季
Sc和Sp的选取原则： Sp ≤ 0 当φ*＝φ 时，线性化方案（59）与源项表达 式取得完全一致
使Sp尽可能接近（d S/dφ ）（最优方案）
如果
dS dφ
≥
0,
取 Sp = 0
如果
dS dφ
<
0,
取
Sp
=
dS dφ
φ =φ *
中国科学院研究生院2010年春季
1
4.2 一维稳态导热
d (k dT ) + S = 0 dx dx
(δx)w
( δx ) e
w W
(δx)e− (δx)e+
e
P ∆x
这些都限制了直接法在求解线性方程组中的应用
中国科学院研究生院2010年春季
求解这组线化的方程组用直接法即使获得了它的真解，那也 仅仅只是我们获得非线性差分方程的近似解的过程中的一个临 时解－－一个有待通过迭代继续进行改进的解，所以我们花费 如此大的气力去取得线化代数方程组的真解是不值得的
求解过程实际上只要给出其一定程度的近似解即可，通过不 断的迭代，实现向非线性方程及其真解的逼近
E
aPTP = aETE + a W TW + b
中国科学院研究生院2010年春季
其中，
aE
=
ke ( δx ) e
aW
=
kw (δx) w
aP = aE + aW − Sp∆x
b = Sc∆x
S = Sc + Sp TP
中国科学院研究生院2010年春季
差分方程的求解
1 非线性差分方程的线性化
aPTP = aETE + a W TW + b