1.2 计算流体力学（CFD）基础

计算流体动力学（Computational Fluid Dynamics, CFD）是近代流体力学、数值数学和计算机科学结合的产物，是一门具有强大生命力的边缘科学。

1.2.1 CFD概述

CFD以电子计算机为工具，应用各种离散化的数学方法，对流体力学的各类问题进行数值实验、计算机模拟和分析研究，以解决各种实际问题。

计算流体力学和相关的计算传热学、计算燃烧学的原理是用数值方法求解非线性联立的质量、能量、组分、动量和自定义的标量的微分方程组，求解结果能预报流动、传热、传质、燃烧等过程的细节，并成为过程装置优化和放大定量设计的有力工具。计算流体力学的基本特征是数值模拟和计算机实验，它从基本物理定理出发，在很大程度上替代了耗资巨大的流体动力学实验设备，在科学研究和工程技术中产生巨大的影响。

计算流体力学是目前国际上一个强有力的研究领域，是进行传热、传质、动量传递及燃烧、多相流和化学反应研究的核心和重要技术，广泛应用于航天设计、汽车设计、生物医学工业、化工处理工业、涡轮机设计、半导体设计、HAVC&R等诸多工程领域，板翅式换热器设计是CFD技术应用的重要领域之一。

CFD在最近20年中得到飞速的发展，除了计算机硬件工业的发展给它提供了坚实的物质基础外，还主要因为无论分析的方法或实验的方法都有较大的限制。例如，由于问题的复杂性，既无法作分析解，也因费用昂贵而无力进行实验确定，而CFD的方法正具有成本低和能模拟较复杂或较理想的过程等优点。

经过一定考核的CFD软件可以拓宽实验研究的范围，减少成本昂贵的实验工作量。在给定的参数下用计算机对现象进行一次数值模拟相当于进行一次数值实验，历史上也曾有过首先由CFD数值模拟发现新现象而后由实验予以证实的例子。

CFD软件一般都能推出多种优化的物理模型，如定常和非定常流动、层流、紊流、不可压缩和可压缩流动、传热、化学反应等。对每一种物理问题的流动特点，都有适合它的数值解法，用户可选择显式或隐式差分格式，以期在计算速度、稳定性和精度等方面达到最佳。

CFD软件之间可以方便地进行数值交换，并采用统一的前、后处理工具，这就省去了科研工作者在计算机方法、编程、前后处理等方面投入的重复、低效的劳动，而可以将主要精力和智慧用于物理问题本身的探索上。

1.2.2 CFD求解力学问题的过程

所有CFD问题的求解过程都可用图1-4表示。如果所求解的问题是瞬态问题，则可将图1-4的过程理解为一个时间步的计算过程，循环这一过程求解下个时间步的解。下面对各求解步骤进行简单介绍。

1．建立控制方程

建立控制方程是求解任何问题前都必须首先进行的。一般来讲，这一步是比较简单的。因为对于一般的流体流动而言，可直接写出其控制方程。假定没有热交换发生，则可直接将连续方程与动量方程作为控制方程使用。一般情况下，需要增加湍流方程。

2．确定边界条件和初始条件

初始条件与边界条件是控制方程有确定解的前提，控制方程与相应的初始条件、边界条件的组合构成对一个物理过程完整的数学描述。

初始条件是所研究对象在过程开始时刻各个求解变量的空间分布情况。对于瞬态问题，必须给定初始条件。对于稳态问题，不需要初始条件。

边界条件是在求解区域的边界上所求解的变量或其导数随地点和时间的变化规律。对于任何问题，都需要给定边界条件。

3．划分计算网格

采用数值方法求解控制方程时，都是想办法将控制方程在空间区域上进行离散，然后求解得到离散方程组。要想在空间域上离散控制方程，必须使用网格。现已发展出多种对各种区域进行离散以生成网格的方法，这些方法统称为网格生成技术。

不同的问题采用不同数值解法时，所需要的网格形式是有一定区别的，但生成网格的方法基本是一致的。目前网格分结构网格和非结构网格两大类。

简单地讲，结构网格在空间上比较规范，如对一个四边形区域，网格往往是成行成列分布的，行线和列线比较明显。而非结构网格在空间分布上没有明显的行线和列线。

对于二维问题，常用的网格单元有三角形和四边形等形式；对于三维问题，常用的网格单元有四面体、六面体、三菱体等形式。在整个计算域上，网格通过节点联系在一起。

目前各种CFD软件都配有专用的网格生成工具，如FLUENT使用Gambit作为前处理软件。多数CFD软件可接收采用其他CAD或CFD/FEM软件产生的网格模型。例如，FLUENT可以接收ANSYS所生成的网格。

4．建立离散方程

对于在求解域内所建立的偏微分方程，理论上是有真解（或称精确解或解析解）的。但由于所处理问题自身的复杂性，一般很难获得方程的真解。因此就需要通过数值方法把计算域内有限数量位置（网格节点或网格中心点）上的因变量值当作基本未知量来处理，从而建立一组关于这些未知量的代数方程组，然后通过求解代数方程组来得到这些节点值，而计算域内其他位置上的值则根据节点位置上的值来确定。

由于所引入的应变量在节点之间的分数假设及推导离散化方程的方法不同，所以形成了有限差分法、有限元法、有限体积法等不同类型的离散化方法。

对于瞬态问题，除了在空间域上的离散外，还要涉及在时间域上的离散。离散后，将要涉及使用何种时间积分方案的问题。

5．离散初始条件和边界条件

前面所给定的初始条件和边界条件是连续性的，如在静止壁面上速度为0，现在需要针对所生成的网格，将连续型的初始条件和边界条件转化为特定节点上的值，如静止壁面上共有90个节点，则这些节点上的速度值应均设为0。

商用CFD软件往往在前处理阶段完成网格划分后，直接在边界上指定初始条件和边界条件，然后由前处理软件自动将这些初始条件和边界条件按离散的方式分配到相应的节点上。

6．给定求解控制参数

在离散空间上建立了离散化的代数方程组，并施加离散化的初始条件和边界条件后，还需要给定流体的物理参数和湍流模型的经验系数等。此外，还要给定迭代计算的控制精度、瞬态问题的时间步长和输出频率等。

7．求解离散方程

进行上述设置后，生成了具有定解条件的代数方程组。对于这些方程组，数学上已有相应的解法，如线性方程组可采用Gauss消去法或Gauss-Seidel迭代法求解，而对于非线性方程组，可采用Newton-Raphson方法。

商用CFD软件往往提供多种不同的解法，以适应不同类型的问题。这部分内容属于求解器设置的范畴。

8．显示计算结果

通过上述求解过程得出了各计算节点上的解后，需要通过适当的手段将整个计算域上的结果表示出来，这时，可采用线值图、矢量图、等值线图、流线图、云图等方式来表示计算结果。

● 线值图是指在二维或三维空间上，将横坐标取为空间长度或时间历程，将纵坐标取为某一物理量，然后用光滑曲线或曲面在坐标系内绘制出某一物理量沿空间或时间的变化情况。

● 矢量图是直接给出二维或三维空间里矢量（如速度）的方向及大小，一般用不同颜色和长度的箭头表示速度矢量。矢量图可以比较容易地让用户发现其中存在的漩涡区。

● 等值线图是用不同颜色的线条表示相等物理量（如温度）的一条线。

● 流线图是用不同颜色的线条表示质点运动轨迹。

● 云图是使用渲染的方式，将流场某个截面上的物理量（如压力或温度）用连续变化的颜色块表示其分布。

1.2.3 CFD数值模拟方法和分类

CFD的数值解法有很多分支，这些方法之间的区别主要在于对控制方程的离散方式。根据离散原理的不同，CFD大体上可以分为有限差分法（FDM）、有限元法（FEM）和有限体积法（FVM）。

1．有限差分法

有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。

有限差分法以Taylor级数展开的方法，把控制方程中的导数用网格节点上的函数值的差商代替，从而创建以网格节点上的值为未知数的代数方程组。

该方法是一种直接将微分问题变为代数问题，从而可以用近似数值解法求解，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

从有限差分格式的精度来划分，有一阶格式、二阶格式和高阶格式；从差分的空间形式来考虑，可分为中心格式和逆风格式；考虑时间因子的影响，差分格式还可分为显格式、隐格式、显隐交替格式等。

目前常见的差分格式主要是上诉几种格式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际情况和柯郎稳定条件决定。

2．有限元法

有限元法（FEM）的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。

采用不同的权函数和插值函数形式，便于构成不同的有限元方法。有限元法最早应用于结构力学，后来随着计算机的发展逐渐用于流体力学的数值模拟。

在有限元法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线性组合来逼近单元中的真解，整个计算域上总体的基函数可以看作由每个单元基函数组成，而整个计算域内的解可以看作由所有单元上的近似解构成。

有限元方法的基本思路和解题步骤可归纳如下。

（1）建立积分方程。根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。

（2）区域单元剖分。根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，并列出自然边界和本质边界的节点序号和相应的边界值。

（3）确定单元基函数。根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，所以在选取基函数时可遵循一定的法则。

（4）单元分析。将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数（即单元中各节点的参数值）的代数方程组（称为单元有限元方程）。

（5）总体合成。在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。

（6）边界条件的处理。一般边界条件有3种形式，分别为本质边界条件（狄里克雷边界条件）、自然边界条件（黎曼边界条件）、混合边界条件（柯西边界条件）。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。

（7）解有限元方程。根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。

3．有限体积法

有限体积法（Finite Volume Method, FVM）又称为控制体积法。其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。

其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值分段分布的剖面。

从积分区域的选取方法看来，有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简而言之，子区域法属于有限体积法的基本方法。

有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积都得到满足，在整个计算区域，自然也就得到满足一样，这是有限体积法吸引人的优点。

某些离散方法，如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。

就离散方法而言，有限体积法可视为有限单元法和有限差分法的中间物，有限单元法必须假定值符号网格点之间的变化规律（即插值函数），并将其作为近似解；有限差分法只考虑网格点上的数值而不考虑其在网格点之间如何变化；有限体积法只寻求节点值，这与有限差分法相类似，但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。

在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。

1.2.4 有限体积法计算区域的离散

从前面的介绍中可以看出，有限体积法是一种分块近似的计算方法，因此其中比较重要的步骤是计算区域的离散和控制方程的离散。

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

● 节点：需要求解的未知物理量的几何位置。

● 控制体积：应用控制方程或守恒定律的最小几何单位。

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

● 网格线：连接相邻两节点面形成的曲线簇。

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

计算区域离散的网格有两类：结构化网格和非结构化网格。

结构化网格（structured grid）的节点排列有序，即给出了一个节点的编号后，立即可以得出其相邻节点的编号，所有内部节点周围的网格数目相同。

结构化网格具有实现容易、生成速度快、网格质量好、数据结构简单化的优点，但不能实现复杂边界区域的离散。

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

1.2.5 有限体积法控制方程的离散

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

对于一维稳态问题，其控制方程如式1-206所示。

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

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

式1-207中，ΔV是控制体积的体积值。当控制体积很微小时，ΔV可以表示为ΔV·A，这里A是控制体积界面的面积。从而有

从式1-208可以看到，对流项和扩散项均已转化为控制体积界面上的值。有限体积法最显著的特点之一就是离散方程中具有明确的物理插值，即界面的物理量要通过插值的方式由节点的物理量来表示。

为了建立所需形式的离散方程，需要找出如何表示式1-208中界面e和w处的ρ、u、Γ、φ为了计算界面上的这些物理参数（包括其导数），需要一个物理参数在节点间的近似分布。和。在有限体积法中规定，ρ、u、Γ、φ和等物理量均是在节点处定义和计算的。因此，

可以想象，线性近似是可以用来计算界面物性值的最直接，也是最简单的方式。这种分布叫作中心差分。如果网格是均匀的，则单个物理参数（以扩散系数Γ为例）的线性插值结果是

(ρuφA)的线性插值结果是

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

对于源项S，它通常是时间和物理量φ的函数。为了简化处理，将S转化为如下线性方式：

式中，SC是常数，SP是随时间和物理量φ变化的项。将式1-210～式1-212代入式1-208，有

整理后得到

记为

a PφP=aWφW+aEφE+b

式1-214中

对于一维问题，控制体积界面 e 和 w 处的面积Ae和Aw均为1，即单位面积。这样ΔV=Δx，式1-215各系数可转化为

根据aPφP=aWφW+aEφE+b，每个节点上都可建立此离散方程，通过求解方程组，就可得到各物理量在各节点处的值。

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

这样，F和D在控制界面上的值分别为

在此基础上，定义一维单元的Peclet数Pe为

式1-219中，Pe表示对流与扩散的强度之比。当Pe为0时，对流扩散演变为纯扩散问题，即流场中没有流动，只有扩散；当Pe>0时，流体沿x方向流动，当Pe数很大时，对流扩散问题演变为纯对流问题。一般在中心差分格式中，有Pe<2的要求。

将式1-218代入式1-216中，有

瞬态问题与稳态问题相似，主要是瞬态项的离散。其一维瞬态问题的通用控制方程为

该方程是一个包含瞬态及源项的对流扩散方程。从左到右，方程中的各项分别是瞬态项、对流项、扩散项及源项。方程中的φ是广义变量，如速度分量、温度、浓度等，Γ为相应于φ的广义扩散系数，S为广义源项。

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

将式1-221在一维计算网格上对时间及控制体积进行积分，有

整理后，有

式1-223中，A是图1-5中控制体积P的界面处的面积。

在处理瞬态项时，假定物理量φ在整个控制体积P上均具有节点处值φP，并用线性插值(φP-φP0)/Δt来表示∂φ/∂t。源项也分解为线性方式S=SC+SPφP。对流项和扩散项的值按中心差分格式通过节点处的值来表示，则有

假定变量φP对时间的积分为

式1-225中，上标0代表t时刻；φP是时刻的值；f为0与1之间的加权因子，当f=0时，变量取老值进行时间积分，当f=1时，变量采用新值进行时间积分。将φP、φE、φW及SC+SPφP进行时间积分，由式1-224可得到

整理后得

扩展F和D的定义，即乘以面积A，有

代入方程1-227，得

同样也向稳态问题引入a P、aW和a E，式1-229变为

式1-230中

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

当0< f <1时，现时刻的未知变量出现在方程的两端，需要解由若干方程组成的方程组才能求出现时刻的变量值，这种方案称为隐式时间积分方案。

● 当f =1时，称为全隐式时间积分方案。

● 当f =1/2时，称为Crank-Nicolson时间积分方案。

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

式1-232中

在三维问题中，计算区域离散如图1-7所示（两个方向的投影）。在二维的基础上增加第三坐标z，控制体积增加的前后界面，分别用t（top）和b（bottom）表示，相应的两个邻点记为T和B。在全隐式时间积分方案下的三维瞬态对流扩散问题的离散方程为

式1-234中

有限体积法常用的离散格式有中心差分格式、一阶迎风格式、混合格式、指数格式、乘方格式、二阶迎风格式、QUICK格式。

各种离散格式对一维、稳态、无源项的对流-扩散问题的通用控制方程（如式1-236所示），均能得到如式1-237所示的形式，对于高阶情况如式1-238所示。

式1-238中，对于一阶情况，a P=aW+a E+(Fe-Fw)，对于二阶情况，a P=aW+a E+WWa EE+a +(Fe-Fw)，其中系数aW和aE取决于所使用的离散格式（高阶还有aWW和aEE）。各种离散格式系数aW和aE的计算公式见表1-5。

1.2.6 CFD常用算法

流场计算的基本过程是在空间上用有限体积法（或其他类似方法）将计算区域离散成许多小的体积单元，在每个体积单元上对离散后的控制方程组进行求解。

其本质是对离散方程进行求解，一般可以分为分离解法（segregated method）和耦合解法（coupled method）两大类，各自又根据实际情况扩展成具体的计算方法。

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

（1）假定初始压力场。

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

（3）利用速度场求解连续方程，使压力场得到修正。

（4）根据需要，求解湍流方程及其他标量方程。

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

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

耦合解法同时求解离散方程组，联立求解出各变量（等），其求解过程如下。

（1）假定初始压力和速度等变量，确定离散方程的系数及常数项等。

（2）联立求解连续方程、动量方程、能量方程。

（3）求解湍流方程及其他标量方程。

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

1.SIMPLE算法

SIMPLE算法就是求解压力耦合方程的半隐方法（Semi-Implicit Method for Pressure Linked Equations），它是Patankar与Spalding在1972年提出的。

在常规离散方法中，压力梯度项的离散会遇到问题。

如图1-8所示，对P控制体积分后，的贡献为Pw-Pe，如w和e为单元中点，则

因此，动量方程将包含相间隔（而非相邻）节点间的压力差。

这样导致求解精度降低，且形成锯齿状的压力场，如图1-9所示。

这类锯齿状压力场对动量方程而言与均匀场相同（奇偶差），因此，高度不均匀的压力场将被动量方程的特殊离散化当作均匀的压力场处理。

连续方程的离散也会遇到问题，如对于一维问题。

对控制体积分后：

ue-uw= 0

即

则

uE-uW= 0（跳开了P点）

这样也会导致奇偶差。

对以上出现的离散问题，用交错网格法能较好地解决，如图1-10所示。在此方法中，将速度变量u、v直接设置在P控制体的边界面上，即P控制体边界面上的u、v不再通过主节点上的值求得，而是直接解得。

SIMPLE算法的基本思想可描述为：对于给定的压力场（它可以是假定的值，或是上一次迭代计算所得到的结果），求解离散形式的动量方程，得出速度场。因为压力场是假定的或不精确的，这样，由此得到的速度场一般不满足连续方程，因此，必须对给定的压力场加以修正。

修正的原则是：与修正后的压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把由动量方程离散形式所规定的压力与速度的关系代入连续方程的离散形式，从而得到压力修正方程，由压力修正方程得出压力修正值，接着根据修正后的压力场，求得新的速度场，然后检查速度场是否收敛，若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算；如此反复，直到获得收敛的解。

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

（1）速度修正方程。

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

根据动量方程的离散方程，有

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

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

将正确的压力场p代入动量离散方程，得到正确的速度场(u, v)。现在，通过离散的动量方程与式1-240和式1-241，并假定源项b不变，有

引入压力修正值与速度修正值的表达式1-242～式1-244，方程1-245和方程1-246可写成

可以看出，由压力修正值p＇可求出速度修正值(u＇,v＇)。式1-247和式1-248还表明，任意一点上速度的修正值由两部分组成：一部分是与该速度在同一方向上的相邻两节点间压力修正值之差，这是产生速度修正值的直接动力；另一部分是由邻点速度的修正值所引起的，这又可以视为四周压力的修正值对所讨论位置上速度改进的间接影响。

为了简化方程1-247和方程1-248的求解过程，在此，引入如下近似处理：略去方程中与速度修正值相关的和。该近似是SIMPLE算法的重要特征。略去后的影响将在后面介绍的SIMPLEC算法中讨论。于是有

以上两式中

将式1-249和式1-250所描述的速度修正值代入方程1-243和方程1-244，有

对于和，存在类似的表达式

以上两式中

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

（2）压力修正方程。

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

针对图1-11所示的标量控制体积，连续方程满足如下离散形式。

将正确的速度值，即式1-252～式1-256，代入连续方程的离散方程1-258，有

整理后，得

式1-260可简化为

式1-261中

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

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

为了求解方程1-261，还必须对压力修正值的边界条件作出说明。实际上，压力修正方程是动量方程和连续方程的派生物，不是基本方程，故其边界条件也与动量方程的边界条件相联系。

在一般的流场计算中，动量方程的边界条件通常有两类。

第一类，已知边界上的压力（速度未知）。

第二类，已知沿边界法向的速度分量。

若已知边界压力，可在该段边界上令p*=，则该段边界上的压力修正值p＇应为零。这类边界条件类似于热传导问题中已知温度的边界条件。

若已知边界上的法向速度，在设计网格时，最好令控制体积的界面与边界相一致，这样，控制体积界面上的速度为已知。

（3）SIMPLE算法的计算步骤。

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

2．其他算法

基于SIMPLE算法的改进算法包括SIMPLEC、SIMPLER和PISO。下面介绍这些改进算法，并对各算法进行对比。

（1）SIMPLER算法。

SIMPLER是英文SIMPLE revised的缩写，是SIMPLE算法的改进版本。它是由SIMPLER算法的创始人之一Patankar完成的。

在SIMPLER算法中，为了确定动量离散方程的系数，一开始就假定了一个速度分布，同时又独立地假定了一个压力分布，两者之间一般是不协调的，从而影响了迭代计算的收敛速度。

实际上，与假定的速度场相协调的压力场是可以通过动量方程求出的，因此不必在初始时刻单独假定一个压力场。

另外，在SIMPLER算法中对压力修正值p＇采用了欠松弛处理，而松弛因子是比较难确定的，因此，速度场的改进与压力场的改进不能同步进行，最终影响收敛速度。于是，Patankar便提出了这样的想法：p＇只用修正速度，压力场的改进则另谋更合适的方法。将上述两方面的思想结合起来，就构成了SIMPLER算法。

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

将离散后的动量方程式重新改写后，有

在SIMPLER算法中，定义伪速度与为

这样，式1-262和式1-263可以变为

以上两式中的系数d，仍沿用式1-251。同样可写出与的表达式。然后将、、与的表达式代入离散后的连续方程1-258，有

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

式1-269中

方程1-269中的系数与前面的压力修正方程中的系数相同，差别仅在于源项b。这里的源项 b 是用伪速度来计算的。因此，离散后的动量方程式，可借助上面得到的压力场来直接求解。这样，可求出速度分量u*和v*。SIMPLER算法流程图如图1-13所示。

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

（2）SIMPLEC算法。

SIMPLEC是英文SIMPLE consistent的缩写，意为协调一致的SIMPLE算法。它也是SIMPLE的改进算法之一。

在SIMPLE算法中，为求解的方便，略去了速度修正方程中的项，从而把速度的修正完全归结为由于压差项的直接作用。这一做法虽然不影响收敛解的值，但加重了修正值p＇的负担，使得整个速度场迭代收敛速度降低。

在略去时，犯了一个“不协调一致”的错误。为了能略去而同时又能使方程基本协调，在方程1-247的等号两端同时减去，有

与其邻点的修正值具有相同的数量级，因而略去所产生的影响远比在方程1-247中不计所产生的影响要小得多。于是有

式1-271中

类似地，有

式1-273中

将式1-273和式1-274代入式1-252和式1-253，得到修正后的速度计算式：

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

由于SIMPLEC算法没有像SIMPLE算法那样将项忽略，因此，得到的压力修正值p＇一般是比较合适的，因此，在SIMPLEC算法中可不再对p＇进行欠松弛处理。但据作者的经验，适当选取一个稍小于1的ap对p＇进行欠松弛处理，对加快迭代过程中解的收敛也是有效的。

（3）PISO算法。

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

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

① 预测步。

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

② 第一步修正。

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

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

就像在SIMPLE算法中一样，将式1-280和式1-281代入连续方程1-258，得到压力修正方程。求解该方程，产生第一个压力修正值p＇。一旦压力修正值已知，可通过式1-280和式1-281获得速度分量u**与v**。

③ 第二步修正。

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

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

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

从式1-284中减去式1-282，从式1-285中减去式1-283，有

式1-286和式1-287中，记号p＇＇是压力的二次修正值。有了该记号，p***可表示为

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

式1-289中，。可写出各系数如下。

现在，求解方程1-289，就可得到二次压力修正值p＇＇，从而可得到二次修正的压力场。

最后，求解方程1-286与方程1-287，得到二次修正的速度场。

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

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

3.SIMPLE系列算法的比较

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

在SIMPLE算法中，压力修正值p＇能够很好地满足速度修正的要求，但压力修正不是十分理想。改进后的SIMPLER算法只用压力修正值p＇来修正速度，另外构建一个更加有效的压力方程来产生“正确”的压力场。

由于在推导SIMPLER算法的离散化压力方程时，没有任何项被忽略，因此所得到的压力场与速度场相适应。

在SIMPLER算法中，正确的速度场将导致正确的压力场，而在SIMPLE算法中则不是这样。因此SIMPLER算法是在很高的效率下正确计算压力场的，这一点在求解动量方程时有明显优势。

虽然SIMPLER算法的计算量比SIMPLE算法高出30%左右，但其较快的收敛速度使得计算时间减少30%～50%。

SIMPLEC算法和PISO算法总体上与SIMPLER算法具有同样的计算效率，相互之间很难区分谁高谁低，对于不同类型的问题每种算法都有自己的优势。

一般来讲，动量方程与标量方程（如温度方程）如果不是耦合在一起的，则PISO算法在收敛性方面显得很好，且效率较高。而在动量方程与标量方程耦合非常密切时，SIMPLEC和SIMPLER算法的效果可能更好些。

1.2.7 计算域网格生成技术

网格分布是流动控制方程数值离散的基础，因此，网格生成技术是CFD成功实现数值模拟的关键前提之一，网格质量的好坏直接影响到计算的敛散性及结果的精度。

网格生成的难度和耗费在整个模拟计算程序中占有较大的比重，“从某种角度看，自动生成绕复杂外形的理想网格甚至比编制一个三维的流场解程序更困难，即使在CFD高度发达的国家，网格生成仍占一个计算任务全部人力与时间的60%”，由此可见网格生成在实现流场解算功能过程中的重要性。

网格生成的实质是物理求解域与计算求解域的转换。一般而言，物理域与计算域间的转换应满足下述基本条件。

（1）生成的网格使物理求解域上的计算节点与计算求解域上的计算节点一一对应，不致于出现物理对应关系不确定的多重映射节点。

（2）生成的网格能够准确反映求解域的复杂几何边界形状变化，能够便于边界条件的处理。

（3）物理求解域上的网格应连续光滑求导，保证控制方程离散过程中一阶甚至多阶偏导数的存在性、连续性。网格中出现的尖点、突跃点都将导致算法发散。

（4）网格的疏密易于控制，能够在气动参数变化剧烈的位置，如击波面、壁面等处加密网格，而在气动参数变化平缓的位置拉疏网格。

计算机技术的发展为计算流体力学步入工程实用阶段提供了可能，如何有效地处理复杂的物面边界，生成高质量的计算网格，是目前计算流体力学一个重要的研究课题。

结构化网格在拓扑结构上相当于矩形域内的均匀网格，其节点定义在每一层的网格线上，因此对于复杂外形物体要生成贴体的结构网格是比较困难的。而非结构化网格节点的空间分布完全是随意的，没有任何结构特性，适应性强，因而适合于处理复杂几何外形，并且由于非结构化网格在其生成过程中都要采用一定的准则进行优化判定，因而生成的网格质量较高。

Wimslow最早在20世纪60年代利用有限面积法用三角形网格对Poisson方程进行了数值求解；而20世纪90年代以后，国外学者Dawes、Hah、Prekwas等采用了非结构化网格进行流场的数值求解，国外的著名商业CFD软件，如FLUENT、Star-CD等在20世纪90年代后都将结构化网格计算方法推广到非结构化网格上。

近年来，国内学者也开始对非结构化网格进行深入的研究与探讨，由此可见非结构化网格已成为目前计算流体力学学科中的一个重要方向。

非结构化网格的生成方法有很多，较常用到的有两种：Delaunay三角化方法和推进阵面法（advancing front method）。

（1）Delaunay三角化方法。

Delaunay三角化方法是按一定的方式在控制体内布置节点。定义一个凸多边形外壳，将所有的点都包含进去，并在外壳上进行三角化的初始化。

将节点逐个加入已有的三角化结构中，根据优化准则破坏原有的三角化结构，并建立新的三角化结构，对有关数据结构进行更新后，继续加点，直到所有的节点都加入其中，三角化过程结束。所生成的三角形网格，可以通过光顺技术进一步提高质量。

常用的一种网格光顺方法称为Laplacian光顺方法，这种光顺方法是通过将节点向这个节点周围的三角形所构成的多边形的形心的移动来实现的。

（2）推进阵面法。

推进阵面法是网格和节点同时生成的一种生成方法，它的基本方法为：根据网格密度控制的需要，在平面上布置一些控制点，给每一个控制点定义一个尺度。根据这些控制点将平面划分成大块的三角形背景网格。

每一个背景网格中的所有点的尺度都可以根据其3个顶点的尺度插值得到。因此相当于布置了一个遍布整个平面的网格尺度函数。根据内边界定义初始阵面，按顺时针方向进行阵面初始化。初始化阵面上的每一条都称为活动边。

由初始阵面上的一条活动边开始推进，根据该活动边的中点所落入的背景网格插值确定该点的尺度，根据该点的尺度及有关的规则确定将要生成的节点位置。

判定该节点是否应被接纳，并根据情况生成新的三角形单元，更新阵面，并沿阵面的方向继续推进生成三角形，直至遇到外边界，网格生成结束。