基于離散等價方程的非結構網格有限差分法

劉君，陳潔，韓芳

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

感受性問題是研究高超聲速邊界層轉捩過程的關鍵，精確刻畫相對來流參數10-4量級的聲波、熵波和渦波擾動經過激波的變化特性對現有的數值模擬是個挑戰。實際上準確模擬激波一直是計算流體力學(CFD)的重要研究內容，早期以激波裝配法為主[1]，自Harten[2]提出TVD格式以來，在超聲速流動應用領域激波捕捉法一枝獨秀，取得了ENO[3]、WENO[4-9]、CNS[10]，WCNS[11]等一系列標志性成果。捕捉法假設激波為空間連續函數，計算必然得到數值解表示的過渡區。這個過渡區對于Euler方程來說不存在，對于Navier-Stokes方程來說常常比理論值大幾個數量級(激波厚度僅為幾個分子自由程)，因此捕捉法面對感受性問題顯得力不從心。幾乎所有高精度捕捉格式需要構建或選擇限制器，通過降階來實現計算穩定，Zhong等[12-14]使用五階WENO格式模擬包含有激波的流場，數值結果表明精度很難超過一階，為此，他們研究感受性時選擇了激波裝配法。

裝配法把激波當作間斷處理，將R-H關系式直接嵌入流場內部。從Euler方程出發模擬相當于精確解；對于Navier-Stokes方程來說，激波厚度與飛行器尺度或網格尺度相比小幾個數量級，可以忽略，當作間斷處理不影響激波上下游參數的精確性。這種理論優勢為發展全場一致高精度算法提供了可能。近年來Zhong等[12-14]采用激波裝配法研究高超聲速邊界層轉捩的感受性問題的工作受到國內湍流研究專家周恒院士的關注。周恒和張涵信[15]指出“對于復雜激波系，用激波裝配法研究擾動通過激波經歷的變化是不現實的”。因此，將激波裝配法應用于高超聲速邊界層轉捩的感受性問題研究，首先需要解決復雜幾何拓撲結構的網格描述問題。

文獻[16]回顧了裝配法發展歷史和綜述了最新進展，本文作者劉君等受邀介紹課題組提出的捕捉/裝配混合算法(Mixed Capturing and Fitting Solver，MCFS)，體現了非結構動網格技術在解決復雜幾何拓撲結構激波裝配方面的優勢。但非結構網格有限體積法很難超過二階精度，用于感受性模擬精度較低，目前高精度格式主要應用于有限差分法。嵌入激波等間斷后流場空間被分割為多個相鄰區域，原則上可以應用基于結構網格的多塊拼接網格技術，但激波相交點等局部區域網格正交性較差且結構網格變形易于扭結，全部采用結構網格描述較為困難。本文的目的是發展基于非結構網格的有限差分法，解決局部復雜幾何形狀區域的網格拼接和激波邊界附近的網格變形難題。

近年來相關學者采用無網格算法的思想，在非結構網格基礎上利用點云擬合局部基函數來構造相應的廣義有限差分法(Generalized Finite Difference，GFD)[17-26]。這類方法在連接局部相鄰點的時候不需要建立局部的貼體坐標系，而本文的算法是基于曲線貼體坐標下的控制方程，本質上是目前常用的有限差分法的推廣，二者構造思路差異明顯。關于GFD的研究進展可以參考國內任玉新團隊的研究[27]，他們在這一領域的研究成果得到了同行認可。

1 離散等價方程及其離散準則

在前期結構網格的研究中，通過理論推導得出貼體坐標系下帶有源項的離散等價方程是有限差分法控制方程的結論，并且提出相應的離散準則[28-29]，這些研究成果是本文研究非結構網格有限差分法的基礎，下面進行簡單介紹。

考慮笛卡爾坐標系(t,x,y,z)下三維守恒型Euler方程：

(1)

從(t,x,y,z)變換到計算空間(τ,ξ,η,ζ)上，得到含源項的形式：

ItU+IxF+IyG+IzH=S

(2)

式中：

E為單位體積流體總能量，即

U、V、W稱為逆變速度，即

在坐標變換函數空間導數連續條件下S=0，因此CFD經典教科書給出方程形式：

(3)

根據離散等價方程的理論推導過程，分析了源項的產生機理，提出了在S≠0條件下離散方程的相容性準則：源項中?/?ξ和?/?η的離散格式應該與對流項的離散格式完全匹配。

考慮本文驗證算例是空間二維的，介紹離散算子時采用如下形式更為合適：

(4)

采用統一算子δ1離散左右兩側，考慮到通量和變換系數呈線性函數，可得離散形式為

代入方程(4)得到半離散形式為

(5)

(6)

如果把方程(6)中換成當地系數，就是經典CFD教科書方程(3)的算子形式,即

(7)

對于變形網格，時間方向導數中離散時要考慮J的變化，按照以上離散準則，如果左側時間導數采用一階顯式格式，源項中J的導數也要匹配：

按照CFD領域習慣，RHS表示空間離散格式計算結果，時間一階顯式格式寫為

(8)

結構網格應用到激波裝配法時遇到分區邊界運動會引起網格變形問題，在研究幾何守恒律的過程中提出以上離散等價方程及其離散準則[28-29]，注意到構造出的有限差分格式表達式中僅用到當地幾何參數，利用該特點可以建立基于非結構網格的有限差分法和激波裝配法。

2 非結構網格有限差分法

在二維情況下，如果從離散點出發只有2條網格線，對亞聲速流動無法應用一階迎風格式；如果從離散點出發有4條網格線，采用一階迎風格式離散正合適；如果從離散點出發有4條以上網格線，選擇正交性較好的4條線采用一階迎風格式離散也沒有問題。目前還沒看到將3條線作為局部曲線坐標的離散算法，這是本文研究重點。

盡管從物理空間(t,x,y)到計算空間(τ,ξ,η)的映射需要點對點的對應關系，但是，變換函數在方程(4)中體現為導數，數值方法計算這些導數只需要相對位置，如圖1所示，1-0-2或1-0-3連成的曲線平移不影響其斜率，因此，把1-0-2映射成計算空間的ξ、把1-0-3映射成η，網格節點位置確定以后，采用中心格式計算如下變換系數的導數不會出現數學奇性：

(9)

圖1 三相鄰節點的坐標變換對應關系Fig.1 Coordinate transformation relationship of three adjacent nodes

圖2 不同相鄰節點數的節點連接關系Fig.2 Node connection relationship of different adjacent nodes number

如圖2所示，原則上采用沿著3條網格線方向拓展模板節點來應用高精度格式，從第1點分支到第4、5點，可以構造二階迎風格式，但此策略又出現選擇組合，影響計算效率，無應用價值。計算空間多維的有限差分法格式大多由一維格式直接推廣而來，圖2中包含5、6、7點的圖形是二維結構網格，構造二階迎風格式只能沿網格線延展，模板中有第5、6點，沒有包含距離更近的第7點，越是高精度格式這種舍近求遠的現象更嚴重，定性分析這是不合理的。本文提出的非結構有限差分法解決了3條網格線基礎上構造一階迎風格式的問題，正在探索如圖2最后一個圖所示的在6條網格線基礎上構造二階迎風格式。

3 演示算例

本文發展非結構網格有限差分法是為了提高激波裝配法的實用性，未來的應用方向主要考慮超聲速流動，因此選擇如下存在理論解的斜激波和激波相交問題來驗證本文算法。

3.1 捕捉法

首先計算激波角β=40°的斜激波。圖3(a)是在笛卡爾坐標系下的直角結構網格上數值模擬得到的結果，稱之為結構網格算例，作為非結構網格計算結果的標準；圖3(b)中每條線長度相等、夾角相等，構成各向同性的六邊形，稱之為標準非結構網格；圖3(c)中線段長度相等、夾角不同，其中兩條網格線夾角為180°，稱之為磚墻非結構網格，考查夾角各向異性影響；圖3(d)網格是在標準非結構網格基礎上引入隨機誤差產生，考查網格質量的影響，稱之為任意非結構網格。以上4種網格的節點數均為1 650。圖3中還給出一階迎風格式計算得到的密度云圖，非結構網格有限差分格式能夠分辨出流場內激波結構，與結構網格的結果基本相符。

圖3 斜激波：4種網格密度云圖Fig.3 Oblique shock wave: density contours of four types of mesh

在固定高度y沿著x方向采用 Tecplot 軟件提取的密度ρ分布曲線和理論值如圖4所示，從圖中定量比較看，不同非結構網格計算的激波位置和結構網格結果一致，過渡區寬度也較為接近，表明精度接近常規的結構網格有限差分。

各算例的收斂曲線如圖5所示，單調下降，沒有出現異常波動，表明從離散等價方程出發的非結構網格離散格式是可以得到收斂解的,圖中t為無量綱時間。

為進一步考查3條網格線的非結構網格有限差分法的算法對網格的適應性和收斂性，選擇含有激波相互干擾結構的流場作為數值模擬算例。如圖6所示，兩道入射激波激波角為β1=41°和β2=30°，相交以后形成兩道激波和一條接觸間斷，流場分為5個區域，其中，Ⅰ、Ⅲ和Ⅴ區采用正交性較好的傘形網格，可以看出，存在網格線相交的奇性點和網格疏密分布不合理這兩個突出的問題，Ⅱ和Ⅳ區采用平行四邊形網格，在激波相交和小夾角情況下難以滿足正交性要求。因此，很難生成高質量的結構網格。

圖4 斜激波：密度隨x變化曲線Fig.4 Oblique shock wave: curves of density change along with x

圖5 斜激波：殘差收斂曲線Fig.5 Oblique shock wave: convergence curves

圖6 異側激波相交流場結構網格空間離散Fig.6 Flow field spatial discrete of opposite side shock wave intersection with structural mesh

同樣采用前面4種形狀的網格進行數值模擬計算，網格分布和流場密度云圖如圖7所示。圖8為固定位置x處沿著y方向的密度分布及其理論值曲線。從圖中可以看出，在不同形狀的非結構網格上計算得到的激波位置和過渡區寬度基本一致，僅在下游強激波一側和常規結構網格有差異，在其他區域非常接近。圖9是各算例的收斂曲線。

圖7 激波相交：激波捕捉法中4種網格密度云圖Fig.7 Shock wave intersection: density contours of four types of mesh

圖8 激波相交:密度隨y變化曲線Fig.8 Shock wave intersection: curves of density change along with y

以上兩個算例驗證了本文算法的有效性。下面通過裝配法展示應用前景。

圖9 激波相交：殘差收斂曲線Fig.9 Shock wave intersection: residual convergence curves

3.2 裝配法

按照GFD格式，在考慮網格點運動情況時，模板函數增加了時間變量，處理起來較為困難，本文從離散等價方程出發建立的相容性準則對時間和空間均可使用。在下面計算中按照方程(8)處理變形網格引起的幾何守恒律問題。

模擬馬赫數Mas=3.0的非定常激波在靜止大氣中運動，初始時刻、無量綱時間t=0.000 63和t=0.001 26的流場如圖10所示，在此基礎上將激波后流場換成β=40°的斜激波波后參數，保持下邊界的流動參數不變，計算從初始正激波收斂到斜激波的過程，初始時刻、無量綱時間t=0.001 09和充分收斂時間t=0.16的流場如圖11所示。以上兩個算例，盡管網格不規則，但是激波陣面保持平穩，運動速度和波后流場參數沒有波動，而且和理論值完全相同，這一結果表明本文算法滿足幾何守恒律，在裝配激波引起的網格變形過程中不會產生誤差。

圖10 運動正激波算例Fig.10 Moving normal shock wave sample

圖11 斜激波算例Fig.11 Oblique shock wave sample

采用激波裝配法數值模擬前面的非對稱激波相交問題裝配法，計算結果如圖12所示，和理論解完全相符。為了考查算法的網格敏感性，針對入射激波角相等(β=40°)的對稱流動，在網格中加入隨機擾動，計算結果如圖13所示，即使上下網格不對稱，計算得到的流場也非常對稱，很好地體現出裝配法精確刻畫激波所產生的效果。

圖12 非對稱激波相交密度云圖Fig.12 Density contours for asymmetric shock wave intersection

圖13 對稱激波相交密度云圖Fig.13 Density contours of symmetric shock wave intersection

4 應用展望

早期的邊界激波裝配法(Boundary Shock-Fitting Method, BSFM)把激波作為計算區域的邊界采用動網格技術進行跟蹤，這種算法適合于位置大致確定的鈍體擾流的脫體激波。為了模擬復雜內部激波結構和在計算過程中激波發生大變形的情況，Moretti[16]提出了在固定背景網格基礎上的浮動激波裝配法(Floating Shock-Fitting Method，FSFM)，這種算法類似于嵌入邊界法(Immersed Boundary Method, IBM)，背景網格常用笛卡爾網格，也可以用非結構網格[1,16]。為了解決多點模板WENO格式應用難題，Zhong等對浮動激波裝配法進行改進，提出激波陣面追蹤法[12-14](Front-Tracking Shock-Fitting Method, FTSFM)，其原理如圖14所示，首先判斷激波面和網格的交點Si±k，再將這些點擬合得到激波面曲線，在此基礎上建立曲線坐標系采用單側WENO模板離散激波上下游，例如，x方向激波上游低壓區Si處的格式模板涉及A～E點，時間推進到下一時刻(圖中虛線)，需要重新判斷激波位置。

圖14 激波陣面追蹤法Fig.14 Front-tracking shock fitting method

下面通過和Zhong等的激波裝配法比較，來說明本文方法的特點：

1) Zhong等的方法基于結構網格，將激波限定在兩個網格點之間，兩個流動參數只能準確描述一種間斷，原則上不能分辨和裝配前面算例中激波相交點處存在多個流動參數的情況，在目前也沒有看到這種方法在類似問題中應用。本文算法采用非結構網格，不限制網格點連接的線段，理論上可以描述任意Riemann結構，在處理激波相交等復雜幾何方面更有優勢。

2) 本文算法的格式也是從曲線坐標體系下方程出發離散，沿著激波切向拓展模板，如圖15(a)中黑線所示，沿著激波法向拓展單側模板，如圖15(b)中綠線所示，和Zhong等的激波陣面追蹤法一樣，也可以應用高精度格式。本文算法在每個網格點上計算量較大，但由于采用動網格裝配跟蹤間斷，計算中不需要每一步辨識和擬合激波，比Zhong等的方法更適合于并行計算。因此，二者的計算效率還有待于針對具體應用問題的比較。

3) 直線激波的上下游參數為常數，一階格式也能計算得到精確解，換言之，對于激波附近的流場，裝配法一階格式的計算結果也比捕捉法任何格式的精度高。激波變彎曲來自于切向參數的非均勻，激波兩側流動參數的法向導數小于切向導數。如圖14所示，Zhong的方法引入當地坐標系的目的是準確表達R-H關系式，WENO格式模板A～E點連成的網格線與激波并沒有正交，因此，在格式離散過程中無法區分激波切向和法向。本文基于非結構網格的算法可以保證網格與激波的正交性，這樣計算中不同方向采用不同格式，有助于提高計算精度。

根據以上分析，定性分析本文算法在超聲速邊界層轉捩過程的感受性問題研究中的應用策略：

圖15 本文算法高精度模板和拼接網格計算Fig.15 High order template and stitching mesh computing of the algorithms in this paper

1) 對于頭部脫體激波，在現有結構網格基礎上先用捕捉法計算，然后采用流場結構辨識算法確定激波大致位置，局部區域采用本文的非結構網格進行激波裝配，激波切向直接采用高精度格式、法向可以從一階格式逐步過渡到高精度格式。

2) 激波前后法向速度間斷、切向速度不變，利用這個特性在模擬飛行環境下自由來流擾動時，可以根據聲波、熵波和渦波的傳播特點進行方向分解，然后精確地施加在激波上游。

3) 對于激波相交等幾何拓撲復雜的流場，局部用本文方法來解決相交點的多值問題和激波兩側的網格正交性問題，其他區域采用常規的結構網格有限差分法。

除了高超聲速的感受性研究的應用前景，本文算法還有其他用處。如圖15(b)所示，對應用多塊拼接網格和加密網格常遇到的相鄰區域網格節點不對應的情況，目前的算法只能作為特殊邊界處理。前面的磚墻非結構網格已演示，本文算法可以采用與內點統一的計算格式。

5 結束語

本文研究是在周恒院士的文章[15]啟發下開展的，基于離散等價方程及其離散準則提出一種新的非結構網格有限差分格式構造算法。從算例演示來看，能夠解決激波裝配法遇到復雜波系時引起的幾何拓撲難題，通過對節點引入擾動，證明了本方法的穩定性和對畸形網格處理的有效性，基本上實現了預期目標，為開展高超聲速邊界層轉捩過程的感受性研究打下基礎。