CFD 基 础（流体力学）

第1章 CFD 基 础

计算流体动力学(computational fluid dynamics，CFD)是流体力学的一个分支，它通过计算机模拟获得某种流体在特定条件下的有关信息，实现了用计算机代替试验装置完成“计算试验”，为工程技术人员提供了实际工况模拟仿真的操作平台，已广泛应用于航空航天、

热能动力、土木水利、汽车工程、铁道、船舶工业、化学工程、流体机械、环境工程等 领域。

本章介绍CFD一些重要的基础知识，帮助读者熟悉CFD的基本理论和基本概念，为计算时设置边界条件、对计算结果进行分析与整理提供参考。

1.1 流体力学的基本概念

1.1.1 流体的连续介质模型

流体质点(fluid particle)：几何尺寸同流动空间相比是极小量，又含有大量分子的微 元体。

连续介质(continuum/continuous medium)：质点连续地充满所占空间的流体或固体。 连续介质模型(continuum/continuous medium model)：把流体视为没有间隙地充满它所占据的整个空间的一种连续介质，且其所有的物理量都是空间坐标和时间的连续函数的一种假设模型：u =u(t,x,y,z)。

1.1.2 流体的性质

1. 惯性

惯性(fluid inertia)指流体不受外力作用时，保持其原有运动状态的属性。惯性与质量有关，质量越大，惯性就越大。单位体积流体的质量称为密度(density)，以r表示，单位为kg/m3。对于均质流体，设其体积为V，质量为m，则其密度为

m?? (1-1)

V对于非均质流体，密度随点而异。若取包含某点在内的体积?V，其中质量?m，则该点密度需要用极限方式表示，即

?m (1-2) ??lim?V?0?V2. 压缩性

作用在流体上的压力变化可引起流体的体积变化或密度变化，这一现象称为流体的可压缩性。压缩性(compressibility)可用体积压缩率k来量度

k??dV/Vd?/? (1-3) ?dpdp式中：p为外部压强。

在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为可压缩流体，例如高速流动的气体。若不考虑流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体，如水、油等。

3. 粘性

粘性(viscosity)指在运动的状态下，流体所产生的抵抗剪切变形的性质。粘性大小由粘度来量度。流体的粘度是由流动流体的内聚力和分子的动量交换所引起的。粘度有动力粘度?和运动粘度?之分。动力粘度由牛顿内摩擦定律导出：

du??? (1-4)

dy式中：?为切应力，Pa；?为动力粘度，Pa ?s；du/dy为流体的剪切变形速率。

运动粘度与动力粘度的关系为

??? (1-5) ?式中：?为运动粘度，m2/s。

在研究流体流动过程中，考虑流体的粘性时，称为粘性流动，相应的流体称为粘性流体；当不考虑流体的粘性时，称为理想流体的流动，相应的流体称为理想流体。

根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛顿流体。牛顿流体严格满足牛顿内摩擦定律且?保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分为塑性流体、假塑性流体、胀塑性流体3种。

塑性流体，如牙膏等，它们有一个保持不产生剪切变形的初始应力?0，只有克服了这个初始应力后，其切应力才与速度梯度成正比，即

du???0?? (1-6)

dy假塑性流体，如泥浆等，其切应力与速度梯度的关系是

?du??????，n?1 (1-7)

?dy?n胀塑性流体，如乳化液等，其切应力与速度梯度的关系是

?du??????，n?1 (1-8)

?dy?n1.1.3 流体力学中的力与压强

1. 质量力

与流体微团质量大小有关并且集中在微团质量中心的力称为质量力(body force)。在重力场中有重力mg；直线运动时，有惯性力ma。质量力是一个矢量，一般用单位质量所具

有的质量力来表示，其形式如下：

f?fxi?fyj?fzk (1-9)

式中：fx，fy，fz为单位质量力在各轴上的投影。

2. 表面力大小与表面面积有关而且分布作用在流体表面上的力称为表面力(surface force)。表面力按其作用方向可以分为两种：一是沿表面内法线方向的压力，称为正压力；另一种是沿表面切向的摩擦力，称为切向力。

对于理想流体的流动，流体质点只受到正压力，没有切向力；对于粘性流体的流动，流体质点所受到的作用力既有正压力，也有切向力。

作用在静止流体上的表面力只有沿表面内法线方向的正压力。单位面积上所受到的表面力称为这一点处的静压强。静压强具有两个特征：①静压强的方向垂直指向作用面； ②流场内一点处静压强的大小与方向无关。

3. 表面张力

在液体表面，界面上液体间的相互作用力称为张力。在液体表面有自动收缩的趋势，收缩的液面存在相互作用的与该处液面相切的拉力，称为液体的表面张力(surface tension)。正是这种力的存在，引起弯曲液面内外出现压强差以及常见的毛细现象等。

试验表明，表面张力大小与液面的截线长度L成正比，即

T??L (1-10)

式中：?为表面张力系数，它表示液面上单位长度截线上的表面张力，其大小由物质种类决定，其单位为N/m。

4. 绝对压强、相对压强及真空度

标准大气压的压强是101325Pa(760mm汞柱)，通常用patm表示。若压强大于大气压，则以该压强为计算基准得到的压强称为相对压强(relative pressure)，也称为表压强，通常用pr表示。若压强小于大气压，则压强低于大气压的值就称为真空度(vacuum)，通常用pv表示。如以压强0Pa为计算的基准，则这个压强就称为绝对压强(absolute pressure)，通常用ps表示。这三者的关系如下：

pr?ps?patm (1-11)

pv?patm?ps (1-12)

在流体力学中，压强都用符号p表示，但一般来说有一个约定：对于液体，压强用相对压强；对于气体，特别是马赫数大于0.1的流动，应视为可压缩流，压强用绝对压强。

压强的单位较多，一般用Pa，也可用bar，还可以用汞柱、水柱，这些单位换算如下：

1Pa=1N/m2

1bar=105Pa

1patm=760mmHg=10.33mH2O=101325Pa

5. 静压、动压和总压

对于静止状态下的流体，只有静压强。对于流动状态的流体，有静压强(static pressure)、动压强(dynamic pressure)、测压管压强(manometric tube pressure)和总压强(total pressure)之

分。下面从伯努利(Bernoulli)方程(也有人称其为伯努里方程)中分析它们的意义。

伯努利方程阐述一条流线上流体质点的机械能守恒，对于理想流体的不可压缩流动其表达式如下：

pv2??z?H (1-13) ?g2g式中：p/?g称为压强水头，也是压能项，为静压强；v2/2g称为速度水头，也是动能项；z称为位置水头，也是重力势能项，这三项之和就是流体质点的总的机械能；H称为总的

水头高。

将式(1-13)两边同时乘以?g，则有

1p??v2??gz??gH (1-14)

21式中：p称为静压强，简称静压；?v2称为动压强，简称动压；?gH称为总压强，简称

2总压。对于不考虑重力的流动，总压就是静压和动压之和。

1.1.4 流体运动的描述

1. 流体运动描述的方法

描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉描述。

拉格朗日(Lagrange)描述也称随体描述，它着眼于流体质点，并将流体质点的物理量认为是随流体质点及时间变化的，即把流体质点的物理量表示为拉格朗日坐标及时间的函数。设拉格朗日坐标为(a,b,c)，以此坐标表示的流体质点的物理量，如矢径、速度、压强等等在任一时刻t的值，便可以写为a、b、c及t的函数。

若以f表示流体质点的某一物理量，其拉格朗日描述的数学表达式为

f?f(a,b,c,t) (1-15)

例如，设时刻t流体质点的矢径即t时刻流体质点的位置以r表示，其拉格朗日描述为

r?r(a,b,c,t) (1-16) 同样，质点的速度的拉格朗日描述是

v?v(a,b,c,t) (1-17) 欧拉描述，也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化，即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q1,q2,q3)，用欧拉坐标表示的各空间点上的流体物理量如速度、压强等，在任一时刻t的值，可写为q1、q2、q3及t的函数。从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个场。因此，欧拉描述实际上描述了一个个物理量的场。

若以f表示流体的一个物理量，其欧拉描述的数学表达式是(设空间坐标取用直角坐标)

f?F(x,y,z,t)?F(r,t) (1-18) 如流体速度的欧拉描述是

v?v(x,y,z,t) (1-19)

2. 拉格朗日描述与欧拉描述之间的关系

拉格朗日描述着眼于流体质点，将物理量视为流体坐标与时间的函数；欧拉描述着眼于空间点，将物理量视为空间坐标与时间的函数。它们可以描述同一物理量，必定互相相关。设表达式f?f(a,b,c,t)表示流体质点(a,b,c)在t时刻的物理量；表达式f?F(x,y,z,t)表示空间点(x,y,z)在时刻t的同一物理量。如果流体质点(a,b,c)在t时刻恰好运动到空间点(x,y,z)上，则应有

?x?x(a,b,c,t)??y?y(a,b,c,t) (1-20) ?z?z(a,b,c,t)?F(x,y,z,t)?f(a,b,c,t) (1-21)

事实上，将式(1-16)代入式(1-21)左端，即有

F(x,y,z,t)?F[x(a,b,c,t),y(a,b,c,t),z(a,b,c,t),t] (1-22)

?f(a,b,c,t)或者反解式(1-16)，得到

?a?a(x,y,z,t)??b?b(x,y,z,t) (1-23) ?c?c(x,y,z,t)?将式(1-23)代入式(1-21)的右端，也应有

f(a,b,c,t)?f[a(x,y,z,t),b(x,y,z,t),c(x,y,z,t),t] (1-24)

?F(x,y,z,t)由此，可以通过拉格朗日描述推出欧拉描述，同样也可以由欧拉描述推出拉格朗日 描述。

3. 随体导数

流体质点物理量随时间的变化率称为随体导数(substantial derivative)，或物质导数、质点导数。

按拉格朗日描述，物理量f表示为f?f(a,b,c,t)，f的随体导数就是跟随质点(a,b,c)的物理量f对时间t的导数?f/?t。例如，速度v(a,b,c,t)是矢径r(a,b,c,t)对时间的偏导数，

?r(a,b,c,t) (1-25) v(a,b,c,t)??t即随体导数就是偏导数。

按欧拉描述，物理量f表示为f?F(x,y,z,t)，但?F/?t并不表示随体导数，它只表示物理量在空间点(x,y,z,t)上的时间变化率。而随体导数必须跟随t时刻位于(x,y,z,t)空间点上的那个流体质点，其物理量f的时间变化率。由于该流体质点是运动的，即x、y、z是变的，若以a、b、c表示该流体质点的拉格朗日坐标，则x、y、z将依式(1-16)变化，从而f =F(x,y,z,t)的变化依连锁法则处理。因此，物理量f =F(x,y,z,t)的随体导数是

DF(x,y,z,t)?DF[x(a,b,c,t),y(a,b,c,t),z(a,b,c,t),t]Dt?F?x?F?y?F?z?F?????x?t?y?t?z?t?t (1-26)

?F?F?F?F?u?v?w??x?y?z?t?F?(v??)F??t式中：D/Dt表示随体导数。

从中可以看出，对于质点物理量的随体导数，欧拉描述与拉格朗日描述大不相同。前者是两者之和，而后者是直接的偏导数。

4. 定常流动与非定常流动

根据流体流动过程以及流动过程中的流体的物理参数是否与时间相关，可将流动分为定常流动(steady flow)与非定常流动(unsteady flow)。

定常流动：流体流动过程中各物理量均与时间无关，这种流动称为定常流动。

非定常流动：流体流动过程中某个或某些物理量与时间有关，则这种流动称为非定常流动。

5. 流线与迹线

常用流线和迹线来描述流体的流动。

迹线(track)：随着时间的变化，空间某一点处的流体质点在流动过程中所留下的痕迹称为迹线。在t =0时刻，位于空间坐标(a,b,c)处的流体质点，其迹线方程为

?dx(a,b,c,t)?udt??dy(a,b,c,t)?vdt (1-27) ?dz(a,b,c,t)?wdt?式中：u、v、w分别为流体质点速度的三个分量；x、y、z为在t时刻此流体质点的空间 位置。

流线(streamline)：在同一个时刻，由不同的无数多个流体质点组成的一条曲线，曲线上每一点处的切线与该质点处流体质点的运动方向平行。流场在某一时刻t的流线方程为

dxdydz (1-28) ??u(x,y,z,t)v(x,y,z,t)w(x,y,z,t)对于定常流动，流线的形状不随时间变化，而且流体质点的迹线与流线重合。在实际流场中除驻点或奇点外，流线不能相交，不能突然转折。

6. 流量与净通量

流量(flux)：单位时间内流过某一控制面的流体体积称为该控制面的流量Q，其单位为m3/s。若单位时间内流过的流体是以质量计算，则称为质量流量Qm；不加说明时“流量”一词概指体积流量。在曲面控制面上有

Q???v?ndA (1-29)

A

净通量(net flux)：在流场中取整个封闭曲面作为控制面A，封闭曲面内的空间称为控制体。流体经一部分控制面流入控制体，同时也有流体经另一部分控制面从控制体中流出，此时流出的流体减去流入的流体，所得出的流量称为流过全部封闭控制面A的净流量(或净通量)，通过式(1-30)计算：

q???v?ndA (1-30)

A对于不可压缩流体来说，流过任意封闭控制面的净通量等于0。 7. 有旋流动与有势流动

由速度分解定理，流体质点的运动可以分解为： (1) 随同其他质点的平动； (2) 自身的旋转运动；

(3) 自身的变形运动(拉伸变形和剪切变形)。

在流动过程中，若流体质点自身做无旋运动(irrotational flow)，则称流动是无旋的，也就是有势的，否则就称流动是有旋流动(rotational flow)。流体质点的旋度是一个矢量，通常用?表示，其大小为

ijk1????? (1-31)

2?x?y?zuvw若?=0，则称流动为无旋流动，否则就是有旋流动。

?与流体的流线或迹线形状无关；粘性流动一般为有旋流动；对于无旋流动，伯努利方程适用于流场中任意两点之间；无旋流动也称为有势流动(potential flow)，即存在一个势函数?(x,y,z,t)，满足：

V?grad? (1-32) 即

u?8. 层流与湍流

流体的流动分为层流流动(laminar flow)和湍流流动(turbulent flow)。从试验的角度来看，层流流动就是流体层与层之间相互没有任何干扰，层与层之间既没有质量的传递也没有动量的传递；而湍流流动中层与层之间相互有干扰，而且干扰的力度还会随着流动而加大，层与层之间既有质量的传递又有动量的传递。

判断流动是层流还是湍流，是看其雷诺数是否超过临界雷诺数。雷诺数的定义如下：

VL (1-34) Re??????? (1-33) ，v?，w??x?y?z?式中：V为截面的平均速度；L为特征长度；?为流体的运动粘度。

对于圆形管内流动，特征长度L取圆管的直径d。一般认为临界雷诺数为2320，即

Re?vd? (1-35)

当Re<2320时，管中是层流；当Re>2320时，管中是湍流。

对于异型管道内的流动，特征长度取水力直径dH，则雷诺数的表达式为

VdRe?H (1-36)

?异型管道水力直径的定义如下：

A (1-37) S式中：A为过流断面的面积；S为过流断面上流体与固体接触的周长。

临界雷诺数根据形状的不同而有所差别。根据试验几种异型管道的临界雷诺数如 表1-1所示。

dH?4表1-1 几种异型管道的临界雷诺数

正方形 正三角形 偏心缝隙 管道截面形状 Re?VdH Va 3?V ? Va? ?(D?d) Rec 对于平板的外部绕流，特征长度取沿流动方向的长度，其临界雷诺数为5×105～3×106。

2070 1930 1000 1.2 CFD基本模型

流体流动所遵循的物理定律，是建立流体运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。

1.2.1 基本控制方程

1. 系统与控制体

在流体力学中，系统是指某一确定流体质点集合的总体。系统以外的环境称为外界。分隔系统与外界的界面，称为系统的边界。系统通常是研究的对象，外界则用来区别于系统。系统将随系统内质点一起运动，系统内的质点始终包含在系统内，系统边界的形状和所围空间的大小可随运动而变化。系统与外界无质量交换，但可以有力的相互作用，及能

量(热和功)交换。

控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。

2. 质量守恒方程(连续性方程)

在流场中，流体通过控制面A1流入控制体，同时也会通过另一部分控制面A2流出控制体，在这期间控制体内部的流体质量也会发生变化。按照质量守恒定律，流入的质量与流出的质量之差，应该等于控制体内部流体质量的增量，由此可导出流体流动连续性方程的积分形式为

??dxdydz????v?ndA?0 (1-38) ?t???VA式中：V表示控制体，A表示控制面。等式左边第一项表示控制体V内部质量的增量；第

二项表示通过控制表面流入控制体的净通量。

根据数学中的奥-高公式，在直角坐标系下可将其化为微分形式：

???(?u)?(?v)?(?w)?u?v?w?0 (1-39) ?t?x?y?z对于不可压缩均质流体，密度为常数，则有

?u?v?w???0 (1-40) ?x?y?z对于圆柱坐标系，其形式为

???vr?(?vr)?(?v?)?(?vz)?????0 (1-41) ?tr?rr???z对于不可压缩均质流体，密度为常数，则有

vr?vr?v??vz????0 (1-42) r?rr???z3. 动量守恒方程(运动方程)

动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S方程，其微分形式表达如下：

?du?pxx?pyx?pzx???F????bxdt?x?y?z???pxy?pyy?pzy?dv???F??? (1-43) ?bydt?x?y?z??dw?pyz?pzz?p????Fbz?xz??dt?x?y?z??式中：Fbx、Fby、Fbz分别是单位质量流体上的质量力在三个方向上的分量；pyx是流体内应力张量的分量。

动量守恒方程在实际应用中有许多表达形式，其中比较常见的有如下几种。 (1) 可压缩粘性流体的动量守恒方程

?du?p?????u2??u?v?w??????fx?????2???????????x?x??dt???x3??x?y?z?????????u?v??????w?u??????????????????y???y?x???z???x?z????dv?p?????v2??u?v?w????????fy?????2?????????y?y??dt???y3??x?y?z???? (1-44) ??????v?w??????u?v?????????????????z?z?y?x?y?x?????????????w2??u?v?w??????dw??f??p?????2?????????z?dt?z?z????z3??x?y?z??????????w?u??????v?w???????x??z????z????z??zy???x?????????(2) 常粘性流体的动量守恒方程

dv???F?gradp?grad(divv)???2v (1-45) dt3(3) 常密度常粘性流体的动量守恒方程

dv???F?gradp???2v (1-46) dt(4) 无粘性流体的动量守恒方程(欧拉方程)

dv???F?gradp (1-47) dt(5) 静力学方程

?F?gradp (1-48)

?(6) 相对运动方程

在非惯性参考系中的相对运动方程是研究像大气、海洋及旋转系统中流体运动的所必须考虑的。由理论力学得知，绝对速度va为相对速度vr及牵连速度ve之和，即

va?vr?ve (1-49)

其中，ve?v0???r，v0为运动系中的平动速度，?是其转动角速度，r为质点矢径。

而绝对加速度aa为相对加速度ar、牵连加速度ae及科氏加速度ac之和，即

aa?ar?ae?ac (1-50)

其中，ae?dv0d???r???(??r)，ac?2??vr。 dtdt将绝对加速度代入运动方程，即得到流体的相对运动方程

dv?r??Fb?divP?ac?2?vr (1-51) dt

4. 能量守恒方程

将热力学第一定律应用于流体运动，把式(1-51)各项用有关的流体物理量表示出来，即是能量方程。如式(1-52)所示。

??????T(?E)?[ui(?E?p)]?k?hJ?u(?)??effj?j?jijeff??Sh (1-52) ?t?xi?xi??xi?j?ui2式中：E?h??；keff是有效热传导系数，keff?k?kt，其中kt是湍流热传导系数，

?2根据所使用的湍流模型来定义；Jj?是组分j的扩散流量；Sh包括了化学反应热以及其他用

p户定义的体积热源项；方程右边的前3项分别描述了热传导、组分扩散和粘性耗散带来的能量输运。

1.2.2 湍流模型

湍流是自然界广泛存在的流动现象。大气、海洋环境的流动，飞行器和船舰的绕流，叶轮机械、化学反应器、核反应器中的流体运动都是湍流。湍流流动的核心特征是其在物理上近乎于无穷多的尺度和数学上强烈的非线性，这使得人们无论是通过理论分析、实验研究还是计算机模拟来彻底认识湍流都非常困难。回顾计算流体力学的发展，特别是活跃的20世纪80年代，不仅提出和发展了一大批高精度、高分辨率的计算格式，从主控方程看相当成功地解决了欧拉方程的数值模拟，可以说欧拉方程数值模拟方法的精度已接近于它有效使用范围的极限；同时还发展了一大批有效的网格生成技术及相应的软件，具体实现了工程计算所需要的复杂外形的计算网格；且随着计算机的发展，无论从计算时间还是从计算费用考虑，欧拉方程都已能适用于各种实践所需。在此基础上，20世纪80年代还进行了求解可压缩雷诺平均方程及其三维定态粘流流动的模拟。20世纪90年代又开始一个非定常粘流流场模拟的新局面，这里所说的粘流流场具有高雷诺数、非定常、不稳定、剧烈分离流动的特点，显然需要继续探求更高精度的计算方法和更实用可靠的网格生成技术。但更为重要的关键性的决策将是，研究湍流机理，建立相应的模式，并进行适当的模拟仍是解决湍流问题的重要途径。

1. 湍流模型分类

湍流流动模型很多，但大致可以归纳为以下3类。

第一类是湍流输运系数模型，即将速度脉动的二阶关联量表示成平均速度梯度与湍流粘性系数的乘积，用笛卡儿张量表示为

??u?uj?2??ui?uj???t?i???k?ij (1-53) ???x?3?xi??j模型的任务就是给出计算湍流粘性系数?t的方法。根据建立模型所需要的微分方程的数目，可以分为零方程模型(代数方程模型)、单方程模型和双方程模型。

第二类是抛弃了湍流输运系数的概念，直接建立湍流应力和其他二阶关联量的输运 方程。

第三类是大涡模拟。前两类是以湍流的统计结构为基础，对所有涡旋进行统计平均。

大涡模拟把湍流分成大尺度湍流和小尺度湍流，通过求解三维经过修正的Navier-Stokes方程(纳维-斯托克斯方程，简称N-S方程)，得到大涡旋的运动特性，而对小涡旋运动还采用上述的模型。

实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度要高，应用简单，节省计算时间，同时也具有通用性。

Fluent 提供的湍流模型包括：单方程(Spalart-Allmaras)模型、双方程模型(标准k-?模型、重整化群k-?模型、可实现k-?模型)及雷诺应力模型和大涡模拟，如图1-1所示。

图1-1 湍流模型详解

2. 平均量输运方程

雷诺平均就是把Navier-Stokes方程中的瞬时变量分解成平均量和脉动量两部分。对于速度，有

ui?ui?ui? (1-54) 式中：ui和ui?分别是平均速度和脉动速度(i?1,2,3)。

类似地，对于压力等其他标量，也有

?????? (1-55)

式中：?表示标量，如压力、能量、组分浓度等。

把上面的表达式代入瞬时的连续与动量方程，并取平均(去掉平均速度上的横线)，可以把连续与动量方程写成如下的笛卡儿坐标系下的张量形式：

????(?ui)?0 (1-56) ?t?xi???ui?uj2?ul???????ij???ui?uj????????? (1-57) ???x?x3?x?x???il??j??j?上面两个方程称为雷诺平均的Navier-Stokes(RANS)方程。它们和瞬时Navier-Stokes方

?dui?p???dt?xi?xj

程有相同的形式，只是速度或其他求解变量变成了时间平均量。额外多出来的项??ui?uj?是雷诺应力，表示湍流的影响。

对于密度变化的流动过程，如燃烧问题，需要采用法夫雷(Favre)平均才可以求解。法夫雷平均就是除了压力和密度本身以外，所有变量都用密度加权平均。变量的密度加权平均定义如下：

????/? (1-58) 式中：符号~表示密度加权平均，对应于密度加权平均值的脉动值用???表示，有???????。显然，这种脉动值的简单平均值不为零，但它的密度加权平均值等于零，即????0，?????0。 为了求解方程(1-57)，必须模拟雷诺应力项以使方程封闭。通常的方法是应用Boussinesq假设，认为雷诺应力与平均速度梯度成正比，表达式如下：

??ui?uj?2??u?????uiuj??t?????k??ti??ij (1-59) ???x??xi??j?xi?3?Boussinesq假设被用于单方程模型和k-?双方程模型。这种近似方法好处是与求解湍

流粘性系数有关的计算时间比较少。例如，在Spalart-Allmaras单方程模型中只多求解一个表示湍流粘性的输运方程；在k-?双方程模型中只需多求解湍动能k和耗散率?两个方程，湍流粘性系数用湍动能k和耗散率?的函数来描述。Boussinesq假设的不足之处是假设?t是个各向同性标量，对于一些复杂流动，该条件并不是严格成立，所以具有其应用局限性。

另外的近似方法是求解雷诺应力各分量的输运方程。这也需要额外再求解一个标量方程，通常是耗散率?方程。这就意味着对于二维湍流流动问题，需要多求解4个输运方程，而三维湍流问题需要多求解7个方程，需要较多的计算时间，要求更高的计算机内存。

在很多情况下基于Boussinesq假设的模型很好用，而且计算量并不是很大。但是，如果湍流场各向异性很明显，如强旋流动以及应力取得的二次流等流动中，求解RSM模型可以得到更好的结果。

3. 常用湍流模型简介

1) 单方程(Spalart-Allmaras)模型

单方程模型求解变量是v，表征出了近壁(粘性影响)区域以外的湍流运动粘性系数。v的输运方程为

??v??dv1????v?????Gv??(???v)?C??Yv (1-60) ???b2???dt?v??x?x?x?j???j???j??式中：Gv是湍流粘性产生项；Yv是由于壁面阻挡与粘性阻尼引起的湍流粘性的减少；?v和Cb2是常数；v是分子运动粘性系数。

湍流粘性系数?t??vfv1，其中，fv1是粘性阻尼函数，定义为fv1?而湍流粘性产生项Gv模拟为Gv?Cb1?Sv，其中S?S??3?3?Cv13，??v。v?v，Cb1和kfv2，fv2?1?21??fv1kd2

1??u?u?是常数，d是计算点到壁面的距离；S?2?ij?ij，?ij??j?i?。在Fluent软件中，

?2???xi?xj?考虑到平均应变率对湍流产生也起到很大作用，S?|?ij|?Cprodmin(0,|Sij|?|?ij|)，其中，?u?1??uCprod=2.0，|?ij|?2?ij?ij，|Sij|?2SijSij，平均应变率Sij??j?i?。

?2???xi?xj?在涡量超过应变率的计算区域计算出来的涡旋粘性系数变小。这适合涡流靠近涡旋中心的区域，那里只有“单纯”的旋转，湍流受到抑止。包含应变张量的影响更能体现旋转对湍流的影响。忽略了平均应变，估计的涡旋粘性系数产生项偏高。

6?1?Cw??v?3湍流粘性系数减少项Yv为Yv?Cw1?fw??，其中，fw?g??，

?g?C6?d??w3??6vg?r?Cw2(r6?r)，r?22，Cw1、Cw2、Cw3是常数，在计算r时用到的S受平均应变

Skd率的影响。

上面的模型常数在Fluent软件中默认值为Cb1?0.1335，Cb2?0.622，?v?2/3，

21/6Cv1?7.1，Cw1?Cb1/k2?(1?Cb2)/?v，Cw2?0.3，Cw3?2.0，k?0.41。

2) 标准k-?模型

标准k-?模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推导得到的，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流动为完全湍流，分子粘性的影响可以忽略。因此，标准k-?模型只适合完全湍流的流动过程模拟。标准k-?模型的湍动能k和耗散率?方程为如下形式：

?t??k?dk???????????Gk?Gb????YM (1-61) ?dt?xi???xk?i?????t????d??????2?? (1-62) ??????C1?(Gk?C3?Gb)?C2???dt?xi???xkk?????i?式中：Gk表示由于平均速度梯度引起的湍动能产生，Gb表示由于浮力影响引起的湍动能产生；YM表示可压缩湍流脉动膨胀对总的耗散率的影响。湍流粘性系数?t??C?率ε的湍流普朗特数分别为?k=1.0，??=1.3。

3) 重整化群k-?模型

重整化群k-?模型是对瞬时的Navier-Stokes方程用重整化群的数学方法推导出来的模型。模型中的常数与标准k-?模型不同，而且方程中也出现了新的函数或者项。其湍动能与耗散率方程与标准k-?模型有相似的形式：

dk???k???(??)?keff??Gk?Gb????YM (1-63) dt?xi??xi?d????????2???R (1-64) ?(???eff)??C1?(Gk?C3?Gb)?C2??dt?xi??xi?kk

k2?。

在Fluent中，作为默认值常数，C1?=1.44， C2?=1.92，C3? ??0.09，湍动能k 与耗散

式中：Gk表示由于平均速度梯度引起的湍动能产生，Gb表示由于浮力影响引起的湍动能产生；这些参数与标准k-?模型中相同。YM表示可压缩湍流脉动膨胀对总的耗散率的影响，

?k和??分别是湍动能k 和耗散率ε的有效湍流普朗特数的倒数。湍流粘性系数计算公式

??2k?v为d??1.72dv，其中，v??eff/?，Cv?100。对于前面方程的积分，可?3????v?1?Cv??以精确到有效雷诺数(涡旋尺度)对湍流输运的影响，这有助于处理低雷诺数和近壁流动问

k2题的模拟。对于高雷诺数，上面方程可以给出：?t??C?，C??0.0845。这个结果非常

?有意思，和标准k-?模型的半经验推导给出的常数C??0.09非常近似。在Fluent中，如果是默认设置，用重整化群k-?模型时是针对的高雷诺数流动问题。如果对低雷诺数问题进行数值模拟，必须进行相应的设置。

4) 可实现k-?模型

可实现k-?模型的湍动能及其耗散率输运方程为

?t??k?dk???????????Gk?Gb????YM (1-65) ?dt?xi???x???k?i??td??????????dt?xi??????????2???CS???C?CC3?Gb (1-66) ??121??xkk?v???i????式中：C1?max?0.43,?，??Sk/?。 ??5??在上述方程中，Gk表示由于平均速度梯度引起的湍动能产生，Gb表示由于浮力影响

引起的湍动能产生；YM表示可压缩湍流脉动膨胀对总的耗散率的影响；C2和C1?是常数；

?k和??分别是湍动能及其耗散率的湍流普朗特数。在Fluent中，作为默认值常数，C1?=1.44，C2=1.9，?k=1.0，??=1.2。

该模型的湍流粘性系数与标准k-?模型相同。不同的是，粘性系数中的C?不是常数，

而是通过公式计算得到C??1A0?AsUK*，其中， U*?SijSij?ΩijΩij，Ωij=Ωij?2?ijk?k，

??ij=?ij?2?ijk?k，?ij表示在角速度?k旋转参考系下的平均旋转张量率。模型常数

SijSjS1kki，A0?4.04，As?6cos?，??arccos(6W)，式中W?S?SijSij，S31??uj?ui?。从这些式子中发现，C?是平均应变率与旋度的函数。在平衡边界层惯Sij??????2??xi?xj?性底层，可以得到C??0.09，与标准k-?模型中采用的常数一样。

该模型适合的流动类型比较广泛，包括有旋均匀剪切流、自由流(射流和混合层)、腔道流动和边界层流动。对以上流动过程模拟结果都比标准k-?模型的结果好，特别是可实现k-?模型对圆口射流和平板射流模拟中，能给出较好的射流扩张角。

双方程模型中，无论是标准k-?模型、重整化群k-?模型还是可实现k-?模型，三个模型有类似的形式，即都有k和?的输运方程，它们的区别在于：①计算湍流粘性的方法

不同；②控制湍流扩散的湍流普朗特数不同；③?方程中的产生项和Gk关系不同。但都包含了相同的表示由于平均速度梯度引起的湍动能产生Gk，表示由于浮力影响引起的湍动能产生Gb；表示可压缩湍流脉动膨胀对总的耗散率的影响YM。

湍动能产生项

Gk???ui?uj??uj?xi (1-67)

Gb??gi?t?TPrt?xi (1-68)

式中：Prt是能量的湍流普特朗数，对于可实现k-?模型，默认设置值为0.85；对于重整化群k-?模型，Prt?1/?，??1/Prt?k/?Cp。热膨胀系数???浮力引起的湍动能产生项变为

Gb??gi1?????，对于理想气体，????T?p?t?? (1-69)

?Prt?xi5) 雷诺应力模型

雷诺应力模型(RSM)是求解雷诺应力张量的各个分量的输运方程。具体形式为

???(?uiuj)?(?Ukuiuj)??[?uiujuk?p(?kjui??ikuj)]??t?xk?xk??xk?Uj?????Ui??uu??uu?ujuk??ik????(giuj??gjui?)? (1-70) ij??x?x?xkkk??????ui?uj??u?ujp???2?i?2??k(ujum?ikm?uium?jkm)???x??x?x?xi?kk?j式中：左边的第二项是对流项Cij，右边第一项是湍流扩散项DijT，第二项是分子扩散项DijL，

第三项是应力产生项Pij，第四项是浮力产生项Gij，第五项是压力应变项?ij，第六项是耗散项?ij，第七项系统旋转产生项Fij。

在式(1-69)中，Cij、DijL、Pij、Fij不需要模拟，而DijT、Gij、?ij、?ij需要模拟以封闭方程。下面简单对几个需要模拟项进行模拟。

DijT可以用Delay和Harlow的梯度扩散模型来模拟，但这个模型会导致数值不稳定，在Fluent中是采用标量湍流扩散模型：

DijT???xk??t?uiuj??? (1-71) ??k?xk???式中：湍流粘性系数用?t??C?k2?来计算，根据Lien和Leschziner，?k?0.82，这和标准

k-?模型中选取1.0有所不同。

压力应变项?ij可以分解为三项，即

?ij??ij,1??ij,2??ijw (1-72)

式中：?ij,1、?ij,2和?ijw分别是慢速项、快速项和壁面反射项，具体表述可以参见文献[2]。

浮力引起的产生项Gij模拟为

?T?T? (1-73) Gij???gi?gj??Prt??x?xji??耗散张量?ij模拟为

?t?23式中：YM?2??Mt2，Mt是马赫数；标量耗散率?用标准k-?模型中采用的耗散率输运方

?ij??ij(???YM) (1-74)

程求解。

6) 大涡模拟

湍流中包含了不同时间与长度尺度的涡旋。最大长度尺度通常为平均流动的特征长度尺度。最小尺度为Komogrov尺度。LES的基本假设是：①动量、能量、质量及其他标量主要由大涡输运；②流动的几何和边界条件决定了大涡的特性，而流动特性主要在大涡中体现；③小尺度涡旋受几何和边界条件影响较小，并且各向同性，大涡模拟(LES)过程中，直接求解大涡，小尺度涡旋模拟，从而使得网格要求比DNS低。

LES的控制方程是对Navier-Stokes方程在波数空间或者物理空间进行过滤得到的。过滤的过程是去掉比过滤宽度或者给定物理宽度小的涡旋，从而得到大涡旋的控制方程：

??ui???u?0 (1-75) ?t?xi?ui?p??ij???(?ui)?(?uiuj)?(?)?? (1-76) ?t?xj?xj?xj?xj?xj式中：?ij为亚网格应力，?ij??uiuj??ui?uj。

很明显，上述方程与雷诺平均方程很相似，只不过大涡模拟中的变量是过滤过的量，

而非时间平均量，并且湍流应力也不同。

1.2.3 初始条件和边界条件

计算流体动力学(CFD)分析中，初始条件和边界条件的正确设置是关键的一步。现有的CFD软件都提供了现成的各种类型的边界条件，这里对有关的初始条件和边界条件作一般讨论。

1. 初始条件

顾名思义，初始条件就是计算初始给定的参数，即t?t0时给出各未知量的函数分布，如 ?u?u(x,y,z,t0)?u0(x,y,z)??v?v(x,y,z,t0)?v0(x,y,z)??w?w(x,y,z,t0)?w0(x,y,z) (1-77) ?p?p(x,y,z,t)?p(x,y,z)00?????(x,y,z,t0)??0(x,y,z)???T?T(x,y,z,t0)?T0(x,y,z)很明显，当流体运动定常时，无初始条件问题。

2. 边界条件

所谓边界条件就是流体力学方程组在求解域的边界上，流体物理量应满足的条件。例如，流体被固壁所限，流体将不应有穿过固壁的速度分量；在水面这个边界上，大气压强认为是常数(一般在距离不大的范围内可如此)；在流体与外界无热传导的边界上，流体与边界之间无温差，如此等。由于各种具体问题不同，边界条件提法千差万别，一般要保持恰当：①保持在物理上是正确的；②要在数学上不多不少，刚好能用来确定积分微分方程中的积分常数，而不是矛盾的或有随意性。

通常流体边界分为流固交界面和流流(液液、液气)交界面，下面分别讨论。 1) 流固分界面边界条件

飞机、船舶在空气及水中运动时的流固分界面，水在岸边及底部的流固分界面，均属这一类。一般而言，流体在固体边界上的速度依流体有无粘性而定。对于粘性流体，流体将粘附于固体表面(无滑移)，即

v|F?v|S (1-78) 式中：v|F是流体速度；v|S是固壁面相应点的速度。式(1-78)表明，在流固边界面上，流体在一点的速度等于固体在该点的速度。对于无粘性流体，流体可沿界面滑移，即有速度的切向分量，但不能离开界面，也就是流体的法向速度分量等于固体的法向速度分量，即

vn|F?v|S (1-79)

另外，也可视所给条件，给出无温差条件：

T|F?T|S (1-80) 式中：T|F是流体温度，T|S是固壁面相应点的温度。

2) 液液分界面边界条件

密度不同的两种液体的分界面就属于这一类。一般而言，对分界面两侧的液体情况经常给出的条件是

v1?v2,T1?T2,p1?p2 (1-81)

对应力及传导热情况给出的条件是

???1?u?u|1??2|2 (1-82) ?n?n?T?TQ?k1|1?k2|2 (1-83)

?n?n3) 液气分界面边界条件

液气分界面最典型的是水与大气的分界面，即自由面。由于自由面本身是运动和变形的，而且其形状常常也是一个需要求解的未知函数，因此就有一个自由面的运动学条件问题。设自由面方程为

F(x,y,z,t)?0 (1-84) 并假定在自由面上的流体质点始终保持在自由面上，则流体质点在自由面上一点的法向速度，应该等于自由面本身在这一点的法向速度。经过一系列推导(参见文献[2])，得到自由液面运动学条件：

?F?v??F?0 (1-85) ?t

如果要考虑液气边界上的表面张力，则在界面两侧，两种介质的压强差与表面张力有如下关系：

?11?p2?p1????? (1-86)

?R1R2?这就是自由面上的动力学条件。当不考虑表面张力时，有

p?pa (1-87) 式中：pa为大气压强。

4) 无限远的条件

流体力学中的很多问题，流体域是无限远的。例如，飞机在空中飞行时，流体是无界的。如果将坐标系取在运动物体上，这时无限远处的边界条件为

当x→?时，

u?u?，p?p? (1-88) 其中下标?表示无穷远处的值。

1.3 CFD模型的离散——有限体积法

1.3.1 CFD模型的数值求解方法概述

从上面的分析看到，CFD模型(控制方程)是一系列偏微分方程组，要得到解析解比较困难，目前，均采用数值方法得到其满足实际需要的近似解。

数值方法求解CFD模型的基本思想是：把原来在空间与时间坐标中连续的物理量的场(如速度场、温度场、浓度场等)，用一系列有限个离散点(称为节点，node)上的值的集合来代替，通过一定的原则建立起这些离散点上变量值之间关系的代数方程(称为离散方程，discretization equation)，求解所建立起来的代数方程以获得所求解变量的近似解。在过去的几十年内已经发展了多种数值解法，其间的主要区别在于区域的离散方式、方程的离散方式及代数方程求解的方法这三个环节上。在CFD求解计算中用得较多的数值方法有：有限差分法(finite difference method，FDM)、有限体积法(finite volume method，FVM)、有限元法(finite element method，FEM)及有限分析法(finite analytic method，FAM)。下面简要介绍，后面将着重介绍有限体积法。

1. 有限差分法

有限差分法是历史上采用最早的数值方法，对简单几何形状中的流动与换热问题也是一种最容易实施的数值方法。其基本点是：将求解区域用与坐标轴平行的一系列网格线的交点所组成的点的集合来代替，在每个节点上，将控制方程中每一个导数用相应的差分表达式来代替，从而在每个节点上形成一个代数方程，每个方程中包括了本节点及其附近一些节点上的未知值，求解这些代数方程就获得了所需的数值解。由于各阶导数的差分表达式可以从Taylor(泰勒)展开式来导出，这种方法又称建立离散方程的Taylor展开法。

有限差分法软件一般研究者自己编写，很少看到商品的有限差分法软件。

2. 有限体积法

在有限体积法中将所计算的区域划分成一系列控制体积，每个控制体积都有一个节点作代表，通过将守恒型的控制方程对控制体积作积分来导出离散方程。在导出过程中，需要对界面上的被求函数本身及其一阶导数的构成作出假定，这种构成的方式就是有限体积法中的离散格式。用有限体积法导出的离散方程可以保证具有守恒特性，而且离散方程系数的物理意义明确，是目前流动与传热问题的数值计算中应用最广泛的一种方法。

Phoenics是最早投入市场的有限体积法软件，Fluent、STAR-CD和CFX都是常用的有限体积法软件，它们在流动、传热传质、燃烧和辐射等方面应用广泛。

3. 有限元法

在有限元法中把计算区域划分成一系列单元体(在二维情况下，单元体多为三角形或四边形)，在每个单元体上取数个点作为节点，然后通过对控制方程做积分来获得离散方程。它与有限体积法区别主要在于如下两点。

(1) 要选定一个形状函数(最简单的是线性函数)，并通过单元体中节点上的被求变量之值来表示该形状函数，在积分之前将该形状函数代入到控制方程中去。这一形状函数在建立离散方程及求解后结果的处理上都要应用。

(2) 控制方程在积分之前要乘上一个权函数，要求在整个计算区域上控制方程余量(即代入形状函数后使控制方程等号两端不相等的差值)的加权平均值等于零，从而得出一组关于节点上的被求变量的代数方程组。

有限元法的最大优点是对不规则区域的适应性好。但计算的工作量一般较有限体积法大，而且在求解流动与换热问题时，对流项的离散处理方法及不可压流体原始变量法求解方面没有有限体积法成熟。

Ansys、Sysweld和北京飞箭公司的FEPG(finite element programs generator)等有限元软件比较流行。

4. 有限分析法

有限分析法是由美籍华裔科学家陈景仁教授在1981年提出的。在这种方法中，也像有限差分法那样，用一系列网格线将区域离散，所不同的是每个节点与相邻的4个网格(二维)问题组成计算单元，即一个计算单元由一个中心节点与8个邻点组成。在计算单元中把控制方程中的非线性项(如Navier-Stokes方程中的对流项)局部线性化(即认为流速已知)，并对该单元上未知函数的变化型线作出假设，把所选定型线表达式中的系数和常数项用单元边界节点上未知的变量值来表示，这样该单元内的被求问题就转化为第一类边界条件下的一个定解问题，可以找出其分析解；然后利用这一分析解，得出该单元中点及边界上8个邻点上未知值间的代数方程，此即为单元中点的离散方程。有限分析法中的系数不像有限体积法中那样有明确的物理意义，对不规则区域的适应性也较差。

1.3.2 有限体积法

从上面的简介看到，有限体积法是一种分块近似的计算方法，其中比较重要步骤是计算区域的离散和控制方程的离散。

1. 计算区域的离散化

所谓区域的离散化(domain discretization)实质上就是用一组有限个离散的点来代替原来的连续空间。一般的实施过程是：把所计算的区域划分成许多个互不重叠的子区域(sub-domain)，确定每个子区域中的节点位置及该节点所代表的控制体积。区域离散后，得到以下四种几何要素。

? 节点(node)：需要求解的未知物理量的几何位置。 ? 控制体积(control volume)：应用控制方程或守恒定律的最小几何单位。

? ?

界面(face)：它定义了与各节点相对应的控制体积的界面位置。

网格线(grid line)：连接相邻两节点面形成的曲线簇。

一般把节点看成是控制体积的代表。在离散过程中，将一个控制体积上的物理量定义并存储在该节点处。图1-2给出了一维问题的有限体积法计算网格，图1-3给出了二维问题的有限体积法计算网格。

图1-2 一维的有限体积法网格

图1-3 二维的有限体积法网格

计算区域离散的网格有两类：结构化网格和非结构化网格。节点排列有序，即当给出了一个节点的编号后，立即可以得出其相邻节点的编号，所有内部节点周围的网格数目相同。这种网格称为结构化网格(structured grid)。结构化网格具有实现容易、生成速度快、网格质量好、数据结构简单的优点，但不能实现复杂边界区域的离散。

而非结构化网格的内部节点以一种不规则的方式布置在流场中，各节点周围的网格数目不尽相同。这种网格虽然生成过程比较复杂，但却有极大的适应性，对复杂边界的流场

计算问题特别有效。

2. 控制方程的离散化

前面给出的流体流动问题的控制方程，无论是连续性方程、动量方程，还是能量方程，都可写成如式(1-89)所示的通用形式，

?(?u?)?div(?u?)?div(?grad?)?S (1-89) ?t对于一维稳态问题，其控制方程如式(1-90)所示：

d(?u?)d?d???????S (1-90) dxdx?dx?式中：从左到右各项分别为对流项、扩散项和源项。方程中的?是广义变量，可以为速度、温度或浓度等一些待求的物理量。?是相应于?的广义扩散系数，S是广义源项。变量?在端点A和B的边界值为已知。

有限体积法的关键一步是在控制体积上积分控制方程，在控制体积节点上产生离散的方程。对一维模型方程(1-90)，在图1-2所示的控制体积P上作积分，有

d(?u?)d?d??dV???dV???VSdV (1-91) ??Vdx??Vdx??dx?式中：?V是控制体积的体积值。当控制体积很微小时，?V可以表示为?V?A，这里A是控制体积界面的面积。从而有

d???d???(?u?A)e?(?u?A)w???A????A??S??V (1-92)

dx?e?dx?w?从式(1-92)看到，对流项和扩散项均已转化为控制体积界面上的值。有限体积法最显著的特

点之一就是离散方程中具有明确的物理插值，即界面的物理量要通过插值的方式由节点的物理量来表示。

u、为了建立所需要形式的离散方程，需要找出如何表示式(1-92)中界面e和w处的?、

d?d??、?和。在有限体积法中规定，?、u、?、?和等物理量均是在节点处定义和

dxdx计算的。因此，为了计算界面上的这些物理参数(包括其导数)，需要一个物理参数在节点间的近似分布。可以想象，线性近似是可以用来计算界面物性值的最直接，也是最简单的方式。这种分布叫做中心差分。如果网格是均匀的，则单个物理参数(以扩散系数?为例)的线性插值结果是

?P??E?????e2 (1-93) ????P???Ww??2(?u?A)的线性插值结果是

?P??E?(?u?A)?(?u)Aeee??2 (1-94) ????P?(?u?A)?(?u)AWwww??2

与梯度项相关的扩散通量的线性插值结果是

????E??P?d???A??A???ee??dx?e???(?x)e? (1-95) ???P??W?d?????A??A?ww????dx?w?(?x)w???对于源项S，它通常是时间和物理量?的函数。为了简化处理，将S转化为如下线性 方式：

S?SC?SP?P (1-96)

式中：SC是常数，SP是随时间和物理量?变化的项。将式(1-93)～式(1-96)代入方程(1-92)，有

(?u)eAe?P??E22????P???P??W??eAe?E??A?ww??(?x)e??(?x)w?(?u)wAw?W??P???(SC?SP?P)??V? (1-97)

整理后得

??e??wA?A?S??V???PewP(?x)(?x)ew??????(?u)w?(?u)e???wAw?Aw??W??eAe?Ae??E?SC??V(?x)2(?x)2we????

记为

aP?P?aW?W?aE?E?b (1-98)

式中：

?w(?u)w?a?A?Aw?W(?x)w2w???e(?u)e?aE?Ae?Ae(?x)e2??? (1-99) ?e?w?a?A?A?S??VP?P(?x)e(?x)wew??(?u)e(?u)wAe?Aw?SP??V??aE?aW?22??b?S??VC??对于一维问题，控制体积界面e和w处的面积Ae和Aw均为1，即单位面积。这样?V??x，式(1-99)中各系数可转化为

?w(?u)w?a??W?(?x)w2???e(?u)e??aE? (1-100) (?x)e2???aP?aE?aW?(?u)e?(?u)w?SP??x?22?b?S??xC?方程(1-98)即为方程(1-90)的离散形式，每个节点上都可建立此离散方程，通过求解方

程组，就可得到各物理量在各节点处的值。

为了后续讨论的方便，定义两个新的物理量F和D，其中F表示通过界面上单位面积的对流质量通量(convective mass flux)，简称对流质量流量，D表示界面的扩散传导性(diffusion conductance)。定义表达式如下：

?F??u?? (1-101) ?D???x?这样，F和D在控制界面上的值分别为

?Fw?(?u)w,Fe?(?u)e??w?e (1-102) ?D?,D?e?w(?x)(?x)ew?在此基础上，定义一维单元的Peclet数Pe如下： F?u (1-103) Pe??D?/?x式中：Pe表示对流与扩散的强度之比。当Pe数为0时，对流-扩散演变为纯扩散问题，即流场中没有流动，只有扩散；当Pe>0时，流体沿x方向流动，当Pe数很大时，对流-扩散问题演变为纯对流问题。一般在中心差分格式中，有Pe<2的要求。

将式(1-101)代入方程(1-100)，有

Fw?a?D?Ww?2??a?D?Fe?Ee (1-104) 2??FF?aP?aE?aW?e?w?SP??x22???b?SC??x对于瞬态问题，与稳态问题相似，主要是瞬态项的离散。其一维瞬态问题的通用控制

方程如下：

?(??)?(?u?)???????????S (1-105) ?t?x?x??x?该方程是一个包含瞬态及源项的对流-扩散方程。从左到右，方程中的各项分别是瞬态

项、对流项、扩散项及源项。方程中的?是广义变量，如速度分量、温度、浓度等，?为相应于?的广义扩散系数，S为广义源项。

对于瞬态问题用有限体积法求解时，在将控制方程对控制体积作空间积分的同时，还必须对时间间隔?t作时间积分。对控制体积所作的空间积分与稳态问题相同，这里仅叙述对时间的积分。

将方程(1-105)在一维计算网格上对时间及控制体积进行积分，有

t??tt??t?(??)?(?u?)dVdt??t??V?t?t??V?xdVdt (1-106)

t??tt??t???????t????dVdt??t?SdVdt?x?x???V?V改写后，有

t??t?t??t?(??)?dtdV?[(?u?A)e?(?u?A)w]dt??V????tt?t????tt??tt??t??d???d??????A????A??dt??tS?Vdtdx?e?dx?w??? (1-107)

式中：A是图1-2中控制体积P的界面处的面积。

在处理瞬态项时，假定物理量?在整个控制体积P上均具有节点处值?P，并用线性插

0)/?t来表示??/?t。源项也分解为线性方式S?Sc?SP?P。对流项和扩散项的值值(?P??P按中心差分格式通过节点处的值来表示，则有

t??t????P????E0?(?P??P)?V??t?(?u)eAeP?(?u)wAwW?dt22?? (1-108)

??t??t?t??t??E??P???P??W????t??eAe????wAw???dt??t(Sc?SP?P)?Vdt(?x)(?x)??e?W?????假定变量?P对时间的积分为

?t??tt?Pdt?[f?P?(1?f)?P0]?t (1-109)

式中：上标0代表t时刻；?P是时刻的值；f为0与1之间的加权因子，当f=0时，变量取老值进行时间积分，当f=1时，变量采用新值进行时间积分。

将?P、?E、?W及SC?SP?P采用类似式(1-109)进行时间积分，式(1-108)可写成

???P????E?V??f?(?u)eAeP?(?u)wAwW??t22???0000???W??P?P??E(1?f)?(?u)eAe?(?u)wAw?22??0?(?P??P)?????P???P??W??f??eAe?E??A?ww???(?x)e??(?x)W?????????? (1-110)

0000???E???P????W??P??(1?f)??eAe???A????ww???(?x)e??(?x)W????0[f(SC?SP?P)?(1?f)?(Sc?SP?P)]?V整理后得

????eAe?wAw??(?u)eAe(?u)wAw???V???f??f??fS?V???P??P?2??t2(?x)(?x)????eW?????(?u)wAw?wAw?0????[f?W?(1?f)?W]?2(?x)W????eAe(?u)eAe?0???[f?E?(1?f)?E]?2??(?x)e???eAe(?u)eAe???V??(1?f)??????t(?x)2?e????0?(?u)wAw(?u)wAw?(1?f)???(1?f)S?V??P?SC?VP?22??? (1-111)

同样引入稳态中关于符号F和D的定义，并将原来的定义作一定扩展，即乘以面积

A，有

?Fw?(?u)wAw,Fe?(?u)eAe? (1-112) ?wAw?eAe?D?,D?e?w(?x)(?x)ew?将式(1-112)代入方程(1-111)，得

??V???t?f(De?Dw)????FF?f?e?w??fSP?V??P2??2?F??F??00??w?Dw?[f?W?(1?f)?W]??De?e?[f?E?(1?f)?E]? 2??2????VF?F??0???(1?f)?De?e??(1?f)?Dw?w??(1?f)SP?V??P????t22??????(1-113)

SC?V同样也像稳态问题引入aP、aW和aE，式(1-113)变为

00aP?P?aW[f?W?(1?f)?W]?aE[f?E?(1?f)?E]???VF?F??0???(1?f)?De?e??(1?f)?Dw?w??(1?f)SP?V??P? (1-114) ???t22??????SC?V式中：

?aP?aP0?f(aE?aW)?f(Fe?Fw)?fSP?V??a?D?Fww?W2? ?Fea?D??Ee2??V?0a??P??t?(1-115)

根据f的取值，瞬态问题对时间的积分有几种方案，当f=0时，变量的初值出现在方程(1-114)的右端，从而可直接求出在现时刻的未知变量值，这种方案称为显式时间积分方案。当0

进一步将一维问题扩展为二维与三维问题。在二维问题中，计算区域离散如图1-4所示。发现只是增加第二坐标y，控制体积增加的上下界面，分别用n(north)和s(south)表示，相应的两个邻点记为N和S。在全隐式时间积分方案下的二维瞬态对流-扩散问题的离散方程为

aP?P?aW?W?aE?E?aN?N?aS?S?b (1-116)

式中：

0?aP?ap?aE?aW?aS?aN?(Fe?Fw)?(Fn?Fs)?SP?V??aW?Dw?max(0,Fw)??aE?De?max(0,?Fe)??aS?Ds?max(0,Fs) (1-117) ??aN?Dn?max(0,?Fn)??00?Va??P?P?t?00??b?SC?V?aP??P在三维问题中，计算区域离散如图1-4所示(两个方向的投影)。在二维的基础上增加第

三坐标z，控制体积增加的前后界面，分别用t(top)和b(bottom)表示，相应的两个邻点记为T和B。在全隐式时间积分方案下的三维瞬态对流-扩散问题的离散方程为

aP?P?aW?W?aE?E?aN?N?aS?S?aT?T?aB?B?b (1-118)

式中：

?aP?aP0?aE?aW?aS?aN?(Fe?Fw)?(Fn?Fs)??(Ft?Fb)?SP?V??aW?Dw?max(0,Fw)??aE?De?max(0,?Fe)?a?D?max(0,F)Sss?? (1-119) ?aN?Dn?max(0,?Fn)??aT?Dt?max(0,?Ft)?aB?Db?max(0,Fb)??a0??0?VP?P?t?00??b?SC?V?aP??P

控制体积边界控制体积边界节点节点NnWwsSWwTtyddy ee ee bBdx dxdx dx dzzPEPE控制体积控制体积 图1-4

1.3.3 有限体积法中常用的离散格式

1. 常用的离散格式

有限体积法常用的离散格式有：中心差分格式、一阶迎风格式、混合格式、指数格式、乘方格式、二阶迎风格式、QUICK格式。各种离散格式对一维、稳态、无源项的对流-扩

散问题的通用控制方程(1-120)均能得到式(1-121)的形式。对于高阶情况如式(1-122)所示。

d(?u?)d?d?????? (1-120) dxdx?dx?aP?P?aW?W?aE?E (1-121)

aP?P?aW?W?aWW?WW?aE?E?aEE?EE (1-122)

式中，对于一阶情况，aP?aW?aE?(Fe?Fw)，对于二阶情况，aP?aW?aE?aWW?aEE? (Fe?Fw)，其中系数aW和aE取决于所使用的离散格式(高阶还有aWW和aEE)。为了便于比

较和编程计算，将各种离散格式的系数总结于表1-2中。

表1-2 不同离散格式下系数aW和aE的计算公式

离散格式 中心差分格式 一阶迎风格式 混合格式 Dw?Fw 2系数aW De?Fe 2系数aE Dw?max(Fw,0) De?max(0,?Fe) ??F??max?Fw,?Dw?w?,0? 2????Fwexp(Fw/Dw) exp(Fw/Dw)?1?F???max??Fe,?De?e?,0? 2????Fe exp(Fe/De)?1指数格式

续表 离散格式 乘方格式 系数aW 系数aE Dwmax[0,(1?0.1|Pe|)3]?max[Fw,0]31Dw??Fw??Fe 221aWW???Fw 2Demax[0,(1?0.1|Pe|)3]?max[?Fe,0] 二阶迎风格式 31De?(1??)Fe?(1??)Fw 22aEE?1(1??)Fe 2QUICK 格式 613Dw??wFw??wFe?(1??w)Fw 8881aWW???wFw 8361De??eFe?(1??e)Fe?(1??e)Fw 8881aEE?(1??e)Fe 8 2. 建立离散方程应遵循的原则

在利用有限体积法建立离散方程时，必须遵守以下四条基本原则。 1) 控制体积界面上的连续性原则

当一个面为相邻的两个控制体积所共有时，在这两个控制体积的离散方程中，通过该界面的通量(包括热通量、质量通量、动量通量)的表达式必须相同。即通过某特定界面从一个控制体积所流出的热通量，必须等于通过该界面进入相邻控制体积的热通量，否则，总体平衡就得不到满足。

2) 正系数原则

中心节点系数aP和相邻节点系数anb必须恒为正值。该原则是求得合理解的重要保证。当违背这一原则，结果往往是得到物理上不真实的解。例如，如果相邻节点的系数为负值，就可能出现边界温度的增加引起的相邻节点温度降低。

3) 源项的负斜率线性化原则

源项斜率为负可以保证正系数原则。从式(1-100)中看到，当相邻节点的系数皆为正值，但有源项Sp的存在，中心节点系数aP仍有可能为负。当我们规定Sp≤0，便可以保证aP为正值。

4) 系数aP等于相邻节点系数之和原则

当源项为0时，我们发现中心节点系数等于相邻节点系数之和，而当有源项存在时也应该保证这一原则，如果不能满足这个条件，可以取SP为0。

1.4 流场数值计算算法分析

上面建立了控制方程的离散方程，即建立了可以数值计算的代数方程组。常规解法只能应付已知速度场求温度场分布这类简单的问题，所以对所生成的离散方程进行某种调整，

并且对各未知量(速度、压力、温度等)的求解顺序及方式进行特殊处理。为此，流场数值计算是针对常规解法存在的主要问题进行改善形成的一系列方法集。流场计算的基本过程是在空间上用有限体积法(或其他类似方法)将计算区域离散成许多小的体积单元，在每个体积单元上对离散后的控制方程组进行求解。其本质是对离散方程进行求解，一般可以分为分离解法(segregated method)和耦合解法(coupled method)两大类，各自又根据实际情况扩展成具体的计算方法，如图1-5所示。

图1-5 流场数值计算方法的分类

1. 分离解法

分离解法不直接求解联立方程组，而是顺序地、逐个地求解各变量代数方程组。分离解法中应用广泛的是压力修正法，其求解基本过程如下。

(1) 假定初始压力场。

(2) 利用压力场求解动量方程，得到速度场。

(3) 利用速度场求解连续方程，使压力场得到修正。 (4) 根据需要，求解湍流方程及其他标量方程。 (5) 判断当前时间步上的计算是否收敛。若不收敛，返回第二步，迭代计算；若收敛，重复上述步骤，计算下一时间步的物理量。

压力修正法有很多实现方式，其中，压力耦合方程组的半隐式方法(SIMPLE算法)应用最为广泛，也是各种商用CFD软件普遍采纳的算法。

2. 耦合解法

耦合解法同时求解离散方程组，联立求解出各变量(等)，其求解过程如下。 (1) 假定初始压力和速度等变量，确定离散方程的系数及常数项等。

(2) 联立求解连续方程、动量方程、能量方程。 (3) 求解湍流方程及其他标量方程。

(4) 判断当前时间步上的计算是否收敛。若不收敛，返回第二步，迭代计算；若收敛，重复上述步骤，计算下一时间步的物理量。

1.4.1 SIMPLE算法详解

1. 交错网格

使用交错网格，主要是为了解决在普通网格上离散控制方程时不能检测有问题的压力场的问题，同时，交错网格也是SIMPLE算法实现的基础。

所谓交错网格(staggered grid)，就是将标量(如压力p、温度T和密度等)在正常的网格节点上存储和计算，而将速度的各分量分别在错位后的网格上存储和计算，错位后网格的中心位于原控制体积的界面上。这样对于二维问题，就有三套不同的网格系统，分别用于存储p、u和v；而对于三维问题，就有四套网格系统，分别用于存储p、u、v、w。

二维流动计算的交错网格系统如图1-6所示，主控制体积为求解压力p的控制体积，称为标量控制体积(scale control volume)或p控制体积，控制体积的节点P称为主节点或标量节点(scale node)(如图1-6(a)所示)。速度u在主控制体积的东、西界面e和w上定义和存储，速度v在主控制体积的南、北界面s和n上定义和存储。u和v各自的控制体积则是分别以速度所在位置(界面e和界面n)为中心的，分别称为u控制体积(u-control volume)和v控制体积(v-control volume)，如图1-6(b)和(c)所示。可以看到，u控制体积和v控制体积是与主控制体积不一致的，u控制体积与主控制体积在x方向上有半个网格步长的错位，而v控制体积与主控制体积则在y方向上有半个步长的错位。

图1-6 控制体积

图1-7中所示的均匀网格是向后错位的，因为u的速度uI,J的i位置到标量节点(I,J)的是-1/2 ?xu；同样，v速度vI,J的j位置到标量节点(I,J )的距离是-1/2?yv。当然，也可以选用向后两个错位的速度网格。

在使用了上述交错网格后，生成离散方程的方法和过程与原来的基于普通网格的方法和过程是一样的，只是需要注意所使用的控制体积有所变化。在交错网格中，由于所有标量(如压力、温度、密度等)仍然在主控制体积上存储，因此，以这些标量为因变量的输运

方程的离散过程及离散结果与前面的一样，主要是在交错网格中生成的u和v两个动量方程的离散方程时，积分用的控制体积不再是原来的主控制体积，而是u和v各自的控制体积，同时压力梯度项从源项中分离出来。例如，对u控制体积，该项积分为yj?1xi??p???dxdy?(pI?1,J?pI,J)Ai,J。从而，对于方程(1-98)的离散方程，对u方向的动?yj?xi?1???x?量方程使用u控制体积，可以写出在位置(i,J)处的关于速度ui，J的动量方程的离散形式：

图1-7 交错网格

ai,Jui,J??anbunb?(pI?1,J?pI,J)Ai,J?bi,J (1-123)

式中：Ai,J是u控制体积的东界面或西界面的面积，在二维问题中实际上是?y，即

Ai,J=?y?yj?1?yj (1-124)

式(1-123)中的b为u动量方程的源项部分(不包括压力在内)。对于稳态问题，有

bi,J＝SuC?Vu (1-125) 式(1-125)中的SuC是对源项Su线性化分解的结果，若Su不随速度u而变化，则有SuC≥Su，SuP=0。式(1-125)中的?Vu是u控制体积的体积。式(1-123)中的压力梯度项已经按线性插值的方式进行了离散，线性插值时使用了u控制体积边界上的两个节点的压力差。

在求和记号?anbunb中所包含的E、W、N和S四个邻点是(i-1,J)，(i+1,J)，(i,J+1)和(i,J+1)，它们的位置及主速度在图1-8中标出。图中阴影部分是u控制体积。图1-8中的u控制体积与图1-7中是一致的，这可从节点的编号看出。但是，图1-8中u控制体积的中心也用P来标记，其界面点也用e、w、n和s来标记，这里的标记与图1-7中的同名标记及系数a是由式(1-98)给定的。即

ai,J＝?anb??F?SuP?Vu (1-126)

图1-8 u控制体积及其邻点的速度分量

系数anb取决于所采用的离散格式，在计算式中含有u控制体积界面上的对流质量流量F与扩散传导性D，在采用新编号系统下的计算公式为

F?Fi?1,J1??I,J??I?1,J???I?2,J?Fw?(?u)w?i,J??ui,J?I?1,Jui?1,J? (1-127a)

22?22???I,J???I?1,J?1????I?1,Jui?1,J?I,Jui,J? (1-127b) 22?22?FI,j?FI?1,j1??I,J??I,J?1???I?1,J?1?Fs?(?u)s???vI,j?I?1,JvI?1,j? (1-127c)

22?22?FI,j?1?FI?1,j?11??I,J?1??I,J???I?1,J?Fn?(?u)n???vI,j?1?I?1,J?1vI?1,j?1? (1-127d)

22?22?Fe?(?u)e?Fi?1,J?Fi,JDw?De??I?1,Jxi?xi?1 (1-127e) (1-127f)

(1-127g) (1-127h)

?I,Jxi?1?xiDs?Dn??I?1,J??I,J??I?1,J?1??I,J?14(yJ?yJ?1)?I?1,J?1??I,J?1??I?1,J??I,J4(yJ?1?yJ)采用交错网格对动量方程离散时，涉及不同类别的控制体积，不同的物理量分别在各

自相应的控制体积的节点上定义和存储。例如，密度是在标量控制体积的节点上存储的，如图1-8中的标量节点(I,J)；而速度分量却是在错位后的速度控制体积的节点上存储的，如图1-8中的速度节点(i,J)。这样就会出现这种情况：在速度节点处不存在密度值，而在标量节点处找不到速度值，当在某个确定位置处的某个复合物理量(如式(1-127)中的流通量F)同时需要该处的密度及速度时，要么找不到该处的密度，要么找不到该处的速度。为此，需要在计算过程中通过插值来解决。式(1-127)表明，标量(密度)及速度分量在u控制体积的界面上是不存在的，这时，根据周边的最近邻点的信息，使用二点或四点平均的办法来

处理。

在每次迭代过程中，用于估计上述各表达式的速度分量u和速度分量v是上一次迭代后的数值(在首次迭代时是初始猜测值)。需要特殊说明的是，这些“已知的”速度值u和v也用于计算方程(1-123)中的系数a，但是，它们与方程(1-123)中的待求uj,J和unb是完全不同的。

还需要说明的是，式(1-127)中的线性插值是基于均匀网格而言的，若网格是不均匀的，应该将式(1-127)中的系数2和4等改为相应的网格长度或宽度值的组合。例如，对于不均匀网格上的Fw，按式(1-128)计算：

Fw?(?u)w?xi?xI?1x?xFi,J?I?1i?1Fi?1,Jxi?xi?1xi?xi?1 ??xi?xI?1?xI?xix?x?I,J?iI?1?I?1,J?ui,J??xi?xi?1?xI?xI?1xI?xI?1??xI?1?xi?1?xI?1?xix?x?I?1,J?i?1I?2?I?2,J?ui?1,J?xi?xi?1?xI?1?xI?2xI?1?xI?2? (1-128)

按上述同样的方式，可以写出在新的编号系统中，对于在位置(I,j)处的关于速度vI,j的离散动量方程：

aI,jvI,j??anbvnb?(pI,J?1?pI,J)AI,j?bI,j (1-129)

建立式(1-129)所使用的v控制体积表示在图1-9中。

在求和记号?anbunb中所包含的4个邻点及其主速度也在图1-9中标出。在系数aI,j和 anb中，同样包含在v控制体积界面上的对流质量流量F与扩散传导性D，计算公式如下：

Fw?(?u)w?Fe?(?u)e?Fs?(?v)s?Fi,J?Fi,J?1Fi?1,J???I,J?1?1????I?1,J??I,Jui,J?I?1,J?1ui,J?1? (1-130a) 22?22??Fi?1,J?11??I?1,J??I,J???I?1,J?1???ui?1,J?I,J?1ui?1,J?1? (1-130b)

22?22???I,J?2???I,J?1?1????I,J?1vI,j?1?I,JvI?1,j? (1-130c) 22?22?FI,j?FI,j?11??I,J??I,J?1???I,J?Fn?(?v)n???vI,j?I,J?1vI,j?1? (1-130d)

22?22?FI,j?1?FI,jDw?De?Ds?Dn??I?1,J?1??I,J?1??I?1,J??I,J4(xI?xI?1)?I,J?1??I?1,J?1??I,J??I?1,J4(xI?1?xI) (1-130e) (1-130f)

?I,J?1yj?yj?1 (1-130g) (1-130h)

?I,Jyj?1?xj同样，在每个迭代层次上，用于估计上述各表达式的速度分量u和速度分量v均取上一次迭代后的数值(在首次迭代时是初始猜测值)。

图1-9 v控制体积及其邻点的速度分量

给定一个压力场p，便可针对每个u控制体积和v控制体积写出式(1-123)和式(1-129)所示的动量方程的离散方程，并可以从中求解出速度场。如果压力场是正确的，所得到的速度场将满足连续方程。

2. SIMPLE算法

SIMPLE算法是目前工程上应用最为广泛的一种流场计算方法，它属于压力修正法的一种。

SIMPLE是英文semi-implicit method for pressure-linked equations的缩写，意为“求解压力耦合方程组的半隐式方法”。该方法由Patankar与Spalding于1972年提出，是一种主要用于求解不可压流场的数值方法(也可用于求解可压流动)。它的核心是采用“猜测-修正”的过程，在交错网格的基础上来计算压力场，从而达到求解动量方程(Navier-Stokes方程)的目的。

SIMPLE算法的基本思想可描述如下：对于给定的压力场(它可以是假定的值，或是上一次迭代计算所得到的结果)，求解离散形式的动量方程，得出速度场。因为压力场是假定的或不精确的，这样，由此得到的速度场一般不满足连续方程，因此，必须对给定的压力场加以修正。修正的原则是：与修正后的压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，我们把由动量方程的离散形式所规定的压力与速度的关系代入连续方程的离散形式，从而得到压力修正方程，由压力修正方程得出压力修正值；接着，根据修正后的压力场，求得新的速度场；然后检查速度场是否收敛；若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算；如此反复，直到获得收敛的解。

在上述求解过程中，如何获得压力修正值(即如何构造压力修正方程)，以及如何根据压力修正值确定“正确”的速度(即如何构造速度修正方程)，是SIMPLE算法的两个关键问题。为此，下面先解决这两个问题，然后给出SIMPLE算法的求解步骤。

1) 速度修正方程

考察一直角坐标系下的二维层流稳态问题。设有初始的猜测压力场p*，我们知道，动量方程的离散方程可借助该压力场得以求解，从而求出相应的速度分量u*和v*。

根据动量方程的离散方程(1-123)和(1-129)，有

*aI,JuI*,J??anbunb?(pI*?1,J?pI*,J)AI,J?bI,J (1-131)

*aI,JvI*,J??anbvnb?(pI*,J?1?pI*,J)AI,J?bI,J (1-132)

现在，定义压力修正值p?为正确的压力场p与猜测的压力场p*之差，有

p?p*?p? (1-133)

同样地，定义速度修正值u?和v?，以联系正确的速度场(u,v)与猜测的速度场 (u*,v*)，有

u?u*?u? (1-134)

v?v*?v? (1-135)

将正确的压力场p代入动量离散方程(1-123)与(1-129)，得到正确的速度场(u,v)。现在，

从方程(1-123)和(1-129)中减去方程(1-131)和(1-132)，并假定源项b不变，有

*ai,j(ui,j?ui*,j)??anb(unb?unb)?[(pi?1,j?pi*?1,j)?(pi,j?pi*,l)]Ai,j

*ai,j(vi,j?vi*,j)??anb(vnb?vnb)?[(pi,j?1?pi*,j?1)?(pi,j?pi*,l)]Ai,j

(1-136) (1-137)

引入压力修正值与速度修正值的表达式(1-133)～式(1-135)，方程(1-136)和(1-137)可 写成

??(pi??1,j?pi?,j)Ai,j (1-138) ai,jui?,j??anbunb??(pi?,j?1?pi?,j)Ai,j (1-139) ai,jvi?,j??anbvnb可以看出，由压力修正值p?可求出速度修正值(u?,v?)。式(1-139)还表明，任一点上速度的修正值由两部分组成：一部分是与该速度在同一方向上的相邻两节点间压力修正值之

差，这是产生速度修正值的直接的动力；另一部分是由邻点速度的修正值所引起的，这又可以视为四周压力的修正值对所讨论位置上速度改进的间接影响。

为了简化式(1-138)和(1-139)的求解过程，在此，引入如下近似处理：略去方程中与速

?和?anbvnb?。该近似是SIMPLE算法的重要特征。略去后的影响将度修正值相关的?anbunb在下一节要介绍的SIMPLEC算法中讨论。于是有

ui?,J?di,J(pi??1,J?pi?,J) (1-140)

vI?,j?dI,j(pI?,j?1?pI?,j) (1-141)

以上两种中，

di,J?Ai,Jai,J,dI,j?AI,jaI,j (1-142)

将式(1-140)和式(1-141)所描述的速度修正值，代入式(1-134)和式(1-135)，有

uI,J?uI*,J?dI,J(pI??1,J?pI?,J) (1-143)

vI,J?vI*,j?dI,J(pI?,J?1?pI?,J) (1-144)

对于ui?1,J和vI,j?1，存在类似的表达式：

ui?1,J?ui*?1,J?di?1,J(pi?,J?pi??1,J) (1-145) vI,j?1?vI*,j?1?dI,j?1(pI?,j?pI?,j?1) (1-146)

以上两式中：

di?1,J?Ai?1,Jai?1,J，dI,j?1?AI,j?1aI,j?1 (1-147)

式(1-143)～式(1-147)表明，如果已知压力修正值p?，便可对猜测的速度场(u*,v*)作出相应的速度修正，得到正确的速度场(u,v)。

2) 压力修正方程

在上面的推导中，只考虑了动量方程，其实，如前所述，速度场还受连续方程(1-38)的约束。按本节开头的约定，这里暂不讨论瞬态问题。对于稳态问题，连续方程可写为

?(?u)?(?v)??0 (1-148) ?x?y针对如图1-10所示的标量控制体积，连续方程(1-148)满足如下离散形式：

图1-10 离散连续方程的标量控制体积

[(?uA)i?1,j?(?uA)i,j]?[(?vA)i,j?1?(?vA)i,j]?0 (1-149)

将正确的速度值，即式(1-143)～式(1-147)，代入连续方程的离散方程(1-149)，有

{?i?1,JAi?1,J[ui*?1,J?di?1,J(pI?,J?pI??1,J)]??i,JAi,J[ui*,J?di,J(pI??1,J?pI?,J)]}? (1-150) **{?I,j?1AI,j?1[vI,j?1?dI,j?1(pI?,J?pI?,J?1)]??I,jAI,j[vI,j?dI,j(pI?,J?1?pI?,J)]}?0整理后，得

[(?dA)i?1,J?(?dA)i,J?(?dA)I,j?1?(?dA)I,j]pI?,J?(?dA)i?1,JpI??1,J?(?dA)i,JpI??1,J?(?dA)I,j?1pI?,J?1?(?dA)I,jpI?,J?1? [(?u?A)i,J?(?u?A)i?1,J?(?u?A)I,j?(?u?A)I,j?1](1-151)

该式可简记为

aI,JpI?,J?aI?1,JpI??1,J?aI?1,JpI??1,J?aI,J?1pI?,J?1?aI,J?1pI?,J?1?bI?,J (1-152)

式中：

aI?1,J?(?dA)i?1,J (1-153a) aI?1,J?(?dA)i,J (1-153b) aI,J?1?(?dA)I,j?1 (1-153c)

aI,J?1?(?dA)I,j (1-153d) aI,j?aI?1,J?aI?1,J?aI,J?1?aI,J?1 (1-153e) bI?,J?(?u*A)i,J?(?u*A)i?1,J?(?v*A)I,j?(?v*A)I,j?1 (1-153f)

式(1-152)表示连续方程的离散方程，即压力修正值p'的离散方程。方程中的源项b?是由于不正确的速度场(u*,v*)所导致的“连续性”不平衡量。通过求解方程(1-152)，可得到

空间所有位置的压力修正值p?。

如同处理式(1-127)中的密度一样，式(1-153)中的?是标量控制体积界面上的密度值，同样需要通过插值得到，这是因为密度?是在标量控制体积中的节点(即控制体积的中心)定义和存储的，在标量控制体积界面上不存在可直接引用的值。无论采用何种插值方法，对于交界面所属的两个控制体积，必须采用同样的?值。

为了求解方程(1-152)，还必须对压力修正值的边界条件作出说明。实际上，压力修正方程是动量方程和连续方程的派生物，不是基本方程，故其边界条件也与动量方程的边界条件相联系。在一般的流场计算中，动量方程的边界条件通常有两类：第一，已知边界上的压力(速度未知)；第二，已知沿边界法向的速度分量。若已知边界压力p，可在该段边界上令p*?p，则该段边界上的压力修正值p?应为零。这类边界条件类似于热传导问题中已知温度的边界条件。若已知边界上的法向速度，在设计网格时，最好令控制体积的界面与边界相一致，这样，控制体积界面上的速度为已知。

3) SIMPLE算法的计算步骤

至此，已经得出了求解速度分量和压力所需要的所有方程。根据前面介绍的SIMPLE算法的基本思想，可给出SIMPLE算法的计算流程，如图1-11所示。

图1-11 SIMPLE算法流程图

1.4.2 其他算法介绍

SIMPLE算法自问世以来，在被广泛应用的同时，也以不同方式不断地得到改善和发展，其中最著名的改进算法包括SIMPLEC、SIMPLER和PISO。本节介绍这些改进算法，并对各算法作一对比。

1. SIMPLER算法

SIMPLER是英文SIMPLE revised的缩写，顾名思义是SIMPLE算法的改进版本。它是由SIMPLER算法的创始者之一Patanker完成的。

我们知道，在SIMPLER算法中，为了确定动量离散方程的系数，一开始就假定了一个速度分布，同时又独立地假定了一个压力分布，两者之间一般是不协调的，从而影响了迭代计算的收敛速度。实际上，不必在初始时刻单独假定一个压力场，因为与假定的速度场相协调的压力场是可以通过动量方程求出的。另外，在SIMPLER算法中对压力修正值p?采用了欠松弛处理，而松弛因子是比较难确定的，因此，速度场的改进与压力场的改进不能同步进行，最终影响收敛速度。于是，Patanker便提出了这样的想法：p?只用修正速度，压力场的改进则另谋更合适的方法。将上述两方面的思想结合起来，就构成了SIMPLER算法。

在SIMPLER算法中，经过离散后的连续方程(1-149)用于建立一个压力的离散方程，而不像在SIMPLE算法中用来建立压力修正方程。从而，可直接得到压力，而不需要修正。但是，速度仍需要通过SIMPLE算法中的修正方程即式(1-143)～式(1-147)来修正。

将离散后的动量方程式(1-123)和式(1-129)重新改写后，有

?anbunb?bi,j?Ai,j(p?p) (1-154) ui,j?i?1,ji,jai,jai,jvi,j??anbnbv?bi,jai,j?Ai,jai,j(pi,j?1?pi,j) (1-155)

?与v?如下： 在SIMPLER算法中，定义伪速度uau?bI,J???nbnbu (1-156)

aI,J??v?anbnbv?bI,JaI,J (1-157)

这样，式(1-154)与式(1-155)可写为

?i,j?di,j(pi?1,j?pi,j) (1-158) ui,j?u?i,j?di,j(pi,j?1?pi,j) (1-159) vi,j?v以上两式中的系数d，仍沿用1.4.1节所给出的计算公式。同样可写出ui?1,j与vi,j?1的表达式。然后，将ui,j、vi,j、ui?1,j与vi,j?1的表达式代入离散后的连续方程(1-149)，有

????

i?1,j?i?1,j?di?1,j(pi,j?pi?1,j)]??i,jAi,j[u?i,j?di,j(pi?1,j?pi,j)]??Ai?1,j[u??i,j?1Ai,j?1[ui,j?1?di,j?1(pi,j?pi,j?1)]??i,jAi,j[vi,j?di,j(pi,j?1?pi,j)]??0 (1-160)

整理后，得到离散后的压力方程：

aI,JpI,J?aI?1,JpI?1,J?aI?1,JpI?1,J?aI,J?1pI,J?1?aI,J?1pI,J?1?bI,J (1-161)

式中：

aI?1,J?(?dA)i?1,J (1-162a)

aI?1,j?(?dA)I,J (1-162b) aI,J?1?(?dA)I,j?1 (1-162c) aI,J?1?(?dA)I,j (1-162d)

aI,J?aI?1,J?aI?1,J?aI,J?1?aI,J?1 (1-162e)

?)i,J?(?uA?)i?1,J?(?vA?)I,j?(?vA?)I,j?1 (1-162f) bI,J?(?uA我们注意到，方程(1-161)中的系数与压力修正方程(1-152)中的系数是一样的，差别仅

在于源项b。这里的源项b是用伪速度来计算的。因此，离散后的动量方程式(1-131)和 式(1-132)，可借助上面得到的压力场来直接求解。这样，可求出速度分量u*和v*。SIMPLER算法流程图如图1-12所示。

图1-12 SIMPLER算法流程图

在SIMPLER算法中，初始的压力场与速度场是协调的，且由SIMPLER算法算出的压力场不必作欠松弛处理，迭代计算时比较容易得到收敛解。但在SIMPLER的每一层迭代中，要比SIMPLE算法多解一个关于压力的方程组，一个迭代步内的计算量较大。总体而言，SIMPLER的计算效率要高于SIMPLE算法。

2. SIMPLEC算法

SIMPLEC是英文SIMPLE consistent的缩写，意为协调一致的SIMPLE算法。它也是SIMPLE的改进算法之一，是由Van Doormal和Raithby提出的。

?项，我们知道，在SIMPLE算法中，为求解的方便，略去了速度修正方程中的?anbunb从而把速度的修正完全归结为由于压差项的直接作用。这一作法虽然并不影响收敛解的值，

但加重了修正值p?的负担，使得整个速度场迭代收敛速度降低。实际上，当我们在略去

?anb?时，犯了一个“不协调一致”的错误。为了能略去anbunb?而同时又能使方程基本协unb??ui?,j)?Ai,j(pi??1,j?pi?,j) (1-163) (ai,j??anb)ui?,j??anb(unb调，试在ui?,j方程(1-138)的等号两端同时减去?anbui?,j，有

??ui?,j)所产?具有相同的数量级，可以预期，ui?,j与其邻点的修正值unb因而略去?anb(unb?所产生的影响要小得多。于是有 生的影响远比在方程(1-138)中不计?anbunbui?,j?di,j(pi??1,j?pi?,j) (1-164)

式中：

di,j?Ai,j(ai,j??anb) (1-165)

类似地，有

vi?,j?di,j(pi?,j?1?pi?,j) (1-166)

式中：

di,j?Ai,j(ai,j??anb) (1-167)

将式(1-166)和式(1-167)代入SIMPLE算法中的式(1-143)和式(1-144)，得到修正后的速度计算式：

ui,j?ui*,j?di,j(pi??1,j?pi?,j) (1-168)

vi,j?vi*,j?di,j(pi?,j?1?pi?,j) (1-169)

式(1-168)和式(1-169)在形式上与式(1-143)和式(1-144)一致，只是其中的系数项d的计算公式不同，现在需要按式(1-165)和式(1-167)计算。

这就是SIMPLEC算法。SIMPLEC算法与SIMPLE算法的计算步骤相同，只是速度修正方程中的系数项d的计算公式有所区别。

?项忽略，因此，得到的压力修由于SIMPLEC算法没有像SIMPLE算法那样将?anbunb正值p?一般是比较合适的，因此，在SIMPLEC算法中可不再对p?进行欠松弛处理。但据

作者的试验，适当选取一个稍小于1的?p对p?进行欠松弛处理，对加快迭代过程中解的收敛也是有效的。

3. PISO算法

PISO是pressure implicit with splitting of operators的缩写，意为压力的隐式算子分割算法。PISO算法是Issa于1986年提出的，起初是针对非稳态可压流动的无迭代计算所建立的一种压力速度计算程序，后来在稳态问题的迭代计算中也较广泛地使用了该算法。

PISO算法与SIMPLE、SIMPLEC算法的不同之处在于：SIMPLE和SIMPLEC算法是两步算法，即一步预测(图1-11中的步骤1)和一步修正(图1-11中的步骤2和步骤3)；而PISO算法增加了一个修正步，包含一个预测步和两个修正步，在完成了第一步修正得到(u,v,p)后寻求二次改进值，目的是使它们更好地同时满足动量方程和连续方程。PISO算法由于使用了预测—修正—再修正三步，从而可加快单个迭代步中的收敛速度。现将三个步骤介绍如下。

1) 预测步

使用与SIMPLE算法相同的方法，利用猜测的压力场p*，求解动量离散方程(1-131)与方程式(1-132)，得到速度分量u*与v*。

2) 第一步修正

所得到的速度场(u*,v*)一般不满足连续方程，除非压力场p*是准确的。现引入对SIMPLE的第一个修正步，该修正步给出一个速度场(u**,v**)，使其满足连续方程。此处的修正公式与SIMPLE算法中的式(1-140)和(1-141)完全一致，只不过考虑到在PISO算法还有第二个修正步，因此，使用不同的记法：

p**?p*?p? (1-170) u**?u*?u? (1-171) v**?v*?v? (1-172)

这组公式用于定义修正后的速度u**与v**：

*??ui**,j?ui,j?di,j(pi?1,j?pi,j) (1-173)

*??vi**,j?vi,j?di,j(pi,j?1?pi,j) (1-174)

就像在SIMPLE算法中一样，将式(1-173)与式(1-174)代入连续方程(1-149)，产生与

式(1-151)具有相同系数和源项的压力修正方程。求解该方程，产生第一个压力修正值p?。一旦压力修正值已知，可通过式(1-173)与式(1-174)获得速度分量u**与v**。

3) 第二步修正

为了强化SIMPLE算法的计算，PISO要进行第二步的修正。u**和v**的动量离散方 程是

*****ai,jui**,j??anbunb?(pi?1,j?pi,j)Ai,j?bi,j (1-175)

*****ai,jvi**,j??anbvnb?(pi,j?1?pi,j)Ai,j?bi,j (1-176)

注意这两式实际就是式(1-131)和式(1-132)。为引用方便，给出新的编号。

再次求解动量方程，可以得到两次修正的速度场(u***,v***)：

********ai,jui***,j??anbunb?(pi?1,j?pi,j)Ai,j?bi,j (1-177) ********ai,jvi***,j??anbvnb?(pi,j?1?pi,j)Ai,j?bi,j (1-178)

注意修正步中的求和项是用速度分量u**和v**来计算的。

现在，从式(1-177)中减去式(1-175)，从式(1-178)中减去式(1-176)，有

***anb(unb?unb)?*****??ui,j?ui,j??di,j(pi???1,j?pi,j) (1-179)

ai,jv***i,j?v**i,j?a?nb***(vnb?vnb)ai,j???di,j(pi??,j?1?pi,j) (1-180)

以上两式中，记号p??是压力的二次修正值。有了该记号，p***可表示为

p***?p**?p?? (1-181)

将u***和v***的表达式(1-179)和(1-180)代入连续方程(1-149)，得到二次压力修正方程：

??????????ai,jpi??,j?ai?1,jpi?1,j?ai?1,jpi?1,j?ai,j?1pi,j?1?ai,j?1pi,j?1?bi,j (1-182) 以上两式中，ai,j?ai?1,j?ai?1,j?ai,j?1?ai,j?1。读者可参考建立方程(1-152)同样的过程，写出各系数如下：

ai?1,j?(?dA)i?1,j (1-183a)

ai?1,j?(?dA)i,j (1-183b) ai,j?1?(?dA)i,j?1 (1-183c) ai,j?1?(?dA)i,j (1-183d)

??A?bI?,J????a?i,J??A? ???a?I,??A?***a(u?u)??nbnbnb?a???i?1,J?a?anb***(unb?unb)? (1-183e)

?aj??A?***(v?v)?nbnbnb???a?v(nb*?*vnb *) nbI,?j1下面对源项b?为何是式(1-183e)的形式作一简要分析和解释。

对比建立方程(1-152)的过程，可以发现，式(1-183e)中的各项，是因在u***和v***的表达

******anb(unb?unb)anb(vnb?vnb)??式(1-179)和式(1-180)中存在和项所导致的，而在u和v的表ai,jai,j达式(1-142)和式(1-143)中没有这样的项，因此，式(1-152)不存在类似式(1-183e)中的各项。

但式(1-152)存在另外一个源项，即[(?u*A)i,j?(?u*A)i?1,j?(?v*A)i,j?(?v*A)i,j?1]，这是因速度u和v的表达式(1-143)和式(1-144)中的u*与v*项所导致的。按此推断，在式(1-183e)中也应该存在类似表达式[(?u**A)i,j?(?u**A)i?1,j?(?v**A)i,j?(?v**A)i,j?1]。但是，由于u**和v**满

********足连续方程，因此??(?uA)i,j?(?uA)i?1,j?(?vA)i,j?(?vA)i,j?1??为0。

现在，求解方程(1-182)，就可得到二次压力修正值p??。这样，通过式(1-184)就可得到二次修正的压力场：

p***?p**?p???p*?p??p?? (1-184)

最后，求解方程(1-179)与(1-180)，得到二次修正的速度场。

在瞬态问题的非迭代计算中，压力场p***与速度场(u***,v***)被认为是准确的。对于稳态流动的迭代计算，PISO算法的实施过程如图1-13所示。

PISO算法要两次求解压力修正方程，因此，它需要额外的存储空间来计算二次压力修正方程中的源项。尽管该方法涉及较多的计算，但对比发现，它的计算速度很快，总体效率比较高。Fluent的用户手册[7]推荐，对于瞬态问题，PISO算法有明显的优势；而对于稳

态问题，可能选SIMPLE或SIMPLEC算法更合适。

图1-13 PISO算法流程图

4. SIMPLE系列算法的比较

SIMPLE算法是SIMPLE系列算法的基础，目前在各种CFD软件中均提供这种算法。SIMPLE的各种改进算法，主要是提高了计算的收敛性，从而可缩短计算时间。

在SIMPLE算法中，压力修正值p?能够很好地满足速度修正的要求，但对压力修正不是十分理想。改进后的SIMPLER算法只用压力修正值p?来修正速度，另外构建一个更加有效的压力方程来产生“正确”的压力场。由于在推导SIMPLER算法的离散化压力方程时，没有任何项被忽略，因此所得到的压力场与速度场相适应。在SIMPLER算法中，正确的速度场将导致正确的压力场，而在SIMPLE算法中则不是这样。所以SIMPLER算法是在很高的效率下正确计算压力场的，这一点在求解动量方程时有明显优势。虽然SIMPLER算法的计算量比SIMPLE算法高出30%左右，但其较快的收敛速度使得计算时间减少30%～50%。

SIMPLEC算法和PISO算法总体上与SIMPLER算法具有同样的计算效率，相互之间很难区分谁高谁低，对于不同类型的问题每种算法都有自己的优势。一般来讲，动量方程与标量方程(如温度方程)如果不是耦合在一起的，则PISO算法在收敛性方面显得很健壮，且效率较高。而在动量方程与标量方程耦合非常密切时，SIMPLEC和SIMPLER算法的效果可能更好些。

另外还有瞬态问题的算法、同位网格的算法等，限于篇幅，请读者查阅相关书籍。

本文来源：https://www.bwwdw.com/article/64mw.html