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

1高等计算流体力学讲义（8）第三章不可压缩流动的数值方法§1基本方程及其性质一、基本方程考虑不可压缩NS方程：0u（1）()ptuuuf（2）其中粘性应力为，2S（3）12()TSuu如果粘性系数为常数，u（4）经无量纲化，常粘性系数不可压缩NS方程可以写为：0()ptuuuufu，其中/为运动粘性系数。NS方程也可以写为无量纲化形式01()Reptuuuufu其中已经吸收到p中（p代表/p）。不可压缩方程的边界条件为：固体壁面：walluu，进口条件：inuu，2出口条件：nu0。不可压缩方程中的压力场可以相差任一常数而对速度场无影响，所以压力场只是在相差任意常数的条件下是确定的。为了确定全场压力值，还应指定流场中某一点的压力。二、不可压N－S方程的特点：（1）方程为二阶偏微分方程，二阶项中包含参数（粘性系数）。边界层、分离、湍流…（2）方程是非线性的，表现为对流项()uu。对一维问题，非线性项为uux。假定u的波数为k的Fourier分量为()sinuutkx（5）则：21sin22uuukxx。即振幅由212uu；波数由2kk。也就是说，振幅呈现非线性变化，且可以产生高频成分。粘性的作用，使得解的结构进一步复杂化，考虑模型方程221Reuutx把（5）式带入模型方程，得2(/Re)()ktute可见，雷诺数越大，或频率越低（流动结构的尺度越大），振幅衰减越慢。综上所述：由于非线性的作用，会产生高频的流动结构；在大雷诺数的条件下，这些高频结构有较长的生命周期，并且与衰减缓慢的低频结构相互作用，使得流动表现出复杂的的非线性、多尺度特征。（3）没有压力对时间的偏导数。由于方程表现出很强的椭圆形，不能利用比较成熟的发展型偏微分方程的数值求解的理论和方法，造成数值求解的困难。尤其是没有明显的计算压力的方程，严格地说，必须耦合求解动量方程和连续方程才能求出压力场和速度场。这种做法，计算量是非常巨大的。目前，不可压缩流动的数值方法还远不如可压缩流动的数值方法更成熟。三、耦合求解方案NS方程的一种可能的求解方案是动量方程和连续方程完全耦合求解，例如可以采用下面的方法。1,,1/21/211/2,,,,,()[]2nnijijnnnnnijijhijijijpLtuuDuuGuuf31,0nijDu差分算子的定义为：,()[()][()]ijxyDDDij,()()()ijxyDDGij1/2,1/2,,()()()ijijxijDx,1/2,1/2,()()()ijijyijDy1,,1/2,()()()2ijijij,1,,1/2()()()2ijijij1,,1,,1,,1,22()2()()()2()()()ijijijijijijhijLxy()表示离散变量；i，j分别为x，y方向的单位向量；x，y分别为x，y方向的空间步长。其中对流项显式处理，用二阶精度的Adams-Bashforth格式离散：1/21,,,()1.5()0.5()nnnijijijDuuDuuDuu。然而，这种方法需要处理非常庞大、形状不规则、刚性很强的稀疏矩阵，计算量非常巨大且不易收敛。所以很少使用。实用的求解不可压缩流动的数值方法都需要把连续方程和动量方程在一定程度上进行解耦。这种解耦处理方法可以使求解效率显著提高，但也由此产生了一系列的新问题。§2MAC方法（F.H.HarlowandJ.E.Welch,Numericalcalculationoftime-dependentviscousincompressibleflowoffluidwithfreesurface,Phys.Fluids8,2182(1965).）一、基本方程考虑不可压缩NS方程，0uvxy（1）422221()Reuuupuuuvtxyxxy（2）22221()Revvvpvvuvtxyyxy（3）在上面的方程中，压力场相当于一种约束条件，使得由（2）、（3）求得的速度场满足连续方程（1）。联立求解（1）～（3）可以求得速度和压力，但计算量很大。为了高效率的计算压力，MAC方法中首先要导出压力满足的方程，称为压力波松方程。对上述方程进行变换(2)(3)xy，整理后得：222222221()()()()ReDuuvvppDDuvuvtxxyyxyxyxy（4）其中uvDxy考虑连续方程，有0D。但是，我们暂时保留与D有关的量。（4）式中左侧的第2、3项可以简化为：22()2()()uuvvDDuvxyxyxy，这样，（4）式可以化为：221()2()()ReuuvvDDDpuvDxyxyxyt（5）（5）式称为压力Poisson方程。经上述变换后，我们把N－S方程（1）、（2）、（3）式用（5）、（2）、（3）式代替。注意到：一般来说，（5）、（2）、（3）式与原始的N－S方程并不严格等价。也就是说，满足（5）、（2）、（3）式的解并不能严格保证0D。这就是为什么要在（4）、（5）式中保留与散度相关的量的原因。二、求解步骤下面介绍MAC方法的求解步骤：i）求解压力Poisson方程（已知,nnuv，求np）（5）式的时间导数项用前差离散：1nnDDDtt我们希望1nnDD，所以取10nD。数值实验表明，这种处理方法可以使散度D的绝5对值保持在非常接近于零的水平，从而不影响数值解的有效性。此时（5）式写为：221()2()()RennnnnnnnnnnuuvvDDDpuvDxyxyxyt（6）对于二维问题，其中的Laplace算子可以用标准的5点中心格式离散，即1,,1,,1,,1,2222nnnnnnijijijijijijnijpppppppxx。（6）式中的其他项用中心差分离散。（6）式称为离散压力泊松方程，可以用求解泊松方程的常规数值方法求解。可知，一旦我们知道了n时刻的速度nu，通过（6）式就可以求出n时刻的压力np。ii）求解动量方程在求出压力后，我们通过求解动量方程计算n＋1时刻的速度1nu：221221()Rennnnnnnnnuupuuuutuvxyxxy（7）221221()Rennnnnnnnnvvpvvvvtuvxyyxy（8）一般，压力和粘性项用中心差分离散，而对流项则可以用中心差分或者迎风差分离散。三、边界条件固壁边界条件：速度满足()wallwallifuuu0u0为了求解压力波松方程，还应补充压力的边界条件，在固壁处，压力边界条件由（6）式简化，得：1Repu，在求解Poisson方程时，我们使用上式的法向分量作为边界条件，即1Repnnu6图1交错网格四、空间离散和交错网格后面我们还将介绍，在求解不可压缩N－S方程时，当压力和速度的存储位置相同时（称为非交错网格），压力场有可能出现棋盘状非物理振荡，这种现象称为“奇偶失连”现象。为了避免非物理振荡，MAC方法在交错网格上求解N－S方程。图（1）给出了物理量的存储位置。在数值计算中，压力泊松方程（（6）式）在(i,j)点求解；x方向动量方程（（7）式）在(1/2,ij)求解；y方向动量方程（（8）式）在（,1/2ij）处求解。下面，我们以下x-方向动量方程为例，简述空间离散的方法。此时，（7）式具体化为：1/2,1/2,1/2,1/2,1/2,11/2,1/2,221/2,1/2,22()()1()Re()()xnynxnijijijnnijijnnijijxnynijijDuDupuvxyxuutuuxy（9）其中1/2,xij和21/2,xij分别为一阶和二阶的中心差分算子，即：1/2,1,,21/2,3/2,1/2,1/2,21/2,1/2,11/2,1/2,122xnnnijijijxnnnnijijijijynnnnijijijijpppuuuuuuuu。,xyDD可以为中心差分或者迎风差分算子，例如：中心差分：1/2,3/2,1/2,1()()2xnnnijijijDuuu1/2,iju1/2,iju,1/2ijv,1/2ijv,ijp7一阶迎风：1/2,1/2,1/2,1/2,3/2,1/2,0()nnnijijijxnijnnijijuuifuDuuuotherwise二阶迎风：1/2,1/2,3/2,1/2,1/2,1/2,3/2,5/2,1.520.50()1.520.5nnnnijijijijxnijnnnijijijuuuifuDuuuuotherwise注意，采用交错网格时，有些计算中需要的物理量在网格点上没有定义。例如，在（7）式时，需要用到1/2,ijv等。这些物理量要通过插值确定，如在均匀网格上，1/2,,1/2,1/21,1/21,1/21()4ijijijijijvvvvv另外，在边界上，常通过构造虚拟网格处理边界条件，具体思路可参见图2。如果边界为固壁，则我们取虚拟网格点上1/2,11/2,ijijuu。MAC方法通过压力波松方程，构造了计算压力场的方法。这个方法实际上是一种压力和速度之间的解耦算法。采用这一方法，计算效率比联立求解连续和动量方程高，但是也存在一些问题。最主要的是：连续方程不能完全满足。通过在压力波松方程中保留uvDxy项，可以使速度场散度保持较小的水平，但不能从本质上解决问题。所以，MAC方法的计算精度不是很高。MAC方法（MarkerandCell）实际上是一种处理有自由面的流动的方法，因为其中涉及不可压缩NS方程的求解，所以我们也用MAC方法称呼相应的N－S方程求解方法。在其他场合，MAC方法通常指处理自由面问题的方法。§3SIMPLE方法图2交错网格下的边界处理v=0v-uu（i，j）8图3交错网格的控制体及变量的存储位置一、SIMPLE方法SIMPLE方法是基于守恒形方程的有限体积方法，守恒型的N－S方程为：0uvxy（1）222221()Reuuuvpuutxyxxy（2）222221()Revuvvpvvtxyyxy（3）前述MAC方法也可采用守恒型方程和有限体积方法离散，请参照下述SIMPLE方法自行把MAC方法改写成有限体积形式。SIMPLE方法与MAC方法一样，在交错网格中离散N－S方程。图3显示了交错网格的控制体及变量的位置。变量指标的编码方法参照图1。在主控制体上积分连续方程（1）式，有：11111/2,1/2,,1/2,1/2()()0nnnnijijijijuuyvvx（4）把x方向的动量方程改写为0uEFtxy（5）其中21ReuEupx（6）主控制体u控制体v控制体(i,j)(i,j-1)(i+1/2,j)(i+1,j)(i-1/2,j)(i,j+1)(i,j+1/2)(i,j-1/2)(i-1,j)puv