东南大学传热学课件第四章导热问题数值解法2

n,i
O • 当 ∆τ 足够小时， (∆τ ) 可略而不计，此时非稳态项的差分格 式可表示为 ( ( ∂t t ni +1) − t ni ) = ∂τ n ,i ∆τ
非稳态项的三种差分格式
• 向前差分
( ( t ni +1) − t ni ) ∂t = ∂τ ∆τ
• 向后差分
∂t t n − t n = ∂τ ∆τ
2
α
=
0.4 × 0.01
2
α
二维非稳态导热内节点差分方程的建立
• 控制方程 • 差分方程
( ( t mi ,+n1) − t mi ,)n
∂ 2t ∂ 2 t ∂t = α 2 + 2 ∂x ∂τ ∂y
∆τ
( ( t mi ,+n1) − t mi ,)n
() () ( ( ( ( t mi +1,n + t mi −1,n − 2t mi ,)n t mi ,)n +1 + t mi ,)n +1 − 2t mi ,)n = α + 2 2 ∆x ∆y
( +1 ( +1 ( ( ) ( ) ( t mi +1,)n + t mi −1,)n − 2t mi ,+n1) t mi ,+n1+1 + t mi ,+n1+1 − 2t mi ,+n1) = α + 2 2 ∆x ∆y
∆τ
二维非稳态导热内节点的差分方程
Φ y+dy
(i +1)
(i +1)
wk.baidu.com
(i +1)
两种差分格式的比较
• 显式差分格式
( ( ( ( ( t ni +1) − t ni ) t ni+)1 + t ni−)1 − 2t ni ) =α ∆τ ∆x 2
• 隐式差分格式
tn
(i +1)
− tn t n +1 + t n −1 − 2t n =α ∆τ ∆x 2
第三节 非稳态导热问题的数值解
非稳态导热与稳态导热的主要区别在于控 制方程中多了一个非稳态项，而扩散项的 离散方法与稳态导热是一样的。因此，本 节将重点讨论非稳态项的离散方法以及扩 散项离散时所取时间层的不同对计算带来 的影响。
一维非稳态导热问题的离散
τ
n,i+1 n-1,i n,i n,i-1 n+1,i
[
]
n = 1,2,3 ⋅ ⋅⋅, N − 1
( ( ( t Ni +1) = t Ni ) (1 − 2 Fo∆ Bi∆ − 2 Fo∆ ) + 2 Fo∆ t Ni )−1 + 2 Fo∆ Bi∆ t f
( ( t2i ) = t−i1)
方程组的求解
• 利用上述方程组，从初始温度 t0 出发，即可依次 求得第二时间层、第三时间层直到 I 时间层上的 温度分布。至于空间步长 ∆τ及时间步长 ∆x 的选取， 原则上步长越小，计算结果越接近于精确解，但 是需要的计算机内存及计算时间则大大增加。此 ∆ 外， τ 与 ∆x 的关系还受到显式差分格式稳定性的 影响。 • 下面，我们从离散方程的结构来分析，说明稳定 性限制的物理意义，再通过数值计算实例予以说 明。
N-2
N-1
N
x
无限大平板换热边界上节点方程的推导
• 从左侧面进入元体的热量 • 从右侧面进入元体的热量
( ( t Ni )−1 − t Ni ) Φx = λ ∆x
( Φ x + dx = h t f − t Ni )
(
)
( ( ∆x t Ni +1) − t Ni ) • 元体自身热力学能的改变量 ∆E = ρc 2 ∆τ
显式差分格式稳定性条件
• 内节点差分方程稳定性条件
1 1 − 2 Fo∆ ≥ 0 Fo∆ = 2 ≤ ∆x 2
α∆τ
• 一维非稳态导热，换热边界上节点差分方 程稳定性条件
1 − 2 Fo∆ Bi∆ − 2 Fo∆ ≥ 0
1 Fo∆ ≤ 2(Bi∆ + 1)
具体计算实例
• 题目：厚 2δ
= 0.06m的无限大平板受对称的冷却， 初始温度 t0 = 100 ℃。在初始瞬间，平板突然被置 于温度 t∞ = 0℃的流体中。已知平板的导热系 数 λ = 40W/ (m ⋅ K ), h = 1000W/ (m 2 K )。试用数值法求解
(i )
(i +1)
(i +1)
(i +1)
两种差分格式的区别
• 格式的形式不同 • 计算工作量不同 显式格式计算工作量小，隐式格式计算工作量大 • 限制条件不同 显式格式对时间步长和空间步长有相互制约的要 求，但隐式格式对时间步长及空间步长之间的关 系没有任何要求。
建立非稳态导热问题 节点方程的热平衡法
2δ
x
控制方程
解：由于问题的对称性，只要研究一半即可，此时，该问 题的控制方程为
∂t ∂ 2t = a 2 ( 0 < x < δ ,τ > 0 ) ∂x ∂τ t ( x,0 ) = t0 ( 0 ≤ x ≤ δ ) ∂t ( x,τ ) =0 ∂x x =0 ∂t ( x,τ ) h t ( x,τ ) − t∞ = −λ ∂x x =δ
m,n+1
λ
m+1,n
∆y + λ ∆y ∆ ∆ 144 x 44 2 3 144 x 44 2 3
左侧面导入元体的热量 右侧面导入元体的热量
() ( t mi −1,n − t mi ,)n
• 化简结果 t n
(i +1)
=
α ∆τ
∆x
2
[t
n +1
(i )
+ t n −1
(i )
]
2α∆τ + 1 − ∆x 2
(i ) t n
一维非稳态导热内节点方程的建立
• 控制方程
∂t ∂τ
(i +1)
n ,i
∂ 2t =α 2 ∂x
(i )
n ,i
• 差分方程
tn
− tn t n +1 + t n −1 − 2t n 2t =α 2 ∆τ ∆x
其温度分布。取 Fo∆ = 1。 • 解：区域离散化，取 ∆x = 0.01m • 则 Bi = h∆x = 1000 × 0.01 = 0.25
∆
λ
40
• 采用如图所示的离散方法，计算结果列于下表
计算区域离散图
n=1 2 3 4
0.01
0.01
0.01
x
差分方程
上述问题的差分方程为
( ( ( ( t ni +1) = Fo∆ t ni+)1 + t ni−)1 + (1 − 2 Fo∆ )t ni )
一维非问题导热问题的差分方程
• 内节点差分方程
( ( ( ( t ni +1) = Fo∆ t ni+)1 + t ni−)1 + (1 − 2 Fo∆ )t ni )
[
]
• 换热边界上的差分方程
( ( ( t Ni +1) = t Ni ) (1 − 2 Fo∆ Bi∆ − 2 Fo∆ ) + 2 Fo∆t Ni )−1 + 2 Fo∆ Bi∆t f
• 将研究区域离散化 • 对各节点所代表的元体建立能量平衡关系式 • 对非稳态导热问题该能量平衡关系式为 从各个方向进入元体的热量之和等于该元体热力 学能的变化量 • 整理化简，得到各节点的差分方程
无限大平板换热边界上节点方程的建立
• 左图示出了一无限大 平板的右侧面的一部 分，其右侧面受到周 围流体的冷却，表面 传热系数为h，流体温 度为 t f • 边界节点为N • 节点N 代表宽度为∆x / 2 的元体
[
]
n = 1,2,3
( ( ( t Ni +1) = t Ni ) (1 − 2 Fo∆ Bi∆ − 2 Fo∆ ) + 2 Fo∆ t Ni )−1 + 2 Fo∆ Bi∆ t f
N = 4, N − 1 = 3
( ( t2i ) = t−i1)
计算结果
t n
0
i
0
1
2
3
4
5
6
7
100
100
100
n
非稳态项的离散
• 如果将函数t 在节点（n,i+1）对点（n,i）作泰勒展开，可 有 ∂t ∆τ 2 ∂ 2t ( ( t ni +1) = t ni ) + ∆τ + + ⋅⋅⋅ 2 ∂τ n ,i 2 ∂τ n ,i • 于是有 ∂t ∂τ
( ( t ni +1) − t ni ) = + O (∆τ ) ∆τ
显式差分格式稳定性分析
由内部节点差分方程可见，在节点n上，i+1时刻 的温度是在该点i 时刻温度的基础上考虑了左右相 邻两点温度的影响后得出的。现在，假设相邻两 点的温度不变，那么合理的情况是:i时刻节点n的 温度越高，则其相继时刻（i+1时刻）的温度也越 高；反之，i时刻节点n的温度越低，则其相继时 刻的温度也越低。所以，在差分方程中要满足这 ( ( tni +1)与 tni )前面的系 种合理性的条件，则差分方程中 ( i +1 ) 数必须保持同方向变化。由于 tn 的系数大于零， 因此 tn(i ) 前面的系数也必须大于零 。
出现解的不稳定性的原因
• 不满足解的收敛条件 • 根据内节点的收敛条件 Fo∆ ≤ 1 / 2 = 0.5 • 根据换热边界节点的收敛条件
1 1 = = 1 / 2.5 = 0.4 Fo∆ ≤ 2(Bi∆ + 1) 2(0.25 + 1)
• 所以此时最达的时间间隔应为
∆τ = Fo∆ min ∆x
∂t ∂τ =
n ,i
(i )
(i −1)
• 向中心差分
tn
(i +1)
− tn 2∆τ
(i −1)
一维非稳态导热内节点方程的建立
• 控制方程 ∂t ∂τ
n ,i
∂ 2t =α 2 ∂x
n ,i
( ( ( ( t ni +1) − t ni ) t ni+)1 + t ni−)1 − 2t ni ) 2t ( • 差分方程 =α 2 ∆τ ∆x
• 根据能量守恒定律 Φ x + Φ x + dx = ∆E
( ( • 整理后得到 t Ni +1) = t Ni ) 1 − 2h∆τ − 2α∆2τ ρc∆x ∆x
2α∆τ (i ) 2h∆τ + ∆x 2 t N −1 + ρc∆x t f
差分方程的进一步演化
• 考察方程中的
区域离散化
• 将所研究平板的一半N等分， 共有N+1个节点，其中节 点1在平板中心截面上，节 点N在平板右侧面上，如图 所示 • 两个节点之间的距离为 ∆x • 节点-1与节点2换热情况对 称，固有相同的温度 • 时间步长取 ∆τ
-1
n=1 2
3
N-1
N
x
差分方程
上述问题的差分方程为
( ( ( ( t ni +1) = Fo∆ t ni+)1 + t ni−)1 + (1 − 2 Fo∆ )t ni )
100
60
148
-109.6
550
1
100
100
100
80
104
19.2
220.2
-328.9
2
100
100
80
84
63.2
91.4
0.9
220
3
100
80
64
67.2
50.6
73.1
0.72
176
对计算结果的说明
从上表可以看出，从 i = 3 这一时刻起出现了这样 的情况：各点温度随时间作忽高忽低的波动，并 且波动幅度越来越大；某点温度越高反而使相继 时刻的温度越低。这种现象是违背热力学第二定 律的。因为这意味着，在该时间间隔中，从某一 时刻起热量将自动地由低温向高温传递。数值计 算中出现的这种计算结果忽高忽低的波动现象， 数学上称为不稳定性。这个例题表明，在数值计 算中避免出现不稳定性是十分重要的。
x
• x为空间坐标，将计算区域划分 为（N-1）等份，得到N个空间 节点；两节点之间的距离为∆x 称为空间步长； • τ为时间坐标，将时间坐标上的 计算区域划分为（i-1）等份， 得到 i 个时间节点。 • 从一个时间层到下一个时间层 的间隔为∆τ，称为时间步长。 • 空间网格与时间网格的交点， 如（n,i），代表了时间—空间 区域中一个节点的位置，相应 的温度记为t (i ) 。
2 h∆ τ ρc∆x
h∆τ λ ∆τ h∆x α∆τ h∆x = = = Fo∆ ⋅ Bi∆ 2 2 ρc∆x ρc ∆x λ ∆x λ
• 定义 Fo∆ = • 定义 Bi∆ =
α ∆τ
∆x 2
称为网格傅立叶数
h∆x
λ
称为网格毕渥数
• 差分方程可演化为
( ( ( t Ni +1) = t Ni ) (1 − 2 Fo∆ Bi∆ − 2 Fo∆ ) + 2 Fo∆t Ni )−1 + 2 Fo∆ Bi∆t f
• 上述方程都是用显式差分格式表示的
数值解求解一维非稳态导热的实例
物理模型：设有一块 厚度为2 δ 的无限大平 壁，初始温度为 t0 。 在初始瞬间将它放置 于温度为 t f 的流体中， 流体与板面间的表面 传热系数h为常数。试 用数值解法确定在非 稳态导热过程中板内 的温度分布。
h, t f
h, t f