基于BVD原理的高保真空間重構(gòu)方法

肖 鋒

(東京工業(yè)大學(xué), 日本 東京 152-8550)

0 引 言

可壓縮流體的運(yùn)動(dòng)變化規(guī)律是空氣動(dòng)力學(xué)的主要研究?jī)?nèi)容。伴隨高馬赫數(shù)流動(dòng)產(chǎn)生的激波等不連續(xù)現(xiàn)象使得可壓縮流動(dòng)的流場(chǎng)結(jié)構(gòu)更為復(fù)雜，并同時(shí)包含連續(xù)與間斷解，給理論和實(shí)驗(yàn)研究帶來(lái)很多本質(zhì)性的困難。在過(guò)去的數(shù)十年中，由于計(jì)算機(jī)硬件的飛速發(fā)展，以及相關(guān)領(lǐng)域?qū)嶋H應(yīng)用的巨大需求，計(jì)算流體力學(xué)作為最活躍的科學(xué)研究領(lǐng)域之一，取得了重大進(jìn)展， 并已成為空氣動(dòng)力學(xué)非常重要的研究手段。

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

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

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

本文對(duì)BVD格式研發(fā)的已有工作做一個(gè)簡(jiǎn)要綜述，內(nèi)容包括BVD的基本思想介紹、幾個(gè)具有實(shí)用意義BVD算法、有代表性的標(biāo)準(zhǔn)算例結(jié)果，以及簡(jiǎn)短總結(jié)和展望。

1 BVD的基本思想

1.1 Godunov方法概要

我們考慮雙曲標(biāo)量守恒律:

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

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

1.2 BVD方法的基本思想

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

由數(shù)值耗散項(xiàng)的表達(dá)式不難看出，減小數(shù)值耗散可從兩方面入手，即:

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

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

為了在空間重構(gòu)過(guò)程中有效地減小數(shù)值耗散，我們提出以下設(shè)計(jì)空間重構(gòu)格式的一般性原理。

如果在開(kāi)發(fā)空間重構(gòu)算法時(shí)能有意識(shí)地根據(jù)BVD原理，盡可能地減小邊界變差值，則可以有效抑制格式的數(shù)值耗散，改進(jìn)數(shù)值計(jì)算結(jié)果。我們注意到，文獻(xiàn)[47]通過(guò)邊界變差最小化的變分約束關(guān)系構(gòu)造了緊致有限體積方法。

2 BVD混合格式的構(gòu)建

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

2.1 BVD許容重構(gòu)函數(shù)集

流體力學(xué)問(wèn)題解的空間分布極為復(fù)雜，既包括各種尺度的渦和波動(dòng)等相對(duì)光滑的流場(chǎng)結(jié)構(gòu)，也包括激波、接觸間斷、不同流體間的物質(zhì)界面等不連續(xù)結(jié)構(gòu)。很難用一類插值函數(shù)準(zhǔn)確描述各種不同的流場(chǎng)結(jié)構(gòu)。在以往的研究中，各種重構(gòu)函數(shù)和重構(gòu)方法被相繼開(kāi)發(fā)出來(lái)，并用于針對(duì)不同應(yīng)用問(wèn)題的高分辨格式。我們把一些有代表性的重構(gòu)函數(shù)分為以下三類BVD許容函數(shù)。

2.1.1 第一類：常系數(shù)多項(xiàng)式

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

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

(5)

為書寫簡(jiǎn)便，我們?cè)诖瞬捎玫染W(wǎng)格距，Δxi=Δx。由此可以直接得到:

2.1.2 第二類：非線性系數(shù)多項(xiàng)式

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

如果采用WENO類算法，則可將插值函數(shù)寫成非線性加權(quán)系數(shù)多項(xiàng)式的形式。以下給出5階WENO格式計(jì)算網(wǎng)格單元邊界左右兩側(cè)變量值的算式。

其中,

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

2.1.3 第三類:非多項(xiàng)式類插值函數(shù)

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

由此，可以將通過(guò)THINC重構(gòu)得到的網(wǎng)格單元左右兩端的值寫為:

這里，

B=eθβ(2C-1)

在實(shí)際應(yīng)用中，可以通過(guò)調(diào)節(jié)陡度參數(shù)β的值，使THINC重構(gòu)函數(shù)適用于各類數(shù)值解分布。數(shù)值頻散/耗散分析表明，如果β取值于1和1.3之間，則THINC重構(gòu)可以表示一類MUSCL格式[37]。取較大β值，則可以很好地?cái)M合躍階間斷分布，并有效地抑制數(shù)值擴(kuò)散。

對(duì)于多維情形，可將THINC重構(gòu)函數(shù)表示為:

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

從以上三類許容函數(shù)中選取適當(dāng)?shù)暮蜓a(bǔ)函數(shù)構(gòu)成BVD許容函數(shù)集，記為

ξ=1,2,…,Ξ

這里Ξ表示BVD許容函數(shù)集里所包含的候補(bǔ)函數(shù)的個(gè)數(shù)，通常Ξ≥2。

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

2.2 BVD選擇算法

我們對(duì)各類BVD許容函數(shù)用于連續(xù)和間斷分布重構(gòu)時(shí)產(chǎn)生的BV或TBV值進(jìn)行了分析及數(shù)值實(shí)驗(yàn)，得到以下主要結(jié)論：

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

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

(3) 對(duì)間斷進(jìn)行插值時(shí)，THINC函數(shù)產(chǎn)生的BV只出現(xiàn)在空間較窄區(qū)域，而基于多項(xiàng)式的重構(gòu)函數(shù)則會(huì)導(dǎo)致較寬區(qū)域內(nèi)出現(xiàn)明顯的BV。

由此，我們提出了各種選擇算法。下面列舉幾個(gè)有較強(qiáng)實(shí)用意義的BVD算法。

2.2.1 BVD算法Ⅰ

BVD許容函數(shù)集包括兩個(gè)候補(bǔ)函數(shù)，即

2.2.2 BVD算法Ⅱ

采用包含兩個(gè)候補(bǔ)函數(shù)的BVD許容函數(shù)集，通過(guò)以下步驟確定重構(gòu)函數(shù)。

2.2.3 BVD算法Ⅲ

針對(duì)包含兩個(gè)候補(bǔ)函數(shù)的BVD許容函數(shù)集，通過(guò)以下步驟確定重構(gòu)函數(shù)。

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

針對(duì)多個(gè)候補(bǔ)函數(shù)構(gòu)成的BVD許容函數(shù)集(Ξ≥3)，算法III可寫成更一般的形式:

文獻(xiàn)[43]采用三個(gè)候補(bǔ)函數(shù)建立了非結(jié)構(gòu)網(wǎng)格的BVD格式，能很好地捕捉間斷和渦結(jié)構(gòu)。

2.2.4 BVD算法Ⅳ

當(dāng)BVD許容函數(shù)集包括3個(gè)以上候補(bǔ)函數(shù)時(shí)，可通過(guò)多步BVD算法確定重構(gòu)函數(shù)。其一般步驟如下:

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

{

}

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

(2) 對(duì)ξ=1，…,m-1循環(huán)計(jì)算

和

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

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

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

(a) 頻散特征

(b) 耗散特征

此外，文獻(xiàn)[44]還提出了以4次常系數(shù)多項(xiàng)式和采用可變?chǔ)轮档腡HINC作為候補(bǔ)函數(shù)的BVD格式。其中，網(wǎng)格單元Ωi的THINC函數(shù)的β值通過(guò)以下算式給出:

j=i-1,i,i+1

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

2.3 多維BVD格式的構(gòu)建

上述構(gòu)建一維BVD格式的思路原則上都可以推廣到多維情形。對(duì)于結(jié)構(gòu)網(wǎng)格可以采用方向分離方法，通過(guò)一維格式進(jìn)行空間重構(gòu)。對(duì)于非結(jié)構(gòu)網(wǎng)格，文獻(xiàn)[43]用重構(gòu)函數(shù)沿網(wǎng)格單元各邊界的積分值計(jì)算BV。針對(duì)網(wǎng)格單元Ωi和與其相鄰的單元Ωij，沿單元邊界Γij的BV值由下式計(jì)算，

需要指出，BVD算法在同一網(wǎng)格單元上要對(duì)所有候補(bǔ)函數(shù)進(jìn)行重構(gòu)，并儲(chǔ)存相應(yīng)的BV或TBV值，會(huì)增加計(jì)算量和存儲(chǔ)量。然而，由于各網(wǎng)格單元間的計(jì)算過(guò)程相互獨(dú)立，BVD算法不會(huì)給并行處理帶來(lái)不利影響。

3 數(shù)值結(jié)果

我們利用標(biāo)量守恒律和歐拉方程的各種標(biāo)準(zhǔn)算例，對(duì)BVD算法進(jìn)行了系統(tǒng)的檢驗(yàn)，并與其他現(xiàn)有方法進(jìn)行了比較。

3.1 精度檢驗(yàn)

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

對(duì)于二維歐拉方程的密度擾動(dòng)傳播算例,初始條件設(shè)為:

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

表1給出各階PnTm-BVD格式的誤差及收斂率，均能達(dá)到相對(duì)應(yīng)的常系數(shù)多項(xiàng)式的階數(shù)。其誤差的絕對(duì)值與采用常系數(shù)多項(xiàng)式作為重構(gòu)函數(shù)的結(jié)果基本一致。

表1 PnTm-BVD格式的數(shù)值誤差及精度檢驗(yàn)

3.2 Jiang-Shu平流實(shí)驗(yàn)

作為檢驗(yàn)數(shù)值格式是否能夠同時(shí)捕捉連續(xù)和間斷解的算例，文獻(xiàn)[13]設(shè)計(jì)了一個(gè)既包括光滑分布又包括間斷解的平流初始條件。

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

圖2 Jiang-Shu平流標(biāo)準(zhǔn)算例在t=2000(1×106步)時(shí)的計(jì)算結(jié)果. 從左至右分別為五階WENO-JS[13]格式、WENO-Z[15]格式和P4T2-BVD格式[38]的計(jì)算結(jié)果，網(wǎng)格單元數(shù)為400.

3.3 一維歐拉方程的數(shù)值實(shí)驗(yàn)

對(duì)于一維歐拉方程，我們利用各類標(biāo)準(zhǔn)算例對(duì)BVD格式進(jìn)行了檢驗(yàn)，在抑制數(shù)值振蕩的同時(shí)，對(duì)間斷和光滑解均能得到高保真的計(jì)算結(jié)果。以下給出幾個(gè)具體算例的結(jié)果。

3.3.1 雙激波相互作用[64]

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

圖3 Two interacting blast waves標(biāo)準(zhǔn)算例在t=0.038時(shí)刻的P4T2-BVD和P4Tβv-BVD格式的計(jì)算結(jié)果(密度), 網(wǎng)格單元數(shù)為400.

圖4 Two interacting blast waves標(biāo)準(zhǔn)算例在t=0.038時(shí)刻的P8T3-BVD(左)和P10T3-BVD(右)格式的計(jì)算結(jié)果(密度),網(wǎng)格單元數(shù)為400.

3.3.2 激波與密度擾動(dòng)相互作用[65]

為了進(jìn)一步檢驗(yàn)BVD格式對(duì)高頻周期波的捕捉能力，我們計(jì)算了激波與密度擾動(dòng)相互作用的算例。圖5顯示不同階數(shù)BVD格式的計(jì)算結(jié)果。可以看到，不同階數(shù)的BVD格式都能抑制數(shù)值振蕩，并在光滑區(qū)正確地選擇多項(xiàng)式作為重構(gòu)函數(shù)。隨著多項(xiàng)式階數(shù)的提高，BVD能有效減小數(shù)值耗散，更精確地捕捉高頻波動(dòng)[40]。

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

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

3.3.3 Le Blanc 激波管問(wèn)題[66]

作為一個(gè)較困難的標(biāo)準(zhǔn)算例，Le Blanc激波管問(wèn)題的解包括非常強(qiáng)的激波和稀疏波，經(jīng)常用來(lái)檢驗(yàn)算法的魯棒性[66]。圖6給出了P4T2-BVD的計(jì)算結(jié)果(密度)。沒(méi)有出現(xiàn)明顯的數(shù)值振蕩。并能很好再現(xiàn)稀疏波、接觸間斷和激波[46]。數(shù)值實(shí)驗(yàn)表明，隨著網(wǎng)格分辨率的提高，右行激波能夠收斂到相應(yīng)精確解的位置。

圖6 Le Blanc一維激波管標(biāo)準(zhǔn)算例在t=6時(shí)刻的P4T2-BVD格式的計(jì)算結(jié)果(密度), 網(wǎng)格單元數(shù)為800.

3.3.4 剛性爆轟波[67]

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

圖7 剛性爆轟波標(biāo)準(zhǔn)算例在t=π/5時(shí)刻的P4T2-BVD格式的計(jì)算結(jié)果(密度),網(wǎng)格單元數(shù)為200.

3.4 二維歐拉方程的數(shù)值實(shí)驗(yàn)

3.4.1 雙馬赫反射[64]

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

3.4.2 爆轟波擾流問(wèn)題[68]

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

(a) P4T2-BVD

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

采用THINC函數(shù)作為候補(bǔ)重構(gòu)函數(shù)的BVD格式可以有效控制數(shù)值耗散，從而避免接觸間斷和物質(zhì)界面在計(jì)算過(guò)程中被平滑擴(kuò)散的現(xiàn)象。文獻(xiàn)[36]采用MUSCL-THINC-BVD格式對(duì)各類兩項(xiàng)可壓縮流動(dòng)進(jìn)行了模擬，取得了良好的計(jì)算結(jié)果。

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

(a) 密度

(b) 溫度

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

4 結(jié) 論

本文概述了BVD方法的基本思想，設(shè)計(jì)相關(guān)格式的基本思路，以及一些具有很強(qiáng)實(shí)用價(jià)值的BVD格式。并通過(guò)單相和兩相可壓縮流動(dòng)的一些典型算例驗(yàn)證了BVD格式的特點(diǎn)和優(yōu)勢(shì)。

BVD方法把減小數(shù)值耗散作為設(shè)計(jì)空間重構(gòu)算法的指導(dǎo)思想，具有明確的物理含義和普遍性，是一條設(shè)計(jì)和改進(jìn)數(shù)值格式的新途徑。針對(duì)研究對(duì)象，選用合適的BVD許容函數(shù)作為空間重構(gòu)的候補(bǔ)選項(xiàng)，并設(shè)計(jì)相應(yīng)的BVD算法可以構(gòu)建一類全新的數(shù)值格式。這些格式與基于多項(xiàng)式重構(gòu)和非線性限制修正的傳統(tǒng)方法有本質(zhì)區(qū)別。數(shù)值實(shí)驗(yàn)結(jié)果表明，BVD格式可以得到明顯優(yōu)于已有算法的計(jì)算結(jié)果，能同時(shí)對(duì)各種尺度的連續(xù)和間斷結(jié)構(gòu)都獲得高保真的數(shù)值解。

采用THINC或類似的Sigmoid函數(shù)作為候補(bǔ)重構(gòu)函數(shù)可以明顯改善基于歐拉網(wǎng)格的數(shù)值方法對(duì)接觸間斷和物質(zhì)界面在計(jì)算過(guò)程中被平滑擴(kuò)散的問(wèn)題，對(duì)解決含有化學(xué)反應(yīng)面或移動(dòng)物質(zhì)界面流動(dòng)問(wèn)題的實(shí)際應(yīng)用具有重要意義。

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

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