1.3　有限体积法

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

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

1.3.1　计算区域离散化与控制方程离散化

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

1.计算区域的离散化

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

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

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

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

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

一般把节点看成是控制体积的代表。在离散过程中，将一个控制体积上的物理量定义并存储在相应节点处。图1-2所示为一维的有限体积法网格，图1-3所示为二维的有限体积法网格。

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

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

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

2.控制方程的离散化

前面给出的关于流体流动问题的控制方程，无论是连续性方程、运动方程，还是能量方程，都可写成

（1-60）

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

对于一维稳态问题，其控制方程为

（1-61）

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

有限体积法的关键一步是在控制体积上建立积分控制方程，以在控制体积节点上产生离散的方程。通过一维模型方程（1-61），在图1-2所示的控制体积P上积分，可得

（1-62）

式中，是控制体积的体积值。当控制体积很小时，可以表示为，A是控制体积界面的面积。从而有

（1-63）

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

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

（1-64）

的线性插值结果为

（1-65）

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

（1-66）

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

（1-67）

式中，是常数，是随时间和物理量变化的项。将式（1-64）～式（1-67）代入方程（1-63），可得

（1-68）

整理后得

记为

（1-69）

式中，

（1-70）

对于一维问题，控制体积界面e和界面w处的面积和均为1，即单位面积，因此，式（1-70）中各系数可转化为

（1-71）

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

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

（1-72）

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

（1-73）

在此基础上，定义一维单元的佩克莱数（Peclet number）为

（1-74）

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

将式（1-72）、（1-73）代入方程（1-71），可得

（1-75）

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

（1-76）

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

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

将式（1-76）在一维计算网格上对时间及控制体积进行积分，可得

（1-77）

改写后，为

（1-78）

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

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

（1-79）

假定变量对时间的积分为

（1-80）

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

将、、及采用类似式（1-80）进行时间积分，式（1-79）可写成

（1-81）

整理后得

（1-82）

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

（1-83）

将式（1-83）代入方程（1-82），得

（1-84）

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

（1-85）

式中，

（1-86）

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

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

（1-87）

（1-88）

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

（1-89）

（1-90）

1.3.2　常用的离散格式与建立离散方程应遵循的原则

1. 常用的离散格式

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

（1-91）

（1-92）

（1-93）

式中，若为一阶迎风格式，；若为二阶迎风格式，。其中，系数和（高阶还有和）取决于所使用的离散格式，为了便于比较和编程计算，将不同离散格式下的系数aW和aE的计算公式总结于表1-2中。

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

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

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

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

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

（2）正系数原则

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

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

源项的斜率为负可以保证遵循正系数原则。从式（1-71）中看到，即使相邻节点的系数皆为正值，但只要有源项的存在，中心节点系数仍有可能为负。规定可以保证为正值。

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

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

1.3.3　离散方程的解法

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

1. 分离解法

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

1）假定初始压力场。

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

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

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

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

在使用分离求解器时，通常可以选择3种压强与速度的关联形式，即 SIMPLE、SIMPLEC 和 PISO。SIMPLE 和 SIMPLEC通常用于定常流动计算，PISO用于非定常流动计算，但是在网格畸变很大时也可以使用PISO格式。

2. 耦合解法

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

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

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

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

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

Fluent采用的离散格式包括一阶迎风格式、指数律格式、二阶迎风格式、QUICK格式、中心差分格式等。Fluent中各流场变量的迭代都由松弛因子控制，因此计算的稳定性与松弛因子紧密相关。在大多数情况下，可以不对松弛因子的默认设置进行修改，因为这些默认值是根据各种算法的特点得出的。而对某些复杂流动，默认设置可能不能满足稳定性要求，计算过程中可能出现振荡、发散等情况，此时需要适当减小松弛因子的值，以保证计算收敛。

在实际计算中，可以先用默认设置进行计算，如果发现残差曲线向上发展，则中断计算，适当调整松弛因子后再继续计算。在修改计算控制参数前，应该先保存当前计算结果，调整参数后，计算需要经过几步调整才能适应新的参数。

在计算发散时，可以考虑将压强、动量、湍流动能和湍流耗散率的松弛因子的默认值分别减小为0.2、0.5、0.5、0.5。在计算格式为SIMPLEC 时，通常没有必要减小松弛因子的值。