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

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

0 引 言

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

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

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

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

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

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

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

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

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

1 流場結(jié)構(gòu)辨識算法

1.1 二維辨識算法流程

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

(d)分片擬合

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

1.2 脫體激波辨識

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

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

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

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

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

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

1.3 內(nèi)部激波及接觸間斷辨識

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

圖4 計算域Fig.4 Calculation domain

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

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

圖5 捕捉法密度云圖及間斷點(diǎn)集分布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 激波及接觸間斷辨識結(jié)果與流場比較Fig.7 Comparison between feature surface extraction result and flow field structure

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

2 三維實(shí)現(xiàn)

2.1 三維辨識算法流程

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

圖8 射線方向示意圖

Fig.8 Ray direction

圖9 單位球

Fig.9 Unit sphere

將單位球的球心作為光源S，球心與單位球上三角形面片的三個頂點(diǎn)ai、aj、ak的連線li、lj、lk構(gòu)成一個子區(qū)域Ap，如圖10(a)所示。若Ap內(nèi)離散點(diǎn)的個數(shù)不小于3，則對Ap內(nèi)所有離散點(diǎn)進(jìn)行平面擬合，擬合出的平面與三條射線的交點(diǎn)分別為bi,p、bj,p和bk,p；若該子區(qū)域內(nèi)離散點(diǎn)的個數(shù)小于3，則跳過該區(qū)域。所有區(qū)域擬合完成后，每條射線li與不同子區(qū)域內(nèi)的擬合平面形成若干交點(diǎn)。將各射線li上的交點(diǎn)取平均值獲得特征點(diǎn)bi，則各子區(qū)域Ap的三條射線邊界上的特征點(diǎn)bi、bj和bk構(gòu)成目標(biāo)三角形面片，如圖10(b)所示。最終，整個離散點(diǎn)集T被擬合成一系列連續(xù)的三角形面片，三角形面片頂點(diǎn)的連接關(guān)系由單位球確定?？蓪ι鲜鋈切蚊嫫瑯?gòu)成的空間曲面重新進(jìn)行網(wǎng)格劃分，獲得可用于裝配法的激波面網(wǎng)格。

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

2.2 三維測試算例

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

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

其中，

式中，i、j、k、l、m為區(qū)間[0,1]內(nèi)的隨機(jī)數(shù)。

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

圖11 橢球隨機(jī)點(diǎn)集圖12 單位球及擬合結(jié)果Fig.11 Random point set of ellipsoidFig.12 Unit sphere and the extracted surface of ellipsoid

圖13 雙橢球隨機(jī)點(diǎn)集圖14 雙橢球面擬合結(jié)果Fig.13 Random point set of double-ellipsoidFig.14 Extracted surface of double-ellipsoid

2.3 三維繞流激波辨識

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

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

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

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

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

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

3 裝配法實(shí)例

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

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

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

圖17 邊界條件及流場結(jié)構(gòu)示意圖[19]Fig.17 Boundary conditions and sketch of flow field structure[19]

(a)捕捉法計算結(jié)果

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

4 結(jié) 論

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

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

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

本文所提方法除了用于激波裝配法外，未來還可應(yīng)用于自適應(yīng)網(wǎng)格算法。根據(jù)不同的流動特征可擬合出相應(yīng)的空間曲面，進(jìn)而可以任意調(diào)整該曲面附近區(qū)域的網(wǎng)格分布，達(dá)到網(wǎng)格根據(jù)流動特征自動調(diào)整的目的。