50 × 0.02 Bi1 = = = 0.01 λ 100
hδ
400 × 0.02 Bi 2 = = =1 λ 8
hδ
第四章 热传导问题数值解法
(i ) N
式中 Fo∆ =
a∆τ 网格傅里叶数 ∆x 2
h∆τ λ ∆τ h∆x = = Fo∆ ⋅ Bi∆ 2 ρc∆x ρc ∆x λ
( ( ( ) t Ni +1) = t Ni ) (1 − 2 Fo∆ ⋅ Bi∆ − 2 Fo∆ ) + 2 Fo∆ t Ni −1 + 2 Fo∆ ⋅ Bi∆ t f
∆τ
从第二式得出
∂t ∂τ
=
n ,i
( ( t ni ) − t ni −1)
∆τ
+ O ( ∆τ ) ≈
( ( t ni ) − t ni −1)
∆τ
difference。 向后差分 back difference。
∂t 二级数相减： 二级数相减： ∂τ
( ( ( ( t ni +1) − t ni −1) t ni +1) − t ni −1) 2 = + O(∆τ ) ≈ 2∆τ 2∆τ
n ,i
( 显式格式
explicit finite difference scheme )
如扩散项用（ +1)时层的值来表示 如扩散项用（i+1)时层的值来表示
( ( ( ( ( tni +1) − tni ) tni++1) − 2tni +1) + tni−+1) 1 =a 1 ∆τ ∆x 2
（隐式格式 implicit finite difference scheme） ）
例题4-3厚 的无限大平板受对称的冷却，初始温度t。 例题 厚2δ=0.06 m 的无限大平板受对称的冷却，初始温度 。=100℃。在初始瞬间， ℃ 在初始瞬间， 平板突然被置于t。 ℃的流体中。 平板突然被置于 。=0℃的流体中。已知平板λ=40 W/(m2·K)，h=1 000W／(m2·K)， ， ／ 试用数值法求解其温度分布。 试用数值法求解其温度分布。取Fo∆=1。 。
例题 t0 = 200 0 C，∞ = 25 t 针肋如右图所示，碳钢=43.2W/(m.K) =43.2W/(m.K)， 针肋如右图所示，碳钢=43.2W/(m.K)， h = 120 W m 2 ⋅ K 求其温度分布及换热量。 求其温度分布及换热量。
0
C
解： P = π d = 0.03141 m
θ 3 −θ 4 ∆x λAc = hP θ 4 ∆x 2
m 2 ∆x 2 (2 + m 2 ∆x 2 )θ 4 = 2θ 3 θ3 − θ 4 = θ4 ⇒ 2
(2 + m (2 + m (2 + m
2 2 2
∆x2 θ ∆x2 ∆x2
) )θ )θ
2 3 4
= θ1 + θ3 =θ
2
+θ4ຫໍສະໝຸດ = 2θ 3得 Φ = λ Ac
θ1 − θ 2
∆x
2 2
=λ
π d 2 θ1 − θ 2
4 ⋅ ∆x
如采用粗网格, ∆x = 15 则类似可以得 :
(2 + m (2 + m
∆x θ 2 = θ1 + θ 3
2
∆x
2
) )θ
3
= 2θ 2
三种情况的计算结果如下: 三种情况的计算结果如下: X 0 175 175 175 10 139.5 139.8 128.13
a∆τ ( i ) 2a∆τ ( = tn+1 + tni−)1 ) + 1 − 2 ( ∆x ∆x 2
− 2t + t ∆x
(i ) n 2
(i ) n −1
∂t ∂x ∂t ∂x
+
n ,i
i i tn +1 − tn = ∆x i i tn − tn −1 = ∆x
−
(i ) tn
误差
15 127.9
20
30
θ θ θ
热量计算： 热量计算：
119.7 113.4 120.13 113.8 114.29
Φ =15.06 W 精确解 Φ =11.94 W 四节点 21% Φ =10.52 W 三节点 30% 如取5 节点, 如取5 节点, 则Φ 的误差为 19%
§ 4-3 非稳态导热问题的数值解法
一维无限大平板非稳态导热显式格式离散方程组及稳定性分析 ①显式格式离散方程组 内节点
( ( ( ( tni +1) = Fo∆ ( tni+)1 + tni−)1 ) + (1 − 2 Fo∆ ) tni ) (1) (tn ＝t0，n=1,2, ,N-1) L
( ( ( ) 边界节点 t Ni +1) = t Ni ) (1 − 2 Fo∆ ⋅ Bi∆ − 2 Fo∆ ) + 2 Fo∆ t Ni −1 + 2 Fo∆ ⋅ Bi∆ t f ( ( (t−i1) = t2i ) )
Θn =
tn − t f t0 − t f
n n
η=
∑ ∆A (t
n =1 n
N
n
− t f )h − tf )× h
∑ ∆A × (t
n =1 n
N
=
0
∑ ∆A ×1× h
n =1 n
n =1 N
∑ ∆A Θ h
n n
N
=
∑ ∆A Θ
n =1 N
N
∑ ∆A
n =1
n
∆An = 2π rn ∆R, n = 2,3...N − 1
( t ni −1)
∂t ( = t ni ) − ∆τ ∂τ
∂t ∂τ =
∆τ 2 ∂ 2 t + 2 ∂τ 2 n ,i
−L
n ,i
( ( t ni +1) − t ni )
从第一式得出
∆τ 向前差分 forward difference
n ,i
( ( t ni +1) − t ni )
+ O(∆τ ) ≈
∆R ∆A1 = 2π r1 2 ∆R ∆AN = 2π rN 2
第四章 热传导问题数值解法
数值结果及分析
网格独立解： 网格独立解： N>=36
第四章 热传导问题数值解法
(m =
h H 3/ 2 ) λ AL
第四章 热传导问题数值解法
Example 4-5 短直肋效率的计算
例题4-5如图所示，一粗而短的肋片的三个表面与温度为 的流体换热 的流体换热， 例题 如图所示，一粗而短的肋片的三个表面与温度为tf的流体换热，且表面 如图所示 传热系均为h。试计算在下表所示的两种条件下肋片的效率， 传热系均为 。试计算在下表所示的两种条件下肋片的效率，并与一维分析解的 结果相比较。（判断肋片可以按一维问题处理的主要依据。 。（判断肋片可以按一维问题处理的主要依据 结果相比较。（判断肋片可以按一维问题处理的主要依据。 ）
θ4
d θ ∆x = 10 = m 2θ 2 dx θ1 节点2 节点2： + θ3 − 2θ 2 = m 2θ 2 2 ∆x 2 2 (2 + m ∆x )θ 2 = θ1 + θ3 同理得节点3 同理得节点3 (2 + m 2 ∆x 2 )θ = θ + θ 3 2 4
2
∆x
节点4 用热力学第一定律， 节点4 用热力学第一定律，导入的热量应等于对流散出 的热量，固有： 的热量，固有：
导热微分方程的建立
由：ρ c ∂t 1 ∂ ∂t 1 ∂ ∂t ∂ ∂t (λ r ) + 2 (λ ) + (λ ) + = ∂τ r ∂r ∂r r ∂φ ∂φ ∂z ∂z K
1 d r dr
(r
dt dr
)+
λ
=0
r = r1 , t = t0 ; dt r = r2 , =0 dr
第四章 热传导问题数值解法
题解 分析：取∆=x=0.01m，则 分析： ，
Bi∆ = 1000 × 0.01 = 0.25 λ 40 1 1 Fo∆ ≤ = ，取 Fo∆＝1的计算结果会振荡。 2(1 + Bi∆ ) 2.50 = h ∆x
计算结果汇总如下表所示。
讨论：从上表可以看出，从i=3时刻起出现了各点温度随时间作忽高忽低的波 动情况，并且波动幅度越来越大；某点温度越高反使其相继时刻的温度越低， (3) (4) t0 > t1(3) t0 < t1(4) 例如 ，但 ，这种现象是荒谬的，它违反了热力学第二 定律（意味着，在该时间间隔中从某一时刻起热量将自动地由低温点向高温 点传递）。数值计算中出现的这种计算结果忽高忽低的波动现象，数学上称 为不稳定性。因此，保证数值计算格式的稳定性是很重要的。
• 多非稳态项 • 扩散项的处理方法与前一样 数学描述 区域离散化 建立节点物理量的代数方程 设立迭代初场 求解代数方程组 解的分析 以一维为例： 以一维为例： 空间坐标 x 1∼N ∆x 空间步长 时间坐标 τ 1∼I ∆ τ 时间步长 （n,i）代表了时间空间区域中的一个 接点位置 t(i)n
在节点（ n,i- 对点（ 将温度函数 t 在节点（n,i+1）和（n,i-1）对点（n,i）作泰勒展开 ∆τ 2 ∂ 2t ∂t ( i +1) (i ) t n = t n + ∆τ + +L 2 ∂τ n ,i 2 ∂τ n ,i
R = R1 , Θ = 1;
dΘ =0 R = R2 , dR