FVM in CFD 學習筆記_第12章_高分辨率格式

學習自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 12 High Resolution Schemeshtml

本章繼續講解對流項格式的發展,討論如何對高階(High Order HO)格式施加有界性來產生高分辨率(High Resolution HR)格式。HR格式是HO廓線和對流有界準則(Convection Boundness Criterion CBC)的結合,以保證解不會出現振盪。將介紹發展HR格式的歸一化變量公式(Normalized Variable Formulation(NVF))和總變差衰減(Total Variation Diminishing(TVD))框架,這倆方法看起來不一樣但實際上幾乎是同樣的。分別展現了在NVF和TVD中用於使HR格式更具形象化的歸一化變量圖(Normalized Variable Diagram (NVD))和Sweby(或 r ψ r-\psi )圖。針對NVF和TVD,具體給出了一些HR格式的函數關係。還有上一章講到的延遲修正(Deferred Correction(DC)),介紹了兩種技術來實如今結構網格和非結構網格的HO和HR格式,即,背風加權因子(Downwind Weighing Factor(DWF))方法和歸一化加權因子(Normalized Weighing Factor(NWF))方法。web

1 歸一化變量公式(Normalized Variable Formulation(NVF))

歸一化變量公式(NVF)是描述和分析高分辨率(HR)格式的一個框架,其由Leonard引入,隨着Gaskell和Lau的簡化對流界限準則(CBC)的引入而聞名,歸一化變量圖(NVD)是發展和分析HO和HR格式的有用工具。算法

NVF是基於局部歸一化變量的面公式,用於構造面 f f 處的變量 ϕ f \phi_f ,該方法用到的是迎風( ϕ C \phi_C )、背風( ϕ D \phi_D )、和遠迎風( ϕ U \phi_U )節點的值,以下圖所示。架構

在這裏插入圖片描述

歸一化變量的表達式爲
ϕ ~ = ϕ ϕ U ϕ D ϕ U \tilde \phi=\frac{\phi-\phi_U}{\phi_D-\phi_U}
根據這個歸一化關係, ϕ f \phi_f 的函數
ϕ f = f ( ϕ U , ϕ C , ϕ D ) \phi_f=f(\phi_U,\phi_C,\phi_D)
轉變爲
ϕ ~ f = f ( ϕ ~ C ) \tilde\phi_f=f(\tilde \phi_C)
之因此沒有了 ϕ U \phi_U ϕ D \phi_D ,是由於他倆的歸一化值變成了
ϕ ~ U = 0 ϕ ~ D = 1 \tilde \phi_U=0 \quad\quad \tilde \phi_D=1
那麼 ϕ C \phi_C 的歸一化值 ϕ ~ C \tilde\phi_C 變成了 ϕ \phi 場光滑性的衡量指標, ϕ ~ C \tilde\phi_C 的值在0到1之間( 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 )表明着單調廓線, ϕ ~ C \tilde\phi_C 的值小於0( ϕ ~ C < 0 \tilde\phi_C<0 )或者大於1( ϕ ~ C > 1 \tilde\phi_C>1 )則代表 C C 處存在極值。此外, ϕ ~ C 0 \tilde\phi_C\approx0 ϕ ~ C 1 \tilde\phi_C\approx1 代表是個梯度突跳(gradient jump)。這些構型都展現在下圖中,圖(a) ϕ ~ C < 0 \tilde\phi_C<0 ,極值;(b) ϕ ~ C > 1 \tilde\phi_C>1 ,極值;(c) ϕ ~ C = 1 \tilde\phi_C=1 ,突跳;(d) ϕ ~ C = 0 \tilde\phi_C=0 ,突跳;(e) 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 ,單調。
在這裏插入圖片描述app

在這裏插入圖片描述

歸一化也可用於把高階HO格式的函數關係轉化爲 ϕ ~ f \tilde\phi_f ϕ ~ C \tilde\phi_C 之間的線性關係,例如,前一章所講到的HO格式的函數關係可轉化成下面的歸一化形式框架

格式 原形式 歸一化形式
迎風Upwind ϕ f = ϕ C \phi_f=\phi_C ϕ ~ f = ϕ ~ C \Rightarrow\tilde\phi_f=\tilde\phi_C
中心差分CD ϕ f = 1 2 ( ϕ C + ϕ D ) \phi_f=\displaystyle\frac{1}{2}(\phi_C+\phi_D) ϕ ~ f = 1 2 ( 1 + ϕ ~ C ) \Rightarrow\tilde\phi_f=\displaystyle\frac{1}{2}(1+\tilde\phi_C)
二階迎風SOU ϕ f = 3 2 ϕ C 1 2 ϕ U \phi_f=\displaystyle\frac{3}{2}\phi_C-\frac{1}{2}\phi_U ϕ ~ f = 3 2 ϕ ~ C \Rightarrow\tilde\phi_f=\displaystyle\frac{3}{2}\tilde\phi_C
FROMM ϕ f = ϕ C + 1 4 ( ϕ D ϕ U ) \phi_f=\phi_C+\displaystyle\frac{1}{4}(\phi_D-\phi_U) ϕ ~ f = ϕ ~ C + 1 4 \Rightarrow\tilde\phi_f=\tilde\phi_C+\displaystyle\frac{1}{4}
QUICK ϕ f = 3 4 ϕ C 1 8 ϕ U + 3 8 ϕ D \phi_f=\displaystyle\frac{3}{4}\phi_C-\frac{1}{8}\phi_U+\frac{3}{8}\phi_D ϕ ~ f = 3 8 + 3 4 ϕ ~ C \Rightarrow\tilde\phi_f=\displaystyle\frac{3}{8}+\displaystyle\frac{3}{4}\tilde\phi_C
背風Downwind ϕ f = ϕ D \phi_f=\phi_D ϕ ~ f = 1 \Rightarrow\tilde\phi_f=1

如此一來,對於基於三個節點值的全部高階格式, ϕ ~ f \tilde\phi_f 總能表示成 ϕ ~ C \tilde\phi_C 的線性函數關係,即, ϕ ~ f = l ϕ ~ C + k \tilde\phi_f=l\tilde\phi_C+k ,其中 l l k k 是由格式所決定的係數。所以,若是把 ϕ ~ f \tilde\phi_f 做爲 ϕ ~ C \tilde\phi_C 的函數繪製在 ( ϕ ~ C , ϕ ~ f ) (\tilde\phi_C,\tilde\phi_f) 平面上,那麼這些格式的函數關係在圖中將呈現爲一條直線,所獲得圖形就是所謂的歸一化變量圖(Normalized Variable Diagram (NVD))。上述格式繪製出的NVD以下。ide

在這裏插入圖片描述

NVD代表,除了一階迎風格式和背風格式,全部的二階和三階格式都經過了點 Q ( 0.5 , 0.75 ) Q(0.5,0.75) (均勻網格狀況下)。實際上,若是一種格式具備二階精度,那麼它必需要經過 Q Q 點。此外,若是其在 Q Q 處的斜率是0.75,那麼這種格式是三階精度的(好比,QUICK格式)。迎風格式是很是擴散的,而背風格式是很是壓縮的(反擴散)。所以,從NVD中能夠形象地推出,若格式的函數關係與迎風格式的很是接近,那麼它是擴散的,若是與背風格式的很是接近,那麼它是壓縮的。svg

前一章所講的高階格式比一階迎風格式的截斷偏差要小不少,還能保持穩定性。然而,這些格式的一個主要缺點是它們的無界性,即,它們傾向於在突跳和尖銳梯度處產生輸運變量的過沖/下衝甚至是振盪。縱然在某些應用中,小的過沖或是振盪是尚可容忍的,然而,在某些情形下它們會致使災難性的後果,好比在湍流計算中(輸運變量多是湍流粘性,算不許就完蛋了)。函數

這些在尖銳梯度處的振盪特性是全部高階線性對流格式的特性。實際上這些格式並非單調的,由於它們會產生局部最大或最小值,即,它們不能保持極值(extrema preserving)。若是一個格式能保持極值,那麼解中的最大值必須不能再增大,且最小值不能再減少(該格式不能產生過沖/下衝)。實際上,Godunov和Ryabenki代表,任何單調的線性數值格式最多隻能是一階精度的,這也就意味着全部高階線性格式不能保持單調性,而若要構造保持單調的格式,就須要使用非線性限制函數(non-linear limiter functions)。有了這樣的理解,發展高階無振盪對流格式的法子就能被分紅兩大組,第一類是,給一階迎風格式添加上反擴散通量來加以限制(顯然,是限制一階迎風格式的假擴散特性),該添加項要使得產生的格式可以求解出尖銳的梯度而不會產生振盪纔好;第二類是,給無界的高階格式引入光滑擴散通量,來衰減非物理振盪。工具

因爲它們的多步驟天性和在平衡兩種通量時的困難,通量混合技術在數值算法中的代價是高昂的。這也就是爲何在本書中,兩種發展高分辨率算法的方法,只講到通量限制方法。第一種是複合流程,其把高階格式與有界低階格式聯合起來,且由一特定判據所控制的開關來切換兩種格式。第二種方法是基於把擴散的一階迎風項上添加一個反擴散通量,由通量限制器實現。此時,獲得的高分辨率格式也稱爲總變差衰減(Total Variational Diminishing(TVD))格式,隨後小節中將會詳細講述。

先講解複合格式方法,用歸一化變量公式(NVF)的架構來說解,並用歸一化變量圖(NVD)使講解更具形象。所以,將首先講解NVF和NVD,NVD的使用將做爲一重要工具,來引出爲確保高階插值算法知足有界性,需定義什麼準則。

2 對流有界準則(The Convection Boundness Criterion(CBC))

數值格式的預期效果應能保持其所描述或近似現象的物理特性,所以有界對流格式應該知足的條件,可經過分析對流的物理特性來更好理解。由於對流是把物理量從上游向下遊輸運的,那麼對於對流的近似應擁有該輸運特性。所以,數值對流格式應該是偏迎風的,否則的話它們缺少對流穩定性。那麼,除了橫跨界面兩側的節點 ϕ C \phi_C ϕ D \phi_D 以外,遠迎風節點的值 ϕ U \phi_U 對於分析對流格式也是很是重要的,而再遠的節點就意義較小了。在前面將的NVF中,值都作了歸一化處理,也考慮了 ϕ U \phi_U 的影響,其很是好使,可有助於肯定數值對流格式在什麼條件下是單調的。Spekreise和Barth和Jaspersen定義了一種單調格式(或者說是有界格式),該格式包含了面周圍的全部鄰居,Leonard和Gaskel和Lau基於他們的單調性定義,只給出了局部座標系統上的鄰近節點的定義形式,即
min ( ϕ C , ϕ D ) ϕ f max ( ϕ C , ϕ D ) \min(\phi_C,\phi_D)\le\phi_f\le\max(\phi_C,\phi_D)
歸一化後,上式條件變爲
min ( ϕ ~ C , 1 ) ϕ f max ( ϕ ~ C , 1 ) \min(\tilde\phi_C,1)\le\phi_f\le\max(\tilde\phi_C,1)
Gaskell和Lau提出的對於隱式穩態流動計算中的對流有界準則(The Convection Boundness Criterion(CBC))代表,一個擁有有界特性的格式,其函數關係應該是連續的,其界限應該是在 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 單調區間內位於下限 ϕ ~ C \tilde\phi_C 和上限 1 1 之間,應該經過點 ( 0 , 0 ) (0,0) ( 1 , 1 ) (1,1) ,應該在 ϕ ~ C > 1 \tilde\phi_C>1 ϕ ~ C < 0 \tilde\phi_C<0 區間的函數關係 f ( ϕ ~ C ) f(\tilde\phi_C) 等於 ϕ ~ C \tilde\phi_C ,上述條件在數學上的表達式爲
ϕ ~ f = { f ( ϕ ~ C ) c o n t i n u o u s f ( ϕ ~ C ) = 1 i f   ϕ ~ C = 1 f ( ϕ ~ C )   w i t h   ϕ ~ C < f ( ϕ ~ C ) < 1 i f   0 < ϕ ~ C < 1 f ( ϕ ~ C ) = 0 i f   ϕ ~ C = 0 f ( ϕ ~ C ) = ϕ ~ C i f   ϕ ~ C < 0   o r   ϕ ~ C > 1 \tilde\phi_f=\begin{cases} f(\tilde\phi_C) & continuous \\ f(\tilde\phi_C)=1 & if~\tilde\phi_C=1 \\ f(\tilde\phi_C)~with~\tilde\phi_C<f(\tilde\phi_C)<1 & if~0<\tilde\phi_C<1 \\ f(\tilde\phi_C)=0 & if~\tilde\phi_C=0 \\ f(\tilde\phi_C)=\tilde\phi_C & if~\tilde\phi_C<0~or~\tilde\phi_C>1 \end{cases}
換句話說,知足上述關係的格式都是有界的格式,上述關係在NVD中的區域即下圖的陰影部分。
在這裏插入圖片描述
在這裏插入圖片描述

結合上面倆圖,就很是容易理解對流有界準則了。當 ϕ C \phi_C 是單調廓線的時候,在單元面處的插值廓線不該產生新的極值,因此其應該由面兩側的節點 ϕ \phi 值所限制,即 min ( ϕ C , ϕ D ) ϕ f max ( ϕ C , ϕ D ) \min(\phi_C,\phi_D)\le\phi_f\le\max(\phi_C,\phi_D) 。當 ϕ C \phi_C 逼近 ϕ D \phi_D 且仍在單調區間的時候, ϕ f \phi_f 也應該趨近於 ϕ D \phi_D 。當 ϕ C \phi_C 等於 ϕ D \phi_D 的時候, ϕ f \phi_f 也應該等於 ϕ D \phi_D ,此時 ( ϕ ~ C , ϕ ~ f ) (\tilde\phi_C,\tilde\phi_f) 經過點 ( 1 , 1 ) (1,1) 。當 ϕ C \phi_C 的值使得 ϕ ~ C > 1 \tilde\phi_C>1 ϕ f \phi_f 應給成迎風值,即, ϕ C \phi_C 。這樣作的效果是,在知足了 ϕ f \phi_f 是由單元面兩側的節點值所限定的前提下,產生了最大可能的出流條件。如此一來,就會把錯誤的振盪衰減掉,由於 ϕ C \phi_C 將趨向於一個更低的值,得益於在這些條件中出流比入流要大的緣故(這句話的意思是說,對於單元 C C 來講,此時其流出的量比流入的量要多,因此 ϕ C \phi_C 會變小,而不會出現振盪)。

所以,若是沒有內在的物理機制(好比源項)來產生極值,那麼極值將會消亡。一樣的機理也對 ϕ ~ C < 0 \tilde\phi_C<0 的狀況成立。隨着在單調區域的 ϕ C \phi_C 逼近 ϕ U \phi_U ϕ f \phi_f 也會等於迎風值 ϕ C \phi_C 直到 ϕ C = ϕ U \phi_C=\phi_U ,代表 ϕ f \phi_f 的廓線要通過點 ( 0 , 0 ) (0,0) 。當 ϕ ~ C < 0 \tilde\phi_C<0 ϕ ~ C > 1 \tilde\phi_C>1 時,對流效應占優,因此這些區域選擇迎風近似就恰如其分了。

3 高分辨率格式

使用NVD來構造有界高階格式,即,高分辨率格式,是相對比較簡單易行的。任何高階格式可使用ad-hoc集合曲線(特殊曲線集,即分段直線所鏈接成的曲線,或者說是多點構造的曲線,兩點間用直線表示)來令其具有有界性。

爲了構造高分辨率格式,在區域 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 的單調廓線應該經過點 ( 0 , 0 ) (0,0) ( 1 , 1 ) (1,1) ,同時須要在NVD中的上三角陰影區域。另外一方面,在非單調區域,即, ϕ ~ C < 0 \tilde\phi_C<0 ϕ ~ C > 1 \tilde\phi_C>1 區域,廓線應該使用迎風廓線。用這種方式構造的一些很是衆所周知的高分辨率格式展現在下圖中(關於MUSCL格式的圖形和後面公式並不對應?鑑於書中所講,應以公式爲準)。

NVD of several HR schemes

爲了提升收斂特性,任何複合高分辨率格式應該避免在其廓線鏈接點處出現較硬的角度,也不該該出現水平或是垂直的廓線。例如,對於SMART格式,它是由QUICK格式構造的,爲了進一步提升其收斂性呢,能夠作些微小的修正。即,把垂直廓線修正成區間 0 < ϕ ~ C < 1 / 6 0<\tilde\phi_C<1/6 ϕ ~ f = 3 ϕ ~ C \tilde\phi_f=3\tilde\phi_C ,對於STOIC格式,則是用 0 < ϕ ~ C < 1 / 5 0<\tilde\phi_C<1/5 區間。同時把SMART和STOIC格式中的水平廓線也做以修正,一種修正方法是把區間 7 / 10 < ϕ ~ C < 1 7/10<\tilde\phi_C<1 所對應的 9 / 10 < ϕ ~ f < 1 9/10<\tilde\phi_f<1 替換成 ϕ ~ f = ϕ ~ C / 3 + 2 / 3 \tilde\phi_f=\tilde\phi_C/3+2/3 。亦可用其它修正方法(好比,能夠對複合廓線在 0.95 < ϕ ~ C < 1 0.95<\tilde\phi_C<1 的水平區間進行修正)。一樣,也能夠對有界CD格式作修正來提升其收斂特性。SMART、STOIC、SUPERBEE格式的修正NVD以下圖所示。

在這裏插入圖片描述

數學上,這些複合高分辨率格式的函數關係爲
MINMOD ϕ ~ f = { 3 2 ϕ ~ C 0 ϕ ~ C 1 2 1 2 ϕ ~ C + 1 2 1 2 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Bounded CD ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 ϕ ~ C e l s e w h e r e OSHER ϕ ~ f = { 3 2 ϕ ~ C 0 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e SMART ϕ ~ f = { 3 4 ϕ ~ C + 3 8 0 ϕ ~ C 5 6 1 5 6 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified SMART ϕ ~ f = { 3 ϕ ~ C 0 ϕ ~ C 1 6 3 4 ϕ ~ C + 3 8 1 6 ϕ ~ C 7 10 1 3 ϕ ~ C + 2 3 7 10 ϕ ~ C 1 ϕ ~ C e l s e w h e r e STOIC ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 2 3 4 ϕ ~ C + 3 8 1 2 ϕ ~ C 5 6 1 5 6 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified STOIC ϕ ~ f = { 3 ϕ ~ C 0 ϕ ~ C 1 5 1 2 ϕ ~ C + 1 2 1 5 ϕ ~ C 1 2 3 4 ϕ ~ C + 3 8 1 2 ϕ ~ C 7 10 1 3 ϕ ~ C + 2 3 7 10 ϕ ~ C 1 ϕ ~ C e l s e w h e r e MUSCL ϕ ~ f = { 2 ϕ ~ C 0 ϕ ~ C 1 4 ϕ ~ C + 1 4 1 4 ϕ ~ C 3 4 1 3 4 ϕ ~ C 1 ϕ ~ C e l s e w h e r e SUPERBEE ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 2 3 2 ϕ ~ C 1 2 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified SUPERBEE ϕ ~ f = { 2 ϕ ~ C 0 ϕ ~ C 1 3 1 2 ϕ ~ C + 1 2 1 3 ϕ ~ C 1 2 3 2 ϕ ~ C 1 2 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e \begin{aligned} \text{MINMOD} \quad &\tilde\phi_f=\begin{cases} \frac{3}{2}\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{1}{2}\tilde\phi_C+\frac{1}{2} & \frac{1}{2}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Bounded CD} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2}\tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{OSHER} \quad &\tilde\phi_f=\begin{cases} \frac{3}{2}\tilde\phi_C & 0\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{SMART} \quad &\tilde\phi_f=\begin{cases} \frac{3}{4}\tilde\phi_C+\frac{3}{8} & 0\le\tilde\phi_C\le\frac{5}{6} \\ 1 & \frac{5}{6}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified SMART} \quad &\tilde\phi_f=\begin{cases} 3\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{6} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{6}\le\tilde\phi_C\le\frac{7}{10} \\ \frac{1}{3}\tilde\phi_C+\frac{2}{3} & \frac{7}{10}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{STOIC} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2}\tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{2}\le\tilde\phi_C\le\frac{5}{6} \\ 1 & \frac{5}{6}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified STOIC} \quad &\tilde\phi_f=\begin{cases} 3\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{5} \\ \frac{1}{2}\tilde\phi_C+\frac{1}{2} & \frac{1}{5}\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{2}\le\tilde\phi_C\le\frac{7}{10} \\ \frac{1}{3}\tilde\phi_C+\frac{2}{3} & \frac{7}{10}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{MUSCL} \quad &\tilde\phi_f=\begin{cases} 2\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{4} \\ \tilde\phi_C+\frac{1}{4} & \frac{1}{4}\le\tilde\phi_C\le\frac{3}{4} \\ 1 & \frac{3}{4}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{SUPERBEE} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2} \tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{2} \tilde\phi_C & \frac{1}{2}\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified SUPERBEE} \quad &\tilde\phi_f=\begin{cases} 2 \tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{3} \\ \frac{1}{2} \tilde\phi_C+\frac{1}{2} & \frac{1}{3} \le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{2} \tilde\phi_C & \frac{1}{2}\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \end{aligned}
根據這種方法還發展了許多其它的格式,好比CLAM、UTOPIA、SHARP、ULTRA-SHARP等,不逐一列舉。

4 TVD框架(總變差衰減框架)

另外一種發展高分辨率對流格式的流行方法是總變差衰減(Total Variation Diminishing(TVD))框架。在數值求解變量 ϕ \phi 的對流偏微分方程時,總變差(TV)定義爲
T V = i ϕ i + 1 ϕ i TV=\sum_i|\phi_{i+1}-\phi_i|
其中 i i 表明空間求解域內的節點標識。若是解中的TV不隨時間變化,那麼該數值方法就被稱爲是總變差衰減的,數學上,等效於
T V ( ϕ t + Δ t ) T V ( ϕ t ) TV(\phi^{t+\Delta t}) \le TV(\phi^t)
在其影響深遠的文章裏,Harten證實了單調格式是TVD的,並且一個TVD格式是保持單調性的。一個保持單調性的格式並不會在求解域內產生任何新的局部極值,即,(原有的)局部最小值不會再減少,(原有的)局部最大值不會再增大。

這兒並不打算給出TVD方法的完整數學推導,而是,簡單地講解下構造TVD格式的方法,這裏的方法是基於Sweby的工做。

回想下前面講過的對流格式,一階迎風格式具備強擴散性,而二階中心差分格式具備高色散性。要是能有一種格式介於迎風和中心差分之間,即,擁有迎風格式的穩定性和中心差分格式的精度,豈不是很好。爲構造這樣的格式,從中心差分格式中尋找靈感
ϕ f = 1 2 ( ϕ D + ϕ C ) C D = ϕ C u p w i n d + 1 2 ( ϕ D ϕ C ) a n t i d i f f u s i v e   f l u x \phi_f=\underbrace{\frac{1}{2}(\phi_D+\phi_C)}_{CD}=\underbrace{\phi_C}_{upwind}+\underbrace{\frac{1}{2}(\phi_D-\phi_C)}_{anti-diffusive~flux}
這裏的 C C D D 依舊錶明着迎風節點和背風節點,而 f f 則表示節點 C C D D 之間的單元面。不難發現,中心差分格式能夠寫成迎風格式和一個通量的組合形式,該通量應該是起到了反擴散效應,從而把本來迎風格式的擴散性整成了中心差分格式的色散性。該通量的效果是使格式有了二階精度,負面效應則是因爲數值擴散效應的減弱而出現了非物理振盪。所以,一個較好的處理方式應該是,把這個反擴散通量的一部分添加到迎風格式中去,從而既能保持二階精度,又不會產生非物理振盪。一種方法是把該通量乘上個限制函數(也叫限制器或通量限制器),在振盪可能發生的區域(如,較大梯度的地方)防止其過分使用(好比直接把它設成0),而在光滑區域發揮其最大貢獻。把這種限制器定義爲 ψ ( r ) \psi(r) ,其中 r r 常取爲兩個連續梯度的比值,這樣, ϕ f \phi_f 的計算爲
ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) w i t h   r f = ϕ C ϕ U ϕ D ϕ C \phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)\quad\quad with~r_f=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
其中 U U 是遠迎風點,而 D D 是背風點, C C 仍爲迎風點。爲了保證反擴散通量的符號特性, ψ ( r f ) \psi(r_f) 須要是非負的,即 ψ ( r f ) 0 \psi(r_f)\ge0 (否則起不到反擴散效果啊)。此外,若是 r f < 0 r_f<0 則代表節點 C C 是局部極值,而 r f = 0 r_f=0 代表有局部的陡峭梯度(有個臺階,突跳狀況),此時應該讓 ψ ( r f ) = 0 \psi(r_f)=0 ,這樣才能不讓反擴散通量起做用,防止出現非物理振盪。那麼,發展TVD格式就變成了找到限制器 ψ ( r f ) \psi(r_f) ,使得數值格式是TVD或者單調的。這些限制器必需要知足哪些條件,才能保證對流格式是保單調的呢?

感受書中這部分(P444-446)的理論推導過程並不正確,但結論毋庸置疑是沒錯的,因爲正確的推導我也沒搞清楚,因此直接把結論給出來吧。那就是這個限制器 ψ ( r ) \psi(r) 必需要知足 r > 0 r>0 時(爲方便計,把下標 f f 去掉了)
ψ ( r ) 2 a n d ψ ( r ) 2 r \psi(r)\le2\quad and \quad \psi(r)\le2r
與前面提到的 r 0 r\le0 ψ ( r ) = 0 \psi(r)=0 聯合起來,給出TVD格式應該知足的條件,即限制器應該知足的條件爲
ψ ( r ) = { m i n ( 2 r , 2 ) r > 0 0 r 0 \psi(r)=\begin{cases} min(2r,2) & r>0 \\ 0 & r\le0 \end{cases}
在這裏插入圖片描述

如上圖所示,這個條件能夠畫在 r ψ r-\psi 圖上,這也被稱爲Sweby圖,圖中的藍色區域即是TVD單調區域。使用這個圖的話,就很容易掌握TVD格式的公式。任何通量限制器的公式只要是能落在TVD單調區域內的都能產生出TVD格式來,Sweby圖跟前面講的NVD是很是相似的。

前面講過的全部格式的限制器也能夠推導出來,還能夠把它們的函數關係繪製在Sweby圖中。顯然,中心差分格式的限制器爲 ψ C D = 1 \psi^{CD}=1 ,而二階迎風格式的限制器爲
ϕ f = ϕ C + 1 2 ψ S O U ( r f ) ( ϕ D ϕ C ) = 3 2 ϕ C 1 2 ϕ U ψ S O U ( r f ) = ϕ C ϕ U ϕ D ϕ C = r f \begin{aligned} \phi_f&=\phi_C+\frac{1}{2}\psi^{SOU}(r_f)(\phi_D-\phi_C)=\frac32\phi_C-\frac12\phi_U \\ \Rightarrow \psi^{SOU}(r_f)&=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}=r_f \end{aligned}
中心差分格式和二階迎風格式的限制器展現在Sweby圖中,以下,CD爲中心差分格式,SOU爲二階迎風格式。
在這裏插入圖片描述

Sweby還注意到,因爲 r 0 r\le0 ψ ( r ) = 0 \psi(r)=0 ,因此在解中有極值的地方就沒有二階精度了。二階迎風和中心差分格式都是二階格式,從圖中可看到,他倆的限制器都通過了點 ( 1 , 1 ) (1,1) 。此外,Van Leer的工做代表,任何二階格式都能寫成中心差分格式和二階迎風格式的加權平均形式。因此呢,若是一個格式有二階精度,那麼它的限制器必需要穿過點 ( 1 , 1 ) (1,1) 纔對,並且,該限制器應該位於中心差分格式和二階迎風格式限制器所圍起來的區域,即上圖中的藍色區域纔對。若是要用上一小節的NVD(歸一化變量圖)來表示該區域,那麼是下面這樣子的。
在這裏插入圖片描述

一樣也能夠求出來其它高階格式的限制器來,即
{ U p w i n d ψ ( r f ) = 0 D o w n w i n d ψ ( r f ) = 2 F R O M M ψ ( r f ) = 1 + r f 2 S O U ψ ( r f ) = r f C D ψ ( r f ) = 1 Q U I C K ψ ( r f ) = 3 + r f 4 \begin{cases} Upwind & \psi(r_f)=0 \\ Downwind & \psi(r_f)=2 \\ FROMM & \psi(r_f)=\dfrac{1+r_f}{2} \\ SOU & \psi(r_f)=r_f \\ CD & \psi(r_f)=1 \\ QUICK & \psi(r_f)=\dfrac{3+r_f}{4} \end{cases}
這些限制器的函數關係以下圖所示,除了迎風格式的限制器外,全部其它格式的限制器都並不是徹底位於單調區間內,這也就是爲何這些格式都是無界的。
在這裏插入圖片描述

經過把這些不一樣格式的限制器 ψ ( r f ) \psi(r_f) 限制到單調區間內,即可把這些高階格式轉化爲高分辨率TVD格式。經過這種方式能夠發展不少種TVD格式,它們的限制器以下圖所示,它們限制器的函數關係以下。

在這裏插入圖片描述

{ S U P E R B E E ψ ( r f ) = max ( 0 , min ( 1 , 2 r f ) , min ( 2 , r f ) ) M I N M O D ψ ( r f ) = max ( 0 , min ( 1 , r f ) ) O S H E R ψ ( r f ) = max ( 0 , min ( 2 , r f ) ) V a n   L e e r ψ ( r f ) = r f + r f 1 + r f M U S C L ψ ( r f ) = max ( 0 , min ( 2 r f , ( r f + 1 ) / 2 , 2 ) ) \begin{cases} SUPERBEE & \psi(r_f)= \max(0,\min(1,2r_f),\min(2,r_f)) \\ MINMOD & \psi(r_f)=\max(0,\min(1,r_f)) \\ OSHER & \psi(r_f)= \max(0,\min(2,r_f)) \\ Van~Leer & \psi(r_f)=\dfrac{r_f+|r_f|}{1+|r_f|} \\ MUSCL & \psi(r_f)= \max(0,\min(2r_f,(r_f+1)/2, 2)) \end{cases}

5 NVF-TVD關係

既然NVF和TVD公式都是用來添加有界性的,只不過是採用不一樣的方法而已,那麼它倆應該是有某種關聯纔對啊。下面先推導出 r f r_f ϕ ~ C \tilde\phi_C 的關係式,而後比較NVF-CBC和TVD-CBC,最後導出通用的轉化關係,可讓任何TVD格式的函數關係寫成NVF框架的形式,或者是反過來,把NVF格式寫成TVD形式。

r f r_f ϕ ~ C \tilde\phi_C 的關係不難推導,直接由它們的定義可得
r f = ϕ C ϕ U ϕ D ϕ C ϕ ~ C = ϕ C ϕ U ϕ D ϕ U } r f = ( ϕ C ϕ U ) / ( ϕ D ϕ U ) ( ϕ D ϕ U + ϕ U ϕ C ) / ( ϕ D ϕ U ) = ϕ ~ C 1 ϕ ~ C ϕ ~ C = r f 1 + r f \left. \begin{aligned} r_f&=\frac{\phi_C-\phi_U}{\phi_D-\phi_C} \\ \tilde\phi_C&=\frac{\phi_C-\phi_U}{\phi_D-\phi_U} \end{aligned} \right\}\Rightarrow r_f=\frac{(\phi_C-\phi_U)/(\phi_D-\phi_U)}{(\phi_D-\phi_U+\phi_U-\phi_C)/(\phi_D-\phi_U)}= \frac{\tilde\phi_C}{1-\tilde\phi_C}\\ \Rightarrow \tilde\phi_C=\frac{r_f}{1+r_f}
使用這個式子,能夠比較兩種框架下的線性格式,好比迎風格式在TVD框架下的限制器爲 ψ ( r f ) = 0 \psi(r_f)=0 ,那麼其應該等效於NVF框架下的迎風格式 ϕ ~ f = ϕ ~ C \tilde\phi_f=\tilde\phi_C ,二者的等效性很容易證實
ψ ( r f ) = 0 ϕ f = ϕ C ϕ ~ f = ϕ ~ C \psi(r_f)=0\Rightarrow\phi_f=\phi_C\Rightarrow\tilde\phi_f=\tilde\phi_C
此外,無論是何種格式,在 r f 0 r_f\le0 時,TVD-CBC中的限制器 ψ ( r f ) = 0 \psi(r_f)=0 變爲純粹的迎風格式,有
r f 0 ϕ ~ C 1 ϕ ~ C 0 { ϕ ~ C 0 ϕ ~ C 1 r_f\le0 \Rightarrow \frac{\tilde\phi_C}{1-\tilde\phi_C}\le0 \Rightarrow \left\{\begin{aligned} \tilde\phi_C\le0 \\ \tilde\phi_C\ge1 \end{aligned}\right.
這跟在NVF-CBC中的迎風格式的使用區間是徹底一致的,即 ϕ ~ C 0 \tilde\phi_C\le0 ϕ ~ C 1 \tilde\phi_C\ge1 時,使用純粹的迎風格式。

還有,在NVF-CBC中,函數關係在區間 0 ϕ ~ C 1 0\le\tilde\phi_C\le1 是單調增長的,而在Sweby圖中該區域則是 0 r f + 0\le r_f \le +\infin ,這兩個區間其實是一樣的,以下所示:
ϕ ~ C 1 r f = ϕ ~ C 1 ϕ ~ C + \tilde\phi_C \rightarrow 1 \Rightarrow r_f=\frac{\tilde\phi_C}{1-\tilde\phi_C} \rightarrow +\infin
更進一步,對於TVD-CBC的條件
ψ ( r f ) 2 \psi(r_f)\le2
等效的NVF-CBC條件爲
ψ ( r f ) = 2 ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) } ϕ f = ϕ C + ( ϕ D ϕ C ) = ϕ D ϕ ~ f = 1 \left. \begin{aligned} &\psi(r_f)=2\\ &\phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C) \end{aligned} \right\} \Rightarrow \phi_f=\phi_C+(\phi_D-\phi_C)=\phi_D \Rightarrow \tilde\phi_f=1

ψ ( r f ) 2 ϕ ~ f 1 \psi(r_f)\le2 \Rightarrow \tilde\phi_f\le1
正是NVF-CBC所應知足的條件。TVD-CBC的限制器 ψ ( r f ) \psi(r_f) 應知足的最後一個條件是
ψ ( r f ) 2 r f \psi(r_f) \le 2r_f
等效的NVF-CBC條件爲
ψ ( r f ) = 2 r f ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) } ϕ f = ϕ C + ϕ C ϕ U ϕ D ϕ C ( ϕ D ϕ C ) = 2 ϕ C ϕ U \left. \begin{aligned} & \psi(r_f) = 2r_f\\ &\phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C) \end{aligned} \right\} \Rightarrow \phi_f=\phi_C+\frac{\phi_C-\phi_U}{\phi_D-\phi_C}(\phi_D-\phi_C)=2\phi_C-\phi_U
歸一化後,爲
ϕ f = 2 ϕ C ϕ U ϕ f ϕ U = 2 ϕ C 2 ϕ U ϕ ~ f = 2 ϕ ~ C \phi_f=2\phi_C-\phi_U \Rightarrow \phi_f-\phi_U=2\phi_C-2\phi_U \Rightarrow \tilde\phi_f=2\tilde\phi_C
這比NVF-CBC要更加嚴格,這是兩種方法的惟一不一樣之處。基於該條件,TVD-CBC和修正的NVF-CBC(用TVD條件表示的NVF)以下圖a所示,修正NVF-CBC的單調區間減少成了迎風線和藍色線所圍的區域;而修正TVD-CBC(用NVF條件表示的TVD)和NVF-CBC以下圖b所示(NVF-CBC中的線 ϕ ~ C = 0 \tilde\phi_C=0 對應TVD-CBC中的線 r f = 0 r_f=0 )。再考慮到二階精度的要求,TVD格式中的二階精度要求必須經過點 ( 1 , 1 ) (1,1) ,即, ψ ( 1 ) = 1 \psi(1)=1 ,那麼其等效的NVF值爲
r f = 1 ϕ ~ C 1 ϕ ~ C = 1 ϕ ~ C = 0.5 ψ ( 1 ) = 1 ϕ f = ϕ C + 1 2 ψ ( 1 ) ( ϕ D ϕ C ) = 1 2 ( ϕ D + ϕ C ) ϕ f ϕ U = 1 2 ( ϕ D ϕ U + ϕ C ϕ U ) ϕ ~ f = 1 2 ( 1 + ϕ ~ C ) = 0.75 \begin{aligned} r_f=1 &\Rightarrow \frac{\tilde\phi_C}{1-\tilde\phi_C}=1 \Rightarrow \tilde\phi_C=0.5 \\ \psi(1)=1 &\Rightarrow \phi_f=\phi_C+\frac{1}{2}\psi(1)(\phi_D-\phi_C)= \frac{1}{2}(\phi_D+\phi_C) \\ & \Rightarrow \phi_f-\phi_U=\frac{1}{2}(\phi_D-\phi_U+\phi_C-\phi_U) \Rightarrow \tilde\phi_f=\frac{1}{2}(1+\tilde\phi_C)=0.75 \end{aligned}
這偏偏就是NVF中的點 ( 0.5 , 0.75 ) (0.5,0.75) 。前面提到過,Van Leer證實了任何二階格式均可以寫成中心差分格式和二階迎風格式的加權平均形式,所以其函數關係應該位於CD和SOU格式的函數關係之間,以下圖c所示。

至此,充分說明了NVF-CBC和TVD-CBC實際上就是一回事,即,異曲同工。
在這裏插入圖片描述

上面的過程也能夠用來把TVD格式轉化爲等效的NVF格式,或者是反過來作轉化。用NVF框架下的格式爲例,面 f f 上的值 ϕ f \phi_f 表示爲
ϕ f = f ( ϕ ~ C ) ( ϕ D ϕ U ) + ϕ U ϕ ~ C = ϕ C ϕ U ϕ D ϕ C \phi_f=f(\tilde\phi_C)(\phi_D-\phi_U)+\phi_U\quad\quad\tilde\phi_C=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
TVD格式的 ϕ f \phi_f 表達式爲
ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) r f = ϕ C ϕ U ϕ D ϕ C \phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)\quad\quad r_f=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
結合上面兩式,得
ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) = f ( ϕ ~ C ) ( ϕ D ϕ U ) + ϕ U \phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)=f(\tilde\phi_C)(\phi_D-\phi_U)+\phi_U

ψ ( r f ) ( ϕ D ϕ C ) ( ϕ D ϕ U ) = 2 f ( ϕ ~ C ) ( ϕ D ϕ U ) ( ϕ D ϕ U ) 2 ( ϕ C ϕ U ) ( ϕ D ϕ U ) = 2 ( f ( ϕ ~ C ) ϕ ~ C ) \psi(r_f)\frac{(\phi_D-\phi_C)}{(\phi_D-\phi_U)}= 2f(\tilde\phi_C)\frac{(\phi_D-\phi_U)}{(\phi_D-\phi_U)}- 2\frac{(\phi_C-\phi_U)}{(\phi_D-\phi_U)} =2(f(\tilde\phi_C)-\tilde\phi_C)
上式左端項能夠寫爲
ψ ( r f ) ( ϕ D ϕ C ) ( ϕ D ϕ U ) = ψ ( r f ) ( ϕ D ϕ U ϕ C + ϕ U ) ( ϕ D ϕ U ) = ψ ( r f ) ( 1 ϕ ~ C ) \psi(r_f)\frac{(\phi_D-\phi_C)}{(\phi_D-\phi_U)}= \psi(r_f)\frac{(\phi_D-\phi_U-\phi_C+\phi_U)}{(\phi_D-\phi_U)}= \psi(r_f)(1-\tilde\phi_C)
因而
ψ ( r f ) ( 1 ϕ ~ C ) = 2 ( f ( ϕ ~ C ) ϕ ~ C ) ψ ( r f ) = 2 f ( ϕ ~ C ) ϕ ~ C 1 ϕ ~ C \psi(r_f)(1-\tilde\phi_C)=2(f(\tilde\phi_C)-\tilde\phi_C) \Rightarrow \psi(r_f)=2\frac{f(\tilde\phi_C)-\tilde\phi_C}{1-\tilde\phi_C}
上式也可寫成
f ( ϕ ~ C ) = ψ ( r f ) + 2 r f 2 ( 1 + r f ) f(\tilde\phi_C)=\frac{\psi(r_f)+2r_f}{2(1+r_f)}
這樣便推得了NVF格式與TVD格式的轉換關係。

舉個例子,迎風格式在NVF框架下的函數關係是 ϕ ~ f = ϕ ~ C \tilde\phi_f=\tilde\phi_C ,則其TVD限制器爲
ϕ ~ f = ϕ ~ C ψ ( r f ) = 2 f ( ϕ ~ C ) ϕ ~ C 1 ϕ ~ C = 2 ϕ ~ C ϕ ~ C 1 ϕ ~ C = 0 \tilde\phi_f=\tilde\phi_C \Rightarrow \psi(r_f)=2\frac{f(\tilde\phi_C)-\tilde\phi_C}{1-\tilde\phi_C}= 2\frac{\tilde\phi_C-\tilde\phi_C}{1-\tilde\phi_C}=0

相關文章
相關標籤/搜索