一種用于間斷裝配的流場結構辨識算法

劉 君,陳澤棟,陳 潔,鄒東陽

(1.大連理工大學 航空航天學院,大連 116024;2.大連理工大學 工程力學系,大連 116024)

0 引 言

流體微團加速到超聲速以后無法感知到達之前遇到的空間物體，只能在物體前通過急劇壓縮的激波才能把速度降下來。因此，激波是可壓縮流動的基本特征。早期計算流體力學(CFD)的算法無法模擬這一現象，直到1954年Lax提出處理包含有間斷流動問題的弱解理論后[1-2]，國內外眾多科學家開始研究在計算過程中能自動分辨出激波的算法，提出TVD、NND、ENO、WENO等具有里程碑意義的激波捕捉格式[3]，解決了航天航空等領域很多工程模擬問題，極大提高了CFD地位，使其發展為一個重要的學科分支。盡管建立這些格式時考慮了激波間斷，但是本質上還是1950年von Neumann提出的把間斷描述成流動參數連續變化過渡區的思路[4]，構造出來的格式存在以下難以克服的理論困境。

數值計算的離散格式是對偏微分方程的近似，按照CFD收斂性理論，數值解隨著網格尺度減小趨向理論解，CFD追求的極限是方程的理論解，原則上說對于方程不能描述的流動現象，數值解和實驗值不具有可比性。海平面環境下馬赫數小于10的激波厚度大約1×10-8m，流體在這么小的空間內受到劇烈壓縮，暫且不說準確模擬所需要的網格尺度及其帶來的計算量，從連續介質假設出發建立的N-S方程是否適用于激波模擬也存在爭議，因為Stoke假設(3λ+2μ=0)不再成立[5]，計算時需要對體積膨脹黏性系數λ進行修正才能得到與實驗相符的激波內部溫度[6]。從事稀薄氣體動力學研究的學者認為激波內部結構涉及到熱力學非平衡，連續介質假設本身就不成立，“從理論上說，Navier-Stokes方程不能用來描述激波結構的[7]”。無黏流體力學把激波表示成間斷，兩側參數和運動速度滿足蘭金-許貢紐(Rankine-Hugoniot,R-H)關系式，這是迄今為止幾乎很少有爭議的理論。所謂間斷在數學上就是在同一空間位置對應至少兩組(激波前后)流體變量，在多維情況下可能更多，例如二維馬赫反射三岔點有四組流體變量。很多基于弱解理論構造的激波捕捉格式考慮了間斷，但是用參數連續變化的思路描述激波注定不可能計算出真正的數學間斷，捕捉到的激波主要是現象模擬。

按照上述觀點，常規N-S方程及其理論解不能用來描述激波內部物理現象，那么基于離散格式的數值解至少在過渡區內沒有物理意義，采用Euler方程產生的后果是無法評價過渡區內數值解的精度(激波前后參數相差很大)，激波捕捉格式對于激波相交等現象更是無能為力。

20世紀60年代Moretti等人將激波裝配方法和Lax-Wendroff格式結合，在CDC 600型計算機上進行，總耗時為360 s，給出了超聲速鈍頭體問題的數值解，2002年Moretti回顧從業36年來CFD領域模擬激波的進展，針對超聲速鈍頭體繞流比較了捕捉法和裝配法計算結果，認為準確模擬激波還需要發展裝配法[8-9]。對于50年前提出的裝配法，文獻[10]的評價具有代表性——“激波裝配法具有計算精度高、物理概念清晰等優點，但是計算過程復雜、計算量大，并且必須事先知道激波的大概位置。一般來說，激波裝配法僅適用于那些激波比較簡單而事先可以估計激波形狀和位置的定常流動問題。”

造成計算復雜的主要原因是求解過程中激波運動和變形。裝配法減少了網格量，造成計算量大的原因是激波及其相交點等空間網格點需要判別、標記或特殊處理，增加代碼編寫復雜性和并行計算難度。隨著CFD網格技術發展，這些難點也得到緩解。國外文獻把最早Moretti提出的稱為邊界激波裝配法 ，僅能處理形狀較為規則和位移不大的脫體激波。后來Richtmyer等人[11]提出了浮動激波裝配方法，在靜止網格基礎上嵌入R-H關系式，允許激波在背景網格上滑動，解決了內部激波裝配難題。近年來Zhong等人將這種方法和高精度有限差分格式相結合有效提高內部光滑流場的計算精度，成功應用于許多高超聲速流動的求解中[12-13]。基于結構網格的邊界激波裝配法和浮動激波裝配法在模擬激波相交和馬赫反射等復雜波系結構時面臨一定困難。Paciorri和Bonfigolioli等人[14-16]采用非結構網格技術和局部網格重構提出了激波陣面追蹤法，記錄激波在背景網格上的位置，對激波附近的背景網格進行挖洞，然后再利用Delanay方法在洞內重構以激波陣面分界的新網格，時間推進過程中根據R-H關系式得到新的激波位置，重新判別、挖洞和填補背景網格。這種算法不需要考慮網格節點之間的連接關系，也就減少了對于特殊點進行處理的麻煩，增加了代碼的普適性，但是在激波邊界節點對應的網格單元只能采用一階插值，像常用的多塊重疊網格技術一樣存在插值引起的精度降低問題。其次，不允許激波在一個時間步內跨越多個網格節點，通常選擇足夠小的時間步長，因而對計算效率有所影響。上述這3類裝配法時間推進過程中需要判別激波位置，增加了計算量，而且激波前后網格數目是變化的，不利于并行計算。2015年劉君等人[17]提出了基于非結構動網格的激波(間斷)裝配方法，從能夠描述網格變形的任意拉格朗日-歐拉坐標系(Arbitrary Lagrangian-Eulerian,ALE)下的控制方程出發，得到激波位置及其兩側的參數后，裝配的激波就變成了一個移動邊界，對于存在內部激波的計算域則變成類似于拼接網格的多個分區，激波運動由上下游參數驅動，采用成熟的動網格技術[18]進行移動和變形。與Paciorri及Bonfigolioli的裝配法相比，在保留非結構網格的優點外，新的裝配法不需要在每一步對計算網格進行挖洞重構，避免了插值誤差，易于并行計算，并且時間步長不受限制。采用這種算法模擬了內部存在正規反射和馬赫反射等復雜激波結構的算例[19-21]，從效果看，較好地解決了傳統裝配法遇到的“計算過程復雜、計算量大”問題。

為了突破“激波裝配法僅適用于那些激波比較簡單而事先可以估計激波形狀和位置的定常流動問題”的限制，劉君等人提出了在采用捕捉法計算的流場基礎上辨識并裝配主要激波結構、計算過程中自動裝配內部細致結構和新生激波的解決思路。在技術途徑上提出了新的數據結構，在激波邊界通量計算時，如果采用Roe、van Leer、AUSM等格式就是捕捉法，如果采用R-H關系式就是裝配法，實現兩種方法之間任意轉換。按照上述思路，為了解決裝配法的應用難題，亟需發展激波、接觸間斷等流場特征結構的辨識算法。

目前，常用的激波辨識算法主要有三種[22]：第一種根據最大變量梯度確定激波面，該方法實現簡單，但需要設置合適的過濾器才能消除錯誤的結果；第二種為法向馬赫數法[24]，該方法認為在激波面上流動的法向馬赫數等于1，該方法與第一種方法相比適用性更強，但也需要設置適當的過濾器才能獲得良好結果，且難以用于探測接觸間斷；第三種為基于特征線理論的方法[25-26]，該方法與前兩種方法相比更精確、更有效，但實施過程復雜且計算量大。

綜上所述，這些激波辨識算法難以在復雜三維問題中廣泛應用，且上述算法辨識出的激波通常為具有一定寬度的激波帶，難以直接用于裝配法。因此，需要發展新的辨識算法來推動激波裝配法實用化進程。

本文提出了一種光源射線法，該方法可將捕捉法中探測到的具有一定寬度的激波帶擬合成一系列離散線段(二維)或三角形面片(三維)，利用曲線/曲面擬合算法進一步將其擬合成連續的曲線/曲面。該曲線/曲面可作為激波裝配法中的激波面，并且可根據計算需要對其進行任意形式的網格劃分。為便于表述，先用二維模型的實際操作過程進行原理說明，然后簡單介紹三維推廣過程，并展示實際工程問題的辨識結果。

1 流場結構辨識算法

1.1 二維辨識算法流程

(a)捕捉法流場 (b)激波單元 (c)激波曲線

(d)分片擬合

(e)確定特征點圖1 光源射線法示意圖Fig.1 Sketch of source-rays method

1.2 脫體激波辨識

針對二維超聲速鈍頭體繞流問題，采用1.1節所述方法進行激波面辨識。算例中，空氣來流馬赫數Ma=5.0、無量綱密度ρ=1、無量綱溫度T=1，來流攻角為0°。鈍頭體中心位置坐標為(0,0)，鈍頭體半徑為0.0254。網格為160×100的四邊形網格，其中沿鈍頭體法向的節點數為100。計算網格和捕捉法流場壓力云圖分別如圖2(a)和圖2(b)所示。

(a)網格 (b)壓力云圖 (c)激波點集 (d)激波面圖2 基于壓力梯度的激波辨識結果Fig.2 Shock wave surface extraction based on pressure gradient

考慮到很多高精度格式的限制器選擇壓力梯度來進行調節，在此也將其作為激波特征參數來探測脫體激波。對于捕捉法得到的壓力場，用格林公式計算出每個網格單元的壓力梯度矢量。作為間斷面，激波厚度理論上處于分子自由程大小這一數量級，其兩側壓力的梯度值為無窮大。但是，在捕捉法數值計算中，激波表現為帶狀分布的過渡區且具有較大壓力梯度。對于不同的問題以及不同的計算網格，激波帶上的壓力梯度值都有所不同。因此，難以通過一個固定的壓力梯度值將激波帶提取出來。本文為了準確獲得激波點集的分布，首先計算出流場中的最大壓力梯度值|p|max，然后針對不同問題設置相應的閾值系數k，壓力梯度閾值|p|thre=k|p|max，所有滿足條件|p|≥|p|thre的單元格心點構成激波點集T。Case1取閾值系數k=0.1，獲得離散點集，如圖2(c)所示。針對上述激波離散點集，采用光源射線法進行曲線擬合。光源S位置固定在(x,y)=(0.1,0)處，相鄰兩條射線間的夾角取常值φi=4.5°。辨識出的結果如圖2(d)所示，從圖中可以看出辨識出的激波曲線與捕捉法流場中的激波點集T吻合度很高，能準確描述激波形狀，對激波曲線重新進行網格離散即可作為激波邊界直接用于激波裝配法。

離散點集的分布取決于所選擇參數和閾值，為了比較不同壓力梯度閾值下所辨識出的激波曲線，在捕捉法流場的基礎上另外補充兩個辨識狀態：Case2取k=0.2，Case3取k=0.3。得到的點集及其最終辨識出來的脫體激波曲線如圖3所示。從3個狀態的辨識結果可以看出，不同的壓力梯度條件下所獲取的離散點集有所不同，導致最終辨識出的激波曲線長度有所區別，但三條激波曲線都能與捕捉法流場中激波的位置和形狀保持一致，將其作為激波裝配法中的激波邊界可快速得到收斂解。得到初始激波位置以后的裝配法計算方法和過程參見文獻[17,19-21]，這里不再贅述。

(a)Case1 (b)Case2 (c)Case3

圖3 不同壓力梯度閾值下的脫體激波辨識結果比較
Fig.3 Comparison of shock surface extraction results under different pressure gradient thresholds

1.3 內部激波及接觸間斷辨識

對波系更為復雜的馬赫反射問題進行數值模擬，計算區域及邊界條件如圖4所示。邊長為LAF=0.4、LAB=0.1、LDE=0.312、LEF=0.85，斜楔BC傾角為14°；邊界AF為超聲速入口、DE為超聲速出口，其余為滑移壁。計算網格為300×218的四邊形網格，其中入口和出口邊界上的網格節點數為300。來流馬赫數Ma=1.9，攻角為0°。

圖4 計算域Fig.4 Calculation domain

由于斜楔BC的偏轉作用，自由來流在點B處產生一道入射激波(Shock1)；激波達到壁面FE處時發生馬赫反射，產生一道反射激波(Shock2)和一條馬赫桿(Shock3)；反射激波在壁面CD處再次發生反射，產生第二道反射激波(Shock4)。流場穩定后的密度云圖以及波系結構如圖5所示。采用激波及接觸間斷探測方法可從流場中獲得主要的四條激波點集和一條接觸間斷點集分布。為了檢驗不同的探測方法對最終結果的影響，激波點集采用密度梯度作為參考量進行探測，接觸間斷點集由密度梯度和壓力梯度共同決定。同樣，先計算出全場最大密度梯度值|ρ|max和最大壓力梯度|p|max。激波探測的密度梯度閾值所有滿足的單元節點構成激波點集，這里取接觸間斷單元由密度梯度和壓力梯度共同決定，滿足密度梯度較大且壓力梯度較小，即其中所有滿足該條件的單元節點構成曲線擬合所需的接觸間斷點集。這里取其中在不同密度梯度閾值下辨識出的接觸間斷長度有一定區別，但位置和形狀基本一致。激波和接觸間斷離散點集分布如圖5所示。

獲得離散點集后，通過在區域內合理布置光源可擬合出每條激波和接觸間斷曲線。首先在Ⅱ區設置一個光源，光源發出等間距射線，采用上文光源射線法可辨識出表征激波1和激波2的曲線，如圖6(a)所示。同理，在Ⅲ區布置光源可辨識出激波4和接觸間斷曲線，如圖6(b)所示。在I區設置光源可辨識出馬赫桿曲線，如圖6(c)所示。所有曲線的最終辨識結果與流場密度云圖的比較如圖7所示?？梢钥闯觯孀R出的曲線與采用捕捉法獲得的波系結構高度吻合，直接用于裝配法有助于快速得到收斂解。得到初始激波位置以后的裝配法計算方法和過程參見文獻[17,19-21]，后文也會有簡單介紹，這里不再贅述。

圖5 捕捉法密度云圖及間斷點集分布Fig.5 Density contours and discontinuity point set distribution under shock-capturing schemes

(a)激波1-2 (b)激波4和接觸間斷 (c)激波3

圖6 激波和接觸間斷辨識過程
Fig.6 Shock wave and contact discontinuity surface extraction process

圖7 激波及接觸間斷辨識結果與流場比較Fig.7 Comparison between feature surface extraction result and flow field structure

值得注意的是，這里光源的位置分布并不唯一，可采用不同的分布方式達到相似的效果，比如可先在Ⅰ區設置光源同時擬合出激波1和激波3曲線，再在Ⅲ區設置光源得到激波2、激波4和接觸間斷曲線。

2 三維實現

2.1 三維辨識算法流程

本文通過引入單位球可確定各射線之間的連接關系，如圖9所示。對半徑為1的球面進行三角形離散，獲得一系列三角形面片，各頂點連接關系確定，且三角形面片滿足連續且不相交條件。

圖8 射線方向示意圖

Fig.8 Ray direction

圖9 單位球

Fig.9 Unit sphere

將單位球的球心作為光源S，球心與單位球上三角形面片的三個頂點ai、aj、ak的連線li、lj、lk構成一個子區域Ap，如圖10(a)所示。若Ap內離散點的個數不小于3，則對Ap內所有離散點進行平面擬合，擬合出的平面與三條射線的交點分別為bi,p、bj,p和bk,p；若該子區域內離散點的個數小于3，則跳過該區域。所有區域擬合完成后，每條射線li與不同子區域內的擬合平面形成若干交點。將各射線li上的交點取平均值獲得特征點bi，則各子區域Ap的三條射線邊界上的特征點bi、bj和bk構成目標三角形面片，如圖10(b)所示。最終，整個離散點集T被擬合成一系列連續的三角形面片，三角形面片頂點的連接關系由單位球確定。可對上述三角形面片構成的空間曲面重新進行網格劃分，獲得可用于裝配法的激波面網格。

(a)分區擬合 (b)三角面片圖10 子區域內三角形面片擬合過程Fig.10 Triangular fitting process in sub-region

2.2 三維測試算例

通過以下兩個算例驗證本文辨識算法在三維問題中的有效性。

在橢球面附近隨機生成100 000個點作為離散點集，隨機點坐標B(x′,y′,z′)滿足：

其中，

式中，i、j、k、l、m為區間[0,1]內的隨機數。

獲得的離散點集如圖11所示。單位球共有736個節點、1468個三角形單元(如圖9)，球心位置坐標為(0,0,0)。單位球位置及橢球面辨識結果如圖12所示。從結果可以看出，辨識出的曲面光滑完整且形狀與離散點分布形狀一致。

圖11 橢球隨機點集圖12 單位球及擬合結果Fig.11 Random point set of ellipsoidFig.12 Unit sphere and the extracted surface of ellipsoid

圖13 雙橢球隨機點集圖14 雙橢球面擬合結果Fig.13 Random point set of double-ellipsoidFig.14 Extracted surface of double-ellipsoid

2.3 三維繞流激波辨識

對三維超聲速球頭繞流進行模擬，來流馬赫數Ma=5.0,無量綱密度ρ=1.0,無量綱溫度T=1.0，來流攻角和側滑角均為0°。球頭的球心位置坐標為(0,0,0)，球頭半徑為0.0254。計算網格共有508 365個的六面體單元，球面網格如圖15(a)所示。圖15(b)為流場壓力云圖，圖15(c)為根據流場壓力梯度獲得的激波帶離散點集。算例中，壓力梯度閾值系數k=0.05離散點為滿足|p|≥|p|thre的所有單元節點，共有23 640個。單位球共包含1734個節點，3464個三角形面片，球心坐標為(2,0,0)。最終辨識出的激波面如圖15(d)所示，除了邊緣附近外的其他區域辨識結果與離散點集吻合較好，將邊緣位置稍作處理即可作為初始激波面用于激波裝配法。對于不同的初始位置的激波，最終收斂解相同，故本例不再展示裝配法的計算過程，計算方法參見文獻[17,19-21]。

(a)球頭表面網格 (b)捕捉法壓力云圖

(c)激波點集 (d)激波面圖15 三維超聲速球頭繞流激波辨識過程Fig.15 Shock wave surface extraction process of 3D blunt body flow

在球頭繞流算例的基礎上，對波系較為復雜的球柱錐組合體外形進行模擬，檢驗特征面提取方法處理三維復雜問題的能力。計算模型如圖16(a)所示，來流馬赫數Ma=4.04，來流攻角為20°。波系結構如圖16(b)所示，在最外層形成一道弓形激波，在錐裙與柱體的結合處形成一個內嵌激波。探測出的激波點集和提取出的激波面分別如圖16(c)和圖16(d)所示。

(a)計算模型 (b)捕捉法壓力云圖

(c)激波點集 (d)激波面圖16 球柱錐組合體繞流激波辨識過程Fig.16 Shock wave surface extraction process of cylinder with a hemispherical nose and a conical flow

3 裝配法實例

利用上述方法獲得初始激波面后，將其作為激波邊界用于激波裝配法，進行流動數值模擬，考查所述激波辨識方法在實際使用中的效果。

對于定常激波裝配法，初始激波位置及形狀會影響計算的收斂過程，但最終的收斂結果都相同。因此，本節直接采用文獻[19]算例的部分結果說明激波裝配法與捕捉法的區別。該算例為具有復雜波系結構的馬赫反射問題，來流馬赫數2.0，偏轉角δ=14°。計算區域為L×L的正方形，邊界上網格節點均勻分布，節點間距為0.01L，內部區域為三角形網格，單元數為22 429。流場波系結構及邊界條件如圖17所示，入射激波IS和反射激波RS、馬赫桿MS交于三波點TP，SS為接觸間斷。

首先采用激波捕捉法獲得初始流場，在此基礎上采用本文辨識算法辨識出激波及接觸間斷的位置和形狀，將其作為裝配法的初始間斷面，結合R-H關系式進行裝配法流場計算。得到初始激波位置以后的計算過程參見文獻[19]，計算結果完全一致。下面給出部分結果，如圖18所示，其中圖18(a)為采用捕捉法的計算結果，圖18(b)為采用裝配法的計算結果。從圖中可以看出，捕捉法流場中激波附近出現等值線的扭曲和震蕩，且激波呈帶狀分布，具有很寬的過渡區。與此相比，采用裝配法所獲得的流場中各區域之間界線分明，激波是一個沒有寬度的界面，激波兩側變量值嚴格滿足R-H關系式。

圖17 邊界條件及流場結構示意圖[19]Fig.17 Boundary conditions and sketch of flow field structure[19]

(a)捕捉法計算結果

(b)裝配法計算結果圖18 裝配法與捕捉法計算結果對比(馬赫數云圖)[19]Fig.18 Comparison between shock-capturing solution and shock-fitting solution[19]

4 結 論

通常對于具有一定“厚度”的空間點集，采用傳統的擬合方法難以將其擬合成形狀一致的光滑曲面，本文提出的光源射線法很好地解決了此類問題。

對于從捕捉法流場中獲取的具有一定“厚度”的激波帶和接觸間斷帶，采用該方法進行空間曲面擬合后可作為激波裝配法中的初始間斷面。尤其對于具有復雜波系結構的流場，該方法提供了一種裝配初始流場的途徑。本文詳細介紹了該算法在二維問題中的實現過程，并將其成功推廣應用于三維問題。算例結果表明，辨識出的激波曲面整體上與捕捉法流場中激波帶分布相吻合，但是三維問題中激波曲面邊緣附近的辨識效果有待進一步提高。對于復雜工程問題，往往存在較為復雜的三維波系結構，間斷面的提取過程更為困難和繁瑣，該方法還需進一步完善。

對于傳統激波裝配法，在初始激波面給定后需要經過若干步迭代才能收斂到準確的形狀和位置。尤其對于三維問題，收斂過程更加漫長，甚至出現難以收斂的情況。采用本文方法可提供較為準確的初始激波，有助于縮短收斂過程。

本文所提方法除了用于激波裝配法外，未來還可應用于自適應網格算法。根據不同的流動特征可擬合出相應的空間曲面，進而可以任意調整該曲面附近區域的網格分布，達到網格根據流動特征自動調整的目的。