第7講-差分方法3

CopyrightbyLiXinliang1知識回顧1）

守恒型方程與守恒型格式守恒型方程：散度型守恒型格式：差量型習(xí)慣寫為僅為記號，與j+1/2點(diǎn)上的值無關(guān)！守恒型方程+守恒型格式=解守恒“解守恒”：數(shù)值解的積分誤差為0（如果邊界準(zhǔn)確）

保證總量（總質(zhì)量、總動量、總能量）嚴(yán)格守恒（無誤差）邊界點(diǎn)早期的CFD:極為重視守恒性；目前CFD:仍很重視守恒性難點(diǎn)——復(fù)雜非線性系統(tǒng)的守恒性很難保證中間項(xiàng)全部消去，只剩兩端CopyrightbyLiXinliang2）

通量分裂——便于使用迎風(fēng)格式方法（1）：逐點(diǎn)分裂（Steger-Warming,VanLeer,L-F）原理：利用了性質(zhì)使得的Jocabian陣特征值純正或純負(fù)優(yōu)點(diǎn)：無需矩陣運(yùn)算，計(jì)算量小，使用方便不足：僅重新組合，沒有做到真正解耦。原因：具體方法：Steger-Warming

L-FVanLeerorCopyrightbyLiXinliang3方法2）特征（投影）分裂在網(wǎng)格基上凍結(jié)系數(shù)…j-2j-1jj+1…在基架點(diǎn)上系數(shù)不變優(yōu)點(diǎn)：特征分解，（局部）解耦——耗散小，數(shù)值振蕩低

缺點(diǎn)：大量矩陣運(yùn)算，計(jì)算量大每計(jì)算一個(gè)點(diǎn)的導(dǎo)數(shù)，要進(jìn)行m次矩陣運(yùn)算（m為網(wǎng)格基點(diǎn)數(shù)）原理描述（非守恒，很少采用；實(shí)際上使用下一頁的方法）

守恒型方式計(jì)算…j-2j-1jj+1…在基架點(diǎn)上系數(shù)不變具體步驟：

假設(shè)已知U，且針對模型方程（線性單波方程）

已構(gòu)造出差分格式（1）1）計(jì)算出教材130頁的公式(6.1.11-6.1.13),式中用到各變量在j+1/2的值（例如）

可使用j,j+1點(diǎn)值的算術(shù)平均（如）或Roe平均；由計(jì)算；方法很多，例如前面介紹的或4CopyrightbyLiXinliang

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

在網(wǎng)格基上計(jì)算…j-2j-1jj+1…計(jì)算fj+1/2用到的點(diǎn)注意，在該網(wǎng)格基上（例如k=j-1,j,j+1）保持不變例如：3)利用已構(gòu)造好的差分格式，計(jì)算通量4）得到總通量

5）計(jì)算差分（j點(diǎn)處）步驟的算法描述（注意：實(shí)際上是兩重循環(huán)）doj=1,Ndok=j-1,j+1(網(wǎng)格基，可以是更多或更少點(diǎn))

enddoenddodoj=1,N

enddo需要多次矩陣運(yùn)算，計(jì)算量大

守恒性好，耗散小，數(shù)值解質(zhì)量好5CopyrightbyLiXinliang

通量分裂+迎風(fēng)差分引入數(shù)值耗散

分裂本身不帶來耗散，但會放大（或減少）差分的耗散舉例：分裂過程耗散如果差分格式無耗散（例如都用中心差分），則通量分裂不帶來耗散。=+向上平移向下平移分裂后的流場越偏離原先流場，則總體耗散越大如使用低精差分度格式，則對分裂形式敏感（推薦使用特征分裂）如使用高精度格式（低耗散），則對分裂形式不敏感（可使用逐點(diǎn)分裂）6CopyrightbyLiXinliang概念澄清耗散放大系數(shù)CopyrightbyLiXinliang7§7.1Roe格式——守恒型格式的范例破壞守恒性后果很嚴(yán)重（?）為了便于使用迎風(fēng)格式、特征分裂解耦，通常把守恒型方程改寫為非守恒型守恒方程+守恒格式=解守恒方程不守恒，即使差分方法守恒，也無法做到解守恒由于a非常數(shù)，無法消去中間項(xiàng)！不再守恒思路：保證特征方向，找回守恒性守恒型方程優(yōu)點(diǎn)：守恒性非守恒方程優(yōu)點(diǎn)：清晰的特征方向兼顧守恒與非守恒方程CopyrightbyLiXinliang81.單方程的Roe格式1階迎風(fēng)（直接從守恒方程出發(fā)）（1）（2）體現(xiàn)了特征方向只有這種表達(dá)式，才能保證（2）與（1）等價(jià)（3）or都無法保證（2）與（1）等價(jià)。簡單的線性平均不行(非線性系統(tǒng)，中點(diǎn)的斜率不等于平均斜率）關(guān)鍵：構(gòu)造Roe格式“平均斜率”，不等于“斜率的平均值”，也不等于中點(diǎn)處的斜率平均斜率CopyrightbyLiXinliang92.方程組的情況（Roe格式的意義）

需滿足如下條件（Uniform特性）單方程的簡單推廣

1）連續(xù)，且2）可通過相似變換對角化，即保證雙曲性3）對于任意UR,UL有單方程的推廣，含義為平均增長率標(biāo)量方程向矩陣方程的簡單推廣，但構(gòu)造很困難。均不滿足Uniform特性平均斜率CopyrightbyLiXinliang103.矩陣的構(gòu)造關(guān)鍵：“向量除以向量”？直接求平均增長率：uf(u)uLuRuRoeRoe點(diǎn)的斜率為平均斜率（根據(jù)拉格朗日中值定理，[UL,UR]區(qū)間內(nèi)肯定存在Roe點(diǎn)）思路1：在UL與UR之間尋找一個(gè)點(diǎn)URoe,該點(diǎn)處的增長率為平均增長率f(u)=u2u二次函數(shù)——Roe點(diǎn)與中點(diǎn)重合標(biāo)量函數(shù)的啟示：Roe點(diǎn)肯定存在（Langrage中值定理）

二次函數(shù)的中點(diǎn)即為Roe點(diǎn)思路2：進(jìn)行坐標(biāo)變換，得到一個(gè)二次（齊）函數(shù)引入如果

是二次（齊）函數(shù)，則其中點(diǎn)即為Roe點(diǎn)重要啟示更準(zhǔn)確地講，應(yīng)當(dāng)是要求為W的線性函數(shù)，

即增長率為線性函數(shù)（中點(diǎn)處的增長率剛好為平均增長率）CopyrightbyLiXinliang11針對Euler方程的具體構(gòu)造引入新變量：則：目的：使得F(w)是W二次齊函數(shù)（增長率為線性函數(shù)）f(U)不是U的二次齊函數(shù)二次齊函數(shù)！中點(diǎn)處的斜率即為平均斜率。Roe點(diǎn)Roe點(diǎn)為：增長率為線性函數(shù)！CopyrightbyLiXinliang12最終：其中如下計(jì)算：平均增長率（矩陣）含義：左、右兩個(gè)狀態(tài)點(diǎn)的某種平均（稱為Roe平均，為密度加權(quán)平均）

該狀態(tài)點(diǎn)對應(yīng)的增長率（矩陣）為平均增長率（矩陣）

實(shí)際上是一種“等效平均”。效果優(yōu)于簡單的算數(shù)（或幾何）平均。

三維情況下，還有其他量（如壓力、溫度、音速等）用這三個(gè)量計(jì)算（5）簡單易記：CopyrightbyLiXinliang13Roe格式的計(jì)算步驟（半離散）已知n時(shí)刻所有網(wǎng)格點(diǎn)上的物理量，對于j點(diǎn)：1）令UR=Uj+1,UL=Uj

（密度、壓力、速度等）

2）采用Roe平均公式（5）計(jì)算Roe平均值3）將Jacobian矩陣進(jìn)行特征分解:

計(jì)算4)計(jì)算5）計(jì)算6）計(jì)算空間導(dǎo)數(shù)7）時(shí)間推進(jìn)，計(jì)算下一時(shí)間步的值。j-1jj+1與前文（第3，4講）的形式相同，僅需把式中的密度、壓力、速度等換成經(jīng)過Roe平均的密度、壓力、速度即可其中：CopyrightbyLiXinliang14可能出現(xiàn)導(dǎo)數(shù)不連續(xù)，可能引起數(shù)值振蕩實(shí)際使用時(shí)可用如下函數(shù)代替——所謂“熵修正”實(shí)際上是在特征值0點(diǎn)周圍增加了耗散Roe格式的優(yōu)點(diǎn)：1）保持守恒性的同時(shí)，嚴(yán)格保證了特征方向2）便于推廣到高精度格式——特征投影分裂中使用Roe平均即可（見本PPT第5頁）。推廣到高階后，雖不再保證嚴(yán)格的特征方向，但仍優(yōu)于采用算數(shù)平均方法。Roe格式的不足：

本身精度只有一階；

推廣到高階后，特征方向無法嚴(yán)格保證；

推廣到二維或三維后，特征方向無法嚴(yán)格保證，出現(xiàn)振蕩。CopyrightbyLiXinliang15作業(yè)7.1對于一維Euler方程：引入新變量：推導(dǎo)出及其Jacobian矩陣的具體表達(dá)式（以W為自變量），并證明對于任意，有：提示：寫出表達(dá)式后，將向量分別代入上式左、右兩端，容易證明相等。要求：推導(dǎo)過程要詳細(xì)，切勿簡單從書本上摘抄。重要的CFD基本功練習(xí)，一定要重視！針對如下Sod激波管問題

用Roe格式計(jì)算其數(shù)值解，畫出t=0.14時(shí)刻密度、速度及壓力的分布；并與精確解進(jìn)行比較（要求數(shù)值解與精確解畫在同一張圖上，便于比較）。

要求：1）空間網(wǎng)格數(shù)100，時(shí)間推進(jìn)格式選用3階Runge-Kutta,時(shí)間步長自選。2）嘗試使用熵修正與不使用熵修正兩種情況（見本PPT15頁）3）歡迎與其他數(shù)值方法得到的結(jié)果對比（最好畫在同一張圖上，便于比較）。16CopyrightbyLiXinliang作業(yè)7.2CopyrightbyLiXinliang17§7.2非物理振蕩及TVD格式1.數(shù)值解中的非物理振蕩

間斷附近非物理振蕩的根源理論1：色散誤差導(dǎo)致各波傳播速度不同（第4講）理論2：物理粘性的錯(cuò)誤計(jì)算思路：物理問題——有粘；物理粘性足以克服本身振蕩

數(shù)值方法錯(cuò)誤計(jì)算了物理粘性——不足以克服振蕩物理問題本身也可能振蕩。但如果錯(cuò)誤計(jì)算物理粘性，則會錯(cuò)誤地加劇（或衰減）振蕩。1）非物理振蕩的原因分析CopyrightbyLiXinliang18數(shù)值實(shí)驗(yàn)

二階中心差分計(jì)算域[0,1],網(wǎng)格點(diǎn)201（Dx=0.005）時(shí)間步長Dx=0.0005

T=0.1時(shí)刻的u分布Re=200Dx=0.005現(xiàn)象：Dx一定時(shí)，減小Reynolds數(shù)可抑制振蕩Reynolds數(shù)一定時(shí)，減小Dx可抑制振蕩暗示是某一特征量Re=2000Dx=0.005Re=2000Dx=0.0005相同CopyrightbyLiXinliang19對流-擴(kuò)散方程的特性：nn+1差分方程：某點(diǎn)的值是上一時(shí)刻周圍幾個(gè)點(diǎn)上值的線性組合物理上要求系數(shù)ak

均非負(fù)含義：某處濃度的增加對下一時(shí)刻周圍濃度的影響為正。j-2j-1jj+1j+2差分方程單調(diào)性（無振蕩）條件：差分方程（1）中的系數(shù)非負(fù)網(wǎng)格Reynolds數(shù)CopyrightbyLiXinliang202）重要概念：網(wǎng)格Reynolds數(shù)以網(wǎng)格尺度度量的Reynolds數(shù)含義：

數(shù)值振蕩——流動尺度為網(wǎng)格尺度

網(wǎng)格Reynolds數(shù)小，該尺度的能量被耗散掉——不發(fā)生振蕩jj+1j-1過于苛刻的條件單方向網(wǎng)格點(diǎn)數(shù)106，三維1018

單純靠物理粘性抑制振蕩，網(wǎng)格間距必須足夠小，通常難以實(shí)現(xiàn)網(wǎng)格足夠小：不會發(fā)生振蕩；網(wǎng)格小于激波的實(shí)際厚度，則不會振蕩網(wǎng)格Reynolds數(shù)足夠小時(shí)，物理粘性發(fā)揮作用，抑制振蕩CopyrightbyLiXinliang213）

人工粘性

物理粘性足夠小才發(fā)揮作用，Reynolds數(shù)很高時(shí)很難做到

思路：人為增加粘性系數(shù)（添加人工粘性）抑制振蕩優(yōu)點(diǎn)：方法簡便，有抑制振蕩效果缺點(diǎn)：改變了物理問題，帶來誤差湍流、分離流等——對粘性敏感：非物理解分離流——對粘性敏感轉(zhuǎn)捩——對粘性敏感很難計(jì)算對粘性敏感的問題改進(jìn)措施：A:局部施加人工粘性B:高階人工粘性VonNeumannMacCormackCopyrightbyLiXinliang224)數(shù)值振蕩的定量描述——總變差對于離散函數(shù){uj}定義總變差:單調(diào)函數(shù)振蕩函數(shù)j=1j=N含義：反映了振蕩的劇烈程度雙曲型守恒方程特點(diǎn)：沿特征線，u不變特征線未相交——總變差不變特征線相交——總變差減小結(jié)論：單個(gè)雙曲型方程，總變差不增（TotalVariationDiminishing：TVD）CopyrightbyLiXinliang232概念：單調(diào)格式、保單調(diào)格式與TVD格式n時(shí)刻：單調(diào)函數(shù)j=1j=Nn+1時(shí)刻：仍是單調(diào)函數(shù)j=1j=N設(shè)n時(shí)刻是單調(diào)的，如果n+1時(shí)刻的解仍保證單調(diào)，則稱該格式為保單調(diào)格式。保單調(diào)格式基本結(jié)論：常系數(shù)的單調(diào)格式只能是一階（Godunov）

單調(diào)格式必是保單調(diào)的；

線性格式，單調(diào)與保單調(diào)等價(jià)格式：如果滿足則稱其為單調(diào)格式。單調(diào)格式：單調(diào)格式保單調(diào)格式：TVD格式總變差不增TVD保單調(diào)單調(diào)CopyrightbyLiXinliang243.TVD格式的理論基礎(chǔ)——Harten定理Harten定理：如果差分格式可寫成如下形式：且則格式（1）是TVD格式（1）可驗(yàn)證：Roe格式是TVD格式保證“系數(shù)非負(fù)”含義：“單調(diào)格式必是TVD格式”CopyrightbyLiXinliang25例7.2.1：考慮線性單波方程：試討論如下Lax-Wendroff格式二階中心人工粘性是否滿足Harten條件單調(diào)格式——只有一階精度對比條件：不滿足Harten條件CopyrightbyLiXinliang26知識回顧：Lax-Wendroff格式Taylor展開，寫出修正方程時(shí)-空二階精度巧妙添加人工粘性，不但克服了不穩(wěn)定性，而且抵消了時(shí)間誤差，提高了時(shí)間精度類似方法：Beam-Warming格式人工粘性二階精度迎風(fēng)差分人工粘性，且提高時(shí)間精度特點(diǎn)：全離散、時(shí)刻耦合CopyrightbyLiXinliang274.構(gòu)建TVD格式思路：對現(xiàn)有格式進(jìn)行改造，使之符合Harten條件通常在Roe、L-W、B-M（或其組合）基礎(chǔ)上改進(jìn)80年代初、這些格式是主流（1）以L-W格式為基礎(chǔ)改造的格式L-W原格式（2階）=1階迎風(fēng)+修正項(xiàng)

新格式=1階迎風(fēng)+限制函數(shù)*修正項(xiàng)引入限制函數(shù)（限制器）1階迎風(fēng)部分修正項(xiàng)CopyrightbyLiXinliang28顯然格式為LW（2階）可驗(yàn)證：格式為B-M（2階）

新格式=1階迎風(fēng)+j*（LW格式-1階迎風(fēng)）新格式：LW,BM均為線性格式，二者組合仍為二階根據(jù)Harten定理，可知時(shí)，可滿足TVD性質(zhì)(2)精度條件Beam-Warming二階精度區(qū)TVD區(qū)二階精度TVD區(qū)（二者交集）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}五個(gè)基架點(diǎn)被分成三個(gè)組1）若高精度逼近，必然利用多個(gè)基架點(diǎn)2）如果該基架點(diǎn)內(nèi)函數(shù)有間斷，會導(dǎo)致振蕩3）間斷不可能處處存在4）把基架點(diǎn)分成多個(gè)組（模板），

每個(gè)模板獨(dú)立計(jì)算j點(diǎn)導(dǎo)數(shù)的逼近?！玫蕉鄠€(gè)差分

5）根據(jù)每個(gè)模板的光滑程度，設(shè)定權(quán)重6）對多個(gè)差分結(jié)果進(jìn)行加權(quán)平均。光滑度越高，權(quán)重越大。如果某模板存在間斷，則權(quán)重趨于0；如果都光滑，則組合成更高階格式。CopyrightbyLiXinliang302.WENO格式的原理描述考慮線性單波方程：注：為了簡便，以非守恒型形式為例講授其思路，實(shí)際使用時(shí)，請采用下一節(jié)介紹的守恒形式（1）確定網(wǎng)格基架點(diǎn)：6個(gè)點(diǎn){j-3,j-2,j-1,j,j+1,j+2}

構(gòu)造出該基架點(diǎn)上的目標(biāo)差分格式計(jì)算這6個(gè)點(diǎn)可構(gòu)造5階迎風(fēng)差分：該格式為WENO的“目標(biāo)”格式，即，光滑區(qū)WENO逼近于該格式。利用Taylor展開，可唯一確定系數(shù)（可利用小程序coeff-schemes.f)實(shí)際上，還可利用分辨率優(yōu)化技術(shù)，可構(gòu)造出新的目標(biāo)格式（降低精度、提高分辨率，見第4講）。目前大量WENO的優(yōu)化版做這種工作。CopyrightbyLiXinliang31將這6個(gè)基架點(diǎn)分割成3個(gè)組（稱為模板）

每個(gè)組獨(dú)立計(jì)算的差分逼近

模板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}利用這三個(gè)模板的基架點(diǎn)，可構(gòu)造出逼近的3階精度差分格式計(jì)算j點(diǎn)的導(dǎo)數(shù)u’,竟然算出了三個(gè)不同的值，怎么辦？ENO方法：選擇最優(yōu)（最光滑）的，舍棄其余兩個(gè)WENO的處理方法：三個(gè)都要，加權(quán)平均它們。利用Taylor展開式，可唯一確定這些系數(shù)）（可利用小程序coeff-schemes.f)也可運(yùn)用優(yōu)化技術(shù)，降低精度、提高分辨率……CopyrightbyLiXinliang32（3）對這3個(gè)差分值進(jìn)行加權(quán)平均，得到總的差分值原則：A.模板內(nèi)函數(shù)越光滑，則權(quán)重越大；模板內(nèi)有間斷時(shí)，權(quán)重趨于0B.三個(gè)模擬內(nèi)函數(shù)都光滑時(shí)，這三個(gè)三階精度的逼近式可組合成一個(gè)五階精度的逼近式。“理想權(quán)重”（3.1）確定理想權(quán)重令：5階精度容易解出：CopyrightbyLiXinliang33（3.2）度量每個(gè)模板內(nèi)函數(shù)的光滑程度

IS越大，表示越不光滑。

光滑區(qū)，不同模板上的IS趨近同一值。具體形式見下一節(jié)。

（3.3）給出實(shí)際權(quán)重構(gòu)造IS方法很多，例如：

：第k個(gè)模板光滑區(qū)逼近O（1）量級間斷區(qū)量級，很大特點(diǎn)：間斷區(qū)權(quán)重很小光滑區(qū)，趨近于理想權(quán)重（3.4）給出最終的差分逼近CopyrightbyLiXinliang343.Jiang&Shu的五階WENO格式守恒型；目前使用的WENO格式均為守恒型針對方程：模板1模板2

模板3構(gòu)造差分格式如下：構(gòu)造方法與前文相同（但注意這里構(gòu)造的是通量，而前文是直接構(gòu)造差分格式）針對整個(gè)網(wǎng)格基，構(gòu)造出5階精度的通量（理想情況下的通量）并構(gòu)造出每個(gè)模板上的通量，計(jì)算出理想權(quán)重。仍利用程序coeff-schemes.f求系數(shù)理想權(quán)重光滑度量因子實(shí)際權(quán)重CopyrightbyLiXinliang35光滑度量因子的計(jì)算（Jiang&Shu）k=1k=2k=3其中：j-2j-1jj+1j+2是使用模板k得到的插值函數(shù)

利用{j-2,j-1,j}點(diǎn)上的值構(gòu)造的插值函數(shù)

，特點(diǎn)：光滑區(qū)趨近同一個(gè)值

非光滑區(qū)值遠(yuǎn)大于光滑區(qū)O(1)j點(diǎn)一階、二階導(dǎo)數(shù)的差分逼近（用模板k計(jì)算）代入CopyrightbyLiXinliang36最終5階WENO格式為正通量情況（a>0）

負(fù)通量情況（a<0）注：正通量差分格式中下標(biāo)“j+k”改成“j-k”即得到負(fù)通量的差分格式(除了第1式不變）注意:是j-1/2而不是j+1/2k=1k=2k=3CopyrightbyLiXinliang374.WENO格式的邊界處理（1）簡易的（降階）處理方法：如果某模板用到邊界外的點(diǎn)，簡單將該模板權(quán)重設(shè)為0即可如果用到邊界點(diǎn)外的點(diǎn)，則該權(quán)重設(shè)為0效果不錯(cuò)，但會邊界降階（推薦）（2）構(gòu)造單邊差分的WENO格式

優(yōu)點(diǎn)：精度高

缺點(diǎn)：穩(wěn)定性不易保證可能會出現(xiàn)負(fù)權(quán)重，造成不穩(wěn)定，可用如下文獻(xiàn)的方法處理：ShiJ,HuCQ,andShuCW,2002,ATechniqueofTreatingNegativeWeightsinWENOSchemes,JournalofComputationalPhysics175,108–127j=12345方法2：特征投影分裂（詳細(xì)步驟見本PPT第5-6頁）計(jì)算

利用上頁的公式計(jì)算正通量及負(fù)通量的WENO通量及

計(jì)算

計(jì)算CopyrightbyLiXinliang38推廣到Euler（或N-S）方程運(yùn)用分裂技術(shù)，可將上述方法推廣到Euler方程方法1：逐點(diǎn)分裂（又稱流通矢量分裂FVS）采用針對正通量（a>0）的方法計(jì)算采用針對負(fù)通量（a>0）的方法計(jì)算可采用Steger-Warming,L-F,VanLeer等分裂，見第4講1)計(jì)算，它們是的函數(shù)（推薦使用Roe平均計(jì)算）效果更好但計(jì)算量較大計(jì)算量小但效果略差有輕微振蕩數(shù)值測試：1維Sod激波管問題，網(wǎng)格點(diǎn)128

差分方法:7階精度經(jīng)典WENO(Jiang&Shu)分裂方式：Steger-WarmingFVS;特征投影分裂（算術(shù)平均/Roe平均）t=0.1時(shí)刻密度及速度分布Steger-WarmingFVS+WENO仍無法避免振蕩特征投影分裂+WENO可避免振蕩