热传导与有限差分

13
研究生课程《高等工程热力学与传热学》电子教案
常见差分方程求解
1) 显式差分方程： 直接求解。 2) 隐式差分方程或者 Crank-Nicholson 差分方程： 一维问题：追赶法； 二维问题：高斯消元法或者迭代法（如高斯－赛德尔法） 。 快速算法－交替方向法（俞昌铭，1981）
14
研究生课程《高等工程热力学与传热学》电子教案
) (
)]
(22)
该差分格式称之为 Crank-Nicholson 格式，它的精度要高于显 式格式和隐式格式，其截断误差为 O[( Δt ) 2 + ( Δx ) 2 ] 。
9
研究生课程《高等工程热力学与传热学》电子教案
Baidu Nhomakorabea常见差分格式的稳定性分析
所谓稳定性问题，是指：若定义第 n 步第 j 个节点的数值求解 误差为 ε n j ，如果定义误差放大倍数 ω 满足
2
(6)
研究生课程《高等工程热力学与传热学》电子教案
三类边界条件
第 1 类边界条件： 给定边界上的温度。 第 2 类边界条件： 给定边界上的法向热流密度。 第 3 类边界条件： 给定外部介质的温度和给定边界上的对流换热面的对流 换热系数。
3
研究生课程《高等工程热力学与传热学》电子教案
偏微分方程的求解
Δt
2
15
T jn +1 − T jn
=a
+1 n +1 +1 T jn + T jn + 1 − 2T j −1
, n + 1 + j = 奇数
(33)
研究生课程《高等工程热力学与传热学》电子教案
跳点差分格式的应用
D2 j=1 Δ x1 Δx Δ x2 j=2 i-1, j i, j i, j+1 i+1, j
⎛ ΔR ⎞ ΔX 1 = ΔX ⎜ + ⎟ 1 ⎜ D1 + D2 ⎟ ⎝ ⎠
(34)
内表面对应单元的弧长为
⎛ ΔR ⎞ ΔX 2 = ΔX ⎜ − ⎟ 1 ⎜ D1 + D2 ⎟ ⎠ ⎝ 2 ， 是 壁 厚 ， D1
T jn 直接求得，故这种差分格式称之为显式差分格式。现在已经知道
其截断误差为 O[Δt + ( Δx ) 2 ]。 如果采用朝后时间差分，则 n+1 时刻 T 对时间 t 的偏导为：
n +1 n ∂T T j − T j = ∂t Δt
(19) (20) (21)
对坐标 x 的二阶偏导用 n+1 时刻的中心差分代替
m =1
10
研究生课程《高等工程热力学与传热学》电子教案
由于微分方程是线性的，级数中每 1 项对其后各时刻的影响可 以叠加，因此可以去掉符号 ∑ 和下标 m，常数 A 也可以不考虑。这 样，经过离散化处理，则 n 时刻节点 j、j+1 和 j-1 处的误差分别为：
⎧ε n j = cos( β iΔx ) ⎪ ⎪ n ⎨ε j +1 = cos[ β ( i + 1) Δx ] = cos( βiΔx ) cos( βΔx ) − sin( βiΔx ) sin( βΔx ) ⎪ n ⎪ ⎩ε j −1 = cos[ β ( i − 1) Δx ] = cos( βiΔx ) cos( βΔx ) + sin( βiΔx ) sin( βΔx )
一、解析解：
稳态问题 分离变量法 虚拟热源法 瞬态问题 近似积分法 格林函数法 拉普拉斯变换法
二、数值解：
有限差分法； 有限元法； 边界元法； 无网格法。
4
研究生课程《高等工程热力学与传热学》电子教案
有限差分法的基本概念
设函数 f(x)可以展开为泰勒展开式：
df ( Δx ) 2 d 2 f ( Δx ) 3 f ( x + Δx ) = f ( x ) + Δx + + dx 2! dx 2 3! df ( Δx ) 2 d 2 f ( Δx ) 3 − f ( x − Δx ) = f ( x ) − Δx + dx 2! dx 2 3! d3 f + ... 3 dx d3 f + ... 3 dx
(30)
由此可见，隐式差分格式无条件稳定。 3) Crank-Nicholson 格式 参考前面的方法，可得隐式差分格式的误差放大倍数为
2 aΔ t 2 β Δx sin ( ) 2 2 ( Δx ) ω= 2 aΔ t 2 β Δx 1+ sin ( ) 2 ( Δx ) 2 1−
(31)
Crank-Nicholson 格式也是无条件稳定的。
4 aΔ t 2 β Δx sin ( ) ≤1 2 2 ( Δx ) aΔ t 1 ≤ ( Δx ) 2 2
(28)
=1−
要让式(28)稳定，那么必须满足
1− ⇒
(29)
12
研究生课程《高等工程热力学与传热学》电子教案
2) 隐式格式 参考前面的方法，可得隐式差分格式的误差放大倍数为
ω=
4aΔt 2 β Δx 1+ sin ( ) 2 2 ( Δx ) 1
将式(25)代入式(26)，可得
11
研究生课程《高等工程热力学与传热学》电子教案
+1 εn = j
⎡ 2 aΔ t ⎤ 2 aΔ t β Δ β Δ + − cos( βiΔx ) cos( i x ) cos( x ) 1 ⎢ 2 2⎥ ( Δx ) ⎣ ( Δx ) ⎦
(27)
将式(27)、(25)代入误差放大倍数公式（式(23)） ，则有
df f ( x + Δx ) − f ( x − Δ x ) 2 = + O (Δx ) dx 2 Δx
[
]
[
]
(12) (13) (14)
此外，还可以有朝前差分公式和朝后差分公式：
有限差分的核心思想就是用各阶差分代替各阶微分。
6
研究生课程《高等工程热力学与传热学》电子教案
瞬态问题的几种常见差分格式
(2) (3)
如果无内热源，则式(2)可以进一步简化为：
∂T k 2 = ∇T ∂t ρc
1
研究生课程《高等工程热力学与传热学》电子教案
正交坐标系中的导热微分方程
直角坐标系下的导热微分方程：
ρc ∂T
∂ 2T ∂ 2T ∂ 2T ⋅ = + + k ∂t ∂ x 2 ∂ y 2 ∂z 2
(4)
柱坐标系下的导热微分方程：
(25)
由于此处讨论的是线性齐次的差分方程，那么对于温度 T 和误 差 ε ，应该各自满足差分方程。于是考虑显式差分格式（即式(18)） ， 则有
+1 εn = j
(Δ x )
aΔt
2
(ε
n j −1
⎡ aΔt ⎤ n 1 2 ε +εn + − j +1 ⎢ 2⎥ j (Δ x ) ⎦ ⎣
)
(26)
+1 εn j ω = n ≤1 εj
(23)
则称该数值解法是稳定的，反之则为不稳定。常见的稳定性分 析方法是傅里叶级数法。一般来说，稳态问题的差分方法不存在稳 定性问题。 1) 显式格式 假设误差函数 ε ( Δx , Δt ) 可以展开成傅里叶级数的形式： ∞ ε ( Δx , Δt ) = ∑ Am cos β m Δx (24)
跳点差分格式
把各节点按 n+j 的奇偶不同分为两组，当时间从 tn 推进到 tn+1 时，先在偶数节点上用显式差分格式：
(Δx ) 再在奇数节点上用隐式差分格式：
Δt
2
T jn +1 − T jn
=a
n n T jn + 1 − 2T j + T j − 1
, n + 1 + j = 偶数
(32)
(Δx ) 注意，由于式(33)中的 T jn++11 、 T jn−+11 位于偶数节点上，实际上已经 由式(32)算出。这样，式(33)形式上是隐式的，算法设计上仍然是显 式的。已经证明，跳点差分格式的截断误差为 O[( Δt ) 2 + ( Δx ) 2 ] ，并且无 条件稳定。
(7) (8) (9)
对以上两式进行改写，可得：
df f ( x + Δx ) − f ( x ) ( Δx ) d 2 f ( Δx ) 2 = − − dx Δx 2 dx 2 6 df f ( x ) − f ( x − Δ x ) ( Δ x ) d 2 f ( Δx ) 2 = + − 2 dx Δx 2 dx 6 d3 f − ... dx 3 d3 f + ... 3 dx
将式(16)、(17)代入式(15)，经过整理，可得：
T jn +1 =
(Δx )
aΔt
2
(T
n j −1
⎡ aΔt ⎤ n 1 2 T + T jn + − +1 ⎢ 2⎥ j (Δx ) ⎦ ⎣
)
7
研究生课程《高等工程热力学与传热学》电子教案
由于 n+1 时刻某节点T jn+1 的温度可以由 n 时刻的温度T jn−1 、 T jn +1 和
(10) 如果将式(9)和式(10)相加，即可得一阶导数的中心差分公 式：
5
研究生课程《高等工程热力学与传热学》电子教案
(11) 如果将式(9)和式(10)相减，即可得二阶导数的中心差分公 式：
d2 f f ( x + Δx ) − 2 f ( x ) + f ( x − Δx ) 2 ( ) O x = + Δ dx 2 (Δx )2 df f ( x + Δx ) − f ( x ) + O (Δ x ) = dx Δx df f ( x ) − f ( x − Δx ) + O (Δ x ) = dx Δx
热传导与有限差分
讲授人：郭蕴华 李格升 武汉理工大学能源与动力工程学院 热能工程系
研究生课程《高等工程热力学与传热学》电子教案
热传导的数学描述
有内热源、各项同性的物体的导热微分方程：
ρc
∂T = ∇ (k∇T ) + q v ∂t
(1)
如果导热系数不随空间位置和温度的变化，则式(1)简化 为：
q k 2 ∂T = ∇ T+ v ρc ρc ∂t
考虑常系数一维抛物型方程：
∂T ∂ 2T =a 2 ∂t ∂x
(15) (16) (17) (18)
如果采用朝前时间差分，则 n 时刻 T 对时间 t 的偏导为：
n +1 n ∂T T j − T j = ∂t Δt
对坐标 x 的二阶偏导用 n 时刻的中心差分代替
n n n ∂ 2T T j +1 − 2T j + T j −1 = ∂x 2 (Δx )2
+1 εn j ω= n εj
⎡ 2aΔt 2aΔt ⎤ Δ Δ + − β i x β x cos( ) cos( ) 1 cos( βiΔx ) ⎢ 2 2⎥ ( Δx ) ⎣ ( Δx ) ⎦ = cos( βiΔx )
=
⎡ 2aΔt 2aΔt ⎤ 2aΔt β Δ x + − = − cos( ) 1 1 [1 − cos( βΔx )] ⎢ 2 2⎥ 2 ( Δx ) ( Δx ) ⎣ ( Δx ) ⎦ 4aΔt 2 β Δx sin ( ) 2 2 ( Δx )
由式(21)被称之为隐式差分格式，其截断误差也是 O[Δt + ( Δx ) 2 ]。
8
研究生课程《高等工程热力学与传热学》电子教案
如果把 T 对 x 的二阶偏导数用式(17)和式(20)的算术平均表示， 则有
T jn +1 − T jn Δt = 2(Δx ) a
2
[(T
n j +1
n +1 n +1 n+1 − 2T jn + T jn + T − 2 T + T −1 j +1 j j −1
Δα
ΔR D1 j=NN
图 1 容器壁的离散化 对水平圆柱体容器的圆筒壁按圆周方向进行离散， 按左右两个半圆各 分成 NN = 18 个等份，每等份对应相同的单元角 Δα ，在壁厚度方向（半径
16
研究生课程《高等工程热力学与传热学》电子教案
方向）划分两等份。整个容器圆筒壁划分成 36 个单元。如图 1 所示（示 意图中仅分为 16 个单元） 。 外表面对应单元的弧长为
ρc ∂T 1 ∂ ⎛ ∂T ⎞ 1 ∂ 2T ∂ 2T = + 2 ⎜r ⎟+ 2 2 k ∂t r ∂r ⎝ ∂r ⎠ r ∂ϕ ∂z
(5)
球坐标系下的导热微分方程：
ρc ∂T
∂ ⎛ ∂T ⎞ 1 ∂ ⎛ 2 ∂T ⎞ 1 1 ∂ 2T ⋅ = 2 ⎜r ⎟+ 2 ⎜ sin θ ⎟+ 2 2 k ∂t r ∂r ⎝ ∂r ⎠ r sin θ ∂θ ⎝ ∂θ ⎠ r sin θ ∂ϕ 2
n +1 n+1 n +1 ∂ 2T T j +1 − 2T j + T j −1 = ∂x 2 (Δx )2
将式(19)、(20)代入式(15)，经过整理，可得：
− ⎛ aΔt n +1 aΔt ⎞ n +1 n +1 n ⎜ ⎟ + + 1 2 T T T T = − j − 1 j j + 1 j ⎜ (Δx )2 (Δx )2 ⎟ (Δx )2 ⎝ ⎠ aΔt