孔令發(fā)，劉 偉，董義道

(國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073)

0 引 言

相比結(jié)構(gòu)化網(wǎng)格和笛卡爾網(wǎng)格，非結(jié)構(gòu)網(wǎng)格因其對復(fù)雜外形的適應(yīng)性、網(wǎng)格生成的自動(dòng)化程度[1-3]以及對運(yùn)動(dòng)邊界、多體分離、流場自適應(yīng)等問題的模擬能力，逐漸成為學(xué)術(shù)研究和應(yīng)用的重點(diǎn)[4]。

基于非結(jié)構(gòu)網(wǎng)格的數(shù)值格式通常包含有限元（Finite Element, FE）[5]、間斷伽遼金（Discontinuous Galerkin, DG）[6-7]和有限體積（Finite Volume, FV）[8-13]等方法，其中，具有二階精度的FV方法因其在魯棒性和計(jì)算開銷較小等方面的優(yōu)勢而被廣泛應(yīng)用于當(dāng)前主流的CFD商業(yè)與開源軟件。這也一定程度上極大促進(jìn)了二階FV格式的發(fā)展，使其趨于成熟。

在二階精度的基礎(chǔ)上，近年來具有高階精度的FV離散格式同樣得到了較快發(fā)展，以滿足高效計(jì)算與高保真數(shù)值模擬的需求[12-15]。Ollivier-Gooch與Jalali等[16-18]系統(tǒng)研究了高階非結(jié)構(gòu)FV格式中的k-exact重構(gòu)[19,20]、曲邊界處理與通量、源項(xiàng)積分等一系列問題。Li[21]在高階非結(jié)構(gòu)數(shù)值格式的框架內(nèi)探索了新的限制器構(gòu)造策略，以及考慮黏性邊界條件的約束k-exact重構(gòu)。Kong等基于k-exact方法，在二階與高階FV離散格式的框架內(nèi)提出了一種新的重構(gòu)模板[22]與單元參考點(diǎn)選擇策略[23]，有效改善了原有格式在各向異性網(wǎng)格上的數(shù)值表現(xiàn)。

相比重構(gòu)與離散算法研究，在非結(jié)構(gòu)FV算法，尤其是高階精度算法框架內(nèi)，針對邊界條件的細(xì)致研究相對較少，然而實(shí)踐表明，采用不同的邊界條件設(shè)置方法可導(dǎo)致流場求解效果的差異[24]。如果邊界條件設(shè)置恰當(dāng)，求解過程的魯棒性以及計(jì)算效果能夠得到明顯改善；反之，若邊界條件設(shè)置不當(dāng)，則有可能收斂至錯(cuò)誤的流場，甚至導(dǎo)致計(jì)算發(fā)散。

對于真實(shí)流動(dòng)，通?？蓪⑦吔鐥l件的施加分為強(qiáng)施加與弱施加兩類。強(qiáng)施加是一種變量施加方法，即通過邊界條件得到邊界位置的變量；而弱施加是一種通量施加方法，即通過邊界條件得到邊界位置的數(shù)值通量。強(qiáng)施加方法一般基于特征理論，認(rèn)為雙曲守恒律方程代表特征波的傳播，在對應(yīng)的邊界位置存在傳入或者傳出計(jì)算域的特征波[25-26]。其中，外傳波的特征一般取決于內(nèi)部流場；而內(nèi)傳波需要基于無黏關(guān)系式，在界面構(gòu)造演化方程[27]。

對基于非結(jié)構(gòu)網(wǎng)格的數(shù)值格式，當(dāng)前較為常用的強(qiáng)施加邊界條件是引入虛擬單元來設(shè)置特征邊界條件[28-29]，該方法能夠準(zhǔn)確反映邊界附近特征波傳播情況，但此過程需要人為判斷特征波傳播方向以得到時(shí)變幅值，對于多維效應(yīng)明顯的流動(dòng)，通常還需要引入橫向修正[29]。此外，對邊界面處特征虛擬單元的變量賦值過程引入了一階誤差，嚴(yán)重破壞了流場求解精度[30-31]。而本文的研究目標(biāo)是探討邊界條件在高階精度有限體積方法中的表現(xiàn)效果，因此不再考慮這種強(qiáng)施加邊界條件在高階精度求解器上的應(yīng)用。

不同于強(qiáng)施加方法，弱施加邊界條件旨在得到邊界面的數(shù)值通量。對于內(nèi)部單元面，可根據(jù)單元面處的變量左右值采用近似黎曼解構(gòu)造數(shù)值通量，而對于邊界面的通量構(gòu)造同樣可采用此方法。Kopriva等[32-33]指出，F(xiàn)UN3D的邊界條件處理也采用了這種方法。但當(dāng)前對亞聲速入口、出口以及遠(yuǎn)場這類邊界條件的主要處理方式是基于一維黎曼不變量來得到邊界面處的右狀態(tài)矢量，這種方法仍需采用特征關(guān)系式，并基于黎曼不變量進(jìn)行繁冗推導(dǎo)。為簡化這類弱施加邊界條件的處理方法，Atkins和Casper[34]針對結(jié)構(gòu)網(wǎng)格上的高精度格式，提出了一種基于有限波模型的精度保持邊界條件施加方法來統(tǒng)一處理包括亞聲速入口、出口以及遠(yuǎn)場等邊界，該方法基于有限波模型得到的精確解，將內(nèi)部流場與外部瞬態(tài)條件有效聯(lián)系起來，得到了高階收斂精度。Dong等[30-31]將這種邊界條件由基于結(jié)構(gòu)網(wǎng)格的有限差分方法推廣至二階精度非結(jié)構(gòu)有限體積方法，并將其與常用的特征虛擬單元邊界條件對比后發(fā)現(xiàn)，采用黎曼邊界條件能夠得到較高的計(jì)算精度與收斂效率，同時(shí)也可顯著提高求解器的計(jì)算穩(wěn)定性。

在此基礎(chǔ)上，本文將這種弱施加邊界條件推廣至高階精度非結(jié)構(gòu)有限體積方法，以檢驗(yàn)其在高階精度求解器中的表現(xiàn)效果。此外，在高階有限體積方法的框架內(nèi)將這種邊界條件推廣至基于制造解方法（MMS）[35-38]的流動(dòng)。這種特殊的數(shù)值流動(dòng)當(dāng)前被廣泛用來檢驗(yàn)算法的有效性，但通常在基于MMS流動(dòng)研究某些問題時(shí)，一般未考慮對其施加邊界條件。較常見的處理方式是所有邊界面處的變量左值由內(nèi)部流場插值得到，而右值均賦值為當(dāng)前位置的精確解。Folkner和Katz[35]以及Wang等[36]在二階精度非結(jié)構(gòu)有限體積方法的框架內(nèi)驗(yàn)證了施加邊界條件的MMS流動(dòng)的數(shù)值表現(xiàn)，經(jīng)檢驗(yàn)，對MMS流動(dòng)施加邊界條件后可達(dá)到二階精度。但在高精度數(shù)值方法中，對此類流動(dòng)施加邊界條件，尤其是黎曼邊界條件的相關(guān)研究還未開展。因此，本文計(jì)劃在高階精度非結(jié)構(gòu)有限體積方法的框架內(nèi)，驗(yàn)證對MMS流動(dòng)和部分真實(shí)流動(dòng)施加黎曼邊界條件的合理性，并對比其與常用無反射邊界條件[35]間的差異，以進(jìn)一步推廣此類弱施加邊界條件在非結(jié)構(gòu)網(wǎng)格數(shù)值格式中的應(yīng)用。

1 非結(jié)構(gòu)有限體積格式

1.1 控制方程

對無黏流動(dòng)探討弱施加黎曼邊界條件，需要求解的控制方程為：

其中，u為待求解變量，F(xiàn)c為無黏通量。該方程可在散度定理的基礎(chǔ)上，通過FV格式離散為：

式中Nf與NG分 別為網(wǎng)格面與通量點(diǎn)的數(shù)量， Φjk為通量點(diǎn)處的數(shù)值通量函數(shù)，本文采用Roe格式[39]來計(jì)算該函數(shù)。 ωk為對應(yīng)通量點(diǎn)的權(quán)系數(shù)[23]，njSj為對應(yīng)單元面的面積矢量。

對于高階方法，一個(gè)通量點(diǎn)無法滿足積分精度要求，因此需要多個(gè)通量點(diǎn)（如圖1所示）。對應(yīng)不同精度下的通量點(diǎn)的數(shù)目可參考文獻(xiàn)[23]。在基于Roe格式構(gòu)造通量點(diǎn)數(shù)值通量時(shí)，需要知道該點(diǎn)的變量左右值，并且變量左右值的獲取依賴于從單元參考點(diǎn)到面通量點(diǎn)的泰勒展開，因此需重構(gòu)控制體單元的變量梯度及其高階導(dǎo)數(shù)。

圖1 數(shù)值通量構(gòu)造示意圖Fig. 1 Schematic diagram of the numerical flux evaluation

1.2 k-exact重構(gòu)算法

對于二階格式，常采用最小二乘（Least Square, LSQ）法與格林高斯（Green Gauss, GG）方法來獲得變量梯度。但對基于積分型控制方程的高階FV格式而言，傳統(tǒng)的LSQ或GG方法均無法滿足高于二階的重構(gòu)精度要求，因此本文采用k-exact方法[16-18]計(jì)算變量的梯度與高階導(dǎo)數(shù)。該方法的詳細(xì)過程本文不再重述，最終需要求解的方程組形式如下：

同時(shí)，權(quán)系數(shù) ωij為參考點(diǎn)間距離的倒數(shù)，以強(qiáng)調(diào)與中心控制體相鄰單元對重構(gòu)過程的影響。

2 弱施加邊界條件

由于強(qiáng)施加方法多依賴于特征虛擬單元來處理相應(yīng)的邊界條件，該方法在對無黏邊界處的虛擬單元賦值過程引入了一階誤差。因此，本文不再討論這種邊界條件強(qiáng)施加方式。

為了在邊界面施加弱邊界條件，根據(jù)近似黎曼公式，需要給定邊界面處的兩個(gè)狀態(tài)，如圖2所示，其中左狀態(tài)由內(nèi)部流場插值得到，右狀態(tài)依賴于相應(yīng)的邊界條件。

圖2 重構(gòu)模板與邊界條件弱施加示意圖Fig. 2 Schematic diagram for the reconstruction stencil and weak-imposition boundary condition

在基于內(nèi)部流場插值時(shí)需要注意避免引入一階誤差。本文研究具有三階精度的有限體積方法，因此首先可通過內(nèi)部流場信息得到邊界面左值，即：

式中WC為邊界單元參考點(diǎn)處的變量值，b為邊界面通量點(diǎn)。因此，通過方程（5）可得到邊界面高斯積分點(diǎn)處的左狀態(tài)矢量。邊界面右值需要根據(jù)相應(yīng)的邊界條件獲得。下面將針對亞聲速入口、出口以及遠(yuǎn)場等不同邊界，給出其邊界條件的具體施加方式。

2.1 無黏壁面邊界

對于無黏壁面，不同弱施加方法的處理過程基本一致。在施加邊界條件時(shí)，需要保證速度的左右狀態(tài)關(guān)于壁面邊界對稱。為滿足這個(gè)條件，右狀態(tài)矢量可定義為：

其中逆變速度由左值計(jì)算得到：

壓力和溫度滿足零梯度條件，即：

2.2 超聲速入口與出口邊界

由于在超聲速條件下特征波傳播方向唯一，因此超聲速入口、出口的處理相對簡單。其中超聲速入口右值直接取來流值，即：

超聲速出口右值與左值相等：

2.3 亞聲速入口與出口邊界

與壁面、超聲速入口和出口的處理方式不同，在處理亞聲速入口、出口以及遠(yuǎn)場時(shí)，本文采用黎曼邊界條件。該方法通過求解對應(yīng)的黎曼問題，在邊界面位置可得到物理上相容的左右狀態(tài)，因此避免了采用特征關(guān)系式與黎曼不變量對這類邊界條件的判斷與推導(dǎo)[35]。

如圖3所示，有限波模型假設(shè)內(nèi)部流場與外部瞬態(tài)狀態(tài)通過兩道由接觸間斷分隔的膨脹波進(jìn)行過渡，圖中標(biāo)有u的虛線代表接觸間斷，另外兩條紅色直線對應(yīng)膨脹波，膨脹波之間的區(qū)域?qū)⑵浞Q為星區(qū)。

現(xiàn)構(gòu)造如下黎曼問題，無窮遠(yuǎn)處的自由來流參數(shù)作為右狀態(tài)矢量，記為WR，而左狀態(tài)矢量WL由內(nèi)部流場插值得到。通過精確求解此黎曼問題，可以獲得星區(qū)的左右狀態(tài)矢量，分別記為與。 星區(qū)的壓力可通過如下關(guān)系式獲得：

其中，c為局部聲速，γ為比熱比。逆變速度的定義為：un=u·n，其中u為笛卡爾坐標(biāo)系下的速度矢量，n為邊界面處的單位法向矢量。在得到星區(qū)壓力的基礎(chǔ)上，可采用等熵關(guān)系式得到星區(qū)的密度：

以及星區(qū)的逆變速度：

在此基礎(chǔ)上，邊界右狀態(tài)可通過星區(qū)狀態(tài)矢量即給定。

圖3 有限波模型示意圖Fig. 3 Schematic diagram of the finite wave model

相比常用的無反射邊界條件[35]，采用黎曼邊界條件處理亞聲速入口、出口與遠(yuǎn)場邊界時(shí)，均可通過精確求解對應(yīng)的黎曼問題來給定邊界右值，無需再基于法向速度與特征關(guān)系來判斷不同邊界，并單獨(dú)處理。因此，這種邊界條件施加方式可有效簡化對亞聲速入口、出口以及遠(yuǎn)場邊界的處理過程。

3 MMS流動(dòng)中的精度驗(yàn)證

當(dāng)前，MMS方法廣泛應(yīng)用于對算法計(jì)算精度的檢驗(yàn)，但在基于該方法驗(yàn)證精度時(shí)，一般未對其施加任何物理邊界條件，即邊界面的右狀態(tài)均以參考點(diǎn)處的精確解賦值。本節(jié)將對MMS流動(dòng)采用黎曼方法處理亞聲速入口、出口邊界，以驗(yàn)證對MMS流動(dòng)施加黎曼邊界條件的合理性。同時(shí)為了驗(yàn)證工作的完整性，還測試了對MMS流動(dòng)施加無黏壁面以及超聲速入口、出口這類具有真實(shí)物理意義邊界條件后的數(shù)值表現(xiàn)，以期在高階精度非結(jié)構(gòu)有限體積框架內(nèi)，為MMS流動(dòng)提供一種便捷合理的邊界處理方式。

測試網(wǎng)格分別采用了六套由疏到密的規(guī)則三角形網(wǎng)格與擾動(dòng)三角形網(wǎng)格。網(wǎng)格示意圖與網(wǎng)格點(diǎn)分布情況見圖4與表1。

圖4 規(guī)則與擾動(dòng)三角形網(wǎng)格示意圖Fig. 4 Schematics of regular and randomly perturbed grids

表1 背景四邊形網(wǎng)格在x與y方向的分布量Table 1 Distribution of background quadrilateral grid in x- and y-directions

當(dāng)來流速度為亞聲速時(shí)，計(jì)算域的左右邊界分別設(shè)置為亞聲速入口與出口；對于超聲速流動(dòng)，左右邊界分別對應(yīng)超聲速入口和出口。此外，設(shè)置壁面邊界時(shí)，考慮到特征波的傳播方向，統(tǒng)一將計(jì)算域上邊界設(shè)置為壁面，示意圖見圖5。

圖5 邊界條件設(shè)置示意圖Fig. 5 Schematic diagram of the corresponding boundary conditions

3.1 無邊界條件

在研究不同邊界條件的影響之前，為了驗(yàn)證求解器的求解精度，首先測試了不添加邊界條件的MMS流動(dòng)計(jì)算效果。在相應(yīng)的邊界高斯積分點(diǎn)處，左值均由內(nèi)部流場插值得到，而右值由MMS制造解給定。其形式如下：

流動(dòng)控制方程如式（1），因此可首先根據(jù)控制方程推導(dǎo)出源項(xiàng)s的統(tǒng)一表達(dá)式為：

基于該解析解，得到的流場如圖6所示。

計(jì)算收斂后，統(tǒng)計(jì)了每套網(wǎng)格誤差的L1、L2以及L∞范數(shù)來充分測試高階求解器的計(jì)算精度。圖7中Slope3表示三階精度標(biāo)準(zhǔn)線。從圖7中可以看出，不添加邊界條件時(shí)，在疏網(wǎng)格上精度一開始并未達(dá)到三階，但隨著網(wǎng)格加密，計(jì)算結(jié)果逐漸趨于設(shè)計(jì)精度。因此，在不施加邊界條件的情況下，首先驗(yàn)證了求解器具有三階精度。接下來對計(jì)算域的不同位置單獨(dú)施加黎曼邊界條件，來驗(yàn)證在高階方法中對此類流動(dòng)施加邊界條件的合理性。

圖6 當(dāng)前制造解流場Fig. 6 Flow fields of current manufactured solutions

圖7 無邊界條件的誤差曲線Fig. 7 Error curves without any boundary condition

3.2 亞聲速入口、出口邊界

對于真實(shí)流動(dòng)的亞聲速入口與出口，可通過精確求解黎曼問題來獲得對應(yīng)的邊界條件。在求解黎曼問題時(shí)，WR通常用無窮遠(yuǎn)處的來流參數(shù)賦值。而MMS流動(dòng)與真實(shí)流動(dòng)不同，一般不包含無窮遠(yuǎn)處的參數(shù)信息。因此嘗試在MMS流動(dòng)中，WR統(tǒng)一用邊界處的制造解提供，并在得到右狀態(tài)的基礎(chǔ)上，通過求解黎曼問題來得到星區(qū)狀態(tài)矢量最終由該矢量來獲得亞聲速入口、出口處的右狀態(tài)參數(shù)。

驗(yàn)證這兩種邊界條件時(shí)，需要保證在對應(yīng)的入口與出口位置速度為亞聲速，因此給定解的形式如下：

對應(yīng)該制造解，得到的流場如圖8所示。

圖8 當(dāng)前制造解流場Fig. 8 Flow fields of current manufactured solutions

對入口與出口條件單獨(dú)驗(yàn)證，在驗(yàn)證其中一個(gè)邊界條件時(shí)，其余邊界右值均使用MMS精確解賦值，這樣有利于研究不同邊界對流動(dòng)求解精度的直接影響。統(tǒng)計(jì)誤差如圖9與圖10所示。

從圖9與圖10反映的結(jié)果可以看出，對MMS流動(dòng)單獨(dú)施加亞聲速入口、出口邊界條件時(shí)，不論是在規(guī)則網(wǎng)格還是擾動(dòng)網(wǎng)格上，均達(dá)到了三階精度。因此，對MMS流動(dòng)的亞聲速入口與出口施加黎曼邊界條件不會(huì)影響流動(dòng)的求解精度。

圖9 施加亞聲速入口邊界條件的誤差曲線Fig. 9 Error curves with the subsonic inlet boundary condition

圖10 施加亞聲速出口邊界條件的誤差曲線Fig. 10 Error curves with the subsonic outlet boundary condition

3.3 無黏壁面邊界

不同的弱施加邊界條件對無黏壁面的處理方式基本一致。壁面法向速度為0，壓力和溫度的右值與左值相等。為保證速度的法向分量為0，驗(yàn)證壁面時(shí)給定解的形式為：

Folkner和Katz等[35]在二階精度有限體積方法中驗(yàn)證了對MMS流動(dòng)施加此類含有真實(shí)物理意義邊界條件的可行性。本節(jié)將其推廣至高階精度非結(jié)構(gòu)有限體積方法，以檢驗(yàn)在高階方法中對MMS流動(dòng)施加此類邊界的合理性?；诜匠蹋?7）得到的流場如圖11所示。

圖11 MMS流場信息Fig. 11 Flow fields of the MMS flow

在驗(yàn)證壁面邊界時(shí)，除壁面外其余邊界右值均由MMS精確解賦值。在規(guī)則與擾動(dòng)網(wǎng)格下的統(tǒng)計(jì)誤差如圖12所示。由圖12可以看出，對MMS流動(dòng)施加相應(yīng)的壁面邊界條件后，不論在規(guī)則網(wǎng)格還是擾動(dòng)網(wǎng)格上，計(jì)算結(jié)果仍可達(dá)到三階精度。由此驗(yàn)證了對MMS流動(dòng)施加壁面邊界條件的可行性，同時(shí)也將文獻(xiàn)[35]中提出的壁面邊界驗(yàn)證方式有效推廣至高階精度非結(jié)構(gòu)有限體積算法。

圖12 施加無黏壁面邊界條件的誤差曲線Fig. 12 Error curves with the inviscid wall boundary condition

3.4 超聲速入口、出口邊界

相比亞聲速入口、出口的邊界條件施加方式，超聲速入口與出口的邊界條件施加相對簡便—入口采用MMS精確解賦值，出口的右值等于左值。為保證入口、出口處為超聲速，給定解的形式為：

圖13為該制造解對應(yīng)的流場云圖。與3.2和3.3節(jié)的驗(yàn)證方式相同，對入口與出口條件單獨(dú)驗(yàn)證。圖14與15分別給出了施加超聲速入口與出口條件得到的誤差。從圖14與圖15反映的結(jié)果來看，單獨(dú)對MMS流動(dòng)施加超聲速入口或出口條件，并未破壞計(jì)算精度，L1、L2與L∞誤差均達(dá)到了三階精度。我們將三類邊界條件在規(guī)則密網(wǎng)格上得到的計(jì)算精度匯總，結(jié)果如表2所示。

從表2可以看出，對MMS流動(dòng)的不同邊界單獨(dú)施加相應(yīng)的邊界條件后，在密網(wǎng)格上L1與L2誤差均達(dá)到了三階精度，L∞誤差略低于但接近三階精度。此外，從誤差統(tǒng)計(jì)曲線可以明顯看出，對MMS流動(dòng)施加合理的邊界條件后，隨著網(wǎng)格加密，求解精度的變化較小，即使在稀疏網(wǎng)格上結(jié)果同樣接近三階精度，并趨于一致收斂。因此，這也從一個(gè)層面論證了在高階方法中對MMS流動(dòng)施加具有真實(shí)物理意義邊界條件的合理性。

綜上，黎曼邊界條件在MMS流動(dòng)中具有較好的數(shù)值表現(xiàn)，尤其在處理亞聲速入口、出口時(shí)，為MMS流動(dòng)的高精度數(shù)值模擬提供了一種簡單有效的邊界條件施加方式。

圖13 MMS流場信息Fig. 13 Flow fields of the MMS flow

圖14 施加超聲速入口邊界條件的誤差曲線Fig. 14 Error curves with the supersonic inlet boundary condition

圖15 施加超聲速出口邊界條件的誤差曲線Fig. 15 Error curves with the supersonic outlet boundary condition

表2 不同弱施加邊界條件在密網(wǎng)格之間的計(jì)算精度Table 2 Accuracy of different weak boundary conditions between vfin and vvfin grids

4 亞聲速無黏圓柱繞流驗(yàn)證

為了測試上述邊界條件在真實(shí)流動(dòng)中的表現(xiàn)，在亞聲速無黏圓柱繞流上檢驗(yàn)其數(shù)值效果。該流動(dòng)是典型的亞聲速算例，邊界條件包括亞聲速遠(yuǎn)場和無黏壁面。由之前的分析可以知，黎曼邊界可有效簡化對亞聲速遠(yuǎn)場的處理過程。為對比黎曼邊界條件在高階求解器中的數(shù)值表現(xiàn)，在處理遠(yuǎn)場邊界時(shí)，還測試了常用的基于一維黎曼不變量的無反射邊界條件[35]。網(wǎng)格采用O型拓?fù)?，壁面附近的局部放大圖如圖16所示。主要來流參數(shù)Ma=0.38，遠(yuǎn)場取20倍圓柱直徑。網(wǎng)格生成可參照文獻(xiàn)[38]。計(jì)算在三套網(wǎng)格上進(jìn)行，密網(wǎng)格徑向與周向分別分布32與128個(gè)網(wǎng)格點(diǎn)。網(wǎng)格徑向分布情況為：

其中α=1.1， 圓柱半徑r0=0.5。對最密網(wǎng)格疏化可得到剩余兩套網(wǎng)格。

圖16 圓柱繞流計(jì)算網(wǎng)格示意圖Fig. 16 Schematic diagram of grids for the cylinder flow simulation

基于兩種不同的弱施加邊界條件，首先測試了由兩種邊界條件得到的流場壓力云圖，如圖17所示。從流場云圖可以看出，由兩種邊界條件得到的流場均具有較好的對稱性。

圖17 不同邊界條件下的圓柱繞流流場Fig. 17 Flow fields of the cylinder flow with different boundary conditions

為進(jìn)一步量化對比兩種邊界條件之間的差異，統(tǒng)計(jì)了沿圓柱表面的總壓恢復(fù)系數(shù)與熵誤差。總壓恢復(fù)系數(shù)的定義如下：

根據(jù)等熵關(guān)系式可得到總壓pt和靜壓p之間的關(guān)系，

由于該算例為無黏亞聲速流動(dòng)，因此理想情況下總壓恢復(fù)系數(shù)應(yīng)當(dāng)精確等于1.0。圖18給出了總壓恢復(fù)系數(shù)沿壁面的分布情況。從圖中可明顯看出，在梳網(wǎng)格上，總壓損失較大，而隨著網(wǎng)格的不斷加密，數(shù)值耗散逐漸降低，總壓恢復(fù)系數(shù)接近理想值。

圖18 圓柱繞流總壓恢復(fù)系數(shù)分布曲線Fig. 18 Total pressure recovery coefficients of the cylinder flow

此外，統(tǒng)計(jì)了由不同邊界條件得到的熵誤差，結(jié)果如圖19所示。從圖19可以看出，采用兩種遠(yuǎn)場邊界條件得到的計(jì)算精度基本一致，但由黎曼邊界得到的誤差值略低于常用的無反射邊界，更利于計(jì)算結(jié)果的準(zhǔn)確性。綜上，在亞聲速圓柱繞流中，黎曼邊界條件的數(shù)值表現(xiàn)優(yōu)于無反射邊界條件。

圖19 熵誤差統(tǒng)計(jì)Fig. 19 Entropy errors

5 Gaussian Pulse非定常熵波驗(yàn)證

本節(jié)引入添加初始高斯脈沖（Gaussian Pulse）的非定常熵波傳播算例，來測試?yán)杪吔鐥l件在出口處的無反射特性。初始擾動(dòng)函數(shù)如下[40]：

式中xc取0.5，計(jì)算域?yàn)閤×y∈[0,1]×[0,1]。并且除密度外，其他變量的初始化[31]如下：

計(jì)算過程中，分別對比了由黎曼邊界條件與無反射邊界條件得到的出口特性。在施加相應(yīng)的邊界條件時(shí)，為防止y方向邊界條件對計(jì)算結(jié)果的影響，計(jì)算域的上下邊界均采用外推邊界。具體的邊界條件施加示意圖如圖20。從圖20可以看出，黎曼邊界條件針對亞聲速入口與出口可通過精確求解對應(yīng)的黎曼問題來統(tǒng)一處理，相比無反射邊界條件更加簡便。

圖20 兩種邊界條件的施加示意圖Fig. 20 Schematic diagram for the imposition of two boundary conditions

為對比兩種邊界在出口處的無反射特性，分別設(shè)置了四套由疏到密的均勻四邊形網(wǎng)格，分布量從60×60至160×160（圖21）。在四套網(wǎng)格下，統(tǒng)計(jì)了添加初始高斯擾動(dòng)的ρ沿x方向的變化情況，同時(shí)參考文獻(xiàn)[40]，分別給出了時(shí)間t= 0.9 s、1.1 s以及1.5 s時(shí)的結(jié)果。統(tǒng)計(jì)結(jié)果如圖22所示。表3、表4為不同時(shí)刻兩種邊界條件得到的L2誤差的具體數(shù)值。

圖21 疏網(wǎng)格與密網(wǎng)格Fig. 21 Coarse and fine grids

圖22 不同網(wǎng)格上密度在三個(gè)時(shí)刻的統(tǒng)計(jì)結(jié)果Fig. 22 Density on different grids at three moments

表3 無反射邊界條件在不同時(shí)刻的L2誤差統(tǒng)計(jì)Table 3 L2 errors of the non-reflective boundary condition at different moments

表4 黎曼邊界條件在不同時(shí)刻的L2誤差統(tǒng)計(jì)Table 4 L2 errors of the Riemann boundary condition at different moments

從圖22中可以看出，隨著網(wǎng)格的不斷加密，兩種邊界條件的密度變化曲線接近，在出口邊界均具有較好的無反射特性。因此說明了對亞聲速非定常擾流施加黎曼邊界條件具有一定的合理性。此外，結(jié)合表3的數(shù)據(jù)可以看出，黎曼邊界條件與無反射邊界條件在四套網(wǎng)格上，不同時(shí)刻的L2誤差接近，并在t=0.9 s與t=1.5 s時(shí)誤差均低于無反射邊界條件。

最后，相比無反射邊界條件，黎曼邊界條件針對亞聲速入口與出口，可通過精確求解黎曼問題來統(tǒng)一處理，有效簡化了邊界條件的處理過程。

6 結(jié) 論

本文對黎曼邊界條件展開了細(xì)致討論，并將其推廣至高精度非結(jié)構(gòu)有限體積方法。現(xiàn)將這種弱施加邊界條件在高階精度數(shù)值方法中的表現(xiàn)情況總結(jié)如下：

1）從方法實(shí)現(xiàn)層面而言，相比傳統(tǒng)無反射邊界條件，黎曼邊界在處理亞聲邊界時(shí)，均可通過精確求解黎曼問題來統(tǒng)一施加，無需判斷特征方向。

2）從數(shù)值結(jié)果來看，在高階精度框架內(nèi)對MMS流動(dòng)單獨(dú)施加黎曼邊界條件，以及壁面、超聲速入口與出口邊界條件后，計(jì)算結(jié)果均可達(dá)到設(shè)計(jì)精度，并趨于一致收斂。

3）在真實(shí)流動(dòng)中，基于黎曼邊界條件和無反射邊界條件得到的流場及總壓恢復(fù)系數(shù)效果接近，并且其熵誤差更低。

4）在添加初始高斯脈沖擾動(dòng)的非定常流動(dòng)中，新的邊界條件在出口處維持了較好的無反射特性。因此，在非定常擾流中施加黎曼邊界條件合理可行。

綜上，黎曼邊界條件在高階精度非結(jié)構(gòu)有限體積方法中的有效性得到了較好地驗(yàn)證，為非結(jié)構(gòu)高精度數(shù)值格式提供了一種便捷有效的亞聲速邊界處理途徑。在接下來的工作中，將進(jìn)一步驗(yàn)證推廣黎曼邊界在黏性流動(dòng)和實(shí)際工程問題中的應(yīng)用。