计算流体力学(中科院力学所)-第12讲-不可压流动

计算流体力学讲义第十二讲不可压缩Navier-Stokes方程的求解李新亮[e-mail]；力学所主楼219；82543801知识点：1讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”讲课录像及讲义上传至网盘拟压缩性方法求解压力Poisson方程法涡流函数法Simple方法CopyrightbyLiXinliang2知识回顾一、代数方程组的求解nnnnnnnnaaaaaaaaaaaaaaaaA321333323122322211131211bAx直接法nnnnnaaaaaaaaaaA0000000333223221131211Gauss消元法nnnnnnnnnnnnnnnnuuuuuuuuuullllllaaaaaaaaaaaaaaaa3332232211312113213231213213333231223222111312111111LU分解法追赶法：jjjjjjjdxcxbxa11jjjjBxAx1CopyrightbyLiXinliang3迭代法),(),(2222yxguyxfyuxu2,,1,1,,1,14jijijijijijifuuuuuJacobi迭代n+1nnnnGauss-Seidel迭代n+1nn+1n+13n+1n+1nn+1n+1n+1nnnLU-ADILU-SGS知识回顾CopyrightbyLiXinliang4知识回顾二、网格生成)(22xyy)(11xyyABCDEFA’B’C’D’E’F’物理空间计算空间1.代数网格生成法2.解椭圆型方程网格生成法CopyrightbyLiXinliang512.1不可压缩Navier-Stokes方程的特点VVVVV2Re10pt密度为常数的不可压缩Navier-Stokes方程组：特点：动量方程与能量方程解耦压力属于约束变量而不是发展变量温度对密度的影响可忽略不计压力不能时间推进求解可压缩不可压缩易处压力可推进求解，易于使用显格式压力方程具有椭圆性，无法推进求解。压力方程收敛性差难处可能出现间断不会出现间断研究重点激波捕捉压力处理)(VVVVCopyrightbyLiXinliang6概念澄清：压力——动力学压力及热力学压力动力学压力——应力的中各向同性部分ijijijpP211p22p12p22pnPpn连续介质微元体的受力平衡：应力的概念nppnpp静止流体或无粘流体中力的平衡——动力学压力的概念热力学压力——由分子动力学性质决定——状态方程热力学压力：分子对固壁的碰撞，产生压力RTp完全气体：动力学压力热力学压力可压缩N-S方程：动力学与热力学耦合；动力学压力=热力学压力不可压缩N-S方程：动力学与热力学解耦由不可压缩条件确定压力（纯动力学概念）1）压力的处理原则CopyrightbyLiXinliang7奇偶失联与交错网格uypyvvxvutvuxpyuvxuutuyvxu22Re1Re10压力项，通常采用中心差分离散yppypxppxpjijijijijiji2,21,1,,,1,1,极端情况：棋盘式压力场高压低压02,021,1,,1,1yppxppjijijiji特点：高压-低压点间隔分布采用中心差分格式计算出：0ypxp流场竟然“保持稳定”“奇偶失联”CopyrightbyLiXinliang8常用措施：交错网格压力p速度uuypyvvxvutvuxpyuvxuutuyvxu22Re1Re10速度v交错网格示意图jip,jiu,2/12/1,jivuxpyuvxuutu2Re1在u的网格点上离散jip,jiu,2/12/1,jivxppxpjijiji,,1,2/1uypyvvxvutv2Re1在v的网格点上离散xppypjijiji,1,2/1,jip,1注：对流项通常采用迎风格式离散xuuxuuxuu后差前差2uuuCopyrightbyLiXinliang92)对流项的处理原则VVVVV2Re10ptVVVVVVVV)()()()()(yvxuuyyvxuuyuvxuu关系式1：)()()(yvxuvyvvxvuyvvxvu关系式2：VωVV221VVωVV2Re1~0pVt兰姆-葛罗米柯等式2/~2Vpp总压表达式优点不足普通型简单，易于迎风有混淆误差守恒型简单，守恒有混淆误差旋度型混淆误差小计算量大螺旋-对称型混淆误差小计算量略大普通型守恒型VV)(VVVω2/)]([VVVVVωCopyrightbyLiXinliang1012.1人工压缩性方法（求解定常方程）VVVVV2Re10ptVVVVV2Re10pttp人工压缩性因子达到定常态0V0流动压缩时（），压力升高流动膨胀时（），压力降低0V0V0Vtp0V0V增大可令压力收敛加快，但会增加方程的刚性（降低时间步长）。人工压缩性因子相当于2ctpctppts202VctpCopyrightbyLiXinliang11对于定常问题，需要迭代到收敛VVVVV2Re10pttp)1,,max(111nnnnnnppvvuu对于非定常问题，需要内迭代（效率较低）VVVVV2Re10ptVVVVV2Re10pttpnknnkkkkkpttppVVVVVV211Re10Step1:得到n时间步的值nVStep2:进行如下内迭代直至收敛Step3:收敛后的V即为1nV内迭代收敛慢，效率较低；通常不使用人工压缩方法解非定常问题。CopyrightbyLiXinliang1212.2求解压力Poisson方法（投影法）1）压力的控制方程VVVVV2Re10pt对动量方程求散度)(2VVpVVVVVV22Re1)(ptpPoisson方程——压力的控制方程无法时间推进需联立求解，通常采用时间分裂法CopyrightbyLiXinliang132）投影法——求解微分型压力Poisson方程0)Re1(02ptVVVVV原理：将时间推进分成三个子步，中间步解出压力可时间推进不能时间推进Step1:预算步0)Re1(2*nntVVVVVStep2:压力修正步001*1nnptVVV*21Vtp求解，得到压力pStep3:最终步0*1ptnVV得到n+1时刻的V以1阶精度时间推进方法为例，实际上可采用更高阶精度时间推进方法：KarniadakisGE,IsraeliM,OrszagSA.1991High-ordersplittingmethodsfortheincompressibleNavier-stokesequations.J.Comp.Phys.97:414-443.CopyrightbyLiXinliang143）投影法——求解离散型压力Poisson方程压力修正步：将离散的动量方程带入离散的连续性方程，得到离散的压力方程001*1nnptVVV)3(0/)(/)()2(0)(/)1(0)(/12/1,12/1,1,2/11,2/11,1,*2/1,12/1,,,1*,2/11,2/1yvvxuuppytvvppxtuunjinjinjinjijijijinjijijijinji交错网格上离散将（1）（2）两式（离散的动量方程）带入（3）式（离散的连续性方程）可得到关于压力p的方程（离散的压力Poisson方程）；该方法可保证（3）式严格满足，因而相容性比方法2）更好CopyrightbyLiXinliang1512.3涡量-流函数方法（二维问题）)3(Re1)2(Re1)1(022uypyvvxvutvuxpyuvxuutuyvxu)3()2(xyxvyu)4(Re12yvxut引入流函数xvyu,)5(2（4）（5）式即涡量-流函数的控制方程计算结束后，如果需要计算压力，则求解如下方程22222yvxvyuxup)3()2(yxCopyrightbyLiXinliang16驱动方腔流动例：求解驱动方腔流动问题描述：如图示边长为L的方腔，上表面流体以常速度U运动，求解里面的流场（假设流动定常）。以涡量-流函数法为例：2Re1yvxut400ReULUL21）离散化jixxjiuuxu,,)(对流项：迎风差分2uuu迎风差分，建议采用高阶的粘性项：采用中心差分2,,1,1222xxjijiji也可采用更高阶的，可借助求差分系数的小程序时间推进：可采用显格式ttnn1CopyrightbyLiXinliang172采用中心差分离散：jijijijijijijiyx,2,1,1,2,,1,122可采用Jocabi，Gauss-Seidel等方法迭代提示：时间推进过程中的中间步无需迭代至收敛，最终（最后一个时间步）收敛即可。2）边界条件UL速度边界条件：上壁面u=1,v=0;其他壁面u=v=0;流函数的边界条件：0边界是一条流线，流线是流函数的等值线涡量的边界条件：由速度给出niniiyyuxuuxyu)(1)(12,1,2,00可用更高阶的格式CopyrightbyLiXinliang1812.4SIMPLE方法基本思想：与（离散型）投影法类似，但速度推进是隐式的；0)Re1(02ptVVVVV0)Re1(**2*ptnVVVVV1）已知预估压力计算速度*p采用隐式离散)1()(......)(*,*,1*,2/32*,2/11*,2/1jijijijijippuauau)2()(......)(*,*1,*2/3,2*2/1,1*2/1,jijijijijippvbvbv2）压力及速度修正'