中科院計算流體力學(xué)最新講義第講差分方法

中科院計算流體力學(xué)最新講義第講差分方法第一頁，共三十八頁，2022年，8月28日1知識回顧1.差分格式的分辨率有效網(wǎng)格點數(shù)：一個波長里面的網(wǎng)格點數(shù)（PPW:PointperWavelength)修正波數(shù)2.群速度控制(GVC)第二頁，共三十八頁，2022年，8月28日3.流通矢量分裂“在二階迎風(fēng)與二階中心格式中選一個”

“選接近一階迎風(fēng)的”GVC2格式第三頁，共三十八頁，2022年，8月28日3§5.1非物理振蕩及TVD格式1.數(shù)值解中的非物理振蕩

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

數(shù)值方法錯誤計算了物理粘性——不足以克服振蕩物理問題本身也可能振蕩。但如果錯誤計算物理粘性，則會錯誤地加?。ɑ蛩p）振蕩。1）非物理振蕩的原因分析理論3：格式不能保單調(diào)第四頁，共三十八頁，2022年，8月28日4數(shù)值實驗

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

T=0.1時刻的u分布Re=200Dx=0.005現(xiàn)象：Dx一定時，減小Reynolds數(shù)可抑制振蕩Reynolds數(shù)一定時，減小Dx可抑制振蕩暗示是某一特征量Re=2000Dx=0.005Re=2000Dx=0.0005相同第五頁，共三十八頁，2022年，8月28日5對流-擴(kuò)散方程的特性：nn+1差分方程：某點的值是上一時刻周圍幾個點上值的線性組合物理上要求系數(shù)ak

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

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

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

單純靠物理粘性抑制振蕩，網(wǎng)格間距必須足夠小，通常難以實現(xiàn)網(wǎng)格足夠?。翰粫l(fā)生振蕩；網(wǎng)格小于激波的實際厚度，則不會振蕩網(wǎng)格Reynolds數(shù)足夠小時，物理粘性發(fā)揮作用，抑制振蕩第七頁，共三十八頁，2022年，8月28日73）

人工粘性

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

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

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

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

新格式=1階迎風(fēng)+限制函數(shù)*修正項1.以二階中心及二階迎風(fēng)格式為基礎(chǔ)的改造2階迎風(fēng)2階中心2個候選格式：思路1：

兩個里面選一個(GVC2,NND2)

思路2：

利用二者的組合改造、新格式二階迎風(fēng)二階中心第十四頁，共三十八頁，2022年，8月28日14顯然格式為2階中心可驗證：格式為2階迎風(fēng)新格式：二者組合仍為二階根據(jù)Harten定理，可知時，可滿足TVD性質(zhì)(2)精度條件二階精度區(qū)TVD區(qū)二階精度TVD區(qū)（二者交集）第十五頁，共三十八頁，2022年，8月28日15限制器（limiter)二階中心二階迎風(fēng)GVC2格式NND2格式GVC2與NND2格式的限制器第十六頁，共三十八頁，2022年，8月28日16常用的限制器Minmod型；Superbee型VanLeer型其他：如VanAlbada第十七頁，共三十八頁，2022年，8月28日17概念：保單調(diào)區(qū)1階迎風(fēng)的預(yù)測值2階中心的修正量2階迎風(fēng)的修正量精度差，但魯棒性好精度高，但有些情況下預(yù)測結(jié)果“不靠譜”作為“標(biāo)桿”檢驗高階修正量是否可用趨勢相反時，不可用；相差超過2倍時，不可用第十八頁，共三十八頁，2022年，8月28日18歷史上，TVD格式是在Roe、L-W、B-M（或其組合）基礎(chǔ)上改進(jìn)80年代初、這些格式是主流2以L-W格式為基礎(chǔ)改造的格式L-W引入限制函數(shù)（限制器）1階迎風(fēng)部分修正項第十九頁，共三十八頁，2022年，8月28日19顯然格式為LW（2階）可驗證：格式為B-M（2階）

新格式=1階迎風(fēng)+j*（LW格式-1階迎風(fēng)）新格式：LW,BM均為線性格式，二者組合仍為二階根據(jù)Harten定理，可知時，可滿足TVD性質(zhì)(2)精度條件Beam-Warming二階精度區(qū)TVD區(qū)二階精度TVD區(qū)（二者交集）第二十頁，共三十八頁，2022年，8月28日20§5.2其他激波捕捉格式簡介

1.Godnov格式思路：利用精確Riemann解難點：精確Riemann解要求間斷兩側(cè)物理量為常數(shù)對策：采用分片常數(shù)代替原先的函數(shù)方法：1）

采用分片常數(shù)代替原先的函數(shù)

即假設(shè)[xj-1/2,xj+1/2]區(qū)間物理量為xj點上的值2）在每個間斷處，精確求解Riemann問題，得到Dt時刻物理量的分布。（假設(shè)Dt足夠小，各間斷之間無相互影響，因此可獨立求解）。3）在區(qū)域[xj-1/2,xj+1/2]上平均，得到Dt時刻xj點上物理量的值。4）重復(fù)1-3，直到指定時刻。j-1jj+1

間斷1間斷2間斷3……第二十一頁，共三十八頁，2022年，8月28日212.NND格式“在二階迎風(fēng)與二階中心格式中選一個”

“選接近一階迎風(fēng)的”如果二階中心與二階迎風(fēng)給出的修正趨勢相反，干脆不修正。Minmod(a,b):a,b符號相同時，取絕對值小的；符號相反時，取0.第二十二頁，共三十八頁，2022年，8月28日k=1/3三階迎風(fēng)

k=0k=-12.MUSCL格式（VanLeer）二階迎風(fēng)Formm格式k=1二階中心1階迎風(fēng)修正部分間斷處都會產(chǎn)生振蕩思路：

振蕩由修正項引起，需要對修正項進(jìn)行限制第二十三頁，共三十八頁，2022年，8月28日23新格式函數(shù)minmod(a,b)

a,b符號相反，則返回0符號相同，則返回絕對值小的1階迎風(fēng)修正部分j-1jj+1j-2j+1/2思路：

修正部分是由j-1,j,j+1三點計算而得，是與的插值舊格式二者的組合：j-1,j計算；j,j+1計算

如果兩個修正方案趨勢（符號）相反，干脆不修正

如果趨勢相同，則?。ń^對值）最小者——避免振蕩兩個修正方案選前者：二階迎風(fēng)

選后者：二階中心

組合之：Formm、3階迎風(fēng)…第二十四頁，共三十八頁，2022年，8月28日24含VanAlbada限制器的3階MUSCL格式推薦方法：VanAlbada限制器光滑區(qū)相差不大間斷區(qū)，

相差很大光滑區(qū)相差不大，s趨近于1間斷區(qū)相差很大，s趨近于0光滑區(qū)，趨近3階迎風(fēng)間斷區(qū)，趨近1階迎風(fēng)理論上并不滿足TVD性質(zhì)；但實際效果很好a<0的情況，把j+k改成j-k就可以了jj第二十五頁，共三十八頁，2022年，8月28日25關(guān)于MUSCL格式的一些注記1.MUSCL格式是一系列格式?？蓮亩A迎風(fēng)、二階中心、Formm,三階迎風(fēng)等格式經(jīng)過限制器改造而得到。限制器有Minmod,VanAlbada等多種形式；MUSCL格式不一定都能滿足TVD特性；NND格式是一種特殊的MUSCL格式（k=1,b=1的情況）MUSCL格式在CFD史上有重要貢獻(xiàn)

MUSCL格式開創(chuàng)了利用基本變量重構(gòu)（而不是利用通量重構(gòu)）的新思路jj+1/2j-1j+1計算方法1：利用計算方法2：利用計算出

再利用計算出實際過程中，為了利用迎風(fēng)技術(shù)，還要進(jìn)行通量分裂第二十六頁，共三十八頁，2022年，8月28日26

§5.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}五個基架點被分成三個組1）若高精度逼近，必然利用多個基架點2）如果該基架點內(nèi)函數(shù)有間斷，會導(dǎo)致振蕩3）間斷不可能處處存在4）把基架點分成多個組（模板），

每個模板獨立計算j點導(dǎo)數(shù)的逼近。——得到多個差分

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

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

每個組獨立計算的差分逼近

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

IS越大，表示越不光滑。

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

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

：第k個模板光滑區(qū)逼近O（1）量級間斷區(qū)量級，很大特點：間斷區(qū)權(quán)重很小光滑區(qū)，趨近于理想權(quán)重（3.4）給出最終的差分逼近各階導(dǎo)數(shù)的差分表達(dá)式都可作為光滑度量因子使用第三十一頁，共三十八頁，2022年，8月28日313.Jiang&Shu的五階WENO格式守恒型；目前使用的WENO格式均為守恒型針對方程：模板1模板2

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

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

，特點：光滑區(qū)趨近同一個值

非光滑區(qū)值遠(yuǎn)大于光滑區(qū)O(1)j點一階、二階導(dǎo)數(shù)的差分逼近（用模板k計算）代入第三十三頁，共三十八頁，2022年，8月28日33光滑因子ISk的進(jìn)一步討論（光滑區(qū)精度分析）同樣：在光滑區(qū)，

趨近于同一個值，相互差距是一個6階小量

是小量，可忽略5階迎風(fēng)格式的通量格式具有5階精度第三十四頁，共三十八頁，2022年，8月28日34最終5階WENO格式為正通量情況（a>0）

負(fù)通量情況（a<0）注：正通量差分格式中下標(biāo)“j+k”改成“j-k”即得到負(fù)通量的差分格式(除了第1式不變）注意:是j-1/2而不是j+1/2第三十五頁，共三十八頁，2022年，8月28日35k=1k=2k=34.WENO格式的邊界處理（1）簡易的（降階）處理方法：如果某模板用到邊界外的點，簡單將該模板權(quán)重設(shè)為0即可如果用到邊界點外的點，則該權(quán)重設(shè)為0效果不錯，但會邊界降階（推薦