中科院计算流体力学最新讲义CFD115讲差分方法.ppt

计算流体力学讲义2011 第五讲 差分方法（3）李新亮；力学所主楼219；82543801,知识点：激波捕捉格式 TVD、WENO、MUSCL、NND,1,Copyright by Li Xinliang,讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”下载地址2：,知识回顾,1.差分格式的分辨率,有效网格点数：一个波长里面的网格点数（PPW:Point per Wavelength),修正波数,2.群速度控制(GVC),Copyright by Li Xinliang,3,3.流通矢量分裂,“在二阶迎风与二阶中心格式中选一个”“选接近一阶迎风的”,GVC2 格式,Copyright by Li Xinliang,4,5.1 非物理振荡及TVD格式,1.数值解中的非物理振荡,间断附近非物理振荡的根源,理论1：色散误差导致各波传播速度不同（第4讲）,理论2：粘性耗散不足,思路：物理问题 有粘；物理粘性足以克服本身振荡 数值方法错误计算了物理粘性不足以克服振荡,物理问题本身也可能振荡。但如果错误计算物理粘性，则会错误地加剧（或衰减）振荡。,1）非物理振荡的原因分析,理论3：格式不能保单调,Copyright by Li Xinliang,5,数值实验,二阶中心差分,计算域0,1,网格点201（Dx=0.005）时间步长Dx=0.0005,T=0.1时刻的u分布,Re=200Dx=0.005,现象：Dx一定时，减小Reynolds数可抑制振荡 Reynolds数一定时，减小Dx可抑制振荡,暗示,是某一特征量,Re=2000Dx=0.005,Re=2000Dx=0.0005,相同,Copyright by Li Xinliang,6,对流-扩散方程的特性：,n,n+1,差分方程：,某点的值是上一时刻周围几个点上值的线性组合,物理上要求系数 ak 均非负,含义：某处浓度的增加对下一时刻周围浓度的影响为正。,j-2 j-1 j j+1 j+2,差分方程单调性（无振荡）条件：差分方程（1）中的系数非负,网格Reynolds数,Copyright by Li Xinliang,7,2）重要概念：网格 Reynolds数 以网格尺度度量的Reynolds数,含义：数值振荡 流动尺度为网格尺度 网格 Reynolds数小，该尺度的能量被耗散掉 不发生振荡,j,j+1,j-1,过于苛刻的条件,单方向网格点数106，三维1018,单纯靠物理粘性抑制振荡，网格间距必须足够小，通常难以实现,网格足够小：不会发生振荡；网格小于激波的实际厚度，则不会振荡,网格Reynolds数足够小时，物理粘性发挥作用，抑制振荡,Copyright by Li Xinliang,8,3）人工粘性,物理粘性 足够小才发挥作用，Reynolds数很高时很难做到,思路：人为增加粘性系数（添加人工粘性）抑制振荡,优点：方法简便，有抑制振荡效果缺点：改变了物理问题，带来误差,湍流、分离流等对粘性敏感：非物理解,分离流 对粘性敏感,转捩对粘性敏感,很难计算对粘性敏感的问题,改进措施：A:局部施加人工粘性 B:高阶人工粘性,Von Neumann,MacCormack,Copyright by Li Xinliang,9,4)数值振荡的定量描述 总变差,对于离散函数uj 定义总变差:,单调函数,振荡函数,j=1,j=N,含义：反映了振荡的剧烈程度,双曲型守恒方程,特点：沿特征线，u不变,特征线未相交总变差不变,特征线相交 总变差减小,结论：单个双曲型方程，总变差不增（Total Variation Diminishing：TVD）,Copyright by Li Xinliang,10,2 概念：单调格式、保单调格式与TVD格式,n时刻：单调函数,j=1,j=N,n+1时刻：仍是单调函数,j=1,j=N,设n时刻 是单调的，如果n+1时刻的解 仍保证单调，则称该格式为保单调格式。,保单调格式,基本结论：常系数的单调格式只能是一阶 单调格式必是保单调的；线性格式，单调与保单调等价,格式：如果满足 则称其为单调格式。,单调格式：,单调格式,保单调格式：,TVD格式,总变差不增,TVD,保单调,单调,Copyright by Li Xinliang,11,3.TVD格式的理论基础 Harten定理,Harten定理：,如果差分格式可写成如下形式：,且,则格式（1）是TVD格式,（1）,可验证：Roe格式是TVD格式,保证“系数非负”,含义：“单调格式必是TVD格式”,Copyright by Li Xinliang,12,例：,考虑线性单波方程：,试讨论如下 Lax-Wendroff 格式,二阶中心,人工粘性,是否满足Harten条件,常系数的单调格式 只有一阶精度,对比条件：,不满足Harten 条件,Copyright by Li Xinliang,13,知识回顾：Lax-Wendroff 格式,Taylor展开，写出修正方程,时-空二阶精度,巧妙添加人工粘性，不但克服了不稳定性，而且抵消了时间误差，提高了时间精度,类似方法：Beam-Warming 格式,人工粘性,二阶精度迎风差分,人工粘性，且提高时间精度,特点：全离散、时刻耦合,Copyright by Li Xinliang,14,4.构建TVD 格式,思路：对现有格式进行改造，使之符合Harten条件,通常在Roe、L-W、B-M（或其组合）基础上改进80年代初、这些格式是主流,原格式（2阶）=1阶迎风+修正项,新格式=1阶迎风+限制函数*修正项,1.以二阶中心及二阶迎风格式为基础的改造,2阶迎风,2阶中心,2个候选格式：,思路1：两个里面选一个(GVC2,NND2)思路2：利用二者的组合,改造、新格式,二阶迎风,二阶中心,Copyright by Li Xinliang,15,显然 格式为2 阶中心可验证：格式为2阶迎风,新格式：,二者组合仍为二阶,根据Harten定理，可知,时，可满足TVD性质,(2)精度条件,二阶精度区,TVD区,二阶精度TVD区（二者交集）,Copyright by Li Xinliang,16,限制器（limiter),二阶中心,二阶迎风,GVC2 格式,NND2 格式,GVC2与NND2格式的限制器,Copyright by Li Xinliang,17,常用的限制器,Minmod型；,Superbee型,Van Leer型,其他：如Van Albada,Copyright by Li Xinliang,18,概念：保单调区,1阶迎风的预测值,2阶中心的修正量,2阶迎风的修正量,精度差，但鲁棒性好,精度高，但有些情况下预测结果“不靠谱”,作为“标杆”检验高阶修正量是否可用,趋势相反时，不可用；相差超过2倍时，不可用,Copyright by Li Xinliang,19,历史上，TVD格式是在Roe、L-W、B-M（或其组合）基础上改进80年代初、这些格式是主流,2 以 L-W格式为基础改造的格式,L-W,引入限制函数（限制器）,1阶迎风部分,修正项,Copyright by Li Xinliang,20,显然 格式为LW（2 阶）可验证：格式为B-M（2阶）,新格式=1阶迎风+j*（LW格式-1阶迎风）,新格式：,LW,BM均为线性格式，二者组合仍为二阶,根据Harten定理，可知,时，可满足TVD性质,(2)精度条件,Beam-Warming,二阶精度区,TVD区,二阶精度TVD区（二者交集）,Copyright by Li Xinliang,21,5.2 其他激波捕捉格式简介,1.Godnov格式,思路：利用精确Riemann解难点：精确Riemann解要求间断两侧物理量为常数 对策：采用分片常数代替原先的函数,方法：1）采用分片常数代替原先的函数 即假设xj-1/2,xj+1/2 区间物理量为xj点上的值 2）在每个间断处，精确求解Riemann问题，得到Dt时刻物理量的分布。（假设Dt足够小，各间断之间无相互影响，因此可独立求解）。3）在区域xj-1/2,xj+1/2 上平均，得到Dt时刻xj点上物理量的值。4）重复1-3，直到指定时刻。,j-1 j j+1,间断1 间断2 间断3,2.NND 格式,“在二阶迎风与二阶中心格式中选一个”“选接近一阶迎风的”如果二阶中心与二阶迎风给出的修正趋势相反，干脆不修正。,Minmod(a,b):a,b 符号相同时，取绝对值小的；符号相反时，取0.,k=1/3 三阶迎风,k=0,k=-1,Copyright by Li Xinliang,23,2.MUSCL格式（Van Leer）,二阶迎风,Formm格式,k=1 二阶中心,1阶迎风,修正部分,间断处都会产生振荡,思路：振荡由修正项引起，需要对修正项进行限制,Copyright by Li Xinliang,24,新格式,函数 minmod(a,b)a,b 符号相反，则返回0 符号相同，则返回绝对值小的,1阶迎风,修正部分,j-1,j,j+1,j-2,j+1/2,思路：修正部分是由 j-1,j,j+1 三点计算而得，是 与 的插值,旧格式,二者的组合：j-1,j 计算；j,j+1计算,如果两个修正方案趋势（符号）相反，干脆不修正 如果趋势相同，则取（绝对值）最小者 避免振荡,两个修正方案,选前者：二阶迎风 选后者：二阶中心 组合之：Formm、3阶迎风,Copyright by Li Xinliang,25,含Van Albada限制器的3阶MUSCL格式,推荐方法：,Van Albada限制器,光滑区 相差不大,间断区，相差很大,光滑区 相差不大，s趋近于1,间断区 相差很大，s趋近于0,光滑区，趋近3阶迎风,间断区，趋近1阶迎风,理论上并不满足TVD性质；但实际效果很好,a0 的情况，把j+k改成j-k就可以了,j,j,Copyright by Li Xinliang,26,关于MUSCL格式的一些注记,1.MUSCL格式是一系列格式。可从二阶迎风、二阶中心、Formm,三阶迎风等格式经过限制器改造而得到。限制器有Minmod,Van Albada等多种形式；MUSCL格式不一定都能满足TVD特性；NND格式是一种特殊的MUSCL格式（k=1,b=1的情况）MUSCL格式在CFD史上有重要贡献,MUSCL格式开创了利用基本变量重构（而不是利用通量重构）的新思路,j,j+1/2,j-1,j+1,计算,方法1：利用 计算,方法2：利用 计算出 再利用 计算出,实际过程中，为了利用迎风技术，还要进行通量分裂,Copyright by Li Xinliang,27,5.3 WENO格式 高精度的激波捕捉法,1.基本思路,j-3,j-2,j-1,j,j+1,j+2,j-3,j-2,j-1,j；j-2,j-1,j,j+1；j-1,j,j+1,j+2,五个基架点被分成三个组,1）若高精度逼近，必然利用多个基架点 2）如果该基架点内函数有间断，会导致振荡 3）间断不可能处处存在 4）把基架点分成多个组（模板），每个模板独立计算j点导数的逼近。得到多个差分 5）根据每个模板的光滑程度，设定权重 6）对多个差分结果进行加权平均。光滑度越高，权重越大。如果某模板存在间断，则权重趋于0；如果都光滑，则组合成更高阶格式。,Copyright by Li Xinliang,28,2.WENO格式的原理描述,考虑线性单波方程：,注：为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采用下一节介绍的守恒形式,（1）确定网格基架点：6个点 j-3,j-2,j-1,j,j+1,j+2 构造出该基架点上的目标差分格式,计算,这6个点可构造5阶迎风差分：,该格式为WENO 的“目标”格式，即，光滑区WENO 逼近于该格式。,利用Taylor展开，可唯一确定系数（可利用小程序coeff-schemes.f),实际上，还可利用分辨率优化技术，可构造出新的目标格式（降低精度、提高分辨率，见第4讲）。目前大量WENO的优化版做这种工作。,Copyright by Li Xinliang,29,将这6个基架点分割成3个组（称为模板）每个组独立计算 的差分逼近,模板1,模板2,模板3,模板1：j-3,j-2,j-1,j 模板2：j-2,j-1,j,j+1模板3：j-1,j,j+1,j+2,利用这三个模板的基架点，可构造出逼近 的3阶精度差分格式,计算j点的导数u,竟然算出了三个不同的值，怎么办？ENO 方法：选择最优（最光滑）的，舍弃其余两个WENO的处理方法：三个都要，加权平均它们。,利用Taylor展开式，可唯一确定这些系数）（可利用小程序coeff-schemes.f),也可运用优化技术，降低精度、提高分辨率,Copyright by Li Xinliang,30,（3）对这3个差分值进行加权平均，得到总的差分值,原则：A.模板内函数越光滑，则权重越大；模板内有间断时，权重趋于0 B.三个模拟内函数都光滑时，这三个三阶精度的逼近式可组合成一个五阶精度的逼近式。,“理想权重”,（3.1）确定理想权重,令：,5阶精度,容易解出：,Copyright by Li Xinliang,31,（3.2）度量每个模板内函数的光滑程度,IS越大，表示越不光滑。光滑区，不同模板上的IS趋近同一值。具体形式见下一节。,（3.3）给出实际权重,构造IS方法很多，例如：,：第k个模板,光滑区逼近 O（1）量级间断区 量级，很大,特点：间断区权重很小 光滑区，趋近于理想权重,（3.4）给出最终的差分逼近,各阶导数的差分表达式都可作为光滑度量因子使用,Copyright by Li Xinliang,32,3.Jiang&Shu 的五阶WENO格式,守恒型；目 前使用的WENO格式均为守恒型,针对方程：,模板1,模板2,模板3,构造差分格式如下：,构造方法与前文相同（但注意这里构造的是通量，而前文是直接构造差分格式）针对整个网格基，构造出5阶精度的通量（理想情况下的通量）并构造出每个模板上的通量，计算出理想权重。,仍利用程序coeff-schemes.f求系数,理想权重,光滑度量因子,实际权重,Copyright by Li Xinliang,33,光滑度量因子的计算（Jiang&Shu）,其中：,j-2 j-1 j j+1 j+2,是使用模板k 得到的插值函数,利用j-2,j-1,j点上的值构造的插值函数,，,特点：光滑区趋近同一个值 非光滑区值远大于光滑区 O(1),j 点一阶、二阶导数的差分逼近（用模板k计算）,代入,Copyright by Li Xinliang,34,光滑因子 ISk 的进一步讨论（光滑区精度分析）,同样：,在光滑区，趋近于同一个值，相互差距是一个6阶小量,是小量，可忽略,5阶迎风格式的通量,格式具有5阶精度,Copyright by Li Xinliang,35,最终5阶 WENO 格式为,正通量情况（a0）,负通量情况（a0）,注：正通量差分格式中下标“j+k”改成“j-k”即得到负通量的差分格式(除了第1式 不变）,注意:是j-1/2而不是j+1/2,Copyright by Li Xinliang,36,4.WENO格式的边界处理,（1）简易的（降阶）处理方法：,如果某模板用到边界外的点，简单将该模板权重设为0即可,如果用到边界点外的点，则该权重设为0,效果不错，但会边界降阶（推荐）,（2）构造单边差分的WENO格式 优点：精度高 缺点：稳定性不易保证可能会出现负权重，造成不稳定，可用如下文献的方法处理：Shi J,Hu CQ,and Shu CW,2002,A Technique of Treating Negative Weights in WENO Schemes,Journal of Computational Physics 175,108127,j=1 2 3 4 5,方法2：特征投影分裂（详细步骤见第4讲或第7讲4-5页）,计算 利用上页的公式计算正通量 及负通量 的WENO通量 及 计算 计算,Copyright by Li Xinliang,37,推广到Euler（或N-S）方程,运用分裂技术，可将上述方法推广到Euler方程,方法1：逐点分裂（又称流通矢量分裂 FVS）,采用针对正通量（a0）的方法计算,采用针对负通量（a0）的方法计算,可采用Steger-Warming,L-F,Van Leer 等分裂，见第4讲,1)计算,，它们是 的函数（推荐使用Roe平均计算）,效果更好但计算量较大,计算量小但效果略差有轻微振荡,Copyright by Li Xinliang,38,5.WENO格式的改进,Martin P 的WENO-SYMBO格式（Martin P,JCP 220,270-289,2006),思路：原WENO格式采用迎风型网格基，耗散偏大 新的WENO 格式采用对称型网格基，耗散小 采用优化方法，提高格式的分辨率,理想权重情况下差分格式的修正波数；(SYMOO:未优化 8阶中心差分；SYMBO 分辨率优化后的4阶格式）,针对如下Sod 激波管问题,用5阶WENO格式计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。要求：空间网格数100，时间推进格式选用3阶Runge-Kutta,时间步长自选。,39,Copyright by Li Xinliang,作业 5.1,针对如下Shu-Osher 激波-密度扰动波干扰问题：,用5阶WENO格式计算其数值解，画出t=0.1时刻密度、速度及压力的分布 要求：1）空间网格数200，时间推进格式选用3阶Runge-Kutta,时间步长自选。2）使用2000个网格点计算，其结果作为“精确解”，与其它结果画在一起，便于比较。,40,Copyright by Li Xinliang,作业 5.2,初值为:,41,Copyright by Li Xinliang,作业 5.3（选作题）,推导7阶精度的WENO格式，给出详细的推导过程及格式的具体表达式,提示：与5阶WENO推导思路相同，但网格基架点扩大了，模板数目也增加了,j-3 j-2 j-1 j j+1 j+2 j+3,正通量(a0)的WENO通量 使用基架点j-3,j-2,j-1,j,j+1,j,j+2,j+3如上图示，将其分割为4个组（模板），每个模板上4个基架点。1）先构建整个网格基（7个点）上7阶迎风格式的通量表达式 2）对于每个模板，构造逼近j点导数的4阶差分格式的通量表达式 3）按照理想权重进行组合可得到理想格式（7阶迎风格式）对照具体表达式，求出理想权重Ck 4）构建WENO通量的加权表达式 5）仿照 本PPT第35页的方法，给出光滑度量因子的具体表达式,可利用小程序求系数coeff-schemes.f减小工作量,