传热流体数值计算

傅立叶定律是导热理论的基础。其向量表达式为:

qgradT (2-1)

2T2T2TT， 222Cxyzt2T2T2TT， a a222Cyztx2式中:q—热流密度，是向量，Kcal/(mh)；gradT—温度梯度，是向量，℃/m ；—导热系数，又称热导率，

Kcal/(mhoC)； 式中的负号表示q的方向始终与gradT相反。

2 导热系数（thermal conductivity）及其影响因素

oKcal/(mhC)）是一个比例常数，在数值上等于每小时每平方米面积上，当物体内温度梯度为1℃/m导热系数（

α称为热扩散率或热扩散系数（thermal diffusivity），单位m2/s.

λ:导热系数，单位W/(m·K)； ρ：密度，单位kg/m3 c：热容，单位J/(kg·K). 思考：如果单元体内有热源：单位体积单位时间的散热量是q 方程怎么变？

4．岩石的热扩散率（导温系数) thermal diffusion coefficient ;thermal diffusivity; thermal degradation 岩石的热扩散率也叫或热扩散系数，表示岩石在加热或冷却时各部分温度趋于一致的能力。它反映岩石的热惯性特征，是一个综合性参数。热扩散率越大的岩石，热能传播温度趋于一致的速度越大，透入的深度也越大。 热扩散系数一般是根据岩石的导热系数(ramuda)、和密度(rou)的测量数据计算得到的。 如果单元体内有热源：单位体积单位时间的散热量是q 则在 时的导热量。 导热系数是指在稳定传热条件下，1m厚的材料，两侧表面的温差为1度（K，°C）,在1秒内，通过1平方米面积传递的热量，用λ表示，单位为瓦/米·度，w/m·k（W/m·K,此处的K可用℃代替）。导热系数为温度梯度1℃/m，单位时间通过每平方米等温面的热传导热流量。单位是：W/(m·K)。 3．热传导微分方程推导 ♥ T 2T2T2TT222dxdydzCdxdydz 左边增加热源散发的热量qdxdydz， txyz方程变为： T在t时刻w界面的温度梯度为

xN 2T2T2TT222dxdydzqdxdydzCdxdydz tyzx2T2T2TqT a222txyz如果是流体除了热传导外，还有流体流动带入和带出微元体的热量 E TTT2Txdxdx 在t时刻e界面的温度梯度为xxxx2单位时间内六面体在x方向流入的热流量为：单位时间内六面体在x方向流出的热流量为：

t Tdydz； xW w s x P n e Δz b Δx Δy 在t时刻

w界面流体速度为U，流体温度为T

单位时间流入微元体的流体质量为：dm1udydz 带入微元体的热量为：uTCdydz e界面流体速度为uT2T2dxdydz；

xx2Tdxdydz 单位时间内六面体在x方向流入的净热量为：2xS uuTdxdydz dx,流体温度为Tdx 单位时间流出微元体的流体质量为：dm2uxxxB 带出微元体的热量为:

图3-1 微分单元体各面上进出流量示意图

uuTdxTdxCdydz xx2Tdxdydz； 单位时间内六面体在y方向流入的净热量为：同理，单位时间内六面体在y方向流入的净热量为：2y2T2T2T2T2dxdydz； 单位时间内流入六面体的总热量为：222dxdydz (3-1) zyzx六面体内介质的质量为:dxdydz。 单位时间六面体内热量的变化量（增加）为:根据热量守恒定律:

uTCdydzuTuTTCdxdydzuCdxdydzCdxdxdydz xxxxTCdxdydz x如果不考虑x方向速度变化，略去高阶微量，则e界面带出微元体的热量为：uTCdydzu单位时间内在x方向流入六面体的净热流量为：uCTdxdydz； xTCdxdydz t同理， y方向：vCTTdxdydz z方向：wCdxdydz yz2T2T2T2T2T2TTT222dxdydzCdxdydz， 222C，

tyztyzxx将上述三项加入到热传导方程（左边，进入微元体的热量）：

1

传热与流体流动的数值计算 传热与流体流动的数值计算

2T2T2TT222dxdydzCdxdydz tyzx2T2T2TTTTT222dxdydzuCdxdydzvCdxdydzwCdxdydzCdxdydz xyzyztx2T2T2TTTTT222uvw （能量方程） yzyztxxR0 R1 0 ΔRNR PI-1 PI PI+1 RI-1 RI PNR 2.2巷壁与风流间的对流换热 运动着的流体与所接触的固体壁面之间的热量传递过程称为对流换热，它是流体(液体或气体)由于宏观相对运动，从

某一区域迁移到温度不同的另一区域时引起热量传递的现象。固体壁面与流体之间存在温度差将产生对流换热，由于实际流体的粘性和壁面摩擦的共同影响，近壁流体分层流动，尤其与壁面直接接触的几何面上，总有一层很薄的流体粘附于表面，该层流体处于静止状态，所以热流通过表面层的传递只能依靠导热。显然，在流体发生热对流的同时，由于流体中温度分布的不均匀，也将伴随产生导热现象。因此，对流换热过程实际上是热对流和热传导的综合作用过程。 牛顿冷却公式

对流换热过程是一个受很多因素影响的复杂过程，如流体的流动状况、流体的物理性质、壁的形状和大小、表面粗糙度等。一般情况下对流换热的计算可采用牛顿冷却公式。根据对流换热定律，可以计算出从壁面某处进入通风风流的显热热流密度：

P1 RI+1 图3-1 巷道围岩内节点划分

根据简化的数学模型，可将巷道围岩划分为一系列等间距 (R)的同心圆，取垂直于长轴的巷道断面角度为，如图3-1所示。

第二章 物理现象的数学描述

控制微分方程：把控制传热、流体流动等有关过程的规律表达成数学形式；详细、完整的推导，阅读标准的教科书；数值解法

方程的形式和意义：我们这里所提到的所有方程都具有一个共同的形式；形式上的一致是构成一个通用解法的基础。

2．1 控制微分方程

2．1—1 微分方程的意义

守恒原理：各个微分方程各自代表着一定的守恒原理．每一个方程以一定的物理量作为它的因变量，方程本身则代表着那些影响该因变量的各个因素之间必定存在着的某种平衡．

这些微分方程的因变量通常具有“比”的性质，即以单位质量为基础来表示各因变量。 这种因变量的例子有：质量分量、速度(即单位质量的动量)以及比焓． 这一类微分方程的各项代表着以单位容积为基础的效果

例： 设想J表示一个典型因变量Φ的流量密度．让我们考虑如图2．1所示的尺寸为dx、dy及dz的控制容积．

qs(TwT) (3) 式中：Tw= 巷道壁面的温度； T= 巷道内风流的平均温度；

= 巷道壁面的换热系数。在围岩与风流的热交换过程中，多半是井巷低温风流流经高温岩壁，井巷壁面向风流放

2oKcal/(mhC)）称为巷壁与风流的换热系数，简称为放热系数。 热，所以矿内常把上式中的对流换热系数（

α T Tw Jx(J在x方向的分量)：进入面积为dydz的流量密度

Jx(Jx/x)dx： 离开与这个面相对曲面上的流量密度

通过该面的整个面积上流出的净流量是：(Jx用同样的方法考虑y与z方向的贡献，有：

圆形巷道（柱体）围岩与风流换热控制方程

地热通过围岩向风流的传热现象与围岩本身的热传导、巷道壁面向风流的对流换热以及壁面上的水分蒸发等因素有关。由于实际情况下围岩的散热是一个很复杂的过程，为了方便本论文的研究，对要研究的物理模型做了简化和假设：

1) 巷道为圆形、无限扩展，围岩岩石均质、各向同性； 2) 不考虑围岩壁面的热辐射作用。

根据上述假设，可得到描述考虑壁面水分蒸发时围岩与风流热质传递的数学方程，如式（3-1）：

/x)dxdydz dxdydz：所讨论区域的容积

T2T1Tta(r2rr) (r0rR;t0)T(r,t)t0T0 (r0rR)  (3-1)

T(r,t)T0 (t0)rRT(TwTa)fLv(mwma) (t0)rrr0式中：R——调热圈半径，m；其他符号的意义同前章所述。

单位容积流出的净流量JxJyJzdivJ （2.1） xyz由于我们的数值方法是通过对一个控制容积进行平衡构成的(就象我们将要在后面看到的那样)，上述divJ的表达方式对我们来说是特别有用的．

以单位容积为基础来表达一项的另一个例子是变化速率项

()/t

如果Φ是某个“比”性质（单位质量）而ρ是密度，那么ρΦ就代表在单位容积内所包含的相应广延性质的大小．于是

2

传热与流体流动的数值计算 传热与流体流动的数值计算

()/t是单位容积内有关性质的变化率．

一个微分方程是这样一些项的组合；其中每一项代表一个以单位容积为基础的效应；而所有的项合在一起则反映着某种平衡或是守恒。

我们现在以几个标准的微分方程为例，并由此找到一个通用的形式．

(u)div(uu)div(gradu)BxVx （2.11） t 其中μ是粘度，p是压力，Bx是沿着x方向的单位容积内的体积力，Vx代表除去以表示的粘性力项之外的其它所有粘性力项。 div(gradu)所2．1—2 化学组分的守恒

令ml代表一种化学组分l的质量分量．当存在有速度场u时，可把ml的守恒表示为： 2．1—5 紊流的时间平均方程 在实际应用中，常常会遇到紊流的问题．就工程中所遇到的实际问题，人们通常关心的是这种流动状态的时间平均特性． 因此人们通过平均运算的方法把不稳态的层流方程转化为对紊流流动的时间平均方程． 在进行这种平均运算的时候，人们假设：紊流中存在有相对平均值的快速而随机的脉动．由平均运算所产生的附加项是所谓的雷诺（Reynolds）应力，紊流热流密度，紊流扩散流量密度等等．紊流模型的任务就是确定用流动的平均性质来表示这些附加量的方法． 许多紊流模型采用紊流粘度或紊流扩散系数的概念来表示紊流应力及流量密度．结果，紊流的时间平均方程就具有了与层流流动方程完全相同的形式。但是，诸如粘度、扩散系数以及导热系数这样一些层流交换系数则需要用相应的有效(即层流加紊流)交换系数取代。 从计算机计算的观点看，按照这种格式描述的紊流即相当于具有一个相当复杂的粘度表达式的层流。 (ml)div(umlJl)Rlt （2.2）

(ml)/t: 单位容积内化学组分l的质量变化率； ρuml:对流流量密度(即一般化流场ρu所携带的流量密度)．

Jl:扩散流量密度，它通常是由ml的梯度引起的． 两部分流量密度(对流与扩散)的散度构成微分方程的第二项． Rl:单位容积的化学组分l的生成率.化学组分的生成系由化学反应所致。当然，依反应实际上是产生还是消毁组分l而定，Rl可能是正的，也可能是负的．对于不参与化学反应的组分，Rl为0．

如果用菲克(Fick)扩散定律来表示扩散流量密度Jl，我们可以写出：

Jllgradml (2.3) 其中，l是扩散系数．将方程(2.3)代入方程(2.2)，可以导出：

(ml)div(uml)div(gradml)Rl (2.4)

t2．1—3 能量方程

(u)div(uu)div(gradu)BxVxtt(u)(v)(w)0

xyz最通用形式表示的能量方程式含有相当数量的各种不同影响因素．因为我们主要关心的是方程的形式而不是其细节，

所以考虑某些限定的情况就足够了。

(v)(v)(v)(v)pvvvuvw(t)(t)(t)对于可忽略粘性耗散作用的稳态的低速流，能量方程可写成

(u)(u)(u)(u)uuupuvw(t)(t)(t)txyzxxyyzzxdiv(uh)div(kgradT)Sk (2.5)

其中h是比焓 k是导热系数 T是温度

txyzxxyyzzySh是容积发热率

(w)(w)(w)(w)wwwp uvw(t)(t)(t)txyzxxyyzzz2.1—6紊流动能方程 现今普遍流行的紊流“双方程模型”把紊流脉动动能k的方程做为其中的方程之一，该方程可以写作： 根据傅里叶(Fourier)热传导定律，项div(kgradT)代表流体中的热传导作用． 对理想气体以及固体与液体，我们可以应用下列关系：

cgradTgradh (2.6) 其中c是定压比热．将该式代入能量方程，就得到：

(k)div(uk)div(Γkgradk)Gt的微分方程与此相类似． （2.12） kdiv(uh)div(gradh)Sk （2.7）

c如果c是常数，则h ~ T关系可以简化为：hcT （2.8） kSkdiv(uT)div(gradT)结果导出：

cc在这种情况下，不管是焓还是温度都可以选作为因变量。 取方程中的速度u为0，就得到了稳态的势传导方程：

（2.9）

其中是k的扩散系数，G是紊流能量的生成率，以及是动能的耗散率．量G是方程中的净源项．对变量2.1—7通用微分方程 ♥ 通过上面简单地列举的一些有关微分方程，可以看出：在这里，我们所感兴超的所有的因变量似乎都服从一个通用的守恒原理．如果用Φ表示固变量，通用的微分方程就是： ()div(u)div(grad)S （2.13） t 其中是扩散系数，S是源项．对于特定意义的Φ，具有特定的量和S（实际上，我们本应采用符号div(kgradT)Sk0 （2.10）

及S2．1—4 动量方程

对于牛顿流体，控制给定方向上的动量守恒的微分方程可以写成类似线性的形式．但是由于必须同时考虑切应力和正应力，加之与流体流动有关的斯托克斯粘性定律要比菲克(Fick)定律或傅里叶定律复杂，因此动量方程就要复杂得多．用u表示x方向的速度，我们可以把相应的动量方程写成下列的形式：

的，然而这样做会导致在随后的文章中应用过多的下标）． 上述通用微分方程中的四项分别是不稳态项、对流项、扩散项以及源项．因变量可以代表各种不同的物理量，如：化学组分的质量分量、焓或温度、速度分量、紊流动能或紊流的长度尺度。与此相应，对于这些变量中的每一个都必须相对应的扩散系数以及源项S赋以适当的意义． 3

传热与流体流动的数值计算 传热与流体流动的数值计算

把div(Γgrad)作为扩散项：

一般说来，因变量Φ是三个空间坐标与时间的函数．于是：

(x,y,z,t) (2.19)

并不是所有的扩散流量密度都是受制于有关变量的梯度的。但是把div(Γgrad)作为扩散项的做法并没有把通用的Φ方程只限于那些以梯度驱动的扩散过程．

凡是不能归入名义的扩散项的因子或项总是可以表示成源项的一部分。如果需要的话，我们甚至于可以把扩散系数取为0．由于大多数的因变量确实需要有一个突出梯度驱动扩散这一特性的扩散项，所以我们把具有这一特性的项以显式的形式写在通用的Φ方程中了．

出现在方程(2．13)中的密度可以通过状态方程与温度和质量分量这样一些变量相关联．这些变量与速度分布都服从通用微分方程．此外，流场还应当满足一个附加的约束条件，即质量守恒或连续性方程．这个方程是：

其中x，y，z以及t都是自变量．

并不是所有的问题都需要考虑所有四个自变量．所涉及的自变量数愈少，需要计算Φ值的位置〔或网格结点)也愈少． 一维问题：当有关的物理量只与一个空间坐标有关； 二维问题：所考虑的问题与二个空间坐标有关； 三维问题：所考虑的问题与三个空间坐标有关； 稳态问题：当所考虑的问题与时间无关时； 不稳态问题：当所考虑的问题与时间有关时。 把空间与时间两个因素放在一起考虑，我们将把一种状态说成是不稳态的一维问题，稳态的三维流动等等．

等温面--等温迁移法♥

象方程(2．19)所表示的那种独立坐标的选择方式并不是唯一的形式．代替把稳态温度分布写成T(x，y，z)，我们可以用另外一种写法：

div(u)0 （2.14） t上面，我们已经用向量的形式写出了有关的通用微分方程(2．13)和连续性方程(2．14)．这些方程的另一种有用的表达方式是直角坐标的张量形式：

zz(x,y,T) (2.20)

()(uj)()S （2.15） (uj)0 （2.16） txjxjxjtxj这里下标j可以取信1、2、3分别代表三个空间坐标．如果在一项内下标j重复出现，就意味着要取三项之和，例如：

(uj)(u1)(u2)(u3) （2.17） xjx1x2x3(Γ)(Γ)(Γ)(Γ) （2.18） xjxjx1x1x2x2x3x3用直角坐标张量形式表达的一个直接好处是，只要简单地把下标j抹掉，就可以由这种形式得到该方程的一维形式．

把任何特定的微分方程改写成通用形式(2．13)的过程，就是把有关因变量的不稳态项、对流项及扩散项转换成共同的标准形式．于是，把扩散项内梯度Φ的系数取为对的表达式；而把方程右端的其余各项之和定义为源项S． 尽管到现在为止我们一直都是把所有的变量当成有因次的量来考虑，但是，以无因次的变量来进行研究则往往比较方便。可以认为任何一个用无因次变量表示的微分方程都具有通用的形式2．13，

这里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的圆柱极坐标系内则是二维的，这是因为：

(r,z) (2.21) 与θ无关。

3．坐标变换可能用来进一步减少自变量的数量，例如：

a．在一块平板上的二维层流边界层给出了速度仅与η有关的相似性状态，其中： cy (2.22) x()div(u)div(Γgrad)S (2.13) (重点) t这时Φ代表无因次的因变量，而和S分别代表扩散系数和源项的无因次形式．在许多情况下，的无因次值可以简化为1，而S可以取值0或1． 热、质传递，流体流动，紊流以及有关的一些现象的所有各有关微分方程都可以看成是通用的Φ方程的一个特殊情况． 我们需要关心的仅仅是方程(2．13)的数值解．在编制计算机程序时，只需要写出一个求解方程(2.13)的通用程序就足够了．们可以对不同意义的Φ，重复使用这个程序；当然需要对相应的和S分别赋以各自合适的表达式，同时也需要给出相应合适的初始条件与边界条件．

通用的Φ方程的概念使我们能够列出一个通用数值方法的公式，并编制通用的计算机程序．

式中c是一个有因次的常数．结果二维的问题便简化为一维的问题了．

b．在半无限固体内的不稳态导热问题具有x和t这样两个自变量．但是，对于某些简单的边界条件，可以把温度描写成只与ξ有关，其中

Cx (2.23) 式中C代表一个适当的有因次常数 t4．改变因变量可能导致自变量数目的减少．例如：

a．在一充分发展的通道流中，温度T与流动方向的坐标x以及横向坐标y有关．但是在具有均匀壁温我们有： (y) (2.24) 式中：Tw的热发展区，

TTw Tb是整体温度，它随x而变化。

TbTwyu,  (2.26)

uc~u~() (2.25) 式中：u~b．平面自由射流是一种二维流，但是我们可以写成：u2．2坐标的性质

到现在为止，我们都在注意因变量．现在我们将回过来考虑自变量（x,y,z,t），并从计算的观点讨论一下它们的性质．

2．2—1 自变量 这里uc代表中心线上的流速，y是横向坐标，δ是射流的特征宽度．uc和δ这两个变量都随流动方向的坐标x而变化． 为了提高计算效率，应当选择合适的坐标系统进行数值计算．

4

传热与流体流动的数值计算 传热与流体流动的数值计算 2.2—3 单向与双向坐标 (重点)♥

【定义】双向的坐标：如果在一个坐标上的一个给定位置处的条件，只受该位置两侧条件变化的影响的话，那么这个坐标就是一个双向的坐标． 单向的坐标：如果在一个坐标上的一个给定位置处的条件，只受该位置一侧条件变化的影响，这样的坐标就是一个单向的坐标．

【例】通常，空间坐标是双向坐标。在一根园棒内的一维稳态导热是—个双向坐标的例子．在棒内任意给定点的温度会受到两端温度变化的影响．

时间坐标则总是单向坐标。在一块固体的不稳态冷却期间，其某一给定瞬时的温度只能受在该瞬时以前所发生的那些条件变化的影响．这是一个常识问题：昨天的事情影响着今天发生的事情，但是明天的条件对于今天会发生什么事情则没有影响．

空间也可以作为单向坐标：在流体流动的作用下，甚至某种空间坐标也能变成非常接近于单向的坐标．如果在一个坐标方向上有很强的单向流动，那么各种重要的影响只能是从上游传播到下游．于是在某一点上的状态主要受其上游条件的影响，而受其下游条件的影响很小．空间坐标的单向特性是一种近似．确实，对流是一种单向的过程，而扩散(这个过程总是存在的)则具有双向的效应．但是，当流量很大时，对流作用远远超过扩散作用，因而空间坐标就近乎是单向的了。

抛物型、椭圆型与双曲型——单向与双向坐标

看起来，通常用来进行微分方程分类的抛物型与椭圆型这两个数学术语相当于我们这里所用的计算方法上的概念——单向与双向坐标．抛物型这一术语表示一种单向的状态，而椭圆型则表示双向的概念。

那么前面提到过的第三种类型，即双曲型，又是什么意思呢？碰巧双曲状态不能准确归入计算方法上的分类．双曲型问题具有一种单向的特性，但是这种单向的特性并不是沿着坐标的方向，而是沿着一些称之为特征线的特殊线的方向。有一些采用特征线的数值方法，但是这些方法只限于双曲问题。另一方面，本书所提出的数

值方法并没有利用双曲问题所具有的这一特殊性质．我们将把双曲问题当成是椭圆问题这一总类中的一个组成部分(即全部是双向坐标)．

计算方法上的含义 上述关于单向和双向坐标的讨论，其动机在于：如果可以用一个单向的坐标来规定一个给定的状态，那么就有可能大大节省计算机的存储量与时间．现在让我们来讨论一个不稳态的二维热传导问题．我们将在计算域内构成一个二维网格结点的阵列。在任一瞬间，将存在一个相应的温度场．在计算机内将对时间上顺序推移的每一瞬时计算这样一个温度场，但是，由于时间是一个单向坐标，因而某一时刻的温度场不受未来时刻温度场的影响．实际上，整个不稳态问题

可以简化为按要求重复进行的一个基本的步骤，即给定t时刻的温度场，求得t＋Δt时刻的温度场．于是计算机内存只要供这两个温度场使用即可；对于所有的各个不同时间步可以一遍又一遍地使用同样的储存空间．

在这样的情况下，从一给定的初始温度场开始，我们就能够“前进”到在时间上顺序推移的各个瞬时．在任何时间步，需要同时处理的未知量只有一个二维的温度数组．这个温度数组与所有未来的温度值无关，而影响这一温度数组的以前瞬时的值则是已知的．于是我们需要求解的就是一个简单得多的方程组，从而可以节省计算机的时间．

3.1数值方法的本质

3.1—1 任务

数值解系与实验数据：

一个微分方程的数值解系由一组构成因变量Φ的分布的数所组成。在这个意义上讲，数值方法有点类似于在实验室中进行的实验。做实验时，仪器的一组读数为我们构成在所研究的域内被测量

的物理量的分布，数值分析与实验室实验两者所获得的数据都只能是一些有限数量的数值。 未知量：各个不同位置上的Φ值

让我们假设用一个关于x的多项式来代表Φ的变化

a0a1xa2x2a3x3amxm （3.1）

并采用数值方法来求得有限数量的系数a0,a1,a2,am,把x的值以及各个a的值代入方程(3．1)，我们就可以

离散化方法

主要任务：我们所感兴趣的现象是受一些微分方程支配的。我们已经把这些微分方程用一个变量为Φ的通用方程全部归纳了．现在我们的主要任务就是推导求解这个通用方程的方法．

()div(u)div(Γgrad)S (2.13)

t假设：为了便于理解，在本章中我们将假设变量Φ仅仅是一个自变量x的函数

计算出在任何位置x的Φ值．

但是如果我们最终的兴趣是得到在各个不同位置上的Φ值，那么这种做法就有点不太方便了。各个a值本身是没有什么特别的意义的．要知道所需要的Φ值，还必须进行前面所述的代入计算过程．

想法：建立一个把一系列给定点上的Φ值作为原始的未知量的方法.实际上求解微分方程的大多数方法都属于这一类，因此我们将把我们的注意力放在这样一类方法上．因而数值方法就是把计算域内有限数量位置(叫做网格结点)上的因变量值当作为基本的未知量来处理．

该方法的任务是

提供一组关于这些未知量的代数方程 并规定求解这组方程的算法． 3.1—2 离散化的概念

离散化方法：把注意力集中在网格结点处的值，用离散的值取代包含在微分方程精确解中的连续信息．于是我们已经离散了Φ的分布，把这一类数值方法叫做离散化方法． (x)w(x)e 离散化方程：关于所选取的网格结点上未知Φ值的代数方程(我们现在就把它们叫做离散化方程)系由支配Φ的微分方程推导而得。 weWPE网格结点之间Φ的分布： x推导过程中，我们必须对网格结点之间Φ如何变化作某种假设．尽x 管可以选择在整个计算域内满足一个简单表达式的关系作为Φ的这种 图3-2 一维问题的网格节点群 “分布”。 采用分段分布是更为实际的方法，即：一定的段仅仅用一个小区域的内部及其边界上的网格结点上的Φ值来描述该区间内Φ的变化；于是一般将计算域分成一定数量的子域或单元．每一个子域可以有一个独立的分布假设．

离散化的概念：这样，连续的计算域已经被离散开了．这种对空间和因变量所作的系统的离散化使得我们有可能用比较容易求解的简单的代数方程取代前面提到过的控制微分方程。

一个离散化方程是连接一组网格结点处Φ值的代数关系式;这样的一个方程系由支配Φ的微分方程推导而得;表示与该微分方程相同的物理信息. 3.1—3 离散化方程的结构

Φ的分布 一定的离散化方程只与少数的几个网格结点有关，这种情况是我们选取分段分布这一特性的结果．因此，在一个网格结点处的Φ值只影响与其紧相邻的一些点上的Φ分布．可以预料到，当网格结点的数目变得很大时，离散化方程的解将趋近于相应微分方程的精确解。这是出自于这样的考虑：当网格结点紧挨在一起时，在相邻点之间的Φ变化就变得很小，因此有关分布假设的实际细节就不那么重要了．

离散化方程不是唯一的：对于一个已知的微分方程，可能的离散化方程决不是唯一的；尽管在网格结点数非常大的极限条件下，预计所有这些可能类型的离散化方程将会给出相同的解．离散化方程的不同形式起因于分布假设以及推导方法的不同．

有限差分和有限元法：可以认为这两种方法是离散化方法的两种可供选择的型式．有限差分法与有限元法之间的区别来自选择分布和推导离散化方程的方法不同．在本书内集中注意的方法具有有限差分法的外形，但是它采用了许多属于典型的有限元方法论所具有的思想．

5

传热与流体流动的数值计算 传热与流体流动的数值计算 3.2推导离散化方程的方法 对于一个已知的微分方程，可以用许多方法推导出所要求的离散化方法。这里我们将简单地介绍几种普通的方法，然后再指出我们所喜爱的一个方法． 3.2—1泰勒级数公式

通过截断泰勒级数来近似表示微分方程的导数构成。对于图3．1中位于结点1与结点3之间中点的结点2 (xddT(k)S0 (3.10) dxdxedTdT(k)e(k)wSdx0 (3.11)

wdxdxx2x1x3x2)，在2周围展开的泰勒级数给出： 分布曲线的假设 需要一个分布假设或是一个内插分式．图3.3中表示了两种简单的分布假设．

阶梯式分布：假设在一个网格结点处的T值代表它周围整个控制容积内的T值（阶梯式分布），对于这样的分布，

斜率

dT/dx在控制容积面(即在w或e)上是不确定的

dT/dx，

分段线性的分布(图3.3b)．这时，在网格结点之间采用线性的内插函数．

离散化方程 如果我们用分段线性分布来计算方程(3.11)中的

在第三项之后截断级数，将二个方程相加及相减： edTdT(k)e(k)wSdx0 (3.11)

wdxdx 把这二个表达式代入微分方程就推得有限差分方程． 这个方法含有这样的假设：Φ的变化多少有点象是x的一个多项式，从而高阶导数是不那么重要的。但是当存在(比方说)指数形式的变化时，这种假设就可能导致人们不希望要的那种公式． 泰勒级数公式的推导是比较直截了当的，但是缺乏弹性并且其中各项的物理意义难以理解。

所得的方程将为：

ke(TETP)kw(TPTW)Sx0 (3.12)

(x)e(x)w其中S为S在整个控制容积内的平均值．离散化方程可缩写成下列形式：

aPTPaETEaWTWb (3.13)

kwkea其中： aE W aPaEaW bSx

(x)w(x)e说明

1. 离散化方程的标准形式: 方程(3．13)

方程的左边: 在中心网格结点上的温度Tp；右侧:由相邻结点上的温度和常数b构成，对二维与三维的情况相邻

结点的数目增加． 一般说来，比较方便的是把方程(3.13)看成具有如下的形式：

3.2—2变分公式 3.2—3加权余数法 3.2—4控制容积公式 (略 见课本P31-34) 3.3 一个说明性的例子 让我们来讨论受下列方程控制的一维稳态热传导问题： aPTPanbTnbb （3.15）

式中下标nb表示一个相邻结点，∑表示对所有的相邻结点求和．

2. 分布假设: 在我们推导方程(3．13)时，我们已经采用了使我们能够估计dT/dx值的最简单的分布假设，当然选用许多其它形式的内插函数本来也是可以的．

3. 此外，我们没有必要对所有的量都采用同样的分布函数，明白这一点是重要的．例如，既没有必要用网格结点之间线性变化的S来计算S，也没有必要由kp和ke之间线性变化的k来计算ke。

4. 即便对于一个确定的变量，也没有必要对方程中所有各项都采用同样的分布函数假设．例如，倘若在方程(3．10)中含有一个单独包含T的附加项，或许对该项使用阶梯分布函数也是可以允许的，而不必坚持如在计算dT/dx时所采

()ddTdiv(u)div(Γgrad)S (k)S0 （3.10） dxdxt其中k是导热系数，T是温度，以及S是单位容积的发热率． 准备 为了推导离散化方程，使用图3.2中所示的网格结点群 (x)w(x)e 集中注意网格结点P

 以网格结点E及W作为它的两个邻点（E表示东侧，即正的x方向，而W表示西侧，或是负的x方向）． weWP 虚线表示控制容积面；就目前来说，它们的准确位置是无关紧要的．字母e与w代表这些面． x 对于所考虑的一维问题，我们将假设在y与z方向为单位厚度．于是，图中所示的控制容积是Δx×1×1。如果我们在整个 图3-2 一维问题的网格节点群 控制容积内积分方程，我们就得到： Ex用的分段线性分布．

指导原则 上述有关选择分布函数的自由度，最终导致不同变型的离散化方程形式．事实是，当网格结点的数目增加时，可以预料所有这些不同形式的方程都会给出相同的解。 要求：即便是采用很粗的网格，解也总应该满足：（1）物理上的真实性；（2）总的平衡。但是，我们将在此提出一个附加的要求，这个要求的加入将使我们能够大大减少我们可以接受的公式(方程)的数目．

物理上的真实性是容易理解的．图3．4所示的变化说明了这一概念．一个真实的变化应当具有与准确变化相同的定

6

传热与流体流动的数值计算 传热与流体流动的数值计算 性倾向．在无内热源的热传导问题中，内部没有一处的温度可以超出边界温度所确定的温度范围之外．当一块热的固体为绕流的流体所冷却时；固体的温度不可能降低到比该流体的温度还低．我们将总是用这样的试验来检验我们的离散化方程．

总平衡的要求意味着对整个计算域应当满足积分守恒．我们将坚持要求热流密度、质量流量以及动量通量必须准确地同相应的源和汇建立平衡，这种平衡不应当只是限于网格结点的数目变得很大时的情况，而是对于任何数目的网格结点都应该得到满足。（1）控制容积公式会使这一总的平衡成为可能，但是如同我们很快就会知道的那样，（2）在计算控制容积界面的热流密度、质量流量以及动量通量时还需要小心处置．

物理上的真实性以及总的平衡这两个约束条件将用来指导我们选择分布假设以及所采用的有关措施．在这些约束条件的基础上，我们将建立起几条基本的法则，应用这些法则我们就可以对现有的一些公式进行鉴别，并进而发展出一些新的公式来．通常需要由数学上的考虑才能作出的一些决断，现在可以直接由物理学上的原理来进行指导了。

源项的处理 一般来说，源项是因变量T本身的函数，因而在构成离散化方程的过程中需要知道这种函数关系。由于离散化方程需要用线性代数的技术来求解，因而我们在形式上只能考虑一种线性的函数关系。在下一章中将讨论对一个已知S~T关系进行“线性化”处理的方法。这里把平均值S表示成下列形式是足够的。

定的控制容积的．

法则2：正系数 大多数的实际问题应该是这样的，即某个网格结点处的因变量值只是通过对流以及扩散的过程才受到相邻网格结点上的值的影响．在其它条件不变的情况下，在一个网格结点处该因变量值的增加应当导致相邻网格结点上该值的增加（而不是减少）

在方程（3.13)中，如果TE的增加必然导致TP增加的话，这就必然要求系数aE与aP具有相同的符号。

挽句话说，对于通用性方程(3．15)，中心结点的系数aP与各相邻结点的系数anb全都必须具有相同的符号．当然，让我们可以把它们全部取为正值或者全部取为负值。让我们决定这样来写我们的离散化方程，使所有的系数均为正值，这样，法则2可以表述如下：

所有的系数(aP以及各相邻给点系数anb)必须总是正的．

aPTPaETEaWTWb (3.13) aPTPanbTnbb （3.15）

说明 由方程(3．14)所给出的系数的定义表明，我们所用的离散化方程[方程(3．13)]遵守正系数法则．但是，正如我们将在后面看到的那样，有许多公式经常违反这一法则．结果往往是得到一个物理上不真实的解．存在负相邻结点系数就会导致这样的情况，在这种情况下，一个边界温度的增加会引起相邻网格结点上的温度的降低．我们将只接受那些确保在所有情况下系数均为正的公式．

SScSPTP (3.16)

式中，Sc代表S的常数部分，而Sp是Tp的系数（显然Sp不代表在结点P所计算的S）．

在方程(3．16)中Tp的出现表明，在表示平均值S时，我们已经假设：Tp值代表整个控制容积内的值，换句话说，已经采用了图3.3中中所示的阶梯式分布(应当注意：在对dT/dx项采用分段线性分布的同时，我们可以自由地对源项采用阶梯式分布)。

应用线性化的源项表达式，离散化方程的样子看起来仍然象方程(3.13)，但是系数的定义[方程(3.24)]要有所改变．新的方程组是：

aEke(x)e

aWkw (x)waPaEaW bSx (3.14)

法则3：源项的负斜率线性化 要是我们研究一下方程(3.18)中系数的定义，就可以发现，即便所有相邻结点的系数均为正，由于SP项的关系，中心结点的系数aP仍可能变为负值．当然，只要SP不为正值，这一危险就完全可以避免．

aEke(x)e

aWkw (x)waPaEaWSPx bScx (3.18)

aPTPaETEaWTWb (3.17)

其中：

于是我们把法则3写成：当源项线性化为

SScSPTP，系数SP必须总是小于或是等于0．

aEke(x)e

aWkw (x)waPaEaWSPx bScx

上述导言性的讨论为我们提供了足够的预备知识，这样，我们就可以来推导我们的离散化方程所应当服从的一些基本法则公式，以确保所得到的解满足物理上真实以及总的平衡这两个要求．这些看来简单的法则具有深远的含义，它们将指导我们推演贯穿本书的那些方法．

3.4 四项基本法则

法则1：在控制容积面上的连续性 当一个面作为两个相邻控制容积的公共面时，在这两个控制容积的离散化方程内必须用相同的表达式来表示通过该面的热流密度、质量流量以及动量通量．

讨论：显然，通过一个特定的面离开一个控制容积的热流通过同一个面进入第二个控制容积的热流密度相同，不然的话，总的平衡就不会满足．尽管这一要求易于理解，但是稍不小心，就可能在一些微妙的地方违反。对于图3.2所示的控制容积，我们或许会用

通过TW、TP及TE的二次曲线来计算界面上的热流密度k注意 这一法则看起来好像有点任意，其实并不然．大多数物理过程确实在源项与因变量之间具有负的斜率关系，实际上，要是SP为正的话，物理状态就可能会变得不稳定了．一个正的SP意味着，当TP增加时，源项也随着增加；如果这时没有有效的散热机构；这可能会反过来导致TP的增加，如此反复进行下去，造成温度飞升的不稳定现象．从计算方法上讲，保持负的SP，使之不致产生不稳定性以及物理上的不真实解，这是至关重要的．在下一章中，我们还要进一步讨论源项线性化的问题．这里要充分注意到：为了使计算成功，负SP的原则是必不可少的．

法则4：相邻结点系数之和 控制微分方程往往只包含有因变量的导数项．

ddT(k)S0 （3.10） dxdx于是，如果T代表因变量，则函数T与T+c(其中c是一个任意常数)两者均满足微分方程．微分方程所具有的这一特性也必定要反映在与之相对应的离散化方程中．因此，当TP以及所有的Tnb都增加同一常数值时，方程(3．15)应当仍然适合．

aPTPanbTnbb （3.15）

aP(TPC)anb(TnbC)b aPCanbC

由这个要求可以得出结论：aP必须等于所有相邻结点的系数之和．因此，我们可以把法则4表述如下： 为了使微分方程在因变量增加一个常数之后也仍然能得到满足，我们要求：aP讨论 很容易理解，方程(3.13，当源项与T无关时)确实满足这一法则．

dT。对下一个控制容积采用同一类的公式意味着：公dx共界面上的梯度dT/dx系由不同的、与正要考察的控制容积有关的分布曲线算得．所得到的dT/dx (因此热流密度)的不连续性绘于图3．5中．

另一种做法也可能导致热流密度的不连续性．这时假定在给定的控制容积的各个表面上，热流密度完全为控制容积中心结点的导热系数kp所控制．这样在考虑P点周围的控制容积时，在界面e处的热流密度(见图3．2)将表示成

anb (3.19)

aPTPaETEaWTWb (3.13)

本法则意味着：中心结点值TP是各相邻结点值Tnb的一个加权平均值．与方程(3.13)不同，方程(3.17，当源项与T有关时)中的系数不服从这一法则．

kp(TPTE)/( x)e。而在把E作为控制容积的中心结点时，热流密度应为ke(TPTE)/( x)e．为了避

免出现这种不连续性，记住这一点是有好处的，即：必须把通过界面上的热流密度看成是属于界面本身，而不是属于一

7

传热与流体流动的数值计算 传热与流体流动的数值计算

aPTPaETEaWTWb (3.17)

当源项与T有关时:不能认为这种情况是对法则4的违背，而应该是本法则对这种情况不适用．当源项与T有关时，T与T＋c两者不能同时满足微分方程．就是在这样的情况下，也不应当把这个法则忘记掉，而应该通过设想方程(3.17)的一个特殊情况来应用这个法则．例如，如果我们取方程(3.17)中的SP为0，那么本法则就变得可以应用了，并且确实是被遵守的。

温度场的确定性：当T以及T+c两者均满足微分方程时，所要求的温度场并不会变成是多值或不确定的．T的值可以由适当的边界条件唯一地确定．遵照法则4就可以确保：例如，倘若边界温度增加一个常数值，那么所有的温度就会准确地增加同一常数值．

理解法则4作用的另一方法是：若是源项不存在，而且所有的相邻结点的温度Tnb都相等，则中心结点的温度

TP必定等于它们．在这样的条件下，只有一个蹩脚的离散化方程才会算出TP≠Tnb的结果来．

第四章 热传导

4．1本章的对象 (略 见课本P46) 4．2 一维稳态热传导 4．2—1 基本方程

作为解释四项基本法则的工具，在3．3节中介绍说明性例子的过程中，我们已经推导了稳态一维导热问题的离散化方程．我们来复习一下那一部分的主要内容．对这个问题的控制微分方程是：

由这个方程可以推导出离散化方程

如何设计一个合适的非均匀网格？？

因为在问题解决之前，T~x的分布是不知道的，我们能够怎样来设计一个合适的非均匀网格呢？ （1）首先，人们通常对于所要得到的解有着某些定性的预计，他们可以从这些预计得到某些指导． （2）其次，可以用初步的粗网格的解来求得T~x 变化的形式，而后可以构成一个合适的非均匀网格．

这就是为什么我们坚持要求我们的方法，即便是在粗网格的条件下也应该给出具有物理意义的解的原因之一。如果这种方法只是对足够细的网格才能给出合理的解，那么一个试探性的由粗网格求得的解将是无用的． 达到给定精确度所需要的网格点数及分布？？

达到给定精确度所需要的网格点数以及这些网格点在计算域内应取的分布方式与所要求解的问题的特性有关．采用仅仅几个网格点进行试探性计算，为弄清有关解的情况提供了一个方便的途径．最后，这恰恰是人们在实验室中相当普遍地采用的一种做法．在那里，人们先进行预备性的实验或是探索性的试验，然后再应用由这些实验(或试验)所得到的数据资料来确定在最终的实验中所应安装的探头的位置与数目。 4．2—3 界面导热系数 界面导热系数：在方程(4．3)中，已经把导热系数ke用来代表通过控制容积面e的k值，类似地，kw代表界面w的k值．当导热系数k是x的函数时，我们往往只知道在网格点W、P、E上的k值．于是我们需要有一个如何由这些网格点值来计算界面导热系数(如ke)的规定。当然，下列讨论并不是针对均匀导热系数的情况而言的．

不均匀导热系数：不均匀导热系数的状况可以由材料的不均匀性(如组合的材料扳)引起．即便对于均匀的材料，由于温度分布的不均匀，加上导热系数随温度而变化的函数关系，也会导致导热系数发生变化．在处理对Φ的通用微分方程时，扩散系数Γ将起着与导热系数相同的作用．在实际问题中，经常会遇到Γ发生明显变化的情况．例如，在紊流的情况下，Γ可以代表紊流粘度，也可以表示紊流导热系数．因此人们迫切希望能有—个适合于不均匀的k或Γ的公式。

求得界面导热系数ke的最直截了当的方法是假设k在P点和E点之间呈线性变化。于是，

kefekP(1fe)kE (4.5)

其中插入因子

ddT(k)S0 dxdxfe，是用图4.l所示的距离所定义的一个比值：

aPTPaETEaWTWb (4.2)

(x)efe(x)e (4.6)

其中：

aEke(x)e

aWkw aPaEaWSPx bScx (4.3) (x)w导热系数的突然变化的情况：我们马上就要说明：在某些情况下这种简单化的做法会导致相当不准确的结果；此外，这种做法不可能精确地处理在组合材料中可能遇到的导热系数的突然变化．在阐明这种替代的办法时，我们认为我们主要关心的不是导热系数在界面e上的局部值．我们的主要目标是要得到一个通过下式描述的界面热流密度qe的良好的表达式：qe网格点P、E及W示于图3．2中，图中还标明了各个不同的距离．控制容积面位于网格点P及其相应的相邻点之间．可以认为这些面的准确位置是任意的．它们的位置可以有许多不同的安排方案。其中的一些将在4．6—1节中讨论．就目前来说，我们将简单地认为e和w相对于网格点P、E和W的位置是已知的．

量Sc及SP由源项的线性化规定，其形式是：Ske(TPTE) (4.7)

(x)eScSPTP (4.4)

至于分布假设，梯度dT/dx已经由T对x的分段线性的变化算得，而对于线性化的源项，假设TP值代表整个控制容积内的值．当然，应当记住：只要不违背四项基本法则，选择其它形式的分布曲线是可能的，也是允许的．这里，选择分布曲线的基本原则是：在四项基本法则的约束范围内宁可采纳简单一些的分布曲线，而只是在那些必要的场合，才采用高级而复杂形式的分布曲线．

一维热传导问题的许多重要方面仍然有待于讨论，我们该转到这些论题上来了． 4．2—2 网格间距 非均匀网格

对于在图3．2中所示的网格点，距离( x)e与( x)w没有必要相等．事实上，人们常常希望采用不均匀的网格间距．这样做可以使我们能够有效地扩大我们计算的功能．一般说来，只有在网格相当细的时候，我们才可能得到精确的解．但是在因变量随x变化相当慢的区域内没有必要采用细的网格．相反，在T~x变化比较陡的区域内则需要细的网格． 不均匀网格的准确度要比均匀网格的差？？

现在在一些人们中间似乎流行着一种错误的概念，这些人以为不均匀网格的准确度要比均匀网格的差些．这样的论断是没有合理依据的．网格间距应当直接与因变量在计算域内的变化方法联系起来．此外，还没有一个通用的法则可以说明相邻网格间距之间最大(或最小)比值究竟应该保持多大．

事实上，在推导离散化方程(4-2)时已经采用了这一表达式，我们想要采用的ke的表达式应当是一个可以由它表达“正确的”qe的式子．

让我们来讨论这样一种情况，围绕着网格点P的控制容积由具有均匀导热系数kp的材料所填满，而围绕着E点的控制容积则由导热系数为kE的材料所填满．对于在P点与E点之间的组合板，根据稳态无内热源一维导热的分析，可以得

TPTE出下列的结果：qe(x)e/kP(x)e/kE (4.8)

1fefe1k() (4.9)

把方程(4．6)一(4．8)合并在一起，就得到：ekPkE当界面l位于P和E之间的中点时，我们有fe＝0.5；于是：

ke0.5(kPkE) (4.10) 或

1112kPkEkekPkE (4.10b)

方程(4．10)说明ke是kP和kE的调和平均值，而不是如方程(4．5)在fe＝0.5时给出的算术平均值．

8

传热与流体流动的数值计算 传热与流体流动的数值计算

将方程(4．9)应用于系数的定义[式(4．3)]中，可以得到aE的下列表达式：

这种最终的不变的状态称之为迭代的收敛。这个收敛的解实际上是非线性方程的正确解，尽管这个解是由线性方程的求解方法得到的．

但是有可能一次次的迭代永远也不会收敛到一个解．T的值可能稳定地漂移或是以一个不断增大的振幅振荡．这一与收敛相对立的过程称之为发散．一个好的数值方法应当使发散的可能性减为最小；如我们在后面将要看到的那样，只要遵守我们的四项基本法则，就会促进收敛；此外，我们还将讨论避免发散的其它一些对策．在这一点上，要充分注意我们的方法不止是限于线性的问题．至少在原则上，任何的非线性问题都可以用刚才概述的选代技术处理． 4．2—5 源项的线性化

当源项S与T有关时，我们采用方程(4．4)所给出的线性形式来表达这一关系．

(x)e(x)eaEkkEP1 (4.11)

对aW可以写出一个类似的表达式．很清楚，aE代表点P和E之间的材料的热导． 这一公式的效能可以由下列两种极限情况很快看出： 1．

1fefe1) (4.9) 可得：ke0 (4.12) 令kE→0，则由方程ke(kPkESScSPTP (4.4)

这样做的原因是：(1)我们名义上的线性结构框架只允许采用一种形式上的线性关系，(2)线性关系的组合要比把S处

理成常数为好．

当S是T的一个非线性函数时，我们必须把它线性化，即规定Sc和Sp的值，而这两个值本身也可能是T的函数。因而在每一次迭代循环，Sc与Sp可能都要根据新的T值重新算出。S的线性化应当是S~T关系的一个良好的表达式。此外，还必须满足非正的

这意味着在一个绝热层的表面上热流密度为零，正象它理应存在的情况那样．相反，算术平均公式在这种情况下均会给出非零的热流密度．

2．

令kP>>

Sp这一基本法则．

T1fefe1kkE，那么：ke() (4.9) keEkPkEfe (4.13)

有许多方法可以用来把给定的S的表达式分解成Sc与SpTp；其中某些方法可以用下面的例子来说明．出现在这些例

子中的数值并不具有特别的意义；符号p用以表示T的估计值或是前一次的迭代值．

[例1] 已知：S＝5—4T．某些可能的线性化是： 1．

这个结果具有两方面的含义：其中的一个含义是显而易见的，而另外的一个含义则比较难于理解．方程(4．13)表明界面的导热系数k，完全与kp无关．这个结果是可以预料到的，因为围绕着P点的材料导热系数高，其热阻与围绕着E点的材料相比可以忽略不计(算术平均公式却保留kp对ke的影响)。另一个含义是：ke不等于kE，而是它的1/fe．稍微转一下弯子就可以明白这一含义是合理的．我们的目的是通过方程(4．7)得到一个正确的qe值．应用方程(4．13)得出：

Sc＝5，Sp＝-4．这是最明显的形式，并且是我们所推荐的．

Sc＝5-4Tp，Sp＝0．这是懒汉的做没他把整个S都写成人Sc，而令Sp＝0．但是，这种做法并不是不Sc=5+7Tp，Sp＝-11．这得出比实际的S~T关系更陡的曲线．这样做的结果特使迭代的收敛速度减慢了．但

kE(TPTE)qe (4.14)

(x)e当kP>> kE时，温度TP将一直扩展到界面e，而温降TP—TE将实际上发生在距离( x)e内．于是正确的热流密度将由方程(4．14)给定．换句话说，在方程(4．13)中的因子fe可以看成是对方程(4．7)中所用的名义距离(2．

切实际的．当S的表达式很复杂时，这样做或许是我们唯一的一种选择．

3．

 x)e是倘若在所研究的问题中还存在着其它的非线性项时，这种减慢的做法实际上可能是受欢迎的．

[例2] 已知S=3+7T．某些可能的线性化是：

1．S=3，S＝7．一般说来，这是不能接受的，因为它使

cp的补偿．

对这个极限情况的讨论表明，这个公式可以适用于导热系数突然变化的情况，而无需在发生突变的邻近区域内采用极细的网格．这样做不仅对组合板的导热计算比较方便，而且还有更加引人入胜的其它一些含义．这些更深的含义已经在帕坦卡(1978)的论文中作了描述．在本书中，我们也将在随后的几章中加以说明．

所推荐的界面导热系数公式(4．9)是建立在无内热源的稳态一维导热状态基础上，其中导热系数在相邻的两个控制容积之间发生阶跃变化．即便在内热源不为零或是导热系数连续变化的场合这种调和平均的表达式也要比算术平均公式好得多。在帕坦卡(1978)的论文中已经针对某些可以求得精确解析解的情况说明了这一点． 4.2—4 非线性

离散化方程(4．2)是一个线性的代数方程

Sp为正．如果所研究的问题不用进行迭代就可以直接求解，那么这种线性化是可以给出正确的解的．但是，倘若由于某种原因(如在方程中其它有一些项是非线性的)而需要采用迭代，则存在正的

cpSp就可能会导致解的发散．

pp2．S＝3+7T，S=0．这是人们在没有自然现成的S可得的情况下，应当采取的措施．

3．S＝3+9T，S= -2．这是一种人为产生的负的S。一般说来，这样做的结果将导致收敛速度减慢．

cpaPTPaETEaWTWb (4.2)

pp[例3] 已知：S=4-5T3．某些可能的线性化是： 1．S＝4-5Tc3p我们将用解线性代数方程的方法来解这样的方程组． 非线性的情况

但是即便在热传导的问题中我们也常常会遇到非线性的情况。导热系数k可能与T有关，或者源项S可能是T的一个非线性函数．这样，离散化方程中的系数本身将与T有关．我们将用迭代的方法来处理这种状况．这个过程包括如下几步：aPTPaETEaWTWb

1．一开始在所有各个网格点上，猜测或估计一组T值． 2．由这些估计的T值，计算出离散化方程中的系数的试探值． 3．解名义上的线性代数方程组，得到一组新的T值． 4．以这些T值作为较好的估计值，返回到第二步并重复这个过程，直到这种进一步的重复计算(称之为迭代)不再引起T值任何有意义的变化时为止．

，S＝0．这是懒汉的做法，这种做法不能很好地利用已知S~T关系这一有利条件．

p2p2．S＝4，S＝-5Tcp．这看起来像是准确的线性化，但是已知的S~T曲线要比这一关系所反映的曲线为陡。

*dS***3*2*SS(TT)45T15T(TT3．推荐的方法：PPPPPP)

dT9

传热与流体流动的数值计算 传热与流体流动的数值计算

Sc410T于是：

*3P式中源项已经按通常的型式线性化．界面热流密度

这一线性化表示，在

(x)iSP15T*2PTp点，所选

qi可以按照方程(4．7)的格式写出．结果是：

qBB择的直线与S~T曲线相切。

4．S＝4+20T3，

cpSp＝-25Tp．这一线性化要地已知的S~T曲

2

ki(TBTI)qB(ScSPTB)x0 （4.16） (x)i这一方程的进一步表达形式与如何给定边界上的热流密度

iI线为陡，它会使收敛速度降低．

这四种可能的线性化与实际的S~T曲线一起示于图4．2中．在这样的图上，正斜率的直线会违反基本法则3．在所有负斜率的直线中，与已知曲线相切的直线通常是最佳的选择．较陡的斜率是可以接受的，但是通常会导致收敛速度减慢．欠陡的直线我们是不希望采纳的，因为它不能体现已知的S随T的下降速度．

就本章所要达到的目的而言，关于源项线性化的讨论到此结束． 4．2—6 边界条件

对于一维问题，让我们来讨论如图4.3所示的网格点组． 在两个边界上各有一个网格点．其余的网格点称为内点．

围绕着每一个内点有一个控制容积．对每一个这样的控制容积可以写出一个象方程(4.2)那样的离散化方程．如果把方程(4.2)看做是关于

qB有关．如果qB值本身是已知的，则所要求的对TB的方程变

成：aBTB其中：

x图4-4 边界附近的半控制容积 aITIb (4.17)

ki(x)i

aIbScxqB aBaISPx (4.18)

如果热流密度于是对

qB系由放热系数h以及环境流体温度Tf规定，那么：qBh(TfTB) (4.19)

ki(x)iTp的一个方程，那么我们就有了对所有内网格点上未知温度所必要的方程．

TB的方程为：aBTBaITIb (4.20)

aPTPaETEaWTWb (4.2)

其中有两个方程包含着边界网格点上的温度．通过处理这些边界温度，就把已知的边界条件引入到我们的数值解法中了．

由于没有必要分别讨论两个边界点，因而我们讨论左边的边界点B，如图4．3所示，该点与第一个内点I相邻．在热传导问题中有三类典型的边界条件，它们是：

1．已知边界温度 2．已知边界热流密度 3．通过放热系数和周围流体的温度来规定边界的热流密度．

如果边界温度已知(即，如果当边界温度未知时，对

式中：aIbScxhTf

aBaISPxh (4.21c)

这样我们就能够构成对所有未知温度的足够数量的方程． 现在我们该来描述求解这些方程的方法了． 4．2—7 线性代数方程的解

一维离散化方程的解可以用标准的高斯(Gauss)消去法得到，由于方程的形式特别简单，消去过程的算法就变得十分方便．有时候，这种算法称之为托马斯(Thomas)算法或是TDMA(三对角矩阵算法)．TDMA的名称来自于这样的事实，即：在写这些方程的系数矩阵时，所有的非零系数均排列在矩阵的三条对角线上．

TB的值是已知的)，这种情况并无任何特别的困难，处理时也不需要外加任何方程．

aPTPaETEaWTWb (4.2)

对方程组AX＝B

TB，我们需要构成一个附加的方程．这样的方程可以通过在图4．3所示的边界附近的“半”

控制容积内对微分方程进行积分来建立(该控制容积只包括位于网格点B一侧的半个容积，因而我们把它叫做半控制容积)．该控制容积的放大图在图4.4中给出．在该控制容积内对方程(4．1)进行积分，

ddT(k)S0 （4.1） dxdx

半控制容积 典型的控制容积

为了表示该算法的方便起见，应用某些不同的符号是必要的．设在图4．3中的网格点记标号为1， 2， 3，…，N，其中1与N代表边界点．离散化方程可以写成为：

BiIx WPEaiTibiTi1ciTi1di (4.22)

式中下标i＝l，2，3，…，N．于是温度i与相邻的温度我们令：

图4-3 内点与边界点的控制容积 TTi1及Ti1有关．考虑到边界点方程的特殊的形式，让

并考虑到热流密度q代表

kdT/dx我们就得到：

c1＝0 及 bN＝0 （4.23）

a1T1b1T2d1 T1f1(T2)

10

qBqi(ScSPTB)x0 （4.15）

传热与流体流动的数值计算 传热与流体流动的数值计算

a2T2b2T3c2T1d2 T2f2(T3) a3T3b3T4c3T2d3 T3f3(T4)

……… TN1fn1(TN) aNTNcNTN1dN 所以温度

这些关系式是递推的关系，因为它们给出了

Pi及Qi与Pi1及Qi1之间的关系．为了开始这个递推过程，我们

注意到：对i=1的方程(4．22)几乎就是方程(4．25)的样子．

T0及TN1不起任何作用．(当边界温度己知时，对这些边界点的方程就只剩下一个无意义的形式．

Ta1=1，b1＝0，c1＝0，以及d1＝T1的已知值)．

aiTibiTi1ciTi1di (4.22)

例如，如果1已知，我们就有：如果1未知，这些条件意味着

c1＝0 及 bN＝0 （4.23）

b1d1P于是P与Q的值由下式给出：1 及 Q1 (4.28)

iia1a1 [值得注意的是，在把1=0代入方程(4．27)之后，

TT1是T2的已知函数．

a1T1b1T2d1 T1f1(T2)

i＝2的方程是

cT1，T2与T3之间的一个关系式． T1，T2与T3之间的关系式就简化为T2与T3之

T2可以用T3的一个关系式来表示．

Pi间的一个关系．换句话说，

bidiciQi1Q (4.27a) iaiciPi1aiciPi1iN (4.27b)

就可以得到这两个式子．]

在P，Q序列的另一端，我们注意到b=0．这导致P=0，因此由方程(4．25，)

ia2T2b2T3c2T1d2 T2f2(T3)

这一代入过程可以一直继续到把步我们实际上也就得到了

NTN表达成TN1的一个关系式．但是由于实际上TN1并不存在，所以到这一

我们得到：TNTiPiTi1Qi (4.25)

TN值．

QN (4.29) 现在我们已经可以开始通过方程(4．25)来回代了．

这使我们可以由此开始“回代”的过程。由

TN得到TN1，由TN1得到TN2，…，由T3得到T2，以及由

算法概要

T2得到T1。这就是TDMA的要点．

a1T1b1T2d1 T1f1(T2) a2T2b2T3c2T1d2 T2f2(T3) a3T3b3T4c3T2d3 T3f3(T4)

………

1．由方程(4．28)计算P，。 1Q1aiTibiTi1ciTi1di (4.22)

P1b1d1Q 及 1a1a1 (4.28)

TN1fn1(TN)

Pi1TiQi1 (4.24)

2．应用递推关系式(4．27)得到对i＝2，3，…N，的

Pi及Qi。

(4.27)

aNTNcNTN1dN

假设在前向代入的过程中，在我们刚刚得到：Ti1以后，我们试图找到一个关系：i将方程（4.24）代入方程（4.22）

就得到：

bidiciQi1PiQi

aiciPi1aiciPi13．令

TPiTi1Qi (4.25)

TN＝QN。

aiTibiTi1ciTi1di (4.22)

4．对i＝N-1，N-2，…，3，2，1应用方程（4.25）得到

TN1，TN2，…，T3，T2，T1。

aiTibiTi1ci(Pi1TiQi1)di (4.26)

TiPiTi1Qi (4.25)

不管什么时候，只要代数方程可以表达成方程(4．22)的形式，那么三对角矩阵的算法就是一个非常有效而又方便的

求解方法．与一般的矩阵方法不同，TDMA需要的计算机贮存量及计算机时间只是正比于N，而不是N2或N3。

4．3 不稳态一维热传导

4．3—1 通用的离散化方程

这个方程可以改写成方程(4．25)的样子．换句话说，系数

Pi与Qi于是可以表示成：

(4.27b)

diciQi1biQPi (4.27a) iaiciPi1aiciPi1

11

传热与流体流动的数值计算 传热与流体流动的数值计算

()tdiv(u)div(Γgrad)S (2.13) cxt(T10PTP)fk(T1)1eET1Pkw(T1PTW))已经知道如何来处理：

(x)e(x)w一维的状态下Φ的通用微分方程中的扩散项与源项． 试探求解不稳态一维热传导微分方程 ♥ )k(T0k00 （4.35）

(1feET0P)w(TPTW)(x)e(x))wcTtx(kTx) （4.30） 在改写这个式子的时候，我们将把上标1去掉，并且记住此后

TTE和TW代表T在t+Δt时刻的新值．结

 考虑不稳态项  将源项暂时撇下不管；对于源项，我们已经没有什么新东西可以谈了． aTfT0PPaEE(1f)T0EaWfTW(1f)T 假设ρc是常数 0因为时间是一个单向坐标，我们由一已知的初始温度分布开始，沿着时间坐标逐步向前求解。 P，W果是：

aP(1f)aE(1f)a0 （4.36）

WTP计算的任务是：已知t时刻T在网格点上的值，求得在t+Δt时刻在相应网格点上的值。 在网格上“老的”(已知的)T值用T0T00akwP，E，T式中：aEkeW表示， (x)

(x) a0cxWP

afa0 (4.37)

ePEfaWaP wt在t+Δt时刻的“新的”(未知的)T值则以T1T1T14．3—2 显式，克兰克—尼科尔森(Crank—Nicolson)模式，以及全隐式模式

P，E，W表示． 对于某个特定的加权因子f的值，离散化方程可以简化为适用于抛物型微分方程的大家所熟悉的模式之一。特别是， 现在对图3．2所示的整个控制容积积分方程(4．30)以推导离散化方程，其中积分的时间间隔为由t到t+Δt． f=0导致显式模式

cTf＝0.5导致克兰克—尼科尔森模式 tx(kTx) （4.30） 以及f=1导致全隐式的模式

我们将简单地讨论这些模式并且最后指出全隐式模式是

我们所喜爱的一种模式．

于是： (x)w(x)ettcettTttT1PdtfTP(1f)T0Pt

wttdtdxttewx(kTx)dxdt（4.31） 为了表达项∂T／∂t，我们将假设：在网格点上的T值代表整个控WwPeE不同的f值可以由图4．5所示的

Tp~t变化关系来说明．

显制容积上的值．于是：

x x式模式实质上是假设老的T0cettTP值代表除了时刻t+Δt之外的整

wttdtdxcx(T1PT0P) （4.32） 图3-2 一维问题的网格节点群 个时间间隔上的TP值．全隐式的模式则假说在时刻t，TP对于k∂T／∂x．参考我们对稳态问题时所采用的措施，我们得到： 的值突然由

T0降到T1tteTttkPP，而在随后的整个时间步上保持为e(TETP)kw(TPTW)twx(kx)dxdtt(x)(x))dt T1，于是在整个时间步期间内，温度为新值T1ewPP所确定．克兰克—尼科尔森模式假设

TP呈线性变化．乍看起来，线

性变化会比其它两种模式更加切合实际．那为什么我们会偏爱全隐式模式呢？答案很快就可以清楚．

对于显式模式(f=0)，方程(4．36)变为(4．38)

cx(T10ttke(TETP)kw(TPTP)t(x)PTW))dt aa0e(x) （4.33） wPTPaEfTE(1f)T0EWfTW(1f)TWa0现在我们需要一个有关

TP(1f)aE(1f)a0）

WT（4.36PP，TE和TW如何随时间由t到t+Δt而变化的关系的假设．可以采取的假设有好几00WTP（4.38）

种，其中的某一些假设可以用下式把它们归纳一般化：tttTfT1PdtP(1f)T0PaT00PPaETEaWTWaPaEat （4.34） 这意味着

T如TT0P与其它未知量(E或TW)无关，而直接可由已知的温度T0P，T0E及W得到．这就是为什

式中f是一个在0与1之间变化的加权因子．对于T

E和TW，应用与上式相类似的公式，我们由方程(4.33)推出： 么把这种模式叫做显式的原因．对于f≠0的任何模式都是隐式的．此时，

TP将与未知量TE和TW有关，并且有必要解一组联立方程．

12

传热与流体流动的数值计算 传热与流体流动的数值计算

但是显式所具有的这一方便性为其本身所受到的一个严重限制所抵消．如果我们记住有关正系数的基本法则(法则

002)，并用它来检查一下方程(4．38)，那么我们就会发现P的系数可能变为负值(我们把P看成是P在时间方向上

0的一个相邻点值)．实际上，要使系数为正，时间步Δt就必须足够小，以使P超过aEaW．对于导热系数均

TTTTP迭代算得新的kP值。

稳态程序的其它环节，如边界条件、源项线性化处理以及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点之间的控制容积面上的热流密度

a匀分布以及Δx=（δx）e=（δx）w这个条件可以写作：

kaEe(x)e

2c(x)kwcx0aW aP (4.37) t2k(x)wt （4.39）

如果违背了这个条件；就可能出现物理上不真实的解，这是因为负的系数意味着：一个较高的T0会导致一个较低

P的Tp值的缘故．方程(4．39)是大家所熟知的有关显式模式稳定性的判定准则．

在我们需要通过减小间距Δx来改进在空间的计算精度时，我们将不得不采用一个更小得多的Δt．

克兰克—尼科尔森模式通常被描述成是无条件稳定的（对于f=0.5）。一个缺乏经验的使用者(用户)往往会就此认为，这就意味着不管采用多大的时间步，总可以取得物理上真实的解的．因此，这样一个使用者，一旦遇到解发生振荡时，就会感到惊讶．数学上的“稳定性’只能确保解的振荡最终将会消失，但是它并不能确保一定会给出物理上似乎合理的解．

对于f=0.5，方程(4．36)中的系数变为a0p(aEaW)/2．在均匀导热系数及均匀网格间距的情况下，这个系

数可以看成是ρcΔx／Δt-k／Δx。这时，一旦时间步不够小，这个系数就可能变为负值，从而有可能产生物理上不真实的结果．

这种如图4．5所示的似乎合理的线性分布是一个只对小的时间间隔有效的好的温度—时间关系表达式．在时间间隔较大时，温度的固有的指数衰减关系相当于开始急剧下降，其斜率很陡，而后则是紧跟着一个平平的尾巴．在隐式模式中所作的假设于是要比克兰克—尼科尔森模式所采用的线性分布更接近于实际，尤其是对于大时间步的情况就更是如此．

qe，我们将假设所得到的qe代表面积为Δy×1的整个表面上的

TTTc(k)(k)S (4.42)

txxyy可以立即转变成离散化方程：

值．通过其它表面的热流密度可以用同样的方式得到．这样，微分方程：

0如果我们要求方程(4．36)中的P的系数务必永不为负，唯一能够确保这一点的f值应为1(当然令f＞l是没有意义

T的)．于是全隐式模式(f＝1)就能满足我们既简单而物理上又满意的要求．鉴于这一原因，在本书中，我们将采用全隐式模式．

必须承认：对于小的时间步，全隐式模式不如克兰克—尼科尔森模式精确，其原因仍可由图4.5看出：对于小的时间间隔，温度—时间曲线是接近于线性的．

探索一种结合这两种模式的优点，而不分担它们各自缺点的模式是诱人的，实际上这一工作已经有人做了，帕坦卡与巴利加(1978)已经描述了一种名叫指数模式的新模式．但是这种模式多少是太复杂了一些． 4.3—3 全隐式离散化方程

这里我们来讨论方程(4．36)的全隐式模式．

aPTPaETEaWTWaNTNaSTSb (4.43) ♥

key式中：aE(x)e

aWkwy(x)waNcxyknxksx0aa S P

(y)nt(y)s0aPaEaWaNaSaPSPxy

bScxyaT00PP

aPTPaEfTE(1f)TaWfTW(1f)T0Pa0E(1f)aE(1f)aWT0W（4.36）

0P 乘积Δx Δy是控制容积的体积。 4．4—2 三维问题的离散化方程

最后，我们再加进两个z方向的相邻点T和B(顶和底)以构成三维的网格图形．可以很容易地看出离散化方程是：

我们引入线性化源项．其结果是：

aPTPaETEaWTWaNTNaSTSaTTTaBTBb (4.45)

keyzkwyz式中：aE aW(x)e(x)waB

aPTPaETEaWTWb (4.40)

式中：aEaNke(x)e

aWkwcx0a

(x)wtPknzxkszxktxy aS aT

(y)n(y)s(y)t000SPx bScxaPTP aPaEaWaPcxyzkbxy0a P

t(y)b00bScxyzaPTP (4.46)

可以看到，当Δt→∞时，方程简化为稳态的离散化方程。 全隐式模式的主要原则是

0aPaEaWaNaSaTaBaPSPxyz (4.46)

TP的新值代表整个时间步上的值。因此如果导热系数kP与温度有关，就应当反复由

13

离散化方程中各个系数的物理意义：相邻系数

aE、aW、aN…、aB的代表点P与相应的相邻点之间的热导．项

传热与流体流动的数值计算 传热与流体流动的数值计算

00aPTP是t时刻控制容积内部所包含的内能(除以Δt)．常数项b则由这一内能项与由SC所造成的在控制容积内的发

热率所组成．中心点系数

a0p是所有的相邻点系数之和(包括“时间相邻”的P的系数

Ta0P)，并包括一项由线性的源项

可以发现，对这个例子，不管在开始计算时初始的估计值是多少，我们都可以得到精确的解．

迭代方法的一个有趣的特征是：在中间的计算步骤上，计算的精度可以是不很高的．由于中间的计算值只是简单地用作为下一次迭代的估计值，因此在中间过程中所作的近似计算．甚至于误差，在计算结束时都将趋于消失．由下面的例子，我们可以进一步加深对这种迭代方法的了解．

方程： T1=0.4T2+0.2 T2=T1+1 (4.49) T1=0T2-1 T2=2.5T1-0.5 (4.50)

所作的贡献．

4．4—3代数方程的解

应当注意：在我们构成离散化方程时，我们只把这些方程写成一种线性的形式而并没有设想一个求解这些方程的特定方法．因此我们在这里可以采用任何一种合适的求解方法．

把方程的推导与它们的求解方法看成是两个分开的步骤是有好处的．在具体选择应用时二者无需互相影响．在计算机程序中，我们可以很方便地把这两个步骤分成两个独立的段，其中任何一段都可以根据需要单独进行修改．

直接解法：至此我们已经由一维的离散化方程直接推广得到了多维的离散化方程．不能很容易地推广到多维问题中去的一段程序是三对角矩阵算法(TDMA)．

求解二维或三维问题的代数方程的直接解法(即那些无需迭代的方法)是非常复杂的，并且需要相当大量的计算机贮存量及时间．

对于一个只需要一次求解代数方程的线性问题，直接法可能是可以接受的，但是对于非线性问题，因为方程必须用最新的系数来重复求解，采用直接法往往是不经济的．因此在进一步讨论时，我们将把直接法排除在外． 达代法：代替直接解法的是解代数方程的达代法．

定义：从一个估计的因变量T的场开始，再利用某种形式的代数方程求得一个改进的场．重复进行这一算法过程最后求得一个充分接近代数方程精确解的解。 优点：迭代方法通常只要求增加很少一点计算机存贮量就够了，对于处理非线性问题，迭代方法特别吸引人．

系数对于一个非线性的问题，在一组固定的系数值下把代数方程的求解过程一直进行到最终收敛是不必要也是不明智的．在这种情况下，在新的一组系数值约定之前应用已经给定的一组系数值对方程只要进行几次迭代就已经足够了。一般说来，似乎在计算系数所付出的代价与求解方程所得要花的时间之间应当存在某种平衡．在计算出一组系数值之后，我们必须进行充分的达代以求解方程，以便从这些系数身上得到足够的好处．但是如果把时间过多地花费在以这些只是暂时性的系数为基础来求解方程上则是不明智的．

迭代方法的种类：有许多求解代数方程的迭代方法．我们只介绍其中的两种方法：第一个方法可以看做是预备性的知识，而第二个方法则是我们将推荐采用的方法．

高斯—塞德尔(Gauss—Seidel)逐点计算法 在所有的迭代法中，最简单的一种方法是高斯—塞德尔法，其中按一定的顺序逐个访问每一个网格点，以计算那里的变量值．在计算机内只需要储存一组T值．开始，这些值代表最初的估计值或是上一次迭代所得到的值．在访问每一个网格结点时，在计算机存贮器中相应的T值交替改变如下，倘若把离散化方程写成

解：

这看起来是令人扫兴的．这里迭代过程已经发散了．更令人惊讶的是方程组(4．50)只是方程组(4．49)的一种简单改写形式，而对于方程组(4．49)，我们在例1中已经得到了收敛的解．

于是，我们可以得出结论，高斯—塞德尔法并不总是可以得到收敛的解的．斯卡巴勒(Scarborough)(1958)已经推导出一个准则公式，当这个公式得以满足时，高斯—塞德尔法就一定可以得到收敛．现在我们不加证明地，把这个斯卡巴勒准则提出来，随后我们将讨论它的含义．

斯卡巴勒准则 高斯—塞德尔法收敛的充分条件是：

说明 (1)这个准则是个充分条件而并不是个必要条件．这就是说有的时候虽然我们违背了这个条件，但迭代仍能收敛．

(2)尽管我们并不主张采用高斯—塞德尔选代法，但是我们似乎还是希望我们的离散化方程应当满足斯卡巴勒准则，以使我们至少有一种迭代方法可以确保收敛。

(3)现在可以知道我们的某些基本法则是满足斯卡巴勒准则的(这些法则系由物理的观点出发，推导而得的)．

例如，负的

SP的存在满足anb/ap1．

我们也可以由我们的正系数要求这—准则说明． 如果某些系数是负的，那么

aPTPanbTnbb (4.47)

式中下标nb代表一个相邻点，于是在被访问的网格点上的

ap(通常等于anb)可能会小于|anb|(因为在某些系数为负时，就会导致

TP由下式算得：

anb|anb|)，从而违反这一准则．

a0PcxytTP*nba*Tnbnbb

000bScxyaPTP aPaEaWaNaSaPSPxy (4.41)

aP (4.48)

*nb(4)当

ap等于anb，以及所有的系数均为正时，对所有的方程我们均得到|anb|/|ap|1．那么至少

nb其中T代表在计算机存贮器中所存在的相邻点的温度值．对于那些在本次迭代过程中已经被访问过的相邻点，T是新鲜的计算；而对于那些尚待访问的相邻点，

有一个方程有可能会使|a|/|ap|小于1的地方究竟在哪儿呢？答案是在域的边界上．对于一个具有确定解的问

nb*是由前一次迭代所得到的值．其实在任何情况下，T*都是取的相Tnbnb题，至少需要在一个边界点上规定温度值．在以该点作为相邻点的离散化方程中就包含有|a思．这是因为在计算

|/|ap|1的意

邻点温度的最新值．当所有的网格点都按这种方式访问过一次之后，就算是完成了一次高斯—塞德尔迭代．

用两个非常简单的例子说明这个方法．方程： T1=0.4T2+0.2 T2=T1+1 (4.49) 解：

|anb|时(为应用斯卡巴勒准则)，应当只对那些未知的相邻点系数求和，而另一方面ap则是

包括边界点系数在内的所有相邻点的系数之和．

高斯—塞德尔法的主要缺点是其收敛速度太慢，特别是当网格点数目很大时尤为明显．收敛速度慢的原因是易于理解的；本方法是以一次迭代传递一个网格间距的速率来传递边界条件所给予的信息的。

14

传热与流体流动的数值计算 传热与流体流动的数值计算 逐行法 现在我们可以很方便地把一维状态的直接法TDMA与高斯—塞德尔法结合起来． ♥

 选择一条网格行(设在y方向选取这样的网格行)，假定沿相邻的行(即所选择行上的点在x与z方向的所有相邻

点的集合)上的T值由其最新值构成．

 用TDMA法求解所选行上的T值．我们将在同一方向的所有各行上进行这种计算．如果想做的话，我们再按

同样的方法在其它方向重复上述的程序．



anbTnbb*TPTPTP (4.54) aP*aPTPaETEaWTWaNTNaSTSb

式中括号内的部分代表由本次迭代所产生的

Tp的变化．这一变化可以通过引进一个松弛因子加以修改，所以

讨论 (1)某一选定的行网格点上的离散化方程中包含着沿两相邻行上网格点(用叉线表示)的温度．如果这些温度用它们

的最新值取代，那么沿选定行网格点(图中用圆点表示)的方程就可以看成是一些一维方程，这些方程可以采用TDMA法来求解。按照这一方法对所有的沿y方向的其它各行进行计算，随后再按类似的方法处理在x方向的所有各行．

(2)逐行法的收敛速度比前面所提到的逐点法快，这是因为不管沿行有多少个网格点，在行两个端点上的边界条件信息都可以立即传递到域的内部．而沿另一方向的信息传递速率则类似于逐点法的传递速率．

(3)通过交替改变应用逐行迭代的行方向，我们就可以把所有边界上的信息迅速传递到中间区域．

(4)几何形状以及其它状态特性往往会导致，例如，y方向的系数远远大于x方向系数(见图4．8)的情况。

anbTnbba*TPTPTP (4.55a) 或 PTPaP*anbTnbb(1)aPTP* (4.55b)

首先，应当注意到当迭代收敛时

Tp变成与T*p相等．方程

ksxkwyknxaSaWaN(y)s (4.44) (x)w(y)n在这种情况下，当对y方向应用TDMA逐行迭代时，收

敛得特别快．这是因为需要代入离散化方程的沿相邻行上的温度估值对于离散化方程的影响不大．

(5)除了应用TDMA所选择的行方向之外，扫描方向(即选择计算行的先后次序)在某些情况下也是重要的．

对于如图4．9所示的边界条件，从左到右的扫描(即，选择域的左边界作为第一行，而后逐渐向右移动行)会把左边界上已知的温度传递到计算域内部；与此相反，由于右边界没有给定温度，自右而左的扫描不可能带进这种有用的信息(同样的考虑方法适用于逐点计算中访问点的先后次序)．当有对流存在是，扫描方向是尤为重要的．十分清楚，由上游向下游扫描的收敛速度要比按相反方向扫描的收敛速度快得多．

其它的一些迭代方法 皮斯曼(Peaceman)和拉什福德(Rachford)(1955)引入了一种常用的、叫做ADI(方向交替的隐式)的逐行求解法．另外一种用于求解多维离散化方程的迭代解法是由斯通(Stone)(1968)所述的强隐式法(SIP)．我们把这些方法的详细研究留给那些感兴趣的读者自己去做．

4．5 超松弛与欠松弛

在代数方程的迭代求解过程中，或是用于处理非线性问题的整体迭代模式中，人们往往希望加快或是减慢前后二次迭代之间因变量值的变化．依因变量的变化究竟是被加速还是被减慢而定，这一过程被称之为超松弛或是欠松弛．

超松弛常常用于和高斯—塞德尔法相结合，这种结合起来的模式是所谓的持续超松弛(SOR)。通常超松弛很少与逐行法结合起来使用．

欠松弛对于非线性的问题是十分有用的工具．在强烈非线性方程组的迭代求解过程中，往往采用欠松弛来避免发散． 引进超松弛或是欠松弛的方法有许多种．这里将描述其中的几种具体做法．我们将针对下列形式的通用离散化方程进行讨论．

(4．55a)意味着T的收敛值确实是满足原始方程(4．52)的．当然任何一种型式的松弛都必须满足这一性质；不管是通过使用任意的松弛因子或是任何类似的手段，所得到的最终收敛解都必须满足原始的离散化方程．

α值的选取、欠松弛与超松弛

当方程(4．55)中的松弛因子α在0与1之间时，它的作用是欠

*T松弛的；即Tp的值更接近于p一些．对一很小的α值，Tp的

变化变得很慢．而当α大于1时，就产生超松弛．

不存在什么选取最佳α值的一般法则．α的最佳值与许多因素有关．诸如所研究问题本身的特性，离散网格点的数目，网格点之间的间距以及所采用的迭代方法等，都是重要的影响因素．通常，可以根据经验以及对所给定的问题所作的试探性计算求得一个合适的α值．

在整个计算期间都保持采用相同的α值是没有必要的．在不同的迭代次数之间可以改变α的值。事实上对每一个网格点选用各自不同的α值。尽管这样做不太方便，也是允许的．

通过惯量进行松弛 进行超松弛或欠松弛的另一种方法是用下面的公式来代替离散化方程(4．52)：

(aPi)TPanbTnbbiTP* (4.56)

式中i是所谓的惯量．对于正的i值，方程(4.56)具有欠松弛的作用，而负的i值则产生超松弛．

同样，这里也不存在获得最佳i值的一般法则；这个最佳值必须由对一特定问题所取得的经验确定。由方程(4．56)我们可以推论：i应当与

aP不相上下，同时i的模值愈大，松弛的作用也愈烈。

有时一个稳态问题的解是通过应用相应的不稳态问题的离散化方程求得的．这样，其中的“时间步”就变成与我们前

aPTPanbTnbb (4.52)

此外，取

*00面的“迭代”的概念相同．“老的”值在这个意义上，方程(4．46h)中的项ppp。

*0与方程(4．66)中p项起着同样的作用．于是，惯量i类似于不稳态公式中的系数p．这样的类比关系为我们提供

T0p就只不过是代表上一次迭代的值

TaTiTaTp*作为前一次迭代所得到的Tp值．

TPapnb采用一个松弛因子 方程(4．52)可以写成

TnbbaP (4.53)

如果我们在方程的右侧加上T*，再减它，我们就有：

了一种确定合理的i值的方法．另一方面，现在可以简单地认为通过不稳态公式解稳态问题的做法是一种特定类型的欠

松弛．其中所选取的时间步愈小，所得的欠松弛也愈强．顺便提一句，一个负的时间步Δt值会产生超松弛．

4．6某些几何上的考虑

4.6—1 控制容积面的位置

到目前为止，还没有专门描述过应在何处放置控制容积面的具体做法问题．因为离散化方程是在一般的条件下推得的，所以这些离散化方程适用于任何特定的控制容积面构成方法．

在许多可能的构成方法中，我们将讨论其中两种不同的替代形式，并讨论它们各自有关的优点．这两种方法将称之为方法A与方法B．为方便起见，描述将针对一个两维的问题，尽管所包含的有关概念也同样适用于一维和三维的情况．

15

传热与流体流动的数值计算 传热与流体流动的数值计算 方法A：控制容积面放在两个网格点之间的中点 最显而易见的构成控制容积的方法是把它们的面放在相邻网格点之间的中点．这种布置方式示于图4．10，其中虚线表示控制容积面．图中故意把网格画得很不均匀；由图中可以看到，这样做所得的一个结果是：一个典型的网格P并不落在包围该点的控制容积的几何中心上．

方法B：网格点放在控制容积的中心 另一种方法(如图4．11所示)是先画控制容积，而后把网格点放在控制容积的几何中心．按照这种布置方式，如果控制容积的尺寸不均匀，则其表面就不会落在两个网格点之间的中点上．

讨论 （1）应当注意对均匀的网格(或是均匀的控制容积尺寸)两种方法趋于一致．因此对两种方法的比较只是在非均匀的网格间距的情况下才是有意义的．

(2)在方法A中把控制容积面放在两个网格点之间的“中点”，这为计算通过该面的热流提供了较高的精度．如在3．4节中所提到的，分段线性分布的斜率与在两个网格点之间的中点所算得的任何一根抛物线的斜率相等．

(3)另一方面，在图4．10中网格点P可能不能在控制容积的几何中心，这一事实却是一个缺点．

于是不能认为在计算源项、导热系数以及类似的其它一些量时

Tp是对控制容积的一个好的代表值．此外，甚至在计

算控制容积面上的热流密度时，方法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与内点I之间，而实际上是通过边界点．如果在边界点周围想象存在一个零厚度的控制容积，那么面i相对于网格点B和I的位置可以认为是与方法B的正常情况相一致的，按照这种布置方式，也就没有必要对邻近边界的控制容积采用特殊的离散化方程，现有的边界条件数据(如已知温度或热流密度)就可以直接用于边界面i上。

4．6—2 其它坐标系

至此我们已经推导出了应用于直角坐标系中一个网格的离散化方程。在本书的余下部分，将继续采用同样的坐标系来处理几乎所有问题．这种坐标系表达方便而且易于理解．但是所提出的这种方法不只限于直角坐标系，它还可以用于任意一种正交坐标系．为了说明在其它坐标系中离散化方程的推导，我们将讨论用r与θ所表示的二维极坐标问题．

与方程(4．42)相应的r θ形式是：

T1T1kTc(rk)()S (4.57)

trrrr在r θ坐标系中的网格与控制容积示于图4．15中．设控制容积在z方向的厚度为1．为了得到这一坐标系下的离散

化方程，我们以r乘方程(4．57)的两边，并在整个控制容积的范围内对r与θ进行积分(这一运算代表体积积分，因为r dr r dθ代表单位厚度的体积元)．按照与4．4—1小节相同的步骤，我们得到离散化方程：

aPTPaETEaWTWaNTNaSTSb(4.58)

式中：

aE

对方法B来说，以规则的控制容积填满计算域、并把边界的网格点放在邻近边界的控制容积的面上，是比较方便的．这样的布置方式示于图4．14中．

0kerre()e

aWkwrrw()w

knrnaN(r)nt

ksrsaS(r)s

aPcV000aaaaaabSVaT EWNSPSPV (4.59) cPP P这里ΔV是控制容积的体况，它等于

0.5(rnrs)   r (应当注意ΔV没有必要等于rp   r，除非

P位于n与S之间的中点)．

上述说明表明由一新的坐标系引入的补充特性主要是几何上的特性．只要所需要求的长度、面积以及体积计算得合适，那就不需要什么新的原则．于是在任何正交坐标系中的离散化方程可以按照同样的方法推导得到．但是如果所用的只是由两个网格点定义的分布曲线的话，那么正交性的要求就是必不可少的．在图4．15中，控制容积面e是与线PE

16

传热与流体流动的数值计算 传热与流体流动的数值计算 相垂直的，正是这一事实才使我们可以只用

TP与TE来计算通过控制容积面的热流密度．对于非正交的网格，就要采

式中

e与w的值要用4．2—3小节中所推荐的方法求得(这一点适用于全书，尽管在前几节中可能并没有重复这

x (5.9)

用较为复杂的离散化方程．

在本书的其余部分，我们将只采用直角坐标系来进行全部的代数推导．但是当几何形状发生明显的变化时，整个处理方法将同样适用于任何一种正交坐标系．

第五章 对流与扩散

5．1 任 务 （略 见课本P91） 5．2 一维稳态对流与扩散 对可能的最简单情况的讨论 讨论只有对流与扩散这两项存在的情况下的一维稳态问题．

控制微分方程是：

些说明)．

为了把方程写得更紧凑，我们定义两个新的符号F与D如下：Fu D因次与物理意义 二者具有同样的因次，F表示对流(或流动)的强度，而D是扩散传导性。应当注意，与D总是保持正值相反，F既可以取正值也可以取负值，依流动方向而定．应用这些新的符号，离散化方程变为：

aPPaEEaWWFeaD式中：Ee2

(5.10)

ddd(u)() （5.4） dxdxdxFaWDww2

FeFwaPDeDwaEaW(FeFw) （5.11）

22式中u代表x方向的速度． 连续性方程：

d(u)0 或 dx讨论 (1)由于连续性u＝常数 （5.5）

Fe=Fw，我们确实得到aPaEaW这一性质．进一步注意到：根据方程(5．11c)，

控制容积和控制容积面的位置

为了推导离散化方程，我们将应用如图5．1所示的三网点群．尽管控制容积面e和w的实际位置并不会影响我们的最终公式，然而假设e位于P和E之间的中点，w位于W与P之间的中点还是比较方便的．

(x)w(x)e w eWPE x

图5.1 一维问题的典型网格节点群

5．2—1 预备性的推导

在图5．1所示的整个控制容积内对方程(5．4)进行积分给出：

离散化方程只是在流场满足连续性时才具有这一性质．

(2)离散化方程(5．10)隐含着Φ分段线性分布的含义．这种形式也就是熟知的中心差分格式，并且是泰勒级数公式的自然结果．

不可接受的离散化方程 (3)考虑一个简单的例子是有启发的，其中

DeDwaPDe1及FeFw4 aEDeFFe1 aWDww3

22FeFDwwaEaW(FeFw)2 22x此外，如果E和

W的值已如，我们就可以由方程

aPPaEEaWW(a)如果因为

(5.10)

2PE3W

E＝200及W＝100，结果是p=50 (b)如果E＝100及W＝200，结果是p＝250

FwFe aWDw系数有时可22ddd(u)() d(u)dxd(d)dx dxdxdxdxdxdxwweep实际上不能落在由其相邻点值所建立起来的范围100---200之外，这些结果显然是不真实的．

(4)实际上也许我们早就可以预料到这些不真实的结果，因为方程 aEDe能会变为负值．当|F|超过2D时，与F是正或是负有关，存在着

(u)e(u)w(dd)e()w (5.6) dxdxaE或aW成为负值的可能．这将违背我们的基本法则

Φ的分布 在上一章中我们已经知道如何由对Φ的一个分段线性分布来表示项d  /d x。很自然对于对流项首

之一，从而产生灾难性的结果．

(5)此外，负系数意味着aP1()(PW) (5.7) 先想到的也是采用同样的分布曲线．其结果是：eEP w2因子1/2 出自于界面位于中点的假设；对不同的界面位置则要采用其它的内插因子．现在，方程(5．6)可以写成：

12aEaW(FeFw)小于∑|anb|，从而违反斯卡巴勒准则．这样，离散化方程的逐点解就可能发散.

这就是为什么所有以前用中心差分格式来求解对流的努力都只能限于低雷诺数(即低的F／D值)的原因． (6)在扩散项为零的情况下(即Γ＝0)， 11(P)w(PW)(u)e(EP)(u)w(PW)eE (5.8) 22(x)e(x)w

aPDe17

FeFDww22中心差分格式导致

aP=0．于是方程(5．10)就变得不适于用逐点法来求解了，并且

传热与流体流动的数值计算 传热与流体流动的数值计算 也不适于采用大多数其它的迭代解法． 由于上述预备性公式推导已经导致一个不可接受的离散化方程，因而必须寻找更好的公式．在下面的几个小节中将描述几个这样的可能方案． 5．2—2 上风方案 众所周知的一帖用于解决前面所遇到的困难的妙药是上风方案，它也叫做上风差分格式、上游差分格式以及供体(施主)室法等． 这个方案首先由库兰特(Courant)、伊萨克逊(Issacson)和里斯(Rees)(1952)提出来，随后又由金特里(Gentry)、马丁(Martin)和戴利(Daly)(1966)，巴拉卡特(Barakat)和克拉克(Clark)(1966)以及朗切尔(Runchal)和沃尔夫茨坦〔Wolfshtein)(1969)等人再次阐明． 预备性公式的弱点 上风方案认为预备性公式的弱点在于假设：界面上的对流性质“槽与管”的模型的基础上[戈斯曼、潘、朗切尔、斯波尔丁和沃尔夫茨坦(1969)]．

图5.2 槽和管的模型

如图5．2所示，可以把控制容积看成是一系列由短管连接起来的搅拌槽．通过管子的流动代表对流，而通过槽壁的传导代表扩散．因为槽是搅拌的，每一个槽中都含有温度均匀的流体．于是可以认为流过每一连接管的流体具有在上游一侧槽内流体的温度．

5．2—3精确解

很幸运，如果Г取作常数，就可以准确地求解控制方程（5．4)，如方程(5．5)所给定，ρu已经是常数．

e是E和P的平均值 (x)w(x)e上风方案提出了一个较好的处理办法．保留扩散项的公式不变，而对流项则按下列假设计算： 界面上的Φ值等于界面上风侧网格点上的Φ值．于是： eP 如果 Fe0 (5.12a) dddd(u)0 或 u＝常数 （5.5） (u)() （5.4）

dxdxdxdx

如果采用一个城0≤x≤L，并利用边界条件：

eE 如果 Fe0 (5.12b) 按类似的方法可以确定w值． Wx wPeE在

x0处，0 （5.16a） 在xL处，L （5.16b）

x方程(5．4)的解是

如果我们定义一个新的运算子的话，条件语句(5．12)可以写得更紧凑．我们特定义『A，B』代表A和B中的大者．于是，上风方案意味着： 图5.1 一维问题的典型网格节点群 0exp(Px/L)1L0exp(P)1 (5.17）

式中P是由下式定义的贝克列(Peclet)数

FeePFe,0EFe,0 (5.13) dd(u)e(u)w()e()w (5.6) dxdx用这个概念代替方程(5．7) PuL (5.18)

11(PW) (5.7) e(EP) w22FwWWFw,0PFw,0 (u)e(u)wFeeFWwPFe,0EFe,0WFw,0PFw,0 dd()e()De(EP)Dw(PW)w dxdx就可把离散化方程写成 aPPaEEaWW式中aE (5.14) DeFe,0 aWDwFW,0 aPDeFe,0DwFw,0aEaW(FeFw) 计论 (1)很显然，方程(5．15)不会产生负的系数．于是，解将总是在物理上真实的；同时，斯卡巴勒准则也将得到满足． (2)上风方案所体现的主导思想的理论基础到底是什么呢？一个关于上风方案的、明白易懂的物理图象是建立在

显然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)所推荐并应用的方案之一．

18

传热与流体流动的数值计算 传热与流体流动的数值计算

5．2—4 指数方案 总流量密度J 考虑一个由对流流量密度ρ u Φ与扩散流量密度 d/dx所组成的总流量密度J是有益的。于是：

看来没有理由在计算指数上花费过多的时间．

实际上需要的方案 我们实际上需要的是一种既具有指数方案定性特征又易于计算的方案．下面我们将介绍两种这样的方案并推荐采用其中的第二个．

5．2—5 混合方案

混合方案是由斯波尔丁(1972)制定的；这一方案也出现在帕坦卡与斯波尔丁(1970)的书中名为“高横向流的修正”的条目下． 无因次形式aE/De与贝克列数之间的曲线关系 为了评价指数方案与混合方案之间的关系，我们将画出系数

dJudx （5.19）

按照这个定义，方程5.4

ddd(u)()dxdxdxaE，或它的无因次形式aE/De与贝克列数之间的

dJ0 （5.20） 变为：dx在图5.1所示的整个控制容积内积分方程（5.20），就得到

曲线关系．由方程(5．26)我们推得：

aEPeDeexp(Pe)1 (5.27)

aE/De随Pe的变化如图5.4所示． 对正的Pe，网格点E是下游的相邻点，看来它的影响随着Pe的增加

而减小。 当

并用距离(δX)e代替L．将 P和E代替0和 L，

由图5．4所示的实线可以看到 对Pe

（5.23）

JeJw0 （5.21）

总流量密度J的分布 现在，精确解(5．17)可以用作为点P与E之间的分布，其中用这个分布关系代入方程(5.19)，就可以给出对Je的表达式

Pe为负值时，E是上游的相邻点并且具有比较大的影响．

aE/De准确变化的某些特殊性质：

(u)e(x)eFePEJeFe(P) （5.22） 式中：Peexp(Pe)1eDeaEa0 ； 对Pe EPe ； 对Pe0 aE1Pe （5.28） DeDeDe2代表这些极限情况的三条直线也示于图5．4中．可以看到这三条直线构成准确曲线的一根包络，并代表着这一准确

曲线的合理近似，混合方案实际上就是由这三条直线所组成，所以

Fe和De由方程(5．9)定义． Fu

应当注意

Dx (5.9)

aEPe； 对2对Pe2 DePe2 aE1PeDe2； 对Pe2

aE0 （5.29） DeJe与点P和E之间界面的位置无关．最后，把方程(5．22)和一个关于Jw的类似关系代入方程(5．21)

PE)Fw(WWP)0 （5.24）

exp(Pe)1exp(Pw)1应用特殊的符号『 』把这些表达式合并成一个紧凑的形式，该符号代表其内所包含的所有量的最大值．于是

F(就得到：epFePea[[F,D,0]] (5.30b) aEDe[[Pe,1,0]] （5.30a） 或 Eee22混合方案的意义可以通过观察下面两点来了解：

(1)在贝克列数为一2≤Pe≤2时，混合方案同中心差分格式一致

(2)在该范围之外，混合方案简化为上风方案，其中已把扩散项置为0，所以混合方案不再具有5．2—3小节末尾所列举的一些缺点．

“混合”这个名称就是代表着中心差分格式与上风方案的组合，但是最好还是把混合方案看成是如图5．4所示的准确曲线的三条近似线．

现在可以把混合方案的对流—扩散离散化方程写成：式中：aE[Fe,De这个方程可以写成标准形式：

aPpaEEaWW (5.25)

式中：aEFwexp(Fw/Dw)Fea W aPaEaW(FeFw) （5.26）

exp(Fw/Dw)1exp(Fe/De)1aPpaEEaWW （5.31）

方案的优点 这些系数表达式定义了指数方案．在应用于一维的稳态问题时，该指数方案保证对任何的贝克列数以及任意数量的网格点均可以得到精确解． 未得到广泛的应用的原因 尽管这种方案具有如此极为理想的性质，但由于下列原因这种方案仍然未得到广泛的应用： (1)指数运算是费时的

(2)对于二维或三维的问题，以及源项不为零的情况，这种方案是不准确的．

FeF,0] aW[Fw,Dew,0] aPaEaW(FeFw)

22应当记住这个公式适用于界面位于两个网格点之间的任何一个任意位置的情况，它不只限于中点界面的情况．

5.2—6幂函数方案 由图5．4，可见在

Pe= 2时，混合方案偏离准确曲线相当大，看来在|P|值一超过2就马上把扩散效应置为0的

做法多少有点过早了．准确曲线的一个更好的近似系由幂函数方案给定，该方案已由帕坦卡(1979a)作了描述．尽管幂函

19

传热与流体流动的数值计算 传热与流体流动的数值计算 数方案较混合方案要复杂一些，但是幂函数表达式在计算时并非是特别费时的，而这些表达式却极好地代表着指数状态．aE的幂函数表达式可以写成如下： 对于Pe10

aEaE510.1PEPe； Pe； 对于10Pe0 DeDeaE510.1PE； 对于PeDe对于0Pe10

10

aE0 （5.33） De把这些方程与方程(5．29)相比较，我们发现：对于

|Pe|10，幂函数方案即与混合方案一致．方程(5．33)的

一种紧凑形式可以写作

aEDe[[0,(10.1FeDe)5]][[0,Fe]] （5.34）

图5.5 两个网格点之间的总流量密度J

可以由表5．1来评判幂函数与准确的指数方案的接近程度； 表5.1由幂函数方案与指数方案所给定的系数值的比较 aE/De Pe -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 Pe 幂函数方案 0.7738 0.5905 0.3277 0.1681 0.07776 0.03125 0 0 aE/De 指数方案 0.7707 0.5820 0.3130 0.1572 0.07463 0.03392 0.00045 4.110-8 A(P)A(P)[[P,0]] （5.42） B(P)A(P)[[P,0]] (5.43)

如果我们现在应用流量关系式(5．37)于界面e和w，并利用方程(5．42)和(5．43)，就可以得到下列通用的对流—扩散公式：

aPpaEEaWW (5.46)

aEDeA(Pe)[[Fe,0]] （5.47a）

式中

0.5 1 2 3 4 5 10 20 aWDwA(Pw)[[Fw,0]] （5.47b）

aPaEaW(FeFw) （5.47c）

现在可以把前面所推得的各种方案看成是仅仅选择不同的函数A(|P|)而已．前面所考虑的各种方案的A(|P|)表达式列于表5．2并用图线示于图5．7中，可以通过与相应的精确解相比较的方法来评判每一种函数的满意程度．

表5.2 各种不同方案（格式）的函数A（|P|）

方案（格式） 中心差分 上风 混合 幂函数 指数（精确解）

5．2—8 各种方案(格式)的结果

在结束一维问题的讨论之前，我们要检查一下对于一定的

对A（|P|）的公式 1-0.5|P| 1 两者之间的差异是太小了，以致很难用图形来进行比较．如前面所述，尽管对许多实际问题来说混合方案用起来一样好，在本书中还是把幂函数方案作为推荐采用的对流—扩散公式．

5．2—7 一个通用化的公式

为了进一步了解对流—扩散公式并构成一个可以把迄今所考虑过的所有方案都包括进去的通用格式，我们现在来探讨一下有关系数的某些一般性质．让我们来讨论图5．5所示的由距离δ分开的网格点i和i十1．我们的兴趣在于表达通过这两个网格点之间的界面上的总流量密度J．利用方程(5．19)，

0,10.5|P| 0,(10.1|P|) 5Jdd*JP Ju （5.19） 我们写出： （5.35）

d(x/)dx式中P是贝克列数，ρuδ／Г．在界面上的Φ值将是Φi与Φi+1的某个加权平均值，J*|P|/[EXP(|P|)-1] BiAi1 (5.37)

J*PiAii1 (5.44) J*Pi1Bii1 (5.45)

E和W值由各个方案所算出的P值．让我们令值

E＝1及W=0(这样做并不失其一般性)．进一步令距离(δx)e与(δx)w相等，于是P将是P（≡ρuδx／Г)的一个函数．对

各个P值由不同的方案所算得的

P值示于图5．8中．(由幂函数方案所得到的结果与指数方案(精确解)两者是太接近了，

它们根本无法分成两条线．) 除了中心差分格式之外，所有方案都给出某种物理上真实的解；与众不同，中心差分格式

则得出某些位于由边界值所规定的0—1范围以外的值．

因为决定数值方案状态的是网格的贝克列数，原则上是有可能通过分细网格(即：采用较小的δx)直到P足够小(＜2)、从而使中心差分格式获得真实的解的．但是在大多数的实际问题中，这种对策需要极细的网格．这样做在经济上往往是

20

传热与流体流动的数值计算 传热与流体流动的数值计算 行不通的．当我们在寻求那种即便在粗的网格下也能给出物理上真实解的方法时，我们无论如何也不可能接受这样的一种约束．

将给出

(PpPp)xyt00JeJwJnJs(SCSPP)xy (5.50)

式中源项已经按常用的方法线性化． 对于不稳态的项，假定

0P与P代表在整个控制容积内的值．“老”的值(即时间步开始时的值)用 0与PP表

示．按照全隐式的做法，所有其它的值(即不带上标的值)都应该看成是“新”值．

量

Je、Jw、Jn和Js是整个控制容积面上的积分总流量，即，Je代表在整个界面e上的Jxdy，依此

(uj)0 (5.1) txj类推．

类似地，我们可以在整个控制容积内积分连续性方程：

(PP)xyFeFwFnFs0 (5.51)

t

5．3 二维问题的离散化方程

现在我们具备了描写相应于通用微分方程(5．2)的离散化方程所需要的全部素材．首先，我们将只推导二维形式的方程，不过同样的方法也适用于三维的情况．

让我们来讨论图5．9所示的控制容积．如果我们应用我们求得总流量密度

式中F、F、F和F是通过控制容积面的质量流量．如果在点e的ρu代表整个界面e上的值，我们就可以写出：

e0wnsFe(u)ey (5.52a)

类似地：

Je的一维的做法，并假设该总流量密度

NxnFw(u)wy Fn(u)nx Fs(u)sx (5.52)

代表面积为Δy × l的整个控制容积面上的值，我们马上就可以写出完整的离散化方程。

5．3—1 推导的细节

那还是在一维的情况下我们就已经知道，只有当满足连续性方程时才能证明

如果以

P乘方程(5．51)，并从方程(5．50)中减去所乘得的结果，

00(PpPp)xyt0JeJWJnJs(SCSPP)xy(5.50) aP等于aEaW．于是我

们有关相邻系数和的基本法则(法则4)只是在我们于推导中应用了连续性方程后才能满足．现在将这一做法说明如下：

方程(5．2)的二维形式可以写作：

WJwwJnPseJe控制容积yE(PP)xyFeFwFnFs0 (5.51)

t我们就得到

JxJy()S （5.48） txy式中扩散)：JJsyxS(PP)0Pxy0t(JnFnP)(JsFsP)(SCSPP)xy根据方程(5．44)和(5．45)：

(JeFeP)(JwFwP) （5.53）

Jx和Jy是由下式定义的总流量密度(对流加

xu （5.49a） x图5.9 二维问题的控制容积

J*PiA(ii1) (5.44) J*Pi1B(ii1) (5.45) JeFePaE(PE) (5.54a) JwFwPaW(WP) (5.54b)

根据

和

 （5.49b）

JyuyJxJy()S (5.48) txyA(P)A(|P|)P,0 (5.42) Pe(u)e(x)eFee (5.23)

De式中u和v代表x和y方向的速度分量．在图5．9所示的控制容积内对方程(5．48)进行积分

Fu,D (5.9) aAD

x21

传热与流体流动的数值计算 传热与流体流动的数值计算 得到：

aEDeA(|Pe|)Fe,0 (5.55a) aWDwA(|Pw|)Fw,0 (5.55b)

A(|P|)0,(10.1|P|)5 (5.60)

可以意识到：就是现在，方程(5．56)中各项系数的物理意义也是容易理解的．相邻点系数。E、

W、N、S代表在四个控制容积面上对流与扩散的影响，它们分别用流量F与传导性D表示．项

这里

De和Dw，象它们的对应变量Fe和Fw一样，包含有面e和w的面积Δ y〔见5．3—2小节中的方程

JnFnP及JsFsP，我们就可以写出离散化方程的最终形式．由于方

aaaa(5．58)〕．把类似的表达式应用于

00aPP是在时刻t控制容积内已

程(5．54)自身的特性；关于相邻点系数和的法则自然得以满足．

当已知的速度与密度场恰恰满足连续性离散化方程时，上述推导与仅仅以方程(5．50)为基础的推导所得到的离散化方程相同。但是当已知流场不满足连续性方程时，两种公式推导就给出不同的结果．出于第三章中提到的原因，我们喜爱满足我们的基本法则的公式推导．

我们怎么可能会遇到一些不满足连续性的流场呢?这种可能性产生于：由于流场常常不是实际给定的，而是通过迭代算得的。这正如在一个热传导问题中，与温度有关的导热系数总是时时更新一样．在达到最终的收敛之前，在迭代的中间阶段的不完善的流场可能并不满足连续性方程。正是由于这一点，我们一直都特别留意要满足法则4．

5．3—2 最终的离散化方程 知的Φ含量除以时间步长．其余各项可作类似的解释．

5．4 三维问题的离散化方程 到此，我们已经达到了我们的目的．我们着手来写一个建立在通用微分方程(5．2)基础上的离散化方程．现在该轮到写三维的问题了(用T与B表示在z方向的“顶部”与“底部”相邻点)：

appaEEaWWaNNaSSaTTaBBb（5.61）

式中

(PP)0Pxy0aEDeA(|Pe|)Fe,0 awDwA(|Pw|)Fw,0 aNDnA(|Pn|)Fn,0

t(JnFnP)(JsFsP)(SCSPP)xy(JeFeP)(JWFWP) aSDsA(|Ps|)Fs,0 aTDtA(|Pt|)Ft,0 aBDbA(|Pb|)Fb,0

ap0Pxyz bSxyza00 cppt00JeFePaE(PE) (5.54a) JwFwPaW(WP) (5.54b)

JnFnPaN(PN) JsFsPaS(SP)

现在可以把二维的离散化方程写成：

apaEaWaNaSaTaBapSPxyz (5.62i)

流量与传导性定义为：

aPPaEEaWWaNNaSSb (5.56)

式中：

0Fe(u)eyz Deeyz Fw(u)wyz Dwwyz

(x)e(x)waEDeA(|Pe|)Fe,0 aWDwA(|Pw|)Fw,0

Pxyt0aNDnA(|Pn|)Fn,0 aSDsA(|Ps|)Fs,0

Fn(v)nzx Dnnzx Fs(v)szx Dsszx

(y)n(y)sapbScxyapp000Ft(w)txy Dttxz Fb(w)bxy Dbbxz

(z)t(z)b贝克列数取为F与D之比，于是，Pe幂函数公式为：

aPaEaWaNaSapSpxy (5.57)

00与P对应于时刻t的已知值，而所有其它的值(P、E、W、N、S等等)都是时刻t +Δ t的P未知值．流量F、Fw、F以及F已经在方程(5．52)中定义．相应的传导性定义为

ens这里

Fe/De，依此类推．对于各种不同的方案，函数A(|P|)列于表5．2中。

A(|P|)0,(10.1|P|)5 (5.64)

wynxsxey Dw Dn Ds (5.58) De(x)e(x)w(y)n(y)sFeFwFnFs,Pw,Pn,Ps贝克列数定义为：PeDeDwDnDs (5.59)

5.5 单向空间坐标 在第二章中我们曾提到可以把坐标分为单向的与双向的两种．在应用计算机计算时，单向坐标具有某些优点． 时间是一个单向坐标，在推导一个按时间向前推进的程序时，我们已经采用了这种单向坐标。 对流—扩散的公式表明一个空间坐标也可能变成为单向的． 5．5—1 使空间坐标成为单向坐标的条件

函数A(|P|)可以按照所要采用的方案由表5．2选取．幂函数方案是推荐的，为此：

22

传热与流体流动的数值计算 传热与流体流动的数值计算

确实，上述的论点系建立在贝克列数足够大的基础上．但是在没有边界条件的其他任何信息的情况下，我们总是可以假设在出流边界处的扩散系数是小的，从而在一个大的贝克列数的条件下进行计算．一个诸如此类的稍稍歪曲实际的假设是我们在对出流边界不知道任何进一步信息而又想取得有意义解的情况下所必须采取的．所得结果的不精确性(如果说存在那么一点点的话)是我们为把计算域与出流边界下游空间隔绝开所付出的代价．

如果由于某种原因忽略出流边界上的扩散变成为严重的问题时，那么我们就可以肯定：那个分析家必定是把出流边界放在一个不合适的位置上了．重新布置边界通常可以使出流的处理成为可以接受的，一个特别坏的出流边界位置的选择是这样的：它的一部分存在有“入流”．这种情况的一个例子示于图5．12中．对于这样的一种坏的边界选择，所得到的解是没有意义的．

现在复习一下处理有关对流—扩散问题中边界条件． 只要没有流体流过计算域的边界，那么流过该边界的边界流量密度就是一纯粹的扩散流，因此可以采用第四章所描述的方法．

对于流体流入计算域的那些边界部分，通常Φ的值是已知的(如果我们不知道流体流所携入的Φ值时，问题是不确定的)．

流体离开计算域的那些边界部分形成出流边界．对这部分边界的处理办法，我们已经在上面作过讨论了． 5．6 假 扩 散 在这一书中，我们将讨论一个在数值分析的专业人员中引起相当大争论的、混淆和误解的题目．存在着某种称之为“假扩散”的情况，这个名词常常被入误解．但是假扩散就其自身的意义上讲，确是反映着大多数对流—扩散公式所存在的一个主要弱点．

5．6—1 关于假扩散的一般观点

在文献中普遍知道这样的一些论点，如：(1)中心差分格式具有二阶的精度，而上风方案只具有一阶的精度；或者(2)上风方案引起严重的假扩散．其言下之意就是中心差分格式比上风方案好．

确实，从泰勒级数展开，我们可以看到中心差分格式具有一个（Δx）3阶的舍入误差，而上风方案具有一个（Δx）2

阶的舍入误差．但是由于在对流—扩散问题中所产生的Φ一x变化是指数的，除了极小的Δ x值(或者最好用相应的贝克列数)之外，在任何情况下，泰勒级数不再是一个好的表达式．在较大的Δ x值时(这是在大多数实际问题中人们所能做到的)，泰勒级数的分析给人以误解，此时，正如我们已经知道的那样，上风方案给出比中心差分格式更为合理的结果．

如果我们把中心差分格式与上风方案的系数[方程(5.11)与(5.15)]作一比较， 中心差分格式的系数

由图5．4或5．6我们已经知道当贝克列数大时，下游相邻点的系数变小．当贝克列数超过10时，幂函数方案将取下游相邻点的系数为0(混合方案在贝克列数大于2时，取下游相邻点的系数为0)．假设在图5．10所示的二维状态中沿正的x方向存在较大的流量．那么对于沿一根y方向的线上的所有网格点P，系数aE均为0．换句话说，P将与

W、而与E无关。这样，由于在任何点上的Φ值将不受x方向下游值的影响，x就成为一个单向的坐标．于N和S 有关，

是在x方向就有可能构成一个前进解的程序．

0(PP)0Pxyt(JnFnP)(JsFsP)(SCSPP)xy

(JeFeP)(JWFWP) JFa()

eePEPE图5.10 具有一个单向空间坐标的问题 （右图）

即便是一个空间坐标就整个计算域而言并不是单向的在实际处理边界条件时，也往往应用其局部的单向特性．这一点将在下面讨论．

5．5—2 出流边界条件 在流动出口的边界上(即流体离开计算域的地方)，人们通常既不知道Φ的值，又不知道它的流量．例如在图5．11所示的流动出口边界上，人们可能不知道温度或热流密度．那么，我们能够用怎样的办法来求解呢?回答却是惊人地简单：在流动出口的边界上，不需要有关边界条件的任何信息．对于所有与出流边界相邻的网格点P，如果贝克列数足够大，系数

NWPSEFeFwFw Fe

apDeDwaEaW(FeFw) aWDwaEDe2222上风方案的系数

aEDeFe,0 aWDwFw,0

apDeFe,0DwFw,0aEaW(FeFw)

FeFeu0,aEDe22FeFe u0,aEDe22

aE将为0，因而系数乘边界值将为0，也就

y没有必要知道边界值．换句话说在出流边界附近的区域，对大的贝克x列数而言，呈现局部的单向状态，因为边界点在计算域的下游，它们对解没有影响．

就可以看到：上风方案相当于用Г+ρuδx/2取代中心差分格式中的Г．换句话说，上风方案似乎把真正的扩散系数扩大了一个虚拟的(因此是假的)扩散系数ρuδx/2，于是，这一人为的扩散系数的引入被人认为是不精确的，是实际情况的一个错误表达式，因此是坏的．此外需要注意这个争辨的麻烦是建立在把中心差分格式假设为精确而又标准的(或者说是把泰勒级数认为是可靠的)基础上，并且用这个参考的框框来观察上风方案．按照这种做法，我们就会发现即便本身就是精确解的指数方案也有某种假扩散．相反，由本章所介绍的理论可以很出结论：所谓的假扩散系数ρuδx/2实际上是在大的贝克列数条件下的一种理想的补充，因为它实际上修正了由中心差分格式可能带来的错误结果．

毫无疑问，在很小的贝克列数时，中心差分格式比上风方案准确．这一点已经在许多图中显示出来了．我们喜爱的一些方案，如指数、混合以及幂因数方案，实际上在很小的贝克列数下与中心差分格式相符．在任何一种方案的情况下，对于很小的贝克列数，假扩散永远也不会是个严重问题．这是因为相比而言，这时实际的扩散相当大．正是在大贝克列数的情况下，假扩散的问题才变得重要了。在这种情况下，中心差分格式几乎不能给出任何合理的解，而我们在前面所讨论的所有其它方案却都给出几乎同样的结果．由于这个原因，我们余下的讨论将集中在很大的贝克列数以及上风方案

23

传热与流体流动的数值计算 传热与流体流动的数值计算 上，但是所得的结论将同样适用于指数、混合以及幂函数方案． 5．6—2 有关假扩散的正确看法 在明白了关于假扩散的一般看法确实给人以误解之后，现在我们回过头来看看到底什么样的情况可以描述成假扩散．首先要认识的是：假扩散是一种多维的现象，在稳态的一维问题中绝对不会有相应形式的假扩散(不稳态的一维问题确实也存在某种形式假扩散的干扰，但是我们将把我们的注意力集中于稳态的问题)． 为了把假扩散的正确意义具体化，让我们来讨论在图5．13中所示的状态．两股速度相等而温度不同的平行流相遇．如果扩散系数Г不等于0，就将形成一温度逐渐由高温到低温变化的混合层，该层的横向厚度就将沿流线方向而增大．另一方面，如果扩散系数Г为0，就不会形成混合层，并且在流线方向将维持温度的不连续性．用于观察假扩散的最好状况就是把真实扩散定为0的状况．如果对Г＝0情况的数值解产生一个逐渐变化的温度分布(这是非零Г的一个特征)，我们就可以得出结论：该数值解方案引起假扩散． pw （5.65）

结果，在每一条水平线上上游的已知值将构成在该线上所有点上的值．这样，上游温度分布的不连续将保存下来．因此，这里没有假扩散的现象发生．

2．在与网格线成45º的方向上的均匀流动 当把网格线布置成与流动方向成45º倾角来求解同样的问题时，情况就大大变化．为了方便起见，让我们采用Δ x＝Δ y的均匀网格．在x与y方向流动速度相等．结果是：上游相邻点的系数(

aW和a)变得相等，而下游相邻点的系数aE和aN变为0．于是我们有：aNaE0 aWFwux

SaSFsuy aPFwFSuxuy uxuy

appaEEaWWaNNaSSb （5.56）

p0.5w0.5s （5.66）

100 图5.13 在扩散存在与不存在时的温度分布 （a）对于Г＝0，中心差分格式会导致10093.7581.2565.625500, (b) 0 10087.568.755034.375aP=0。因此不能采用常用的迭代方法来求解代数方程．如果试图用直接法来求解Fw2 100755031.2518.75方程，那么要不是得不到唯一解，就是产生非常不真实的解． aEDeFe2100507512.56.25 aWDwFeFwapDeDwaEaW(FeFw)热 流22 上风方案的含义 我们现在来试一试对两个不同的网格方向用上风方案来求解图5.13b所示的问题． 1．在x方向的均匀流动 让我们来讨论图5.14所示的问题．流动沿着x方向，并且在左边界上具有已知的阶跃不连续的温度．由于Г＝0以及在y方向没有流动，因而系数aN与aS为0． 下游相邻点的系数aE也为0．于是aP必定等于冷 流00WSW0NPS0E0 图5.15 流动与网格线成45倾角的情况

对于图5．15所示的网格点，可以设左边界温度为100º，底边界温度为0来表示不连续．在内点上所得的解写在每

个网格点的近旁．要是没有假扩散的话，我们本来会在通过左下角的对角线的上方得到值100，而在对角线的下方得值0．相反，所得到的实际的解却的确代表着一个逐渐变化的温度分布，它非常象是在图5．13a中所示的那种情况．

要点：(1)当流场方向与网格点线成一倾斜角，并在与流动方向相垂直的方向上存在有非零的因变量梯度时，就会有假扩散出现．(2)对二维情况的假扩散系数的一个近似表达式已经由德瓦尔戴维斯(De Vahl Davis)和马林森

aW，这就导致： Uxysin2(Mallinson)(1972)给定，这个公式是：假 （5.67）

4(ysin3xcos3)式中假是假扩散系数，U是合速度，θ是速度向量与x方向之间的夹角(在0与90º之间)．由该方程可以很容易地知道：当合流动方向与其中的一组网格线相重时，不存在假扩散；而当流动方向与网格线之间的夹角成45º时，假扩散最为严重．

(3)减小Δx和Δy可以减小假扩散的大小；此外，只要有可能，应当把网格线布置在接近于流动的方向。 (4)因为在许多问题中存在着实际的扩散，应使假扩散与实际值相比足够小．

aEDeFe,0 aWDwFw,0 apDeFe,0DwFw,0aEaW(FeFw) （5.15c） 图5.14 沿着x方向流动的情况

24

传热与流体流动的数值计算 传热与流体流动的数值计算 (5)采用中心差分格式不是解决假扩散的灵丹妙药．如前所述，在贝克列数大时，中心差分格式会给出极不真实的解． (6)产生假扩散的基本原因是由于把流过每一控制容积面的流动处理成局部一维流动这样一种做法．对于在图5．15中插图所示的状态，由斜流对流到网格点P的Φ值实际上来自于角上的网格点SW，但是这一对流却被处理成分别来自网格点W与S的两股分开流动的作用结果．

(7)那些可能给出较小假扩散的方案应当考虑流动的多维特性．在离散化方程中包括进更多的相邻点也许是必要的．尽管已经研究出几个这样的方案〔如雷思比(1976b)〕．并且这些方案已经表明假扩散明显地减小了，但是这些方案明显地要复杂得多，并且迄今所经受的检验也还很不够．为此，我们将不在此讨论这些方案．

(8)有关假扩散的比较详细的讨论已经由雷思比(1976a)给出．

第六章 流场的计算

6．1 制订一个特殊程序的必要性

6．1—1 主要的困难

在第五章中，我们推导了在已知流场的情况下求解对Φ的通用微分方程的公式．但是除了某些特殊的情况之外，流场是不可能已知的； 相反，我们必须由相应的控制方程计算局部的速度场与密度场。

为取得收敛解的麻烦的原因所在．

当希望得到压力值时：这个量已经被如此巧妙地消去了，往往是人们所想要取得的一个重要结果或是为了计算密度和其它流体物性的一个中间结果。此外由涡量求出压力所付的代价会把本方法在节省计算机时间与储存量方面的优点抵销掉了．

然而更为严重的缺点是，本方法不能方便地外推到三维的情况，因为在这种情况下流函数并不存在．可是大多数实际问题都是三维的问题，一个只能用于二维问题的方法要在这方面应用显然要受到极大的限制．

对于三维的问题，一个以涡量为基础的方法采用六个因变量即三个涡量向量分量和三个速度势向量分量[例如，参看阿齐兹[Aziz)和黑勒姆斯(Hellums)的论文(1967)]．于是其复杂性实际上超过直接处理三个速度分量和压力本身．此外，涡量向量与速度势向量所包含的概念要比速度分量与压力的意义难以理解和接受。为使我们保持推导出物理上有意义而又明了的方法的愿望，我们需要寻求一种采用所谓的基本变量(即速度分量和压力)的方法．

本章的主要任务 将连续性方程中的间接信息转换成用于计算压力的一个直接算法．这个方法会带来几个次要的困难，在我们着手解决这个任务之前我们将讨论一下这几个困难．

6．2 某些有关的困难

6.2—1 压力梯度项的表达

(u)div(uu)div(gradu)BxVx （2.11） t(u)pdiv(uu)div(gradu) （2.11） tx速度分量受动量方程控制，而动量方程则是对Φ的通用微分方程的特殊情况(取Φ＝u，Г＝μ，等等)．

于是我们便想得出结论：我们已经研究出了求解动量方程、因而得到速度场的方法．事实上，问题并不如想象的那么简单

困难在哪儿呢?

如果说动量方程的非线性是一个困难的话，我们只要提醒一下我们自己，在处理热传导的问题时，我们是如何用迭代法来处理非线性问题的就够了．特别要指出的是，作为动量方程因变量u的一个函数的对流系数ρ u与作为温度T的一个函数的导热系数，这两者之间是没有什么不同的．从一个估计的速度场开始，我们可以迭代求解动量方程，从而得到速度分量的收敛解．

计算速度场的真正困难在于未知压力场．压力梯度构成动量方程中源项的一部分．然而还没有一个可以用来求得压力场的明显方程．对一个给定的压力场，求解动量方程确实并没有什么特别因难．可是确定压力场的途径则似乎是相当迷惘的．

压力场间接地通过连续性方程规定

(u)pdiv(uu)div(gradu) （2.11） tx对于如图6．1所示的一维情况，如果我们着手构成x方向动量方程的离散化公式。 那么唯一的新问题：dp/dx项对整个控制容积的积分表达式．

pwpe。

（6.1）

该项对离散化方程的贡献是压力降

可以假设压力是分段线性分布的． 如果把控制容积面e和w选在相应网格点之间的中点 我们就有：pwpepWpPpPpEppEW222uj0 （5.1）

txj(u)(v)(w)(u)0 （6.20） txyz这意味着动量方程将包含有两个相间而不是相邻网格点之间的压力差． (x)w意味着：压力取的粗的网 就效果上讲压力取的是较实际所采用的网格为粗的网格上的值．这一点导致解的精度下降． w但是存在着比这一点远为严重得多的另一个含义．这一点可由图W6．2最为清楚地看出，图中假设压力场用网格点处的压力值表示．不能认为这样100x 一个锯齿形的压力场是合乎实际的；但是由于各处相间网格上的压力值相等，在任5(x)ePe1005300273002730027100510051005300273002730027Ex273002730027300当正确的压力场代入动量方程时，所得到的速度场满足连续性方程．

直接求解由动量方程与连续性方程所推得的整个离散化方程组．

因为即便是单个因变量的离散化方程组，我们也都乐意采用迭代法来求解，因而直接求解速度分量以及压力的整套方程似乎已经超出了我们所考虑的问题范围了。

6．1—2 以涡量为基础的方法

由确定压力所带来的困难已经导致人们提出若干从控制方程中消去压力的方法．于是在二维的问题中，通过交叉微分从两个动量方程中消去压力，就可以导出一个涡量输运方程（这一推导的有关要点将在习题6、1中概述）．把涡量的定义与二维稳态情况下的流函数的定义相结合就是众所周知的“流函数／涡量法”，除了其它许多人之外要特别指出的是迈克斯(1963)、弗罗姆(Fromm)和哈洛(Harlow)(1963)、皮尔逊(Pearson)(1965)、巴拉卡特和克拉克(1966)以及朗切尔和沃尔夫茨坦(1969)等人都描述了这一方法，随后，戈斯曼、潘、朗切尔、斯波尔丁和沃尔夫茨坦(1969)等人所写的书把这一方法讲得更为清楚，并使该方法变得更易于接受了．

流因数／涡量法具有某些吸引人的优点．其中，压力消失了，代替处理连续性方程和两个动量方程的问题，我们只需要解两个方程以求得流函数和涡量。某些边界条件可以相当容易地规定：当一个无旋的外流位于计算域的邻近时，可以方便地取边界的涡量为0．但是流函数／涡量法具有某些比较突出的缺点．在壁面上的涡量值很难给定，而这往往成

一网格点P处，所对应的 图6.1 三点网格节点群 pwpe=0．这是一个多么荒唐的结果：动量方程100100510055100对这样一个波形压力场的“感受”竟然与对均匀的压力场的“感受”一样． 波形压力场=均匀的压力场 P=100500100500100500 图6.2 锯齿形压力场 在二维的情况下，可以更为明显地看清这个困难．正如x方向的动量受yx5pwpe的影响一样，y方向的动量则受pspn影响；于是压力pP对两

个方向的动量都没有影响．记住这一点我们可以得出结论在图6．3中所示的压力场(该压力场由布置成棋盘形的四个任图6.3 棋盘形压力场 意的压力值所组成)不会在x或y方向产生压力合力。于是一个高度不均匀的压力场就被动量方程的这种特殊离散形式处 理成为一个均匀的压力场了。由于动量方程对于这样一种不均匀的压力场毫无反应，所以一旦在迭代求解的过程中形成

25

传热与流体流动的数值计算 传热与流体流动的数值计算 这样的压力场，它就会一直保留到解完全收敛时为止。

应当注意在图6．2与6．3中所采用的实际数字并不具有任意特殊的意义，它们只是表示一个可以由任何任意数字构成的图形而己．容易想象三维状态可能会存在一个甚至更为复杂的图形，在这种情况下动量方程也还是把它当作一个均匀的压力场看待．

如果得到一个相当光滑的压力场的解，那么在该解上附加上一个棋盘型的压力场就可以构成任意数量的附加解．由于棋盘场包含着没有压力合力的作用的意思，动量方程将不会受该附加压力场的影响．得出这种荒谬的解的数值方法是不理想的．

6．2—2 连续性方程的表达

当我们试图构成连续性方程的离散化公式时，一个与上述动量方程相类似的困难发生了。对于一维稳态的常密度流场，连续性方程简化为：

而与控制容积面是否在两网格点之间的中点上无关．

很容易看出究竟应该怎样来确定速度分量v和w的位置的．在图6．6中表示的是一个二维的网格点图形，其中u和v的位置在各自的控制容积表面上．相应的三维图形可以由上述一维和二维的图形直接想象出来．

交错网格的直接结果是：无需对有关的速度分量进行任何内插运算，就可以计算出通过控制容积表面的质量流量(在第五章中出现的各个F值)．但是这一特性，尽管为建立通用的对Φ的离散化方程提供了某些方便，却还算不上是交错网格的一个重要优点．

(u)(v)(w)du(u)0 （6.20） 0 （6.2） txyzdx如果我们在图6．1所示的控制容积内积分这一方程，我们就有：

ueuw0 （6.3）

uW0 （6.5）

与前面对p一样，对u采用分段线性分布并取控制容积面于中点位置就导致：

uPuEuWuP0 （6.4） 或：uE22yx图6.5 对u的交错位置

交错网格具有两方面的重要优点．

对于典型的控制容积(如图6．6中的阴影所示)，很清楚，离散化的连续性方程将含有相邻速度分量的差，这样就避免了诸如图6．4所示的那种波形速度场会满足连续性方程的情况．在交错网格的情况下，只有“合理”的速度场才有可能满足连续性方程．

交错网格的第二方面的优点是：两相邻网格点之间的压力差现在成了位于这两个网格点之间的速度分量的自然驱动力．结果如图6．2与图6．3所示的那种压力场不再会被当做均匀的压力场了，并且也不可能成为有可能出现的一些解了．

于是，离散化的连续性方程要求相间点而不是相邻点上的速度相等．其结果，如图6．4所示的那样，完全不合乎实际的速度场却的确满足离散化连续性方程(6．5)．在二维和三维的情况下，对所有的速度分量都可以构成类似的分布图形，这些速度分量满足连续性方程但却不能认为是合理的或有意义的解．

u=100400100400100400x

图6.4 波形的速度场

在推导有关速度分量与压力的数值方法的公式之前，必须解决这些困难．在文献中，可能发现有些方法并没有对这些困难予以特别的注意．其中，通过在边界上所作的某些特殊处理，通过对边界条件的超规定说明，通过相对于一个光滑的初始估值的欠松弛，或是由于交好运都可避免可能出现的不合乎实际的解．但是这些方法中的大多数会把图6．2—6．4所示的那些类型的压力场与速度场接受为满意的解．因此，如果没有什么特殊的对策的话，总是存在着得到这种解的危险．

在我们着手去描述一个解决这些困难的方法之前，注意到在数值分析中造成麻烦的祸根似乎与一阶导数有关这一点是有趣的．二阶导数的“表现”总是好的，它不会产生什么困难．另一方面，在第五章中所遇到的所有复杂问题可能全都是起因于代表对流项的一阶导数；而在这里，动量方程中的压力的一阶导数以及连续性方程中的速度的一阶导数，都会引起相当大约麻烦．

6．3 一种解决困难的妙法——交错的网格

如果我们不一定要在同样的网格点上计算所有的变量的话，前面所述的这些困难就可以迎刃而解．如果愿意的话，我们可以对每一个因变量采用不同的网格．当然，如果这样做不会带来任何好处的话，我们是不会这样做的．然而把速度分量放在与所有其它变量不同的网格上会给我们带来相当大的好处．这个好处就是6．2节中所述的困难将完全消失，

这样一种对速度分量的位移式或“交错式”的网格，是首先由哈洛和韦尔奇(Welch)(1965)在他们的MAC方法中采用的．随后又在哈洛及其同事们所发展的其它方法中得到应用．它构成了卡雷特(Caretto)、柯尔(Curr)和斯波尔丁(1972)的SIVA程序以及帕坦卡和斯波尔丁(1972a)的SIMPLE程序的基础．

采用交错网格时，速度分量是对位于控制容积表面上的点进行计算的．于是，x方向的速度u要在垂直于x方向的表面上进行计算．在图6．5中，

u的位置用短箭头表示

而网格点(以下称主网格点)用小圆圈表示。虚线表示控制容积面。

需要注意，相对于主网格点而言，u位置只在x方向是错位的。换句话说，u的位置落在x方向连接两相邻主网格点的联线上．u的位置是否准确地位于两网格点之间的中点与控制容积的定义方法有关．u的位置必须落在控制容积面上，

yx

图6.6 对u与v的交错位置u;v;Ο其他的变量

26

传热与流体流动的数值计算 传热与流体流动的数值计算 于是在6．2节中所述的困难可以归因于；在相同的网格点上计算所有的变量这样一种做法．采用交错的网格之后，这些困难就完全消除了．

消除这些困难有它自己的代价．一个以交错网格为基础的计算机程序必须提供有关速度分量位置的全部指示值与几何信息，而且必须进行相当烦琐冗长的内插．然而交错网格所带来的好处是值得付出这样的代价的．

6．4 动量方程

我们再次提醒读者，如果压力场已知的话，那么就可以应用在第五章中对通用变量所推得的公式求得动量方程的解．在动量方程中，Φ代表有关的速度分量，Г与S则需附以适当的意义．采用交错网格就使离散化的动量方程与在主网格点处计算的其它变量Φ的离散化方程略有不同．这种不同是由于对动量方程采用交错的控制容积引起的，它只是细节性的，而不是本质的．

x方向动量方程 一个交错的控制容积如图6．7所示．

该控制容积相对于围绕主网格点P的正常控制容积的位置被错开了． 错位只发生在x方向

垂直于x方向的面通过主网格P和E． 这样的布置实现了交错网格的主要优点之一，差值

aeueanbunbb(pppE)Ae anvnanbnbvb(pppN)An （6.7）

NnPyx

图6.8 对v的控制容积

其中(pPpPpE可以用来计算作用在速度u的控制容积上的压力合力。

pN)An为相应的压力合力．对于三维的情况，可以写出一个类似的对速度分量w的方程．



动量方程只有在压力场已知或是按照某种方法估计出来的情况下才能求解，除非采用正确的压力场，所得的速度场将不会满足连续性方程．用u、v及w表示这样的一个以估计的压力场速度场将由下列离散化方程的解求得：

图6.7 对u的控制容积

在计算图6．7所示的u的控制容积面上的扩散系数与质量流量时，需要进行适当的内插，但是这种内插可以采用基本上与第五章中所描述的公式相同的式子进行．所得到的离散化方程可以写成：

p为基础的不完善的速度场．这种“带星号的”

aeueanbunbb(pppE)Ae （6.8）

anvnanbvnbb(pppN)An （6.9）

****(u)pdiv(uu)div(gradu) （2.11） txaeueanbunbb(pppE)Ae （6.6）

相邻点项的数目与问题的维数有关．对于图6．7所示的二维情况，u的四个相邻点示于控制容积的外边；对于

三维的情况，应当包括u的六个相邻点．

相邻系数anb系按控制容积面上的复合对流—扩散效应来计算．

项b的定义与方程(5．57)或(5．62)中的定义方法相同，但是压力梯度不包括在源项量SC与SP中．

压力梯度构成方程(6．6)的最后一项．由于压力场也要到最后计算，把压力项隐含在动量源项中可能会不太方

便．项(pPpE)Ae是作用在u控制容积上的压力合力，Ae是压力差作用的面积．对于二维的问题，Ae将是Δy×l，而在三维的情况下，

anwtanbw**nbb(pppT)At （6.10）

**在这些方程中，速度分量与压力已经加了上标“*”．可以注意到，位置t位于网格点P和T之间的z方向的网格线上．

6．5 压力与速度的修正 我们的目的是寻找到一个改进压力估计值

p的方法以使所算得的带星号的速度场将逐渐地接近于满足连续性方程．

让我们假定正确的压力p由下式得到：p=p*+p’ （6.11） 其中p称之为压力修正． 我们需要知道速度分量是怎样来响应压力的这一变化的． 按照类似的方法，可以引入相应的速度修正

u、v以及w：u=u*+u’ v=v*+v’ w=w*+w’ （6.12）

如果我们从方程(6．6)减去方程(6．8)，我们就得到：

Ae将代表ΔyΔz。

aeue'anbunb'b(pp'pE')Ae （6.13） aeueanbunbb(pppE)Ae （6.6）

aeueanbunbb(pppE)Ae （6.8）

其它方向的动量方程 可按类似的方法处理．图6．8表示对y方向动量方程的控制容积，它在y方向是错位的．对vn的离散化方程可以写成：

27

传热与流体流动的数值计算 传热与流体流动的数值计算

anvnanbvnbb(pppN)An （6.9）

****根据全隐式的做法，速度与密度的新值(即在时刻t+Δ t那些值)被认为在整个的时间步内都有效，老的密度p (即在时刻t的密度)只是在处理项0anwtanbwnbb(pppT)At （6.10）

****/t时才出现．

在此，我们将大胆地决定从方程中摒弃项

anbunb。

按照这些规定，方程(6．20)的积分形式变为：

可以把这样的做法看成是仅仅为了计算上的方便，而对这种做法最好不要太予注意。 在6．7—2小节中我们将要进一步讨论略去项这样做的结果是：

anbunb所产生的影响．

(PP)xyz[(u)e(u)w]yzt （6.21）

[(v)n(v)s]zx[(w)t(w)b]xy0如果我们现在用诸如方程(6．17)--(6．19)那样一些速度修正公式来代替所有的速度分量

0Aaeue'(pp'pE')Ae （6.14） 或 ue'de(pp'pE') （6.15） 其中：dee

ae方程(6．15)将称之为速度修正公式，这个公式也可以写作： ueuede(pp'pE') （6.17）

这个公式说明，相应于压力的修正，带星号的速度ue该如何修正以产生ue． 类似地，其它方向的速度分量修正公式可以写成：

*ueuede(pp'pE') （6.17） vnvndn(pp'pN') （6.18） wtwtdt(pp'pT') （6.19）

在重新整理之后，我们即得到下列对

***p的离散化方程：

vnvndn(pp'pN') （6.18） wtwtdt(pp'pT') （6.19）

这样，我们就完成了获得

6．6 压力修正方程

我们现在来把连续性方程转换为压力修正方程 假设密度ρ与压力之间没有直接的关系．稍后一点我们再讨论这个假设的内在含义． 推导是对三维情况的；一维或二维的形式可以很容易地由此得到．

**aPPP'aEPE'aWPW'aNPN'aSPS'aTPT'aBPB'b （6.22）

其中:

p的离散化方程而需要进行的一切准备工作．现在我们该转到这个任务上来了．

aEedeyz aWwdwyz aNndnzx

aSsdszx aTtdtxy aBbdbxy

0aPaEaWaNaSaTaB

(u)(v)(w)(u)0 （6.20） 连续性方程是：txyz我们将在图6．9所示的阴影部分的控制容积内积分这个方程(为方便起见图中只表示了一个二维的示图)．记住，同

样的控制容积曾被用来推导对通用变量Φ的离散方程．

(P)xyzbP[(u*)w(u*)e]yzt

[(v*)s(v*)n]zx[(w*)b(w*)t]xy界面上的密度 因为通常只计算在主网格点处密度ρ的值，界面上的密度(如式，由网格点处的值计算得到．不管采用什么样的内插公式，都必须保持

e)可以来用任何一种方便的内插公

NvnWuwyx图6.9 对连续性方程的控制容积

为了对项e在其界面所属的两个控制容积内连续(参

PvsSueE看第三章中的基本法则1)。

由方程(6．23h)可见：压力修正方程中的项b实际上是按带星号的速度取值的离散化连续性方程(6．21)左侧的负值．如

0果b值为0，那就意味着带星号的速度值与现时的pp一起确实满足连续性方程，从而不再有必要对压力进行进

/t进行积分，我们将假设密度p代表整个控制容积上的值．

28

此外，假设一个位于控制容积面上的速度分量(如ue)支配着整个面上的质量流量。

一步的修正．于是，项b代表一个“质量源”，该质量源必须由压力修正(通过与其相应的速度修正一起)去消除．

至此，我们已经推导出来为求得速度分量与压力所需要的全部方程．现在我们仍应该来全面地看看整个的求解算法了。

6．7 “SIMPLE”算法

我们现在正要研究的用以进行流场计算的程序已被定名为SIMPLE，其全称是Semi—Implicit Method for Pressure—Linked Equations，意思是解压力耦合方程的半隐式法．稍后一点我们将来讨论这一名称的意义。帕坦卡与斯波尔丁(1972)，卡雷特、戈斯曼、帕坦卡和斯波尔丁(1972)以及帕坦卡(1975)都已对这一程序作了描述． 6．7—1 计算进行的顺序

主要的计算步骤，按照它们执行的先后次序是： ♥

传热与流体流动的数值计算 传热与流体流动的数值计算

(1)．估计压力场

p．



提处理该方程的边界条件的某些要点是适宜的．

通常在边界上有两类条件．(1) 已知边界上的压力(速度是未知的)，(2)规定垂直于边界的速度分量。

给定边界上的压力 如果规定边界上的压力场估计值p等于给定的压力值pgiven，那么在边界上的p值将为0，于是这就类似于在热传导问题中给定温度边界条件的情况．

给定边界上的法向速度 如果这样来设计网格，以便该边界与一个控制容积面相吻合，这种情况看起来就象是图6．10所示的那种状态．速度

(2)．求解动量方程[如方程(6．8)--(6．10)]以得到u、v及w．

aeueanbunbb(pppE)Ae （6.8）

anvnanbvnbb(pppN)An （6.9）

****ue是已知的．

*p方程时，通过边界面上的流量不应当用ue以及相应的修正来表示，而是应当用ue本

anwtanbwnbb(pppT)At （6.10）

****在推导对所示控制容积的身来表示．

于是，在

(3)．解p方程．aPPP'aEPE'aWPW'aNPN'aSPS'aTPT'aBPB'b （6.22）

aEedeyz aWwdwyz aNndnzx aSsdszx aTtdtxy aBbdbxy

p方程中将不出现pE，或者说在p的方程中aE将为0(式6.17 6.23a)．于是也就不需要有关pE边界的任何信息了．

aPaEaWaNaSaTaB

N(PP)xyzb[(u*)w(u*)e]yzt

****[(v)s(v)n]zx[(w)b(w)t]xy

(4)．由方程(6．11)，以p加p计算p． p=p*+p’ （6.11）

(5)．利用速度修正公式(6．17)一(6．19)由带星号的速度值 计算u、v和w．

0WSueE已知

图6.10 连续性方程的边界控制容积

6．8 一个修订的算法“SIMPLER”

SIMPLE算法已经广为应用，并且很受用户欢迎．例如，在第九章中将要介绍的所有的流体流动计算均系用这一算法完成的．但是，为了力求改进其收敛速度，已经制订出一个修订的版本．这个版本叫做SIMPLER，它代表修订的SIMPLE的意思[帕坦卡(1979a)]． 6．8—1 问题的提出

ueuede(pp'pE') （6.17） vnvndn(pp'pN') （6.18）

**wtwtdt(pp'pT') （6.19）

(6)．求解那些通过源项、流体物性等等影响流场的其它一些物理量Φ(如温度、浓度以及紊度等物理量)的离散化方程．(如果一个特定的Φ并不影响流场，那么较好的办法是在已经求得了流场的收敛解之后再来计算这个量．)

(7)．把经过修正的压力p处理成一个新的估计的压力p，返回到第二步，重复全部过程，直至求得收敛的解时为止．

6．7—3 压力修正方程的边界条件

动量方程是通用的Φ方程的特殊情况，因此我们的通用边界条件处理方法也同样适用于这些动量方程的边界条件．但是由于p方程并不是一个基本的方程，提一

*aeue'anbunb'b(pp'pE')Ae （6.13）

在推导p方程的过程中所引入的近似(忽略掉项

anbunb)导致过于夸大了的压力修正，因此欠松弛成为迭代过

程中的基本做法．

由于从速度修正公式中消除了相邻点速度修正的影响，压力修正就完全担负起修正速度的重担，从而导致一个相当强烈的压力修正场．在大多数情况下，作这样的假设是合理的，即：可以认为压力修正方程用来修正速度时是相当好的，而用来修正压力时则相当差．

为了评价这一段论述，让我们来考察一个非常简单的例子：一个已知入口边界速度的一维常密度流动问题．显而易见，本问题中的速度只受连续性的控制，因此在第一次迭代终了所得到的满足连续性的速度场自然就是最终的答案了．但是由于p方程的近似特性，所算得的压力却远远偏离最终的解．尽管在迭代过程中早就得到了正确的速度场，可是收敛的压力场却需要经过许多次迭代才能建立．

如果我们只是用压力修正方程来修正速度，而提供某些其他的方法来获得改进的压力场，我们就构成了一个更为有效的算法．这就是SIMPLER的基础． 6．8—2 压力方程

29

传热与流体流动的数值计算 传热与流体流动的数值计算

一个用以求得压力场的方程可由如下方法推得，首先把动量方程写成：

(u)div(uu)div(gradu)BxVx （2.11） tˆeuanbnbub （6.26）

ae3．计算压力方程(6．30)的系数，求解方程(6．30)以得到压力场．

aeueanbunbb(pppE)Ae （6.6）

ueanbnbaPPPaEPEaWPWaNPNaSPSaTPTaBPBb （6.30）

4．把这个压力场当作p看待，求解动量方程以得到u、v及w。

ubaede(pppE) （6.25）

5．计算质量源b[方程(6．23h)]，随后求解p方程(6．22)。

ˆe，即： 其中de的定义已经由方程(6.16)给定．现在我们这样来定义一个假速度uˆeuanbnbubaPPP'aEPE'aWPW'aNPN'aSPS'aTPT'aBPB'b (6．22)

6．采用方程(6．17)—(6．19)修正速度场，但不修正压力场．

ae （6.26）

ˆe系由相邻点的速度unb组成而与压力无关，方程(6．25)现在就变成： 可以注意到：uueuede(pp'pE') vnvndn(pp'pN') wtwtdt(pp'pT')

7．如果有必要的话，求解对其它各个Φ的离散化方程．

8．返回步骤2，并重复上述步骤直到收敛时为止．

6．8—4 讨论

1．可以很容易看出：对于6．8—1小节中所讨论的一维问题，SIMPLER算法可以马上给出一个收敛的解．一般说来，由于压力修正方程产生合理的速度场，而压力方程可以算出一个给定速度场的直接结果(不带有任何一点近似)，这样收敛到最终的解应当要快得多．

2．在SIMPLE算法中，估计的压力场起着重要的作用．与此相反，SIMPLER算法不采用估计的压力，而是由一已知的速度场求出压力场．

3．如果已知的速度场碰巧就是正确的速度场，那么在SIMPLER中的压力方程就将得出正确的压力场，任何进一步的迭代都是没有必要的．相反，如果采用相同的正确的速度场和一个估计的压力场来开始SIMPLE程序，那么，实际上一开始情况就会变得很槽．采用估计的压力将会导致与已知的正确速度不同的带星号速度．于是p方程的近似将会在第一次迭代的末尾产生不正确的速度与压力场．在这种情况下，尽管事实上我们在一开始就有了一个正确的速度场，然而要达到收敛却仍然需要进行许多次的迭代．

4．由于压力方程与压力修正方程极其相似，在6．7—3小节中对p方程的边界条件的有关讨论也适合于压力方程．此外，在6．7—4小节中所讨论的压力的相对特性，或许在提及压力方程时，早已描述了．

5．尽管已经发现SIMPLER程序比SIMPLE程序收敛得快，但是应当承认一次SIMPLER迭代需要较多的计算机工作．首先，除了要解SIMPLE中的所有方程之外，还必须解压力方程，其次，在SIMPLE中并不存在与SIMPLER相对

***ˆede(pppE) （6.27） ueu类似地，我们可以写出：

vnvndn(pppN) （6.28） wtwtdt(pppT) （6.29）

ˆ、vˆ、wˆ出现在原来u、v、w所处的位置上，而很容易看出这些方程与方程(6.17)—(6.19)之间的相似性．这里uˆ、vˆ、wˆ的新的速度—压力关系来完成6.6节中的推导，就可压力p本身取代了p的位置．根据这个情况要是用包含有u以得到一个有关压力的方程．这个方程可以写成：

(u)(v)(w)(u)0 （6.20） txyzaPPPaEPEaWPWaNPNaSPSaTPTaBPBb （6.30）

其中aE、aW、aN、aS、aT、aB以及aP由方程(6．23a)一(6．23g)给顶，b由下式给出：

(P)xyzbP[(u)w(u)e]yz （6.31） t[(v)s(v)n]zx[(w)b(w)t]xy应当注意到：对b的不同表达式是压力方程(6．30)与压力修正方程(6．22)之间的唯一不同．对b的表达式(6．31)

0ˆ、vˆ、wˆ的工作．然而由于SIMPLER只需要较少的迭代次数就可以达到收敛，因而每次迭代所增加的工应的需要计算u作量已远远被总的计算工作量的节省所补偿．

ˆ、vˆ、wˆ，而对p方程的b则用带星号的速度项来计算． 采用假速度u尽管压力方程与压力修正方程几乎是相同的，但是两者之间存在着一个主要的差别：在推导压力方程时，没有作任

何的近似假设．于是，如果用一个正确的速度场来计算假速度，压力方程将马上给出正确的压力． 6．8—3 “SIMPLER”算法

修订的算法是这样的，其中解压力方程以求得压力场，而求解压力修正方程仅仅是为了修正速度．运算的步骤可以叙述如下：

1．由一估计的速度场开始．

ˆ、vˆ、wˆ． 2．计算动量方程的系数，随后由诸如方程(6．26)那样一些方程通过代入相邻点速度unb的值计算u

30