第7講-差分方法3

1、Copyright by Li Xinliang1知識回顧知識回顧1） 守恒型方程與守恒型格式守恒型方程與守恒型格式守恒型方程：守恒型方程： 散度型散度型0)(xuftu0)(UFUt守恒型格式：守恒型格式： 差量型差量型1jjjxf習慣寫為習慣寫為)(12/12/1jjjffxxf僅為記號，僅為記號，與與j+1/2點點上的值無關上的值無關！守恒型方程守恒型方程+ 守恒型格式守恒型格式=解守恒解守恒“解守恒解守恒”： 數(shù)值解的積分誤差為數(shù)值解的積分誤差為0 （如果邊界準確）（如果邊界準確） 保證總量（總質量、總動量、總能量）嚴格守恒保證總量（總質量、總動量、總能量）嚴格守恒 （無誤差）（無誤差

2、）0)(12/12/1jjjffxtu0)(10)(12/12/12/12/1Njjjjjjjffxutffxtu邊界點早期的早期的CFD:極為重視守恒性極為重視守恒性；目前目前CFD: 仍很重視守恒性仍很重視守恒性難點復雜非線性系統(tǒng)的守恒性很難保證.)()()()(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陣特征值陣特征值純正或純負純正或純負優(yōu)點：優(yōu)點： 無需矩陣運算，計算量小，使用方便無需矩陣運算，計算量小，使用方便不足：不足： 僅重新組合，沒有做到真正解耦。僅重新組合，沒有做到真正解耦。 原因：原因：xxxUAUAf具體方法：具體方法： Steger-Warming L-F Van Leer2)(2/122kkk2*kk2kkkorCopyright by Li Xinliang3方法方法2） 特征（投影）分裂特征（投影）分裂)()(jjjjjjjxxxxxVVSUSSUSSUA1j1jj1j

4、00 xtxtUAUf(U)U0jjjxtUAU在網(wǎng)格基上凍結系數(shù) j-2 j-1 j j+1 在基架點上系數(shù)在基架點上系數(shù) 不變不變jjxUAjAUSVj)(12/12/1jjjxxVVV 優(yōu)點：特征分解，（局部）解耦優(yōu)點：特征分解，（局部）解耦耗散小，數(shù)值振蕩低耗散小，數(shù)值振蕩低 缺點：大量矩陣運算，計算量大缺點：大量矩陣運算，計算量大 每計算一個點的導數(shù)，要進行每計算一個點的導數(shù)，要進行m次矩次矩陣運算陣運算 （m為網(wǎng)格基點數(shù)）為網(wǎng)格基點數(shù)）原理描述原理描述 （非守恒，很少采用；（非守恒，很少采用； 實際上使用下一頁的方法）實際上使用下一頁的方法） 守恒型方式守恒型方式計算計算x)(Uf

5、xxjjj/ )()(2/12/1ffUf j-2 j-1 j j+1 在基架點上系數(shù) 不變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） 在網(wǎng)格基上計算在網(wǎng)格基上計算) 1, 1.(2/12/1jjjkegkjjkUSV j-2 j-1 j j+1 計算fj+1/2用到的點注意，在該網(wǎng)格基上（例如注意，在該網(wǎng)格基上（例如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步驟的算法描述步驟的算法描述 （注意：（注意： 實際上是兩重循環(huán)）實際上是兩重循環(huán)） do j=1,N do k=j-1,j+1 (網(wǎng)格基，可以

8、是更多或更少點) enddo enddo do j=1,N enddo),(),(112/1112/1jjjjjjjjffVV,VVVV,VV)(2/12/112/12/1jjjjSVVfxxjjj/ )()(2/12/1ffUf 需要多次矩陣運算，需要多次矩陣運算，計算量大計算量大 守恒性好，耗散小，守恒性好，耗散小，數(shù)值解質量好數(shù)值解質量好5Copyright by Li XinliangkjjkUSV2/12/1 通量分裂通量分裂+迎風差分迎風差分 引入數(shù)值耗散引入數(shù)值耗散 分裂本身不帶來耗散，分裂本身不帶來耗散， 但會放大（或減少）差分的耗散但會放大（或減少）差分的耗散舉例：0)(xU

9、ftUfff分裂過程),(21),(21UffUfffffxxxUftUxxx00耗散耗散如果差分格式無耗散（例如都用中心差分），則通量分裂不帶來耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂后的流場越偏離原先流場，則分裂后的流場越偏離原先流場，則總體耗散越大總體耗散越大如使用低精差分度格式，如使用低精差分度格式， 則對分裂形式敏感則對分裂形式敏感 （推薦使用特征分裂）（推薦使用特征分裂）如使用高精度格式（低耗散），則對分裂形式不敏感如使用高精度格式（低耗散），則對分裂形式不敏感 （可使用逐點分裂）（可使用逐點分裂）6Copyright by Li Xinliang概念澄清概

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

11、utjjjjjj由于由于a非常數(shù)，無法消非常數(shù)，無法消去中間項去中間項 ！.)()()()(2/532/732/322/522/112/312/12/1uauauauauauauuajjjj不再守恒思路：思路： 保證特征方向，找回守恒性保證特征方向，找回守恒性守恒型方程優(yōu)點：守恒型方程優(yōu)點： 守恒性守恒性非守恒方程優(yōu)點：非守恒方程優(yōu)點： 清晰的特征方向清晰的特征方向兼顧守恒與非守恒方程Copyright by Li Xinliang81. 單方程的單方程的Roe格式格式0)(xuftu0)(假設,0)(ufuaxuuatu0/ )(0/ )(11axffaxffxfnjnjnjnjnj1階迎風

12、（直接從守恒方程出發(fā)）)(212112/12/11njnjnjnnjnjuuafffj)(12/12/1njnjnjffxxf0當)(0當,)()(),(111112/1jjnjnjnjnjnjnjnjjjnjuuuauuuuufufuuaa（1）（2）體現(xiàn)了特征方向體現(xiàn)了特征方向只有這種表達式，只有這種表達式，才能保證才能保證 （2）與（）與（1）等價）等價（3）2/ )(12/1jjjaaa)2/ )(12/1jjjuuaaor都無法保證 （2） 與（1）等價。簡單的線性平均不行 (非線性系統(tǒng)，中點的斜率不等于平均斜率）關鍵： 構造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） 連續(xù)，且連續(xù)，且 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經(jīng)常記為)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點的斜率為平均斜率點的斜率為平均斜率（根據(jù)拉格朗日中值定理，（根據(jù)拉格朗日中值定理，UL,UR區(qū)區(qū)間內(nèi)肯定存在間內(nèi)肯定存在Roe點）點）思路思路1： 在在UL與與UR之間尋找一個點之間尋找一個點URoe, 該點該點處的增長率為平均增長率處的增長率為平均增

16、長率f(u)=u2u二次函數(shù)二次函數(shù) Roe點與點與中點重合中點重合標量函數(shù)的啟示：標量函數(shù)的啟示： Roe點肯定存在（點肯定存在（Langrage 中值定理）中值定理） 二次函數(shù)的中點即為二次函數(shù)的中點即為Roe點點思路思路2： 進行坐標變換，得到一進行坐標變換，得到一個二次（齊）函數(shù)個二次（齊）函數(shù)F(W)F(U(W)F(U)U(W)U引入如果如果 是二次（齊）函數(shù)，則其是二次（齊）函數(shù)，則其中點中點 即為即為Roe點點重要啟示F(W)/2W(WWRL更準確地講，應當是要求更準確地講，應當是要求 為為W的線性函數(shù)，的線性函數(shù)， 即增長率為線性函數(shù)即增長率為線性函數(shù) （中點處的增長率剛好為平

17、均增長率）（中點處的增長率剛好為平均增長率）WF(W)Copyright by Li Xinliang1122312121211wwwwwwU(W)UHuwww12/1321WpEH針對針對Euler方程的具體構造方程的具體構造引入新變量：則：目的：目的： 使得使得F(w)是是W二次齊函數(shù)二次齊函數(shù) （增長（增長率為線性函數(shù)）率為線性函數(shù)）0 xtf(U)UTTpEupuuEu)(,()(,),(2UfUf(U)不是U的二次齊函數(shù)23223121211)(wwwwwwwWf二次齊函數(shù)！二次齊函數(shù)！)W)(WWC(ffLRLR)/2W(WWLR中點處的斜率即為平均斜率。中點處的斜率即為平均斜率。

18、Roe點Roe點為：)WU(U 23122312011210wwwwwwwwwf(W)C(W)增長率為線增長率為線性函數(shù)！性函數(shù)！Copyright by Li Xinliang12最終：)UA()U,(UALR)2)(1(22uHc其中其中 如下計算：如下計算：U平均增長率（矩陣）平均增長率（矩陣）)/()()/()(2/ )(2RLRRLLRLRRLLRLHHHuuu含義：含義：左、右兩個狀態(tài)點的某種平均左、右兩個狀態(tài)點的某種平均 （稱為（稱為Roe平均，為密度加權平均）平均，為密度加權平均） 該狀態(tài)點對應的增長率（矩陣）為平均增長率（矩陣）該狀態(tài)點對應的增長率（矩陣）為平均增長率（矩陣）

19、 實際上是一種實際上是一種“等效平均等效平均”。 效果優(yōu)于簡單的算數(shù)（或幾何）平均。效果優(yōu)于簡單的算數(shù)（或幾何）平均。 )/()()/()(RLRRLLRLRRLLwwwvvv三維情況下，還有其他量（如壓力、溫度、音速等）用這三個量計算其他量（如壓力、溫度、音速等）用這三個量計算pEH)21(12uHp（5）2/ )(RL簡單易記：Copyright by Li Xinliang13 Roe 格式的計算步驟格式的計算步驟 （半離散）（半離散）0)(xtUfU已知已知n時刻所有網(wǎng)格點上的物理量，對于時刻所有網(wǎng)格點上的物理量，對于j點：點： 1） 令令UR=Uj+1,UL=Uj （密度、壓力、速度

20、等）（密度、壓力、速度等） 2） 采用采用Roe平均公式（平均公式（5）計算）計算Roe平均值平均值 3） 將將Jacobian矩陣矩陣 進行特征分解進行特征分解: 計算計算 4) 計算計算 5）計算）計算 6） 計算空間導數(shù)計算空間導數(shù) 7）時間推進，計算下一時間步的值。）時間推進，計算下一時間步的值。 j-1 j j+1U)UA(SS)UA(1與前文（第與前文（第3，4講）的形式相同，講）的形式相同，僅需把式中的密度、壓力、速度等換僅需把式中的密度、壓力、速度等換成經(jīng)過成經(jīng)過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可能出現(xiàn)導數(shù)不連續(xù)，可能出現(xiàn)導數(shù)不連續(xù)， 可能引起數(shù)值振蕩可能引起數(shù)值振蕩實際使用時實際使用時 可用如下函數(shù)代替可用如下函數(shù)代替 所謂所謂“熵修正熵修正”當當2/ )()(22fk)(f實際上是在特征值實際上是在特征值0點周圍增加了耗散點周圍增加了耗散Roe 格式的優(yōu)點：格式的優(yōu)點： 1） 保持守恒性的同時，嚴格保證了特征方向保持守恒性的同時，嚴格保證了特征方向 2） 便于推廣到高精度格式便于推廣到高精度

22、格式 特征投影分裂中使用特征投影分裂中使用Roe平均即可平均即可 （見（見本本PPT 第第5頁）。推廣到高階后，雖不再保證嚴格的特征方向，頁）。推廣到高階后，雖不再保證嚴格的特征方向， 但仍優(yōu)于采但仍優(yōu)于采用算數(shù)平均方法。用算數(shù)平均方法。 Roe 格式的不足：格式的不足： 本身精度只有一階；本身精度只有一階； 推廣到高階后，特征方向無法嚴格保證推廣到高階后，特征方向無法嚴格保證 ； 推廣到二維或三維后，特征方向無法嚴格保證，出現(xiàn)振蕩。推廣到二維或三維后，特征方向無法嚴格保證，出現(xiàn)振蕩。Copyright by Li Xinliang15作業(yè)作業(yè) 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格式格式計算其數(shù)值解，畫出計算其數(shù)值解，畫出t=0.14時刻密度、速度及壓力的分布；并與精確解進時刻密度、速度及壓力的分布；并與精確解進行比較（要求數(shù)值解與精確解畫在同一張圖上，便于比較）。行比較（要求數(shù)值解與精確解畫在同一張圖上，便于比較）。 要求：要求： 1） 空間網(wǎng)格數(shù)空間網(wǎng)格數(shù)100

25、， 時間推進格式選用時間推進格式選用3階階Runge-Kutta,時間步長自選。時間步長自選。 2） 嘗試使用熵修正與不使用熵修正兩種情況（見本嘗試使用熵修正與不使用熵修正兩種情況（見本PPT 15頁）頁） 3） 歡迎與其他數(shù)值方法得到的結果對比（最好畫在同一張圖上，便于比較）。歡迎與其他數(shù)值方法得到的結果對比（最好畫在同一張圖上，便于比較）。 16Copyright by Li Xinliang作業(yè)作業(yè) 7.2Copyright by Li Xinliang17 7.2 非物理振蕩及非物理振蕩及TVD格式格式1. 數(shù)值解中的非物理振蕩數(shù)值解中的非物理振蕩 間斷附近間斷附近非物理振非物理振蕩的

26、根源蕩的根源理論理論1：色散誤差導致各波傳：色散誤差導致各波傳播速度不同播速度不同 （第（第4講）講）理論理論2：物理粘性的錯誤計算：物理粘性的錯誤計算思路：思路： 物理問題物理問題 有粘；有粘； 物理粘性足以克服本身振蕩物理粘性足以克服本身振蕩 數(shù)值方法錯誤計算了物理粘性數(shù)值方法錯誤計算了物理粘性不足以克服振蕩不足以克服振蕩物理問題本身也可能振蕩。但如物理問題本身也可能振蕩。但如果錯誤計算物理粘性，則會錯誤果錯誤計算物理粘性，則會錯誤地加?。ɑ蛩p）振蕩。地加?。ɑ蛩p）振蕩。1） 非物理振蕩的原因分析非物理振蕩的原因分析Copyright by Li Xinliang18數(shù)值實驗數(shù)值實驗)

27、2(Re12112111njnjnjnjnjnjnjuuuxxuutuu 二階中心差分二階中心差分5 . 0 05 . 0 1)0 ,(xxxu當當 Re122xuxutu計算域計算域0,1,網(wǎng)格點網(wǎng)格點201 （ x=0.005 ）時間步長時間步長 x=0.0005 T=0.1時刻的u分布Re=200 x=0.005 現(xiàn)象：現(xiàn)象： x一定時，減小一定時，減小Reynolds數(shù)可抑制振蕩數(shù)可抑制振蕩 Reynolds數(shù)一定時，減小數(shù)一定時，減小 x可抑制振蕩可抑制振蕩暗示xRe是某一特征量Re=2000 x=0.005 Re=2000 x=0.0005 相同Copyright by Li Xi

28、nliang19 Re122xuxutu對流對流-擴散方程的特性：擴散方程的特性：nn+1) 1 (1knkjknjuau差分方程：差分方程：某點的值是上一時刻周圍幾個點上值的線性組合物理上要求系數(shù)物理上要求系數(shù) ak 均非負均非負含義： 某處濃度的增加對下一時刻周圍濃度的影響為正。j-2 j-1 j j+1 j+2差分方程單調性（無振蕩）條件：差分方程單調性（無振蕩）條件： 差分方程差分方程 （1）中的系數(shù)非負）中的系數(shù)非負)2(Re12112111njnjnjnjnjnjnjuuuxxuutuunjnjnjnjuxuxuxu111)21Re1()Re21 ()21Re1(xt /2ReRe

29、xx網(wǎng)格網(wǎng)格Reynolds數(shù)數(shù)xxtRe211xtnjnjnjnjuauauau2211221.Copyright by Li Xinliang20 xxReRe2） 重要概念：重要概念： 網(wǎng)格網(wǎng)格 Reynolds數(shù)數(shù) 以網(wǎng)格尺度度量的以網(wǎng)格尺度度量的Reynolds數(shù)數(shù)含義：含義： 數(shù)值振蕩數(shù)值振蕩 流動尺度為網(wǎng)格尺度流動尺度為網(wǎng)格尺度 網(wǎng)格網(wǎng)格 Reynolds數(shù)小，該尺度的能量數(shù)小，該尺度的能量被耗散掉被耗散掉 不發(fā)生振蕩不發(fā)生振蕩jj+1j-12ReRexx過于苛刻的條件過于苛刻的條件6102Re單方向網(wǎng)格點數(shù)單方向網(wǎng)格點數(shù)106， 三維三維1018 單純靠物理粘性抑制振蕩，網(wǎng)格間

30、距必須足夠小，通常難以實現(xiàn)單純靠物理粘性抑制振蕩，網(wǎng)格間距必須足夠小，通常難以實現(xiàn)網(wǎng)格足夠?。翰粫l(fā)生振蕩；網(wǎng)格足夠?。翰粫l(fā)生振蕩； 網(wǎng)格小于激波的實際厚度，則不會振蕩網(wǎng)格小于激波的實際厚度，則不會振蕩網(wǎng)格網(wǎng)格Reynolds數(shù)足夠小時，物理粘性發(fā)揮作用，抑制振蕩數(shù)足夠小時，物理粘性發(fā)揮作用，抑制振蕩Copyright by Li Xinliang213） 人工粘性人工粘性 物理粘性物理粘性 足夠小才發(fā)揮作用，足夠小才發(fā)揮作用， Reynolds數(shù)很高時很難做到數(shù)很高時很難做到 xxReRe思路：思路： 人為增加粘性系數(shù)人為增加粘性系數(shù) （添加人工粘性）（添加人工粘性） 抑制振蕩抑制振蕩優(yōu)點

31、：方法簡便，有抑制振蕩效果優(yōu)點：方法簡便，有抑制振蕩效果缺點：改變了物理問題，帶來誤差缺點：改變了物理問題，帶來誤差湍流、分離流等對粘性敏感： 非物理解分離流 對粘性敏感轉捩對粘性敏感很難計算對粘性敏感的問題改進措施：改進措施： A: 局部施加人工粘性局部施加人工粘性 B: 高階人工粘性高階人工粘性xUxUxIxb21)(Von NeumannxUxppcuxxa2232|MacCormackCopyright by Li Xinliang224) 數(shù)值振蕩的定量描述數(shù)值振蕩的定量描述 總變差總變差對于離散函數(shù)對于離散函數(shù)uj 定義總變差定義總變差:jjjuuTV1單調函數(shù)振蕩函數(shù)j=1j=N

32、NuuTV1NuuTV1含義：含義： 反映了振蕩的劇烈程度反映了振蕩的劇烈程度0)(xuftu雙曲型守恒方程雙曲型守恒方程特點： 沿特征線 ， u不變adtdx/ufa特征線未相交總變差不變特征線相交 總變差減小結論：結論： 單個雙曲型方程，總變差不增單個雙曲型方程，總變差不增（Total Variation Diminishing： TVD）Copyright by Li Xinliang232 概念：概念： 單調格式、保單調格式與單調格式、保單調格式與TVD格式格式n時刻： 單調函數(shù)j=1j=Nn+1時刻： 仍是單調函數(shù)j=1j=N設設n時刻時刻 是單調的，如果是單調的，如果n+1時刻的解

33、時刻的解 仍保證單調，則稱該格式為仍保證單調，則稱該格式為保單調格式保單調格式。保單調格式保單調格式 nju1njuknkjknjuau10ka基本結論：基本結論： 常系數(shù)的單調格式只能是一階常系數(shù)的單調格式只能是一階 （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)(保證保證“系數(shù)非負系數(shù)非負”含義：含義：“單調格式必是單調格式必是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時時-空二階精度空二階精度巧妙添加人工粘性，不但克服了不穩(wěn)定性，而且抵消了時間誤差，提高了時間巧妙添加人工粘性，不但克服了不穩(wěn)定性，而且抵消了時間誤差，提高了時間精度精度類似方法：類似方法： Beam-Warming 格式格式211222xuuuxunjnjnjnj人工粘性)2(22341222

38、121njnjnjnjnjnjnjnjuuuxtaxuuuatuu二階精度迎風差分二階精度迎風差分人工粘性，且提高時間精度特點：特點： 全離全離散、時刻耦合散、時刻耦合Copyright by Li Xinliang274. 構建構建TVD 格式格式思路：思路： 對現(xiàn)有格式進行改造，使之符合對現(xiàn)有格式進行改造，使之符合Harten條件條件通常在通常在Roe、L-W、B-M （或其組合）基礎上改進（或其組合）基礎上改進80年代初、這些格式是主流年代初、這些格式是主流0 xuatu)2(221122111njnjnjnjnjnjnjuuuxtaxuuatuu（1） 以以 L-W格式為基礎改造的格式

39、格式為基礎改造的格式L-W原格式（原格式（2階）階） = 1階迎風階迎風+ 修正項修正項 新格式新格式= 1階迎風階迎風+ 限制函數(shù)限制函數(shù)*修正項修正項)( xO )2(2)()(21112111njnjnjnjnjnjnjuuuauuauuxt)(2)1 ()(111jjjjnjnjuuaauuauu1jjjuuu)()(1112121njnjnjnjnjnjnjnjuuDuuCuu引入限制函數(shù)引入限制函數(shù) （限制器）（限制器）)(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均為線性格式，均為線性格式，二者組合仍為二階二者組合仍為二階根據(jù)Harten定理，可知2)()(1jjjrrr時，可滿足TVD性質(2) 精度條件精度條件Beam-Warming二階精度區(qū)二階精度區(qū)TVD區(qū)區(qū)二階精度二階精度TVD區(qū)（二區(qū)（二者

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） 如果該基架點內(nèi)函數(shù)有間斷，會導致振蕩如果該基架點內(nèi)函數(shù)有間斷，會導致振蕩 3） 間斷不可能處處存在間斷不可

42、能處處存在 4） 把基架點分成多個組（模板），把基架點分成多個組（模板）， 每個模板每個模板獨立獨立計算計算j點導數(shù)的逼近。點導數(shù)的逼近。 得到得到多個差分多個差分 5）根據(jù)每個模板的光滑程度，設定權重）根據(jù)每個模板的光滑程度，設定權重 6） 對多個差分結果進行加權平均對多個差分結果進行加權平均 。光滑度越高，。光滑度越高，權重越大。如果某模板存在間斷，則權重趨于權重越大。如果某模板存在間斷，則權重趨于0；如果都光滑，則組合成更高階格式。如果都光滑，則組合成更高階格式。Copyright by Li Xinliang302. WENO格式的原理描述格式的原理描述) 1 (0 xuaxu0a考慮

43、線性單波方程：考慮線性單波方程：注：注： 為了簡便，以非守恒型形式為例講授其思路，實際使用時，請采為了簡便，以非守恒型形式為例講授其思路，實際使用時，請采用下一節(jié)介紹的守恒形式用下一節(jié)介紹的守恒形式（1） 確定網(wǎng)格基架點：確定網(wǎng)格基架點： 6個點個點 j-3 , j-2,j-1,j,j+1,j+2 構造出該基架點上的目標差分格式構造出該基架點上的目標差分格式jxu計算這這6個點可構造個點可構造5階迎風差分：階迎風差分：)2(26154132231jjjjjjjuauauauauauau該格式為該格式為WENO 的的“目標目標”格式，格式，即，即， 光滑區(qū)光滑區(qū)WENO 逼近于該格式。逼近于該格

44、式。利用利用Taylor展開，可唯一確定系數(shù)展開，可唯一確定系數(shù) （可利用小程序（可利用小程序coeff-schemes.f )實際上，還可利用分辨率優(yōu)化技術，可構造出新的目標實際上，還可利用分辨率優(yōu)化技術，可構造出新的目標格式（降低精度、提高分辨率，見第格式（降低精度、提高分辨率，見第4講）。目前大量講）。目前大量WENO的優(yōu)化版做這種工作。的優(yōu)化版做這種工作。)2()60/()3302060152(21123xuuuuuuujjjjjjjCopyright by Li Xinliang31將這將這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點的導數(shù)點的導數(shù)u, 竟然算出了三個不同的值，怎么辦？竟然算出了三個不同的值，怎么辦

46、？ENO 方法：方法： 選擇最優(yōu)（最光滑）的，舍棄其余兩個選擇最優(yōu)（最光滑）的，舍棄其余兩個WENO的處理方法：的處理方法： 三個都要，加權平均它們。三個都要，加權平均它們。ju利用Taylor展開式，可唯一確定這些系數(shù)）（可利用小程序coeff-schemes.f )也可運用優(yōu)化技術，降低精度、提高分辨率)6/()111892(123)1(xuuuuujjjjjCopyright by Li Xinliang32（3） 對這對這3個差分值進行加權平均，得到總的差分值個差分值進行加權平均，得到總的差分值)2(26154132231jjjjjjjuauauauauauaujjjjjuauauau

47、au)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. 模板內(nèi)函數(shù)越光滑，則權重越大；模板內(nèi)函數(shù)越光滑，則權重越大； 模板內(nèi)有模板內(nèi)有間斷時，權重趨于間斷時，權重趨于0 B. 三個模擬內(nèi)函數(shù)都光滑時，這三個模擬內(nèi)函數(shù)都光滑時，這三個三階精度三個三階精度的逼近式可組合成的逼近式可組合成一個五階精度一個五階精度的逼近式。的逼近式。“理想權重”（3.1） 確定理想權重確定理想權重)3(3)2(2)1(1jj

48、jjuCuCuCu令：5階精度容易解出：10/3,10/6,10/1321CCCCopyright by Li Xinliang33（3.2） 度量每個模板內(nèi)函數(shù)的光滑程度度量每個模板內(nèi)函數(shù)的光滑程度 )()()(kjkufIS IS越大，表示越不光滑。越大，表示越不光滑。 光滑區(qū)，不同模板上的光滑區(qū)，不同模板上的IS趨近同一值。趨近同一值。 具體形式見下一節(jié)。具體形式見下一節(jié)。 （3.3） 給出實際權重給出實際權重321kk)3(3)2(2)1(1jjjjuuuupkkkISC)( 構造構造IS方法很多，方法很多， 例如：例如： )(kju ：第k個模板)()254(12221232)1(x

49、OxuuuuuxISjjjjj)()22(1222112)2(xOxuuuuxISjjjj光滑區(qū)逼近光滑區(qū)逼近 O（1）量級）量級間斷區(qū)間斷區(qū) 量級，很大量級，很大jxu2221x特點：特點： 間斷區(qū)權重很小間斷區(qū)權重很小 光滑區(qū)，趨近于理想權重光滑區(qū)，趨近于理想權重（3.4） 給出最終的差分逼近給出最終的差分逼近Copyright by Li Xinliang343. Jiang & Shu 的五階的五階WENO格式格式) 1 (0,)(,0)(aauufxufxu 守恒型；目守恒型；目 前使用的前使用的WENO格式均為守恒型格式均為守恒型針對方程：模板模板1模板模板2 模板模板3構

50、造差分格式如下：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構造方法與前文相同構造方法與前文相同 （但注意這里構造的是通量，（但注意這里構造的是通量，而前文是直接構造差分格式）而前文是直接構造差分格式）針對整個網(wǎng)格基，構造出針對整個網(wǎng)格基，構造出5階精度的通量（理想情階精度的通量（理想情況下的通量）況下的通量）并構造出每個模板上的通量，計算出理想權重。并

51、構造出每個模板上的通量，計算出理想權重。)32747132(60121122/1jjjjjjffffff321kk3 , 2 , 1,)(kISCpkkk10/3,10/6,10/1321CCC仍利用程序仍利用程序coeff-schemes.f求系數(shù)求系數(shù)理想權重理想權重光滑度量因子實際權重Copyright by Li Xinliang35光滑度量因子的計算光滑度量因子的計算 （Jiang & Shu）21212)(2/12/1lklllxxkdxqxxISjj2/12/12/12/122232)()(jjjjxxkxxkkdxxqxxdxxqxxIS)(2)()(21)()(kjj

52、kjjjkfxxfxxfxq k=1 k=2 k=3其中：其中：j-2 j-1 j j+1 j+2 是使用模板是使用模板k 得到的插值函數(shù)得到的插值函數(shù) )(kjx)(xqk利用利用j-2,j-1,j點上的點上的值構造的插值函數(shù)值構造的插值函數(shù) )(1xq)(1xqx2122121)2(1213)34(41jjjjjjffffffIS2112112)2(1213)(41jjjjjfffffIS2212213)2(1213)43(41jjjjjjffffffIS，2)(22)()(1213)(kjkjkfxf xIS 特點：特點： 光滑區(qū)趨近同一個值光滑區(qū)趨近同一個值 非光滑區(qū)值遠大于光滑區(qū)非光

53、滑區(qū)值遠大于光滑區(qū) O(1)222)(1213)(jjkfxf xIS j 點一階、二階導數(shù)的差點一階、二階導數(shù)的差分逼近（用模板分逼近（用模板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(1213

54、)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/1WENOjWENOjjj

55、VVSf)2, 1, 1, 2(jjjjjkk1/2j1/2jkUSV 1) 計算計算1/2j1/2j1S,S1/2j)A(USS1/2j1/2j1/2j11/2j，它們是，它們是 的函數(shù)（推薦使用的函數(shù)（推薦使用Roe平均計算平均計算 ） 2/1jUVVWENOj2/1VWENOj2/1V效果更好效果更好但計算量較大但計算量較大計算量小計算量小但效果略差但效果略差有輕微振蕩有輕微振蕩2/1jUx0.60.620.640.660.60.22x00.80.81arithmatic averageRoe averageSteger-Warming FVS數(shù)值測試：數(shù)值測試： 1維維Sod激波管問題，網(wǎng)格點激波管問題，網(wǎng)格點128 差分方法差分方法: 7階精度經(jīng)典階精度經(jīng)典WENO (Jiang & Shu)分裂方式：分裂方式： Steger-Warming FVS; 特征投影分裂（算術平均特征投影分裂（算術平均/Roe平均）平均）t=0.1 時刻密度及速度分布時刻密度及速度分布xu00.80.81arithmatic