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

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

0 引 言

相比結構化網格和笛卡爾網格，非結構網格因其對復雜外形的適應性、網格生成的自動化程度[1-3]以及對運動邊界、多體分離、流場自適應等問題的模擬能力，逐漸成為學術研究和應用的重點[4]。

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

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

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

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

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

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

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

1 非結構有限體積格式

1.1 控制方程

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

其中，u為待求解變量，Fc為無黏通量。該方程可在散度定理的基礎上，通過FV格式離散為：

式中Nf與NG分 別為網格面與通量點的數量， Φjk為通量點處的數值通量函數，本文采用Roe格式[39]來計算該函數。 ωk為對應通量點的權系數[23]，njSj為對應單元面的面積矢量。

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

圖1 數值通量構造示意圖Fig. 1 Schematic diagram of the numerical flux evaluation

1.2 k-exact重構算法

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

同時，權系數 ωij為參考點間距離的倒數，以強調與中心控制體相鄰單元對重構過程的影響。

2 弱施加邊界條件

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

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

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

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

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

2.1 無黏壁面邊界

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

其中逆變速度由左值計算得到：

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

2.2 超聲速入口與出口邊界

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

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

2.3 亞聲速入口與出口邊界

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

如圖3所示，有限波模型假設內部流場與外部瞬態(tài)狀態(tài)通過兩道由接觸間斷分隔的膨脹波進行過渡，圖中標有u的虛線代表接觸間斷，另外兩條紅色直線對應膨脹波，膨脹波之間的區(qū)域將其稱為星區(qū)。

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

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

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

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

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

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

3 MMS流動中的精度驗證

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

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

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

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

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

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

3.1 無邊界條件

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

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

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

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

圖6 當前制造解流場Fig. 6 Flow fields of current manufactured solutions

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

3.2 亞聲速入口、出口邊界

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

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

對應該制造解，得到的流場如圖8所示。

圖8 當前制造解流場Fig. 8 Flow fields of current manufactured solutions

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

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

圖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，驗證壁面時給定解的形式為：

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

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

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

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

3.4 超聲速入口、出口邊界

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

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

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

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

圖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 不同弱施加邊界條件在密網格之間的計算精度Table 2 Accuracy of different weak boundary conditions between vfin and vvfin grids

4 亞聲速無黏圓柱繞流驗證

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

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

圖16 圓柱繞流計算網格示意圖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

為進一步量化對比兩種邊界條件之間的差異，統(tǒng)計了沿圓柱表面的總壓恢復系數與熵誤差?？倝夯謴拖禂档亩x如下：

根據等熵關系式可得到總壓pt和靜壓p之間的關系，

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

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

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

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

5 Gaussian Pulse非定常熵波驗證

本節(jié)引入添加初始高斯脈沖（Gaussian Pulse）的非定常熵波傳播算例，來測試黎曼邊界條件在出口處的無反射特性。初始擾動函數如下[40]：

式中xc取0.5，計算域為x×y∈[0,1]×[0,1]。并且除密度外，其他變量的初始化[31]如下：

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

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

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

圖21 疏網格與密網格Fig. 21 Coarse and fine grids

圖22 不同網格上密度在三個時刻的統(tǒng)計結果Fig. 22 Density on different grids at three moments

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

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

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

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

6 結 論

本文對黎曼邊界條件展開了細致討論，并將其推廣至高精度非結構有限體積方法?，F將這種弱施加邊界條件在高階精度數值方法中的表現情況總結如下：

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

2）從數值結果來看，在高階精度框架內對MMS流動單獨施加黎曼邊界條件，以及壁面、超聲速入口與出口邊界條件后，計算結果均可達到設計精度，并趨于一致收斂。

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

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

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