1.4 流场数值计算算法分析

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

1．分离解法

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

（1）假定初始压力场。

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

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

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

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

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

2．耦合解法

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

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

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

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

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

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-7中所示的均匀网格是向前错位的，因为u的速度u的I位置到标量节点（I，J）的距离是－1/2；同样，v速度v的j位置到标量的距离是－1/2。当然，也可以选用向后错位的速度网格。

在使用了上述交错网格后，生成离散方程的方法和过程与原来的基于普通网格的方法和过程是一样的，只是需要注意所使用的控制体积有所变化。在交错网格中，由于所有标量（如压力、温度、密度等）仍然在主控制体积上存储，因此，以这些标量为因变量的输运方程的离散过程及离散结果与前面的一样，主要是在交错网格中生成的u和v两个动量方程的离散方程时，积分用的控制体积不再是原来的主控制体积，而是u和v各自的控制体积，同时压力梯度项从源项中分离出来。例如，对u控制体积，该项积分为。从而，对于方程（1-104）的离散方程，对u方向的动量方程使用u控制体积，可以写出在位置（i，J）处的关于速度ui，J的动量方程的离散形式：

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

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

式（1-131）中的S是对源项S线性化分解为的结果，若不随速度u而变化，则有式（1-131）中的V是u控制体积的体积。式（1-129）中的压力梯度项已经按线性插值的方式进行了离散，线性插值时使用了u控制体积边界上的两个节点的压力差。

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

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

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

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

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

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

建立式（1-135）所使用的v控制体积表示在图1-9中。

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

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

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

2．SIMPLE算法

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

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

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

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

1）速度修正方程

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

根据动量方程的离散方程（1-129）和方程（1-135），有

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

同样地，我们定义速度修正值u′和v′，以联系正确的速度场（u，v）与猜测的速度场（u*，v*），有

将正确的压力场p代入动量离散方程（1-129）与方程（1-135），得到正确的速度场（u，v）。现在，从方程（1-129）和方程（1-135）中减去方程（1-137）和方程（1-138），并假定源项b不变，有：

引入压力修正值与速度修正值的表达式（1-139）、式（1-140）和式（1-141），方程（1-142）和方程（1-143）可写成

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

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

其中，

将式（1-146）和式（1-147）所描述的速度修正值，代入式（1-140）和式（1-141），有

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

其中，

式（1-149）～式（1-153）表明，如果已知压力修正值p′，便可对猜测的速度场（u*，v*）作出相应的速度修正，得到正确的速度场（u，v）。

2）压力修正方程

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

针对图1-10所示的标量控制体积，连续方程（1-154）满足如下离散形式：

将正确的速度值，即式（1-149）～式（1-153），代入连续方程的离散方程（1-155），有

整理后，得

该式可简记为

其中，

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

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

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

3）SIMPLE算法的计算步骤

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

1.4.2 其他算法介绍

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

1．SIMPLER算法

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

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

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

将离散后的动量方程（1-129）和（1-135）重新改写后，有

在SIMPLER算法中，定义伪速度与如下：

这样，式（1-160）与式（1-161）可写为

以上两式中的系数d，仍沿用1.4.1节所给出的计算公式。同样可写出ui＋1，J与vI，j＋1的表达式。然后，将ui，J、vI，j、ui＋1，J与vI，j＋1的表达式代入离散后的连续方程（1-155），有

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

其中，

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

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

2．SIMPLEC算法

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

在SIMPLE算法中，为求解的方便，略去了速度修正方程中的项，从而把速度的修正完全归结为由压差项的直接作用。这一做法虽然并不影响收敛解的值，但加重了修正值p′的负担，使得整个速度场迭代收敛速度降低。实际上，当在略去时，犯了一个“不协调一致”的错误。为了能略去而同时又能使方程基本协调，试在方程（1-144）的等号两端同时减去，有

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

其中，

类似地，有

其中，

将式（1-172）和式（1-173）代入SIMPLE算法中的式（1-149）和式（1-150），得到修正后的速度计算式：

式（1-174）和式（1-175）在形式上与式（1-149）和式（1-150）一致，只是其中的系数项d的计算公式不同，现在需要按式（1-171）和式（1-173）计算。

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

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

3．PISO算法

PISO是英文pressure implicit with spltting 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-137）与（1-138），得到速度分量u*与v*。

2）第一步修正

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

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

就像在SIMPLE算法中一样，将式（1-179）与式（1-180）代入连续方程（1-155），产生与式（1-157）具有相同系数和源项的压力修正方程。求解该方程，产生第一个压力修正值p′。一旦压力修正值已知，可通过方程（1-179）与式（1-180）获得速度分量u**与v**。

3）第二步修正

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

注意这两式实际就是式（1-137）和式（1-138）。为引用方便，给出新的编号。

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

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

现在，从式（1-183）中减去式（1-181），从式（1-184）中减去式（1-182），有

其中，记号p″是压力的二次修正值。有了该记号，p***可表示为

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

其中，ai，j＝ai＋1，j＋ai－1，j＋ai，j＋1＋ai，j－1。读者可参考建立方程（1-158）同样的过程，写出各系数如下：

下面对源项b′为何是式（1-189e）的形式，作一简要分析和解释。

对比建立方程（1-158）的过程，可以发现，式（1-189e）中的各项，是因在u***和v***的表达式（1-185）和式（1-186）中存在和项所导致的，而在u和v的表达式（1-149）和式（1-150）中没有这样的项，因此，式（1-158）不存在类似式（1-189e）中的各项。但式（1-158）存在另外一个源项，即[（ρu*A）i，J－（ρu*A）i＋1，J＋（ρv*A）I，j－（ρv*A）I，j＋1]，这是因速度u和v的表达式（1-149）和式（1-150）中的u*与v*项所导致的。按此推断，在式（1-189e）中也应该存在类似表达式[（ρu**A）i，J－（ρu**A）i＋1，J＋（ρv**A）I，j－（ρv**A）I，j＋1]。但是，由于u**和v**满足连续方程，因此[（ρu**A）i，J－（ρu**A）i＋1，J＋（ρv**A）I，j－（ρv**A）I，j＋1]为0。

现在，求解方程（1-188），就可得到二次压力修正值p″。这样，通过下式就可得到二次修正的压力场：

最后，求解方程（1-185）与式（1-186），得到二次修正的速度场。

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

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

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的效果可能更好些。

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