计算流体力学课件：差分方法3

差分方法（3）知识点：激波捕捉格式——TVD、WENO、MUSCL、NND

1知识回顾1.差分格式的分辨率有效网格点数：一个波长里面的网格点数（PPW:PointperWavelength)修正波数2.群速度控制(GVC)33.流通矢量分裂“在二阶迎风与二阶中心格式中选一个”

“选接近一阶迎风的”GVC2格式4§5.1非物理振荡及TVD格式1.数值解中的非物理振荡

间断附近非物理振荡的根源理论1：色散误差导致各波传播速度不同（第4讲）理论2：粘性耗散不足思路：物理问题——有粘；物理粘性足以克服本身振荡

数值方法错误计算了物理粘性——不足以克服振荡物理问题本身也可能振荡。但如果错误计算物理粘性，则会错误地加剧（或衰减）振荡。1）非物理振荡的原因分析理论3：格式不能保单调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.005Re=2000Dx=0.0005相同6对流-扩散方程的特性：nn+1（线性）差分方程：某点的值是上一时刻周围几个点上值的线性组合物理上要求系数ak

均非负含义：某处浓度的增加对下一时刻周围浓度的影响为正。j-2j-1jj+1j+2差分方程单调性（无振荡）条件：差分方程（1）中的系数非负网格Reynolds数72）重要概念：网格Reynolds数以网格尺度度量的Reynolds数含义：

数值振荡——流动尺度为网格尺度

网格Reynolds数小，该尺度的能量被耗散掉——不发生振荡jj+1j-1过于苛刻的条件单方向网格点数106，三维1018

单纯靠物理粘性抑制振荡，网格间距必须足够小，通常难以实现网格足够小：不会发生振荡；网格小于激波的实际厚度，则不会振荡网格Reynolds数足够小时，物理粘性发挥作用，抑制振荡83）

人工粘性

物理粘性足够小才发挥作用，Reynolds数很高时很难做到

思路：人为增加粘性系数（添加人工粘性）抑制振荡优点：方法简便，有抑制振荡效果缺点：改变了物理问题，带来误差湍流、分离流等——对粘性敏感：非物理解分离流——对粘性敏感转捩——对粘性敏感很难计算对粘性敏感的问题改进措施：A:局部施加人工粘性B:高阶人工粘性VonNeumannMacCormack人工粘性系数大梯度区，加大人工粘性光滑区为二阶小量人工粘性项94)数值振荡的定量描述——总变差对于离散函数{uj}定义总变差:单调函数振荡函数j=1j=N含义：反映了振荡的剧烈程度双曲型守恒方程特点：沿特征线，u不变特征线未相交——总变差不变特征线相交——总变差减小结论：单个双曲型方程，总变差不增（TotalVariationDiminishing：TVD）102概念：单调格式、保单调格式与TVD格式n时刻：单调函数j=1j=Nn+1时刻：仍是单调函数j=1j=N设n时刻是单调的，如果n+1时刻的解仍保证单调，则称该格式为保单调格式。保单调格式基本结论：常系数的单调格式只能是一阶

单调格式必是保单调的；

线性格式，单调与保单调等价格式：如果满足则称其为单调格式。单调格式：单调格式保单调格式：TVD格式总变差不增TVD保单调单调113.TVD格式的理论基础——Harten定理Harten定理：如果差分格式可写成如下形式：且则格式（1）是TVD格式（1）可验证：Roe格式是TVD格式保证“系数非负”含义：“单调格式必是TVD格式”12例7.2.1：考虑线性单波方程：试讨论如下Lax-Wendroff格式二阶中心人工粘性是否满足Harten条件常系数的单调格式——只有一阶精度对比条件：不满足Harten条件13知识回顾：Lax-Wendroff格式Taylor展开，写出修正方程时-空二阶精度巧妙添加人工粘性，不但克服了不稳定性，而且抵消了时间误差，提高了时间精度类似方法：Beam-Warming格式人工粘性二阶精度迎风差分人工粘性，且提高时间精度特点：全离散、时刻耦合144.构建TVD格式思路：对现有格式进行改造，使之符合Harten条件通常在Roe、L-W、B-M（或其组合）基础上改进80年代初、这些格式是主流原格式（2阶）=1阶迎风+修正项

新格式=1阶迎风+限制函数*修正项1.以二阶中心及二阶迎风格式为基础的改造2阶迎风2阶中心2个候选格式：思路1：

两个里面选一个(GVC2,NND2)

思路2：

利用二者的组合改造、新格式二阶迎风二阶中心15显然格式为2阶中心可验证：格式为2阶迎风新格式：二者组合仍为二阶根据Harten定理，可知时，可满足TVD性质(2)精度条件二阶精度区TVD区二阶精度TVD区（二者交集）且16限制器（limiter)二阶中心二阶迎风GVC2格式NND2格式GVC2与NND2格式的限制器17常用的限制器Minmod型；Superbee型VanLeer型其他：如VanAlbada18概念：保单调区1阶迎风的预测值2阶中心的修正量2阶迎风的修正量精度差，但鲁棒性好精度高，但有些情况下预测结果“不靠谱”作为“标杆”检验高阶修正量是否可用趋势相反时，不可用；相差超过2倍时，不可用19历史上，TVD格式是在Roe、L-W、B-M（或其组合）基础上改进80年代初、这些格式是主流2以L-W格式为基础改造的格式L-W引入限制函数（限制器）1阶迎风部分修正项20显然格式为LW（2阶）可验证：格式为B-M（2阶）

新格式=1阶迎风+j*（LW格式-1阶迎风）新格式：LW,BM均为线性格式，二者组合仍为二阶根据Harten定理，可知时，可满足TVD性质(2)精度条件Beam-Warming二阶精度区TVD区二阶精度TVD区（二者交集）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-1jj+1

间断1间断2间断3……2.NND格式“在二阶迎风与二阶中心格式中选一个”

“选接近一阶迎风的”如果二阶中心与二阶迎风给出的修正趋势相反，干脆不修正。Minmod(a,b):a,b符号相同时，取绝对值小的；符号相反时，取0.k=1/3三阶迎风

k=0k=-1232.MUSCL格式（VanLeer）二阶迎风Formm格式k=1二阶中心1阶迎风修正部分间断处都会产生振荡思路：

振荡由修正项引起，需要对修正项进行限制24新格式函数minmod(a,b)

a,b符号相反，则返回0符号相同，则返回绝对值小的1阶迎风修正部分j-1jj+1j-2j+1/2思路：

修正部分是由j-1,j,j+1三点计算而得，是与的插值旧格式二者的组合：j-1,j计算；j,j+1计算

如果两个修正方案趋势（符号）相反，干脆不修正

如果趋势相同，则取（绝对值）最小者——避免振荡两个修正方案选前者：二阶迎风

选后者：二阶中心

组合之：Formm、3阶迎风…25含VanAlbada限制器的3阶MUSCL格式推荐方法：VanAlbada限制器光滑区相差不大间断区，

相差很大光滑区相差不大，s趋近于1间断区相差很大，s趋近于0光滑区，趋近3阶迎风间断区，趋近1阶迎风理论上并不满足TVD性质；但实际效果很好a<0的情况，把j+k改成j-k就可以了jj26关于MUSCL格式的一些注记1.MUSCL格式是一系列格式。可从二阶迎风、二阶中心、Formm,三阶迎风等格式经过限制器改造而得到。限制器有Minmod,VanAlbada等多种形式；MUSCL格式不一定都能满足TVD特性；NND格式是一种特殊的MUSCL格式（k=1,b=1的情况）MUSCL格式在CFD史上有重要贡献

MUSCL格式开创了利用基本变量重构（而不是利用通量重构）的新思路jj+1/2j-1j+1计算方法1：利用计算方法2：利用计算出

再利用计算出实际过程中，为了利用迎风技术，还要进行通量分裂27

§5.3WENO格式——高精度的激波捕捉法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；如果都光滑，则组合成更高阶格式。282.WENO格式的原理描述考虑线性单波方程：注：为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采用下一节介绍的守恒形式（1）确定网格基架点：6个点{j-3,j-2,j-1,j,j+1,j+2}

构造出该基架点上的目标差分格式计算这6个点可构造5阶迎风差分：该格式为WENO的“目标”格式，即，光滑区WENO逼近于该格式。利用Taylor展开，可唯一确定系数（可利用小程序coeff-schemes.f)实际上，还可利用分辨率优化技术，可构造出新的目标格式（降低精度、提高分辨率，见第4讲）。目前大量WENO的优化版做这种工作。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)也可运用优化技术，降低精度、提高分辨率……30（3）对这3个差分值进行加权平均，得到总的差分值原则：A.模板内函数越光滑，则权重越大；模板内有间断时，权重趋于0B.三个模拟内函数都光滑时，这三个三阶精度的逼近式可组合成一个五阶精度的逼近式。“理想权重”（3.1）确定理想权重令：5阶精度容易解出：31（3.2）度量每个模板内函数的光滑程度

IS越大，表示越不光滑。

光滑区，不同模板上的IS趋近同一值。具体形式见下一节。

（3.3）给出实际权重构造IS方法很多，例如：

：第k个模板光滑区逼近O（1）量级间断区量级，很大特点：间断区权重很小光滑区，趋近于理想权重（3.4）给出最终的差分逼近各阶导数的差分表达式都可作为光滑度量因子使用323.Jiang&Shu的五阶WENO格式守恒型；目前使用的WENO格式均为守恒型针对方程：模板1模板2

模板3构造差分格式如下：构造方法与前文相同（但注意这里构造的是通量，而前文是直接构造差分格式）针对整个网格基，构造出5阶精度的通量（理想情况下的通量）并构造出每个模板上的通量，计算出理想权重。仍利用程序coeff-schemes.f求系数理想权重光滑度量因子实际权重33光滑度量因子的计算（Jiang&Shu）k=1k=2k=3其中：j-2j-1jj+1j+2是使用模板k得到的插值函数

利用{j-2,j-1,j}点上的值构造的插值函数

，特点：光滑区趋近同一个值

非光滑区值远大于光滑区O(1)j点一阶、二阶导数的差分逼近（用模板k计算）代入34光滑因子ISk的进一步讨论（光滑区精度分析）同样：在光滑区，

趋近于同一个值，相互差距是一个6阶小量

是小量，可忽略5阶迎风格式的通量格式具有5阶精度35最终5阶WENO格式为正通量情况（a>0）

负通量情况（a<0）注：正通量差分格式中下标“j+k”改成“j-k”即得到负通量的差分格式(除了第1式不变）注意:是j-1/2而不是j+1/2k=1k=2k=3364.WENO格式的边界处理（1）简易的（降阶）处理方法：如果某模板用到边界外的点，简单将该模板权重设为0即可如果用到边界点外的点，则该权重设为0效果不错，但会边界降阶（推荐）（2）构造单边差分的WENO格式

优点：精度高

缺点：稳定性不易保证可能会出现负权重，造成不稳定，可用如下文献的方法处理：ShiJ,HuCQ,andShuCW,2002,ATechniqueofTreatingNegativeWeightsinWENOSchemes,JournalofComputationalPhysics175,108–127j=12345方法2：特征投影分裂（详细步骤见第4讲或第7讲4-5页）计算

利用上页的公式计算正通量及负通量的WENO通量及

计算

计算37推广到Euler（或N-S）方程运用分裂技术，可将上述方法推广到Euler方程方法1：逐点分裂（又称流通矢量分裂FVS）采用针对正通量（a>0）的方法计算采用针对负通量（a>0）的方法计算可采用Steger-Warming,L-F,VanLeer等分裂，见第4讲1)计算，它们是的函数（推荐使用Roe平均计算）效果更好但计算量较大计算量小但效果略差有轻微振荡385.WENO格式的改进MartinP的WENO-SYMBO格式（MartinP,JCP220,270-289,2006)思路：

原WENO格式采用迎风型网格基，耗散偏大

新的WENO格式采用对称型网格基，耗散小

采用优化方法，提高格式的分辨率理想权重情况下差分格式的修正波数；(SYMOO:未优化——8阶中心差分；SYMBO分辨率优化后的