基于BVD原理的高保真空間重構方法

肖 鋒

(東京工業大學, 日本 東京 152-8550)

0 引 言

可壓縮流體的運動變化規律是空氣動力學的主要研究內容。伴隨高馬赫數流動產生的激波等不連續現象使得可壓縮流動的流場結構更為復雜，并同時包含連續與間斷解，給理論和實驗研究帶來很多本質性的困難。在過去的數十年中，由于計算機硬件的飛速發展，以及相關領域實際應用的巨大需求，計算流體力學作為最活躍的科學研究領域之一，取得了重大進展， 并已成為空氣動力學非常重要的研究手段。

針對描述可壓縮流體的歐拉方程，至今已提出了各種數值算法。其中，以具有嚴格守恒特性的Godunov格式[1]為基礎開發的各類高階守恒格式已成為求解可壓縮流體運動的主流算法?；贕odunov格式的高階守恒格式一般包括重構和演化兩個重要步驟。其中，演化是根據重構得到的離散變量，利用控制方程的物理性質計算相應時間步長內的數值通量，進而預報下一時刻的物理量變化。關于歐拉方程及其黏性擴展情形下的數值通量計算，可參閱有關專著及論文[2-5]。本文集中討論如何進行物理變量的空間重構。

作為構造高階Godunov方法的關鍵，空間重構一直是關注的焦點。至今已提出了針對不同空間離散框架的各種重構方法[6]，包括有限體積或有限差分格式[7-19]、緊致格式[20-22]等基于單個單元自由度的方法， 也有間斷伽列金方法[23-24]、間斷伽列金/有限體積混合方法[25-27]、通量重構方法[28-30]、多矩有限體積方法[31-34]等包含多個局地自由度的高精度方法，并已被用于求解各類實際問題。這些方法通常采用高階多項式進行重構，在用于計算解較為光滑的問題時，能夠得到較高收斂率和較精確的數值結果。然而，當求解的物理問題含有激波或多相復雜介質界面等強間斷時，必須使用非線性限制映射或限制器來抑制虛假的數值振蕩。其中，使用WENO(weighted essentially non-oscillatory)[12-13]思想構造的高分辨格式既可有效地抑制數值振蕩又能對光滑解獲得較高的收斂精度，在可壓縮流體的數值模擬中得到廣泛應用?，F有算法中設計非線性限制器的基本做法是，在解較光滑的區域盡量保持高階多項式的性質，而在解出現間斷處采用更平直或較低階的重構函數。這類做法通常在抑制數值振蕩與控制數值耗散之間，在計算魯棒性和求解精度等方面無法兼顧，很難對光滑和間斷解同時給出令人滿意的計算結果。在實際應用中主要表現為：(1)存在顯著的數值耗散，會抹平接觸間斷和較小尺度的流場結構；(2)多數非線性權函數無法使重構函數在光滑情形下完全恢復到原來的常系數高階多項式插值函數，且依賴于一些人為確定的參數，給實際應用帶來不便。作為相關領域長期以來一個未能很好解決的問題，如何構造能同時精確捕捉光滑和間斷解的高保真空間重構方法，一直是可壓縮流體數值方法研究的熱點。

為了更好地模擬包含間斷解的可壓縮復雜流動，我們近年提出了以減小空間重構在網格單元邊界上的變差為原則設計數值格式的思想，即，BVD(Boundary Variation Diminishing)原理。根據此原理，我們發展了一類可用于求解單相和多相可壓縮流體的高保真數值算法[35-46]。這些算法能有效克服其他現有高分辨格式存在的問題，可以同時精確捕捉光滑和間斷解。各種標準算例結果表明，BVD格式在模擬各種具有強間斷的可壓縮流體運動方面，與其他的現有格式相比具有明顯的優勢。 利用BVD原理進行空間重構可望成為構建可壓縮流體新型高保真數值方法的有效途徑。

本文對BVD格式研發的已有工作做一個簡要綜述，內容包括BVD的基本思想介紹、幾個具有實用意義BVD算法、有代表性的標準算例結果，以及簡短總結和展望。

1 BVD的基本思想

1.1 Godunov方法概要

我們考慮雙曲標量守恒律:

通過以上步驟(I)和(II)獲得邊界數值通量后，一般可以通過線方法，采用常微分方程的數值解法對半離散方程(2)進行時間積分,得到下一時刻的數值解。

至今的高分辨數值方法研究，基本集中于開發具備更好數值性能的空間重構(I)和黎曼算子(II)。

1.2 BVD方法的基本思想

從近似黎曼算子的基本形式(3)可以看出，數值通量一般由中心格式(右端第一項)和數值耗散/數值黏性(右端第二項)兩部分組成。如果在構造數值格式時，能夠減小數值耗散項，將有效降低格式的數值耗散，從而改善計算結果。

由數值耗散項的表達式不難看出，減小數值耗散可從兩方面入手，即:

(B) 減小空間重構得到的邊界變差值

各類黎曼算子的研究開發在一定意義上可以看作是為實現目標(A)所做的努力。圍繞在保證計算穩定的條件下，如何控制減小數值耗散系數，提出了各類近似黎曼算子。這方面的研究可參閱相關專著[2]。

為了在空間重構過程中有效地減小數值耗散，我們提出以下設計空間重構格式的一般性原理。

如果在開發空間重構算法時能有意識地根據BVD原理，盡可能地減小邊界變差值，則可以有效抑制格式的數值耗散，改進數值計算結果。我們注意到，文獻[47]通過邊界變差最小化的變分約束關系構造了緊致有限體積方法。

2 BVD混合格式的構建

基于BVD原理，我們提出了構建高保真數值格式的基本框架，并發展了一類基于BVD原理的混合格式。BVD格式主要包括BVD許容重構函數集和BVD選擇算法兩個部分。

2.1 BVD許容重構函數集

流體力學問題解的空間分布極為復雜，既包括各種尺度的渦和波動等相對光滑的流場結構，也包括激波、接觸間斷、不同流體間的物質界面等不連續結構。很難用一類插值函數準確描述各種不同的流場結構。在以往的研究中，各種重構函數和重構方法被相繼開發出來，并用于針對不同應用問題的高分辨格式。我們把一些有代表性的重構函數分為以下三類BVD許容函數。

2.1.1 第一類：常系數多項式

這類重構函數用于逼近光滑解。包括各階泰勒多項式、緊致格式[48-49]以及各種頻散/耗散優化格式[50-54]等。作為例子，我們具體給出針對網格單元Ωi的四階多項式:

其中系數ck(k=0,1,…,4)通過以下約束條件求出:

(5)

為書寫簡便，我們在此采用等網格距，Δxi=Δx。由此可以直接得到:

2.1.2 第二類：非線性系數多項式

在高分辨激波捕捉格式的發展過程中，為避免高階多項式在間斷解附近出現的非物理振蕩，提出了許多非線性限制投影算法。例如，采用TVD限制器的2階MUSCL格式[7]，可將網格單元邊界左右兩側的變量值寫為：

如果采用WENO類算法，則可將插值函數寫成非線性加權系數多項式的形式。以下給出5階WENO格式計算網格單元邊界左右兩側變量值的算式。

其中,

關于非線性權重系數ωj(j=0,1,2)的計算方法，至今已提出了各種方案[12-19,24]。

2.1.3 第三類:非多項式類插值函數

采用某些具有單調特性的非多項式類插值函數能夠避免空間重構中的數值振蕩。例如，采用雙曲函數、有理函數或對數函數[55-57]可以有效地抑制數值振蕩。文獻[58-61]使用雙曲正切函數擬合躍階型間斷，并發展了一類可用于捕捉移動界面的代數VOF(Volume of Fluid)方法[62]，稱為THINC(Tangent of Hyperbola for INterface Capturing)方法。和其他Sigmoid函數一樣，THINC重構函數具有嚴格的單調性，并能很好地模擬激波、接觸間斷以及多相流中的物質界面等具有躍階分布的變量。

由此，可以將通過THINC重構得到的網格單元左右兩端的值寫為:

這里，

B=eθβ(2C-1)

在實際應用中，可以通過調節陡度參數β的值，使THINC重構函數適用于各類數值解分布。數值頻散/耗散分析表明，如果β取值于1和1.3之間，則THINC重構可以表示一類MUSCL格式[37]。取較大β值，則可以很好地擬合躍階間斷分布，并有效地抑制數值擴散。

對于多維情形，可將THINC重構函數表示為:

其中，Pi(x,y,z)+di=0為界面方程。di通過與式(9)類似的多維體積分約束關系計算。Pi(x,y,z)可采用高階多項式，以更準確地描述界面的幾何形狀。關于多維高階THINC重構可參閱相關文獻[60-61]。

從以上三類許容函數中選取適當的候補函數構成BVD許容函數集，記為

ξ=1,2,…,Ξ

這里Ξ表示BVD許容函數集里所包含的候補函數的個數，通常Ξ≥2。

一般根據流場結構的主要特征以及希望改善的數值解性質選取候補函數。例如，針對以間斷解為主的應用問題，可以選取非線性限制多項式(第二類)和THINC函數(第三類)構成BVD許容函數集[35-36,43]。或選取采用不同β值的THINC函數構成BVD許容函數集[37]。針對以聲場或旋渦為主的應用問題，為了充分利用高階常系數多項式的良好性質，可以避開非線性限制多項式，只選取高階常系數多項式(第一類)和THINC函數(第三類)構成BVD許容函數集[38-42,44-46]。研究表明，通過構造合適的BVD算法，即使不采用非線性限制多項式，也能有效地抑制數值振蕩，如后敘的PnTm-BVD格式。

2.2 BVD選擇算法

我們對各類BVD許容函數用于連續和間斷分布重構時產生的BV或TBV值進行了分析及數值實驗，得到以下主要結論：

(1) 對于光滑分布，采用第一類許容重構函數(常系數多項式)能夠得到更小的BV或TBV值。且多項式階數越高，邊界變差值越小。與Weierstrass多項式逼近定理一致。

(2) 對于躍階型間斷解，在第二類許容重構函數(帶限制器的非線性系數多項式)中，高階格式的TBV值小于低階格式的TBV值。數值耗散較大的格式產生較大的TBV值。與基于非線性系數多項式重構相比，采用較大β值的THINC重構可以獲得更小的TBV值。

(3) 對間斷進行插值時，THINC函數產生的BV只出現在空間較窄區域，而基于多項式的重構函數則會導致較寬區域內出現明顯的BV。

由此，我們提出了各種選擇算法。下面列舉幾個有較強實用意義的BVD算法。

2.2.1 BVD算法Ⅰ

BVD許容函數集包括兩個候補函數，即

2.2.2 BVD算法Ⅱ

采用包含兩個候補函數的BVD許容函數集，通過以下步驟確定重構函數。

2.2.3 BVD算法Ⅲ

針對包含兩個候補函數的BVD許容函數集，通過以下步驟確定重構函數。

(2) 選取使得邊界總變差最小的函數作為最終的重構函數，

針對多個候補函數構成的BVD許容函數集(Ξ≥3)，算法III可寫成更一般的形式:

文獻[43]采用三個候補函數建立了非結構網格的BVD格式，能很好地捕捉間斷和渦結構。

2.2.4 BVD算法Ⅳ

當BVD許容函數集包括3個以上候補函數時，可通過多步BVD算法確定重構函數。其一般步驟如下:

for (ξ=2;ξ≤Ξ;ξ++)

{

}

其中，BVDξ(·,·)表示基于BVD原理，針對兩個候補函數的選擇算法，如上述的算法Ⅰ至算法Ⅲ。

(2) 對ξ=1，…,m-1循環計算

和

j=i-1,i,i+1。

文獻[38-40]分別給出了P4T2-BVD、P6T3-BVD、P8T3-BVD、P10T3-BVD、P12T3-BVD和P14T3-BVD等高階BVD格式。這些格式對光滑解可以得到對應的常系數多項式的最高收斂階數，亦即，PnTm-BVD格式具有n+1階。同時，能有效避免數值振蕩。與其他現有的高階方法不同，這類BVD格式完全不需要使用傳統的非線性限制映射來消除數值振蕩。

使用ADR(Approximate Dispersion Relation)方法[63]可以分析數值格式的頻散和耗散性質。圖1顯示各階PnTm-BVD格式的數值頻散和耗散特征。為了比較，這里還顯示了WENO-JS[13]、WENO-M[14]、WENO-Z[15]以及基于各階格式對應的常系數多項式的線性格式的頻散/耗散關系。對于所有波數，PnTm-BVD可以回復其對應高階線性格式的頻散/耗散關系，而采用非線性限制多項式的高階格式則在高波數區明顯偏離線性格式。

(a) 頻散特征

(b) 耗散特征

此外，文獻[44]還提出了以4次常系數多項式和采用可變β值的THINC作為候補函數的BVD格式。其中，網格單元Ωi的THINC函數的β值通過以下算式給出:

j=i-1,i,i+1

該格式稱為P4Tβv-BVD格式。

2.3 多維BVD格式的構建

上述構建一維BVD格式的思路原則上都可以推廣到多維情形。對于結構網格可以采用方向分離方法，通過一維格式進行空間重構。對于非結構網格，文獻[43]用重構函數沿網格單元各邊界的積分值計算BV。針對網格單元Ωi和與其相鄰的單元Ωij，沿單元邊界Γij的BV值由下式計算，

需要指出，BVD算法在同一網格單元上要對所有候補函數進行重構，并儲存相應的BV或TBV值，會增加計算量和存儲量。然而，由于各網格單元間的計算過程相互獨立，BVD算法不會給并行處理帶來不利影響。

3 數值結果

我們利用標量守恒律和歐拉方程的各種標準算例，對BVD算法進行了系統的檢驗，并與其他現有方法進行了比較。

3.1 精度檢驗

對于光滑解，BVD格式一般都能從候補函數中選擇高階重構，從而獲得高收斂率。以下列舉采用不同階數常系數多項式作為候補函數的BVD格式的收斂檢驗算例。文獻[38-40]分別針對平流和歐拉方程進行了網格收斂率檢驗。對于一維平流方程，選用以下具有一定光滑度的初始條件，

對于二維歐拉方程的密度擾動傳播算例,初始條件設為:

(x,y)∈[-1,1]×[-1,1]

表1給出各階PnTm-BVD格式的誤差及收斂率，均能達到相對應的常系數多項式的階數。其誤差的絕對值與采用常系數多項式作為重構函數的結果基本一致。

表1 PnTm-BVD格式的數值誤差及精度檢驗

3.2 Jiang-Shu平流實驗

作為檢驗數值格式是否能夠同時捕捉連續和間斷解的算例，文獻[13]設計了一個既包括光滑分布又包括間斷解的平流初始條件。

圖2分別顯示了經過一百萬步后WENO-JS[13]，WENO-Z[15]和P4T2-BVD格式[38]的計算結果。可以看出，由于數值誤差的不斷積累，WENO-JS和WENO-Z格式已無法保持初始分布的特征。而P4T2-BVD格式則能很好保持初始分布，對連續和間斷解都能得到高保真的計算結果。P4T2-BVD格式既能有效地避免數值振蕩，又可以消除數值耗散，使解的結構在長期積分后仍能可靠再現。這是現有其他方法很難做到的。

圖2 Jiang-Shu平流標準算例在t=2000(1×106步)時的計算結果. 從左至右分別為五階WENO-JS[13]格式、WENO-Z[15]格式和P4T2-BVD格式[38]的計算結果，網格單元數為400.

3.3 一維歐拉方程的數值實驗

對于一維歐拉方程，我們利用各類標準算例對BVD格式進行了檢驗，在抑制數值振蕩的同時，對間斷和光滑解均能得到高保真的計算結果。以下給出幾個具體算例的結果。

3.3.1 雙激波相互作用[64]

本算例通過兩個激波的相互作用及邊界反射產生復雜的流場結構，包括激波、接觸間斷和膨脹波，被廣泛用于檢驗數值方法是否產生數值振蕩，以及對各類間斷和連續解的捕捉能力。圖3給出了P4T2-BVD以及P4Tβv-BVD格式的計算結果?？梢钥吹紹VD格式能夠精確地得到各類特征的數值解。在圖4中我們進一步顯示9階和11階BVD格式的計算結果，高階BVD格式也能有效抑制數值振蕩，保持良好的魯棒性，并得到高保真的數值計算結果。在上述結果中，BVD格式能夠有效地控制數值耗散，可使最左端的接觸間斷厚度保持在3至4個網格之內。這是現有基于多項式重構及非線性限制映射的高分辨格式難以實現的。

圖3 Two interacting blast waves標準算例在t=0.038時刻的P4T2-BVD和P4Tβv-BVD格式的計算結果(密度), 網格單元數為400.

圖4 Two interacting blast waves標準算例在t=0.038時刻的P8T3-BVD(左)和P10T3-BVD(右)格式的計算結果(密度),網格單元數為400.

3.3.2 激波與密度擾動相互作用[65]

為了進一步檢驗BVD格式對高頻周期波的捕捉能力，我們計算了激波與密度擾動相互作用的算例。圖5顯示不同階數BVD格式的計算結果?？梢钥吹剑煌A數的BVD格式都能抑制數值振蕩，并在光滑區正確地選擇多項式作為重構函數。隨著多項式階數的提高，BVD能有效減小數值耗散，更精確地捕捉高頻波動[40]。

(a) P4T2-BVD[38]和P6T3-BVD格式[38]的計算結果

(b) P8T3-BVD[38]和P10T3-BVD格式[38]的計算結果

3.3.3 Le Blanc 激波管問題[66]

作為一個較困難的標準算例，Le Blanc激波管問題的解包括非常強的激波和稀疏波，經常用來檢驗算法的魯棒性[66]。圖6給出了P4T2-BVD的計算結果(密度)。沒有出現明顯的數值振蕩。并能很好再現稀疏波、接觸間斷和激波[46]。數值實驗表明，隨著網格分辨率的提高，右行激波能夠收斂到相應精確解的位置。

圖6 Le Blanc一維激波管標準算例在t=6時刻的P4T2-BVD格式的計算結果(密度), 網格單元數為800.

3.3.4 剛性爆轟波[67]

與現有的數值方法相比，BVD格式能夠消除數值耗散，精確地捕捉不連續界面。圖7顯示了包括剛性化學反應的爆轟波計算結果。P4T2-BVD格式能有效控制反應面的溫度躍階厚度，從而保證反應面始終以正確的速度傳播[42,46]。

圖7 剛性爆轟波標準算例在t=π/5時刻的P4T2-BVD格式的計算結果(密度),網格單元數為200.

3.4 二維歐拉方程的數值實驗

3.4.1 雙馬赫反射[64]

作為檢驗二維算法的標準算例，雙馬赫反射問題自提出以來，在相關文獻中得到廣泛應用。在本算例中，沿兩個反射激波間的滑移線會出現開爾文-亥姆霍茲不穩定，進而產生渦列。在數值計算中，這些渦列是否出現和發展取決于數值耗散的大小。圖8分別給出了P4-T2-BVD和P14-T3-BVD格式的計算結果。由于BVD原理旨在減小格式的數值耗散，以此為基礎設計的格式都能更多地保留流場的細微結構，捕捉到更小尺度的流體運動。與其他五階格式相比，P4-T2-BVD格式的結果顯示了充分發展的渦列。隨著精度提高，十五階的P14-T3-BVD格式可以捕捉到更多旋渦結構[44]。

3.4.2 爆轟波擾流問題[68]

二維爆轟波的90°凸角繞流在繞角下流區域會產生密度和壓力的低值區，極易在數值結果中出現負值，導致計算失敗。文獻[46]在P4-T2-BVD格式中引入了MOOD(Multi dimensional Optimal Order Detection)正定修正算法[69]，從而保證了熱力學變量數值解的正定性。圖9顯示二維爆轟波90°凸角繞流的密度和溫度分布，對于反應區和各種流場結構均獲得較精確的計算結果。

(a) P4T2-BVD

3.4.3 空氣激波與水柱相互作用[70]

采用THINC函數作為候補重構函數的BVD格式可以有效控制數值耗散，從而避免接觸間斷和物質界面在計算過程中被平滑擴散的現象。文獻[36]采用MUSCL-THINC-BVD格式對各類兩項可壓縮流動進行了模擬，取得了良好的計算結果。

圖10給出了空氣與水柱相互作用的計算結果。BVD格式能有效地保持水/氣界面的厚度，并對各種流場結構都能得到高保真的計算結果。與文獻[70]及其他現有算法相比，BVD格式能在抑制數值振蕩的同時，獲得間斷解的高保真計算結果，顯示出在模擬可壓縮多相介質等具有顯著間斷解結構特征應用問題時的明顯優勢和潛力。文獻[43]將基于MUSCL和THINC的BVD格式擴展到非結構網格，并計算了包含移動界面的可壓縮多相流算例，獲得了良好的數值結果。

(a) 密度

(b) 溫度

圖10 空氣與水柱相互作用的兩相流算例在t=2.15時刻MUSCL-THINC-BVD格式的計算結果(密度梯度), 網格單元數為2000×500.

4 結 論

本文概述了BVD方法的基本思想，設計相關格式的基本思路，以及一些具有很強實用價值的BVD格式。并通過單相和兩相可壓縮流動的一些典型算例驗證了BVD格式的特點和優勢。

BVD方法把減小數值耗散作為設計空間重構算法的指導思想，具有明確的物理含義和普遍性，是一條設計和改進數值格式的新途徑。針對研究對象，選用合適的BVD許容函數作為空間重構的候補選項，并設計相應的BVD算法可以構建一類全新的數值格式。這些格式與基于多項式重構和非線性限制修正的傳統方法有本質區別。數值實驗結果表明，BVD格式可以得到明顯優于已有算法的計算結果，能同時對各種尺度的連續和間斷結構都獲得高保真的數值解。

采用THINC或類似的Sigmoid函數作為候補重構函數可以明顯改善基于歐拉網格的數值方法對接觸間斷和物質界面在計算過程中被平滑擴散的問題，對解決含有化學反應面或移動物質界面流動問題的實際應用具有重要意義。

相關研究工作仍在繼續，例如，可以通過BVD判據，在光滑區引入中心格式進一步控制數值耗散[71]，以及將BVD格式擴展到間斷伽列金的數值框架[72]。希望通過不斷努力和持續發展，能夠提出具有更好性能的BVD格式，并成為解決實際應用問題的有力工具。

致謝:鄧希、謝彬、孫紫堯、程李東、金鵬、陳春剛、Shimizu、Inaba、Tann、Wakimura等東京工業大學工學院肖研究室的學生以及合作者姜振華、Loubère、Abe等對本項研究工作做出了重要貢獻，特此致謝。