罗大雷-导热问题的数值解法
4第四章导热问题的数值解法
基本概念：网格线、节点、步长、控制容 积
(m,n) N
n
y
y
x x
M m
（b）
（3）建立节点物理量的代数方程（离散方程） 节点上物理量的代数方程称离散方程。其过程 如下： • 首先划分各节点的类型；
2. 数值解(numerical method): 用某种方式把微分方程化为关 于各个离散点(节点)的代数方程,通过解代数方程获得问题近 似解的方法。
连续——离散(任意情况)
一、 数值解法的实质
对物理问题进行数值解法的基本思路可以概括为：把 原来在时间、空间坐标系中连续的物理量的场，如导热 物体的温度场等，用有限个离散点上的值的集合来代替， 通过求解按一定方法建立起来的关于这些值的代数方程， 来获得离散点上被求物理量的值。该方法称为数值解法。
（5） 求解代数方程组
如前图所示，除 m=1 的左边界上各节点的温度已知外， 其余(M-1)N个节点均需建立离散方程，共有(M-1)N个方 程，则构成一个封闭的代数方程组。实际工程问题代数方 程的个数在103-106数量级，只有利用现代计算机才能迅 速获得所需要的解。 1）常物性、无内热源（或具有均匀的内热源）的导热 代数方程一经建立，其中各项系数在整个求解过程中不再 变化——线性代数方程组；
这些离散点上被求物理量值的集合称为该物理量的 数值解。
二、 物理问题的数值求解过程
建立控制方程及定解条件 确定节点（区域离散化）
设立温度场的迭代初值
建立节点物理量的代数方程
求解代数方程
改进初场
是否收敛 否
是 解的分析
第4章 热传导问题的数值解法
w
A t
x
Δy 1 tm1,n
Δx
tm,n
e
A
t x
Δy 1 tm1,n
Δx
tm,n
s
A
t y
Δx 1 tm,n1
Δy
tm,n
n
A
t y
Δx
1
tm,n1 Δy
tm,n
热平衡时：
w e s n 0
如果 Δx Δy ：
1 tm,n 4 (tm1,n tm1,n
Δx 1!
t
x m,n
Δx 2 2t 2! x2
m,n
tm1,n
tm,n
Δx 1!
t x
m,n
Δx 2 2!
2t x 2
m,n
tm1,n tm,n
Δx 1!
t
x m,n
Δx 2 2t 2! x2
m,n
两式相加：
tm1,n tm1,n 2tm,n
Δx
2 2t x 2
③ 内部角点
2
hx
3t
m,n
2
tm1,n
tm,n1
tm1,n
tm,n1
3x 2 m ,n 2
2hx
tf
hΔx Δx
1 h BiΔ ——网格 Bi 数
4.3.2 处理不规则区域的阶梯型逼近法
如果边界为曲线，或倾斜边界， 用阶梯形折线模拟真实边界。
4.3.3 求解代数方程的迭代法
求解代数方程组： 直接解法——高斯消元法，矩阵求逆； 迭代法——逐次逼近法。
) 1
2FoΔBiΔt
f
式中，网格毕渥数 BiΔ hΔx ，网格傅里叶数 FoΔ aΔ Δx2
No.08 1013 4 导热问题的数值解法
Φ +Φ +Φ +Φ = 0 下 左 右 上
(m,n+1)
∆y
(m-1,n) (m, n) (m+1,n)
(m-1,n)
(m+1,n)
∆y
(m,n-1)
(m,n-1)
y o
∆x
∆x
15
x
以二维、稳态、有内热源的导热问题为例，此时： 以二维、稳态、有内热源的导热问题为例，此时：
Φi = Φ +Φ +Φ ＋ 右 +Φg = 0 下 左 Φ 上
其节点方程为：
ti +1, j − 2ti , j + ti −1, j
∆x 2
& ti , j +1 − 2ti , j + ti , j −1 Φv ,i , j + + =0 2 ∆y λ
13
(2)热平衡法 (控制容积平衡法)
基本思想： 基本思想 对每个有限大小的控制容积应用能量守恒，从而获得
温度场的代数方程组。它从基本物理现象和基本定律出发，不必事 先建立控制方程，依据能量守恒和Fourier导热定律即可。 能量守恒： 流入控制体的总热流量＋控制体内热源生成热= 控制体内能的增量 流入控制体的总热流量＋控制体内热源生成热
即：
Φ i + Φ g = Φ st
单位： 单位 [ W ]
注意：
横坐 标节 点编 号 N
(m,n)
n
（m，n）
∆y
y x
纵坐 标节 点编 号
∆x
m
M
8
N
(m,n)
网格( 网格 grid )划分 划分
网格划分方法: 网格划分方法 方法1： 方法 ： 先确定节点， 先确定节点，后定界面 方法2： 方法 ： 先确定界面，后定节点 先确定界面， 均分网格: 均分网格
第4章 导热问题的数值解法(含控制容积)
ti−1, j
二方程相加， 二方程相加，得：
ti+1, j − 2ti, j +ti−1, j ∂2t 2 = + 0(∆x2 ) ∂x ∆x2 i, j ti, j+1 − 2ti, j +ti, j−1 ∂2t 2 = + 0(∆y2 ) ∂y ∆y2 i, j
ti, j −ti−1, j ∂2t ∆x ∂3t ∆x2 ∂t + 2 − 3 +...... = ∂x 2 ∂x ∆x ! ∂x i, j i, j ! i, j 3 = ti, j −ti−1, j ∆x + 0(∆x)
2012-5-9 2
§4-1 导热问题数值求解的基本思想 及内节点离散方程的建立
2012-5-9
3
一、数值解法的基本思想 用导热问题所涉及的空间和时间区域内有限个离散 称为节点 节点) 点(称为节点)的温度近似值来代替物体内实际连续的温 度分布， 度分布 ， 将连续温度分布函数的求解问题转化为各节 点温度值的求解问题， 点温度值的求解问题 ， 将导热微分方程的求解问题转 化为节点温度代数方程的求解问题。 化为节点温度代数方程的求解问题。 数值解法的基本内容与步骤: 数值解法的基本内容与步骤:
第4章 导热问题的数值解法共30页
若取上面式右边的前三项，并将式①和式③相加 移项整理即得二阶导数的中心差分：
2t tm 1 ,n2 tm ,ntm 1 ,no( x2)
x2m ,n
x2
截断误差
同样可得：
未明确写出的级数余项中
的Δx的最低阶数为2
2t tm ,n 12tm ,ntm ,n 1o( y2)
y2m ,n
y2
28.05.2020 - 8 -
(3) 实验法: 是传热学的基本研究方法，a 适应性不好； b 费用昂贵
数值解法：有限差分法（finite-difference）、 有限元法（finite-element） 、 边界元法（boundary- element）、 分子动力学模拟（MD）
28.05.2020 - 2 -
第4章 导热问题的数值解法——§4-1 导热问题数值求解的基本思
2 例题条件
y
h3t f
W
t0
t2 x 2
t2 y 2
0
x 0, t t0
x H,
t x
h2 (t
tf)
h2t f
y 0,
t y
h1 (t
tf)
yW ,
t y
h3 (t
tf)
h1t f
Hx
二维矩形域内稳态无内热源，
常物性的导热问题
28.05.2020 - 4 -
第4章 导热问题的数值解法———§4-1 导热问题数值求解的基本思想
第4章 导热问题的数值解法———§4-1 导热问题数值求解的基本思想
以二维、稳态、有内热源的导热问题为例 此时：
Φ 上 Φ 下 Φ 左 ＋ Φ 右 Φ v 0 左Ad dxtyd dxt
导热问题的数值解法
第四章 导热问题的数值解法
4
§4-1 导热问题数值求解的基本思想 及内部节点离散方程的建立
1物 理 问 题 的 数 值 求 解 过 程
建立控制方程及定解条件
确定节点（区域离散化）
设立温度场的迭代初值
建立节点物理量的代数方程
求解代数方程
改进初场
是否收敛 否
是 解的分析
第四章 导热问题的数值解法
5
2 例题条件
tm1,n
tm,n
t x
x
m,n
2t x2
m,n
x2 2!
3t x3
m,n
x3 3!
第四章 导热问题的数值解法
9
若取上面式右边的前三项，并将式①和式③相加 移项整理即得二阶导数的中心差分：
2t
tm1,n 2tm,n tm1,n o(x2 )
x2 m,n
x 2
同样可得：
2t tm,n1 2tm,n tm,n1 o(y2 )
代数方程组的求解方法：直接解法、迭代解法
第四章 导热问题的数值解法
24
直接解法：通过有限次运算获得代数方程精确解; 矩阵求逆、高斯消元法
缺点：所需内存较大、方程数目多时不便、不适用于非线性 问题（若物性为温度的函数，节点温度差分方程中的系数不 再是常数，而是温度的函数。这些系数在计算过程中要相应 地不断更新）
25
例如：根据第 k 次迭代的数值 可以求得节点温度：
t1(k)、t2(k)....tn(k)
t1(k1) a11t1(k) a12t2(k) ...... a1ntn(k) b1(k)
在计算后面的节点温度时应按下式（采用最新值）
t2(k1) a21t1(k1) a22t2(k) ...... a2ntn(k) b2(k) t3(k1) a31t1(k1) a32t2(k1) ...... a3ntn(k) b3(k) ....................................................... tn(k1) an1t1(k 1) an2t2(k 1) ...... ann1tn(k11) anntn(k ) bn(k )
导热问题的数值解
a 1 f 2 x 2
隐式差分格式是无条件稳定的
边界节点的稳定性条件：
1 f 2b 2
其中：b h x
谢谢！
数值求解过程
3.1 区域的离散化
3.2 建立节点方程
控制方程
（无内热源）
用泰勒级数法导得温度节点方程式：
法二：元体热平衡法
本质：傅里叶导热定律和能量守恒定律的体现。
QE QW QN QS qV xy 0
3.3 设立迭代初场
代数方程组的求解方法有直接解法与迭代解法， 传热问题的有限差分法中主要采用迭代法。采 用迭代法求解时，需对被求的温度场预先设定 一个解，这个解称为初场，并在求解过程中不 断改进。
ti
k 1
a a k k k (ti 1 ti 1 ) (1 2 2 )ti 2 x x
a 引入f 2 x
得 : ti
k 1
f (ti 1 ti 1 ) (1 2 f )ti
k k
k
为了加快计算的进程而调整x和 的大小 时，必须遵守使式中 ti k 的系数大于或至少 等于零。即
有限差分形式
显示格式：
计算工作量小，但受时间及空间步长的限制。 隐式格式：
ft
p 1 i 1
(1 2t )t
p 1 i
ft
p 1 i 1
t
பைடு நூலகம்
p i
不受时间及空间的步长影响，计算工作量大。
5. 差分格式的稳定性
稳定性的表述：如果初始条件和边界条件有微 小的变化，最后的解是否也只 有微小的变化，由此来判断解 的稳定与否。 傅里叶级数法
第五章--导热问题的数值解
5-2 稳态导热的数值分析
把能量守恒关系应用于每个元体，在稳态导热的情况下，从所有 相邻的元体导入的热量和该元体本身的发热量之代数和应为零。 这样，对于图5-3 所示的内部元体，有以下的关系：
QE QW QN QS qV xy 0
( 5-2-13 )
在根据傅里叶定律计算相邻两个元体的导热量时，假设两个节点
间的温度是线性分布的。这实际上也是采用了以一阶差商近似偏
导数的概念。例如：
QE
ti1, j ti, j
x
y
，QW
ti1, j ti, j
x
y
代入式(5-2-13)并整理，得内部元体的节点方程为
ti1, j
2ti, j ti1, j (x)2
ti,
j 1
2ti, j (y)2
ti,
2ti, j ti1, j (x)2
ti, j1
2ti, j (y)2
ti,
j 1
qV ,i,
j
0
(5-2-4)
5-2 稳态导热的数值分析
如果采用正方形的网格，即△x＝△y 式(5-2-4 ) 简化为
，且无内热源(qV＝0)
，则
1
ti, j 4 ti1, j ti1, j ti, j1 ti, j1
5-1 导数的有限差分近似表达式
有限差分的数学基础是用差商代替微商，即用有限差分代替导数。 若f(x)是连续函数，则它的导数定义为
dt
f (x x) f (x) f
lim
lim
dx x0
x
x0 x
(5-1-1)
在这里，df /dx称为微商(导数)，△f/△x 称为有限差商。微商是有 限差商当△x趋于零时的极限。在△x没有达到零以前，△f/△x只 是df/dx的近似，而两者的差值就是用差商代替微商的偏差。
传热学—第4章 热传导问题的数值解法
⎧a11t1 + a12 t2 + a13t3 = b1 ⎪ ⎨a21t1 + a22 t2 + a23t3 = b2 ⎪a t + a t + a t = b 33 3 3 ⎩ 31 1 32 2
假定初场
⎧ (1) ⎪t1 = ⎪ ⎪ Jacobi ⎨t(1) = 2 ⎪ ⎪ (1) ⎪t3 = ⎩
4.1.1 4 1 1 基本思想 把原来在时间、空间坐标系中连续的物理量的场， 用有限个离散点上的值的集合来代替，通过求解按 定方 建 起来 关 值 代数方程 来获 一定方法建立起来的关于这些值的代数方程，来获 得离散点上被求物理量的值。 这些离散点上被求物理量值的集合称为该物理量 的数值解。
4.1.1 基本思想
λ Δy
Δx = Δy 时： tm −1,n
+ tm+1,n + tm,n+1 + tm,n−1 − 4tm,n = 0
tm ,n
1 = ( tm−1,n + tm+1,n + tm,n+1 + tm ,n−1 ) 4
与Taylor级数法相比，热平衡法物理意义明显。
4.3.1 边界节点离散方程的建立
4-2 内部节点离散方程的建立
4.2.1 4 2 1 Taylor级数展开法
4-2 内部节点离散方程的建立 内部节点离散方程的建
∂ 2t ∂x 2
=
m ,n
tm+1 n − 2tm ,n + tm −1 n 1, 1, Δx 2
控制方程
∂ 2t ∂ 2t + =0 ∂x 2 ∂y 2
∂ 2t ∂y 2
传热学第五章导热问题数值解法
2 h∆y
=
h∆y 2 2 + λ t m ,n+1 + t m ,n−1 + 2(t m−1,n + Bi ⋅ t f ) 2(2 + Bi )
9
λ
tf
3）方程组的求解 ） 对应每个未知量（一个节点温度） 一条方程 一条方程（ 对应每个未知量（一个节点温度）→一条方程（一 个节点方程） 方程组有唯一解 个节点方程）→方程组有唯一解 解法：（ ：（1） 解法：（ ）矩阵法 （2）迭代法：高斯 赛德尔迭代 ）迭代法：高斯-赛德尔迭代
λ∆y
t m−1,n − t m ,n ∆x ∆x t m ,n−1 − t m ,n +λ + ∆y ⋅ h(t f − t m ,n ) = 0 2 ∆y ∆y ∆x t m ,n+1 − t m ,n +λ ∆y 2
如取正方形网络
∆x = ∆y
上式简化为： 上式简化为：
t m ,n =
2t m−1,n + t m ,n+1 + t m ,n−1 +
对于非稳态导热问题，除了空间上进行网格划分外， 对于非稳态导热问题，除了空间上进行网格划分外， 还要把时间分割成许多间隔。 还要把时间分割成许多间隔。
4
2）有限元法 ）
把整个求解域离散成为有限个子域， 把整个求解域离散成为有限个子域，每一子域内运 用变分法， 用变分法，即利用与原问题中微分方程相等价的变 分原理来进行推导，从而使原问题的微分方程组退 分原理来进行推导， 化到代数联立方程组，得到数值解。 化到代数联立方程组，得到数值解。 有限元法和差分法都是常用的数值计算方法， 有限元法和差分法都是常用的数值计算方法，差分 法计算模型对于不规则的几何形状难以应用。 法计算模型对于不规则的几何形状难以应用。有限 元法能够很好地适应复杂的几何形状、 元法能够很好地适应复杂的几何形状、复杂的材料 特性和复杂的边界条件。 特性和复杂的边界条件。
第4章-热传导问题的数值解法(1)
y 2
y 2
m,n
将差分表达式代入控制方程
2t x 2
2t y 2
0
得:
tm1,n
tm1,n x2
2tm,n
tm,n1
tm,n1 y 2
2tm,n
0
如果 x y 则有:
tm,n
1 4
tm1,n tm1,n tm,n1 tm,n1
第4章 热传导问题的数值解法
4.1 导热问题数值求解的基本思想 4.2 内节点离散方程的建立方法 4.3 边界节点离散方程的建立及代数方程的求解 4.4 非稳态导热问题的数值解法
4.1 导热问题数值求解的基本思想
4.1.1 基本思想 4.1.2 导热问题数值求解的基本步骤
返回
4.1.1 基本思想
2、边界上的外部角点
边界节点D代表的区域为1/4个普通元体大小 的面积。对该外部节点元体应用能量平衡
y
2
tm1,n tm,n x
x tm,n1 tm,n 2 y
xy 4
.
m,n
x y 2
qw
0
如 x y ，则有：
在均分网格中，一、二阶导数常见的差分表达式如下表所示：
返回
4.2.2 热平衡法（热力学第一定律）
n
热平衡法不是在控制方程的基础上进行离
散，而是直接对元体应用热力学第一定律
和傅里叶定律，从而得到该节点温度的离 w
e
散方程。
二维稳态常导热系数无内热源的稳态导热
问题，对元体(m,n)列出能量守恒方程：
在空间-时间坐标系中对所研究的空间区域 和时间区域进行离散
导热问题的数值解法一
T5 = Tr
2 = 0 = λ 1 W / m ⋅K) (
一维稳态常物性导热问题
有限容积积分(内点)
∆x x i-1 i-1/2 i i+1/2 i+1
∫
xi +1/ 2
xi −1/ 2
d dT λ dx dx
dT 0⇒λ dx = dx
对所求变量给出预估值，然后根据离 散方程对其不断改进，直至求得收敛解。
一维稳态常物性导热问题
Jacobi迭代法： 任意点上未知量的更新都用上一轮迭代中 所获得的值来计算
1 (I ) (I ) Ti aiTi + c T = + 1 i i −1 bi
( I +ቤተ መጻሕፍቲ ባይዱ)
(
)
Gauss-Seidel迭代法： 每一步计算取邻点的最新值
界面导热系数的计算方式
λ= i +1/ 2 λ= i −1/ 2
1 ( λi + λi +1 ) 2 1 ( λi + λi −1 ) 2
算术平均
哪种好？
∆x x i-1 i i-1/2 i+1/2 i+1
λi +1/ 2 =
2 1
λi λi −1/ 2 =
1
+
1
λi +1
1
2
调和平均
λi
+
λi −1
Ti +1 − Ti ⇒ = 1 Ti +1 − Ti 2 ⇒λ = i +1/ 2 1 1 1 1 + + 2λi +1 2λi λi +1 λi
11导热问题的数值解法
1 ti ti1 ti1 2
?
i1 i 5
解：
2 1000 W m ℃ L 0 .5 m A 0 .01 m
tW1 100 ℃
i 1 , 2 , , 4 , 5
tW2 500 ℃
x 0 . 1 m
i 1
tw t 1 1 Q A left 1 x2
1 k k 1 k 1 k t t t t t i , j i 1 , j i 1 , j i , j 1 i , j 1 4
k 0 , 1 , 2 , 3 ,
k 1 k max t t i ,j i ,j
稳态导热的有限差分方法 以常物性、无内热源矩形区
生成网格 将求解区域划 分为有限个网 格单元
j n
y
域的二维稳态导热为例
差分方程的建立－差商代替微商
t x , y i , j ti, j
j y
i, j1 i1 , j
i 1 , j i, j 1
y
t i x , j y
tW1 100 ℃
i 1 , 2 , , 4 , 5
tW2 500 ℃
x 0 . 1 m
1 t1 t2 2tw1 i 1 3 1 ti ti1 ti1 i 2 , , 4 2 1 t5 t4 2 tw2 i 5 3
t t 1 2
t 3
解
ti,kj1 ti,kj max k ti, j
特点：计算第 k+1次的值时， 全部使用第k次 迭代的值，收 敛速度比较慢
差分方程组的求解 高斯－赛德尔迭代：
t
第四章导热问题的数值求解
第四章导热问题的数值求解随着计算机的普及应用和性能的不断改善，以及相关的数值计算方法的发展和应用程序的开发，传热学数值计算方法作为数值求解传热问题的有效工具也得到了相应的发展，利用计算机求解传热学问题愈来愈受到人们的普遍重视，而且在计算复杂传热问题中显示出它的优越性，因而成为传热学的一个重要的分支。
数值传热的相关内容也很自然地成为工程类学生学习传热学课程的不可缺少的部分。
为了使学生能简要地掌握传热学数值计算的基本方法，在这里我们以导热问题为例对传热学数值计算方法做一个简单的介绍。
4-1导热问题数值解概述在第二章和第三章中我们对较为简单的导热问题，如一维、二维简单几何形状和边界条件的稳态导热和非稳态导热、以及一些特殊导热问题，象通过肋片的导热和忽略内阻的集总导热系统，进行了分析求解。
然而对于一些更为复杂的导热问题，如复杂的几何形状和边界条件以及物性变化较大的情况，分析求解往往很复杂或者根本不可能。
此时求解问题的唯一途径是利用数值分析的办法获得数值解。
数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再求解。
这里要介绍的是后一种方法。
如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。
作为一本入门的教材，这里仅向读者简要地介绍用有限差分析方法从微分方程确立代数方程的处理过程。
有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量（如温度、压力、速度和热流等），用有限个离散点上的数值集合来近似表示。
有限差分的数学基础是用差商代替微商（导数），而几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。
在图4-1中可以看出有限差分表示的温度场与真实温度场的区别。
图中用T0、T1、T2…表示连续的温度场T；Δx为步长，它将区域的x方向划分为有限个数的区域，Δx0、Δx1、Δx2…，它们可以相等，也可以不相等。
当Δx相等时，T1处的真实变化率a可以用平均变化率b、c或d来表示，其中b、c和d分别表示三种不同差分格式下的温度随时间的变化率，即：T图4-1温度场的有限差分表示b 为向后差分格式xx x T x T dx dT ∆∆--≈)()(111； c 为向前差分格式xx T x T dx dT ∆-∆≈)()(111＋；d 为中心差分格式xx x T x x T dx dT ∆∆--∆+≈)()(111。
传热学-第四章-热传导问题的数值解法
36
n=N
w (m-1,n)
(m,n+1)
n e (m+1,n)
s (m,n)
(m,n-1)
y
n=1
m=1
m
x
m=M
Monday, March 30, 2020
37
而获得离散点上被求物理量的值；并称之为数值解 a. 在很大程度上弥补了分析法的缺点；适应性强，特别对于复杂问题更显其 优越性； b. 与实验法相比成本低 数值解法： 有限差分法（finite-difference）、
有限元法（finite-element） 、 边界元法（boundary- element）、 分子动力学模拟（MD）
局限性： 简单几何形状及边界条件
稳态问题：直接积分法 非稳态问题：分离变量法 解析解（analytical solution）
工程实际中面临的大部分问题几何形状和边界条件要复杂的多，由于数学上 的困难还不能给出解析解，导致目前解析解只能作为某些简单问题的参照依 据，不能解决实际问题。
Monday, March 30, 2020
y
n=1
m=1
m
x
m=M
Monday, March 30, 2020
16
1.边界节点离散方程的建立： (1) 平直边界上的节点
qw
(m,n+1)
(m-1,n)
(m,n)
(m,n-1)
qw
y x
Monday, March 30, 2020
17
(2) 内部角点
qw
(m-1,n)
(m,n+1) (m,n)
4
(2) 实验法: 是传热学的基本研究方法： a 偏向于机理研究； b.受场地，燃料动力源等因素的影响，无法完全复现研究对象，具有时间、
4 导热问题的数值解(2014)
t(k)
i
+
Fo
t + t (k +1)
i −1
(k +1) i+1
( ) ( ) ( ) t(k+1)
i
=
1 1+ 2Fo
t(k)
i
+
Fo 1+ 2Fo
t + t (k +1)
i −1
(k +1) i +1
边界节点的离散化方程（隐式）
ρc
t (k +1)
i
−
t(k)
i
∆τ
∆x 2
=
−λ
ti(k
边界节点
∆U = ∆ΦL + ∆ΦR
∆Φ L ⇒
i −1
∆U ⇐ ∆Φ R = −h(ti(k ) − t∞ )
i
x
边界节点的离散化方程（显式）
ρc
t (k +1)
i
− ti(k )
∆τ
∆x 2
=
−λ
t(k)
i
−
t(k)
i −1
∆x
+
h(t∞
− ti(k ) )
t (k +1)
i
=
t(k)
2
t (k )
i −1
+
t (k )
i +1
−
2ti(k )
( ) t(k+1)
( ) ( ) i
=
1 −
2a∆τ ∆x 2
ti(k )
+
a∆τ ∆x 2
t + t (k ) (k )
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
导热问题数值解法初次研究对物理物体的数值求解的基本思想可以概括为：把原有的时间、空间坐标系中连续的物理量的场，如导热问题的温度场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上的值。
这些离散点上的被求解物理量的值的集合称为该物理量的数值解。
物理模型在四个输气的管道中间有一个各边长为10厘米的薄铁片，求导热其达到稳态后，这块铁片的温度分布。
四个输气管道里的气体温度是恒值分别为100℃、200℃、500℃、1000℃。
因此，可以看成是二维矩形域内的稳态、无内热源、常物性的导热问题。
建立数学模型描写物理问题的微分方程称为控制方程，导热微分方程为：22220t t xy∂∂+=∂∂ （1）其四个边界分别为第一类边界条件，1234t =1005002001000===℃、t ℃、t ℃、t ℃。
区域离散化用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。
相邻两节点的距离称为步长，记为x ∆、y ∆。
本模型x 、y 方向是各自均分的，各自为100个子区域。
节点的位置以该节点在两个方向上的标号m 、n 来表示。
每一个节点都可以看成是以它为中心的一个小区域的代表，由相邻的两节点连接的中垂线构成。
为叙述方便，我们把节点所代表的小区域称为元体。
数学模型离散化它的建立是数值求解过程中的重要环节，主要有泰勒级数展开法及热平衡法两种，取节点（m ，n ）及其临点为例。
泰勒级数展开法以节点（m ，n ）处的二阶偏导数为例用这种方法来导出其差分表达式。
对节点(1,)m n +及(1,)m n -分别写出函数t 对（m ，n ）的泰勒级数展开式：2233441,,,,,,2342624m n m n m nm nm nm nt x t x t xtt t xx xxx +∂∆∂∆∂∆∂=+∆++++∂∂∂∂ （2）2233441,,,,,,2342624m n m n m n m nm nm nt x t x t xtt t xxxxx-∂∆∂∆∂∆∂=-∆+-++∂∂∂∂ （3）将式（2）、（3）相加得24421,1,,,,24212m n m n m n m nm nt xtt t t x xx+-∂∆∂+=+∆++∂∂ （4） 将（4）式改写成2,2m n t x∂∂的表达式，有21,,1,2,222()m n m n m n m nt t t t O x xx+--+∂=+∆∂∆ （5）这是用三个离散点上的值来计算二阶导数2,2m nt x∂∂的严格的表达式，其中符号2()O x ∆表示未明确写出的级数余项中x ∆的最低阶数为2。
在进行数值计算时，我们希望得出用三个相邻节点上的值表示的二阶导数的近似表达式。
于是，略去式(5)中的余项后，得21,,1,,222m n m n m n m n t t t t xx+--+∂=∂∆ （6-1）这就是二阶导数的差分表达式，称为中心差分。
同理可有2,1,,1,222m n m n mn m n t t t t xx+--+∂=∂∆ （6-2）将式（6-1）、（6-2）代入式（1）得离散方程1,,1,,1,,122220m n mn m n mn mn mn t t t t t t xy+-+--+-++=∆∆ （7）如果，x y ∆=∆，则式（7）即为下式，,1,1,,1,11()4m n m n m n mn mn t t t t t +-+-=+++ （8）热平衡法采用这种方法时，对每个节点所代表的元体可用傅里叶导热定律直接写出其能量守恒形式的表达式。
此时把节点看成是元体的代表。
通过元体的界面所传导的热流量可以对有关的两个节点应用傅里叶定律写出。
例如，从节点(1,)m n -通过界面w （节点(1,)m n -西面的界面）传导到节点(,)m n 的热流量可写为1,,m n m nw t tk yxφ--=∆∆ （9-1）类似的可以写出通过其他三个界面e 、n 、s 而传到给节点的(,)m n 的热量。
1,,m n m ne t tk yxφ+-=∆∆ （9-2）,1,m n m nn t tk xy φ+-=∆∆ （9-3）,1,m n m ns t tk xyφ--=∆∆ （9-4）对于所研究的问题，元体(,)m n 的能量守恒方程为0w e n s φφφφ+++= （10）将（9-1）、（9-2）、（9-3）、（9-4）分别代入得1,,1,,,1,,1,0m n m nm n m nm n m nm n m nt t t t t t t t k yk yk xk xxxyy-++-----∆+∆+∆+∆=∆∆∆∆ （11）（11）式中都已热量导入元体(,)m n 的方向为正，对其进一步两边同除以k x y ∆∆可得出式1,,1,,1,,122220m n mn m n mn mn mn t t t t t t xy+-+--+-++=∆∆即为式（7）。
如果，x y ∆=∆，则式（7）即为下式，,1,1,,1,11()4m n m n m n mn mn t t t t t +-+-=+++即为式（8）。
有上述推导过程可见，用热平衡法导出式（11）的思路和过程与建立导热微分方程的思路和过程完全一致，所不同的只是建立导热微分方程所讨论的是一个微元体，而此处的为一个有限大小的元体。
建立迭代初场和给出边界条件四个输气管道里的气体温度是恒值分别为100℃、200℃、500℃、1000℃，因此可设四条边的温度分别为1234t =1005002001000===℃、t ℃、t ℃、t ℃。
此次大作业采用迭代法，因此需要对被求的温度场预先假定一个值，设为四条边最低温度100=0t ℃，此值在迭代中不断更新。
求解离散化方程采用FORTRAN 语言求解离散化方程，原程序如下： program Heattransfer_Problem_kk12000 implicit noneinteger:: m_x=101 !x 方向的网格点数!integer:: n_y=101 !y方向的网格点数!integer:: k_k=12000 !迭代次数!integer m !x方向的网格代号!integer n !y方向的网格代号!integer k !迭代次数!real,allocatable:: t(:,:) !声明可变大小的数组!!定义初值函数!real,external:: fallocate(t(m_x,n_y))!赋初值（第一时间层）!do m=1,m_x !x方向循环!do n=1,n_y !y方向循环!t(m,n)=f(m,n)end doend do!计算温度场的结果!do k=1,k_k !迭代次数增加!do m=2,m_x-1 !x方向循环!do n=2,n_y-1 !y方向循环!!计算温度场的函数，输入一点周围四点温度值，返回这一点下一迭代层的值!t(m,n)=(t(m+1,n)+t(m-1,n)+t(m,n+1)+t(m,n-1))/4end doend do!if(max(abs(t(m,n,k)-t(m,n,k-1)),0.1*0.1*0.1)<=0.1*0.1*0.1) exitend do!输出温度场结果，结果的格式符合Tecplot的格式!open(unit=10,file='temperature_distributing.txt')write(10,*) '"m","n","k","t"'do n=1,n_y !y方向循环!do m=1,m_x !x方向循环!write(10,*)m,n,k,t(m,n) !输出t(m,n,k),m,n,k !end doend do!主程序结束!stopend!定义初值函数体!function f(x,y)implicit nonereal:: finteger:: xinteger:: yreal,parameter :: t1=100real,parameter :: t2=500real,parameter :: t3=200real,parameter :: t4=1000!给左边赋t1!if (x==1 .and. y>1 .and. y<101) thenf=t1!给右边赋t3!else if (x==101 .and. y>1 .and. y<101) thenf=t3!给下边赋t2!else if (y==1) thenf=t2!给上边赋t4!else if (y==101) thenf=t4!给内部点赋假定值t1!elsef=t1end ifreturnend解的分析温度场分布见“导热问题求解结果.txt”文件。
在程序中，分别令迭代次数为10、25、50、100、200、500、1000、1500、2000、5000、7500、10000、15000、12000、11000，从结果可以看出，迭代100步时，边界的温度影响传到了全场，迭代12000步收敛了。
取m=50,n=50的节点温度为例，迭代次数为10、25、50、100、200、500、1000、1500、2000、5000、7500、10000、15000、12000、11000，温度分别为100.0000℃(10)、 100.0000℃(25)、 100.0000℃(50)、 100.0044℃(100)、 100.5426℃(200)、 131.3818℃(500)、 235.3835℃(1000)、 315.1268℃(1500)、 365.4749℃(2000)、440.8724℃(5000)、 444.0441℃(7500)、 444.9626℃(10000)、 444.9761℃(15000)、444.9761℃(12000)、 444.9761℃(11000)。
此次求解具有良好的收敛性，温度场分布符合物理规律。