1、The Elements of Computational Fluid Dynamics第七章 不可压缩流动的数值方法7.1 基本方程7.2 涡量-流函数方法7.3 SIMPLE方法7.1 基本方程0()/01()Re-vvvvfpvtppuuuufput 运动黏性系数，代表。压力场可以相差任一常数而对速度场无影响，为了确定全场压力值，应指定流场中某一点的压力。：不包含压力的时间导数项单一，表现出椭圆 抛物组合型的特介质不可压缩流动的控制方程无量纲化形式控制方程的特点点。不存在状态方程，压力不具有热力学意义；压力对速度场加以限制，使连续性方程得到满足。压力场的波动具有无穷大的传播速度，瞬间传遍
2、全场，以使不可压缩条件在任何时间、任何地点满足。1,1/21/211/2,1,Navier-Stokes-Navier-Stokes1()2Re0nni ji jnnnnni ji jhi ji ji jni juuD uuGpL uuftDu目前，不可压缩流动的数值方法不如可压缩流动的数值方法成熟。动量求解不可压缩方程的困难在于处理“压力 速度”耦合问题。不可压缩方程的一种可能的求解方案：方程和连续方程完全耦合求解。其中，,1/2,1/2,1/2,1/2,1,11/2,1/21,(*)(*)(*);(*)(*)(*)(*)(*)(*)(*)(*);(*)(*)(*)(*)(*)(*);(*)
3、22(*)2(*)(*)(*)i jxyi jxyijiji ji jxi jyi jiji ji ji jiji jiji jihi jDD iDjGDiDjDDxyL1,1,1221/21,(*)2(*)(*)Adams-Bashforth()1.5()0.5()ji ji ji jnnni ji ji jxyD uuD uuD uu，如：用时间方向具有二阶精度的如果对流项采用显散式处理格式离-得到一个非常庞大、形状不规则、刚性很强的稀疏线性方程组。计算量非常大，不易收敛。，得到一个非线性方程组，其求解更如果对流项采用隐式处理因此，完全耦合求解方法在实用中很少使用。实用中，为了处理压力 速
4、度耦合问题，需要把连续方程和动量方程在一定程度上进为困难。行解耦。7.2 涡量流函数方法Navier-StokesNavier-Stokes计算二维不可压缩方程的有效算法。基本思想：对方程进行变换，写成涡量流函数的形式。7.2.1 基本方程222222220(7.2.1)1(7.2.2)Re1(7.2.3)Reuvxyuuupuuuvtxyxxyvvvpvvuvtxyyxy 22(7.2.2)(7.2.3)1(7.2.4)Re,(7.2.5)(7.2.6)(7.2.4)(7.2.6)(7.2.5)xytxyvuyxuvvuxyx 引入涡量方程，消去方程中的压力项引入流函数，使连续方程自然满涡量
5、对 求偏导，对 求偏导，推导出涡量满足的方程为满足：通过计算涡量；通过计算流函数；通过计算速度分量。如果需要计算压求解过程力，把足222Poisson2Poissonxyyuu vvpxyxy 方向的动量方程对 求导数，把 方向的动量方程对 求导数，二者求和后，利用连续性方程，可以得到压力的方程，求解压力的方程，可以得到压力的分布。7.2.2 差分格式21,1,1,1,1,1,1,1,122,1(7.2.4)Re1.FTCS22221Re()()2.-txynnnnnni ji jijiji ji jnni ji jnnnnnniji jiji ji ji jiuvuvtxyxy数值求解涡量满
6、足的方程（对流 扩散方程如果采用格式如果采用迎风格式）：122,2222,2222122Re()()()()()()(),()22nnnnji ji ji jnnxxxi jxi jnnnni jyi jyxi jyi jnnyi jyi juutxxvvyyxy 其中，当01时，一阶迎风格式；当时，二阶迎风格式。21,11,222121,22223.Crank-Nicolson1(7.2.4)Re4412Re()()()()Poisstxynnnni ji ji ji jnnnnxi jxi jyi jyi jnnnnxi jyi jxi jyi juvuvtxxxyxy 数值求解流函如果数
7、满足的采用格式21111111,1,1,11,22(7.2.6)2o2()()nnnnnnniji jiji ji ji jni jxy 二阶精度的差分格式为方程7.2.3 边界条件7.1图：二维槽道流动的计算域(1)0,0(2)0()()()0ywallininletininletininletjjMuiuuyuuyvvy固体壁面边界：边界条件为进口边界：边界条件为分量形式（假定进口流动方向是水平的）(3)xiMx出口边界：当槽道的长度远大于其高度时，边界条件可以认为是充分发展条件，即物理量（除压力）沿 方向的导数等于零。0,0,0,0,1,0,1,02,1(1)00()02(-,1)Myy
8、iyyi MinyxywallwallyyywallwallwalliiiiijjMu dyconstvujuyi 固壁边界。固体壁在涡量 流函数方法中，边界条件的处理：面是流线，在的固壁上，可以指定在的固壁上，可以指定在的固壁上，离散形式其中，虚拟网格点处的是未知的利用,0,1,02,1,1,00202iiiiiiyuy 0,0,00,1,0,1,0,10,0,10,221,1,0,1,0,0,10,0,10,22(2)()22022()2jyjinynnnnnnjjjjjjjjjjnnnnnjjjjjjuy dyxyvxxy 进口边界。流函数：可以用，通过适当的数值积分方法计算涡量：利用2
9、2,1,2,1,(3)0,02,xxxxxoutletoutletoutletMjMjMjMjMjvxxx出口边界。假定流动是充分发展的，有当采用一阶近似时1PoissonRe1Repup nnu 为了求解压力方程，还应补充压力的边界条件。在固壁处，利用法向分量作为边界条件，即n+1i,j1,1.-(1)0(2)1(3)(4)Poisson(5)(6)Poisson(7)ni jNavierStokesnnn流程采用涡量 流函数方法计算不可压缩方程的过程：初始化：给定初始速度场，初始化涡量和流函数，置时间指标计算开始，置求解涡量方程的差分格式，得到求解流函数的方程的差分格式，得到利用速度分量与
10、流函数之间的关系，计算速度场如果需要，求解压力方程计算压力场如果达到(2)规定的时间或者已经收敛到定常解，结束计算；否则转到。7.2.4 求解方法1111111,1,1,11,222.PoissonPoissonPoisson22()()nnnnnniji jiji ji ji jni jxy 方程的解法方程是椭圆型方程，代表物理量在空间的平衡分布过程，必须全场联立求解。在全场联立流函数的二维方程的差分格式，得到一个具有稀疏五对角系数矩阵的线性方程组。1111121,1,1,1,1,11,(0),(),(1)()(,1,1,Poisson1()4(1)Jacobi1(4nnnnnni jiji
11、ji ji ji jnnni ji ji jki jkki jijijxyhhk 在实际计算中，常采用迭代法求解。假定，二维方程的差分格式简化为的迭代初值第 次迭代的值记为几种迭代格式：迭代：)()()2,1,1,(1)(1)(1),1,1(1)(),1,)(2)Gauss-Seidel(,),0,0,JacobiJacobi1(4kkki ji ji jxykkki jiji jkki jijhi j iMjM迭代：迭代过程中，要对每一个网格点进行扫描，如果扫描的次序为则在迭代中，计算时，与均已经得到。因此，迭代可以改进为(1)()(1)21,1,1,)kkkiji ji ji jh充分利用
12、最新的迭代结果，可以加快收敛速度。(1)(1)(1)()()21,1,1,1,(1)(1)(1)()(,1,11,1,(3)111()444111(444kkkkkiji jiji ji ji jkkkkki ji ji jijijihjj线迭代：沿 方向的线迭代：在每个常数的网格线上求解一个三对角线方程组，增强了求解过程的耦合程度，可以加快收敛速度。沿 方向的线迭代：)2,(1)(1)(1),(1)(1)(),)(4)(1)011i jkkki ji ji jkkki ji ji jh 松弛迭代：松弛是迭代法中的常用技术。当用某种迭代方法计算出后，把视为中间结果，记作，本次迭代的最终值通过下
13、式计算松弛因子。当时，称为亚松弛，有利于迭代法的收敛，消除或推迟发散现象；称为超松弛，选择合适的松弛因子，有利于提高收敛速度，改进迭代法的计算效率。松弛因子一般存在一个上限，大于上限时，迭代过程不可能收敛。-Navier-Stokes-注：涡量 流函数方法在计算二维不可压缩方程中取得了成功。但是，当把涡量 流函数方法推广到求解三维问题时出现了困难。因为三维流动的流函数无法直接定义，需要引入多个流函数，其物理意义不如二维问题明确。另，根据涡量 流函数方法的思路，近年来，提出了计算三维不可压缩流动的涡量 速度方法，引入假想矢量势的涡量 矢量势方法等，但其应用还不够普遍。7.3 SIMPLE方法SI
14、MPLE(Semi-Implicit Method for Pressure-Linked Equations)SIMPLE方法：求解压力关联方程的半隐式方法讨论方法的有限体积形式。7.3.1 交错网格和非交错网格Navier-Stokes,(1/2,)(,1/2)SIMPLEu viji j非交错网格：所有物理量均存储于控制体的中心。直接在非交错网格上求解不可压缩方程，会出现数值解压力分布的异常。交错网格：压力、密度等标量因变量存储于控制体的中心，称为主控制体；速度分量分别存储在以和为中心的控制体中，称为动量控制体。不会出现压力分布的异常。下面介绍交错网格上的方法。交错网格7.3.2 动量方
15、程的离散22Navier-Stokes01Re1Reuvxyuuuvuputxyxvuvvvpvtxyy 守恒型的方程首先，假定压力是已知的，介绍动量方程的离散方案。1.控制方程及有限体积格式x考虑 方向的动量方程1/2,1,1/2,1/21/2,1/211/2,1/2,11111,1/2,1/21/2,1/2(1/2,)(7.3.4)()()0Euler()()0(7.3.ijiji jijijnnijijnnnniji jijijuijux yEEyFFxtuux yEEyFFxt 在动量控制体上积分，得时间方向采用隐式格式，8)20(7.3.4)1,Re1ReuEFtxyuEupxuFu
16、vy22111/2,1/2,111/2,1/2,1,1,1/2,1,1Re22xxxxxnnnxnnijiji jnnijijnx i ji jni jijnx i jE FECDpuCuDpxCuCu uuuuuuCuu uC 2.通量的离散其中，为对流项，为扩散项，压力：非线性项，采用隐式格式时，要进行线性化处理，令离散方案：中心差分 一阶迎风对流项,1,1/2,11,1/2,3/2,1,11,1/2,3/2,if0otherwise1.50.5if01.50.5otherwisei jni jijnni jijiji jnx i jnni jijijuu uuuuuCuuu 二阶迎风11
17、1,1/2,1/2,1111,1/2,1/2,1,max(0,),min(0,)2nnnx i ji jiji jiji ji ji ji jnnnnx i ji jiji jiji ji jni ji jCuuuuCuuu当采用一阶迎风格式离散时，其中，采用其它格式离散的对流项表示为一阶迎风格式加上修正项中心差分：,111/2,1/2,11111,1/2,3/2,1/2,3/2,111/2,1/2,1,20.5()()1Rei jnniji jijnnnnni ji jijiji jijijnnijijnjxx iuuuuuDuuDpuux二阶迎风：中心差分：取主控制扩散项压力项体中心的值1
18、1111,1/2,1/2,11ReRennnnni ji jiji jiji ji jEEuup通量 的通用离散形式为1,ReyyyyFFCDuCvuDy 对通量 的离散：其中，1/2,1/2,1/21,1/2111/2,1/2,11,1/2,1/21/2,1/211/2,1/21/2,1/2,1/21,1/2,1/211/2,1/21/2,11()22if0onniji jijnnijijny ijijnijijijny ijnijijyvvvuuCvvuvCCvu对进行线性化处理：令，中心差分：一阶迎风：111/2,1/21/2,1/2,11/2,1/21,1/2,1/2111/2,1/2
19、1/2,11/2,2therwise(1.50.5)if0(1.50.5)otherwisennijijijijny ijnnijijijvuuvCvuu二阶迎风：111/2,11/2,1,1/2,1/2()1Rennijijny ijuuDy 扩散项的离散：11111/2,1/21/2,1/21/2,1/2,1/21/2,11/2,1/21/2,1/21/2,1/21/2,1/21/2,1/211ReRemax(0,),min(0,)nnnnijijijijijijijijijijFFuuvv通量 的通用离散形式为1/2,1/2,1/2,1/2,3/2,3/2,1/2,11/2,11/2,1
20、1/2,11,1/2,1/2,1/2,1,1,3.()()2()Reijijijijijijijijijijniji jijijijijijuauauauauaux yy ppuatx yayt 动量方程的离散形式动量方程的离散形式1/2,1/21/2,1/21/2,3/2,1,2()Re11()()ReReijijiji jijijxayay ，1/2,11/2,1/21/2,11/2,1/21/2,1,1/2,1/21/2,1/2,1/2,1/21,1/21,1/21,1/21,1/2,1/211()()ReRe()()(ijijijijijiji jijiji ji jijijijiji
21、 jaxaxayxvbvbvbvb ，动量方程的离散形式,1/2,3/2,3/2,1,1/2,1/2)()i ji ji jni ji ji ji jvbvx yx ppvbt 7.3.3 SIMPLE方法11111/2,1/2,1/2,1/2SIMPLE-()()0nnnnijiji ji juuyvvx 方法：提供了解决压力 速度耦合问题的一种“半隐式”方案。确定压力场，使速度场满足连续方程的离散形式，基本思想：通过估计的压力场计算速度场的估计值；通过把离散的动量方程代入离散的连续性方程计算速度和压力的校正量。1/2,1/2,1/2,1/2,3/2,3/2,1/2,11/2,11/2,11
22、/2,11,1/2,1/2,1/2,1/21,1/21,1/21,1/21,1()ijijijijijijijijijijniji jijiji ji jijijijijauauauauaux yy ppuatbvbvbv/2,1/2,1/2,3/2,3/2,1,1/2,1/2()1i ji ji ji jni ji ji ji jbvbvx yx ppvbtnuv 通过上式计算时刻速度分量的估计值 和。2.1,(*)u vu vpnpppuuuvvv 压力修正和速度修正，一般不满足离散的连续方程式，因此，要对，和 进行修正。设时刻的数值解与估计值之间满足下列关系：校正量1.,nppppuv在
23、假定压力场下估算速度给出压力的估计值，一般可取代入和动量方程的离散形式，1/2,1/2,1,1/2,1/2,1,1/2,1/2,1,1/2,1/2,1/2,1,1/2SIMPLE()()()()ijijiji ji ji ji ji jijijiji jiji ji ji ji ji jauy ppbvx ppyuuppaxvvppb 方法在发展的早期主要用于求解定常问题，近似认为，速度修正方程为,1,1,1,1,1,1,1,1,(7.3.31)(7.3.31)11i ji jijijijiji ji ji ji ji jcpcpcpcpcpcpnuvpppnp代入离散的连续方程，得到压力修正
24、量满足的方程，用迭代法求解上式，得到压力修正量，代入得到时刻的 和，由计算时刻的。1/2,1/2,1/2,1/2,3/2,3/2,1/2,11/2,11/2,11/2,11,1/2,1/21,1/21,1/21,1/21,1/2,1/2,()ijijijijijijijijijijiji ji ji jijijijiji ji jauauauauauy ppbvbvbvbv校正量满足下列方程：1/2,3/2,3/2,1,()i ji ji ji jbvx pp3.边界条件在求解动量方程时，需要用到速度边界条件。不可压缩流动最常见的边界条件：给定速度或速度梯度。0,0uv假定边界是固体壁面，满足
25、,1/2(1/2)0i jvjv：动量控制体的中心在壁面动量方程的边界处理上方处，可以取法1/2,11/2,1/2,11/2,1/2,1/2(1/2,1)(1/2,)02ijijijijijijuuijuuuuu：引入虚拟网格点，如，取可以用于处方程的离散。在边界附近假定 沿竖直方向线性分布，动量方程的边界处理方法则满足物理边界条件。,1/2,1/2,1/2,1,1/2,1/2,1/2,1,1,1,1,1,1,1,1,1(1)0()0i ji ji ji ji ji ji ji ji ji ji ji jijijijiji ji ji ji jvxvvppbvvppcpcpcpcpcp压力修正
26、方程的边已知流动的法向速度，。根据速度修正方程代入压力界修正量满足的方程：条件：,1,1,1,1,1,1,1,1()(2)i ji ji ji jijijijiji ji ji ji jcccpcpcpcpcpu得消去了虚拟网格点上的未知量，得到一个封闭的方程组。给定边界处的压力。边界上的压力修正量为零，仿照动量方程的边界处理方法，引入虚拟网格点，在虚拟网格点上给定合适的压力修正量，使边界上压力修正量为零的条件得到满足。1.1,puvnuuuvvvppp由速度初值计算离散动量方程中各项的系数；2.给定压力场的估计值；3.数值求解离散的动量方程，计算速度分量的估计值 和；4.求解压力修正量满足的方程，计算压力修正值；5.求解速度修正量满足的方程，计算速度修正值；6.计算时刻的速度和压力，7.判断计算是否收敛；是，退出；否，算法流程：转到1。
