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），确定每个子区域中的节点位置及该节点所代表的控制体积（control volume）。区域离散后，得到以下四种几何要素：

节点（node）：需要求解的未知物理量的几何位置；

控制体积（control volume）：应用控制方程或守恒定律的最小几何单位；

界面（face）：它定义了与各节点相对应的控制体积的界面位置；

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

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

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

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

2．控制方程的离散化

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

式中，div表示散度（计算方法如式（1-39）所示）；grad表示梯度（计算方法如式（1-38）所示）。

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

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

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

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

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

为了建立所需要形式的离散方程，需要找出如何表示式（1-98）中界面e和w处的ρ、u、Γ、ϕ和。在有限体积法中规定，ρ、u、Γ、ϕ和等物理量均是在节点处定义和计算的。因此，为了计算界面上的这些物理参数（包括其导数），需要一个物理参数在节点间的近似分布。可以想象，线性近似是可以用来计算界面物性值的最直接，也是最简单的方式，这种分布叫做中心差分。如果网格是均匀的，则单个物理参数（以扩散系数Γ为例）的线性插值结果是

（ρuϕA）的线性插值结果是

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

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

式中，SC是常数；SP是随时间和物理量ϕ变化的项。将式（1-100）～式（1-102）代入方程（1-98），有

整理后得

记为

式中，

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

方程（1-104）即为方程（1-96）的离散形式，每个节点上都可建立此离散方程，通过求解方程组，就可得到各物理量在各节点处的值。

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

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

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

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

将式（1-108）代入方程（1-106）有

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

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

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

将方程（1-111）在一维计算网格上对时间及控制体积进行积分有

改写后，有

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

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

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

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

将ϕP、ϕE、ϕW及SC＋SPϕP采用类似式（1-115）进行时间积分，式（1-114）可写成：

整理后得

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

将式（1-118）代入方程（1-117），得

同样也像稳态问题引入ap、aw和aE，上式变为

式中，

根据f的取值，瞬态问题对时间的积分有几种方案，当f＝0时，变量的初值出现在方程（1-120）的右端，从而可直接求出在现时刻的未知变量值，这种方案称为显式时间积分方案。当0＜f＜1时，有现时刻的未知变量出现在方程（1-120）的两端，需要解若干个方程组成的方程组才能求出现时刻的变量值，这种方案称为隐式时间积分方案。当f＝1时，就成为了全隐式时间积分方案了。当f＝1/2时，称为Crank-Nicolson时间积分方案。

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

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

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

1．常用的离散格式

有限体积法常用的离散格式有：中心差分格式、一阶迎风格式、混合格式、指数格式、乘方格式、二阶迎风格式、QUICK格式。各种离散格式对一维、稳态、无源项的对流-扩散问题的通用控制方程（1-126）均能得到式（1-127）的形式。对于高阶情况如式（1-128）所示。

式中，对于一阶情况，aP＝aW＋aE＋（Fe－Fw），对于二阶情况，aP＝aW＋aE＋aWW＋aEE＋（Fe－Fw），其中系数aW和aE取决于所使用的离散格式（高阶还有aWW和aEE），为了便于比较和编程计算，将各种离散格式的系数总结于表1-2中。

2．建立离散方程应遵循的原则

在利用有限体积法建立离散方程时，必须遵守以下四条基本原则：

1）控制体积界面上的连续性原则

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

2）正系数原则

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

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

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

4）系数ap等于相邻节点系数之和原则

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