第7讲-差分方法3

1、Copyright by Li Xinliang1知识回顾知识回顾1） 守恒型方程与守恒型格式守恒型方程与守恒型格式守恒型方程：守恒型方程： 散度型散度型0)(xuftu0)(UFUt守恒型格式：守恒型格式： 差量型差量型1jjjxf习惯写为习惯写为)(12/12/1jjjffxxf仅为记号，仅为记号，与与j+1/2点点上的值无关上的值无关！守恒型方程守恒型方程+ 守恒型格式守恒型格式=解守恒解守恒“解守恒解守恒”： 数值解的积分误差为数值解的积分误差为0 （如果边界准确）（如果边界准确） 保证总量（总质量、总动量、总能量）严格守恒保证总量（总质量、总动量、总能量）严格守恒 （无误差）（无误差

2、）0)(12/12/1jjjffxtu0)(10)(12/12/12/12/1Njjjjjjjffxutffxtu边界点早期的早期的CFD:极为重视守恒性极为重视守恒性；目前目前CFD: 仍很重视守恒性仍很重视守恒性难点复杂非线性系统的守恒性很难保证.)()()()(2/52/72/32/52/12/32/12/1ffffffffjjj中间项全部消去，只剩两端Copyright by Li Xinliang2） 通量分裂通量分裂 便于使用迎风格式便于使用迎风格式0 xtf(U)U方法（方法（1）：）： 逐点分裂逐点分裂 （Steger-Warming, Van Leer, L-F）)(0UfA

3、AUxUtSUSUAffff1,原理： 利用了性质xxUAAU)(使得使得UAf的的Jocabian阵特征值阵特征值纯正或纯负纯正或纯负优点：优点： 无需矩阵运算，计算量小，使用方便无需矩阵运算，计算量小，使用方便不足：不足： 仅重新组合，没有做到真正解耦。仅重新组合，没有做到真正解耦。 原因：原因：xxxUAUAf具体方法：具体方法： Steger-Warming L-F Van Leer2)(2/122kkk2*kk2kkkorCopyright by Li Xinliang3方法方法2） 特征（投影）分裂特征（投影）分裂)()(jjjjjjjxxxxxVVSUSSUSSUA1j1jj1j

4、00 xtxtUAUf(U)U0jjjxtUAU在网格基上冻结系数 j-2 j-1 j j+1 在基架点上系数在基架点上系数 不变不变jjxUAjAUSVj)(12/12/1jjjxxVVV 优点：特征分解，（局部）解耦优点：特征分解，（局部）解耦耗散小，数值振荡低耗散小，数值振荡低 缺点：大量矩阵运算，计算量大缺点：大量矩阵运算，计算量大 每计算一个点的导数，要进行每计算一个点的导数，要进行m次矩次矩阵运算阵运算 （m为网格基点数）为网格基点数）原理描述原理描述 （非守恒，很少采用；（非守恒，很少采用； 实际上使用下一页的方法）实际上使用下一页的方法） 守恒型方式守恒型方式计算计算x)(Uf

5、xxjjj/ )()(2/12/1ffUf j-2 j-1 j j+1 在基架点上系数 不变jjxUAjA)(2/12/112/12/1jjjjSVVfUSV2/1j具体步骤：具体步骤： 假设已知假设已知 U， 且针对模型方程（线性单波方程）且针对模型方程（线性单波方程） 已构造出差分格式已构造出差分格式xvvvxvvvjjjjjj/ )(;/ )(2/12/12/12/10 xvatv（1）1） 计算出计算出2/12/112/1,jjjSS教材教材130页的公式页的公式(6.1.11-6.1.13), 式中用到各变量在式中用到各变量在j+1/2的值（例如的值（例如 ） 可使用可使用j, j+

6、1 点值的算术平均点值的算术平均 （如（如 ） 或或Roe平均平均 ； 由由 计算；方法很多，例如前面介绍的计算；方法很多，例如前面介绍的 或或 2/1ju2/ )(12/1jjjuuu4Copyright by Li Xinliang2kkk2*kk 均可2/1j2/1j2kkk2*kk推荐使用推荐使用Roe-平均！平均！2） 在网格基上计算在网格基上计算) 1, 1.(2/12/1jjjkegkjjkUSV j-2 j-1 j j+1 计算fj+1/2用到的点注意，在该网格基上（例如注意，在该网格基上（例如k=j-1,j,j+1） 保持不变保持不变2/12/1,jjSjjjjjjjjjjj

7、jjuuuuwhenuuuuuuwhenuuu1111112/12/ )(2/ )3(例如： 1211121212/12/ )(2/ )3(jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuu3) 利用已构造好的差分格式，计利用已构造好的差分格式，计算通量算通量4） 得到总通量得到总通量2/12/1,jjVV)(2/12/112/12/1jjjjVVSf 5） 计算差分计算差分 （j点处）点处）xxjjj/ )()(2/12/1ffUf步骤的算法描述步骤的算法描述 （注意：（注意： 实际上是两重循环）实际上是两重循环） do j=1,N do k=j-1,j+1 (网格基，可以

8、是更多或更少点) enddo enddo do j=1,N enddo),(),(112/1112/1jjjjjjjjffVV,VVVV,VV)(2/12/112/12/1jjjjSVVfxxjjj/ )()(2/12/1ffUf 需要多次矩阵运算，需要多次矩阵运算，计算量大计算量大 守恒性好，耗散小，守恒性好，耗散小，数值解质量好数值解质量好5Copyright by Li XinliangkjjkUSV2/12/1 通量分裂通量分裂+迎风差分迎风差分 引入数值耗散引入数值耗散 分裂本身不带来耗散，分裂本身不带来耗散， 但会放大（或减少）差分的耗散但会放大（或减少）差分的耗散举例：0)(xU

9、ftUfff分裂过程),(21),(21UffUfffffxxxUftUxxx00耗散耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂后的流场越偏离原先流场，则分裂后的流场越偏离原先流场，则总体耗散越大总体耗散越大如使用低精差分度格式，如使用低精差分度格式， 则对分裂形式敏感则对分裂形式敏感 （推荐使用特征分裂）（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）（可使用逐点分裂）6Copyright by Li Xinliang概念澄清概

10、念澄清耗散放耗散放大系数大系数Copyright by Li Xinliang7 7.1 Roe格式格式 守恒型格式的范例守恒型格式的范例0)(xuftu)(,0)(ufuaxuuatu破坏守恒性破坏守恒性后果很严重（后果很严重（?）为了便于使用迎风格式、特征分裂解耦，为了便于使用迎风格式、特征分裂解耦， 通常把守恒型方程改写为非守恒型通常把守恒型方程改写为非守恒型守恒方程守恒方程+守恒格式守恒格式=解守恒解守恒方程不守恒，方程不守恒， 即使差分方法守恒，也无法即使差分方法守恒，也无法做到解守恒做到解守恒)(12/12/1jjjuuxxu02/12/1xuuatujjjj02/12/1xuua

11、utjjjjjj由于由于a非常数，无法消非常数，无法消去中间项去中间项 ！.)()()()(2/532/732/322/522/112/312/12/1uauauauauauauuajjjj不再守恒思路：思路： 保证特征方向，找回守恒性保证特征方向，找回守恒性守恒型方程优点：守恒型方程优点： 守恒性守恒性非守恒方程优点：非守恒方程优点： 清晰的特征方向清晰的特征方向兼顾守恒与非守恒方程Copyright by Li Xinliang81. 单方程的单方程的Roe格式格式0)(xuftu0)(假设,0)(ufuaxuuatu0/ )(0/ )(11axffaxffxfnjnjnjnjnj1阶迎风

12、（直接从守恒方程出发）)(212112/12/11njnjnjnnjnjuuafffj)(12/12/1njnjnjffxxf0当)(0当,)()(),(111112/1jjnjnjnjnjnjnjnjjjnjuuuauuuuufufuuaa（1）（2）体现了特征方向体现了特征方向只有这种表达式，只有这种表达式，才能保证才能保证 （2）与（）与（1）等价）等价（3）2/ )(12/1jjjaaa)2/ )(12/1jjjuuaaor都无法保证 （2） 与（1）等价。简单的线性平均不行 (非线性系统，中点的斜率不等于平均斜率）关键： 构造2/1jaRoe 格式“平均斜率平均斜率”，不等于，不等于

13、“斜率的平均值斜率的平均值”，也不等于中点处的斜率，也不等于中点处的斜率njnjnjnjuuufuf11)()(平均平均斜率斜率Copyright by Li Xinliang92. 方程组的情况方程组的情况 （Roe 格式的意义）格式的意义）0)(xtUfUUf(U)Af(U)AU, 0 xt)(12/12/1jjjxxfff)(12/1jjjU,Uff)U(U)U,(UA21)f(U)f(U21)U,(UfLRLRLRLRSS)U,(UA1LRSS)U,(UA1LR 需满足如下条件需满足如下条件 （Uniform 特性）特性）单方程的简单推广单方程的简单推广 1） 连续，且连续，且 2）

14、可通过相似变换对角化，即可通过相似变换对角化，即)(UAU)(U,ASSA1保证双曲性3）对于任意）对于任意UR,UL有有)(,()()()()(),(111111jjjjjjjjjjjjuuuuaufufuuufufuua)U)(UU,(UA)f(U)f(ULRLRLR单方程的推广，含义单方程的推广，含义为为平均增长率平均增长率标量方程向矩阵方程的简单推广，标量方程向矩阵方程的简单推广，但但 构造很困难。构造很困难。)U,(UALR)(21),(),()(21),(VUAVUAVAUAVUA均不满足Uniform特性)()(1LRjjU,UfU,Uf经常记为)U,(UALR)U,(UALR)

15、U,(UALR)U,(UALR平均斜率平均斜率Copyright by Li Xinliang103. 矩阵矩阵 的构造的构造)U,(UALR)U)(UU,(UA)f(U)f(ULRLRLR关键：关键：错误！/ )()U(U)f(U)f(U)U,(UALRLRLR“向量除以向量向量除以向量” ？直接求平均增长率：uf(u)uLuRuRoeRoe点的斜率为平均斜率点的斜率为平均斜率（根据拉格朗日中值定理，（根据拉格朗日中值定理，UL,UR区区间内肯定存在间内肯定存在Roe点）点）思路思路1： 在在UL与与UR之间寻找一个点之间寻找一个点URoe, 该点该点处的增长率为平均增长率处的增长率为平均增

16、长率f(u)=u2u二次函数二次函数 Roe点与点与中点重合中点重合标量函数的启示：标量函数的启示： Roe点肯定存在（点肯定存在（Langrage 中值定理）中值定理） 二次函数的中点即为二次函数的中点即为Roe点点思路思路2： 进行坐标变换，得到一进行坐标变换，得到一个二次（齐）函数个二次（齐）函数F(W)F(U(W)F(U)U(W)U引入如果如果 是二次（齐）函数，则其是二次（齐）函数，则其中点中点 即为即为Roe点点重要启示F(W)/2W(WWRL更准确地讲，应当是要求更准确地讲，应当是要求 为为W的线性函数，的线性函数， 即增长率为线性函数即增长率为线性函数 （中点处的增长率刚好为平

17、均增长率）（中点处的增长率刚好为平均增长率）WF(W)Copyright by Li Xinliang1122312121211wwwwwwU(W)UHuwww12/1321WpEH针对针对Euler方程的具体构造方程的具体构造引入新变量：则：目的：目的： 使得使得F(w)是是W二次齐函数二次齐函数 （增长（增长率为线性函数）率为线性函数）0 xtf(U)UTTpEupuuEu)(,()(,),(2UfUf(U)不是U的二次齐函数23223121211)(wwwwwwwWf二次齐函数！二次齐函数！)W)(WWC(ffLRLR)/2W(WWLR中点处的斜率即为平均斜率。中点处的斜率即为平均斜率。

18、Roe点Roe点为：)WU(U 23122312011210wwwwwwwwwf(W)C(W)增长率为线增长率为线性函数！性函数！Copyright by Li Xinliang12最终：)UA()U,(UALR)2)(1(22uHc其中其中 如下计算：如下计算：U平均增长率（矩阵）平均增长率（矩阵）)/()()/()(2/ )(2RLRRLLRLRRLLRLHHHuuu含义：含义：左、右两个状态点的某种平均左、右两个状态点的某种平均 （称为（称为Roe平均，为密度加权平均）平均，为密度加权平均） 该状态点对应的增长率（矩阵）为平均增长率（矩阵）该状态点对应的增长率（矩阵）为平均增长率（矩阵）

19、 实际上是一种实际上是一种“等效平均等效平均”。 效果优于简单的算数（或几何）平均。效果优于简单的算数（或几何）平均。 )/()()/()(RLRRLLRLRRLLwwwvvv三维情况下，还有其他量（如压力、温度、音速等）用这三个量计算其他量（如压力、温度、音速等）用这三个量计算pEH)21(12uHp（5）2/ )(RL简单易记：Copyright by Li Xinliang13 Roe 格式的计算步骤格式的计算步骤 （半离散）（半离散）0)(xtUfU已知已知n时刻所有网格点上的物理量，对于时刻所有网格点上的物理量，对于j点：点： 1） 令令UR=Uj+1,UL=Uj （密度、压力、速度

20、等）（密度、压力、速度等） 2） 采用采用Roe平均公式（平均公式（5）计算）计算Roe平均值平均值 3） 将将Jacobian矩阵矩阵 进行特征分解进行特征分解: 计算计算 4) 计算计算 5）计算）计算 6） 计算空间导数计算空间导数 7）时间推进，计算下一时间步的值。）时间推进，计算下一时间步的值。 j-1 j j+1U)UA(SS)UA(1与前文（第与前文（第3，4讲）的形式相同，讲）的形式相同，仅需把式中的密度、压力、速度等换仅需把式中的密度、压力、速度等换成经过成经过Roe平均的密度、压力、速度平均的密度、压力、速度即可即可SS1,SS)U,(UA1LR)U(U)U,(UA21)f

21、(U)f(U21)U,(UffLRLRLRLR2/1j)(1)(1/2j1/2jffUfxxj其中：),(321diagCopyright by Li Xinliang14),(321diag可能出现导数不连续，可能出现导数不连续， 可能引起数值振荡可能引起数值振荡实际使用时实际使用时 可用如下函数代替可用如下函数代替 所谓所谓“熵修正熵修正”当当2/ )()(22fk)(f实际上是在特征值实际上是在特征值0点周围增加了耗散点周围增加了耗散Roe 格式的优点：格式的优点： 1） 保持守恒性的同时，严格保证了特征方向保持守恒性的同时，严格保证了特征方向 2） 便于推广到高精度格式便于推广到高精度

22、格式 特征投影分裂中使用特征投影分裂中使用Roe平均即可平均即可 （见（见本本PPT 第第5页）。推广到高阶后，虽不再保证严格的特征方向，页）。推广到高阶后，虽不再保证严格的特征方向， 但仍优于采但仍优于采用算数平均方法。用算数平均方法。 Roe 格式的不足：格式的不足： 本身精度只有一阶；本身精度只有一阶； 推广到高阶后，特征方向无法严格保证推广到高阶后，特征方向无法严格保证 ； 推广到二维或三维后，特征方向无法严格保证，出现振荡。推广到二维或三维后，特征方向无法严格保证，出现振荡。Copyright by Li Xinliang15作业作业 7.1 0 xtf(U)U对于一维对于一维Eul

23、er方程：方程：TTpEupuuEu)(,()(,),(2UfU引入新变量：引入新变量：Huwww12/1321W推导出推导出 及其及其Jacobian矩阵矩阵 的具体表达式的具体表达式（以（以W为自变量），并证明对于任意为自变量），并证明对于任意 ，有：，有： f(U(W)f(w) Wf(w)C(W)W)(W2WWC()f(W)f(WLRLRLRLRW,W提示：提示： 写出写出 表达式后，将向量表达式后，将向量 分别代入上式左、右两端，容分别代入上式左、右两端，容易证明相等。易证明相等。 要求：推导过程要详细，切勿简单从书本上摘抄。要求：推导过程要详细，切勿简单从书本上摘抄。C(W)f(w)

24、,LRW,W重要的重要的CFD基本功练习，一定要重视！基本功练习，一定要重视！pEH针对如下针对如下Sod 激波管问题激波管问题0)()(0)()(0)(2xpuEutExputuxut5 . 01 . 0 ,125. 0 , 05 . 01 , 1 , 0),(:0 xxput 1 , 0 x 用用Roe格式格式计算其数值解，画出计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。行比较（要求数值解与精确解画在同一张图上，便于比较）。 要求：要求： 1） 空间网格数空间网格数100

25、， 时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。 2） 尝试使用熵修正与不使用熵修正两种情况（见本尝试使用熵修正与不使用熵修正两种情况（见本PPT 15页）页） 3） 欢迎与其他数值方法得到的结果对比（最好画在同一张图上，便于比较）。欢迎与其他数值方法得到的结果对比（最好画在同一张图上，便于比较）。 16Copyright by Li Xinliang作业作业 7.2Copyright by Li Xinliang17 7.2 非物理振荡及非物理振荡及TVD格式格式1. 数值解中的非物理振荡数值解中的非物理振荡 间断附近间断附近非物理振非物理振荡的

26、根源荡的根源理论理论1：色散误差导致各波传：色散误差导致各波传播速度不同播速度不同 （第（第4讲）讲）理论理论2：物理粘性的错误计算：物理粘性的错误计算思路：思路： 物理问题物理问题 有粘；有粘； 物理粘性足以克服本身振荡物理粘性足以克服本身振荡 数值方法错误计算了物理粘性数值方法错误计算了物理粘性不足以克服振荡不足以克服振荡物理问题本身也可能振荡。但如物理问题本身也可能振荡。但如果错误计算物理粘性，则会错误果错误计算物理粘性，则会错误地加剧（或衰减）振荡。地加剧（或衰减）振荡。1） 非物理振荡的原因分析非物理振荡的原因分析Copyright by Li Xinliang18数值实验数值实验)

27、2(Re12112111njnjnjnjnjnjnjuuuxxuutuu 二阶中心差分二阶中心差分5 . 0 05 . 0 1)0 ,(xxxu当当 Re122xuxutu计算域计算域0,1,网格点网格点201 （ x=0.005 ）时间步长时间步长 x=0.0005 T=0.1时刻的u分布Re=200 x=0.005 现象：现象： x一定时，减小一定时，减小Reynolds数可抑制振荡数可抑制振荡 Reynolds数一定时，减小数一定时，减小 x可抑制振荡可抑制振荡暗示xRe是某一特征量Re=2000 x=0.005 Re=2000 x=0.0005 相同Copyright by Li Xi

28、nliang19 Re122xuxutu对流对流-扩散方程的特性：扩散方程的特性：nn+1) 1 (1knkjknjuau差分方程：差分方程：某点的值是上一时刻周围几个点上值的线性组合物理上要求系数物理上要求系数 ak 均非负均非负含义： 某处浓度的增加对下一时刻周围浓度的影响为正。j-2 j-1 j j+1 j+2差分方程单调性（无振荡）条件：差分方程单调性（无振荡）条件： 差分方程差分方程 （1）中的系数非负）中的系数非负)2(Re12112111njnjnjnjnjnjnjuuuxxuutuunjnjnjnjuxuxuxu111)21Re1()Re21 ()21Re1(xt /2ReRe

29、xx网格网格Reynolds数数xxtRe211xtnjnjnjnjuauauau2211221.Copyright by Li Xinliang20 xxReRe2） 重要概念：重要概念： 网格网格 Reynolds数数 以网格尺度度量的以网格尺度度量的Reynolds数数含义：含义： 数值振荡数值振荡 流动尺度为网格尺度流动尺度为网格尺度 网格网格 Reynolds数小，该尺度的能量数小，该尺度的能量被耗散掉被耗散掉 不发生振荡不发生振荡jj+1j-12ReRexx过于苛刻的条件过于苛刻的条件6102Re单方向网格点数单方向网格点数106， 三维三维1018 单纯靠物理粘性抑制振荡，网格间

30、距必须足够小，通常难以实现单纯靠物理粘性抑制振荡，网格间距必须足够小，通常难以实现网格足够小：不会发生振荡；网格足够小：不会发生振荡； 网格小于激波的实际厚度，则不会振荡网格小于激波的实际厚度，则不会振荡网格网格Reynolds数足够小时，物理粘性发挥作用，抑制振荡数足够小时，物理粘性发挥作用，抑制振荡Copyright by Li Xinliang213） 人工粘性人工粘性 物理粘性物理粘性 足够小才发挥作用，足够小才发挥作用， Reynolds数很高时很难做到数很高时很难做到 xxReRe思路：思路： 人为增加粘性系数人为增加粘性系数 （添加人工粘性）（添加人工粘性） 抑制振荡抑制振荡优点

31、：方法简便，有抑制振荡效果优点：方法简便，有抑制振荡效果缺点：改变了物理问题，带来误差缺点：改变了物理问题，带来误差湍流、分离流等对粘性敏感： 非物理解分离流 对粘性敏感转捩对粘性敏感很难计算对粘性敏感的问题改进措施：改进措施： A: 局部施加人工粘性局部施加人工粘性 B: 高阶人工粘性高阶人工粘性xUxUxIxb21)(Von NeumannxUxppcuxxa2232|MacCormackCopyright by Li Xinliang224) 数值振荡的定量描述数值振荡的定量描述 总变差总变差对于离散函数对于离散函数uj 定义总变差定义总变差:jjjuuTV1单调函数振荡函数j=1j=N

32、NuuTV1NuuTV1含义：含义： 反映了振荡的剧烈程度反映了振荡的剧烈程度0)(xuftu双曲型守恒方程双曲型守恒方程特点： 沿特征线 ， u不变adtdx/ufa特征线未相交总变差不变特征线相交 总变差减小结论：结论： 单个双曲型方程，总变差不增单个双曲型方程，总变差不增（Total Variation Diminishing： TVD）Copyright by Li Xinliang232 概念：概念： 单调格式、保单调格式与单调格式、保单调格式与TVD格式格式n时刻： 单调函数j=1j=Nn+1时刻： 仍是单调函数j=1j=N设设n时刻时刻 是单调的，如果是单调的，如果n+1时刻的解

33、时刻的解 仍保证单调，则称该格式为仍保证单调，则称该格式为保单调格式保单调格式。保单调格式保单调格式 nju1njuknkjknjuau10ka基本结论：基本结论： 常系数的单调格式只能是一阶常系数的单调格式只能是一阶 （Godunov） 单调格式必是保单调的；单调格式必是保单调的； 线性格式，单调与保单调等价线性格式，单调与保单调等价格式：格式：如果满足如果满足 则称其为单调格式则称其为单调格式。 ).,(, 11nqjnpjnpjnjuuuGu单调格式：单调格式：kjuG 单调格式单调格式保单调格式：保单调格式：TVD格式格式 总变差不增总变差不增nnTVTV1jjjuuTV1TVD保单调

34、保单调单调单调Copyright by Li Xinliang243. TVD格式的理论基础格式的理论基础 Harten定理定理Harten定理：定理：如果差分格式可写成如下形式：如果差分格式可写成如下形式：)()(1112121njnjnjnjnjnjnjnjuuDuuCuu01kknkjknjauau且且10 , 021212121njnjnjnjD, C DC则格式（则格式（1）是）是TVD格式格式（1）)()(1112121njnjnjnjnjnjnjnjuuDuuCuunjnjnjnjnjnjnjnjuDCuCuDu)1 (21212121111可验证：可验证： Roe格式是格式是T

35、VD格式格式)(12/12/1njnjnjffxxf)(212112/12/11njnjnjnnjnjuuafffj0当)(0当,)()(),(111112/1jjnjnjnjnjnjnjnjjjnjuuuauuuuufufuuaaauuf)(保证保证“系数非负系数非负”含义：含义：“单调格式必是单调格式必是TVD格式格式”Copyright by Li Xinliang25例例7.2.1：auuf)(0)(xuftu考虑线性单波方程考虑线性单波方程：0a)2(221122111njnjnjnjnjnjnjuuuxtaxuuatuu试讨论如下试讨论如下 Lax-Wendroff 格式格式二阶中

36、心人工粘性是否满足是否满足Harten条件条件单调格式 只有一阶精度)()(1112121njnjnjnjnjnjnjnjuuDuuCuu)2(221122111njnjnjnjnjnjnjuuuxtaxuuatuu)(21),(21222/1222/1aaDaaCnjnjxt /10 , 021212121njnjnjnjD, C DC对比条件：对比条件：211222xuuuxunjnjnjnj不满足不满足Harten 条件条件Copyright by Li Xinliang26知识回顾：知识回顾： Lax-Wendroff 格式格式)2(221122111njnjnjnjnjnjnjuuu

37、xtaxuuatuuTaylor展开，展开，写出修正方程写出修正方程),()24121()61(6121332222xtOxuutaxuuatutuuxxxxxxxxxxtttttt0 xuatuxxttxtuauauu2),(66124332222xtOxuatuxutaauuxxxtttxxxxxt时时-空二阶精度空二阶精度巧妙添加人工粘性，不但克服了不稳定性，而且抵消了时间误差，提高了时间巧妙添加人工粘性，不但克服了不稳定性，而且抵消了时间误差，提高了时间精度精度类似方法：类似方法： Beam-Warming 格式格式211222xuuuxunjnjnjnj人工粘性)2(22341222

38、121njnjnjnjnjnjnjnjuuuxtaxuuuatuu二阶精度迎风差分二阶精度迎风差分人工粘性，且提高时间精度特点：特点： 全离全离散、时刻耦合散、时刻耦合Copyright by Li Xinliang274. 构建构建TVD 格式格式思路：思路： 对现有格式进行改造，使之符合对现有格式进行改造，使之符合Harten条件条件通常在通常在Roe、L-W、B-M （或其组合）基础上改进（或其组合）基础上改进80年代初、这些格式是主流年代初、这些格式是主流0 xuatu)2(221122111njnjnjnjnjnjnjuuuxtaxuuatuu（1） 以以 L-W格式为基础改造的格式

39、格式为基础改造的格式L-W原格式（原格式（2阶）阶） = 1阶迎风阶迎风+ 修正项修正项 新格式新格式= 1阶迎风阶迎风+ 限制函数限制函数*修正项修正项)( xO )2(2)()(21112111njnjnjnjnjnjnjuuuauuauuxt)(2)1 ()(111jjjjnjnjuuaauuauu1jjjuuu)()(1112121njnjnjnjnjnjnjnjuuDuuCuu引入限制函数引入限制函数 （限制器）（限制器）)(2)1 ()(111jjjjnjnjuuaauuauuj1阶迎风部分修正项jjjjjjjuuuurr11 ),(0aCopyright by Li Xinlia

40、ng28)(2)1 ()(111jjjjnjnjuuaauuauujjjjjjjjuuuurr11 ),(显然显然 格式为格式为LW （2 阶）阶）可验证：可验证： 格式为格式为B-M （2阶）阶） 1 新格式新格式= 1阶迎风阶迎风+ *（LW格式格式-1阶迎风）阶迎风）新格式：r)(2)1 ()(111jjjjnjnjuuaauuauurLW, BM均为线性格式，均为线性格式，二者组合仍为二阶二者组合仍为二阶根据Harten定理，可知2)()(1jjjrrr时，可满足TVD性质(2) 精度条件精度条件Beam-Warming二阶精度区二阶精度区TVD区区二阶精度二阶精度TVD区（二区（二者

41、交集）者交集）Copyright by Li Xinliang29 7.3 WENO格式格式 高精度的激波捕捉法高精度的激波捕捉法 1. 基本思路基本思路0 xuaxu0axu26154132231jjjjjjjfafafafafafaf j-3, j-2,j-1,j,j+1,j+2j-3,j-2,j-1,j； j-2,j-1,j,j+1； j-1,j,j+1,j+2五个基架点被分成三个组五个基架点被分成三个组1） 若高精度逼近若高精度逼近 ， 必然利用多个基架点必然利用多个基架点 2） 如果该基架点内函数有间断，会导致振荡如果该基架点内函数有间断，会导致振荡 3） 间断不可能处处存在间断不可

42、能处处存在 4） 把基架点分成多个组（模板），把基架点分成多个组（模板）， 每个模板每个模板独立独立计算计算j点导数的逼近。点导数的逼近。 得到得到多个差分多个差分 5）根据每个模板的光滑程度，设定权重）根据每个模板的光滑程度，设定权重 6） 对多个差分结果进行加权平均对多个差分结果进行加权平均 。光滑度越高，。光滑度越高，权重越大。如果某模板存在间断，则权重趋于权重越大。如果某模板存在间断，则权重趋于0；如果都光滑，则组合成更高阶格式。如果都光滑，则组合成更高阶格式。Copyright by Li Xinliang302. WENO格式的原理描述格式的原理描述) 1 (0 xuaxu0a考虑

43、线性单波方程：考虑线性单波方程：注：注： 为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采用下一节介绍的守恒形式用下一节介绍的守恒形式（1） 确定网格基架点：确定网格基架点： 6个点个点 j-3 , j-2,j-1,j,j+1,j+2 构造出该基架点上的目标差分格式构造出该基架点上的目标差分格式jxu计算这这6个点可构造个点可构造5阶迎风差分：阶迎风差分：)2(26154132231jjjjjjjuauauauauauau该格式为该格式为WENO 的的“目标目标”格式，格式，即，即， 光滑区光滑区WENO 逼近于该格式。逼近于该格

44、式。利用利用Taylor展开，可唯一确定系数展开，可唯一确定系数 （可利用小程序（可利用小程序coeff-schemes.f )实际上，还可利用分辨率优化技术，可构造出新的目标实际上，还可利用分辨率优化技术，可构造出新的目标格式（降低精度、提高分辨率，见第格式（降低精度、提高分辨率，见第4讲）。目前大量讲）。目前大量WENO的优化版做这种工作。的优化版做这种工作。)2()60/()3302060152(21123xuuuuuuujjjjjjjCopyright by Li Xinliang31(2) 将这将这6个基架点分割成个基架点分割成3个组（称为模板）个组（称为模板） 每个组独立计算每个组

45、独立计算 的差分逼近的差分逼近 模板模板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阶精度差分格式阶精度差分格式jujjjjjuauauauau)1(41)1(32)1(23)1(1)1(1)2(4)2(31)2(22)2(1)2(jjjjjuauauauau2)3(41)3(3)3(21)3(1)3(jjjjjuauauauau计算计算j点的导数点的导数u, 竟然算出了三个不同的值，怎么办？竟然算出了三个不同的值

46、，怎么办？ENO 方法：方法： 选择最优（最光滑）的，舍弃其余两个选择最优（最光滑）的，舍弃其余两个WENO的处理方法：的处理方法： 三个都要，加权平均它们。三个都要，加权平均它们。ju利用Taylor展开式，可唯一确定这些系数）（可利用小程序coeff-schemes.f )也可运用优化技术，降低精度、提高分辨率)6/()111892(123)1(xuuuuujjjjjCopyright by Li Xinliang32（3） 对这对这3个差分值进行加权平均，得到总的差分值个差分值进行加权平均，得到总的差分值)2(26154132231jjjjjjjuauauauauauaujjjjjuau

47、auauau)1(41)1(32)1(23)1(1)1(1)2(4)2(31)2(22)2(1)2(jjjjjuauauauau2)3(41)3(3)3(21)3(1)3(jjjjjuauauauau)3(3)2(2)1(1jjjjuuuu原则：原则： A. 模板内函数越光滑，则权重越大；模板内函数越光滑，则权重越大； 模板内有模板内有间断时，权重趋于间断时，权重趋于0 B. 三个模拟内函数都光滑时，这三个模拟内函数都光滑时，这三个三阶精度三个三阶精度的逼近式可组合成的逼近式可组合成一个五阶精度一个五阶精度的逼近式。的逼近式。“理想权重”（3.1） 确定理想权重确定理想权重)3(3)2(2)1

48、(1jjjjuCuCuCu令：5阶精度容易解出：10/3,10/6,10/1321CCCCopyright by Li Xinliang33（3.2） 度量每个模板内函数的光滑程度度量每个模板内函数的光滑程度 )()()(kjkufIS IS越大，表示越不光滑。越大，表示越不光滑。 光滑区，不同模板上的光滑区，不同模板上的IS趋近同一值。趋近同一值。 具体形式见下一节。具体形式见下一节。 （3.3） 给出实际权重给出实际权重321kk)3(3)2(2)1(1jjjjuuuupkkkISC)( 构造构造IS方法很多，方法很多， 例如：例如： )(kju ：第k个模板)()254(12221232

49、)1(xOxuuuuuxISjjjjj)()22(1222112)2(xOxuuuuxISjjjj光滑区逼近光滑区逼近 O（1）量级）量级间断区间断区 量级，很大量级，很大jxu2221x特点：特点： 间断区权重很小间断区权重很小 光滑区，趋近于理想权重光滑区，趋近于理想权重（3.4） 给出最终的差分逼近给出最终的差分逼近Copyright by Li Xinliang343. Jiang & Shu 的五阶的五阶WENO格式格式) 1 (0,)(,0)(aauufxufxu 守恒型；目守恒型；目 前使用的前使用的WENO格式均为守恒型格式均为守恒型针对方程：模板模板1模板模板2 模板

50、模板3构造差分格式如下：xffxfWENOWENOjjWENOj/ )(2/12/1)3(2/13)2(2/12)1 (2/112/1jjjWENOffffjjjjjfffH6/116/73/112)1(2/111)2(2/13/16/56/1jjjjfffH21)3(2/16/16/53/1jjjjfffH构造方法与前文相同构造方法与前文相同 （但注意这里构造的是通量，（但注意这里构造的是通量，而前文是直接构造差分格式）而前文是直接构造差分格式）针对整个网格基，构造出针对整个网格基，构造出5阶精度的通量（理想情阶精度的通量（理想情况下的通量）况下的通量）并构造出每个模板上的通量，计算出理想权

51、重。并构造出每个模板上的通量，计算出理想权重。)32747132(60121122/1jjjjjjffffff321kk3 , 2 , 1,)(kISCpkkk10/3,10/6,10/1321CCC仍利用程序仍利用程序coeff-schemes.f求系数求系数理想权重理想权重光滑度量因子实际权重Copyright by Li Xinliang35光滑度量因子的计算光滑度量因子的计算 （Jiang & Shu）21212)(2/12/1lklllxxkdxqxxISjj2/12/12/12/122232)()(jjjjxxkxxkkdxxqxxdxxqxxIS)(2)()(21)()(

52、kjjkjjjkfxxfxxfxq k=1 k=2 k=3其中：其中：j-2 j-1 j j+1 j+2 是使用模板是使用模板k 得到的插值函数得到的插值函数 )(kjx)(xqk利用利用j-2,j-1,j点上的点上的值构造的插值函数值构造的插值函数 )(1xq)(1xqx2122121)2(1213)34(41jjjjjjffffffIS2112112)2(1213)(41jjjjjfffffIS2212213)2(1213)43(41jjjjjjffffffIS，2)(22)()(1213)(kjkjkfxf xIS 特点：特点： 光滑区趋近同一个值光滑区趋近同一个值 非光滑区值远大于光滑

53、区非光滑区值远大于光滑区 O(1)222)(1213)(jjkfxf xIS j 点一阶、二阶导数的差点一阶、二阶导数的差分逼近（用模板分逼近（用模板k计算）计算）)(2xO 代入Copyright by Li Xinliang36最终最终5阶阶 WENO 格式为格式为) 1 (0,)(,0)(aauufxufxuxfffWENOWENOjWENOjj/ )(2/12/1)3(2/13)2(2/12)1 (2/112/1jjjWENOffffj2122121)2(1213)34(41jjjjjjffffffIS2112112)2(1213)(41jjjjjfffffIS2212213)2(12

54、13)43(41jjjjjjffffffIS321kk610, 2p10/3,10/6,10/1321CCC3 , 2 , 1,)(kISCpkkk0,)(,0)(aauufxufxu 正通量情况（正通量情况（ a0）)3(2/13)2(2/12) 1 (2/112/1jjjWENOffffj 负通量情况负通量情况 （a0）的方法计算）的方法计算xf采用针对负通量（采用针对负通量（a0）的方法计算）的方法计算可采用可采用Steger-Warming, L-F, Van Leer 等分裂，等分裂， 见第见第4讲讲)(12/12/1jjjxxfff)(2/12/112/12/1WENOjWENOj

55、jjVVSf)2, 1, 1, 2(jjjjjkk1/2j1/2jkUSV 1) 计算计算1/2j1/2j1S,S1/2j)A(USS1/2j1/2j1/2j11/2j，它们是，它们是 的函数（推荐使用的函数（推荐使用Roe平均计算平均计算 ） 2/1jUVVWENOj2/1VWENOj2/1V效果更好效果更好但计算量较大但计算量较大计算量小计算量小但效果略差但效果略差有轻微振荡有轻微振荡2/1jUx0.60.620.640.660.60.22x00.80.81arithmatic averageRoe averageSteger-Warming FVS数值测试：数值测试： 1维维Sod激波管问题，网格点激波管问题，网格点128 差分方法差分方法: 7阶精度经典阶精度经典WENO (Jiang & Shu)分裂方式：分裂方式： Steger-Warming FVS; 特征投影分裂（算术平均特征投影分裂（算术平均/Roe平均）平均）t=0.1 时刻密度及速度分布时刻密度及速度分布xu00.80.81arithmatic