第7讲-差分方法3

CopyrightbyLiXinliang1知识回顾1）

守恒型方程与守恒型格式守恒型方程：散度型守恒型格式：差量型习惯写为仅为记号，与j+1/2点上的值无关！守恒型方程+守恒型格式=解守恒“解守恒”：数值解的积分误差为0（如果边界准确）

保证总量（总质量、总动量、总能量）严格守恒（无误差）边界点早期的CFD:极为重视守恒性；目前CFD:仍很重视守恒性难点——复杂非线性系统的守恒性很难保证中间项全部消去，只剩两端CopyrightbyLiXinliang2）

通量分裂——便于使用迎风格式方法（1）：逐点分裂（Steger-Warming,VanLeer,L-F）原理：利用了性质使得的Jocabian阵特征值纯正或纯负优点：无需矩阵运算，计算量小，使用方便不足：仅重新组合，没有做到真正解耦。原因：具体方法：Steger-Warming

L-FVanLeerorCopyrightbyLiXinliang3方法2）特征（投影）分裂在网格基上冻结系数…j-2j-1jj+1…在基架点上系数不变优点：特征分解，（局部）解耦——耗散小，数值振荡低

缺点：大量矩阵运算，计算量大每计算一个点的导数，要进行m次矩阵运算（m为网格基点数）原理描述（非守恒，很少采用；实际上使用下一页的方法）

守恒型方式计算…j-2j-1jj+1…在基架点上系数不变具体步骤：

假设已知U，且针对模型方程（线性单波方程）

已构造出差分格式（1）1）计算出教材130页的公式(6.1.11-6.1.13),式中用到各变量在j+1/2的值（例如）

可使用j,j+1点值的算术平均（如）或Roe平均；由计算；方法很多，例如前面介绍的或4CopyrightbyLiXinliang

均可推荐使用Roe-平均！2）

在网格基上计算…j-2j-1jj+1…计算fj+1/2用到的点注意，在该网格基上（例如k=j-1,j,j+1）保持不变例如：3)利用已构造好的差分格式，计算通量4）得到总通量

5）计算差分（j点处）步骤的算法描述（注意：实际上是两重循环）doj=1,Ndok=j-1,j+1(网格基，可以是更多或更少点)

enddoenddodoj=1,N

enddo需要多次矩阵运算，计算量大

守恒性好，耗散小，数值解质量好5CopyrightbyLiXinliang

通量分裂+迎风差分引入数值耗散

分裂本身不带来耗散，但会放大（或减少）差分的耗散举例：分裂过程耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。=+向上平移向下平移分裂后的流场越偏离原先流场，则总体耗散越大如使用低精差分度格式，则对分裂形式敏感（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感（可使用逐点分裂）6CopyrightbyLiXinliang概念澄清耗散放大系数CopyrightbyLiXinliang7§7.1Roe格式——守恒型格式的范例破坏守恒性后果很严重（?）为了便于使用迎风格式、特征分裂解耦，通常把守恒型方程改写为非守恒型守恒方程+守恒格式=解守恒方程不守恒，即使差分方法守恒，也无法做到解守恒由于a非常数，无法消去中间项！不再守恒思路：保证特征方向，找回守恒性守恒型方程优点：守恒性非守恒方程优点：清晰的特征方向兼顾守恒与非守恒方程CopyrightbyLiXinliang81.单方程的Roe格式1阶迎风（直接从守恒方程出发）（1）（2）体现了特征方向只有这种表达式，才能保证（2）与（1）等价（3）or都无法保证（2）与（1）等价。简单的线性平均不行(非线性系统，中点的斜率不等于平均斜率）关键：构造Roe格式“平均斜率”，不等于“斜率的平均值”，也不等于中点处的斜率平均斜率CopyrightbyLiXinliang92.方程组的情况（Roe格式的意义）

需满足如下条件（Uniform特性）单方程的简单推广

1）连续，且2）可通过相似变换对角化，即保证双曲性3）对于任意UR,UL有单方程的推广，含义为平均增长率标量方程向矩阵方程的简单推广，但构造很困难。均不满足Uniform特性平均斜率CopyrightbyLiXinliang103.矩阵的构造关键：“向量除以向量”？直接求平均增长率：uf(u)uLuRuRoeRoe点的斜率为平均斜率（根据拉格朗日中值定理，[UL,UR]区间内肯定存在Roe点）思路1：在UL与UR之间寻找一个点URoe,该点处的增长率为平均增长率f(u)=u2u二次函数——Roe点与中点重合标量函数的启示：Roe点肯定存在（Langrage中值定理）

二次函数的中点即为Roe点思路2：进行坐标变换，得到一个二次（齐）函数引入如果

是二次（齐）函数，则其中点即为Roe点重要启示更准确地讲，应当是要求为W的线性函数，

即增长率为线性函数（中点处的增长率刚好为平均增长率）CopyrightbyLiXinliang11针对Euler方程的具体构造引入新变量：则：目的：使得F(w)是W二次齐函数（增长率为线性函数）f(U)不是U的二次齐函数二次齐函数！中点处的斜率即为平均斜率。Roe点Roe点为：增长率为线性函数！CopyrightbyLiXinliang12最终：其中如下计算：平均增长率（矩阵）含义：左、右两个状态点的某种平均（称为Roe平均，为密度加权平均）

该状态点对应的增长率（矩阵）为平均增长率（矩阵）

实际上是一种“等效平均”。效果优于简单的算数（或几何）平均。

三维情况下，还有其他量（如压力、温度、音速等）用这三个量计算（5）简单易记：CopyrightbyLiXinliang13Roe格式的计算步骤（半离散）已知n时刻所有网格点上的物理量，对于j点：1）令UR=Uj+1,UL=Uj

（密度、压力、速度等）

2）采用Roe平均公式（5）计算Roe平均值3）将Jacobian矩阵进行特征分解:

计算4)计算5）计算6）计算空间导数7）时间推进，计算下一时间步的值。j-1jj+1与前文（第3，4讲）的形式相同，仅需把式中的密度、压力、速度等换成经过Roe平均的密度、压力、速度即可其中：CopyrightbyLiXinliang14可能出现导数不连续，可能引起数值振荡实际使用时可用如下函数代替——所谓“熵修正”实际上是在特征值0点周围增加了耗散Roe格式的优点：1）保持守恒性的同时，严格保证了特征方向2）便于推广到高精度格式——特征投影分裂中使用Roe平均即可（见本PPT第5页）。推广到高阶后，虽不再保证严格的特征方向，但仍优于采用算数平均方法。Roe格式的不足：

本身精度只有一阶；

推广到高阶后，特征方向无法严格保证；

推广到二维或三维后，特征方向无法严格保证，出现振荡。CopyrightbyLiXinliang15作业7.1对于一维Euler方程：引入新变量：推导出及其Jacobian矩阵的具体表达式（以W为自变量），并证明对于任意，有：提示：写出表达式后，将向量分别代入上式左、右两端，容易证明相等。要求：推导过程要详细，切勿简单从书本上摘抄。重要的CFD基本功练习，一定要重视！针对如下Sod激波管问题

用Roe格式计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。

要求：1）空间网格数100，时间推进格式选用3阶Runge-Kutta,时间步长自选。2）尝试使用熵修正与不使用熵修正两种情况（见本PPT15页）3）欢迎与其他数值方法得到的结果对比（最好画在同一张图上，便于比较）。16CopyrightbyLiXinliang作业7.2CopyrightbyLiXinliang17§7.2非物理振荡及TVD格式1.数值解中的非物理振荡

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

数值方法错误计算了物理粘性——不足以克服振荡物理问题本身也可能振荡。但如果错误计算物理粘性，则会错误地加剧（或衰减）振荡。1）非物理振荡的原因分析CopyrightbyLiXinliang18数值实验

二阶中心差分计算域[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相同CopyrightbyLiXinliang19对流-扩散方程的特性：nn+1差分方程：某点的值是上一时刻周围几个点上值的线性组合物理上要求系数ak

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

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

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

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

人工粘性

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

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

单调格式必是保单调的；

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

新格式=1阶迎风+限制函数*修正项引入限制函数（限制器）1阶迎风部分修正项CopyrightbyLiXinliang28显然格式为LW（2阶）可验证：格式为B-M（2阶）

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

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

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

IS越大，表示越不光滑。

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

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

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

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

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

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

非光滑区值远大于光滑区O(1)j点一阶、二阶导数的差分逼近（用模板k计算）代入CopyrightbyLiXinliang36最终5阶WENO格式为正通量情况（a>0）

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

优点：精度高

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

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

计算

计算CopyrightbyLiXinliang38推广到Euler（或N-S）方程运用分裂技术，可将上述方法推广到Euler方程方法1：逐点分裂（又称流通矢量分裂FVS）采用针对正通量（a>0）的方法计算采用针对负通量（a>0）的方法计算可采用Steger-Warming,L-F,VanLeer等分裂，见第4讲1)计算，它们是的函数（推荐使用Roe平均计算）效果更好但计算量较大计算量小但效果略差有轻微振荡数值测试：1维Sod激波管问题，网格点128

差分方法:7阶精度经典WENO(Jiang&Shu)分裂方式：Steger-WarmingFVS;特征投影分裂（算术平均/Roe平均）t=0.1时刻密度及速度分布Steger-WarmingFVS+WENO仍无法避免振荡特征投影分裂+WENO可避免振荡