计算传热学第4讲扩散方程的数值解讲解

计算传热学 第4讲
扩散方程的数值解
Numerical Solution of Diffusion Equations
主要内容
一维稳态问题的数值解 一维非稳态问题 多维非稳态问题的离散化 差分方程的求解
主要目的
掌握用数值方法求解传热问题的整体步骤 数值方法的计算机实现 边界条件的处理
所有方法都适用
– 把跃变界面作为边界
4.1 一维稳态导热问题的数值解
把跃变界面作为边界
– 可以考虑接触热阻rc (W·m2)/K – 满足流的唯一性原则，
qi


dT dx
i


dT dx
i
Ti
Ti
Ti

T i

qirc
4.1 一维稳态导热问题的数值解
e
(x)1
2
O[(x)12 ]
(13)
边界条件的处理


dT dx
x0

qB
(8)
所以，
dT T2 T1 S (x)1
(14)
dx x0 (x)1 e 2
将之代入式（8）


dT dx
x0

1
dT dx
x0

1[
边界条件的处理


dT dx
x0

qB
(8)
代入式（16），整理后得到，
dT dx
x0

(x)2
1 (1

Lx
)
[T3
(Lx
1)2T2
(2
Lx )LxT1]
(17)
代入方程（8），
1 (x)2 (1
Lx
)
[T3

(Lx
1) 2 T2

(x)e
(x)-e
Ww P
(x)-w
x
eE (x)+e
图 1 一维问题空间区域的离散化
4.1 一维稳态导热问题的数值解
e
e
SA(x)dx w
w (Sc S pT ) A(x)dx
(Sc A)P x (S p A)P xTP
(6)
在上面的积分过程中，我们假定：
qB
1
e
法
(x)1
dT T2 T1
(9)
dx x0 (x)1



dT dx
x0

1
T2 T1
(x)1

qB
2
3
(x)2
边界条件的处理
整理后得到：
T1
T2

(x)1 1
qB
(10)
特点：
最简单的处理方法
只有一阶精度
与控制方程的精度不匹配
阅读与作业
阅读要求：陶文铨《数值传热学》第4章 作业：P124 题4－1；P125 题4－7 完成课外作业第一题和第二题
4.1 一维稳态导热问题的数值解
控制方程：
1 A( x)
d dx
A( x)
dT dx


S

0
(1)
其中，A(x)是面积函数。定义如下：
直角坐标系： A(x)＝1（无限大平板导热问题）
– 线性插值法（算术平均） – 调和平均法 – 待求变量插值 – Kirchhoff变换法
4.1 一维稳态导热问题的数值解
4.1.5 跃变界面的处理
– 把跃变界面作为控制面
调和平均法
– 把跃变界面作为节点
算术平均法 Kirchhoff变换法 待求变量插值法
– 把跃变界面放置在其它位置
– 待求变量T在控制容积P上为常数 – 整个控制容积的A(x)为常数，且等于P点的值。
4.1 一维稳态导热问题的数值解
将（5）和（6）代入方程（4），

A x
e
(TE

TP
)


A x

w
(TP

TW
)

(Sc A)P x (S p A)P xTP 0
dT dx
x0

dT dx
2

d 2T dx2
(x)1 O[(x)12 ]
2
(16)
按二阶精度的差商公式
dT dx
2

1
(x)2 (1 Lx ) [T3
L2xT1
(L2x
1)T2 ]
d 2T dx2
2

2 1 Lx
T3
(1 Lx )T2
(x)1(x)2
LxT1
边界条件的处理
– 元体能量平衡法：
在研究边界节点所代表的控 制容积（元体）的能量平衡
– 流入CV的能量＋内热源发 出的热量＝流出CV的能量
– 流入CV的能量
– 通过边界流入的热量qB – 通过控制面流入的热量q1 – 内热源发出的热量＝ ½(x)1S – 流出CV的能量＝0
x=0
½(x)1 e
边界条件的处理
x=0 (x)1
(x)2
qB 1
2
3
(x)2
( x)3
说明：
– 尽管它与一阶Taylor级数展开法的结果形式上相同， 但它却是二阶精度的！
– 请大家证明这一结论。
采用内节点法划分网格时，即使在均匀网络的前 提下，第1个近边界节点也不是等步长的。
从图中可以清楚地看出这一点 即使 (x)2= (x)3
T2 T1
(x)1

S
e
(x)1 ]
2

qB
边界条件的处理
整理后得到，
T1
T2

(x)1 e

e 1
qB

1 2

(x)1S
(15)
特点
二阶精度 不具有一般性 推导繁琐
边界条件的处理
x=0
e
qB 1
e
(x)1
2
3
(x)2
二阶精度的Taylor级数展开1
qB 1
2
3
e
(x)1
(x)2
q1

e
T2 T1
(x)1
边界条件的处理
代入能量守恒关系，
qB

e
T2 T1
(x)1

1 2
(x)1 S

0
整理后得到，
特点
T1
T2

(x)1 e
[qB


1 2
(x)1 S]
(11)
灵活，便于处理各种复杂的边界条件

e (x)
2

S p (x)2 T2


e (x)2
T3

Sc (x)2

qB
(24)
或者写成，
aPT2 aET3 Sc' (x)2
(25)
其中，
边界条件的处理
其中
aP

e (x)2

S p (x)2
Sc'

Sc

qB (x)2
Sc Sc,ad
( x)1也不等于 ( x)2 所以要对第一个内部节点给予特别注意。
x=0 (x)1
(x)2
qB 1
2
3
边界条件的处理
注意：
(x)2
( x)3
例如，对于直角坐标系，对C点于VW2节（的点节左点2控控，1制）制面重面w合e！与，节即
aPT2 aWT1 aET3 b2 与左边(界20重) 合！
边界条件的处理
代入与第1个近边界节点的差分方程（21），

1 (x)1

e (x)2

S p (x)2 T2


1 (x)1
T2

(x)1 1
qB



e (x)2
T3

Sc
(x)2
(23)
整理后得到，
边界条件的处理
整理后得到，
aE

e (x)2
(25)
Additional
source term!
Sc,ad

qB (x)2
(26)
边界条件的处理
对于第3类边界条件，也可以做类似的处理，但 是这时，
qB (Tf T1)
(27)
S Sc' Sp' T Sc' Sc Sc,ad
S
' p
Sp
A x e

(x)e (x)e 2(x)e
SP AP
(8a) (8b)
aP

aW
aE

(S p A)P 2
2x
(x)w
(x)
w
(x)w

(x)e (x)e (x)e

(8c)
4.1 一维稳态导热问题的数值解
4.1.4交界面参数的计算
柱坐标系： A(x)＝x（极坐标系中的一维问题，无限 长圆筒壁导热问题）
球坐标系：A(x)＝ x2 （通过球壁的导热）
变截面肋片： A(x)
4.1 一维稳态导热问题的数值解
4.1.1求解区域的离散化
(x)w (x)+w
(x)e (x)-e
W
w
P
(x)-w
x
图 1 一维问题空间区域的离散化
(12)
e
qB 1
2
3
e
(x)1
(x)2
边界条件的处理
x=0
e
qB 1
e
(x)1
2
3
(x)2
– 另一方面，参照附图，
dT dx
x0

dT dx
e

d 2T dx2
e
(x)1
2
O[(x)12 ]

dT dx
e

S
e
(x)1
2
O[(x)12 ]

T2 T1
(x)1

S
(7)
整理后得到，
aPTP aWTW aETE bP
(8)
其中，
4.1 一维稳态导热问题的数值解
aPTP aWTW aETE bP
(8)
其中，
aW
A x w
(8a)
aE


A

x e
(8b)
aP aW aE (S p A)P x
e
E
(x)+e
4.1 一维稳态导热问题的数值解
4.1.2 源项的线性化
– 在通常情况下，S＝S（T） – 线性化：
S＝Sc＋SpT – 其中，按负斜率源项原则，
Sp＝Sp(T*) 0
（2） （3）
4.1 一维稳态导热问题的数值解
4.1.3 控制方程的离散化
– 将方程（1）两边通乘A（x），并对x从w到e积分：
二阶精度，与内部节点精度等级匹配（请证明！！）
推导过程较繁
边界条件的处理
控制方程法
– 在直角坐标的条件下，方程（1）变为，
d dT S 0
dx dx
– 假定在节点1～2之间的导热系数为常数，且恒等于e，
则有，对于点e，
x=0
d 2T dx2
e
S
e
0
4.1.6 边界条件的处理
– 直角坐标
x=0
– 左边界，第二类边界条件 qB


dT dx
x0

qB
(8)
注意：qB的正方向与x轴的正方向一致！！
边界条件的处理
– 网格是用外节点法划分的
边界上出现半个控制容积
x=0
e
边界节点的差分方程可以用
下述方法推出： 一阶精度的Taylor级数展开
d dx
A( x)
dT dx


SA(
x)

0
e w
d dx
A(
x)
dT dx
dx

e
SA(x)dx 0
w
(4)
(x)w
(x)e
(x)+w
(x)-e
Ww P
(x)-w
x
eE (x)+e
图 1 一维问题空间区域的离散化
4.1 一维稳态导热问题的数值解
(8c)
bP (Sc A)P x
(8d)
下标：大写字母表示在节点处取值，小写字母表示在 相应的控制面处取值
4.1 一维稳态导热问题的数值解
可能的改进方案：对源项积分时采用线性分布
aW

A

x w

(x)w (x)w 2(x)w
SP AP
aE




1 (x)1
T1


e (x)2
T3

Sc
(x)2
(21)
边界条件的处理
附加源项法（Additional source term method)
– 以内节点法为例 – 由方程（19）解出边界节点上的待求变量T1，
T1

T2

(x)1 1
qB
(22)
代入与第1个近边界节点的差分方程（21），
(2

Lx
) LxT1 ]

qB
整理后得到，
(2

Lx
)LxT1

(Lx
1) 2 T2

T3

(x)2
(1
1
Lx
)
qB
(18)
边界条件的处理
二阶精度 具有一般性 增加计算工作量 一般很少采用
T1

T2

(x)1 e

e 1
qB

1 2

(x)1S
aW
x w

1 (x)1
aE

x
e

e (x)2
aP aW aE (S p A)2 (x)2
bP (Sc A)2 (x)2
边界条件的处理
或写成，

1 (x)1

e (x)2

S p (x)2 T2
e w
d dx
A(
x)
dT dx
dx

A(
x)
dT dx

e

A(
x)
dT dx

w

(A)e

dT dx
e

(A)w

dT dx
w

(A)e
TE TP
(x)e
(A)w
TP TW
(x)w
(5)
(x)w
(x)+w