1 傅立叶定律
傅立叶定律是导热理论的基础。其向量表达式为:
q g r a d T λ=-? (2-1)
式中:q —热流密度，是向量，2
/()Kcal m h ；gradT —温度梯度，是向量，℃/m ；λ—导热系数，又称热导率，
/()Kcal mh C ； 式中的负号表示q 的方向始终与gradT 相反。
2 导热系数（thermal conductivity ）及其影响因素
导热系数λ（/()Kcal mh C ）是一个比例常数，在数值上等于每小时每平方米面积上，当物体内温度梯度为1℃/m 时的导热量。 导热系数是指在稳定传热条件下，1m 厚的材料，两侧表面的温差为1度（K ，°C ）,在1秒内，通过1平方米面积传递的热量，用λ表示，单位为瓦/米·度，w/m·k （W/m·K,此处的K 可用℃代替）。导热系数为温度梯度1℃/m ，单位时间通过每平方米等温面的热传导热流量。单位是：W/(m·K)。 3．热传导微分方程推导 ? 在t 时刻w 界面的温度梯度为x
T
??
在t 时刻e 界面的温度梯度为dx x T x T dx x x T
x T 2
2??+??=????+?? 单位时间内六面体在x 方向流入的热流量为：dydz x
T
??-λ； 单位时间内六面体在x 方向流出的热流量为：
dydz dx x T x T ??
??????+??-22λ；
单位时间内六面体在x 方向流入的净热量为：dxdydz x
T
22??λ 图3-1 微分单元体各面上进出流量示意图
同理，单位时间内六面体在y 方向流入的净热量为：dxdydz y
T
22??λ； 单位时间内六面体在y 方向流入的净热量为：dxdydz z T 22??λ； 单位时间内流入六面体的总热量为：dxdydz z T y T x
T ???
?????+??+??222222λ (3-1) 六面体内介质的质量为:dxdydz ρ。 单位时间六面体内热量的变化量（增加）为:Cdxdydz t
T
ρ?? 根据热量守恒定律:
Cdxdydz t T
dxdydz z T y T x
T ρλ??=????????+??+??222222， C t T z T y T x T ρλ??=????????+??+??222222，
t T
z T y T x T C ??=
????????+??+??222222ρλ， t T z T y T x
T a ??=
????????+??+??222222， C a ρλ
= α称为热扩散率或热扩散系数（thermal diffusivity ），单位m 2/s.
λ:导热系数，单位W/(m·K)； ρ：密度，单位kg/m 3 c ：热容，单位J/(kg·K). 思考：如果单元体内有热源：单位体积单位时间的散热量是q 方程怎么变？
4．岩石的热扩散率（导温系数) thermal diffusion coefficient ;thermal diffusivity; thermal degradation 岩石的热扩散率也叫或热扩散系数，表示岩石在加热或冷却时各部分温度趋于一致的能力。它反映岩石的热惯性特征，是一个综合性参数。热扩散率越大的岩石，热能传播温度趋于一致的速度越大，透入的深度也越大。 热扩散系数一般是根据岩石的导热系数(ramuda)、和密度(rou)的测量数据计算得到的。 如果单元体内有热源：单位体积单位时间的散热量是q 则在
Cdxdydz t T dxdydz z T y T x T ρλ??=??
??????+??+??222222 左边增加热源散发的热量qdxdydz ，
方程变为：
Cdxdydz t T qdxdydz dxdydz z T y T x
T ρλ??=+????????+??+??222222
t T
q z T y T x T a ??=
+??
??????+??+??ρ222222 如果是流体除了热传导外，还有流体流动带入和带出微元体的热量
在t 时刻
w 界面流体速度为U ，流体温度为T
单位时间流入微元体的流体质量为：udydz dm ρ=1 带入微元体的热量为：uTCdydz ρ e 界面流体速度为dx x u u ??+
,流体温度为dx x T T ??+ 单位时间流出微元体的流体质量为：dydz dx x u u dm ?????
?
??+
=ρ2 带出微元体的热量为:
Cdydz dx x T T dx x u u ??
??????+??????
??+
ρ dxdydz x
T
dx x u C Cdxdydz x T u TCdxdydz x u uTCdydz ????+??+??+ρρρ
ρ 如果不考虑x 方向速度变化，略去高阶微量，则e 界面带出微元体的热量为：Cdxdydz x
T
u
uTCdydz ??+ρρ 单位时间内在x 方向流入六面体的净热流量为：dxdydz x
T
uC
??-ρ； 同理， y 方向：dxdydz y T vC
??-ρ z 方向：dxdydz z
T wC ??-ρ 将上述三项加入到热传导方程（左边，进入微元体的热量）：
P
Δz
Δx
Δy
e
s x
n
b t
w
W
E
N
S
T
B
Cdxdydz t T dxdydz z T y T x
T ρλ??=????????+??+??222222
Cdxdydz t T dxdydz z T
wC dxdydz y T vC dxdydz x T uC dxdydz z T y T x
T ρρρρλ??=????????+??-+??-????????+??+??222222
???
?????+??+??+??=????????+??+??z T w y T v x T u t T z T y T x
T ρρρρλ222222 （能量方程）
2.2巷壁与风流间的对流换热
运动着的流体与所接触的固体壁面之间的热量传递过程称为对流换热，它是流体(液体或气体)由于宏观相对运动，从
某一区域迁移到温度不同的另一区域时引起热量传递的现象。固体壁面与流体之间存在温度差将产生对流换热，由于实际流体的粘性和壁面摩擦的共同影响，近壁流体分层流动，尤其与壁面直接接触的几何面上，总有一层很薄的流体粘附于表面，该层流体处于静止状态，所以热流通过表面层的传递只能依靠导热。显然，在流体发生热对流的同时，由于流体中温度分布的不均匀，也将伴随产生导热现象。因此，对流换热过程实际上是热对流和热传导的综合作用过程。 牛顿冷却公式
对流换热过程是一个受很多因素影响的复杂过程，如流体的流动状况、流体的物理性质、壁的形状和大小、表面粗糙度等。一般情况下对流换热的计算可采用牛顿冷却公式。根据对流换热定律，可以计算出从壁面某处进入通风风流的显热热流密度：
)
(T T q w s -=α (3) 式中：T w = 巷道壁面的温度； T= 巷道内风流的平均温度；
α= 巷道壁面的换热系数。在围岩与风流的热交换过程中，多半是井巷低温风流流经高温岩壁，井巷壁面向风流放
热，所以矿内常把上式中的对流换热系数α（
2
/()Kcal m h C ）称为巷壁与风流的换热系数，简称为放热系数。
圆形巷道（柱体）围岩与风流换热控制方程
地热通过围岩向风流的传热现象与围岩本身的热传导、巷道壁面向风流的对流换热以及壁面上的水分蒸发等因素有关。由于实际情况下围岩的散热是一个很复杂的过程，为了方便本论文的研究，对要研究的物理模型做了简化和假设：
1) 巷道为圆形、无限扩展，围岩岩石均质、各向同性； 2) 不考虑围岩壁面的热辐射作用。
根据上述假设，可得到描述考虑壁面水分蒸发时围岩与风流热质传递的数学方程，如式（3-1）：
020200001() (;0)(,) ()(,) (0)()() (0)t r R w a v w a r r T T T
a r r R t t r r r T r t T r r R T r t T t T
T T f L m m t r λασ===????=+?<<>??????
?=<≤?
??=≥????=-+-≥???
(3-1)
式中：R ——调热圈半径，m ；其他符号的意义同前章所述。
根据简化的数学模型，可将巷道围岩划分为一系列等间距 (R ?)的同心圆，取垂直于长轴的巷道断面角度为θ?，如图3-1所示。
第二章 物理现象的数学描述
控制微分方程：把控制传热、流体流动等有关过程的规律表达成数学形式；详细、完整的推导，阅读标准的教科书；数值解法
方程的形式和意义：我们这里所提到的所有方程都具有一个共同的形式；形式上的一致是构成一个通用解法的基础。
2．1 控制微分方程
2．1—1 微分方程的意义
守恒原理：各个微分方程各自代表着一定的守恒原理．每一个方程以一定的物理量作为它的因变量，方程本身则代表着那些影响该因变量的各个因素之间必定存在着的某种平衡．
这些微分方程的因变量通常具有“比”的性质，即以单位质量为基础来表示各因变量。 这种因变量的例子有：质量分量、速度(即单位质量的动量)以及比焓． 这一类微分方程的各项代表着以单位容积为基础的效果
例： 设想J 表示一个典型因变量Φ的流量密度．让我们考虑如图2．1所示的尺寸为dx 、dy 及dz 的控制容积．
x J (J 在x 方向的分量)：进入面积为dydz 的流量密度
dx x J J x x )/(??+： 离开与这个面相对曲面上的流量密度
通过该面的整个面积上流出的净流量是：dxdydz x J x )/(?? dxdydz ：所讨论区域的容积
用同样的方法考虑y 与z 方向的贡献，有：
divJ z
J y J x J z
y x =?+?+?=
单位容积流出的净流量 （2.1） 由于我们的数值方法是通过对一个控制容积进行平衡构成的(就象我们将要在后面看到的那样)，上述divJ 的表达方式对我们来说是特别有用的．
以单位容积为基础来表达一项的另一个例子是变化速率项
t ??/)(ρφ
如果Φ是某个“比”性质（单位质量）而ρ是密度，那么ρΦ就代表在单位容积内所包含的相应广延性质的大小．于是
T
T w
图3-1 巷道围岩内节点划分
R N
R
P 1
P NR
P I+1 P I P I-1 R I+
1 R I
R I-1 R 1
R 0
Δ
α
t ??/)(ρφ是单位容积内有关性质的变化率．
一个微分方程是这样一些项的组合；其中每一项代表一个以单位容积为基础的效应；而所有的项合在一起则反映着某种平衡或是守恒。
我们现在以几个标准的微分方程为例，并由此找到一个通用的形式．
2．1—2 化学组分的守恒
令m l 代表一种化学组分l 的质量分量．当存在有速度场u 时，可把m l 的守恒表示为： l
l l l R J m t
m =++??)div()
(u ρρ （2.2）
t m l ??/)(ρ: 单位容积内化学组分l 的质量变化率； ρu m l :对流流量密度(即一般化流场ρu 所携带的流量密度)．
J l :扩散流量密度，它通常是由m l 的梯度引起的． 两部分流量密度(对流与扩散)的散度构成微分方程的第二项． R l :单位容积的化学组分l 的生成率.化学组分的生成系由化学反应所致。当然，依反应实际上是产生还是消毁组分l 而定，R l 可能是正的，也可能是负的．对于不参与化学反应的组分，R l 为0．
如果用菲克(Fick)扩散定律来表示扩散流量密度J l ，我们可以写出：
l l l m J grad Γ-= (2.3) 其中，l Γ是扩散系数．将方程(2.3)代入方程(2.2)，可以导出：
l l l l R m m t
m +Γ=+??)grad div()div()
(u ρρ (2.4)
2．1—3 能量方程
最通用形式表示的能量方程式含有相当数量的各种不同影响因素．因为我们主要关心的是方程的形式而不是其细节，
所以考虑某些限定的情况就足够了。
对于可忽略粘性耗散作用的稳态的低速流，能量方程可写成
k S T k h +=)grad div()div(u ρ (2.5)
其中h 是比焓 k 是导热系数 T 是温度
h S 是容积发热率
根据傅里叶(Fourier)热传导定律，项)(kgradT div 代表流体中的热传导作用． 对理想气体以及固体与液体，我们可以应用下列关系：
h T c grad grad = (2.6) 其中c 是定压比热．将该式代入能量方程，就得到：
k S h c
k
h +=)grad div()div(u ρ （2.7）
如果c 是常数，则h ~ T 关系可以简化为：cT h = （2.8） 结果导出：
c
S T c k T k
+=)g r a d d i v ()d i v (u ρ （2.9）
在这种情况下，不管是焓还是温度都可以选作为因变量。 取方程中的速度u 为0，就得到了稳态的势传导方程：
0)grad div(=+k S T k （2.10）
2．1—4 动量方程
对于牛顿流体，控制给定方向上的动量守恒的微分方程可以写成类似线性的形式．但是由于必须同时考虑切应力和正应力，加之与流体流动有关的斯托克斯粘性定律要比菲克(Fick)定律或傅里叶定律复杂，因此动量方程就要复杂得多．用u 表示x 方向的速度，我们可以把相应的动量方程写成下列的形式：
x x V B u u t
u ++=+??)grad div()div()
(μρρu （2.11） 其中μ是粘度，p 是压力，Bx 是沿着x 方向的单位容积内的体积力，Vx 代表除去以)grad div(u μ所
表示的粘性力项之外的其它所有粘性力项。
2．1—5 紊流的时间平均方程
在实际应用中，常常会遇到紊流的问题．就工程中所遇到的实际问题，人们通常关心的是这种流动状态的时间平均特性．
因此人们通过平均运算的方法把不稳态的层流方程转化为对紊流流动的时间平均方程．
在进行这种平均运算的时候，人们假设：紊流中存在有相对平均值的快速而随机的脉动．由平均运算所产生的附加项是所谓的雷诺（Reynolds ）应力，紊流热流密度，紊流扩散流量密度等等．紊流模型的任务就是确定用流动的平均性质来表示这些附加量的方法．
许多紊流模型采用紊流粘度或紊流扩散系数的概念来表示紊流应力及流量密度．结果，紊流的时间平均方程就具有了与层流流动方程完全相同的形式。但是，诸如粘度、扩散系数以及导热系数这样一些层流交换系数则需要用相应的有效(即层流加紊流)交换系数取代。
从计算机计算的观点看，按照这种格式描述的紊流即相当于具有一个相当复杂的粘度表达式的层流。
x x V B u u t u ++=+??)grad div()div()
(μρρu 0)
()()(=??+??+??+
??z
w y v x u t
ρρρρ
x
p z u z y u y x u x z u w y u v x u u t u t t t ??-??+??+??+??+??+??=??+??+??+??)()()()()()()(μμμμμμρρρρy
p z v
z y v y x v x z v w y v v x v u t v t t t ??-??+??+??+??+??+??=??+??+??+??)()()()()()()(μμμμμμρρρρz
p
z w z y w y x w x z w w y w v x w u t w t t t ??-
??+??+??+??+??+??=??+??+??+??)()()()()()()(μμμμμμρρρρ 2.1—6紊流动能方程
现今普遍流行的紊流“双方程模型”把紊流脉动动能k 的方程做为其中的方程之一，该方程可以写作：
ρερρ-+=+??G k Γk t
k k )grad div()div()
(u （2.12）
其中Γ是k 的扩散系数，G 是紊流能量的生成率，以及ε是动能的耗散率．量ρε-G 是方程中的净源项．对变量
ε的微分方程与此相类似．
2.1—7通用微分方程 ?
通过上面简单地列举的一些有关微分方程，可以看出：在这里，我们所感兴超的所有的因变量似乎都服从一个通用的守恒原理．如果用Φ表示固变量，通用的微分方程就是：
S t
+Γ=+??)grad div()div()
(φφρρφu （2.13） 其中Γ是扩散系数，S 是源项．对于特定意义的Φ，具有特定的量Γ和S （实际上，我们本应采用符号
φ
φS 及Γ的，然而这样做会导致在随后的文章中应用过多的下标）．
上述通用微分方程中的四项分别是不稳态项、对流项、扩散项以及源项．因变量可以代表各种不同的物理量，如：化学组分的质量分量、焓或温度、速度分量、紊流动能或紊流的长度尺度。与此相应，对于这些变量中的每一个都必须相对应的扩散系数Γ以及源项S 赋以适当的意义．
把)gra div(φd Γ
作为扩散项：
并不是所有的扩散流量密度都是受制于有关变量的梯度的。但是把)gra div(φd Γ作为扩散项的做法并没有把通用的Φ方程只限于那些以梯度驱动的扩散过程．
凡是不能归入名义的扩散项的因子或项总是可以表示成源项的一部分。如果需要的话，我们甚至于可以把扩散系数取为0．由于大多数的因变量确实需要有一个突出梯度驱动扩散这一特性的扩散项，所以我们把具有这一特性的项以显式的形式写在通用的Φ方程中了．
出现在方程(2．13)中的密度可以通过状态方程与温度和质量分量这样一些变量相关联．这些变量与速度分布都服从通用微分方程．此外，流场还应当满足一个附加的约束条件，即质量守恒或连续性方程．这个方程是：
0)div(=+??u ρρ
t
（2.14） 上面，我们已经用向量的形式写出了有关的通用微分方程(2．13)和连续性方程(2．14)．这些方程的另一种有用的表达方式是直角坐标的张量形式：
S x x u x t j j j j +??Γ??=??+??)()()(φ
φρρφ （2.15） 0)(=??+??j j
u x t ρρ （2.16） 这里下标j 可以取信1、2、3分别代表三个空间坐标．如果在一项内下标j 重复出现，就意味着要取三项之和，例如：
)()()()(33
2211u x u x u x u x j j ρρρρ??+??+??=?? （2.17） )()()()(3
32211x Γx x Γx x Γx x Γx j j ????+????+????=????φφφφ （2.18） 用直角坐标张量形式表达的一个直接好处是，只要简单地把下标j 抹掉，就可以由这种形式得到该方程的一维形式．
把任何特定的微分方程改写成通用形式(2．13)的过程，就是把有关因变量的不稳态项、对流项及扩散项转换成共同的标准形式．于是，把扩散项内梯度Φ的系数取为对Γ的表达式；而把方程右端的其余各项之和定义为源项S ． 尽管到现在为止我们一直都是把所有的变量当成有因次的量来考虑，但是，以无因次的变量来进行研究则往往比较方便。可以认为任何一个用无因次变量表示的微分方程都具有通用的形式2．13，
S Γt
+=+??)grad div()div()
(φφρρφu (2.13) (重点) 这时Φ代表无因次的因变量，而Γ和S 分别代表扩散系数和源项的无因次形式．在许多情况下，Γ的无因次值可以简化为1，而S 可以取值0或1． 热、质传递，流体流动，紊流以及有关的一些现象的所有各有关微分方程都可以看成是通用的Φ方程的一个特殊情况． 我们需要关心的仅仅是方程(2．13)的数值解．在编制计算机程序时，只需要写出一个求解方程(2.13)的通用程序就足够了．们可以对不同意义的Φ，重复使用这个程序；当然需要对相应的Γ和S 分别赋以各自合适的表达式，同时也需要给出相应合适的初始条件与边界条件．
通用的Φ方程的概念使我们能够列出一个通用数值方法的公式，并编制通用的计算机程序．
2．2坐标的性质
到现在为止，我们都在注意因变量．现在我们将回过来考虑自变量（x,y,z,t ），并从计算的观点讨论一下它们的性质．
2．2—1 自变量
一般说来，因变量Φ是三个空间坐标与时间的函数．于是：
),,,(t z y x φφ= (2.19)
其中x ，y ，z 以及t 都是自变量．
并不是所有的问题都需要考虑所有四个自变量．所涉及的自变量数愈少，需要计算Φ值的位置〔或网格结点)也愈少． 一维问题：当有关的物理量只与一个空间坐标有关； 二维问题：所考虑的问题与二个空间坐标有关； 三维问题：所考虑的问题与三个空间坐标有关； 稳态问题：当所考虑的问题与时间无关时； 不稳态问题：当所考虑的问题与时间有关时。 把空间与时间两个因素放在一起考虑，我们将把一种状态说成是不稳态的一维问题，稳态的三维流动等等．
等温面--等温迁移法?
象方程(2．19)所表示的那种独立坐标的选择方式并不是唯一的形式．代替把稳态温度分布写成T(x ，y ，z)，我们可以用另外一种写法：
),,(T y x z z = (2.20)
这里z 变为因变量，它代表在位置(x ，y)相应于温度T 的等温面高度．
一种采用该表达形式的方法已经由迪克斯(Dix)和西切克(Cizek)(1970)以及克兰克(Crank)及其合作者[克兰克和费尔(Phahle)(1973)；克兰克和格普塔(Gupta)(1975)；克兰克和克劳利（Crowley)（1978）]研究出来了．这种方法叫做等温迁移法．但是这种方法只限于对坐标是单调函数的温度场；对于更为一般的温度场，一定的T 、x 和y 可能对应有几个高度z ；这样，从计算的目的看，z 就不适合作为因变量了．
2.2—2 坐标的合适选择
恰当明智地选择坐标系统有时可以减少所得要的自变量数目．由于网格结点的数目一般都与自变量的数目有关，因而以较少的自变量进行研究可以大大节省计算时间．
我们用几个特殊的例子来说明坐标的选择究竟是怎样影响自变量的数目的．
1．飞机周围的流体流动：在一个静止的坐标系上看以恒定速度飞行的飞机周围的流体流动是不稳态的；但是相对于固定在飞机上的移动坐标系而言，流动则是稳态的．
2．在一圆管内的轴对称流动：直角坐标系内是三维的，但是在r 、θ、z 的圆柱极坐标系内则是二维的，这是因为：
),(z r φφ= (2.21) 与θ无关。
3．坐标变换可能用来进一步减少自变量的数量，例如：
a ．在一块平板上的二维层流边界层给出了速度仅与η有关的相似性状态，其中： x
cy
=
η (2.22) 式中c 是一个有因次的常数．结果二维的问题便简化为一维的问题了．
b ．在半无限固体内的不稳态导热问题具有x 和t 这样两个自变量．但是，对于某些简单的边界条件，可以把温度描写成只与ξ有关，其中
t
Cx
=
ξ (2.23) 式中C 代表一个适当的有因次常数 4．改变因变量可能导致自变量数目的减少．例如：
a ．在一充分发展的通道流中，温度T 与流动方向的坐标x 以及横向坐标y 有关．但是在具有均匀壁温w T 的热发展区，
我们有： )(y θθ= (2.24) 式中：w
b w
T T T T --≡
θ T b 是整体温度，它随x 而变化。
b ．平面自由射流是一种二维流，但是我们可以写成：)(~~ηu u
= (2.25) 式中：c
u u
u ≡~, δηy ≡ (2.26)
这里u c 代表中心线上的流速，y 是横向坐标，δ是射流的特征宽度．u c 和δ这两个变量都随流动方向的坐标x 而变化． 为了提高计算效率，应当选择合适的坐标系统进行数值计算．
图3-2 一维问题的网格节点群
W P E x
w e w x )(δe x )(δx ?2.2—3 单向与双向坐标 (重点)?
【定义】双向的坐标：如果在一个坐标上的一个给定位置处的条件，只受该位置两侧条件变化的影响的话，那么这个坐标就是一个双向的坐标． 单向的坐标：如果在一个坐标上的一个给定位置处的条件，只受该位置一侧条件变化的影响，这样的坐标就是一个单向的坐标．
【例】通常，空间坐标是双向坐标。在一根园棒内的一维稳态导热是—个双向坐标的例子．在棒内任意给定点的温度会受到两端温度变化的影响．
时间坐标则总是单向坐标。在一块固体的不稳态冷却期间，其某一给定瞬时的温度只能受在该瞬时以前所发生的那些条件变化的影响．这是一个常识问题：昨天的事情影响着今天发生的事情，但是明天的条件对于今天会发生什么事情则没有影响．
空间也可以作为单向坐标：在流体流动的作用下，甚至某种空间坐标也能变成非常接近于单向的坐标．如果在一个坐标方向上有很强的单向流动，那么各种重要的影响只能是从上游传播到下游．于是在某一点上的状态主要受其上游条件的影响，而受其下游条件的影响很小．空间坐标的单向特性是一种近似．确实，对流是一种单向的过程
，
而扩散(这个过程总是存在的)则具有双向的效应．但是，当流量很大时，对流作用远远超过扩散作用，因而空间坐标就近乎是单向的了。
抛物型、椭圆型与双曲型——单向与双向坐标
看起来，通常用来进行微分方程分类的抛物型与椭圆型这两个数学术语相当于我们这里所用的计算方法上的概念——单向与双向坐标．抛物型这一术语表示一种单向的状态，而椭圆型则表示双向的概念。
那么前面提到过的第三种类型，即双曲型，又是什么意思呢？碰巧双曲状态不能准确归入计算方法上的分类．双曲型问题具有一种单向的特性，但是这种单向的特性并不是沿着坐标的方向，而是沿着一些称之为特征线的特殊线的方向。有一些采用特征线的数值方法，但是这些方法只限于双曲问题。另一方面，本书所提出的数
值方法并没有利用双曲问题所具有的这一特殊性质．我们将把双曲问题当成是椭圆问题这一总类中的一个组成部分(即全部是双向坐标)．
计算方法上的含义 上述关于单向和双向坐标的讨论，其动机在于：如果可以用一个单向的坐标来规定一个给定的状态，那么就有可能大大节省计算机的存储量与时间．现在让我们来讨论一个不稳态的二维热传导问题．我们将在计算域内构成一个二维网格结点的阵列。在任一瞬间，将存在一个相应的温度场．在计算机内将对时间上顺序推移的每一瞬时计算这样一个温度场，但是，由于时间是一个单向坐标，因而某一时刻的温度场不受未来时刻温度场的影响．实际上，整个不稳态问题
可以简化为按要求重复进行的一个基本的步骤，即给定t 时刻的温度场，求得t ＋Δt 时刻的温度场．于是计算机内存只要供这两个温度场使用即可；对于所有的各个不同时间步可以一遍又一遍地使用同样的储存空间．
在这样的情况下，从一给定的初始温度场开始，我们就能够“前进”到在时间上顺序推移的各个瞬时．在任何时间步，需要同时处理的未知量只有一个二维的温度数组．这个温度数组与所有未来的温度值无关，而影响这一温度数组的以前瞬时的值则是已知的．于是我们需要求解的就是一个简单得多的方程组，从而可以节省计算机的时间．
离散化方法
主要任务：我们所感兴趣的现象是受一些微分方程支配的。我们已经把这些微分方程用一个变量为Φ的通用方程全部归纳了．现在我们的主要任务就是推导求解这个通用方程的方法．
S Γt
+=+??)grad div()div()
(φφρρφu (2.13)
假设：为了便于理解，在本章中我们将假设变量Φ仅仅是一个自变量x 的函数
3.1数值方法的本质
3.1—1 任务
数值解系与实验数据：
一个微分方程的数值解系由一组构成因变量Φ的分布的数所组成。在这个意义上讲，数值方法有点类似于在实验室中进行的实验。做实验时，仪器的一组读数为我们构成在所研究的域内被测量
的物理量的分布，数值分析与实验室实验两者所获得的数据都只能是一些有限数量的数值。 未知量：各个不同位置上的Φ值
让我们假设用一个关于x 的多项式来代表Φ的变化
m m x a x a x a x a a +++++= 332210φ （3.1）
并采用数值方法来求得有限数量的系数,,,,210m a a a a 把x 的值以及各个a 的值代入方程(3．1)，我们就可以计算出在任何位置x 的Φ值．
但是如果我们最终的兴趣是得到在各个不同位置上的Φ值，那么这种做法就有点不太方便了。各个a 值本身是没有什么特别的意义的．要知道所需要的Φ值，还必须进行前面所述的代入计算过程．
想法：建立一个把一系列给定点上的Φ值作为原始的未知量的方法.实际上求解微分方程的大多数方法都属于这一类，因此我们将把我们的注意力放在这样一类方法上．因而数值方法就是把计算域内有限数量位置(叫做网格结点)上的因变量值当作为基本的未知量来处理．
该方法的任务是
提供一组关于这些未知量的代数方程 并规定求解这组方程的算法． 3.1—2 离散化的概念
离散化方法：把注意力集中在网格结点处的值，用离散的值取代包含在微分方程精确解中的连续信息．于是我们已经离散了Φ的分布，把这一类数值方法叫做离散化方法． 离散化方程：关于所选取的网格结点上未知Φ值的代数方程(我们现在就把它们叫做离散化方程)系由支配Φ的微分方程推导而得。
网格结点之间Φ的分布：
推导过程中，我们必须对网格结点之间Φ如何变化作某种假设．尽管可以选择在整个计算域内满足一个简单表达式的关系作为Φ的这种
“分布”。
采用分段分布是更为实际的方法，即：一定的段仅仅用一个小区域的内部及其边界上的网格结点上的Φ值来描述该区间内Φ的变化；于是一般将计算域分成一定数量的子域或单元．每一个子域可以有一个独立的分布假设．
离散化的概念：这样，连续的计算域已经被离散开了．这种对空间和因变量所作的系统的离散化使得我们有可能用比较容易求解的简单的代数方程取代前面提到过的控制微分方程。
一个离散化方程是连接一组网格结点处Φ值的代数关系式;这样的一个方程系由支配Φ的微分方程推导而得;表示与该微分方程相同的物理信息. 3.1—3 离散化方程的结构
Φ的分布 一定的离散化方程只与少数的几个网格结点有关，这种情况是我们选取分段分布这一特性的结果．因此，在一个网格结点处的Φ值只影响与其紧相邻的一些点上的Φ分布．可以预料到，当网格结点的数目变得很大时，离散化方程的解将趋近于相应微分方程的精确解。这是出自于这样的考虑：当网格结点紧挨在一起时，在相邻点之间的Φ变化就变得很小，因此有关分布假设的实际细节就不那么重要了．
离散化方程不是唯一的：对于一个已知的微分方程，可能的离散化方程决不是唯一的；尽管在网格结点数非常大的极限条件下，预计所有这些可能类型的离散化方程将会给出相同的解．离散化方程的不同形式起因于分布假设以及推导方法的不同．
有限差分和有限元法：可以认为这两种方法是离散化方法的两种可供选择的型式．有限差分法与有限元法之间的区别来自选择分布和推导离散化方程的方法不同．在本书内集中注意的方法具有有限差分法的外形，但是它采用了许多属于典型的有限元方法论所具有的思想．
图3-2 一维问题的网格节点群 W P E
x
w e w x )(δe x )(δx
?3.2推导离散化方程的方法
对于一个已知的微分方程，可以用许多方法推导出所要求的离散化方法。这里我们将简单地介绍几种普通的方法，然后再指出我们所喜爱的一个方法． 3.2—1泰勒级数公式
通过截断泰勒级数来近似表示微分方程的导数构成。对于图3．1中位于结点1与结点3之间中点的结点2 (2312x x x x x
-=-=?)，在2周围展开的泰勒级数给出：
在第三项之后截断级数，将二个方程相加及相减：
把这二个表达式代入微分方程就推得有限差分方程．
这个方法含有这样的假设：Φ的变化多少有点象是x 的一个多项式，从而高阶导数是不那么重要的。但是当存在(比方说)指数形式的变化时，这种假设就可能导致人们不希望要的那种公式．
泰勒级数公式的推导是比较直截了当的，但是缺乏弹性并且其中各项的物理意义难以理解。
3.2—2变分公式 3.2—3加权余数法 3.2—4控制容积公式 (略 见课本P31-34)
3.3 一个说明性的例子
让我们来讨论受下列方程控制的一维稳态热传导问题：
S Γt
+=+??)grad div()div()(φφρρφu 0)(=+S dx dT
k dx d （3.10）
其中k 是导热系数，T 是温度，以及S 是单位容积的发热率．
准备 为了推导离散化方程，使用图3.2中所示的网格结点群 ● 集中注意网格结点P
● 以网格结点E 及W 作为它的两个邻点（E 表示东侧，即正
的x 方向，而W 表示西侧，或是负的x 方向）．
● 虚线表示控制容积面；就目前来说，它们的准确位置是无关
紧要的．字母e 与w 代表这些面．
对于所考虑的一维问题，我们将假设在y 与z 方向为单位
厚度．于是，图中所示的控制容积是Δx×1×1。如果我们在整个控制容积内积分方程，我们就得到：
0)(=+S dx
dT k dx d (3.10) 0)()(=+-?dx S dx
dT
k dx dT k e w w e (3.11)
分布曲线的假设 需要一个分布假设或是一个内插分式．图3.3中表示了两种简单的分布假设．
阶梯式分布：假设在一个网格结点处的T 值代表它周围整个控制容积内的T 值（阶梯式分布），对于这样的分布，
斜率
dx dT /在控制容积面(即在w 或e)上是不确定的
分段线性的分布(图3.3b)．这时，在网格结点之间采用线性的内插函数．
离散化方程 如果我们用分段线性分布来计算方程(3.11)中的
dx dT /，
0)()(=+-?dx S dx
dT
k dx dT k e w w e (3.11)
所得的方程将为：
0)()
()()(=?+---x S x T T k x T T k w
W P w e P E e δδ (3.12)
其中S 为S 在整个控制容积内的平均值．离散化方程可缩写成下列形式：
b T a T a T a W W E E P P ++= (3.13)
其中： e
e
E x k a )(δ= w w W x k a )(δ= W E P a a a += x S b ?=
说明
1. 离散化方程的标准形式: 方程(3．13)
方程的左边: 在中心网格结点上的温度Tp ；右侧:由相邻结点上的温度和常数b 构成，对二维与三维的情况相邻
结点的数目增加． 一般说来，比较方便的是把方程(3.13)看成具有如下的形式：
b T a T a nb nb P P +=∑ （3.15）
式中下标nb 表示一个相邻结点，∑表示对所有的相邻结点求和．
2. 分布假设: 在我们推导方程(3．13)时，我们已经采用了使我们能够估计dx dT /值的最简单的分布假设，当然选用许多其它形式的内插函数本来也是可以的．
3. 此外，我们没有必要对所有的量都采用同样的分布函数，明白这一点是重要的．例如，既没有必要用网格结点之间线性变化的S 来计算S ，也没有必要由kp 和ke 之间线性变化的k 来计算ke 。
4. 即便对于一个确定的变量，也没有必要对方程中所有各项都采用同样的分布函数假设．例如，倘若在方程(3．10)中含有一个单独包含T 的附加项，或许对该项使用阶梯分布函数也是可以允许的，而不必坚持如在计算dx dT /时所采
用的分段线性分布．
指导原则 上述有关选择分布函数的自由度，最终导致不同变型的离散化方程形式．事实是，当网格结点的数目增加时，可以预料所有这些不同形式的方程都会给出相同的解。 要求：即便是采用很粗的网格，解也总应该满足：（1）物理上的真实性；（2）总的平衡。但是，我们将在此提出一个附加的要求，这个要求的加入将使我们能够大大减少我们可以接受的公式(方程)的数目．
物理上的真实性是容易理解的．图3．4所示的变化说明了这一概念．一个真实的变化应当具有与准确变化相同的定
性倾向．在无内热源的热传导问题中，内部没有一处的温度可以超出边界温度所确定的温度范围之外．当一块热的固体为绕流的流体所冷却时；固体的温度不可能降低到比该流体的温度还低．我们将总是用这样的试验来检验我们的离散化方程．
总平衡的要求意味着对整个计算域应当满足积分守恒．我们将坚持要求热流密度、质量流量以及动量通量必须准确地同相应的源和汇建立平衡，这种平衡不应当只是限于网格结点的数目变得很大时的情况，而是对于任何数目的网格结点都应该得到满足。（1）控制容积公式会使这一总的平衡成为可能，但是如同我们很快就会知道的那样，（2）在计算控制容积界面的热流密度、质量流量以及动量通量时还需要小心处置．
物理上的真实性以及总的平衡这两个约束条件将用来指导我们选择分布假设以及所采用的有关措施．在这些约束条件的基础上，我们将建立起几条基本的法则，应用这些法则我们就可以对现有的一些公式进行鉴别，并进而发展出一些新的公式来．通常需要由数学上的考虑才能作出的一些决断，现在可以直接由物理学上的原理来进行指导了。
源项的处理 一般来说，源项是因变量T 本身的函数，因而在构成离散化方程的过程中需要知道这种函数关系。由于离散化方程需要用线性代数的技术来求解，因而我们在形式上只能考虑一种线性的函数关系。在下一章中将讨论对一个已知S~T 关系进行“线性化”处理的方法。这里把平均值S 表示成下列形式是足够的。
P P c T S S S += (3.16)
式中，Sc 代表S 的常数部分，而Sp 是Tp 的系数（显然Sp 不代表在结点P 所计算的S ）．
在方程(3．16)中Tp 的出现表明，在表示平均值S 时，我们已经假设：Tp 值代表整个控制容积内的值，换句话说，已经采用了图3.3中中所示的阶梯式分布(应当注意：在对dx dT /项采用分段线性分布的同时，我们可以自由地对源项采用阶梯式分布)。
应用线性化的源项表达式，离散化方程的样子看起来仍然象方程(3.13)，但是系数的定义[方程(3.24)]要有所改变．新的方程组是：
b T a T a T a W W E E P P ++= (3.17)
其中：
e
e E x k a )(δ=
w
w
W x k a )(δ=
x S a a a P W E P ?-+= x S b c ?=
上述导言性的讨论为我们提供了足够的预备知识，这样，我们就可以来推导我们的离散化方程所应当服从的一些基本法则公式，以确保所得到的解满足物理上真实以及总的平衡这两个要求．这些看来简单的法则具有深远的含义，它们将指导我们推演贯穿本书的那些方法．
3.4 四项基本法则
法则1：在控制容积面上的连续性 当一个面作为两个相邻控制容积的公共面时，在这两个控制容积的离散化方程内必须用相同的表达式来表示通过该面的热流密度、质量流量以及动量通量．
讨论：显然，通过一个特定的面离开一个控制容积的热流通过同一个面进入第二个控制容积的热流密度相同，不然的话，总的平衡就不会满足．尽管这一要求易于理解，但是稍不小心，就可能在一些微妙的地方违反。对于图3.2所示的控制容积，我们或许会用
通过TW 、TP 及TE 的二次曲线来计算界面上的热流密度dx
dT
k
。对下一个控制容积采用同一类的公式意味着：公共界面上的梯度dx dT /系由不同的、与正要考察的控制容积有关的分布曲线算得．所得到的dx dT / (因此热流密度)的不连续性绘于图3．5中．
另一种做法也可能导致热流密度的不连续性．这时假定在给定的控制容积的各个表面上，热流密度完全为控制容积中心结点的导热系数kp 所控制．这样在考虑P 点周围的控制容积时，在界面e 处的热流密度(见图3．2)将表示成
e E P p x T T k ) /()(δ-。而在把E 作为控制容积的中心结点时，热流密度应为e E P e x T T k ) /()(δ-．为了避
免出现这种不连续性，记住这一点是有好处的，即：必须把通过界面上的热流密度看成是属于界面本身，而不是属于一
定的控制容积的．
法则2：正系数 大多数的实际问题应该是这样的，即某个网格结点处的因变量值只是通过对流以及扩散的过程才受到相邻网格结点上的值的影响．在其它条件不变的情况下，在一个网格结点处该因变量值的增加应当导致相邻网格结点上该值的增加（而不是减少）
在方程（3.13)中，如果T E 的增加必然导致T P 增加的话，这就必然要求系数a E 与a P 具有相同的符号。
挽句话说，对于通用性方程(3．15)，中心结点的系数a P 与各相邻结点的系数a nb 全都必须具有相同的符号．当然，让我们可以把它们全部取为正值或者全部取为负值。让我们决定这样来写我们的离散化方程，使所有的系数均为正值，这样，法则2可以表述如下：
所有的系数(a P 以及各相邻给点系数a nb )必须总是正的．
b T a T a T a W W E E P P ++= (3.13) b T a T a nb nb P P +=∑ （3.15）
说明 由方程(3．14)所给出的系数的定义表明，我们所用的离散化方程[方程(3．13)]遵守正系数法则．但是，正如我们将在后面看到的那样，有许多公式经常违反这一法则．结果往往是得到一个物理上不真实的解．存在负相邻结点系数就会导致这样的情况，在这种情况下，一个边界温度的增加会引起相邻网格结点上的温度的降低．我们将只接受那些确保在所有情况下系数均为正的公式．
e
e E x k a )(δ=
w
w
W x k a )(δ=
W E P a a a += x S b ?= (3.14)
法则3：源项的负斜率线性化 要是我们研究一下方程(3.18)中系数的定义，就可以发现，即便所有相邻结点的系数均为正，由于S P 项的关系，中心结点的系数a P 仍可能变为负值．当然，只要S P 不为正值，这一危险就完全可以避免．
e
e E x k a )(δ=
w
w
W x k a )(δ=
x S a a a P W E P ?-+= x S b c ?= (3.18)
于是我们把法则3写成：当源项线性化为
P P c T S S S +=，系数S P 必须总是小于或是等于0．
注意 这一法则看起来好像有点任意，其实并不然．大多数物理过程确实在源项与因变量之间具有负的斜率关系，实际上，要是S P 为正的话，物理状态就可能会变得不稳定了．一个正的S P 意味着，当T P 增加时，源项也随着增加；如果这时没有有效的散热机构；这可能会反过来导致T P 的增加，如此反复进行下去，造成温度飞升的不稳定现象．从计算方法上讲，保持负的S P ，使之不致产生不稳定性以及物理上的不真实解，这是至关重要的．在下一章中，我们还要进一步讨论源项线性化的问题．这里要充分注意到：为了使计算成功，负S P 的原则是必不可少的．
法则4：相邻结点系数之和 控制微分方程往往只包含有因变量的导数项．
0)(=+S dx
dT k dx d （3.10） 于是，如果T 代表因变量，则函数T 与T+c (其中c 是一个任意常数)两者均满足微分方程．微分方程所具有的这一特性也必定要反映在与之相对应的离散化方程中．因此，当T P 以及所有的T nb 都增加同一常数值时，方程(3．15)应当仍然适合．
b T a T a nb nb P P +=∑ （3.15）
b C T a C T a nb nb P P ++=+∑)()( ∑=C a C a nb P
由这个要求可以得出结论：a P 必须等于所有相邻结点的系数之和．因此，我们可以把法则4表述如下： 为了使微分方程在因变量增加一个常数之后也仍然能得到满足，我们要求：∑=nb P a a (3.19)
讨论 很容易理解，方程(3.13，当源项与T 无关时)确实满足这一法则．
b T a T a T a W W E E P P ++= (3.13)
本法则意味着：中心结点值T P 是各相邻结点值T nb 的一个加权平均值．与方程(3.13)不同，方程(3.17，当源项与T 有关时)中的系数不服从这一法则．
b T a T a T a W W E E P P ++= (3.17)
当源项与T 有关时:不能认为这种情况是对法则4的违背，而应该是本法则对这种情况不适用．当源项与T 有关时，T 与T ＋c 两者不能同时满足微分方程．就是在这样的情况下，也不应当把这个法则忘记掉，而应该通过设想方程(3.17)的一个特殊情况来应用这个法则．例如，如果我们取方程(3.17)中的S P 为0，那么本法则就变得可以应用了，并且确实是被遵守的。
温度场的确定性：当T 以及T+c 两者均满足微分方程时，所要求的温度场并不会变成是多值或不确定的．T 的值可以由适当的边界条件唯一地确定．遵照法则4就可以确保：例如，倘若边界温度增加一个常数值，那么所有的温度就会准确地增加同一常数值．
理解法则4作用的另一方法是：若是源项不存在，而且所有的相邻结点的温度T nb 都相等，则中心结点的温度
T P 必定等于它们．在这样的条件下，只有一个蹩脚的离散化方程才会算出T P ≠T nb 的结果来．
第四章 热传导
4．1本章的对象 (略 见课本P46) 4．2 一维稳态热传导 4．2—1 基本方程
作为解释四项基本法则的工具，在3．3节中介绍说明性例子的过程中，我们已经推导了稳态一维导热问题的离散化方程．我们来复习一下那一部分的主要内容．对这个问题的控制微分方程是：
0)(=+S dx
dT k dx d 由这个方程可以推导出离散化方程
b T a T a T a W W E E P P ++= (4.2)
其中：
e
e E x k a )(δ=
w
w
W x k a )(δ=
x S a a a P W E P ?-+= x S b c ?= (4.3) 网格点P 、E 及W 示于图3．2中，图中还标明了各个不同的距离．控制容积面位于网格点P 及其相应的相邻点之间．可以认为这些面的准确位置是任意的．它们的位置可以有许多不同的安排方案。其中的一些将在4．6—1节中讨论．就目前来说，我们将简单地认为e 和w 相对于网格点P 、E 和W 的位置是已知的．
量Sc 及S P 由源项的线性化规定，其形式是：P P c T S S S
+= (4.4)
至于分布假设，梯度dx dT /已经由T 对x 的分段线性的变化算得，而对于线性化的源项，假设T P 值代表整个控制容积内的值．当然，应当记住：只要不违背四项基本法则，选择其它形式的分布曲线是可能的，也是允许的．这里，选择分布曲线的基本原则是：在四项基本法则的约束范围内宁可采纳简单一些的分布曲线，而只是在那些必要的场合，才采用高级而复杂形式的分布曲线．
一维热传导问题的许多重要方面仍然有待于讨论，我们该转到这些论题上来了． 4．2—2 网格间距 非均匀网格
对于在图3．2中所示的网格点，距离w e x x ) () (δδ与没有必要相等．事实上，人们常常希望采用不均匀的网格间距．这样做可以使我们能够有效地扩大我们计算的功能．一般说来，只有在网格相当细的时候，我们才可能得到精确的解．但是在因变量随x 变化相当慢的区域内没有必要采用细的网格．相反，在T~x 变化比较陡的区域内则需要细的网格． 不均匀网格的准确度要比均匀网格的差？？
现在在一些人们中间似乎流行着一种错误的概念，这些人以为不均匀网格的准确度要比均匀网格的差些．这样的论断是没有合理依据的．网格间距应当直接与因变量在计算域内的变化方法联系起来．此外，还没有一个通用的法则可以说明相邻网格间距之间最大(或最小)比值究竟应该保持多大．
如何设计一个合适的非均匀网格？？
因为在问题解决之前，T~x 的分布是不知道的，我们能够怎样来设计一个合适的非均匀网格呢？ （1）首先，人们通常对于所要得到的解有着某些定性的预计，他们可以从这些预计得到某些指导． （2）其次，可以用初步的粗网格的解来求得T~x 变化的形式，而后可以构成一个合适的非均匀网格．
这就是为什么我们坚持要求我们的方法，即便是在粗网格的条件下也应该给出具有物理意义的解的原因之一。如果这种方法只是对足够细的网格才能给出合理的解，那么一个试探性的由粗网格求得的解将是无用的． 达到给定精确度所需要的网格点数及分布？？
达到给定精确度所需要的网格点数以及这些网格点在计算域内应取的分布方式与所要求解的问题的特性有关．采用仅仅几个网格点进行试探性计算，为弄清有关解的情况提供了一个方便的途径．最后，这恰恰是人们在实验室中相当普遍地采用的一种做法．在那里，人们先进行预备性的实验或是探索性的试验，然后再应用由这些实验(或试验)所得到的数据资料来确定在最终的实验中所应安装的探头的位置与数目。 4．2—3 界面导热系数 界面导热系数：在方程(4．3)中，已经把导热系数k e 用来代表通过控制容积面e 的k 值，类似地，k w 代表界面w 的k 值．当导热系数k 是x 的函数时，我们往往只知道在网格点W 、P 、E 上的k 值．于是我们需要有一个如何由这些网格点值来计算界面导热系数(如k e )的规定。当然，下列讨论并不是针对均匀导热系数的情况而言的．
不均匀导热系数：不均匀导热系数的状况可以由材料的不均匀性(如组合的材料扳)引起．即便对于均匀的材料，由于温度分布的不均匀，加上导热系数随温度而变化的函数关系，也会导致导热系数发生变化．在处理对Φ的通用微分方程时，扩散系数Γ将起着与导热系数相同的作用．在实际问题中，经常会遇到Γ发生明显变化的情况．例如，在紊流的情况下，Γ可以代表紊流粘度，也可以表示紊流导热系数．因此人们迫切希望能有—个适合于不均匀的k 或Γ的公式。
求得界面导热系数k e 的最直截了当的方法是假设k 在P 点和E 点之间呈线性变化。于是，
E e P e e k f k f k )1(-+= (4.5)
其中插入因子
e f ，是用图4.l 所示的距离所定义的一个比值：
e
e e x x
f )()(δδ+
=
(4.6)
导热系数的突然变化的情况：我们马上就要说明：在某些情况下这种简单化的做法会导致相当不准确的结果；此外，这种做法不可能精确地处理在组合材料中可能遇到的导热系数的突然变化．在阐明这种替代的办法时，我们认为我们主要关心的不是导热系数在界面e 上的局部值．我们的主要目标是要得到一个通过下式描述的界面热流密度q e 的良好的表达式：e
E P e e x T T k q )()
(δ-=
(4.7)
事实上，在推导离散化方程(4-2)时已经采用了这一表达式，我们想要采用的k e 的表达式应当是一个可以由它表达“正确的”q e 的式子．
让我们来讨论这样一种情况，围绕着网格点P 的控制容积由具有均匀导热系数k p 的材料所填满，而围绕着E 点的控制容积则由导热系数为k E 的材料所填满．对于在P 点与E 点之间的组合板，根据稳态无内热源一维导热的分析，可以得
出下列的结果：E
e P e E
P e k x k x T T q /)(/)(+-+-=
δδ (4.8)
把方程(4．6)一(4．8)合并在一起，就得到：1
)1(-+-=E
e P e e
k f k f k (4.9)
当界面l 位于P 和E 之间的中点时，我们有f e ＝0.5；于是：
)(5.01
11
---+=E P k k k e (4.10) 或
E
P E
P e k k k k k +=
2 (4.10b)
方程(4．10)说明e k 是P k 和E k 的调和平均值，而不是如方程(4．5)在e f ＝0.5时给出的算术平均值．
将方程(4．9)应用于系数的定义[式(4．3)]中，可以得到E a 的下列表达式：
1
)()(-+-??
????+=E e P e E k x k x a δδ (4.11)
对W a 可以写出一个类似的表达式．很清楚，E a 代表点P 和E 之间的材料的热导． 这一公式的效能可以由下列两种极限情况很快看出： 1．
令E k →0，则由方程1
)1(-+-=E
e P e e k
f k f k (4.9) 可得：0→e k (4.12)
这意味着在一个绝热层的表面上热流密度为零，正象它理应存在的情况那样．相反，算术平均公式在这种情况下均会给出非零的热流密度．
2．
令P k >>
E k ，那么：1)1(-+-=E e P e e k f k f k (4.9) e
E
e f k
k → (4.13)
这个结果具有两方面的含义：其中的一个含义是显而易见的，而另外的一个含义则比较难于理解．方程(4．13)表明界面的导热系数k ，完全与k p 无关．这个结果是可以预料到的，因为围绕着P 点的材料导热系数高，其热阻与围绕着E 点的材料相比可以忽略不计(算术平均公式却保留k p 对k e 的影响)。另一个含义是：k e 不等于k E ，而是它的1/f e ．稍微转一下弯子就可以明白这一含义是合理的．我们的目的是通过方程(4．7)得到一个正确的q e 值．应用方程(4．13)得出：
+
-=e E P E e x T T k q )()(δ (4.14)
当P k >> E k 时，温度P T 将一直扩展到界面e ，而温降P T —E T 将实际上发生在距离+e x ) (δ内．于是正确的热流密度将由方程(4．14)给定．换句话说，在方程(4．13)中的因子f e 可以看成是对方程(4．7)中所用的名义距离e
x ) (
δ的补偿．
对这个极限情况的讨论表明，这个公式可以适用于导热系数突然变化的情况，而无需在发生突变的邻近区域内采用极细的网格．这样做不仅对组合板的导热计算比较方便，而且还有更加引人入胜的其它一些含义．这些更深的含义已经在帕坦卡(1978)的论文中作了描述．在本书中，我们也将在随后的几章中加以说明．
所推荐的界面导热系数公式(4．9)是建立在无内热源的稳态一维导热状态基础上，其中导热系数在相邻的两个控制容积之间发生阶跃变化．即便在内热源不为零或是导热系数连续变化的场合这种调和平均的表达式也要比算术平均公式好得多。在帕坦卡(1978)的论文中已经针对某些可以求得精确解析解的情况说明了这一点． 4.2—4 非线性
离散化方程(4．2)是一个线性的代数方程
b T a T a T a W W E E P P ++= (4.2)
我们将用解线性代数方程的方法来解这样的方程组． 非线性的情况
但是即便在热传导的问题中我们也常常会遇到非线性的情况。导热系数k 可能与T 有关，或者源项S 可能是T 的一个非线性函数．这样，离散化方程中的系数本身将与T 有关．我们将用迭代的方法来处理这种状况．这个过程包括如下几步：b T a T a T a W W E E P P ++=
1．一开始在所有各个网格点上，猜测或估计一组T 值． 2．由这些估计的T 值，计算出离散化方程中的系数的试探值． 3．解名义上的线性代数方程组，得到一组新的T 值． 4．以这些T 值作为较好的估计值，返回到第二步并重复这个过程，直到这种进一步的重复计算(称之为迭代)不再引起T 值任何有意义的变化时为止．
这种最终的不变的状态称之为迭代的收敛。这个收敛的解实际上是非线性方程的正确解，尽管这个解是由线性方程的求解方法得到的．
但是有可能一次次的迭代永远也不会收敛到一个解．T 的值可能稳定地漂移或是以一个不断增大的振幅振荡．这一与收敛相对立的过程称之为发散．一个好的数值方法应当使发散的可能性减为最小；如我们在后面将要看到的那样，只要遵守我们的四项基本法则，就会促进收敛；此外，我们还将讨论避免发散的其它一些对策．在这一点上，要充分注意我们的方法不止是限于线性的问题．至少在原则上，任何的非线性问题都可以用刚才概述的选代技术处理． 4．2—5 源项的线性化
当源项S 与T 有关时，我们采用方程(4．4)所给出的线性形式来表达这一关系．
P P c T S S S += (4.4)
这样做的原因是：(1)我们名义上的线性结构框架只允许采用一种形式上的线性关系，(2)线性关系的组合要比把S 处
理成常数为好．
当S 是T 的一个非线性函数时，我们必须把它线性化，即规定S c 和S p 的值，而这两个值本身也可能是T 的函数。因而在每一次迭代循环，S c 与S p 可能都要根据新的T 值重新算出。S 的线性化应当是S~T 关系的一个良好的表达式。此外，还必须满足非正的
p S 这一基本法则．
有许多方法可以用来把给定的S 的表达式分解成S c 与S p T p ；其中某些方法可以用下面的例子来说明．出现在这些例
子中的数值并不具有特别的意义；符号*
p T
用以表示T 的估计值或是前一次的迭代值．
[例1] 已知：S ＝5—4T ．某些可能的线性化是： 1．
c S ＝5，p S ＝-4．这是最明显的形式，并且是我们所推荐的．
2．
c S ＝5-4*
p T ，p S ＝0．这是懒汉的做没他把整个S 都写成人c S ，而令p S ＝0．但是，这种做法并不是不
切实际的．当S 的表达式很复杂时，这样做或许是我们唯一的一种选择．
3．
c S =5+7*p T ，p S ＝-11．这得出比实际的S~T 关系更陡的曲线．这样做的结果特使迭代的收敛速度减慢了．但
是倘若在所研究的问题中还存在着其它的非线性项时，这种减慢的做法实际上可能是受欢迎的．
[例2] 已知S=3+7T ．某些可能的线性化是：
1．c
S =3，p
S ＝7．一般说来，这是不能接受的，因为它使
p S 为正．如果所研究的问题不用进行迭代就可以直接
求解，那么这种线性化是可以给出正确的解的．但是，倘若由于某种原因(如在方程中其它有一些项是非线性的)而需要采用迭代，则存在正的
p S 就可能会导致解的发散．
2．c
S ＝3+7*p
T ，p
S =0．这是人们在没有自然现成的p
S 可得的情况下，应当采取的措施．
3．c
S ＝3+9*p
T ，p S = -2．这是一种人为产生的负的p
S 。一般说来，这样做的结果将导致收敛速度减慢．
[例3] 已知：S=4-5T 3．某些可能的线性化是： 1．c
S ＝4-5*p
T
3，p S ＝0．这是懒汉的做法，这种做法不能很好地利用已知S~T 关系这一有利条件．
2．c
S ＝4，p
S ＝-5*p
T
2．这看起来像是准确的线性化，但是已知的S~T 曲线要比这一关系所反映的曲线为陡。
3．推荐的方法：)(1554)(*
2*3***
*P P P P P P T T T T T T dT dS S S ---=-??
? ??+=
x
?B
q i
x )(δB
i
I
图4-4 边界附近的半控制容积
于是：
2*3*15104P
P P
c T
S T
S -=+= 这一线性化表示，在
*
p
T 点，所选
择的直线与S~T 曲线相切。
4．c S ＝4+20*p
T 3，
p S ＝-25*p T 2
．这一线性化要地已知的S~T 曲
线为陡，它会使收敛速度降低．
这四种可能的线性化与实际的S~T 曲线一起示于图4．2中．在这样的图上，正斜率的直线会违反基本法则3．在所有负斜率的直线中，与已知曲线相切的直线通常是最佳的选择．较陡的斜率是可以接受的，但是通常会导致收敛速度减慢．欠陡的直线我们是不希望采纳的，因为它不能体现已知的S 随T 的下降速度．
就本章所要达到的目的而言，关于源项线性化的讨论到此结束． 4．2—6 边界条件
对于一维问题，让我们来讨论如图4.3所示的网格点组． 在两个边界上各有一个网格点．其余的网格点称为内点．
围绕着每一个内点有一个控制容积．对每一个这样的控制容积可以写出一个象方程(4.2)那样的离散化方程．如果把方程(4.2)看做是关于
p T 的一个方程，那么我们就有了对所有内网格点上未知温度所必要的方程．
b T a T a T a W W E E P P ++= (4.2)
其中有两个方程包含着边界网格点上的温度．通过处理这些边界温度，就把已知的边界条件引入到我们的数值解法中了．
由于没有必要分别讨论两个边界点，因而我们讨论左边的边界点B ，如图4．3所示，该点与第一个内点I 相邻．在热传导问题中有三类典型的边界条件，它们是：
1．已知边界温度 2．已知边界热流密度 3．通过放热系数和周围流体的温度来规定边界的热流密度．
如果边界温度已知(即，如果B T 的值是已知的)，这种情况并无任何特别的困难，处理时也不需要外加任何方程．
当边界温度未知时，对
B T ，我们需要构成一个附加的方程．这样的方程可以通过在图4．3所示的边界附近的“半”
控制容积内对微分方程进行积分来建立(该控制容积只包括位于网格点B 一侧的半个容积，因而我们把它叫做半控制容积)．该控制容积的放大图在图4.4中给出．在该控制容积内对方程(4．1)进行积分，
0)(=+S dx
dT k dx d （4.1）
并考虑到热流密度q 代表
dx kdT /-我们就得到：
0)(=?++-x T S S q q B P c i B （4.15）
式中源项已经按通常的型式线性化．界面热流密度
i q 可以按照方程(4．7)的格式写出．结果是：
0)()()
(=?++--x T S S x T T k q B P c i
I B i B δ （4.16）
这一方程的进一步表达形式与如何给定边界上的热流密度
B q 有关．如果B q 值本身是已知的，则所要求的对B T 的方程变
成：b T a T a I I B B += (4.17)
其中：
i
i I x k a )(δ=
B c q x S b +?= x S a a P I B ?-= (4.18)
如果热流密度B q 系由放热系数h 以及环境流体温度f T 规定，那么：)(B f B T T h q -= (4.19)
于是对
B T 的方程为：b T a T a I I B B += (4.20)
式中：i
i I x k a )(δ=
f
c hT x S b +?=
h x S a a P I B +?-= (4.21c)
这样我们就能够构成对所有未知温度的足够数量的方程． 现在我们该来描述求解这些方程的方法了． 4．2—7 线性代数方程的解
一维离散化方程的解可以用标准的高斯(Gauss)消去法得到，由于方程的形式特别简单，消去过程的算法就变得十分方便．有时候，这种算法称之为托马斯(Thomas)算法或是TDMA (三对角矩阵算法)．TDMA 的名称来自于这样的事实，即：在写这些方程的系数矩阵时，所有的非零系数均排列在矩阵的三条对角线上．
b T a T a T a W W E E P P ++=
(4.2)
对方程组AX ＝B
为了表示该算法的方便起见，应用某些不同的符号是必要的．设在图4．3中的网格点记标号为1， 2， 3，…，N ，其中1与N 代表边界点．离散化方程可以写成为：
i i i i i i i d T c T b T a ++=-+11 (4.22)
式中下标i ＝l ，2，3，…，N ．于是温度i T
与相邻的温度1+i T 及1-i T 有关．考虑到边界点方程的特殊的形式，让
我们令：
1c ＝0 及 N b ＝0 （4.23）
12111d T b T a += )(211T f T =
图
4-3 内点与边界点的控制容积
B
i I
W
P E
半控制容积
典型的控制容积
x
2123222d T c T b T a ++= )(322T f T = 3234333d T c T b T a ++= )(433T f T =
……… )(11N n N T f T --= N N N N N d T c T a +=-1 所以温度
0T 及1+N T 不起任何作用．(当边界温度己知时，对这些边界点的方程就只剩下一个无意义的形式．
例如，如果1T
已知，我们就有：1a =1，1b ＝0，1c ＝0，以及1d ＝1T 的已知值)．
如果1T
未知，这些条件意味着
1T 是2T 的已知函数．
12111d T b T a += )(211T f T =
i ＝2的方程是
1T ，2T 与3T 之间的一个关系式． 1T ，2T 与3T 之间的关系式就简化为2T 与3T 之
间的一个关系．换句话说，
2T 可以用3T 的一个关系式来表示．
2123222d T c T b T a ++= )(322T f T =
这一代入过程可以一直继续到把N T 表达成1+N T 的一个关系式．但是由于实际上1+N T 并不存在，所以到这一
步我们实际上也就得到了
N T 值．
这使我们可以由此开始“回代”的过程。由
N T 得到1-N T ，由1-N T 得到2-N T ，…，由3T 得到2T ，以及由
2T 得到1T 。这就是TDMA 的要点．
12111d T b T a += )(211T f T = 2123222d T c T b T a ++= )(322T f T = 3234333d T c T b T a ++= )(433T f T =
………
)(11N n N T f T --=
N N N N N d T c T a +=-1
假设在前向代入的过程中，在我们刚刚得到：111---+=i i i i Q T P T (4.24)
以后，我们试图找到一个关系：i i i i Q T P T +=+1 (4.25)
将方程（4.24）代入方程（4.22）
i i i i i i i d T c T b T a ++=-+11 (4.22)
就得到：
i i i i i i i i i d Q T P c T b T a +++=--+)(111 (4.26)
这个方程可以改写成方程(4．25)的样子．换句话说，系数
i P 与i Q 于是可以表示成：
1
--=i i i i
i P c a b P (4.27a) 11---+=i i i i i i i
P c a Q c d Q (4.27b)
这些关系式是递推的关系，因为它们给出了
i P 及i Q 与1-i P 及1-i Q 之间的关系．为了开始这个递推过程，我们
注意到：对i=1的方程(4．22)几乎就是方程(4．25)的样子．
i i i i i i i d T c T b T a ++=-+11 (4.22)
1c ＝0 及 N b ＝0 （4.23）
于是i P 与i Q 的值由下式给出：111a b P = 及 1
1
1a d Q = (4.28)
[值得注意的是，在把1c
=0代入方程(4．27)之后，
1--=
i i i i i P c a b P (4.27a) 1
1
---+=i i i i i i i
P c a Q c d Q (4.27b)
就可以得到这两个式子．]
在i
P ，i
Q 序列的另一端，我们注意到N
b =0．这导致N
P =0，因此由方程(4．25，)
i i i i Q T P T +=+1 (4.25)
我们得到：N N
Q T = (4.29) 现在我们已经可以开始通过方程(4．25)来回代了．
算法概要
i i i i i i i d T c T b T a ++=-+11 (4.22)
1．由方程(4．28)计算1P ，1
Q 。 111a b P =
及 1
1
1
a d Q = (4.28)
2．应用递推关系式(4．27)得到对i ＝2，3，…N ，的
i P 及i Q 。
1--=i i i i i P c a b P
1
1---+=i i i i i i i P c a Q c d Q (4.27)
3．令
N T ＝N Q 。
4．对i ＝N-1，N-2，…，3，2，1应用方程（4.25）得到
1-N T ，2-N T ，…，3T ，2T ，1T 。
i i i i Q T P T +=+1 (4.25)
不管什么时候，只要代数方程可以表达成方程(4．22)的形式，那么三对角矩阵的算法就是一个非常有效而又方便的
求解方法．与一般的矩阵方法不同，TDMA 需要的计算机贮存量及计算机时间只是正比于N ，而不是N 2或N 3。
4．3 不稳态一维热传导
4．3—1 通用的离散化方程
图3-2 一维问题的网格节点群 W
P
E
x
w
e
w
x )(δe
x )(δx
?S Γt
+=+??)grad div()div()
(φφρρφu (2.13) 已经知道如何来处理：
一维的状态下Φ的通用微分方程中的扩散项与源项． 试探求解不稳态一维热传导微分方程 ?
)(x
T k x t T c ????=??ρ （4.30）
● 考虑不稳态项 ● 将源项暂时撇下不管；对于源项，我们已经没有什么新东西可以谈了． ● 假设ρc 是常数
因为时间是一个单向坐标，我们由一已知的初始温度分布开始，沿着时间坐标逐步向前求解。 计算的任务是：已知t 时刻T 在网格点上的值，求得在t+Δt 时刻在相应网格点上的值。
在网格上“老的”(已知的)T 值用0P T ，0
E T ，0W T
表示，
在t+Δt 时刻的“新的”(未知的)T 值则以1P T ，1
E T ，1W T
表示．
现在对图3．2所示的整个控制容积积分方程(4．30)以推导离散化方程，其中积分的时间间隔为由t 到t+Δt ．
)(x
T
k x t T c ????=??ρ （4.30）
于是：
dxdt x
T
k x dtdx t T c t t t e w e w t
t t
???
?
?+?+????=??)(ρ（4.31）
为了表达项?T ／?t ，我们将假设：在网格点上的T 值代表整个控
制容积上的值．于是：
)(01
P P e w t
t t
T T x c dtdx t
T c -?=???
?
?+ρρ （4.32） 对于k?T ／?x ．参考我们对稳态问题时所采用的措施，我们得到：
dt x T T k x T T k dxdt x T k x t t t w W P w e P
E e t
t t
e
w ??
??+?+??????---=????))()()()()(δδ
dt x T T k x T T k T T x c t
t t
w W P w e P E e P
P ?
?+??
?
?
??---=-?))()()()()(01δδρ （4.33） 现在我们需要一个有关
P T ，E T 和W T 如何随时间由t 到t+Δt 而变化的关系的假设．可以采取的假设有好几
种，其中的某一些假设可以用下式把它们归纳一般化：[]
t T f fT dt T P P t
t t P ?-+=?
?+01
)1( （4.34）
式中f 是一个在0与1之间变化的加权因子．对于
E T 和W T ，应用与上式相类似的公式，我们由方程(4.33)推出：
??
?
?????-----
???
?????---=-??))()()()()1())()()()()(00001
11101
w w e e w w e e P P x T T k x T T k f x T T k x T T k f T T t
x c W P P E W P P E δδδδρ （4.35）
在改写这个式子的时候，我们将把上标1去掉，并且记住此后
P T ，E T 和W T 代表T 在t+Δt 时刻的新值．结
果是：
[][
]
[]
00
0)1()1()1()1(P W E P W
W W E E E P P T a f a f a
T f fT a T f fT a T a ----+-++-+= （4.36）
式中：e
e E x k a )(δ=
w w W x k a )(δ= t
x c a P ??=ρ0
P W E P a fa fa a ++= (4.37)
4．3—2 显式，克兰克—尼科尔森(Crank —Nicolson)模式，以及全隐式模式
对于某个特定的加权因子f 的值，离散化方程可以简化为适用于抛物型微分方程的大家所熟悉的模式之一。特别是， f=0导致显式模式
f ＝0.5导致克兰克—尼科尔森模式 以及f=1导致全隐式的模式
我们将简单地讨论这些模式并且最后指出全隐式模式是我们所喜爱的一种模式．
[]
t T f fT dt T P P t
t t
P ?-+=?
?+01
)1(
不同的f 值可以由图4．5所示的
p T ~t 变化关系来说明．
显式模式实质上是假设老的0
P T
值代表除了时刻t+Δt 之外的整
个时间间隔上的P T 值．全隐式的模式则假说在时刻t ，P
T 的值突然由
0P T 降到1
P
T ，而在随后的整个时间步上保持为1P T
，于是在整个时间步期间内，温度为新值1
P T
所确定．克兰克—尼科尔森模式假设
P T 呈线性变化．乍看起来，线
性变化会比其它两种模式更加切合实际．那为什么我们会偏爱全隐式模式呢？答案很快就可以清楚．
对于显式模式(f=0)，方程(4．36)变为(4．38)
[][
]
[]
000
0)1()1()1()1(P
W E P
W
W W E E E P P T
a f a f a
T f fT a T f fT a T a ----+-++-+=（4.36）
[]
00
00P W E P W W E E P P T a a a T a T a T a --++=（4.38）
这意味着
P T 与其它未知量(如E T 或W T )无关，而直接可由已知的温度0P T ，0E T 及0
W T 得到．这就是为什
么把这种模式叫做显式的原因．
对于f≠0的任何模式都是隐式的．此时，
P T 将与未知量E T 和W T 有关，并且有必要解一组联立方程．
但是显式所具有的这一方便性为其本身所受到的一个严重限制所抵消．如果我们记住有关正系数的基本法则(法则
2)，并用它来检查一下方程(4．38)，那么我们就会发现0P T 的系数可能变为负值(我们把0
P T 看成是P T
在时间方向上
的一个相邻点值)．实际上，要使系数为正，时间步Δt 就必须足够小，以使0
P a
超过W E a a +．对于导热系数均
匀分布以及Δx=（δx ）e=（δx ）w 这个条件可以写作：
e
e
E x k a )(δ=
w w W x k a )(δ= t
x c a P ??=ρ0
(4.37) k x c t 2)(2?<
?ρ （4.39）
如果违背了这个条件；就可能出现物理上不真实的解，这是因为负的系数意味着：一个较高的0P
T 会导致一个较低
的T p 值的缘故．方程(4．39)是大家所熟知的有关显式模式稳定性的判定准则．
在我们需要通过减小间距Δx 来改进在空间的计算精度时，我们将不得不采用一个更小得多的Δt ．
克兰克—尼科尔森模式通常被描述成是无条件稳定的（对于f=0.5）。一个缺乏经验的使用者(用户)往往会就此认为，这就意味着不管采用多大的时间步，总可以取得物理上真实的解的．因此，这样一个使用者，一旦遇到解发生振荡时，就会感到惊讶．数学上的“稳定性’只能确保解的振荡最终将会消失，但是它并不能确保一定会给出物理上似乎合理的解．
对于f=0.5，方程(4．36)中的系数变为2/)(0
W E p
a a a +-．在均匀导热系数及均匀网格间距的情况下，这个系
数可以看成是ρcΔx ／Δt -k ／Δx 。这时，一旦时间步不够小，这个系数就可能变为负值，从而有可能产生物理上不真实的结果．
这种如图4．5所示的似乎合理的线性分布是一个只对小的时间间隔有效的好的温度—时间关系表达式．在时间间隔较大时，温度的固有的指数衰减关系相当于开始急剧下降，其斜率很陡，而后则是紧跟着一个平平的尾巴．在隐式模式中所作的假设于是要比克兰克—尼科尔森模式所采用的线性分布更接近于实际，尤其是对于大时间步的情况就更是如此．
如果我们要求方程(4．36)中的0
P T
的系数务必永不为负，唯一能够确保这一点的f 值应为1(当然令f ＞l 是没有意义
的)．于是全隐式模式(f ＝1)就能满足我们既简单而物理上又满意的要求．鉴于这一原因，在本书中，我们将采用全隐式模式．
必须承认：对于小的时间步，全隐式模式不如克兰克—尼科尔森模式精确，其原因仍可由图4.5看出：对于小的时间间隔，温度—时间曲线是接近于线性的．
探索一种结合这两种模式的优点，而不分担它们各自缺点的模式是诱人的，实际上这一工作已经有人做了，帕坦卡与巴利加(1978)已经描述了一种名叫指数模式的新模式．但是这种模式多少是太复杂了一些． 4.3—3 全隐式离散化方程
这里我们来讨论方程(4．36)的全隐式模式．
[][
]
[]
0000)1()1()1()1(P
W E P
W
W W E
E E P P T
a f a f a
T
f fT a T f fT a T a ----+-++-+=（4.36）
我们引入线性化源项．其结果是：
b T a T a T a W W E E P P ++= (4.40)
式中：e
e E x k a )(δ=
w w W x k a )(δ=
t
x c a P
??=ρ0
0P P c T a x S b +?= x S a a a a P P
W E P ?-++=0 可以看到，当Δt→∞时，方程简化为稳态的离散化方程。 全隐式模式的主要原则是
P T 的新值代表整个时间步上的值。因此如果导热系数P k 与温度有关，就应当反复由
P T 迭代算得新的P k 值。
稳态程序的其它环节，如边界条件、源项线性化处理以及TDMA 也都完全适用于不稳态的问题。 我们有关一维问题的详细讨论可以推广到三维与三维问题的条件． 4．4二维与三维问题
4．4—1 二维问题的离散化方程
二维网格的一小部分示于图4.6．
对网格点P 来说：
● 点E 与W 是它沿x 方向的相邻点
● 而N 与S(代表北方和南方)则是它沿y 方向的相邻点 ● 围绕着P 的控制容积用虚线表示 ● 该控制容积沿z 方向的厚度假设为1
在这里把图3.2中对距离Δx ，(δx)e 等引入的符号名称推广到二维的情况。
控制容积面相对于网格点的实际位置仍然是自由的．把控制容积面布置在两个网格点之间的中点，这是一种很显然的可能做法这里我们将推导可以适用于任何一种控制容积布置方式的离散化方程。
我们已经知道了怎样来计算P 与E 点之间的控制容积面上的热流密度
e q ，我们将假设所得到的e q 代表面积为Δy×1的整个表面上的
值．通过其它表面的热流密度可以用同样的方式得到．这样，微分方程：
S y
T k y x T k x t T c +????+????=??)()(ρ (4.42)
可以立即转变成离散化方程：
b T a T a T a T a T a S S N N W W E E P P ++++= (4.43) ?
式中：e
e E x y
k a )(δ?=
w w W x y k a )(δ?=
n n N y x k a )(δ?= s
s S y x k a )(δ?= t y x c a P ???=ρ0
0P P c T
a y x S
b +??=
y x S a a a a a a P P S N W E P ??-++++=0
乘积Δx Δy 是控制容积的体积。 4．4—2 三维问题的离散化方程
最后，我们再加进两个z 方向的相邻点T 和B(顶和底)以构成三维的网格图形．可以很容易地看出离散化方程是：
b T a T a T a T a T a T a T a B B T T S S N N W W E E P P ++++++= (4.45)
式中：e e E x z y k a )(δ??=
w
w W x z
y k a )(δ??=
n n N y x z k a )(δ??=
s s S y x z k a )(δ??= t
t T y y
x k a )(δ??=
b
b B y y x k a )(δ??=
t z y x c a P ????=ρ0
0P P c T a z y x S b +???= (4.46)
z y x S a a a a a a a a P P B T S N W E P ???-++++++=0
(4.46)
离散化方程中各个系数的物理意义：相邻系数
E a 、W a 、N a …、B a
的代表点P 与相应的相邻点之间的热导．项
00P P T a 是t 时刻控制容积内部所包含的内能(除以Δt)．常数项b 则由这一内能项与由C S 所造成的在控制容积内的发
热率所组成．中心点系数
p a
是所有的相邻点系数之和(包括“时间相邻”的0
P T 的系数
P a
)，并包括一项由线性的源项
所作的贡献．
4．4—3代数方程的解
应当注意：在我们构成离散化方程时，我们只把这些方程写成一种线性的形式而并没有设想一个求解这些方程的特定方法．因此我们在这里可以采用任何一种合适的求解方法．
把方程的推导与它们的求解方法看成是两个分开的步骤是有好处的．在具体选择应用时二者无需互相影响．在计算机程序中，我们可以很方便地把这两个步骤分成两个独立的段，其中任何一段都可以根据需要单独进行修改．
直接解法：至此我们已经由一维的离散化方程直接推广得到了多维的离散化方程．不能很容易地推广到多维问题中去的一段程序是三对角矩阵算法(TDMA)．
求解二维或三维问题的代数方程的直接解法(即那些无需迭代的方法)是非常复杂的，并且需要相当大量的计算机贮存量及时间．
对于一个只需要一次求解代数方程的线性问题，直接法可能是可以接受的，但是对于非线性问题，因为方程必须用最新的系数来重复求解，采用直接法往往是不经济的．因此在进一步讨论时，我们将把直接法排除在外． 达代法：代替直接解法的是解代数方程的达代法．
定义：从一个估计的因变量T 的场开始，再利用某种形式的代数方程求得一个改进的场．重复进行这一算法过程最后求得一个充分接近代数方程精确解的解。
优点：迭代方法通常只要求增加很少一点计算机存贮量就够了，对于处理非线性问题，迭代方法特别吸引人．
系数对于一个非线性的问题，在一组固定的系数值下把代数方程的求解过程一直进行到最终收敛是不必要也是不明智的．在这种情况下，在新的一组系数值约定之前应用已经给定的一组系数值对方程只要进行几次迭代就已经足够了。一般说来，似乎在计算系数所付出的代价与求解方程所得要花的时间之间应当存在某种平衡．在计算出一组系数值之后，我们必须进行充分的达代以求解方程，以便从这些系数身上得到足够的好处．但是如果把时间过多地花费在以这些只是暂时性的系数为基础来求解方程上则是不明智的．
迭代方法的种类：有许多求解代数方程的迭代方法．我们只介绍其中的两种方法：第一个方法可以看做是预备性的知识，而第二个方法则是我们将推荐采用的方法．
高斯—塞德尔(Gauss —Seidel)逐点计算法 在所有的迭代法中，最简单的一种方法是高斯—塞德尔法，其中按一定的顺序逐个访问每一个网格点，以计算那里的变量值．在计算机内只需要储存一组T 值．开始，这些值代表最初的估计值或是上一次迭代所得到的值．在访问每一个网格结点时，在计算机存贮器中相应的T 值交替改变如下，倘若把离散化方程写成
b T a T a nb nb P P +=∑ (4.47)
式中下标nb 代表一个相邻点，于是在被访问的网格点上的
P T 由下式算得：
P
nb P
a b
T a T nb +=
∑*
(4.48)
其中*
nb
T 代表在计算机存贮器中所存在的相邻点的温度值．对于那些在本次迭代过程中已经被访问过的相邻点，*
nb
T 是新鲜的计算；而对于那些尚待访问的相邻点，
*
nb T 是由前一次迭代所得到的值．其实在任何情况下，*nb
T 都是取的相邻点温度的最新值．当所有的网格点都按这种方式访问过一次之后，就算是完成了一次高斯—
塞德尔迭代．
用两个非常简单的例子说明这个方法．方程： T 1=0.4T 2+0.2 T 2=T 1+1 (4.49) 解：
可以发现，对这个例子，不管在开始计算时初始的估计值是多少，我们都可以得到精确的解．
迭代方法的一个有趣的特征是：在中间的计算步骤上，计算的精度可以是不很高的．由于中间的计算值只是简单地用作为下一次迭代的估计值，因此在中间过程中所作的近似计算．甚至于误差，在计算结束时都将趋于消失．由下面的例子，我们可以进一步加深对这种迭代方法的了解．
方程： T 1=0.4T 2+0.2 T 2=T 1+1 (4.49)
T 1=0T 2-1 T 2=2.5T 1-0.5 (4.50)
解：
这看起来是令人扫兴的．这里迭代过程已经发散了．更令人惊讶的是方程组(4．50)只是方程组(4．49)
的一种简单改写形式，而对于方程组(4．49)，我们在例1中已经得到了收敛的解．
于是，我们可以得出结论，高斯—塞德尔法并不总是可以得到收敛的解的．斯卡巴勒(Scarborough)(1958)已经推导出一个准则公式，当这个公式得以满足时，高斯—塞德尔法就一定可以得到收敛．现在我们不加证明地，把这个斯卡巴勒准则提出来，随后我们将讨论它的含义．
斯卡巴勒准则 高斯—塞德尔法收敛的充分条件是：
说明 (1)这个准则是个充分条件而并不是个必要条件．这就是说有的时候虽然我们违背了这个条件，但迭代仍能收敛．
(2)尽管我们并不主张采用高斯—塞德尔选代法，但是我们似乎还是希望我们的离散化方程应当满足斯卡巴勒准则，以使我们至少有一种迭代方法可以确保收敛。
(3)现在可以知道我们的某些基本法则是满足斯卡巴勒准则的(这些法则系由物理的观点出发，推导而得的)．
例如，负的
P S 的存在满足1/<∑p nb a a ．
我们也可以由我们的正系数要求这—准则说明． 如果某些系数是负的，那么
p a (通常等于∑nb a )可能会小于||∑nb a (因为在某些系数为负时，就会导致
||∑<∑nb nb a a )，从而违反这一准则．
t
y
x c a P
???=ρ0
0P P c T a y x S b +??= y x S a a a a a a P P S N W E P ??-++++=0 (4.41)
(4)当
p a 等于∑nb a ，以及所有的系数均为正时，对所有的方程我们均得到1||/||=∑p nb a a ．那么至少
有一个方程有可能会使∑||/||p nb
a a
小于1的地方究竟在哪儿呢？答案是在域的边界上．对于一个具有确定解的问
题，至少需要在一个边界点上规定温度值．在以该点作为相邻点的离散化方程中就包含有1||/||<∑p nb
a a 的意
思．这是因为在计算
||∑nb a 时(为应用斯卡巴勒准则)，应当只对那些未知的相邻点系数求和，而另一方面p a 则是
包括边界点系数在内的所有相邻点的系数之和．
高斯—塞德尔法的主要缺点是其收敛速度太慢，特别是当网格点数目很大时尤为明显．收敛速度慢的原因是易于理解的；本方法是以一次迭代传递一个网格间距的速率来传递边界条件所给予的信息的。
逐行法 现在我们可以很方便地把一维状态的直接法TDMA 与高斯—塞德尔法结合起来． ?
● 选择一条网格行(设在y 方向选取这样的网格行)，假定沿相邻的行(即所选择行上的点在x 与z 方向的所有相邻
点的集合)上的T 值由其最新值构成．
● 用TDMA 法求解所选行上的T 值．我们将在同一方向的所有各行上进行这种计算．如果想做的话，我们再按
同样的方法在其它方向重复上述的程序．
●
b T a T a T a T a T a S S N N W W E E P P ++++=
讨论 (1)某一选定的行网格点上的离散化方程中包含着沿两相邻行上网格点(用叉线表示)的温度．如果这些温度用它们
的最新值取代，那么沿选定行网格点(图中用圆点表示)的方程就可以看成是一些一维方程，这些方程可以采用TDMA 法来求解。按照这一方法对所有的沿y 方向的其它各行进行计算，随后再按类似的方法处理在x 方向的所有各行．
(2)逐行法的收敛速度比前面所提到的逐点法快，这是因为不管沿行有多少个网格点，在行两个端点上的边界条件信息都可以立即传递到域的内部．而沿另一方向的信息传递速率则类似于逐点法的传递速率．
(3)通过交替改变应用逐行迭代的行方向，我们就可以把所有边界上的信息迅速传递到中间区域．
(4)几何形状以及其它状态特性往往会导致，例如，y 方向的系数远远大于x 方向系数(见图4．8)的情况。
w w W x y k a )(δ?=n
n N y x
k a )(δ?=
s s S y x k a )(δ?= (4.44) 在这种情况下，当对y 方向应用TDMA 逐行迭代时，收
敛得特别快．这是因为需要代入离散化方程的沿相邻行上的温度估值对于离散化方程的影响不大．
(5)除了应用TDMA 所选择的行方向之外，扫描方向(即选择计算行的先后次序)在某些情况下也是重要的．
对于如图4．9所示的边界条件，从左到右的扫描(即，选择域的左边界作为第一行，而后逐渐向右移动行)会把左边界上已知的温度传递到计算域内部；与此相反，由于右边界没有给定温度，自右而左的扫描不可能带进这种有用的信息(同样的考虑方法适用于逐点计算中访问点的先后次序)．当有对流存在是，扫描方向是尤为重要的．十分清楚，由上游向下游扫描的收敛速度要比按相反方向扫描的收敛速度快得多．
其它的一些迭代方法 皮斯曼(Peaceman)和拉什福德(Rachford)(1955)引入了一种常用的、叫做ADI(方向交替的隐式)的逐行求解法．另外一种用于求解多维离散化方程的迭代解法是由斯通(Stone)(1968)所述的强隐式法(SIP)．我们把这些方法的详细研究留给那些感兴趣的读者自己去做．
4．5 超松弛与欠松弛
在代数方程的迭代求解过程中，或是用于处理非线性问题的整体迭代模式中，人们往往希望加快或是减慢前后二次迭代之间因变量值的变化．依因变量的变化究竟是被加速还是被减慢而定，这一过程被称之为超松弛或是欠松弛．
超松弛常常用于和高斯—塞德尔法相结合，这种结合起来的模式是所谓的持续超松弛(SOR)。通常超松弛很少与逐行法结合起来使用．
欠松弛对于非线性的问题是十分有用的工具．在强烈非线性方程组的迭代求解过程中，往往采用欠松弛来避免发散． 引进超松弛或是欠松弛的方法有许多种．这里将描述其中的几种具体做法．我们将针对下列形式的通用离散化方程进行讨论．
b T a T a nb nb P P +=∑ (4.52)
此外，取
*p T 作为前一次迭代所得到的p T 值．
采用一个松弛因子
方程(4．52)可以写成
P
nb
P
a b
T a T nb +=
∑ (4.53)
如果我们在方程的右侧加上*p
T ，再减它，我们就有：
???
? ??-++=∑**
P nb P T a b T a T T P nb P (4.54) 式中括号内的部分代表由本次迭代所产生的
p T 的变化．这一变化可以通过引进一个松弛因子加以修改，所以
???
? ??-++=∑**
P nb P T a b T a T T P nb P α (4.55a) 或 *)
1(P nb T a b T a T a P
nb P P α
αα-++=∑ (4.55b)
首先，应当注意到当迭代收敛时
p T 变成与*
p T
相等．方程
(4．55a)意味着T 的收敛值确实是满足原始方程(4．52)的．当然任何一种型式的松弛都必须满足这一性质；不管是通过使用任意的松弛因子或是任何类似的手段，所得到的最终收敛解都必须满足原始的离散化方程．
α值的选取、欠松弛与超松弛
当方程(4．55)中的松弛因子α在0与1之间时，它的作用是欠
松弛的；即p T 的值更接近于*p T 一些．对一很小的α值，p T 的
变化变得很慢．而当α大于1时，就产生超松弛．
不存在什么选取最佳α值的一般法则．α的最佳值与许多因素有关．诸如所研究问题本身的特性，离散网格点的数目，网格点之间的间距以及所采用的迭代方法等，都是重要的影响因素．通常，可以根据经验以及对所给定的问题所作的试探性计算求得一个合适的α值．
在整个计算期间都保持采用相同的α值是没有必要的．在不同的迭代次数之间可以改变α的值。事实上对每一个网格点选用各自不同的α值。尽管这样做不太方便，也是允许的．
通过惯量进行松弛 进行超松弛或欠松弛的另一种方法是用下面的公式来代替离散化方程(4．52)：
*)(P iT b T a T i a nb nb P P ++=+∑ (4.56)
式中i 是所谓的惯量．对于正的i 值，方程(4.56)具有欠松弛的作用，而负的i 值则产生超松弛．
同样，这里也不存在获得最佳i 值的一般法则；这个最佳值必须由对一特定问题所取得的经验确定。由方程(4．56)我们可以推论：i 应当与
P a 不相上下，同时i 的模值愈大，松弛的作用也愈烈。
有时一个稳态问题的解是通过应用相应的不稳态问题的离散化方程求得的．这样，其中的“时间步”就变成与我们前
面的“迭代”的概念相同．“老的”值0
p T
就只不过是代表上一次迭代的值
*p T 。
在这个意义上，方程(4．46h)中的项0
0p p T
a 与方程(4．66)中*p iT 项起着同样的作用．于是，惯量i 类似于不稳态公式中的系数0
p a
．这样的类比关系为我们提供
了一种确定合理的i 值的方法．另一方面，现在可以简单地认为通过不稳态公式解稳态问题的做法是一种特定类型的欠
松弛．其中所选取的时间步愈小，所得的欠松弛也愈强．顺便提一句，一个负的时间步Δt 值会产生超松弛．
4．6某些几何上的考虑
4.6—
1
控制容积面的位置
到目前为止，还没有专门描述过应在何处放置控制容积面的具体做法问题．因为离散化方程是在一般的条件下推得的，所以这些离散化方程适用于任何特定的控制容积面构成方法．
在许多可能的构成方法中，我们将讨论其中两种不同的替代形式，并讨论它们各自有关的优点．这两种方法将称之为方法A 与方法B ．为方便起见，描述将针对一个两维的问题，尽管所包含的有关概念也同样适用于一维和三维的情况．
方法A ：控制容积面放在两个网格点之间的中点 最显而易见的构成控制容积的方法是把它们的面放在相邻网格点之间的中点．这种布置方式示于图4．10，其中虚线表示控制容积面．图中故意把网格画得很不均匀；由图中可以看到，这样做所得的一个结果是：一个典型的网格P 并不落在包围该点的控制容积的几何中心上．
方法B ：网格点放在控制容积的中心 另一种方法(如图4．11所示)是先画控制容积，而后把网格点放在控制容积的几何中心．按照这种布置方式，如果控制容积的尺寸不均匀，则其表面就不会落在两个网格点之间的中点上．
讨论 （1）应当注意对均匀的网格(或是均匀的控制容积尺寸)两种方法趋于一致．因此对两种方法的比较只是在非均匀的网格间距的情况下才是有意义的．
(2)在方法A 中把控制容积面放在两个网格点之间的“中点”，这为计算通过该面的热流提供了较高的精度．如在3．4节中所提到的，分段线性分布的斜率与在两个网格点之间的中点所算得的任何一根抛物线的斜率相等．
(3)另一方面，在图4．10中网格点P 可能不能在控制容积的几何中心，这一事实却是一个缺点．
于是不能认为在计算源项、导热系数以及类似的其它一些量时
p
T 是对控制容积的一个好的代表值．此外，甚至在计
算控制容积面上的热流密度时，方法A 也不是无可非议的．例如图4．10中的e 点并不位于它所在的控制容积面的中点．这样，把e 点的热流密度看成是代表整个控制容积面上的值的做法也会造成某些误差．
(4)方法B 就没有这些缺点．因为根据定义，P 点位于控制容积的中心，同时象e 这样一些点也位于它们自己所代表的面的中心(见图4．11)．但是控制容积面并不位于相邻网格点之间的中点．因此，与方法A 不同，方法B 并不具有前面所提到过的抛物线的偶然性质所具有的好处．
(5)或许，方法B 的决定性的优点还在于它的方便性．因为控制容积是迄今为止所提出的离散化方法的基本单位，先画控制容积的边界而后决定网格点的位置是比较方便的，例如对一组合的固体，我们可以把控制容积面放在材料性质发生不连续的地方(见图4．12)．
类似地，边界条件的不连续性也易于处理．如果边界的一部分是绝热的，其余的部分是等温的，我们便可以这样来设计控制容积，以避免在一个控制容积面内存在不连续性；这种情况也示于图4．12中．在方法A 中，因为我们必须首
先规定网格点的位置，要想把控制容积面布置在所希望放的位置就要困难得多；
(6)需要对计算域的边界附近的控制容积的设计作些补充考虑．如在图4.13中所表示的那样，方法A 在边界的网格点周围造成“半”控制容积（在4
．2
—6小节中引入)．
对方法B 来说，以规则的控制容积填满计算域、并把边界的网格点放在邻近边界的控制容积的面上，是比较方便的．这样的布置方式示于图4．14中．
一个典型的边界面不是放在边界点B 与内点I 之间，而实际上是通过边界点．如果在边界点周围想象存在一个零厚度的控制容积，那么面i 相对于网格点B 和I 的位置可以认为是与方法B 的正常情况相一致的，按照这种布置方式，也就没有必要对邻近边界的控制容积采用特殊的离散化方程，现有的边界条件数据(如已知温度或热流密度)就可以直接用于边界面i 上。
4．6—2 其它坐标系
至此我们已经推导出了应用于直角坐标系中一个网格的离散化方程。在本书的余下部分，将继续采用同样的坐标系来处理几乎所有问题．这种坐标系表达方便而且易于理解．但是所提出的这种方法不只限于直角坐标系，它还可以用于任意一种正交坐标系．为了说明在其它坐标系中离散化方程的推导，我们将讨论用r 与θ所表示的二维极坐标问题．
与方程(4．42)相应的r θ形式是：
S T r k r r T rk r t T c +????+??=??)(1)(1θ
θρ (4.57)
在r θ坐标系中的网格与控制容积示于图4．15中．设控制容积在z 方向的厚度为1．为了得到这一坐标系下的离散
化方程，我们以r 乘方程(4．57)的两边，并在整个控制容积的范围内对r 与θ进行积分(这一运算代表体积积分，因为r dr r dθ代表单位厚度的体积元)．按照与4．4—1小节相同的步骤，我们得到离散化方程：
b
T a T a T a T a T a S S N N W W E E P P ++++=
(4.58)
式中：
e
e e E r r k a )(δθ?=
w
w w W r r k a )(δθ?=
n
n n N r r k a )(δθ?=
s
s s S r r k a )(δθ?=
t
V
c a P ??=
ρ0
0P P c T a V S b +?= V S a a a a a a P P S N W E P
?-++++=0 (4.59) 这里ΔV 是控制容积的体况，它等于
r r r s n )(5.0??+θ (应当注意ΔV 没有必要等于r r p
??θ，除非
P 位于n 与S 之间的中点)．
上述说明表明由一新的坐标系引入的补充特性主要是几何上的特性．只要所需要求的长度、面积以及体积计算得合适，那就不需要什么新的原则．于是在任何正交坐标系中的离散化方程可以按照同样的方法推导得到．但是如果所用的只是由两个网格点定义的分布曲线的话，那么正交性的要求就是必不可少的．在图4．15中，控制容积面e 是与线PE
相垂直的，正是这一事实才使我们可以只用
P T 与E T 来计算通过控制容积面的热流密度．对于非正交的网格，就要采
用较为复杂的离散化方程．
在本书的其余部分，我们将只采用直角坐标系来进行全部的代数推导．但是当几何形状发生明显的变化时，整个处理方法将同样适用于任何一种正交坐标系．
第五章 对流与扩散
5．1 任 务 （略 见课本P91） 5．2 一维稳态对流与扩散 对可能的最简单情况的讨论
讨论只有对流与扩散这两项存在的情况下的一维稳态问题．
控制微分方程是：
)()(dx
d dx d u dx d φφρΓ= （5.4） 式中u 代表x 方向的速度． 连续性方程：
0)(=u dx
d
ρ 或 ＝常数u ρ （5.5）
控制容积和控制容积面的位置
为了推导离散化方程，我们将应用如图5．1所示的三网点群．尽管控制容积面e 和w 的实际位置并不会影响我们的最终公式，然而假设e 位于P 和E 之间的中点，w 位于W 与P 之间的中点还是比较方便的．
5．2—1 预备性的推导
在图5．1所示的整个控制容积内对方程(5．4)进行积分给出：
)()(dx d dx d u dx d φ
φρΓ= ??Γ=e
w e
w
dx dx d dx d dx u dx d )()(φφρ w e w e dx
d dx d u u )()()()(φ
φφρφρΓ-Γ
=- (5.6) Φ的分布 在上一章中我们已经知道如何由对Φ的一个分段线性分布来表示项 x /d φd Γ。很自然对于对流项首
先想到的也是采用同样的分布曲线．其结果是：)(2
1
P E e φφφ+= )(2
1
W P w
φφφ+= (5.7) 因子1/2 出自于界面位于中点的假设；对不同的界面位置则要采用其它的内插因子．现在，方程(5．6)可以写成：
w
W P w e P E e W P w P E e x x u u )()()()()()(21)()(21δφφδφφφφρφφρ-Γ--Γ=+-+ (5.8) 式中
e Γ与w Γ的值要用4．2—3小节中所推荐的方法求得(这一点适用于全书，尽管在前几节中可能并没有重复这
些说明)．
为了把方程写得更紧凑，我们定义两个新的符号F 与D 如下：u F ρ≡ x
D δΓ≡
(5.9)
因次与物理意义
二者具有同样的因次，F 表示对流(或流动)的强度，而D 是扩散传导性。应当注意，与D 总是保持正值相反，F 既可以取正值也可以取负值，依流动方向而定．应用这些新的符号，离散化方程变为：
W
W E E P P a a a φφφ+= (5.10)
式中：2
e
e E F D a -
=
2
w
w W F
D a +=
)(2
2w e W E w
w e e P F F a a F D F D a -++=-++= （5.11）
讨论 (1)由于连续性
e F =w F ，我们确实得到W E P a a a +=这一性质．进一步注意到：根据方程(5．11c)，
离散化方程只是在流场满足连续性时才具有这一性质．
(2)离散化方程(5．10)隐含着Φ分段线性分布的含义．这种形式也就是熟知的中心差分格式，并且是泰勒级数公式的自然结果．
不可接受的离散化方程
(3)考虑一个简单的例子是有启发的，其中
1==w
e D D 及4==w e
F F 12
-=-
=e
e E F D a 32=+=w w W F D a
2)(2
2=-++=-++
=w e W E w w e e P F F a a F
D F D a 此外，如果
E φ和
W φ
的值已如，我们就可以由方程
W
W E E P P a a a φφφ+= (5.10)
W E P φφφ32+-=
(a)如果E φ＝200及W φ＝100，结果是p φ=50 (b)如果E φ＝100及W φ＝200，结果是p φ＝250
因为
p φ实际上不能落在由其相邻点值所建立起来的范围100---200之外，这些结果显然是不真实的．
(4)实际上也许我们早就可以预料到这些不真实的结果，因为方程 2
e
e E F D a -= 2w w W F D a +=系数有时可
能会变为负值．当|F|超过2D 时，与F 是正或是负有关，存在着
E a 或W a 成为负值的可能．这将违背我们的基本法则
之一，从而产生灾难性的结果．
(5)此外，负系数意味着)(w e W E P
F F a a a -++=小于∑|a nb |，从而违反斯卡巴勒准则．这样，离散化方程的
逐点解就可能发散.
这就是为什么所有以前用中心差分格式来求解对流的努力都只能限于低雷诺数(即低的F ／D 值)的原因． (6)在扩散项为零的情况下(即Γ＝0)，
2
2w w e e P F
D F D a -++
=中心差分格式导致
P a =0．于是方程(5．10)就变得不适于用逐点法来求解了，并且
图5.1 一维问题的典型网格节点群 W P E
x w e w x )(δe
x )(δx
?
图5.1 一维问题的典型网格节点群
W
P
E
x
w
e
w
x )(δe
x )(δx
?也不适于采用大多数其它的迭代解法．
由于上述预备性公式推导已经导致一个不可接受的离散化方程，因而必须寻找更好的公式．在下面的几个小节中将描述几个这样的可能方案．
5．2—2 上风方案
众所周知的一帖用于解决前面所遇到的困难的妙药是上风方案，它也叫做上风差分格式、上游差分格式以及供体(施主)室法等．
这个方案首先由库兰特(Courant)、伊萨克逊(Issacson)和里斯(Rees)(1952)提出来，随后又由金特里(Gentry)、马丁(Martin)和戴利(Daly)(1966)，巴拉卡特(Barakat)和克拉克(Clark)(1966)以及朗切尔(Runchal)和沃尔夫茨坦〔Wolfshtein)(1969)等人再次阐明．
预备性公式的弱点
上风方案认为预备性公式的弱点在于假设：界面上的对流性质
e φ是E φ和P φ的平均值
上风方案提出了一个较好的处理办法．保留扩散项的公式不变，而对流项则按下列假设计算： 界面上的Φ值等于界面上风侧网格点上的Φ值．于是：
P e φφ= 如果 0>e F (5.12a)
E e φφ= 如果 0 F (5.12b) 按类似的方法可以确定 w φ值． 如果我们定义一个新的运算子的话，条件语句(5．12)可以写得更紧凑．我们特定义『A ，B 』代表A 和B 中的大者．于是，上风方案意味着： [][]0,0,e E e P e e F F F --=φφφ (5.13) w e w e dx d dx d u u )()()()(φ φφρφρΓ-Γ=- (5.6) 用这个概念代替方程(5．7) )(2 1 P E e φφφ+= )(21W P w φφφ+= (5.7) [][]0,0,w P w W W w F F F --=φφφ [][][][]0,0,0,0,)()(w P w W e E e P w W e e w e F F F F F F u u -+---=-=-φφφφφφφρφρ w W P w P E e e D D dx d dx d )()()()(φφφφφ φ---=Γ-Γ 就可把离散化方程写成 W W E E P P a a a φφφ+= (5.14) 式中[] 0,e e E F D a -+= []0,W w W F D a += [][])(0,0,w e W E w w e e P F F a a F D F D a -++=-+++= 计论 (1)很显然，方程(5．15)不会产生负的系数．于是，解将总是在物理上真实的；同时，斯卡巴勒准则也将得到 满足． (2)上风方案所体现的主导思想的理论基础到底是什么呢？一个关于上风方案的、明白易懂的物理图象是建立在 “槽与管”的模型的基础上[戈斯曼、潘、朗切尔、斯波尔丁和沃尔夫茨坦(1969)]． 图5.2 槽和管的模型 如图5．2所示，可以把控制容积看成是一系列由短管连接起来的搅拌槽．通过管子的流动代表对流，而通过槽壁的传导代表扩散．因为槽是搅拌的，每一个槽中都含有温度均匀的流体．于是可以认为流过每一连接管的流体具有在上游一侧槽内流体的温度． 5．2—3精确解 很幸运，如果Г取作常数，就可以准确地求解控制方程（5．4)，如方程(5．5)所给定，ρu 已经是常数． )()(dx d dx d u dx d φ φρΓ= （5.4） 0)(=u dx d ρ 或 ＝常数u ρ （5.5） 如果采用一个城0≤x≤L ，并利用边界条件： 在 0=x 处，0φφ= （5.16a ） 在L x =处，L φφ= （5.16b ） 方程(5．4)的解是 1)ex p(1 )/ex p(00--= --P L Px L φφφφ (5.17） 式中P 是由下式定义的贝克列(Peclet)数 Γ ≡ uL P ρ (5.18) 显然P 是对流与扩散的强度之比。 精确解的性质可以由图5．3了解．在该图上，绘出了不同的贝可列数时的Φ~x 的变化。 在贝可列数为0的极限条件下，我们的问题成了纯扩散（或热传导）的问题，并且Φ~x 的变化是线性的。 当沿正x 方向流动时（即P 值为正），域内的Φ值看来受上游的影响要大些。对于大的正P 值，该域的大部分范围内Φ值非常接近于上游的值Φ0． 当P 为负值时，图象正好相反．如果流体沿负x 方向流动，则ΦL 成为上游的值，它支配着域内的Φ值．当P 为相当大的负值时，在该域的大部分范围内Φ值非常接近于ΦL 。 含意 为了构成离散化方程，我们现在可以从图5.3得到指导，该图描述了网格点之间的适宜的Φ---x 分布．(1)由图可以很容易地明白为什么我们的原来的推导不能给出一个满意的公式．除非 |P| 值非常小之外，Φ--x 曲线均偏离线性很远．(2)当|P|值很大时，x ＝L ／2(界面)处的Φ值几乎等于上风边界上的Φ值．这精确地符合上风方案的假设， 上风方案的缺点 （1）但是在上风方案中把这种假设用于所有的|P|，而不只是用于大的|P|值。 (2)当|P|值很大时，在x ＝L ／2处，d Φ/d x 几乎为0。这样扩散也就几乎不存在．上风方案总是由一线性的Φ--x 分布计算扩散项，从而在大的|P|值条件下过高地估计了扩散值。 以精确解为基础得到离散化方程---指数方案 要是直接由图5．3所示的精确解为基础得离散化方程，那么所得到的方案本来是不会有这样的缺陷的。让我们来着手推导这样的方案，我们将把这一方案称作为指数方案．指数方案是建立在由斯波尔丁(1972)首先提出的公式基础之上的，并且是雷思比(Raithby)和托兰斯(Torrance)(1盯4)所推荐并应用的方案之一． 5．2—4 指数方案 总流量密度J 考虑一个由对流流量密度ρ u Φ与扩散流量密度dx d /φΓ-所组成的总流量密度J 是有益的。于是： dx d u J φ φρΓ -= （5.19） 按照这个定义，方程5.4 )()(dx d dx d u dx d φφρΓ=变为：0=dx dJ （5.20） 在图5.1所示的整个控制容积内积分方程（5.20），就得到 0=-w e J J （5.21） 总流量密度J 的分布 现在，精确解(5．17)可以用作为点P 与E 之间的分布，其中用P φ和E φ代替0φ和L φ， 并用距离(δX)e 代替L ．将这个分布关系代入方程(5.19)，就可以给出对e J 的表达式 )1)ex p((--+=e E P P e e P F J φφφ （5.22） 式中：e e e e e e D F x u P =Γ= )()(δρ （5.23） e F 和e D 由方程(5．9)定义． u F ρ≡ x D δΓ≡ (5.9) 应当注意 e J 与点P 和E 之间界面的位置无关．最后，把方程(5．22)和一个关于w J 的类似关系代入方程(5．21) 就得到：0)1 )exp(()1)exp((=--+---+w P W W e E P p e P F P F w φ φφφφφ （5.24） 这个方程可以写成标准形式： W W E E p P a a a φφφ+= (5.25) 式中：1 )/exp(-= e e e E D F F a 1)/e x p ()/e x p (-=w w w w w W D F D F F a )(w e W E P F F a a a -++= （5.26） 方案的优点 这些系数表达式定义了指数方案．在应用于一维的稳态问题时，该指数方案保证对任何的贝克列数以及任意数量的网格点均可以得到精确解． 未得到广泛的应用的原因 尽管这种方案具有如此极为理想的性质，但由于下列原因这种方案仍然未得到广泛的应用： (1)指数运算是费时的 (2)对于二维或三维的问题，以及源项不为零的情况，这种方案是不准确的． 看来没有理由在计算指数上花费过多的时间． 实际上需要的方案 我们实际上需要的是一种既具有指数方案定性特征又易于计算的方案．下面我们将介绍两种这样的方案并推荐采用其中的第二个． 5．2—5 混合方案 混合方案是由斯波尔丁(1972)制定的；这一方案也出现在帕坦卡与斯波尔丁(1970)的书中名为“高横向流的修正”的条目下． 无因次形式a E /D e 与贝克列数之间的曲线关系 为了评价指数方案与混合方案之间的关系，我们将画出系数 E a ，或它的无因次形式e E D a /与贝克列数之间的 曲线关系．由方程(5．26)我们推得： 1)exp(-= e e e E P P D a (5.27) e E D a /随e P 的变化如图5.4所示． 对正的e P ，网格点E 是下游的相邻点，看来它的影响随着e P 的增加 而减小。 当 e P 为负值时，E 是上游的相邻点并且具有比较大的影响． 由图5．4所示的实线可以看到e E D a /准确变化的某些特殊性质： 对∞→e P 0→e E D a ； 对-∞→e P e e E P D a -→ ； 对0=e P 2 1e e E P D a -= （5.28） 代表这些极限情况的三条直线也示于图5．4中．可以看到这三条直线构成准确曲线的一根包络，并代表着这一准确 曲线的合理近似，混合方案实际上就是由这三条直线所组成，所以 对2- E P D a -=； 对22≤≤-e P 2 1e e E P D a -=； 对2>e P 0=e E D a （5.29） 应用特殊的符号『 』把这些表达式合并成一个紧凑的形式，该符号代表其内所包含的所有量的最大值．于是 ]]0,21,[[e e e E P P D a --= （5.30a ） 或 ]]0,2 ,[[e e e E F D F a --= (5.30b) 混合方案的意义可以通过观察下面两点来了解： (1)在贝克列数为一2≤Pe ≤2时，混合方案同中心差分格式一致 (2)在该范围之外，混合方案简化为上风方案，其中已把扩散项置为0，所以混合方案不再具有5．2—3小节末尾所列举的一些缺点． “混合”这个名称就是代表着中心差分格式与上风方案的组合，但是最好还是把混合方案看成是如图5．4所示的准确曲线的三条近似线． 现在可以把混合方案的对流—扩散离散化方程写成：W W E E p P a a a Φ+Φ=Φ （5.31） 式中：]0,2 ,[e e e E F D F a - -= ]0,2,[w e w W F D F a += )(w e W E P F F a a a -++= 应当记住这个公式适用于界面位于两个网格点之间的任何一个任意位置的情况，它不只限于中点界面的情况． 5.2—6幂函数方案 由图5．4，可见在 e P =± 2时，混合方案偏离准确曲线相当大，看来在|P|值一超过2就马上把扩散效应置为0的 做法多少有点过早了．准确曲线的一个更好的近似系由幂函数方案给定，该方案已由帕坦卡(1979a)作了描述．尽管幂函 数方案较混合方案要复杂一些，但是幂函数表达式在计算时并非是特别费时的，而这些表达式却极好地代表着指数状态．E a 的幂函数表达式可以写成如下： 对于10- P e e E P D a -=； 对于010<≤-e P ()e E e E P P D a -+=5 1.01； 对于100≤≤e P ()5 1.01E e E P D a -=； 对于10>e P 0=e E D a （5.33） 把这些方程与方程(5．29)相比较，我们发现：对于 10||>e P ，幂函数方案即与混合方案一致．方程(5．33)的 一种紧凑形式可以写作 ]],0[[]])1.01(,0[[5e e e e E F D F D a -+- = （5.34） 可以由表5．1来评判幂函数与准确的指数方案的接近程度； 表5.1由幂函数方案与指数方案所给定的系数值的比较 P e a E /D e P e a E /D e 幂函数方案 指数方案 幂函数方案 指数方案 -20 -10 -5 -4 -3 -2 -1 -0.5 0 20.00 10.00 5.031 4.078 3.168 2.328 1.590 1.274 1 20.00 10.00 5.034 4.075 3.157 2.313 1.582 1.271 1 0.5 1 2 3 4 5 10 20 0.7738 0.5905 0.3277 0.1681 0.07776 0.03125 0 0 0.7707 0.5820 0.3130 0.1572 0.07463 0.03392 0.00045 4.1?10-8 两者之间的差异是太小了，以致很难用图形来进行比较．如前面所述，尽管对许多实际问题来说混合方案用起来一样好，在本书中还是把幂函数方案作为推荐采用的对流—扩散公式． 5．2—7 一个通用化的公式 为了进一步了解对流—扩散公式并构成一个可以把迄今所考虑过的所有方案都包括进去的通用格式，我们现在来探讨一下有关系数的某些一般性质．让我们来讨论图5．5所示的由距离δ分开的网格点i 和i 十1．我们的兴趣在于表达通过这两个网格点之间的界面上的总流量密度J ．利用方程(5．19)， dx d u J ΦΓ-Φ=ρ （5.19） 我们写出：)/(* δδx d d P J J Φ-Φ=Γ≡ （5.35） 式中P 是贝克列数，ρuδ／Г．在界面上的Φ值将是Φi 与Φi+1的某个加权平均值，1* +Φ-Φ=i i A B J (5.37) ()1*+Φ-Φ=Φ-i i i A P J (5.44) ()11*++Φ-Φ=Φ-i i i B P J (5.45) 图5.5 两个网格点之间的总流量密度J ]]0,[[)()(P P A P A -+= （5.42） ]]0,[[)()(P P A P B += (5.43) 如果我们现在应用流量关系式(5．37)于界面e 和w ，并利用方程(5．42)和(5．43)，就可以得到下列通用的对流—扩散公式： W W E E p P a a a Φ+Φ=Φ (5.46) 式中 ]]0,[[)(e e e E F P A D a -+= （5.47a ） ]]0,[[)(w w w W F P A D a += （5.47b ） )(w e W E P F F a a a -++= （5.47c ） 现在可以把前面所推得的各种方案看成是仅仅选择不同的函数A(|P|)而已．前面所考虑的各种方案的A(|P|)表达式列于表5．2并用图线示于图5．7中，可以通过与相应的精确解相比较的方法来评判每一种函数的满意程度． 表5.2 各种不同方案（格式）的函数A （|P|） 方案（格式） 对A （|P|）的公式 中心差分 上风 混合 幂函数 指数（精确解） 1-0.5|P| 1 [][]||5.01,0P - [][]5 |)|1.01(,0P - |P|/[EXP(|P|)-1] 5．2—8 各种方案(格式)的结果 在结束一维问题的讨论之前，我们要检查一下对于一定的 E φ和W φ值由各个方案所算出的P φ值．让我们令值 E φ＝1及W φ=0(这样做并不失其一般性)．进一步令距离(δx)e 与(δx)w 相等，于是P φ将是P （≡ρuδx ／Г)的一个函数．对 各个P 值由不同的方案所算得的 P φ值示于图5．8中．(由幂函数方案所得到的结果与指数方案(精确解)两者是太接近了， 它们根本无法分成两条线．) 除了中心差分格式之外，所有方案都给出某种物理上真实的解；与众不同，中心差分格式 则得出某些位于由边界值所规定的0—1范围以外的值． 因为决定数值方案状态的是网格的贝克列数，原则上是有可能通过分细网格(即：采用较小的δx)直到P 足够小(＜2)、从而使中心差分格式获得真实的解的．但是在大多数的实际问题中，这种对策需要极细的网格．这样做在经济上往往是 第四章循环流化床锅炉炉内传热计算 循环流化床锅炉炉膛中的传热是一个复杂的过程，传热系数的计算精度直接影响了受热面设计时的布置数量，从而影响锅炉的实际出力、蒸汽参数和燃烧温度。正确计算燃烧室受热面传热系数是循环流化床锅炉设计的关键之一，也是区别于煤粉炉的重要方面。 随着循环流化床燃烧技术的日益成熟，有关循环流化床锅炉的炉膛传热计算思想和方法的研究也在迅速发展。许多著名的循环流化床制造公司和研究部门在此方面也做了大量的工作，有的已经形成商业化产品使用的设计导则。 但由于技术保密的原因，目前国内外还没有公开的可以用于工程使用的循环流化床锅炉炉膛传热计算方法，因此对它的研究具有重要的学术价值和实践意义。 清华大学对CFB锅炉炉膛传热作了深入的研究，长江动力公司、华中理工大学、浙江大学等单位也对CFB锅炉炉膛中的传热过程进行了有益的探索。根据已公开发表的文献报导，考虑工程上的方便和可行，本章根椐清华大学提出的方法，进一步分析整理，作为我们研究的基础。为了了解CFB锅炉传热计算发展过程，也参看了巴苏的传热理论和计算方法，浙江大学和华中理工大学的传热计算与巴苏的相近似。 4.1 清华的传热理论及计算方法 4.1.1 循环流化床传热分析 CFB锅炉与煤粉锅炉的显著不同是CFB锅炉中的物料(包括煤灰、脱硫添加剂等)浓度C p 大大高于煤粉炉，而且炉内各处的浓度也不一样，它对炉内传热起着重要作用。为此首先需要计算出炉膛出口处的物料浓度C p，此处浓度可由外循环倍率求出。而炉膛不同高度的物料浓度则由内循环流率决定，它沿炉膛高度是逐渐变化的，底部高、上部低。近壁区贴壁下降流的温度比中心区温度低的趋势，使边壁下降流减少了辐射换热系数；水平截面方向上的横向搅混形成良好的近壁区物料与中心区物料的质交换，同时近壁区与中心区的对流和辐射的热交换使截面方向的温度趋于一致，综合作用的结果近壁区物料向壁面的辐射加强，总辐射换热系数明显提高。在计算水冷壁、双面水冷壁、屏式过热器和屏式再热器时需采用不同的计算式。物料浓度C p对辐射传热和对流传热都有显著影响。燃烧室的平均温度是床对受热面换热系数的另一个重要影响因素。床温的升高增加了烟气辐射换热并提高烟气的导热系数。虽然粒径的减小会提高颗粒对受热面的对流换热系数，在循环流化床锅炉条件下，燃烧室内部的物料颗粒粒径变化较小，在较小范围内的粒径变化时换热系数的变化不大，在进行满负荷传热计算时可以忽略，但在低负荷传热计算时，应该考虑小的颗粒有提高传热系数的能力。 炉内受热面的结构尺寸，如鳍片的净宽度、厚度等，对平均换热系数的影响也是非常明显的。鳍片宽度对物料颗粒的团聚产生影响；另一方面，宽度与扩展受热面的利用系数有关。根 数值计算大作业 一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。物体的导热系数λ为1.0w/m·K。边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K; 要求： 1、写出问题的数学描述； 2、写出内部节点和边界节点的差分方程； 3、给出求解方法； 4、编写计算程序(自选程序语言)； 5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图； 6、就一个工况下（自选）对不同网格数下的计算结果进行讨论； 7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论； 8、对4个不同表面传热系数的计算结果进行分析和讨论。 9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。（自选项） 1、写出问题的数学描述 设H=0.1m 微分方程 22220t t x y ??+=?? x=0，0 y=H ，0 《传热学》试题库 第一章概论 一、名词解释 1．热流量：单位时间内所传递的热量 2．热流密度：单位传热面上的热流量 3．导热：当物体内有温度差或两个不同温度的物体接触时，在物体各部分之间不发生相对位移的情况下，物质微粒(分子、原子或自由电子)的热运动传递了热量，这种现象被称为热传导，简称导热。 4．对流传热：流体流过固体壁时的热传递过程，就是热对流和导热联合用的热量传递过程，称为表面对流传热，简称对流传热。 5．辐射传热：物体不断向周围空间发出热辐射能，并被周围物体吸收。同时，物体也不断接收周围物体辐射给它的热能。这样，物体发出和接收过程的综合结果产生了物体间通过热辐射而进行的热量传递，称为表面辐射传热，简称辐射传热。 6．总传热过程：热量从温度较高的流体经过固体壁传递给另一侧温度较低流体的过程，称为总传热过程，简称传热过程。 7．对流传热系数：单位时间内单位传热面当流体温度与壁面温度差为1K是的对流传热量，单位为W／(m2·K)。对流传热系数表示对流传热能力的大小。 8．辐射传热系数：单位时间内单位传热面当流体温度与壁面温度差为1K是的辐射传热量，单位为W／(m2·K)。辐射传热系数表示辐射传热能力的大小。 9．复合传热系数：单位时间内单位传热面当流体温度与壁面温度差为1K是的复合传热量，单位为W／(m2·K)。复合传热系数表示复合传热能力的大小。 10．总传热系数：总传热过程中热量传递能力的大小。数值上表示传热温差为1K时，单位传热面积在单位时间内的传热量。 四、简答题 1．试述三种热量传递基本方式的差别，并各举1～2个实际例子说明。 （提示：从三种热量传递基本方式的定义及特点来区分这三种热传递方式） 2．请说明在传热设备中，水垢、灰垢的存在对传热过程会产生什么影响?如何防止? （提示：从传热过程各个环节的热阻的角度，分析水垢、灰垢对换热设备传热能力与壁面的影响情况）3. 试比较导热系数、对流传热系数和总传热系数的差别，它们各自的单位是什么? （提示：写出三个系数的定义并比较，单位分别为W/(m·K)，W/(m2·K)，W/(m2·K)） 4．在分析传热过程时引入热阻的概念有何好处?引入热路欧姆定律有何意义? （提示：分析热阻与温压的关系，热路图在传热过程分析中的作用。） 5．结合你的工作实践，举一个传热过程的实例，分析它是由哪些基本热量传递方式组成的。 （提示：学会分析实际传热问题，如水冷式内燃机等） 6．在空调房间内，夏季与冬季室内温度都保持在22℃左右，夏季人们可以穿短袖衬衣，而冬季则要穿毛线衣。试用传热学知识解释这一现象。 （提示：从分析不同季节时墙体的传热过程和壁温，以及人体与墙表面的热交换过程来解释这一现象（主 绪论 1．设有一大平壁，面积为A ，它的一侧为温度为1f t 的热流体，另一侧为温度为2f t 的冷流体；两侧对流换热表面传热系数分别为12h h 及；壁面温度分别为12w w t t 和；壁的材料的导热系数为λ，厚度为δ，传热过程处于稳态。写出单位面积上的传热量的计算式（不考虑辐射换热）。 第一章~第四章 导热部分 1． 厚度为δ的大平壁，无内热源，λ为常数，平壁两侧表面分别维持均匀稳定的温度1w t 和2w t 。写出这一稳态导热过程的完整数学描写。 2． 上题，若λ与温度有关，其余不变，写出这一稳态导热过程的完整数学描写。 3． 1题，若壁两侧壁面均给出第三类边界条件，即已知：X =0处，流体的温度为1f t ，对流换热表面传热系数为1h ；X=δ处，流体的温度为2f t ，对流换热表面传热系数为2h ，且1f t >2f t 。写出这一稳态导热过程的完整数学描写。 4． 对第1题求解，得出温度()t f x =的关系式，并进一步写出导热量q 的关系式。 5． 厚为δ的无限大平壁，无内热源，稳态导热时，壁内温度分布情况如图所示。说明①②③三种情况下，材料导热系数0(1)bt λλ=+中，b 何时为正、为负、为零？ 6． 用一平底壶烧开水，壶底与水接触面的温度为111oC ，通过壶底的热流密度为424002 /w m ，如在壶底结一层水垢厚3mm ，1/w m C λ=??，此时水垢层与水接触面上的温度和通过的热流密度均不变，计算： (1) 水垢层与壶底接触面上的温度； (2) 单位面积上的导热热阻。 7． 人对冷热的感觉以皮肤表面的热损失作为衡量依据。设人体脂肪层的厚度为3mm ，其内表面温度为36oC 且保持不变，冬季的某一天，气温为–15oC ，无风条件下，裸露的皮肤外表面与空气的表面传热系数为252/()w m k ?，某一风速时，表面传热系数为652/()w m k ?，人体脂肪层的导热系数为0.2/w m k ?，确定： (1) 要使无风天的感觉与某一风速、气温–15oC 时感觉一样，则无风天气温是多少？ (2) 在同样是–15oC 的气温时，无风天和某一风速时的刮风天，人皮肤单位面积上的热损失之比是多少？（按大平壁处理）。 第四章 复习题 1、 试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。 2、 试说明用热平衡法建立节点温度离散方程的基本思想。 3、 推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散方程的过程十分相似， 为什么前者得到的是精确描述，而后者解出的确实近似解。 4、 第三类边界条件边界节点的离散那方程，也可用将第三类边界条件表达式中的一阶导数 用差分公式表示来建立。试比较这样建立起来的离散方程与用热平衡建立起来的离散方程的异同与优劣。 5．对绝热边界条件的数值处理本章采用了哪些方法？试分析比较之． 6．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？ 7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解时是否因为初场的假设不合适而造成？ 8．有人对一阶导数()()()2 21,253x t t t x t i n i n i n i n ?-+-≈ ??++ 你能否判断这一表达式是否正确，为什么？ 一般性数值计算 4-1、采用计算机进行数值计算不仅是求解偏微分方程的有力工具，而且对一些复杂的经验公式及用无穷级数表示的分析解，也常用计算机来获得数值结果。试用数值方法对Bi=0.1,1,10的三种情况计算下列特征方程的根:)6,2,1( =n n μ 3,2,1,tan == n Bi n n μμ 并用计算机查明，当2 .02≥=δτ a Fo 时用式（3-19）表示的级数的第一项代替整个级数（计 算中用前六项之和来替代）可能引起的误差。 解：Bi n n =μμtan ，不同Bi 下前六个根如下表所示： Bi μ 1 μ2 μ3 μ 4 μ 5 μ 6 0.1 0.3111 3.1731 6.2991 9.4354 12.5743 15.7143 1.0 0.8603 3.4256 6.4373 9.5293 12.6453 15.7713 10 1.4289 4.3058 7.2281 10.2003 13.2142 16.2594 Fo=0.2及0.24时计算结果的对比列于下表： Fo=0.2 δ=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.94879 0.62945 0.11866 前六和的值 0.95142 0.64339 0.12248 比值 0.99724 0.97833 0.96881 Fo=0.2 0=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.99662 0.96514 0.83889 前六项和的值 0.994 0.95064 0.82925 比值 1.002 1.01525 1.01163 Fo=0.24 δ=x 1、Jacobi 迭代 在Jacobi 迭代法中任一点上未知值的更新是用上一轮迭代中所获得的各邻 点之值来计算的，即 kk k k l l n l k n k a b T a T /)(1)1()(+=∑≠=- k=1,2,...,L 1×M 1 这里带括号的上角标表示迭代轮数。所谓一轮是指把求解区域中每一节点之值都更新一次的运算环节。显然，采用Jacobi 迭代式，迭代前进的方向（又称扫描方向）并不影响迭代收敛速度。这种迭代法收敛速度很慢，一般较少采用。但对强烈的非线性问题，如果两个层次的迭代之间未知量的变化过大，容易引起非线性问题迭代的发散。在规定每一层次计算的迭代轮次数的情况下，有利于Jacobi 迭代有利于非线性问题迭代的收敛。 2、Gauss-Seidel 迭代 在这种迭代法中，每一种计算总是取邻点的最新值来进行。如果每一轮迭代按T 的下角标由小到大的方式进行，则可表示为： kk k M L k l n l kl k l l n l kl n k a b T a T a T /)(1 11 ) 1(1 1) ()(++ =∑∑?+=--≠= 此时迭代计算进行的方向（即扫描方向）会影响到收敛速度，这是与边界条件的影响传入到区域内部的快慢有关的。 3、例题： 一矩形薄板几何尺寸如图所示，薄板左侧的边界温度T L =100K ，右侧温度T R =300K ，上侧温度T T =200K ，下侧温度T B =200K ，其余各面绝热，求板上个节点的温度。要求节点数目可以变化，写出程序。 解析： ⑴列出描述问题的微分方程和定解条件。 22 220t t x y ??+=??；对于离散化的问题，其微分方程根据热平衡原理得到： 《传热学》复习题 一、判断题 1．稳态导热没有初始条件。（） 2．面积为A的平壁导热热阻是面积为1的平壁导热热阻的A倍。（） 3．复合平壁各种不同材料的导热系数相差不是很大时可以当做一维导热问题来处理（） 4．肋片应该加在换热系数较小的那一端。（） 5．当管道外径大于临界绝缘直径时，覆盖保温层才起到减少热损失的作用。（） 6．所谓集总参数法就是忽略物体的内部热阻的近视处理方法。（） 7．影响温度波衰减的主要因素有物体的热扩散系数，波动周期和深度。（） 8．普朗特准则反映了流体物性对换热的影响。（） 9. 傅里叶定律既适用于稳态导热过程，也适用于非稳态导热过程。（） 10.相同的流动和换热壁面条件下，导热系数较大的流体，对流换热系数就较小。（） 11、导热微分方程是导热普遍规律的数学描写，它对任意形状物体内部和边界都适用。( ) 12、给出了边界面上的绝热条件相当于给出了第二类边界条件。 ( ) 13、温度不高于350℃，导热系数不小于0.12w/(m.k)的材料称为保温材料。 ( ) 14、在相同的进出口温度下，逆流比顺流的传热平均温差大。 ( ) 15、接触面的粗糙度是影响接触热阻的主要因素。 ( ) 16、非稳态导热温度对时间导数的向前差分叫做隐式格式，是无条件稳定的。 ( ) 17、边界层理论中，主流区沿着垂直于流体流动的方向的速度梯度零。 ( ) 18、无限大平壁冷却时，若Bi→∞，则可以采用集总参数法。 ( ) 19、加速凝结液的排出有利于增强凝结换热。 ( ) 20、普朗特准则反映了流体物性对换热的影响。( ) 二、填空题 1．流体横向冲刷n排外径为d的管束时，定性尺寸是。 2．热扩散率（导温系数）是材料指标，大小等于。 3．一个半径为R的半球形空腔，空腔表面对外界的辐射角系数为。 4．某表面的辐射特性，除了与方向无关外，还与波长无关，表面叫做表面。 5．物体表面的发射率是ε，面积是A，则表面的辐射表面热阻是。 6．影响膜状冷凝换热的热阻主要是。 传热学计算练习题 1．某平壁燃烧炉是由一层耐火砖与一层普通砖砌成，两层的厚度均为100mm ，其导热系数分别为0.9W/（m·℃）及0.7W/（m·℃）。待操作稳定后，测得炉膛的内表面温度为700℃，外表面温度为130℃。为了减少燃烧炉的热损失，在普通砖外表面增加一层厚度为40mm 、导热系数为0.06W/（m·℃）的保温材料。操作稳定后，又测得炉内表面温度为740℃，外表面温度为90℃。设两层砖的导热系数不变，试计算加保温层后炉壁的热损失比原来的减少百分之几？(%5.68) 2．在外径为140mm 的蒸气管道外包扎保温材料，以减少热损失。蒸气管外壁温度为390℃，保温层外表面温度不大于40℃。保温材料的λ与t 的关系为λ=0.1+0.0002t （t 的单位为℃，λ的单位为W/（m·℃））。若要求每米管长的热损失Q/L 不大于450W/m ，试求保温层的厚度以及保温层中温度分布(b= 71mm)( t=－501lnr －942)。 3．有一列管式换热器，由38根φ25mm×2.5mm 的无缝钢管组成。苯在管内流动，由20℃被加热至80℃，苯的流量为8.32kg/s 。外壳中通入水蒸气进行加热。试求管壁对苯的传热系数(1272 W/（m 2·℃）)。当苯的流量提高一倍，传热系数有何变化(2215 W/（m 2·℃）)。 4．在预热器内将压强为101.3kPa 的空气从10℃加热到50℃。预热器由一束长度为1.5m ，直径为φ86×1.5mm 的错列直立钢管所组成。空气在管外垂直流过，沿流动方向共有15行(对流传热核准系数为1.02)，每行有管子20列，行间与列间管子的中心距为110mm 。空气通过管间最狭处的流速为8m/s 。管内有饱和蒸气冷凝。试求管壁对空气的平均对流传热系数(56W/（m 2·℃）)。注：(空气流过15排管束时，对流传热核准系数为1.02) 5．热空气在冷却管管外流过，α2=90W/（m 2·℃），冷却水在管内流过， α1=1000W/（m 2·℃）。冷却管外径d o =16mm ，壁厚b=1.5mm ，管壁的λ=40W/（m·℃）。试求： ①总传热系数K o ；(80.8W/（m 2·℃）) ②管外对流传热系数α2增加一倍，总传热系数有何变化？(增加了82.4%) ③管内对流传热系数α1增加一倍，总传热系数有何变化？(增加了6%) 6．有一碳钢制造的套管换热器，内管直径为φ89mm×3.5mm ，流量为2000kg/h 的苯在内管中从80℃冷却到50℃。冷却水在环隙从15℃升到35℃。苯的对流传热系数αh =230W/（m 2·K ），水的对流传热系数αc =290W/（m 2·K ）。忽略污垢热阻。试求：①冷却水消耗量；(1335 kg/h)②并流和逆流操作时所需传热面积(并流6.81 m 2，逆流5.83 m 2)；③如果逆流操作时所采用的传热面积与并流时的相同，计算冷却水出口温度与消耗量(46.6℃，846 kg/h)，假设总传热系数随温度的变化忽略不计。 7．有一台运转中的单程逆流列管式换热器，热空气在管程由120℃降至80℃，其对流传热系数α1=50W/（m 2·K ）。壳程的冷却水从15℃升至90℃，其对流传热系数α2=2000W/（m 2·K ），管壁热阻及污垢热阻皆可不计。当冷却水量增加一倍时，试求①水和空气的出口温度t'2和T'2，忽略流体物性参数随温度的变化；(t'2=61.9℃，T '2=69.9℃)②传热速率Q'比原来增加了多少？(25%) 8．为了得到热水，0.361 MPa (t s =140℃) 的水蒸气在管外凝结（如图3所示），其表面传热系数29500W/(m K) o h 。冷却水在盘管内流动，流速为0.8m/s ，黄铜管外径为18mm ，壁厚为1.5mm ， 1管道总传热系数 管道总传热系数是热油管道设计和运行管理中的重要参数。在热油管道稳态运行方案的工艺计算中，温降和压降的计算至关重要，而管道总传热系数是影响温降计算的关键因素，同时它也通过温降影响压降的计算结果。1.1 利用管道周围埋设介质热物性计算K 值管道总传热系数K 指油流与周围介质温差为1℃时，单位时间内通过管道单位传热表面所传递的热量，它表示油流至周围介质散热的强弱。当考虑结蜡 层的热阻对管道散热的影响时，根据热量平衡方程可得如下计算表达式： （1-1）1112ln 111ln 22i i n e n w i L L D D D KD D D D ααλλ-+???? ?????=+++????????∑式中：——总传热系数，W /(m 2·℃)；K ——计算直径，m ；（对于保温管路取保温层内外径的平均值，对于e D 无保温埋地管路可取沥青层外径）；——管道内直径，m ；n D ——管道最外层直径，m ；w D ——油流与管内壁放热系数，W/(m 2·℃)；1α ——管外壁与周围介质的放热系数，W/(m 2·℃)；2α ——第层相应的导热系数，W/(m·℃)；i λi ，——管道第层的内外直径，m ，其中；i D 1i D +i 1,2,3...i n =——结蜡后的管内径，m 。L D 为计算总传热系数，需分别计算内部放热系数、自管壁至管道最外径K 1α的导热热阻、管道外壁或最大外围至周围环境的放热系数。 2α（1）内部放热系数的确定1α放热强度决定于原油的物理性质及流动状态，可用与放热准数、自然1αu N 对流准数和流体物理性质准数间的数学关系式来表示[47]。r G r P 在层流状态（Re<2000），当时：500Pr 、室内一根水平放置的无限长的蒸汽管道，其保温层外径d=583 mm，外表面 实测平均温度及空气温度分别为，此时空气与管道外 表面间的自然对流换热的表面传热系数h=3.42 W /(m2 K),墙壁的温度近似取为 室内空气的温度，保温层外表面的发射率 问：(1)此管道外壁的换热必须考虑哪些热量传递方式； (2)计算每米长度管道外壁的总散热量。(12分) 解： (1)此管道外壁的换热有辐射换热和自然对流换热两种方式。 (2)把管道每米长度上的散热量记为qi 当仅考虑自然对流时，单位长度上的自然对流散热 q i,c =二d h t =二dh (j - t f ) = 3.14 0.583 3.42 (48 - 23 ) 二156 .5(W / m) 近似地取墙壁的温度为室内空气温度，于是每米长度管道外表面与室内物体及墙壁 之间的辐射为： q i厂d (T； -T；) = 3.14 0.583 5.67 10》0.9 [(48 273)4-(23 273)4] = 274.7(W /m) 总的散热量为q i = q i,c +q i,r = 156.5 +274.7 = 431.2(W/m) 2、如图所示的墙壁，其导热系数为50W/(m- K),厚度为50mm在稳态情况下的 墙壁内的一维温度分布为：t=200-2000x 2，式中t的单位为°C, x单位为m 试 求： t (1) 墙壁两侧表面的热流密度； (2) 墙壁内单位体积的内热源生成的热量 2 t =200 —2000x 解：(1)由傅立叶定律: ① dt W q ' (―4000x) = 4000二x A dx 所以墙壁两侧的热流密度： q x _. =4000 50 0.05 =10000 (1)由导热微分方程 茫?生=0得: dx 扎 3、一根直径为1mm 勺铜导线，每米的电阻为2.22 10 。导线外包有厚度为 0.5mm 导热系数为0.15W/(m ? K)的绝缘层。限定绝缘层的最高温度为 65°C,绝 缘层的外表面温度受环境影响，假设为40°C 。试确定该导线的最大允许电流为多 少？ 解：(1)以长度为L 的导线为例，导线通电后生成的热量为I 2RL ,其中的一部分 热量用于导线的升温，其热量为心务中：一部分热量通过绝热层的 导热传到大气中，其热量为：门二 1 , d In 2 L d 1 根据能量守恒定律知：l 2RL -门 述二厶E = I 2RL -门 即 E = — L dT m = I 2RL - t w1 _tw2 4 di 1 , d 2 In 2 L d 1 q v 、d 2t ——' 2 dx =-(7000)= 4000 50 二 200000 W/m 3 t w1 - t w2 。 2 q x 卫=4000.： 0 = 0 四、非线笥问题迭代式解法的收敛性 每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量的变化不太大←多数情形下非线性问题迭代式解法是可以收敛的）。 使相邻两层次间未知量变化不太大的措施： 1、欠松弛迭代 常用逐次欠弛线迭法（SLUR ）：一组临时系数下逐线迭代求解+对所得的解施以欠松弛，再用欠松弛后的解去计算新的系数，常数，以进入下一层次的迭代。 实施：常把欠松弛处理纳入迭代过程，而不是在一个层次迭代完成后再行欠松弛。 )( ) ()()1(n p p n n n p n p t a b bt a t t -∑+=+ω )()1() 1()( n p p n n n p p t a b b t b a t a ω ωω -+++∑=+ ∑+=+')1('b b bt a t a n n n p p )('))(1(',n p p p p t a b b a a ωωω-+==，用交替方向线迭代法求解这一方程，就实现了SLUR 的迭代求解。为一般化起见，上式中b t n 上没有标以迭代层次的符号（J ，GS 时不相同）。 2、采用拟非稳态法 前面已指出，稳态问题的迭代解法与非稳态问题的步进法十分相似。对于非线性稳态问题，从代数方程的一组临时系数进入到另一组临时系数亦好象非稳态问题前进了一个时间层，非稳态问题的物理特性：系数热惯性越大（↑??=τρ/v c a o p ），温度变化越慢，仿此，对稳态非线性问题，可在离散方程中加入拟非稳态项，以减小未知量托两个层次间的变化，即 由 )()1()1()()(n p o p n n n p o p p n n n n p p n t a b b bt a t a V S b a b b bt a t V S b a ++∑=+?-∑?+∑=?-∑++ o p p n n p o p n n n p a V S b a t a b b bt a t +?-∑++∑= +) ()1( 一直进行到b t t n p ,收敛，虚拟时间步τ?的大小通过计算实践确定。 3、采用Jacobi 点迭代法 中止迭代的判据（该层次迭代）除前述变化率判据外，还可以规定迭代的轮数，例如规定进行4-6次ADI 线迭代就结束该层次上的计算。此时，用收敛速度低的丁迭代也就起到了欠松弛的作用。 五、迭代法的收敛速度 1、收敛速度 对给定的代数方程组（包括是临时系数的情形），采用不同的迭代方法求解时，使一定的初始误差缩小成α倍所需要的迭代轮数K 是不相的。1<α 4.3 热工设计 4.3.1 本系统用于外墙外保温时的保温层设计厚度，应根据《河南省公共建筑节能设计标准》（DBJ41/075-2006）、《河南省居住建筑节能设计标准（寒冷地区）》（DBJ41/062-2005）、《河南省居住建筑节能设计标准（夏热冬冷地区）》（DBJ41/071-2006）规定的外墙传热系数限值，通过热工计算确定。 4.3.2 ZCK无机复合保温板用于外墙外保温时，其导热系数（λ）、蓄热系数（S）设计计算值和修正系数按下表取值。 表4.3.2 ZCK无机复合保温板λ、S、修正系数 4.3.3 热工计算示例，以采用60mm保温板为例。 示例一：200mm混凝土剪力墙外贴60mm保温板，计算如下： Ra=R内+R1+R2+R3+R4+R外=0.11+0.0215+0.1149+1.1429+0.005+0.04=1.4343 Ka=1/R=1/1.4333=0.70W/(m2.K) 其中：R内为内表面换热阻，0.11m2.K/W； R1为水泥砂浆层热阻，0.02/0.81=0.0215 m2.K/W； R2为混凝土剪力墙层热阻，0.2/1.74=0.1149 m2.K/W； R3为保温板层热阻，0.06/（0.05*1.05）=1.1429 m2.K/W； R4为抗裂砂浆层热阻，0.005/0.93=0.005 m2.K/W； R外为外表面换热阻，0.04m2.K/W； 示例二：200mm加气混凝土砌块外贴60mm保温板，计算如下： Rb=R内+R1+R2+R3+R4+R外=0.11+0.0215+0.80+1.1429+0.005+0.04=2.1194 Kb=1/R=1/2.1194=0.47W/(m2.K) 其中：R内为内表面换热阻，0.11m2.K/W； R1为水泥砂浆层热阻，0.02/0.81=0.0215 m2.K/W； R2为加气混凝土砌块层热阻，0.2/（0.20*1.25）=0.80 m2.K/W； R3为保温板层热阻，0.06/（0.05*1.05）=1.1429 m2.K/W； R4为抗裂砂浆层热阻，0.005/0.93=0.005 m2.K/W； 传热学基础试题 一、选择题 1．对于燃气加热炉：高温烟气→内炉壁→外炉壁→空气的传热过程次序为 A.复合换热、导热、对流换热 B.对流换热、复合换热、导热 C.导热、对流换热、复合换热 D.复合换热、对流换热、导热 2．温度对辐射换热的影响（ ）对对流换热的影响。 A.等于 B.大于 C.小于 D.可能大于、小于 3．对流换热系数为1000W/（m 2·K ）、温度为77℃的水流经27℃的壁面，其对流换热的热流密度为（ ） A.8×104W/m 2 B.6×104 W/m 2 C.7×104 W/m 2 D.5×104 W/m 2 4．在无内热源、物性为常数且温度只沿径向变化的一维圆筒壁(t 1 >t 2,r 1 A. 0.5 B. 0.65 C. 0.15 D. 0.35 7.准则方程式Nu=f(Gr,Pr)反映了( )的变化规律。 A.强制对流换热 B.凝结对流换热 C.自然对流换热 D.核态沸腾换热 8．当采用加肋片的方法增强传热时，将肋片加在（）会最有效。 A. 换热系数较大一侧 B. 热流体一侧 C. 换热系数较小一侧 D. 冷流体一侧 9. 某热力管道采用两种导热系数不同的保温材料进行保温，为了达到较好的保温效果，应将( )材料放在内层。 A. 导热系数较大的材料 B. 导热系数较小的材料 C. 任选一种均可 D. 不能确定 10．下列各种方法中，属于削弱传热的方法是( ) A.增加流体流速 B.管内加插入物增加流体扰动 C. 设置肋片 D.采用导热系数较小的材料使导热热阻增加 11.由炉膛火焰向水冷壁传热的主要方式是( ) A.热辐射 B.热对流 C.导热 D.都不是 12.准则方程式Nu=f(Gr,Pr)反映了( )的变化规律。 A.强制对流换热 B.凝结对流换热 C.自然对流换热 D.核态沸腾换热 13.判断管内紊流强制对流是否需要进行入口效应修正的依据是( ) A.l/d≥70 B.Re≥104 C.l/d<50 D.l/d<104 14.下列各种方法中，属于削弱传热的方法是( ) A.增加流体流度 B.设置肋片 C.管内加插入物增加流体扰动 D.采用导热系数较小的材料使导热热阻增加 15.冷热流体的温度给定，换热器热流体侧结垢会使传热壁面的温度( ) A.增加 B.减小 C.不变 D.有时增加，有时减小 1 傅立叶定律 傅立叶定律是导热理论的基础。其向量表达式为: q gradT λ=-? (2-1) 式中:q —热流密度，是向量，2 /()Kcal m h ；gradT —温度梯度，是向量，℃/m ；λ—导热系数，又称热导率， /()Kcal mh C o ； 式中的负号表示q 的方向始终与gradT 相反。 2 导热系数（thermal conductivity ）及其影响因素 导热系数λ（ /()Kcal mh C o ）是一个比例常数，在数值上等于每小时每平方米面积上，当物体内温度梯度为1℃/m 时的导热量。 导热系数是指在稳定传热条件下，1m 厚的材料，两侧表面的温差为1度（K ，°C ）,在1秒内，通过1平方米面积传递的热量，用λ表示，单位为瓦/米·度，w/m·k （W/m·K,此处的K 可用℃代替）。导热系数为温度梯度1℃/m ，单位时间通过每平方米等温面的热传导热流量。单位是：W/(m·K)。 3．热传导微分方程推导 ? 在t 时刻w 界面的温度梯度为x T ?? 在t 时刻e 界面的温度梯度为dx x T x T dx x x T x T 2 2??+??=????+?? 单位时间内六面体在x 方向流入的热流量为：dydz x T ??-λ； 单位时间内六面体在x 方向流出的热流量为： dydz dx x T x T ?? ??????+??-22λ； 单位时间内六面体在x 方向流入的净热量为：dxdydz x T 22??λ 图3-1 微分单元体各面上进出流量示意图 同理，单位时间内六面体在y 方向流入的净热量为：dxdydz y T 22??λ； 单位时间内六面体在y 方向流入的净热量为：dxdydz z T 22??λ； 单位时间内流入六面体的总热量为：dxdydz z T y T x T ??? ?????+??+??222222λ (3-1) 六面体内介质的质量为:dxdydz ρ。 单位时间六面体内热量的变化量（增加）为:Cdxdydz t T ρ?? 根据热量守恒定律: Cdxdydz t T dxdydz z T y T x T ρλ??=????????+??+??222222， C t T z T y T x T ρλ??=????????+??+??222222， t T z T y T x T C ??= ????????+??+??222222ρλ， t T z T y T x T a ??= ????????+??+??222222， C a ρλ = α称为热扩散率或热扩散系数（thermal diffusivity ），单位m 2/s. λ:导热系数，单位W/(m·K)； ρ：密度，单位kg/m 3 c ：热容，单位J/(kg·K). 思考：如果单元体内有热源：单位体积单位时间的散热量是q 方程怎么变？ 4．岩石的热扩散率（导温系数) thermal diffusion coefficient ;thermal diffusivity; thermal degradation 岩石的热扩散率也叫或热扩散系数，表示岩石在加热或冷却时各部分温度趋于一致的能力。它反映岩石的热惯性特征，是一个综合性参数。热扩散率越大的岩石，热能传播温度趋于一致的速度越大，透入的深度也越大。 在t 时刻 w 界面流体速度为U ，流体温度为T 单位时间流入微元体的流体质量为：udydz dm ρ=1 带入微元体的热量为：uTCdydz ρ e 界面流体速度为dx x u u ??+ ,流体温度为dx x T T ??+ 单位时间流出微元体的流体质量为：dydz dx x u u dm ????? ? ??+ =ρ2 带出微元体的热量为: Cdydz dx x T T dx x u u ?? ??????+?????? ??+ ρ dxdydz x T dx x u C Cdxdydz x T u TCdxdydz x u uTCdydz ????+??+??+ρρρ ρ 如果不考虑x 方向速度变化，略去高阶微量，则e 界面带出微元体的热量为：Cdxdydz x T u uTCdydz ??+ρρ 单位时间内在x 方向流入六面体的净热流量为：dxdydz x T uC ??-ρ； 同理， y 方向：dxdydz y T vC ??-ρ z 方向：dxdydz z T wC ??-ρ 导热系数、传热系数、热阻值概念及热工计算方法 导热系数λ[W/(m.k)]： 导热系数是指在稳定传热条件下，1m厚的材料，两侧表面的温差为1度(K,℃),在1小时内，通过1平方米面积传递的热量，单位为瓦/米?度(W/m?K,此处的K可用℃代替)。导热系数可通过保温材料的检测报告中获得或通过热阻计算。 传热系数K [W/(㎡?K)]： 传热系数以往称总传热系数。国家现行标准规范统一定名为传热系数。传热系数K 值，是指在稳定传热条件下，围护结构两侧空气温差为1度(K，℃)，1小时内通过1平方米面积传递的热量，单位是瓦/平方米?度(W/㎡?K，此处K可用℃代替)。传热系数可通过保温材料的检测报告中获得。 热阻值R(m.k/w)： 热阻指的是当有热量在物体上传输时，在物体两端温度差与热源的功率之间的比值。单位为开尔文每瓦特（K/W）或摄氏度每瓦特（℃/W）。 传热阻： 传热阻以往称总热阻，现统一定名为传热阻。传热阻R0是传热系数K的倒数，即R0=1/K，单位是平方米*度/瓦（㎡*K/W）围护结构的传热系数K值愈小，或传热阻R0值愈大，保温性能愈好。 （节能）热工计算： 1、围护结构热阻的计算 单层结构热阻： R=δ/λ 式中：δ—材料层厚度(m)；λ—材料导热系数[W/(m.k)] 多层结构热阻： R=R1+R2+----Rn=δ1/λ1+δ2/λ2+----+δn/λn 式中: R1、R2、---Rn—各层材料热阻(m.k/w) δ1、δ2、---δn—各层材料厚度(m) λ1、λ2、---λn—各层材料导热系数[W/(m.k)] 2、围护结构的传热阻 R0=Ri+R+Re 式中: Ri —内表面换热阻(m.k/w)(一般取0.11) Re —外表面换热阻(m.k/w)(一般取0.04) R —围护结构热阻(m.k/w) 3、围护结构传热系数计算 K=1/ R0 1、对流传热总是概括地着眼于壁面和流体主体之间的热传递，也就是将边界层的（热传导）和边界层外的（对流传热）合并考虑，并命名为给热。 2、在工程计算中，对两侧温度分别为 t1,t2 的固体，通常采用平均导热系数进行热传导计算。平均导热系数的两种表示方法是或。答案；λ = 3、图 3-2 表示固定管板式换热器的两块管板。由图可知，此换热器为或。体的走向为 管程，管程流 1 1 4 2 2 3 3 5 图 3-2 3-18 附图答案：4；2 → 4 → 1 → 5 → 3；3 → 5 → 1 → 4 → 2 4、4.黑体的表面温度从 300℃升至 600℃，其辐射能力增大到原来的（5.39）倍. 答案: 5.39 分析: 斯蒂芬-波尔兹曼定律表明黑体的辐射能力与绝对温度的 4 次方成正比, ? 600 + 273 ? 摄氏温度,即 ? ? =5.39。 ? 300 + 273 ? 5、 3-24 用 0.1Mpa 的饱和水蒸气在套管换热器中加热空气。空气走管内， 20℃升至 60℃，由则管内壁的温度约为（100℃） 6、热油和水在一套管换热器中换热，水由 20℃升至 75℃。若冷流体为最小值流体，传热效率 0.65，则油的入口温度为 (104℃)。 7、因次分析法基础是 (因次的一致性)，又称因次的和谐性。 8、粘度的物理意义是促使流体产生单位速度梯度的(剪应力) 9、如果管内流体流量增大 1 倍以后，仍处于滞流状态，则流动阻力增大到原来的(2 倍) 10、在滞流区，若总流量不变，规格相同的两根管子串联时的压降为并联时4 倍。 11、流体沿壁面流动时，在边界层内垂直于流动方向上存在着显著的（速度梯度），即使（粘度）很小，（内摩擦应力）仍然很大，不容忽视。 12、雷诺数的物理意义实际上就是与阻力有关的两个作用力的比值，即流体流动时的（惯性力）与（粘性力）之比。 13、滞流与湍流的本质区别是（滞流无径向运动，湍流有径向运动） 二、问答题：问答题： 1、工业上常使用饱和蒸汽做为加热介质而不用过热蒸汽，为什么？答：使用饱和蒸汽做为加热介质的方法在工业上已得到广泛的应用。这是因为饱和蒸汽与低于其温度的壁面接触后，冷凝为液体，释放出大量的潜在热量。虽然蒸汽凝结后生成的凝液覆盖着壁面，使后续蒸汽放出的潜热只能通过先前形成的液膜传到壁面，但因气相不存在热阻，冷凝传热的全部热阻只集中在液膜，由于冷凝给热系数很大，加上其温度恒定的特点，所以在工业上得到日益广泛的应用。如要加热介质是过热蒸汽，特别是壁温高于蒸汽相应的饱和温度时，壁面上就不会发生冷凝现象，蒸汽和壁面之间发生的只是通常的对流传热。此时，热阻将集中在靠近壁面的滞流内层中，而蒸气的导热系数又很小，故过热蒸汽的对流传热系数远小于蒸汽的冷凝给热系数，这就大大限制了过热蒸汽的工业应用。 2、下图所示的两个 U 形管压差计中,同一水平面上的两点 A、或 C、的压强是否相等？ B D P1 A P2 p 水 B C 空气 C 水银图 1－1 D 水 P1 1－1 附图 P2 A B D p h1 。 答：在图 1—1 所示的倒 U 形管压差计顶部划出一微小空气柱。空气柱静止不动，说明两侧的压强相等，设为 P。由流体静力学基本方程式： p A = p + ρ空气 gh1 + ρ水 gh1 p B = p + ρ空气 gh1 + ρ空气 gh 1 Q ρ水 > ρ空气 p C = p + ρ空气 gh1 ∴ p A> pB 即 A、B 两点压强不等。而传热系数计算方法
传热学数值计算大作业2014011673
传热学试题库含答案
传热学补充题
传热学第四版课后题答案第四章
计算传热学数值模拟
传热学习题及参考答案
传热学练习题(学生)
管道总传热系数计算18
传热学计算例题
热物理过程的数值模拟-计算传热学3.(DOC)
玻璃的传热系数计算
传热学基础试题及答案
传热流体数值计算
导热系数、传热系数、热阻值概念及热工计算方法简述实用版)
流体力学与传热学
相关主题
文本预览