基于離散等價方程的滑移壁面邊界約束算法

魏雁昕，劉 君

(大連理工大學 航空航天學院，大連 116024)

0 引 言

計算流體力學（Computational Fluid Dynamics, CFD）的實質是求解偏微分方程的初邊值問題。盡管實際物理流動千變萬化，但是在很多問題的計算中原始方程是相同的，數值模擬得到的流場差異源自空間形狀或初場不同。CFD 應用于定常流動情況時，收斂的數值解與初始條件無關，在流動方程已經確定的情況下，決定流場特性的主要因素是邊界的求解方法。從CFD 發展的早期開始，許多學者針對有限差分法的邊界條件進行了廣泛的研究[1-4]。針對外流場模擬的空氣動力學問題，兩類邊界條件起著至關重要的作用：一類稱之為無反射邊界條件，用于計算亞聲速遠場邊界；另一類則為無穿透邊界條件，用于處理壁面邊界。對于遠場邊界，通常使用特征變量方法得到遠場特征邊界條件[5]或者基于特征線相容關系得到遠場Riemann 邊界條件[6]。遠場邊界距離物體較遠時（如對于翼型計算，最小距離通常取10 倍弦長），使用上述兩種方法求解亞聲速遠場邊界均可得到滿意的計算結果，因此不作為本文的研究內容。而物面邊界的計算決定了近壁面附近的流場特性，直接影響物體氣動力/熱預測的準確性。

本文針對有限差分法的滑移壁面邊界條件開展研究，以二維歐拉方程為例，在壁面處存在4 個未知的流動變量，只有法向無穿透條件是準確已知的，即vn=0，其余3 個物理量無法直接給出。CFD 早期的邊界處理是基于雙曲方程的特征理論，將邊界約束帶入原始方程，把邊界條件的影響作為源項處理，推導出針對此類邊界的修正方程，通過求解修正后的控制方程得到邊界狀態[3]。因為修正后的控制方程難以寫成守恒形式，邊界節點和內點的離散格式不相容，降低了計算的穩定性，也不利于軟件的工程化應用，所以此類方法使用較少。隨著CFD 的發展，目前壁面邊界條件通常有兩種不同的施加方式：在每個計算步內直接修正邊界節點的值，該方法稱為強邊界條件；或者在邊界節點處建立差分方程，通過計算邊界節點通量更新邊界節點狀態，該方法稱為弱邊界條件。強邊界條件利用邊界附近內部流場節點信息，通過外推插值等方式獲得邊界面的節點物理量。常用的強邊界條件包括壓力外推邊界和特征邊界[7-8]。弱邊界條件則需要在節點上建立差分方程，本文依據離散方式的不同將弱邊界條件分為兩類：一類是基于流場內側信息，在法向方向構造單側差分格式進行求解[9]，該方法沒有對法向速度進行限制，因此無法嚴格保持物面無穿透的強制條件；另一類主流處理方法是將計算區域進行拓展，在物體內部設置虛擬節點，通過補充模板的方式將邊界節點作為內點進行求解。由于虛擬節點不是實際存在的，因此該方法的關鍵在于如何提高虛擬節點參數的合理性。最為常用的是基于“反射原理”提出的鏡像法（symmetry technique, ST）[7]，該方法將內部節點參數按照鏡面對稱的方式構造虛擬節點參數。此外，Dadone 等[10]在壁面附近引入旋渦流動假設模型，提出一種通過熵和總焓對稱分布假設構造虛擬節點的曲率修正對稱技術（curvature-corrected symmetry technique, CCST）。該方法使用局部動量方程獲得虛擬節點壓力，在此過程中引入了物體曲率的影響，可以更加真實地反映壁面附近的壓力和速度分布，CCST 基于格心型有限體積方法提出[10-11]，逐漸被推廣至高階譜方法[12]。近年來，在貼體網格基礎上，CCST 方法也被應用于浸沒邊界法的邊界計算中，稱之為虛擬體單元法（ghost body-cell method, GBCM）[13-14]，不過，目前還未見該方法應用于貼體網格有限差分法的相關報告。

調研文獻發現，當前研究中所有使用虛擬節點的物面邊界方法在確定虛擬節點過程中僅考慮了法向相鄰節點的影響，即僅使用流場內側法向節點信息構造虛擬節點參數，并沒有考慮壁面切向相鄰節點的影響。本文基于上述認識，在壁面節點上建立滿足邊界條件的約束方程，將虛擬節點參數表示為計算模板中法向和切向相鄰節點的函數，通過求解約束方程得到滿足邊界約束的虛擬節點參數。該方法中虛擬節點受壁面切向節點的影響，可以有效體現邊界曲率對流動的影響，算例的對比驗證表明：該方法具有良好的計算精度，可以有效降低近邊界附近的數值誤差。

1 傳統滑移壁面求解方法

1.1 壓力外推邊界條件

壓力外推是CFD 中最為常用的一類強邊界條件，文中用PI-n表 示，其中n為插值多項式階數。該方法在壁面處滿足法向動量方程，建立(?p/?n)|wall=0的假設模型，使用單側迎風格式對上述壓力梯度進行離散，通過多項式插值得到邊界節點壓力值。壁面法向速度強制為零，切向速度使用內側參數外插得到，再依據等熵條件確定壁面節點密度[7,15]。以二維問題中MUSCL（ monotone upstream-centered schemes for conservation laws）五點格式為例，壁面速度由式（1）確定，這種處理方式本質上是通過內點外推邊界值，在內側區域流場梯度變化較大的情況下，高階的外插算法容易導致邊界值“過沖”[10]。

式中：下標 (i,0)表 示壁面邊界節點信息， (i,1)和(i,2)則 為流場內側節點參數， n表示壁面節點處的單位法向矢量。此外，當使用高階格式，如五階WENO（weighted essentially non-oscillatory）等時，在邊界區域常采用不同于內部節點的格式進行計算。一種處理方法是基于單邊差分的邊界格式，在邊界內側節點處進行降階處理，使用三階迎風WENO/ENO 計算內側節點數值通量，邊界節點同樣使用內部節點壓力外推確定，該方法的具體內容參見文獻[15]。

1.2 虛擬節點邊界條件

弱邊界條件需要在邊界節點處建立差分方程，通過計算節點通量的方式對壁面參數進行更新。由于缺少計算模板，通常在計算區域外設置虛擬節點，采用補充模板的方式將邊界節點當作內部節點處理，邊界約束通過虛擬節點參數體現。最為常用的方法為ST，根據鏡面對稱原則確定虛擬節點的參數，如式（2～3）所示：

其中： ρgj、pgj和vgj分別為第j個虛擬節點的密度、壓力和速度矢量，j為垂直壁面方向的節點編號。使用反射技術構造虛擬節點的處理方式，壁面無穿透的條件是弱施加的，當式（2）和式（3）應用于無曲率平直壁面時，可以嚴格滿足法向無穿透條件，在物體存在曲率的條件下，鏡像對稱的處理方式會引入較大誤差，并且無法保證壁面節點的法向速度恒等于零。

2 滑移壁面邊界約束算法

調研文獻認識到，在現有基于虛擬節點的弱邊界方法中，虛擬節點的變量是直接由內側法向節點通過鏡像或者外推插值確定的。當該種方法應用于包含曲率的壁面邊界時，在計算邊界節點通量過程中，切向速度和法向速度不能完全解耦，即切向方向通量分裂后會產生非零的法向動量通量。因此，即便在計算初始時刻強制壁面法向速度vn=0，經過時間推進后也會導致vn≠0。本文方法以嚴格滿足壁面無穿透條件為出發點，增加法向速度隨時間導數 ?vn/?t=0的強制約束，保證離散方程在時間推進過程中嚴格滿足vn=0的邊界條件。按照上述思路構造的邊界算法兼顧弱邊界和強邊界條件的優點：與弱邊界算法一樣，邊界和內部節點使用統一的離散格式；同時與強邊界算法一樣，壁面數值解可以嚴格滿足vn=0的強制條件。

需要說明的是，有限差分法應用于包含曲率的壁面邊界的貼體網格時需要進行坐標變換，將控制方程由笛卡爾坐標系變換至計算坐標系后再進行格式離散。在此過程中容易引入坐標變換誘導誤差，文獻[16]中使用多種類型網格進行自由流保持測試，表明該誘導誤差廣泛存在于差分格式計算中。本文重點在于求解存在彎曲的滑移壁面邊界，因此，在對邊界條件進行對比研究前需要先消除坐標變換引入的誘導誤差。本文采用劉君提出的離散等價方程及其離散 準 則（discrete equivalence equation and its discrete rule，DEER）[17]，以坐標變換后得到的包含非零源項的控制方程為原始方程，按照與差分格式相同的準則離散對流項和坐標變換源項。公式推導和算例驗證表明：DEER 方法可以在任意擾動網格上實現自由流保持，消除了坐標變換引入的誘導誤差。DEER 方法的詳細推導參見文獻[17-20]。

笛卡爾坐標 (t,x,y)下二維守恒型Euler 方程：

上式中守恒變量和對流通量表示為：

求解過程中需要將方程從笛卡爾坐標系 (t,x,y)變換到計算空間 ( τ,ξ,η)，因此實際應用差分格式離散的是計算空間上包含源項S的控制方程：

本文不涉及網格運動，因此僅給出靜止網格條件下表達形式：

按照DEER 算法要求，對流項和坐標變換產生的源項需要使用相同差分格式離散，將式（15）合并可以得到滑移壁面的邊界約束方程：

1）對于無黏流動，假設虛擬節點 ?η與內部節點η存在總焓相等條件：

3）根據內部節點 η=1的法向速度方向，分別使用兩種模型對虛擬節點的切向速度進行計算，如圖1 所示，其中黑色箭頭方向為 η=1節點處的速度方向。

圖1 切向速度假設模型示意圖Fig.1 Sketch of wall tangential velocity

圖2 DEEB 方法求解流程圖Fig.2 Flow chart of DEEB method

（1）將 η=1節點處法向速度方向指向流場內部的流動定義為切向速度模型一，使用多項式插值得到虛擬節點的切向速度：

其中s為計算模板中法向方向已知節點個數，cj為多項式插值系數。對于一階迎風格式，s=2，則虛擬節點的切向速度表示為：

（2）將 η=1節點處法向速度方向指向壁面的流動定義為切向速度模型二，此時在壁面附近可能存在激波，高階的外推插值可能會導致過沖，在這種情況下，將切向速度進行鏡像處理，表示為：

需要說明的是，式（19）和式（21）是在壁面存在曲率條件下設立的假設模型，對于平直壁面，無論內點η=1處法向速度方向正負，都直接使用式（21）確定虛擬節點切向速度，此時本文方法等價于傳統的鏡像方法。本文邊界算法在曲線坐標系下的原始方程是離散等價方程，將壁面無穿透條件表示為時間導數恒為零的約束方程，因此，稱其為“基于離散等價方程的約束邊界算法（discrete equivalent equation with boundary constraint algorithm，DEEB）”。

在二維問題中，每個虛擬節點需要確定四個參數，使用一階迎風格式時，通過式（16～19）或式（21）即可獲得虛擬節點的參數。由于約束方程式（16）與離散格式和通量分裂格式相關，無法給出統一的顯式表達式，本文使用弦截法求解該非線性方程。對于二階或者高階格式，為了確定更多的虛擬節點參數（二階需要兩個虛擬節點，五階WENO 需要三個虛擬節點），需要補充壓力遞推條件。在無黏流動中，壁面邊界本質上是一條貼壁流動的流線，可以利用計算模板中的內部節點多項式插值得到壁面壓力的左值，同樣利用虛擬節點可以插值得到壁面壓力的右值：

上文介紹的DEEB 方法主要是基于網格與壁面正交的條件，然而在許多問題的模擬中，網格難以保持和壁面正交化。對于非正交網格，可以利用離散等價方程計算中僅使用當地點坐標變換系數的特點[16-17]，在壁面節點處建立局部坐標系，其中 ξ方向保持不變，構造新的 η*方向與壁面正交，如圖3 所示。以空間一階格式為例，在該局部坐標系中， η*軸與j=1的原始網格線相交于圖中藍色節點處，使用反距離加權插值[21]得到交點u?的物理量。一階迎風格式在二維條件下使用交點周圍的5 個相鄰節點進行插值，插值模板表表示為：

圖3 非正交網格局部坐標系（一階迎風格式）Fig.3 Local coordinate system of non-orthogonal grids (first-order upwind)

其中： Δlm表 示模板中第m個節點到交點的距離，對于高階格式，只需要擴大插值模板即可。需要說明的是，模板中心點的選取準則是與 η*正交的原始網格線上距離交點最近的網格節點，即圖3 中節點 (n,1)。綜上所述，在壁面節點處構造正交的局部坐標系，然后使用式（16～19）或式（21）確定虛擬節點（空心藍色節點）參數，后續計算流程與正交網格下完全相同。

當DEEB 方法應用于五階WENO 等高精度格式時，拉格朗日插值多項式即使在光滑解條件下也可能存在穩定性問題[22]。本文使用Tan 等提出的WENO型插值方法，通過引入非線性權的方式確定插值模板中每個節點的貢獻，再使用高階插值得到局部正交坐標系中對應位置節點參數，該插值方法具體細節可參見文獻[22]。

雖然通過求解約束方程得到的虛擬節點參數可以提高邊界計算的精度，但是在這一過程中我們不希望引入過多額外的計算消耗，幸運的是，DEEB 方法僅在壁面節點處激活，且壁面節點相對于整個計算節點數量總是相對較少的。此外，隨著壁面附近流場區域逐漸收斂，由于虛擬節點參數保持恒定，此時約束方程的求解停止，通過分析可知，DEEB 方法的使用不會消耗過多的計算成本。

3 算例驗證

本小節以MUSCL3 格式[7]和五階WENO-Z 格式[23]為例，對本文所提方法進行算例驗證，以此考查DEEB 方法應用于不同離散格式的準確性和適用性。算例模擬中近似黎曼求解器采用Vanleer[24]格式，使用多步Rung-Kutta[25]法進行時間離散。需要說明的是，本文DEEB 方法中的約束方程是耦合在時間推進過程中的，因此可以使用現有的限制函數計算壁面附近存在激波的流動問題。由于本節算例目的在于對不同邊界計算方法進行對比，為了排除限制器引入的數值耗散，因此選取了不存在激波的光滑流場進行數值模擬。

3.1 MUSCL 格式計算結果

3.1.1 二維亞聲速圓柱繞流

首先使用馬赫數0.38 條件下圓柱繞流問題驗證DEEB 方法針對固定曲率問題的模擬能力，計算無量綱自由來流為 (ρ,u,v,p)=(1.4,0.38,0,1.0)，圓柱半徑r=0.5，計算區域取10 倍圓柱直徑，外邊界統一使用黎曼遠場邊界條件。計算區域使用4 套網格進行離散，網格周向和徑向節點個數分別為 40×10、 80×20、160×40和 3 20×80，將4 套網格定義為Mesh1～Mesh4。周向節點均勻分布，法向壁面附近局部加密。采用延拓網格對周向分布的周期邊界進行處理，排除該邊界對計算的影響。

圖4 不同方法計算得到的馬赫數等值線圖Fig.4 Mach number contours obtained by different methods

圖5 不同方法計算得到的熵增等值線圖Fig.5 Entropy contours obtained by different methods

圖6 中統計了不同方法的熵增誤差 L1范數隨網格細化的收斂曲線，N為網格節點總數，從圖中可知，PI-2 方法和ST 方法由于邊界處理引入了較大誤差，導致全場熵增收斂精度僅為 1.360 0和 1.090 0，本文DEEB 方法通過約束方程提高虛擬節點參數的合理性，使用與內部節點相同的離散格式進行求解，雖然測試結果表明DEEB 方法并未達到MUSCL 格式的二階設計精度，但是計算得到的誤差絕對值和收斂精度均優于兩種傳統求解方法。DEEB 方法沒有達到格式設計精度的原因可能是在邊界約束方程求解中使用切向速度模型（式（19、21））和壓力遞推假設（式（22））作為計算補充條件，在此過程中可能引入了模型誤差，這也是下一步研究中需要解決的首要問題。

圖6 不同方法的熵增誤差L1 范數對比Fig.6 L1 errors of entropy obtained by different methods

對于無黏亞聲速圓柱繞流，攻角為0°條件下，理論上流場完全對稱，因此常使用阻力系數評價計算對稱性。圖7 給出了4 套網格條件下使用不同方法計算得到的阻力系數，由曲線可知，隨著網格加密，阻力系數逐漸收斂至網格無關解，DEEB 在Mesh3 和Mesh4 條件下阻力系數相差微小，CD值約為 0.000 4，展示了其網絡收斂性較其他方法的優勢。

圖7 亞聲速圓柱繞流阻力系數Fig.7 Grid indenpence test of the drag coefficient of a subsonic cylinder

另一方面，阻力系數為積分量，反映的是圓柱整體壓力分布情況，而當局部正壓力誤差與局部負壓力誤差累加后得到的阻力值反而更小，因此認為阻力系數不能嚴格反映上下游壓力對稱性。為此，本文通過考查上下游物面網格對稱點壓力差進行對稱性評價。壓力差定義為：

如圖8 所示，其中pdown為下游物面點（橘黃色節點A?） 壓力，pup為 上游物面點（綠色節點A）壓力。本文計算中以上游駐點為起始點，定義順時針方向夾角為θ ，統計 θ ∈(90°，270°)物面上（藍色壁面）所有節點的δp值。

圖8 壓力對稱測點示意圖Fig.8 Sketch of symmetric pressure measuring points

壁面上下游節點壓力差分布曲線如圖9 所示，θ=180°對應下游駐點。從圖中可知，DEEB 方法壓力差明顯小于PI-2 和ST 方法。注意到ST 方法在駐點處不對稱性最為顯著，最大壓力差δpST,max=0.019 3，而DEEB 方法得到了最大壓力差僅為δpDEEB,max=0.000 4。此外，圖10 給出了收斂流場中壁面法向速度（使用來流速度進行無量綱化）分布曲線，DEEB 方法將法向動量通量為0（式（16））作為約束加入離散方程中，因此，時間推進過程中法向速度恒等于0；ST 方法存在非零的壁面法向速度，無法嚴格保證物面無穿透條件；PI-2 方法由于使用壓力外推的強邊界處理方式，在每一時間步強制壁面法向速度為0，所以與DEEB 方法得到的法向速度曲線重合。

圖9 壁面上下游節點壓力差分布曲線Fig.9 Pressure difference between symmetric measuring points

圖10 壁面法向速度分布曲線Fig.10 Azimuthal variation of the surface wall-nomal velocity

圖11 為Mesh3 對應的局部網格圖，其中紅色圓環所示節點為DEEB 算法激活區域，從圖中可知，約束算法僅在壁面節點處使用。表1 給出了該網格條件下三種方法使用Rung-Kutta 法顯式計算50 000 步所消耗的CPU 總時間Tall，其中以PI-2 方法的CPU時間為基準得到歸一化時間T?。可以看出，ST 方法和DEEB 方法由于使用虛擬節點，在邊界處求解差分方程的計算成本略高于壓力外推的邊界方法。在本算例中節點總數為6 400，其中壁面邊界節點數為160，這也是DEEB 方法求解邊界約束方程需要的節點數量，僅占總數的 2.5%，因此DEEB 方法不會引入過多的計算時間，計算總時間僅為PI-2 方法的1.055倍。

表1 二維亞聲速圓柱繞流CPU 計算總時間Table 1 CPU time for simulating subsonic inviscid flows arounda 2D cylinder

圖11 DEEB 算法激活區域（Mesh3）Fig.11 Activated region for DEEB method (mesh3)

3.1.2 二維平滑凸起通道內亞聲速流動第二個測試算例為馬赫數0.50 的自由流在平滑凸起通道內的流動[26]，來流攻角0°，用于考查DEEB 方法在任意曲率滑移壁面的模擬效果。計算網格量 200×50， 其中壁面方向設置 200個節點，豎直方向節點總數為 50。邊界條件設置如圖12所示，下邊界為光滑凸壁面，邊界形狀由函數y=0.625exp(?25x2)定義，其中 ?2 ≤x≤2。圖13 為局部區域網格圖，可以看出，下壁面凸起處網格與壁面非正交，因此可以考查本文方法在非正交網格條件下的計算效果。圖14給出了3 種方法得到的下壁面法向速度曲線，DEEB 方法在非正交網格條件下依然可以滿足法向無穿透條件，而ST 方法在凸起處存在較大的法向速度，無法保證壁面法向速度為零。圖15 顯示了3 種方法模擬得到的馬赫數等值線圖，除本文DEEB 方法外，其余兩種方法得到的等值線在凸起下游位置處發生拐折。

圖12 二維平滑凸起通道內亞聲速流動計算邊界Fig.12 Sketch of the computational domain and boundary conditions

圖13 局部區域網格圖Fig.13 Zoom-in view of local grids

圖14 下壁面法向速度曲線Fig.14 Wall-normal velocity at the slip wall

圖15 不同方法模擬得到的馬赫數等值線圖Fig.15 Mach number contours around a smooth bump obtained by different methods

圖16 給出了對應的熵增云圖，從圖中可以看出DEEB 方法再次得到了最小的壁面熵增，最大值僅為ΔsDEEB,max=0.000 38； PI-2 方法最大熵增為ΔsPI-2,max=0.000 73，且在凸起上游位置也存在明顯熵增，表明在流場變化較為劇烈的邊界處，壓力外推差值會引入較大 誤 差；ST 最 大 熵 增 ΔsST,max=0.002 2。本 例 和3.1.1 節測試結果均表明，在曲率滑移壁面條件下使用鏡像對稱的方式構造虛擬節點會嚴重降低計算精度。DEEB 方法通過求解約束方程得到了嚴格滿足邊界條件的虛擬節點參數，可以有效降低近邊界區域的計算誤差。

圖16 不同方法模擬得到的熵增云圖Fig.16 Entropy contours around a smooth bump obtained by different methods

3.2 五階WENO 格式計算結果

本節使用五階WENO-Z 格式進行離散，邊界附近分別使用邊界格式降階處理和虛擬節點兩種方法與DEEB 方法進行對比計算。其中，使用邊界格式的處理方法對壁面節點處壓力使用三階外推，密度依據等熵關系確定[15]，速度滿足法向無穿透條件，該方法表示為PI-3；補充虛擬節點的計算方法依舊使用ST 表示。

3.2.1 二維亞聲速橢圓繞流

使用二維橢圓繞流對本文方法在高精度格式中的計算性能進行測試，橢圓由方程x2/4+y2/9=1確定，為了保證全場亞聲速流動，自由來流馬赫數Ma=0.28， 沿橢圓周向分布 1 60 個節點，沿徑向分布50個節點，壁面附近網格如圖17 所示。同樣使用馬赫數和熵增作為不同方法計算能力的評價指標，從圖18 可知，DEEB 方法獲得了對稱性最好的馬赫數輪廓，PI-3 方法在下游駐點流線上產生了非物理凸起，而ST 方法誤差最大，下游位置難以形成對稱分布流場。這是因為五階WENO 格式需要構造三個虛擬節點，對于存在曲率的壁面邊界，使用鏡像法確定虛擬節點參數引入的誤差較大，進而對每一個包含虛擬節點的子模板均產生影響，最終導致了較大的數值誤差。圖19 給出了橢圓后緣位置處速度矢量分布，可以看出ST 方法產生了明顯的非物理渦結構，進一步說明了在高精度格式中，對存在曲率的壁面邊界應盡可能避免使用鏡像法進行邊界計算。從熵增等值線（圖20）中可以看出，DEEB 方法僅在曲率變化最大處存在微弱熵增，有效降低了邊界區域的數值耗散。PI-3 方法因為對邊界內側節點進行了降階處理，同時使用外推插值的方法確定壁面壓力和速度，雖然使用了等熵假設計算壁面密度，但是熵增還是遠大于DEEB 方法的模擬結果。圖21 的壁面法向無量綱速度曲線顯示了DEEB 方法應用于WENO 格式時依然可以嚴格滿足法向無穿透條件。

圖17 橢圓繞流局部網格圖Fig.17 Grids nearby an ellipticesh grids

圖18 橢圓繞流馬赫數等值線Fig.18 Mach number contours around an elliptic

圖19 橢圓后緣繞流速度矢量圖Fig.19 Velocity vectors at the trailing edge

圖20 橢圓后緣繞流熵增等值線Fig.20 Entropy contours around the elliptic

圖21 橢圓表面法向速度曲線Fig.21 Wall-normal velocity on the elliptic surface

最后，對不同方法計算得到的阻力系數（以橢圓長軸為參考面積）進行定量對比（CD,PI?2=5.56×10?4，CD,ST=9.74×10?2，CD,DEEB=1.20×10?4），DEEB 方法計算所得阻力系數最接近理論值，而ST 方法誤差最為顯著，這也與馬赫數等值線和熵增等值線結果相吻合，再次表明在存在曲率條件下，傳統鏡像法存在精度嚴重降低的問題。

3.2.2 Prandtl-Meyer 流動

計算模型和網格分布如圖22 所示，下壁面是Ma∞=1.357的自由來流經過20°折轉角形成的流線，網格線與壁面保持正交，流向和法向網格點個數為60×60，在第一道膨脹波前設置10 層均勻網格用于避免入口邊界處理的影響，上邊界和超聲速出口邊界采用外推處理。上側邊界距離壁面較遠，保證膨脹波在上邊界產生的誤差不會影響到下壁面，無量綱來流參數為 (ρ,u,v,p)=(5.4, 20/9, 0, 31/3)。

圖22 計算網格和邊界條件設置Fig.22 Grids and boundary conditions

圖23 和圖24 給出了不同方法模擬得到的流場等值線圖，圖中彩色等值線為計算結果，黑色等值線為相同標尺下對應的理論值。從圖中可以看出，ST方法得到的計算誤差最大，在壁面附近等值線存在顯著的非物理彎折，而本文DEEB 方法得到的結果與理論值吻合度最高，在壁面處也可以保持等值線的平直。

圖23 不同方法計算得到的速度u 等值線圖Fig.23 u contours obtained by different methods

圖24 不同方法計算得到的壓力等值線圖Fig.24 Pressure contours obtained by different methods

以圖25 為例對PI-3 方法的誤差產生機理進行簡要分析，在這一方法中，壁面壓力和切向速度利用內部法向節點三階外推插值獲得，雖然壁面節點（紅色節點）在第一道膨脹波位置處受到膨脹波的影響，但是此時內部節點（藍色節點）位于膨脹波波前，使用內部節點外插得到的壁面節點仍然保持均勻流參數，而這與理論流場顯然是不相符的，因此在膨脹波初始位置處壓力外推插值的方法存在誤差極值。

圖25 PI-3 方法外推插值誤差分析Fig.25 Extrapolation error of PI-3 method

同樣的，使用圖26 所示流場對ST 方法的誤差機理進行分析，圖中彩色等值線為壓力理論值，壁面節點計算時需要補充三個虛擬節點（j=?1,?2,?3）以進行通量求解，按照傳統的對稱反射技術（即式（2～3））可得：

圖26 ST 方法虛擬節點構造誤差分析Fig.26 Ghost node construction error of ST method

Prandtl-Meyer 流動的一個顯著特點是在來流參數一定的前提下，同一偏折角度對應的特征線上所有變量值均為定值，因此將圖26 中特征線（理論等值線）延長至物體內部后可以對虛擬節點參數的合理性進行判定。圖中結果顯示j=?1節點位置處對應的理論值介于p=7.0和p=7.5兩條等值線之間，即滿足7＜p?1＜7.5條件。而使用傳統對稱方法獲得的虛擬節點壓力p?1=p1＞7.5， 同理j=?2和j=?3兩個虛擬節點存在相同的問題，且距離壁面越遠的虛擬節點誤差越大，這也解釋了在五階WENO 格式中使用ST 方法的計算效果遠差于MUSCL 格式的原因。

此外，通過提取沿下壁面的壓力誤差曲線來定量對比不同方法的計算效果，其中，壓力誤差定義為：

其中：pi,j為 壓力的數值解，pei,j是 對應節點的壓力理論值，圖27 顯示了不同方法壁面處壓力誤差曲線，可以看出所有方法在第一道膨脹波附近都存在較大誤差。分析認為此處流場由均勻流開始膨脹，雖然流場參數保持連續，但是參數的導數存在間斷，因此會導致誤差極值。本文DEEB 方法在此處的最大誤差僅為0.001 6，遠小于PI-3 方法和ST 方法的壓力誤差。此外，圖28 給出了壁面節點的法向速度曲線，DEEB 方法雖然使用了虛擬節點，但是因為約束方程的作用，依舊可以保證法向無穿透的強制約束條件。上述等值線圖和壓力曲線分析表明：DEEB 方法可以直接應用于五階WENO 格式，通過求解約束方程獲得的虛擬節點參數更為合理，邊界節點使用與內部相同的離散格式，可以在滿足邊界約束的前提下，有效提高壁面附近流場計算的準確性。

圖27 壁面壓力誤差曲線Fig.27 Relative error of the wall pressure

圖28 壁面法向速度曲線Fig.28 Surface wall-normal velocity

4 結 論

本文以壁面無穿透條件為約束準則，在構造虛擬節點過程中引入了壁面相鄰切向節點的影響，通過求解約束方程獲得了合理的虛擬節點參數。在此基礎上，提出了一種面向有限差分方法的無黏流動邊界算法，該方法既能保證數值解嚴格滿足無穿透條件，又能在邊界點使用與內點完全相同的離散格式。通過亞聲速和超聲速算例模擬驗證了本文方法的準確性：

1）首先，針對MUSCL 格式開展數值模擬，通過亞聲速圓柱繞流算例與傳統計算方法進行對比，發現DEEB 方法能夠使在圓柱下游處馬赫數和壓力等物理量對稱性得到明顯改善，近邊界區域的熵增也顯著減小。

2）其次，使用平滑凸起通道內的亞聲速流動算例考查了本文方法在任意曲率壁面邊界條件下的模擬能力。

3）最后，通過橢圓繞流和存在解析解的Prandtl-Meyer 流動對本文方法在WENO 格式下的適用性進行了驗證。數值結果表明：DEEB 方法可以有效改善傳統虛擬節點方法誤差較大、無法滿足壁面無穿透的缺點，降低邊界計算引入的數值誤差。

本文工作為有限差分格式的物面邊界求解提供了一種新的思路，當前研究僅針對無黏流動和壁面處僅存在法向速度約束的情況，后續將對邊界約束方程進行推廣，開展針對黏性物面的邊界算法研究，相關工作目前正在進行中。