計(jì)算流體力學(xué)講義

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

；力學(xué)所主樓219；82543801知識(shí)點(diǎn)：離散誤差的Fourier分析;

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

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

1講義、課件上傳至

(流體中文網(wǎng)）->“流體論壇”->“CFD基礎(chǔ)理論

”CopyrightbyLiXinliang§3.3差分格式的進(jìn)一步分析1.耗散與色散誤差2CopyrightbyLiXinliang精確解1階迎風(fēng)2階迎風(fēng)數(shù)值實(shí)驗(yàn)

時(shí)間推進(jìn)：3步TVD型Runge-Kutta,且時(shí)間步長(zhǎng)足夠小（誤差忽略）空間離散：1階及2階迎風(fēng)格式（20個(gè)網(wǎng)格點(diǎn)）實(shí)驗(yàn)觀察到的現(xiàn)象——兩類誤差：

振幅誤差

相位誤差（波速誤差）

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

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

其他格式……假設(shè)對(duì)于：有隱含假設(shè):線性差分格式

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

理想情況：

的誤差導(dǎo)致解的幅值誤差——耗散誤差

的誤差導(dǎo)致解傳播速度的誤差——色散誤差假設(shè)對(duì)于：有

反映了一個(gè)波內(nèi)的點(diǎn)數(shù)。PPW（波內(nèi)的點(diǎn)數(shù)）=CopyrightbyLiXinliang5耗散、色散誤差分別由修正波數(shù)的實(shí)部和虛部決定。關(guān)鍵參數(shù)：修正波數(shù)含義：反應(yīng)波數(shù)（譜）空間內(nèi)差分的誤差任意函數(shù)：定義：求導(dǎo)數(shù)，精確解差分解Fourier分析的任務(wù)計(jì)算出，并考差其與的逼近程度。考察格式分辨率（resolution）的重要指標(biāo)

精度：反映時(shí)的情況

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

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

方法1.理論計(jì)算

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

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

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

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

要求：精度不限；

網(wǎng)格基架點(diǎn)數(shù)不超過6個(gè)；

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

（即kr,ki曲線近可能接近準(zhǔn)確解）

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

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

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

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

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

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

快格式(FST)：

慢格式(SLW):

混合格式(MXD):

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

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

時(shí)間推進(jìn)：3階Runge-Kutta

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

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

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

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

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

4）

a<0時(shí)，同樣思路構(gòu)造

（利用對(duì)稱性，僅需把下標(biāo)j+k換成j-k即可）

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

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

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

目的：確定波傳播方向，便于使用迎風(fēng)差分方法：1）逐點(diǎn)分裂2）嚴(yán)格特征分裂

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

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

計(jì)算特征值

分裂特征值，計(jì)算

帶入左式，計(jì)算正、負(fù)流通矢量

計(jì)算計(jì)算設(shè)置，并注意對(duì)于曲線坐標(biāo)系僅需令三維問題同樣處理

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

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

優(yōu)點(diǎn)：無需矩陣運(yùn)算，計(jì)算量小

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

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

缺點(diǎn)：由于未解耦，各變量的誤差會(huì)相互傳遞27CopyrightbyLiXinliang概念澄清：

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

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

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

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

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

均可本人