高等计算流体力学讲义(8)

第三章 不可压缩流动的数值方法

§1 基本方程及其性质

一、基本方程

    考虑不可压缩NS方程：

        （1）

        （2）

其中粘性应力为，

        （3）

如果粘性系数为常数，

        （4）

经无量纲化，常粘性系数不可压缩NS方程可以写为：

    ，

其中为运动粘性系数。NS方程也可以写为无量纲化形式

其中已经吸收到p中（代表）。不可压缩方程的边界条件为： 

固体壁面：，

进口条件：，

出口条件：。

不可压缩方程中的压力场可以相差任一常数而对速度场无影响，所以压力场只是在相差任意常数的条件下是确定的。为了确定全场压力值，还应指定流场中某一点的压力。

二、不可压N－S方程的特点：

（1）方程为二阶偏微分方程，二阶项中包含参数（粘性系数）。

边界层、分离、湍流…

（2）方程是非线性的，表现为对流项。

对一维问题，非线性项为。假定u的波数为k的Fourier分量为

                                                           （5）

则：。即振幅由；波数由。也就是说，振幅呈现非线性变化，且可以产生高频成分。粘性的作用，使得解的结构进一步复杂化，考虑模型方程

把（5）式带入模型方程，得

可见，雷诺数越大，或频率越低（流动结构的尺度越大），振幅衰减越慢。

综上所述：由于非线性的作用，会产生高频的流动结构；在大雷诺数的条件下，这些高频结构有较长的生命周期，并且与衰减缓慢的低频结构相互作用，使得流动表现出复杂的的非线性、多尺度特征。

（3）没有压力对时间的偏导数。

由于方程表现出很强的椭圆形，不能利用比较成熟的发展型偏微分方程的数值求解的理论和方法，造成数值求解的困难。尤其是没有明显的计算压力的方程，严格地说，必须耦合求解动量方程和连续方程才能求出压力场和速度场。这种做法，计算量是非常巨大的。目前，不可压缩流动的数值方法还远不如可压缩流动的数值方法更成熟。

三、耦合求解方案

NS方程的一种可能的求解方案是动量方程和连续方程完全耦合求解，例如可以采用下面的方法。

差分算子的定义为：

表示离散变量；i，j分别为x，y方向的单位向量；，分别为x，y方向的空间步长。其中对流项显式处理，用二阶精度的Adams-Bashforth格式离散：

    。

    然而，这种方法需要处理非常庞大、形状不规则、刚性很强的稀疏矩阵，计算量非常巨大且不易收敛。所以很少使用。实用的求解不可压缩流动的数值方法都需要把连续方程和动量方程在一定程度上进行解耦。这种解耦处理方法可以使求解效率显著提高，但也由此产生了一系列的新问题。

§2 MAC方法

（F. H. Harlow and J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8, 2182 (1965).）

一、基本方程    

考虑不可压缩NS方程，

                             （1）

              （2）

              （3）

在上面的方程中，压力场相当于一种约束条件，使得由（2）、（3）求得的速度场满足连续方程（1）。联立求解（1）～（3）可以求得速度和压力，但计算量很大。为了高效率的计算压力，MAC方法中首先要导出压力满足的方程，称为压力波松方程。对上述方程进行变换，整理后得：

        （4）

其中

考虑连续方程，有。但是，我们暂时保留与D有关的量。（4）式中左侧的第2、3项可以简化为：

    ，

这样，（4）式可以化为：

（5）

（5）式称为压力Poisson方程。经上述变换后，我们把N－S方程（1）、（2）、（3）式用（5）、（2）、（3）式代替。注意到：一般来说，（5）、（2）、（3）式与原始的N－S方程并不严格等价。也就是说，满足（5）、（2）、（3）式的解并不能严格保证。这就是为什么要在（4）、（5）式中保留与散度相关的量的原因。

二、求解步骤    

下面介绍MAC方法的求解步骤：

i）求解压力Poisson方程（已知，求）

（5）式的时间导数项用前差离散：

我们希望，所以取。数值实验表明，这种处理方法可以使散度的绝对值保持在非常接近于零的水平，从而不影响数值解的有效性。此时（5）式写为：

        （6）

对于二维问题，其中的Laplace算子可以用标准的5点中心格式离散，即

    。

（6）式中的其他项用中心差分离散。（6）式称为离散压力泊松方程，可以用求解泊松方程的常规数值方法求解。可知，一旦我们知道了n时刻的速度，通过（6）式就可以求出n时刻的压力。

ii）求解动量方程

在求出压力后，我们通过求解动量方程计算n＋1时刻的速度：

        （7）

        （8）

一般，压力和粘性项用中心差分离散，而对流项则可以用中心差分或者迎风差分离散。

三、边界条件

固壁边界条件：速度满足

为了求解压力波松方程，还应补充压力的边界条件，在固壁处，压力边界条件由（6）式简化，得：

    ，

在求解Poisson方程时，我们使用上式的法向分量作为边界条件，即

图1 交错网格

四、空间离散和交错网格

后面我们还将介绍，在求解不可压缩N－S方程时，当压力和速度的存储位置相同时（称为非交错网格），压力场有可能出现棋盘状非物理振荡，这种现象称为“奇偶失连”现象。为了避免非物理振荡，MAC方法在交错网格上求解N－S方程。图（1）给出了物理量的存储位置。在数值计算中，压力泊松方程（（6）式）在(i,j)点求解；x方向动量方程（（7）式）在()求解；y方向动量方程（（8）式）在（）处求解。下面，我们以下x-方向动量方程为例，简述空间离散的方法。此时，（7）式具体化为：

（9）

其中和分别为一阶和二阶的中心差分算子，即：

    。

可以为中心差分或者迎风差分算子，例如：

中心差分： 

一阶迎风： 

二阶迎风： 

注意，采用交错网格时，有些计算中需要的物理量在网格点上没有定义。例如，在（7）式时，需要用到等。这些物理量要通过插值确定，如在均匀网格上，

另外，在边界上，常通过构造虚拟网格处理边界条件，具体思路可参见图2。如果边界为固壁，则我们取虚拟网格点上。

    MAC方法通过压力波松方程，构造了计算压力场的方法。这个方法实际上是一种压力和速度之间的解耦算法。采用这一方法，计算效率比联立求解连续和动量方程高，但是也存在一些问题。最主要的是：连续方程不能完全满足。通过在压力波松方程中保留项，可以使速度场散度保持较小的水平，但不能从本质上解决问题。所以，MAC方法的计算精度不是很高。

    MAC方法（Marker and Cell）实际上是一种处理有自由面的流动的方法，因为其中涉及不可压缩NS方程的求解，所以我们也用MAC方法称呼相应的N－S方程求解方法。在其他场合，MAC方法通常指处理自由面问题的方法。

§3 SIMPLE方法

图3 交错网格的控制体及变量的存储位置

一、 SIMPLE方法

SIMPLE方法是基于守恒形方程的有限体积方法，守恒型的N－S方程为：

                                               （1）

                     （2）

                     （3） 

前述MAC方法也可采用守恒型方程和有限体积方法离散，请参照下述SIMPLE方法自行把MAC方法改写成有限体积形式。

SIMPLE方法与MAC方法一样，在交错网格中离散N－S方程。图3显示了交错网格的控制体及变量的位置。变量指标的编码方法参照图1。

在主控制体上积分连续方程（1）式，有：

                                        （4）    

把x方向的动量方程改写为

        （5）

其中

        （6）

        （7）

在U动量控制体(i+1/2,j)上积分（4）式,得

        （8）

（8）式中E、F的具体算法见下一小节。为了保证计算的稳定性，（8）式在时间方向常采用隐式离散方案，如Euler隐式格式：

        （9）

把E、F的表达式代入（9）式，整理得，

       （10）

其中系数等为由n时刻物理量决定的系数，是数值离散中的已知量（如体积力）的合并形式，称为源项。（约定，无上标的物理量u，p均为n＋1时刻的值。）

    在V控制体上积分y方向的动量方程，同样可以得到y方向速度满足的类似（10）的关系：

        （11）

在已知精确压力时，用迭代法求解（10）（11）式可以得到n＋1时刻的速度u，v。由于n＋1时刻的压力未知，我们首先给出压力的估计值，

代入（10）（11），记n＋1时刻的速度估计值为，有：

        （12）

        （13）

求解（12）（13）得。设n+1时刻的数值解与估计值之间满足下列关系：

    ，    （14）

其中为压力和速度的校正量。易知

        （15）

        （16）

最初的SIMPLE方法主要用于求解定常问题，其中的一个关键步骤是略去（15）（16）式右侧第一个括号中的各项（为什么可以这样化简？），从而有

        （17）

（17）称为速度修正方程。把（17）和同理得到的的速度修正方程代入（4）式，可以得到压力修正量满足的方程：

        （18）

（18）式称为压力修正方程。用迭代法求解（18）式，可得压力修正量。从而可根据（17）（14）式计算出n＋1时刻的。

    从上面的叙述可以看出，SIMPLE方法的基本思想是：通过估计的压力场计算速度场的估计值；通过把离散的动量方程带入离散的连续方程计算速度和压力的校正量。SIMPLE方法和MAC方法一样，也是一种速度、压力的解耦算法。上面所述方法，由于引入了（15）、（16）式到（17）式的近似，严格地说只适用于定常流动的计算。为了计算非定常流动，一般需要引入内迭代，即：从n时刻的流场出发，执行一次上述求解过程，把计算得到的当作n＋1时刻的第一次内迭代值，记为。然后把当作n＋1时刻压力的改进的估计值，从n时刻的流场出发，再执行一次上述求解过程，把把计算得到的当作n＋1时刻的第二次内迭代值，记为。重复这个过程k次，直至为止。其中是我们给定的压力修正量的收敛指标。

二、通量E，F的计算方法

以x-方向的动量方程为例，通量E可以写为

        （19）

其中为对流项，为扩散项，p是压力。对流项是非线性项，为了使用隐式格式，对流项要进行线性化处理，即

对流项的离散有多种方案，如：

中心差分

        （20）

一阶迎风

        （21）

二阶迎风

        （22）

扩散项用中心插分离散，

        （23）

在交错网格上，压力项可以直接取主控制体中心点的值。

F也可以写为对流项和扩散项的和，即

        （24）

当采用一阶迎风格式离散时，

        （25）

当采用二阶迎风格式离散时，

扩散项的离散公式为：

        （26）

注意到，当对流项采用二阶迎风格式或者其他高阶格式离散时，（i+1/2,j）处的动量方程包括诸如之类的项，因而不能直接写为（10）式的形式。这是的处理方法是：在计算定常流动时，这些项可取n时刻的值，合并到中；在计算非定常流动时，这些项取前一次内迭代的值，同样合并到中。这种处理方法，可以使离散方程在对流项的各种离散方案中，均可写成（10）式的形式。

SIMPLE方法有很多变形和改进形式，如SIMPLEC，SIMPLEST等。

§4 非交错网格上的SIMPLE方法，动量内插技术

一、动量内插方法

基于交错网格的计算方法如MAC和SIMPLE等的优点是可以有效的避免奇偶失连的问题，其缺点是编程复杂，且比较难于推广到任意形状控制体的情况。为了克服这些缺点，C. M. Rhie 和W. L. Chow （AIAA J. 21, 1525 ，1983)发展了基于动量内插技术的非交错网格上的有限体积方法。与交错网格不同，非交错网格只有一套控制体，所有物理量均存储在网格中心。参见图4。

考虑方程

                                               （1）

                     （2）

                     （3） 

在控制体上（i，j）积分连续方程（1）式，有：

                                   （4）

在同样控制体上积分x方向动量方程，有：

        （5）

对E，F采用与SIMPLE方法类似的差分方法，整理后有

        （6）

其中。对y方向动量方程积分，类似地有

        （7）

在已知压力时，用迭代法求解（6）（7）式可以得到n＋1时刻的速度u，v。由于n＋1时刻的压力未知，我们首先给出压力的估计值，

代入（6）（7），记n＋1时刻的速度估计值为，有：

        （8）

        （9）

求解（8）（9）得，其中的压力项采用中心差分离散。设

        （10）

易知

        （11）

        （12）

同交错网格上的SIMPLE方法一样，略去（11）（12）式右侧第一个括号中的各项，从而有

        （13）

在计算出的情况下，由（13）式可算得。下面介绍的算法。

    交错网格上SIMPLE方法通过把速度修正方程代入连续方程得到压力修正方程，从而求出。在这个步骤中，利用了交错网格中主控制体的边界分别是U动量控制体和V动量控制体中心这一特点。在非交错网格中，只有一套控制体，为了把动量方程代入连续方程，必须进行某种插值处理。其中最直接的方法是：

    。    （14）

把（13）、（14）式代入（4）式，可得压力修正量的Poisson方程。然而，通过这种方法得到的压力场具有明显的奇偶失连现象。 我们下面将要介绍的动量插值技术是一种可以有效的避免奇偶失连问题的方法。在（i，j）点，（6）式可以具体化为：

        （14）

其中

同理,在(i+1,j)点,

        （15）

其中

参考（4）式,其中的可由插值决定:

        （16）

动量插值的关键步骤是把（16）式近似为

        （17）

对于估计的压力值,速度的估计值满足

        （18）

采用与SIMPLE方法同样的近似,可得(i,j)控制体的(i+1/2,j)边界处的速度修正方程:

        （19）

同理推导(i-1/2,j),(i,j+1/2),(i,j-1/2)边界处的速度修正方程,并代入（4）式,可得(i,j)处的压力修正方程:

        （20）

总结整个求解步骤。

二、非交错网格奇偶失连现象初步分析

为什么在非交错网格上会出现奇偶失连的现象呢？我们给出一个直观的解释。假定在如下图的网格上，n时刻流动速度为零，压力分布如图所示：

并由此得到压力修正方程。通过这种方法求解N－S方程，我们可以验证：在n＋1时刻速度仍为零，这显然与物理过程不符。这种解法中，压力和速度没有实现正确的耦合，从而造成速度场对上述“棋盘格”状的压力分布不敏感，从而可能造成奇偶失连现象。

    S. W. Armfield (Computers and Fluids, 20:1-17, 1991.)认为，离散的压力修正方程在非交错网格上不具有椭圆性，是产生奇偶失连现象的主要原因。他定义了一种所谓椭圆测度，记为。当时，称数值方法具有强椭圆性；当时，称数值方法具有弱椭圆性；当时，数值方法不具有椭圆性。按照这一理论，我们可以证明，在交错网格上的SIMPLE方法；非交错网格上，压力修正方程由（13）、（14）式得到时，；而采用动量插值的方法推导压力修正方程时，在一定假设下，可以证明。大量的研究表明，椭圆测度与奇偶失连现象有密切关系。一般而言，时，只要很小的就可以避免奇偶失连问题；但当雷诺数很高或者流场存在奇异性时，则需要较大的才能克服奇偶失连问题；当时，奇偶失连现象一般不会发生。

§5投影方法(Projection Method)

    投影方法亦称格式(Splitting Scheme)，或分数步方法(Fractional Step Method)，是一种基于原始变量公式来计算粘性不可压缩流动的有效方法。它最早由Chorin在六十年代末提出来，类似的思想也出现在Temam的文章中。Chorin在他的一系列文章中，对投影法进行了详尽地阐述，并认为投影法在时间方向上至多只有一阶精度。然而二十多年的实践证明，投影法的实际应用比Chorin所预期的要好得多，出现了多样化的投影形式和一些高阶精度的格式。投影方法与SIMPLE方法由很多类似之处，但是计算非定常问题的效率要高于SIMPLE方法。最近几年，由于湍流DNS和LES的需要，投影方法受到了高度重视。

一、投影法的基本原理及其求解步骤

考虑一个n维有界域，不可压缩流的Navier-Stokes方程可以写成：

                             (1)

                                             (2)

边界条件：

                                           (3)

投影方法，与MAC方法和SIMPLE方法一样，也是一种速度和压力的解耦算法，其根据是Hodge分解定理：任意的向量（场）均可唯一的分解为一个散度为零的向量和一个旋度为零的向量的和，即：

                                         (4)

其中满足。仿照矢量分析的叫法，我们把这个分解过程叫做投影。设为向量投影到散度为零平面的投影算子，则有。

因此，由投影公式，我们可以把(1.1)式写成：

式中已不含压力项，可以用数值方法求解。

   下面介绍Chorin构造的投影格式，先将(1)式和(2)式写成时间离散的半离散形式：

                   (5)

                                           (6)

其中的空间导数项可采用适当的方法离散，如对流项采用迎风或中心差分，而粘性项和散度项采用中心格式。边界条件为：

                                       （7）

Chorin采用算子的方法，把（6）式分为两个算子的和：

这样，就引入了一个中间速度场。投影方法的第一步（预测步）是求解中间速度场：

                       (9)

中间速度场不必满足散度为零的条件。严格地说，其边界条件是未知的，我们必须人为地构造的近似边界条件，Chorin取：

    。                                     (10)

投影方法的第二步称为投影过程。考虑（8b）式，或者由Hodge分解定理，可以表示为：

                                    (11)

其中

                                         （12） 

显然，这一步是一个标准的投影过程：向量被分解为散度为零的向量和旋度为零的向量的和。由（11）和（12）式，得：

                                             (13)

即投影过程和一个压力波松方程相对应。投影过程的边界条件可由（7）、（10）式和（11）式提供：

                                          (14)

由(12)和(14)式计算得到，然后将和代入到(11)式中得到。

二、二阶精度投影方法

Chorin的投影格式在时间方向上只有一阶精度。从投影方法的实施过程中我们可以看出：投影方法人为地对动量方程的离散形式（5）进行算子，引入了中间速度场，从而需要给出的边界条件。显然，的边界条件对格式的精度有直接影响。由于采用了算子和人工边界条件，投影方法在时间方向的精度难以提高，而在空间方向达到高精度相对容易。所以，我们主要讨论改进时间精度的方法，以构造出时间方向二阶精度的投影格式。和Chorin的做法类似，我们首先给出NS方程的二阶半离散格式：

         (15a)

                                   （15b）

边界条件：

式中：表示二阶对流项在时间处的二阶近似项，通常用Adams－Barshforth方法显式计算，即：

    。

我们把（15a）式写成算子的形式：

是n＋1/2时刻压力的估计值。把（16b）带入（16a），并与（15a）比较，可得：

这样，投影方法的实施可以分为以下三步：

Step1求解中间速度场：

                  (17)

                                        (18)

（18）式为的边界条件，其具体形式待定。

Step2 投影：

                                    (19)

。                                          

的边界条件为：

    ，

其形式待定。在求得后，n+1时刻的速度由（19）式确定。                  

Step3 计算压力：

。                       (20)

前面已经提及，中间速度场的边界条件对投影方法的精度有直接的影响。另一方面，如果是的一个较好近似，则也将是的良好近似，这样，中间速度场的边界条件为（10）式即可获得较好效果。因此，二阶精度的投影格式可以通过选取合适的或者获得。Brown等（Journal of Computational Physics 168, 4–499 (2001)）推荐了两种二阶精度的投影格式。

① 

其中分别为法向量和切向量。

②

    如果上述投影过程逼近（15）到至少二阶精度，则这些投影格式也具有二阶精度。但是，直接证明这一点并不容易。Brown等人采用正则模态分析（Normal Mode Analysis）的方法证明了这两种投影格式具有二阶精度。正则模态分析方法在分析不可压NS方程的精度和性能中有广泛的应用，但不应看作是严格的证明。最近，Shen Jie等证明，对于一般的区域，这种投影格式速度可以达到二阶精度，而压力的精度为3/2阶。最近，我们导出了投影方法具有时间方向二阶和三阶精度的条件，从理论上严格证明了投影方法在适当的人工边界条件下，可以达到高阶精度。

§6 人工压缩方法

人工压缩方法也叫拟压缩性方法，其基本思想是在连续方程中引入与压力有关的人工压缩项，使修改后的方程与可压缩流动方程具有相同的类型。这样，求解可压缩流动的数值方法就可以被用来计算不可压缩问题。

考虑可压缩的连续方程

        （1）

流体的压缩性一般用压缩性系数来衡量。

                                                              （2）

对于等熵完全气体，此时。由（2）式，（1）式可以改写为：

        （3）

音速为，所以（3）式也可写为

        （4）

对于不可压缩流动，，所以连续方程简化为。人工压缩方法假定对不可压流动，（4）式中的音速为人为给定的有限量，定义，（4）式可写成

        （5）

其中是人工压缩系数。显然，上式只适用于定常流动，计算收敛时，对定常流动有，从而连续方程可以得到满足。人工压缩方法也可以发展到计算非定常不可压流动，但要同时采用双时间步的时间推进方法。

    采用拟压缩性近似，并用代替，流动的支配方程为：

        （6）

（6）式在形式上和数学性质上与可压缩的N－S方程类似，可以应用计算可压缩流动的数值方法求解。

将不可压守恒型方程组写成矩阵形式：

其中：

，， ， 

人工压缩方法实际上是对矩阵进行预处理，作变换，得：

其中

。

处理过的方程能够利用时间推进的方法进行计算。

令

则矩阵A、B的特征值分别为：

如果取水中的声速，显然每个矩阵的条件数都很大，因此偏微分方程组是刚性的。所以实际应用中的值通过试算来确定。