CFD 基 础（流体力学）

图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)