一種用于間斷裝配的流場(chǎng)結(jié)構(gòu)辨識(shí)算法

劉 君,陳澤棟,陳 潔,鄒東陽(yáng)

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

0 引 言

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

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

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

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

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

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

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

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

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

1 流場(chǎng)結(jié)構(gòu)辨識(shí)算法

1.1 二維辨識(shí)算法流程

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

(d)分片擬合

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

1.2 脫體激波辨識(shí)

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

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

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

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

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

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

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

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

圖4 計(jì)算域Fig.4 Calculation domain

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

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

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

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

圖6 激波和接觸間斷辨識(shí)過(guò)程
Fig.6 Shock wave and contact discontinuity surface extraction process

圖7 激波及接觸間斷辨識(shí)結(jié)果與流場(chǎng)比較Fig.7 Comparison between feature surface extraction result and flow field structure

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

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

2.1 三維辨識(shí)算法流程

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

圖8 射線方向示意圖

Fig.8 Ray direction

圖9 單位球

Fig.9 Unit sphere

將單位球的球心作為光源S，球心與單位球上三角形面片的三個(gè)頂點(diǎn)ai、aj、ak的連線li、lj、lk構(gòu)成一個(gè)子區(qū)域Ap，如圖10(a)所示。若Ap內(nèi)離散點(diǎn)的個(gè)數(shù)不小于3，則對(duì)Ap內(nèi)所有離散點(diǎn)進(jìn)行平面擬合，擬合出的平面與三條射線的交點(diǎn)分別為bi,p、bj,p和bk,p；若該子區(qū)域內(nèi)離散點(diǎn)的個(gè)數(shù)小于3，則跳過(guò)該區(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)所示。最終，整個(gè)離散點(diǎn)集T被擬合成一系列連續(xù)的三角形面片，三角形面片頂點(diǎn)的連接關(guān)系由單位球確定?？蓪?duì)上述三角形面片構(gòu)成的空間曲面重新進(jìn)行網(wǎng)格劃分，獲得可用于裝配法的激波面網(wǎng)格。

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

2.2 三維測(cè)試算例

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

在橢球面附近隨機(jī)生成100 000個(gè)點(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個(gè)節(jié)點(diǎn)、1468個(gè)三角形單元(如圖9)，球心位置坐標(biāo)為(0,0,0)。單位球位置及橢球面辨識(shí)結(jié)果如圖12所示。從結(jié)果可以看出，辨識(shí)出的曲面光滑完整且形狀與離散點(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 三維繞流激波辨識(shí)

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

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

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

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

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

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

3 裝配法實(shí)例

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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