第4講-差分方法2

計算流體力學講義

第四講有限差分法（2）李新亮lixl@

；力學所主樓219；82543801知識點：離散誤差的Fourier分析;

間斷周圍數(shù)值振蕩的原因；GVC格式；

模型方程向N-S方程的推廣；

1講義、課件上傳至

(流體中文網）->“流體論壇”->“CFD基礎理論

”CopyrightbyLiXinliang§3.3差分格式的進一步分析1.耗散與色散誤差2CopyrightbyLiXinliang精確解1階迎風2階迎風數(shù)值實驗

時間推進：3步TVD型Runge-Kutta,且時間步長足夠小（誤差忽略）空間離散：1階及2階迎風格式（20個網格點）實驗觀察到的現(xiàn)象——兩類誤差：

振幅誤差

相位誤差（波速誤差）

CopyrightbyLiXinliang3對以上“實驗現(xiàn)象”進行理論分析半離散分析：假設時間推進是精確的，僅分析空間離散帶來的誤差（難度小、常用）全離散分析：同時分析時、空離散的誤差（難度大）考查問題：實際上就是普通三角函數(shù)，采用復數(shù)形式僅僅是為了理論推導方便。用實數(shù)形式sin(kx),cos(kx)推導形式上略顯繁瑣。

精確解：差分格式：（1）

其他格式……假設對于：有隱含假設:線性差分格式

非線性系統(tǒng)作用于單波，會產生多個諧波（2）差分沒有誤差CopyrightbyLiXinliang4令：（1）式化為：“半離散化”：空間導數(shù)差分計算，時間方程（常微）精確計算如果，無誤差分析（修正波數(shù)）與誤差的關系

理想情況：

的誤差導致解的幅值誤差——耗散誤差

的誤差導致解傳播速度的誤差——色散誤差假設對于：有

反映了一個波內的點數(shù)。PPW（波內的點數(shù)）=CopyrightbyLiXinliang5耗散、色散誤差分別由修正波數(shù)的實部和虛部決定。關鍵參數(shù)：修正波數(shù)含義：反應波數(shù)（譜）空間內差分的誤差任意函數(shù)：定義：求導數(shù)，精確解差分解Fourier分析的任務計算出，并考差其與的逼近程度?？疾旄袷椒直媛剩╮esolution）的重要指標

精度：反映時的情況

分辨率：網格點數(shù)很少（例如波里面只有6個點）時的性能對于多尺度問題，分辨率更重要。犧牲精度，提高分辨率

優(yōu)秀的差分格式，1個波長里面6個點即可精度分辨率CopyrightbyLiXinliang6如何計算修正波數(shù)？定義：

方法1.理論計算

根據(jù)差分具體表達式及定義計算例1：令則：于是：1階迎風例2：2階迎風CopyrightbyLiXinliang7方法2：數(shù)值計算定義：Step1）選取計算域[0,2p],計算網格（例如64，128）Step2）給定波數(shù)k,生成函數(shù)值Step3)調用差分子程序，得到導數(shù)值Step4)通過Fourier反變換，得到譜：假設已有求差分的子程序（黑箱，已知是線性的）線性黑箱強調：研究CFD本身，不能只使用理論手段，還要用數(shù)值手段根據(jù)修正波數(shù)的定義，有Step5)改變k的值，重復2-5，得到對于的依賴關系。畫圖非線性情況會產生高次諧波，造成step4中隱含的假設無法成立

將Fourier分析手段拓展到非線性系統(tǒng)需要研究的課題隱含條件：只有波數(shù)為k的那個譜不為0（線性系統(tǒng)）CopyrightbyLiXinliang8中心差分格式的色散特性0：精確解；1：4階普通2：6階普通；3：4階緊致4：6階緊致；5：6階超緊致迎風差分格式的色散特性0：精確解，1：2階迎風2：5階迎風偏心3：3階迎風緊致4：5階迎風緊致每個波長里面2個網格點，譜方法的分辨率，差分法分辨率的極限（只有無窮階精度才能達到）20階超緊致格式——接近譜方法CopyrightbyLiXinliang9不同差分格式的色散誤差曲線結論：要求分辨率相同的情況下，

采用高階格式可放寬空間網格步長，從而減少計算量重要方向：高分辨率差分格式0：精確解1：2階迎風2:3階迎風3：3階迎風緊致4：5階迎風緊致指定誤差要求的情況下，不同差分格式能模擬的最大a(a越大，所需網格越少）

作業(yè)題1：構造高分辨率差分格式，并進行理論分析及數(shù)值實驗針對單波方程：對于空間導數(shù)，構造出一種不超過6點格式；并進行Fourier誤差分析，畫出kr,ki的曲線。

要求：精度不限；

網格基架點數(shù)不超過6個；

能夠分辨的波數(shù)范圍盡量寬；

（即kr,ki曲線近可能接近準確解）

給出差分的具體表達式，畫出kr,ki的曲線；

說明構造格式的階數(shù)，并采用本PPT第5頁的方法給出的精度驗證；

形如：……另外，進行如下數(shù)值驗證：空間采用20個網格點，采用新構造的差分格式離散；時間推進采用3步Runge-Kutta方法，時間步長可足夠小（例如0.01）。給出t=20,50兩個時刻的數(shù)值解，與精確解比較（畫圖），并給出數(shù)值解的L2模誤差。10CopyrightbyLiXinliang提示：1.如不使用優(yōu)化技術，則格式構造方法簡單，Taylor展開后解代數(shù)方程組即可。2.建議嘗試使用優(yōu)化技術

例：假設格式形式如下如果要求其有5階精度，則通過Taylor展開可得到6個方程，6個系數(shù)可直接解出。我們要求其有4階精度（當然3階，2階也可），于是Taylor展開只能提供5個方程。6個未知數(shù)（a1-a6），5個方程；有1個自由參數(shù)。調整這個自由參數(shù)，使得kr,ki曲線最為理想。

如何調整？1)可以人工調整，觀察kr,ki曲線，選取滿意的。2）可自動調整，設立一個優(yōu)化目標函數(shù)。例如

調整自由參數(shù)，使得該目標函數(shù)取最大值。思路：犧牲精度，提高分辨率11CopyrightbyLiXinliang附錄：部分差分格式…j-2j-1jj+1…表中的迎風差分格式均針對a>0當a<0時，需把下標的“j+k”換成“j-k”（例如把j+2換成j-2,把j-1換成j+1);并在表達式前加上“-”號。例：迎風偏斜格式：上游的基架點更多些（或上游權重更大）12CopyrightbyLiXinliang§3.4數(shù)值解的群速度及間斷處數(shù)值振蕩來源對于：有修正波數(shù)數(shù)值解色散誤差：數(shù)值解傳播的速度與精確解不一致數(shù)值解傳播偏快數(shù)值解傳播偏慢0：精確解；1：2階迎風；2：5階迎風偏心3：3階迎風緊致；4：5階迎風緊致

快格式(FST)：

慢格式(SLW):

混合格式(MXD):

特點：波數(shù)越高，誤差越嚴重1.色散誤差與群速度13CopyrightbyLiXinliangt=0.5時刻的精確解及數(shù)值解空間離散：五種不同格式；時間推進：3階Runge-Kutta【數(shù)值實驗】波的傳播問題觀察現(xiàn)象：1）高波數(shù)成分誤差嚴重，低波數(shù)成分誤差不明顯；2）二階Pade格式的解傳播速度快于精確解，其余格式偏慢；3）迎風型格式有耗散，尤其是二階迎風格式；

概念：群速度——波包傳播的速度14CopyrightbyLiXinliang2.間斷附近數(shù)值振蕩的來源【數(shù)值實驗】間斷的傳播計算域[0,1];計算網格點100

時間推進：3階Runge-Kutta

空間離散：1）二階中心差分

2階迎風及2階中心格式的色散特性2）二階迎風差分15CopyrightbyLiXinliang過激波數(shù)值振蕩的根源——色散誤差導致群速度不一致快格式慢格式波前振蕩波后振蕩=+++…群速度控制的基本思路（群速度控制GVC:Fu&Ma）：

間斷前、后分別采用快格式和慢格式，可有效抑制振蕩Zhuang&Zhang:抑制波動原則

示意圖：間斷的Fourier分解好思路16CopyrightbyLiXinliang利用GVC的思想構造可計算間斷的差分方法1）間斷的前后判據(jù)

簡易方法：則j點在間斷左側j-1jj+1則j點在間斷右側原理：越靠近間斷，振蕩越劇烈（a>0時，右側為“前”）2）根據(jù)GVC的思想構造格式間斷前：快格式；間斷后：慢格式；格式GVC23）改寫成為守恒型非線性情況，通常守恒型效果更好NND格式GVC2a17CopyrightbyLiXinliang

4）

a<0時，同樣思路構造

（利用對稱性，僅需把下標j+k換成j-k即可）

采用GVC2a(NND2a)格式的計算結果——消除振蕩使用NND2a（守恒形式）；NND2（普通形式）及1階迎風格式的計算結果將j換成j+118CopyrightbyLiXinliang作業(yè)題2：構造更高分辨率的GVC格式對于空間導數(shù)，構造出一種不超過6點的GVC格式。要求：a.精度不限；b.網格基架點數(shù)不超過6個；c.求解模型方程

計算結果間斷盡量保持“銳利”;計算結果振蕩盡量小。振蕩的定量判據(jù)：總變差（TotalVariation）：

間斷“銳利”的定量判據(jù)：間斷區(qū)內的點數(shù)？（自行設計）給出差分格式的表達式、色散/耗散分析（ki,kr曲線）；給出模型方程t=0.2的結果（空間100個網格點，計算域[0,1]，時間推進可采用3階Runge-Kutta方法）；與精確解及NND2a進行比較（畫在同一張圖上）建議：利用優(yōu)化方法19CopyrightbyLiXinliang§3.5從模型方程推廣到N-S方程（Euler方程）格式F+格式F-（教科書第6章）對流項：信息（波）從上游傳至下游——上游信息更重要——迎風差分擴散項：信息從中心向周圍擴散——不區(qū)分上、下游——中心差分迎風差分優(yōu)點：有效利用信息傳播的方向，增強穩(wěn)定性微分與差分方程的影響域N-S方程：單波方程：單波方程——一個波，容易判斷波傳播方向N-S對流項（Euler）——方程組：多波問題，復雜雙曲方程組的原則——特征分解，找到獨立傳播的波常系數(shù)矩陣A的情況——完全解耦，獨立求解變系數(shù)矩陣A的情況——局部討論20CopyrightbyLiXinliang1.Jacobian系數(shù)矩陣及其性質重要性質特點：A可以像常數(shù)一樣，和求導運算交換21CopyrightbyLiXinliang2.對流項的分裂

目的：確定波傳播方向，便于使用迎風差分方法：1）逐點分裂2）嚴格特征分裂

1）逐點分裂利用性質=+優(yōu)點：耗散小缺點：導數(shù)間斷方式A:特點：不必進行矩陣運算，計算量小Steger-Warming分裂A:Steger-Warming分裂22CopyrightbyLiXinliangSteger-Warming具體步驟（以一維為例）已知1）計算2）計算3）計算4）帶入（1）式得到5）利用不同的迎風格式，分別計算

（1）（后差，前差）6）計算7）時間推進23CopyrightbyLiXinliang二維問題的steger-Warming分裂令：則：具體使用步驟，以計算為例令

計算特征值

分裂特征值，計算

帶入左式，計算正、負流通矢量

計算計算設置，并注意對于曲線坐標系僅需令三維問題同樣處理

二維、三維具體公式見傅德薰等《計算空氣動力學》4.7節(jié)（158-162）書中公式有一定的排版錯誤，使用前務必重新仔細推導！24CopyrightbyLiXinliangB:L-F分裂特點：正特征值負特征值=+缺點：耗散偏大局部L-F分裂，每個點上計算全局L-F分裂，全局（一維）上計算足夠大數(shù)學性質（光滑性）最好，但耗散偏大常數(shù)與迎風格式結合，等價于人工粘性例如，可取25CopyrightbyLiXinliang方式很多=+S-W:L-F:=+VanLeer:=+26CopyrightbyLiXinliang

分裂后失去了A的性質（可以像常數(shù)一樣與求導交換）逐點分裂：

優(yōu)點：無需矩陣運算，計算量小

缺點：分裂后改變了特征方向，耗散大利用了性質一般情況下：變系數(shù)，

不能與導數(shù)交換實質：沒有做到解耦；只是把原變量重新組合，組合后波的傳播方向的保證f+向正向傳播，f-向負向傳播

缺點：由于未解耦，各變量的誤差會相互傳遞27CopyrightbyLiXinliang概念澄清：

流通矢量分裂本身不帶來耗散，但其會影響到差分的耗散;舉例：分裂過程耗散如果差分格式無耗散（例如都用中心差分），則通量分裂不帶來耗散。=+向上平移向下平移分裂差分格式耗散分裂后的流場越偏離原先流場，則總體耗散越大精確滿足，不引入誤差！如使用低精差分度格式，則對分裂形式敏感（推薦使用特征分裂）如使用高精度格式（低耗散），則對分裂形式不敏感（可使用逐點分裂）28CopyrightbyLiXinliang2.嚴格特征分裂——基架點上凍結系數(shù)常系數(shù)方程組：完全解耦變系數(shù)情況——局部凍結系數(shù)…j-2j-1jj+1…在基架點上系數(shù)不變計算：在差分基架點上Aj

不變，可按常矩陣處理局部凍結系數(shù)分別采用后差和前差優(yōu)點：嚴格保證（局部）特征方向，數(shù)值解質量好；缺點：大量矩陣運算，計算量大。29CopyrightbyLiXinliang通常寫成守恒型差分，計算…j-2j-1jj+1…在基架點上系數(shù)不變具體步驟：

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

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

可使用j,j+1點值的算術平均（如）或Roe平均（教材6.4節(jié)）；由計算；方法很多，例如前面介紹的或30CopyrightbyLiXinliang

均可本人