常思源，白曉征，劉君

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

激波的探測(cè)與可視化是計(jì)算流體力學(xué)(CFD)領(lǐng)域中一個(gè)非常富有挑戰(zhàn)的問題。在某些流場(chǎng)參數(shù)的等值線云圖中，通常認(rèn)為激波間斷位于大量等值線聚集處，但對(duì)于存在接觸間斷、膨脹波和旋渦等復(fù)雜波系干擾的流場(chǎng)，該方法很容易對(duì)激波產(chǎn)生誤判，因此有必要發(fā)展更加精準(zhǔn)的激波探測(cè)(Shock Wave Detection)技術(shù)。

至今，很多學(xué)者[1-14]基于激波捕捉(Shock-Capturing)法計(jì)算出的流場(chǎng)提出了各種各樣激波探測(cè)的方法。1985年，Buning和Steger[1]首次提出將沿激波法向(近似用壓強(qiáng)梯度方向代替)馬赫數(shù)等于1的等值面視為激波面；該方法隨后被Darmofal[2]稱為基于正則馬赫數(shù)(Normal Mach Number)的激波探測(cè)法，并逐漸在流場(chǎng)可視化中被廣泛應(yīng)用；Liou等[3]在此基礎(chǔ)上介紹了3種濾波技術(shù)，以剔除辨識(shí)出的偽激波；隨后，Lovely和Haimes[4]又成功將其推廣到非定常流動(dòng)中用以探測(cè)運(yùn)動(dòng)激波。1993年，Pagendarm和Seitz[5]詳細(xì)地介紹了一種基于流場(chǎng)密度方向?qū)?shù)(Directional Derivative)來(lái)辨識(shí)激波的方法，即將密度二階方向?qū)?shù)為0的等值面視為激波陣面，并結(jié)合密度的一階方向?qū)?shù)(密度梯度在當(dāng)?shù)厮俣仁噶可系耐队?來(lái)過濾噪點(diǎn)；其后，Ma等[6]的數(shù)值算例表明基于正則馬赫數(shù)過濾噪點(diǎn)可能會(huì)獲得更好的結(jié)果。1994年，Van Rosendale[7]提出了一種針對(duì)二維非結(jié)構(gòu)網(wǎng)格的激波擬合算法，通過比較網(wǎng)格節(jié)點(diǎn)之間的密度梯度來(lái)移動(dòng)網(wǎng)格使其邊緣與激波線對(duì)齊；Ma等[6]認(rèn)為該算法中局部加權(quán)密度梯度的策略可以用于識(shí)別出激波附近的網(wǎng)格節(jié)點(diǎn)，但可能需要結(jié)合其他流場(chǎng)特征來(lái)實(shí)現(xiàn)更準(zhǔn)確的辨識(shí)。2011年，Kanamori和Suzuki[8]巧妙地將特征線理論(Theory of Characteristics)運(yùn)用于二維激波探測(cè)中，并成功推廣到非定常[8]和三維流場(chǎng)[9]計(jì)算中，相比前述方法，該算法雖然較為復(fù)雜，但不需要人為設(shè)置一些經(jīng)驗(yàn)性的閾值參數(shù)進(jìn)行濾波就可以獲得更加良好的結(jié)果。2017年，Akhlaghi等[10]通過分析數(shù)值紋影圖內(nèi)某些流動(dòng)參數(shù)的高斯分布(Gaussian Distribution)，提供了一種新的激波辨識(shí)方法，對(duì)連續(xù)流和稀薄流都比較適用。近年來(lái)還有學(xué)者[11-12]嘗試將深度學(xué)習(xí)(Deep Learning)融入大規(guī)模流場(chǎng)的激波識(shí)別中，取得了不錯(cuò)的結(jié)果。

國(guó)內(nèi)針對(duì)激波探測(cè)的研究報(bào)道相對(duì)較少，馬千里等[13]基于壓強(qiáng)梯度計(jì)算正則馬赫數(shù)，并利用激波物理特征，結(jié)合光線投射法提出了一種基于兩級(jí)采樣的多激波提取方法。2013年，Wu等[14]在其綜述中總結(jié)了以往激波辨識(shí)方法，并指出了對(duì)運(yùn)動(dòng)激波和弱激波的精確識(shí)別是非常有挑戰(zhàn)的研究方向。

整體來(lái)看，以上這些激波探測(cè)法具有一個(gè)共同的特征，即它們都是基于“當(dāng)?shù)貥?biāo)準(zhǔn)(Local Criteria)”進(jìn)行辨識(shí)的；也就是說(shuō)，激波位置的識(shí)別僅通過分析局部一個(gè)網(wǎng)格節(jié)點(diǎn)/單元(最多考慮一層相鄰網(wǎng)格)上的流動(dòng)參數(shù)來(lái)實(shí)現(xiàn)。Paciorri和Bonfiglioli[15]認(rèn)為這些探測(cè)技術(shù)的當(dāng)?shù)匦再|(zhì)正是其局限性和缺陷的根源。一方面，這些方法往往將探測(cè)到的一系列散點(diǎn)[6]或小線段[8-9]視為激波線/面，但由于數(shù)值耗散和振蕩，這些激波線/面處經(jīng)常會(huì)出現(xiàn)噪聲、偽激波分支或空隙[14]。另一方面，所有這些算法都不能判斷出流場(chǎng)中是否存在激波干擾點(diǎn)，即不能清晰地給出流場(chǎng)中激波干擾的模式，如λ激波、異側(cè)激波相交、激波-壁面反射等。

近年來(lái)，有一批學(xué)者[16-19]重新開啟了激波裝配(Shock-Fitting)法的研究與應(yīng)用，相比當(dāng)前主流的激波捕捉法，該方法必須顯式地處理激波，因此在開始計(jì)算前，通常都需要從激波捕捉流場(chǎng)獲取相關(guān)信息，如合理的初始流動(dòng)參數(shù)、激波的近似位置、激波的干擾結(jié)構(gòu)等。然而，由于上述激波探測(cè)法的缺陷，其很難直接用于激波裝配的初始化中[15]，最終造成這些激波裝配算法難免需要依據(jù)先驗(yàn)知識(shí)進(jìn)行人工干預(yù)來(lái)設(shè)置初始激波及激波干擾點(diǎn)的位置。

針對(duì)上述問題，Paciorri和Bonfiglioli[15]于2012年介紹了一種二維激波模式識(shí)別技術(shù)，首先根據(jù)傳統(tǒng)基于當(dāng)?shù)孛芏榷A方向?qū)?shù)的方法[6]提取出激波點(diǎn)云，隨后利用霍夫變換(Hough Transform)搜尋這些散點(diǎn)中的特征直線并去除噪點(diǎn)，再采用最小二乘(Least-Squares)法擬合得到各條激波曲線，最后通過一套模糊邏輯(Fuzzy Logic)算法識(shí)別出激波干擾的模式。然而，該方法實(shí)施起來(lái)比較繁瑣，對(duì)于流場(chǎng)中長(zhǎng)度較短的激波分支容易產(chǎn)生誤判，且很難適用于非定常流場(chǎng)。

本文提出了一種面向全局(Global)的二維激波模式識(shí)別算法，對(duì)探測(cè)出的激波帶網(wǎng)格單元進(jìn)一步處理，通過聚類分析(Cluster Analysis)方法將激波單元?jiǎng)澐殖稍S多簇；分析各個(gè)簇的相鄰關(guān)系就可以清晰地判斷出激波干擾的模式，最終基于Bézier曲線擬合算法優(yōu)化的激波線和數(shù)值捕捉的激波相比在形狀和位置上都較為吻合。該算法不僅邏輯簡(jiǎn)單，且不受網(wǎng)格類型的限制，同時(shí)在非定常流動(dòng)中也具有較高的可靠性。

1 二維激波模式識(shí)別算法

結(jié)合一個(gè)定常無(wú)黏激波反射算例(幾何模型和流場(chǎng)如圖1所示，其中p表示壓強(qiáng)，下標(biāo)∞表示來(lái)流狀態(tài))，從3個(gè)主要方面介紹了該二維激波模式識(shí)別算法的步驟和特點(diǎn)。在該算例中，來(lái)流馬赫數(shù)Ma∞=1.5，楔角θ為10°，入口和出口高度分別為2.17L和2L。入射斜激波S1在上壁面發(fā)生馬赫反射產(chǎn)生了強(qiáng)度較強(qiáng)的馬赫干擾S2，而反射激波S3經(jīng)膨脹波異側(cè)干擾后強(qiáng)度下降，最后在下壁面產(chǎn)生了正規(guī)反射，次反射激波S4與出口邊界相交。此外，雖然該流場(chǎng)是基于均勻的非結(jié)構(gòu)三角形網(wǎng)格(網(wǎng)格尺寸Δ=0.015L)計(jì)算的，但是該識(shí)別算法同樣適用于結(jié)構(gòu)網(wǎng)格、混合網(wǎng)格和非均勻網(wǎng)格等，在2.2節(jié)有相關(guān)描述。

1.1 辨識(shí)激波單元

本文算法是針對(duì)網(wǎng)格單元而非網(wǎng)格散點(diǎn)實(shí)現(xiàn)的，因此首先需要從圖1所示流場(chǎng)提取出表征激波位置的網(wǎng)格單元，即圖2所示激波單元。考慮到引言所述若干激波探測(cè)方法的適用性和魯棒性，本文基于流場(chǎng)單元格心的壓強(qiáng)梯度，設(shè)計(jì)了一套針對(duì)激波的辨識(shí)準(zhǔn)則，具體步驟為

圖2 定常激波反射：初始辨識(shí)的激波單元

(1)

由于考慮了流向，δp的正、負(fù)號(hào)分別對(duì)應(yīng)流體處于壓縮和膨脹的狀態(tài)，而激波處通常表現(xiàn)出很強(qiáng)的壓縮性，因此首先要過濾掉δp<0的單元。

2) 由于數(shù)值計(jì)算難免會(huì)存在耗散和誤差，會(huì)造成一些非激波區(qū)域的δp值為一個(gè)較小的正量，故過濾掉滿足式(2)的單元：

δp<ξ1δpmax

(2)

式中：δpmax為流向壓強(qiáng)梯度的全場(chǎng)最大值；ξ1為全局濾波系數(shù)，當(dāng)流場(chǎng)中激波的強(qiáng)度比較懸殊時(shí)，該系數(shù)過大可能會(huì)濾掉某些弱激波而導(dǎo)致部分激波帶缺失，因此經(jīng)大量數(shù)值試驗(yàn)的調(diào)試，取ξ1=0.01。流場(chǎng)中除了激波間斷外可能還存在接觸間斷，由于數(shù)值捕捉出的接觸間斷兩側(cè)壓強(qiáng)近似相等，因此基于式(2)流向壓強(qiáng)梯度閾值法可以有效地過濾掉接觸間斷。

3) 然而幾乎不可能只通過設(shè)置滿意的ξ1，既能完整地辨識(shí)多激波又能去除全部噪聲。此外，因全局濾波強(qiáng)度不大，經(jīng)第2步辨識(shí)的激波帶可能會(huì)出現(xiàn)一些偽激波單元，為了在一定程度上提高辨識(shí)方法的普適性，考慮局部特性對(duì)滿足式(3)的單元進(jìn)行二次去噪：

(3)

定常激波反射算例經(jīng)上述方法辨識(shí)出的具有一定厚度的“激波帶”如圖3所示，受激波強(qiáng)弱影響，激波帶的厚度并不均勻，通常有2～5層激波單元。由于數(shù)值誤差的存在，激波帶往往存在很多“毛刺”和“空穴”，為了使后續(xù)聚類分析更加準(zhǔn)確，從兩方面對(duì)原始辨識(shí)的激波帶進(jìn)行了光順優(yōu)化。首先一方面標(biāo)記“空穴單元”為激波單元，另一方面剔除“毛刺單元”。所謂空穴單元是指其所有節(jié)點(diǎn)均落在激波帶上的初始非激波單元(圖3中綠色單元)，而毛刺單元是指共面相鄰激波單元的個(gè)數(shù)小于2的初始激波單元(圖3中藍(lán)色單元)。

圖3 定常激波反射：激波帶的優(yōu)化

1.2 激波單元的聚類分析

聚類分析是一種十分重要的統(tǒng)計(jì)數(shù)據(jù)分析方法，在數(shù)據(jù)挖掘、模式識(shí)別、圖像分析和機(jī)器學(xué)習(xí)等許多領(lǐng)域受到廣泛應(yīng)用；它的目標(biāo)是將數(shù)據(jù)集分成許多簇，使得同一簇內(nèi)的數(shù)據(jù)點(diǎn)相似度盡可能大，而不同簇間的數(shù)據(jù)點(diǎn)相似度盡可能小[20]。本文采用經(jīng)典的K-means聚類(K-Means Clustering)[21]算法對(duì)優(yōu)化后的激波帶進(jìn)行聚類處理，下面舉例簡(jiǎn)要介紹該算法對(duì)圖4(a)數(shù)據(jù)點(diǎn)的聚類流程：

圖4 K-means聚類算法示意圖

1) 選擇K個(gè)初始“聚類中心”。在圖4(b)中，初始化選擇了兩個(gè)聚類中心，即圖中紅、藍(lán)“×”形點(diǎn)。

2) 對(duì)任意一個(gè)數(shù)據(jù)點(diǎn)，求其到K個(gè)聚類中心的“距離”(取歐氏距離)，將該數(shù)據(jù)點(diǎn)歸到距離其最近的聚類中心所在的“簇”。如圖4(c)所示，所有數(shù)據(jù)點(diǎn)在經(jīng)過第一輪遍歷后被分為紅、藍(lán)標(biāo)識(shí)的兩個(gè)簇。

3) 利用均值等方法更新K個(gè)簇的中心值。此時(shí)對(duì)當(dāng)前標(biāo)記為紅色和藍(lán)色的數(shù)據(jù)點(diǎn)分別求出其新的質(zhì)心，作為新的聚類中心，如圖4(d)所示，紅色和藍(lán)色聚類中心的位置發(fā)生了變動(dòng)。

4) 重復(fù)步驟2)～3)，若所有聚類中心在迭代更新后位置都保持不變，則迭代結(jié)束；否則繼續(xù)迭代。最終得到的兩個(gè)簇及其聚類中心如圖4(f)所示。

K-means聚類算法操作簡(jiǎn)單，易于實(shí)現(xiàn)，且十分高效，一般僅需5～10次迭代即可收斂。但是研究表明，該算法對(duì)初始聚類中心的位置選取比較敏感，不同的隨機(jī)初始聚類中心得到的聚類結(jié)果可能差異顯著。因此，在對(duì)激波單元聚類時(shí)，合理的初始聚類中心位置顯得至關(guān)重要，其既不能偏離激波帶過遠(yuǎn)，也不能分布的過密或過稀。下面介紹下針對(duì)激波帶初始聚類中心的選取方法。

首先任意選取一激波單元作為搜索起點(diǎn)，以該單元格心為圓心，n倍當(dāng)?shù)鼐W(wǎng)格尺寸為半徑r作“搜索圓”(經(jīng)大量算例評(píng)估，一般n取6～7即可獲得合理的結(jié)果)，如圖5所示；隨后將格心點(diǎn)落在該圓區(qū)域內(nèi)且未被搜尋到的激波單元?dú)w為一簇，計(jì)算該簇中所有激波單元格心坐標(biāo)的平均值，作為新的搜索圓圓心坐標(biāo)；反復(fù)迭代直到所有激波單元均被搜尋到，則所有搜索圓的圓心即為初始聚類中心。

圖5 定常激波反射：確定初始聚類中心

取所有激波單元的格心作為數(shù)據(jù)點(diǎn)集，對(duì)所有激波單元進(jìn)行K-means聚類；經(jīng)過若干次迭代后最終所有激波單元被劃分到不同的簇，如圖6所示，為方便顯示，緊鄰的簇用不同顏色加以區(qū)分。

圖6 定常激波反射：K-means聚類結(jié)果

為了后續(xù)便于分析激波拓?fù)浣Y(jié)構(gòu)及其干擾模式，根據(jù)每個(gè)簇的相鄰特征將簇進(jìn)一步分成4類：

1) 終止簇。若流場(chǎng)內(nèi)部某簇僅有1個(gè)相鄰簇，則該簇為終止簇；其一般位于流場(chǎng)內(nèi)部激波起始/終止處，注意該簇不與邊界相鄰。

2) 普通簇。若流場(chǎng)內(nèi)部某簇僅有2個(gè)相鄰簇，則該簇為普通簇；其一般呈“串狀”普遍分布于激波帶中，是數(shù)量最多的簇。

3) 分叉簇。若流場(chǎng)內(nèi)部某簇有3個(gè)及以上相鄰簇，則該簇為分叉簇；其一般分布于厚度較大的激波帶處，如三波點(diǎn)、四波點(diǎn)和弱激波處等。

4) 邊界簇。若某簇直接與流場(chǎng)邊界相鄰，則該簇為邊界簇；其通常存在于激波-壁面反射處、壁面拐角處以及流場(chǎng)進(jìn)出口邊界處等。

其中，終止簇、分叉簇和邊界簇統(tǒng)稱關(guān)鍵簇，其在流場(chǎng)中個(gè)數(shù)雖少，卻往往位于激波拓?fù)浠蛐螤畎l(fā)生變化的關(guān)鍵位置處。

1.3 激波的模式識(shí)別與擬合

如何自動(dòng)準(zhǔn)確地辨識(shí)出流場(chǎng)中激波干擾點(diǎn)的位置是一個(gè)難題。下面通過對(duì)激波帶K-means聚類得到的各個(gè)簇進(jìn)行“近鄰分析”，給出了一種激波模式識(shí)別和激波線擬合的策略：

1) 判斷是否存在需要進(jìn)行合并的簇。具體分兩方面，若某分叉簇與其他分叉簇相鄰，則將這些分叉簇合并；若某邊界簇與其他邊界簇(注意它們必須對(duì)應(yīng)同一個(gè)邊界)相鄰，則將這些邊界簇合并。

2) 若進(jìn)行了簇合并的操作，需要重新對(duì)簇進(jìn)行分類，并計(jì)算新簇的聚類中心。如圖6(a)無(wú)需進(jìn)行簇合并，而圖6(b)經(jīng)簇合并后變成圖7，即原4個(gè)相鄰分叉簇合并后形成1個(gè)新的普通簇，一同緊靠下壁面的2個(gè)相鄰邊界簇合并后形成1個(gè)新的邊界簇。

圖7 定常激波反射：簇合并后

3) 識(shí)別關(guān)鍵激波點(diǎn)。規(guī)定若某分叉簇與周邊3個(gè)簇(見圖6(a))或4個(gè)簇(見圖23)相鄰，則該分叉簇對(duì)應(yīng)的聚類中心即為三波點(diǎn)或四波點(diǎn)；若某分叉簇僅與1個(gè)簇相鄰，則其聚類中心即為激波終止/起始點(diǎn)(見圖25)；若某邊界簇與2個(gè)簇相鄰(圖7下邊界處)，則取該邊界簇對(duì)應(yīng)的邊界中心點(diǎn)作為激波-壁面反射點(diǎn)；同理，若某邊界簇僅與1個(gè)簇相鄰(圖7右邊界處)，則邊界中心點(diǎn)即為激波-邊界干擾點(diǎn)。通過以上近鄰分析，便可識(shí)別出激波-激波干擾或激波-邊界干擾等模式。

4) 記錄各條激波所包含的簇。從第3)步得到的某一關(guān)鍵簇開始搜索，依次記錄下該分支上兩兩相鄰的所有普通簇，直至搜索到另一關(guān)鍵簇停止，此時(shí)該一系列簇即對(duì)應(yīng)一條激波；繼續(xù)搜索確定下一激波分支上的簇，直到遍歷完流場(chǎng)中所有的簇。

5) 擬合激波線。依次連接激波包含的一系列簇所對(duì)應(yīng)的聚類中心，便可得到各條激波線；然而此時(shí)得到的激波線可能比較曲折(圖7中的白色實(shí)線)，為此，本文進(jìn)一步采用著名的Bézier曲線擬合算法[22]獲取更加光滑的激波線，下面對(duì)其作簡(jiǎn)要介紹。

當(dāng)給定一系列有序散點(diǎn)P0、P1、…、Pn時(shí)，則其n階Bézier曲線的一般參數(shù)化形式為

(4)

如4個(gè)點(diǎn)對(duì)應(yīng)的三階Bézier曲線表達(dá)式為

B3(λ)=P0(1-λ)3+3P1λ(1-λ)2+

3P2λ2(1-λ)+P3λ3λ∈[0,1]

(5)

式中：Pi稱為Bézier曲線的控制點(diǎn)，且該曲線開始于點(diǎn)P0(對(duì)應(yīng)λ=0)并結(jié)束于點(diǎn)Pn(對(duì)應(yīng)λ=1)。值得注意的是，Bézier曲線是有方向性的，控制點(diǎn)的次序不同會(huì)擬合得到不同的曲線；且Bézier曲線是直線的充分必要條件是所有控制點(diǎn)共線。

因此，根據(jù)Bézier曲線的特點(diǎn)，若第4)步記錄的某條激波包含n+1個(gè)簇，則將兩個(gè)關(guān)鍵簇的聚類中心分別作為起始點(diǎn)P0和終止點(diǎn)Pn，而將其余依次緊鄰的普通簇的聚類中心作為控制點(diǎn)P1～Pn-1，進(jìn)行n階Bézier曲線擬合便得到更加光滑的激波線，如圖8所示。

圖8 定常激波反射：擬合的激波線

圖9對(duì)比了最終擬合的激波線與流場(chǎng)壓強(qiáng)云圖，4條激波線、三波點(diǎn)、激波-壁面反射點(diǎn)和激波-邊界干擾點(diǎn)的位置都比較準(zhǔn)確，顯示出該算法具有良好的有效性和準(zhǔn)確性。

圖9 定常激波反射：激波線擬合結(jié)果與云圖對(duì)比

1.4 算法小結(jié)

該二維激波模式識(shí)別算法的總體流程如圖10所示，主要包含激波辨識(shí)、聚類分析和識(shí)別擬合3方面內(nèi)容；此三者層層遞進(jìn)，激波辨識(shí)是基礎(chǔ)，聚類分析是核心，而識(shí)別擬合是關(guān)鍵。

圖10 算法流程圖

從激波捕捉流場(chǎng)出發(fā)，通過流向壓強(qiáng)梯度閾值法進(jìn)行二級(jí)濾波辨識(shí)提取出具有一定厚度的激波帶，到最終準(zhǔn)確識(shí)別出激波-激波干擾點(diǎn)和激波-邊界干擾點(diǎn)等關(guān)鍵位置，且擬合出的各條激波線的形狀和位置都具有較高的精度?？傊撍惴ǖ淖畲髢?yōu)勢(shì)是基本不需要人工參與，且算法簡(jiǎn)單高效，可以準(zhǔn)確實(shí)現(xiàn)從“激波帶”到“激波線”的擬合識(shí)別。

2 算例驗(yàn)證

通過4個(gè)超聲速流動(dòng)算例從多方面考察了該算法的準(zhǔn)確性、適用性和可靠性。以下激波捕捉流場(chǎng)均基于格心型有限體積法[23]計(jì)算，空間對(duì)流項(xiàng)通量采用van Leer格式，時(shí)間推進(jìn)采用顯式4步Runge-Kutta格式，保證在光滑流場(chǎng)中時(shí)空均具有二階精度，全場(chǎng)均不考慮黏性。

2.1 定常斜激波

首先，采用兩種條件下的定常斜激波算例，驗(yàn)證本算法最終擬合出的激波線的準(zhǔn)確性。計(jì)算區(qū)域?yàn)閤∈[0,L]，y∈[0,0.9L]，在x=0.1L處存在楔角，均勻超聲速來(lái)流經(jīng)過楔面會(huì)產(chǎn)生一道筆直的斜激波，上邊界為超聲速出口，激波不會(huì)反射。為了產(chǎn)生相對(duì)較弱和較強(qiáng)的斜激波，楔面角θw分別取3°和20°，對(duì)應(yīng)的來(lái)流馬赫數(shù)Ma∞分別為1.469和1.959。計(jì)算網(wǎng)格尺寸均為0.02L，全場(chǎng)分別有5 257個(gè)和4 472個(gè)三角形網(wǎng)格。

兩種狀態(tài)下的斜激波經(jīng)辨識(shí)和聚類分析的結(jié)果如圖11所示，弱激波辨識(shí)出的激波帶較寬，最終一共得到25個(gè)聚類中心；而強(qiáng)激波辨識(shí)出的激波帶更加銳利，最終一共得到了22個(gè)聚類中心。將這些聚類中心作為控制點(diǎn)按順序代入式(4)擬合得到Bézier曲線，并與相應(yīng)的理論值進(jìn)行對(duì)比(圖12)，可見擬合值具有較高的位置精度。表1進(jìn)一步對(duì)比了激波角，可見無(wú)論是弱激波還是強(qiáng)激波，最終擬合出的激波線的激波角度和理論值的相對(duì)偏差都不到1%。

圖11 定常斜激波：聚類結(jié)果

圖12 定常斜激波：激波線擬合值與理論值對(duì)比

表1 定常斜激波：激波角擬合值和理論值對(duì)比

2.2 定常激波反射

仍考慮圖1所示的定常激波反射流動(dòng)，測(cè)試本算法對(duì)不同類型網(wǎng)格的可靠性。首先分別基于全場(chǎng)均勻的四邊形結(jié)構(gòu)網(wǎng)格和三角形/四邊形混合網(wǎng)格重新進(jìn)行數(shù)值計(jì)算，平均網(wǎng)格尺寸均為0.015L。 圖13給出了這兩種類型網(wǎng)格下局部流場(chǎng)壓強(qiáng)等值線云圖(圖例與圖1一致)。

圖13 定常激波反射：不同種類網(wǎng)格下的壓強(qiáng)云圖

根據(jù)1.1節(jié)所述的激波探測(cè)方法，采用相同的濾波系數(shù)，分別獲得滿足條件的激波帶。在對(duì)激波單元進(jìn)行K-means聚類分析時(shí)，結(jié)構(gòu)網(wǎng)格算例全場(chǎng)一共確定了64個(gè)初始聚類中心，而混合網(wǎng)格算例共搜尋到了82個(gè)，經(jīng)若干步后得到收斂的簇及其聚類中心。此時(shí)，根據(jù)簇的類別和相鄰簇的個(gè)數(shù)判斷某些簇是否需要進(jìn)行合并，圖14和圖15分別給出了兩種網(wǎng)格下壁面馬赫反射處和正規(guī)反射處聚類合并后的結(jié)果。

圖14(a)和圖15(a)中各存在一個(gè)分叉簇，因其與周圍3個(gè)激波帶分支中的3個(gè)普通簇相鄰，可以判定該分叉簇的聚類中心位于流場(chǎng)三波點(diǎn)附近。此外，反射激波S3在受到膨脹波干擾后，強(qiáng)度顯著減弱；觀察圖14(b)和圖15(b)壁面激波反射處的激波帶，其厚度要明顯大于三波點(diǎn)附近激波帶的厚度，相應(yīng)地，各簇包含的單元個(gè)數(shù)也明顯增加，但并沒有存在錯(cuò)誤的分叉簇，且下壁面和出口邊界處的兩個(gè)邊界簇都被準(zhǔn)確地識(shí)別出來(lái)。值得注意的是，由于激波-壁面反射點(diǎn)比較靠近出口邊界，圖14(b)中存在兩個(gè)邊界簇直接相鄰的情況，但由于它們鄰近的邊界不同，因此在簇合并操作中并未錯(cuò)誤地將它們合并。

圖14 定常激波反射：基于結(jié)構(gòu)網(wǎng)格的聚類結(jié)果

圖15 定常激波反射：基于混合網(wǎng)格的聚類結(jié)果

以上采用的都是均勻網(wǎng)格，下面進(jìn)一步考察了本算法在非均勻網(wǎng)格中的適用性。如圖16所示，在非結(jié)構(gòu)網(wǎng)格的基礎(chǔ)上，沿入射激波的法向進(jìn)行加密，第一層網(wǎng)格高度為0.002L，向激波兩側(cè)各推進(jìn)5層四邊形網(wǎng)格，網(wǎng)格高度增長(zhǎng)率取1.2，隨后將激波右側(cè)四邊形網(wǎng)格沿對(duì)角線剖分成三角形網(wǎng)格。圖17顯示了對(duì)辨識(shí)出的激波帶進(jìn)行聚類分析后的結(jié)果，雖然激波帶厚度差異明顯，但并不影響分叉簇的位置識(shí)別。

圖16 定常激波反射：非均勻網(wǎng)格下的壓強(qiáng)云圖

圖17 定常激波反射：基于非均勻網(wǎng)格的聚類結(jié)果

最后對(duì)比了以上3種網(wǎng)格下最終的擬合結(jié)果(圖18)，4條擬合激波線的形狀和位置都與數(shù)值捕捉的激波吻合得很好；從局部放大圖看出，三者擬合出的三波點(diǎn)和邊界干擾點(diǎn)位置稍有不同，這是由于流場(chǎng)數(shù)值誤差和初始激波帶辨識(shí)差異造成的，難以完全根除，但偏差都在一個(gè)網(wǎng)格尺度左右，仍顯示出較好的一致性。

2.3 同側(cè)和異側(cè)激波干擾

對(duì)波系較為復(fù)雜的多激波干擾問題進(jìn)行了數(shù)值驗(yàn)證，幾何模型和流場(chǎng)壓強(qiáng)等值線如圖19所示，全場(chǎng)采用均勻的非結(jié)構(gòu)Delaunay三角形網(wǎng)格離散，平均網(wǎng)格尺寸Δ1=0.015L。Ma∞=3.0的高超聲速均勻氣流在經(jīng)過收縮管道時(shí)，經(jīng)上壁面20°楔角壓縮產(chǎn)生了一道斜激波S1，與此同時(shí)，先后經(jīng)下壁面10°和20°楔角二次壓縮產(chǎn)生了兩道斜激波S2和S3，隨后S2和S3發(fā)生了同側(cè)正規(guī)激波相交并匯聚成了一道新激波S4，最終激波S1與S4發(fā)生了異側(cè)正規(guī)相交并產(chǎn)生了兩條透射激波S5和S6。全場(chǎng)一共產(chǎn)生了六道激波，并存在一個(gè)三波點(diǎn)和一個(gè)四波點(diǎn)結(jié)構(gòu)。

圖19 激波干擾：幾何模型與壓強(qiáng)云圖

根據(jù)流向壓強(qiáng)梯度過濾得到的激波帶分布如圖20所示，在兩個(gè)激波干擾處聚集了大量的激波單元。通過進(jìn)一步增補(bǔ)原激波帶外緣兩側(cè)的空穴單元，優(yōu)化后的激波帶更加光順。隨后對(duì)優(yōu)化的激波帶進(jìn)行K-means聚類，結(jié)果如圖21(a)所示，在每條激波帶分支中都分布著大量呈串狀依次排列的普通簇，而在分支交匯處卻聚集著若干分叉簇，這些分叉簇都至少與周圍3個(gè)簇直接相鄰，因此需要進(jìn)行簇合并操作。

圖20 激波干擾：辨識(shí)出的激波帶

在將某些簇合并后(圖21(b))，清晰地存在兩個(gè)較大的分叉簇，分別與周邊3和4個(gè)普通簇相鄰，因此其聚類中心可以被視為三波點(diǎn)和四波點(diǎn)位置。最終分別將各個(gè)分支中的聚類中心進(jìn)行擬合得到了所有6條激波線，結(jié)果如圖22所示。擬合的激波線、激波相交點(diǎn)均落在數(shù)值捕捉激波的過渡區(qū)內(nèi)，顯示出該算法良好的精度。

圖22 激波干擾：擬合的激波線

實(shí)際上，激波(尤其是激波-激波干擾點(diǎn))的探測(cè)精度在很大程度上取決于流場(chǎng)的分辨率。也就是說(shuō)，捕捉得到的激波過渡區(qū)寬度越小，越有利于精確探測(cè)到激波及其干擾點(diǎn)的位置。下面考察了網(wǎng)格的疏密程度對(duì)激波探測(cè)的影響。

在其他相同計(jì)算條件下，僅僅改變流場(chǎng)網(wǎng)格的尺寸，分別設(shè)計(jì)疏(Δ2=0.02L)和密(Δ3=0.01L)兩套非結(jié)構(gòu)網(wǎng)格，重新進(jìn)行流場(chǎng)計(jì)算。經(jīng)激波辨識(shí)、聚類分析和簇合并后的結(jié)果如圖23所示；網(wǎng)格越稀，劃分的簇越少，且聚類中心兩兩相距越遠(yuǎn)。圖24對(duì)比了3種網(wǎng)格尺寸下最終擬合出的激波線，可以發(fā)現(xiàn)，三者大部分基本完全吻合，但激波干擾點(diǎn)處差異較大。對(duì)于同側(cè)激波干擾結(jié)構(gòu)，由于三道激波的斜率都比較接近，因此當(dāng)網(wǎng)格較稀時(shí)，數(shù)值捕捉激波的過渡區(qū)較寬，會(huì)造成兩條入射激波帶提前交匯(對(duì)比圖23中兩圖可見)，分叉簇的位置會(huì)發(fā)生較大改變，最終造成了辨識(shí)出的三波點(diǎn)位置出現(xiàn)較大差異。另一方面，對(duì)于異側(cè)激波干擾產(chǎn)生的四波點(diǎn)結(jié)構(gòu)，雖然稀網(wǎng)格對(duì)應(yīng)的分叉簇也較大，但其聚類中心的位置和密網(wǎng)格的結(jié)果偏差并不顯著。

圖23 激波干擾：不同網(wǎng)格尺寸下的聚類結(jié)果

圖24 激波干擾：3種網(wǎng)格尺寸下擬合的激波線對(duì)比

該算例表明，本文算法對(duì)同側(cè)激波干擾點(diǎn)的識(shí)別精度網(wǎng)格依賴性相對(duì)較強(qiáng)，即網(wǎng)格疏密對(duì)結(jié)果的影響較大；而對(duì)于異側(cè)激波干擾點(diǎn)，包括上文激波反射算例中的三波點(diǎn)結(jié)構(gòu)，其識(shí)別精度受網(wǎng)格尺度的影響較?。欢鴮?duì)于干擾點(diǎn)以外的激波線，擬合精度受網(wǎng)格尺度的影響最小，始終吻合得比較好。

2.4 前向臺(tái)階非定常流動(dòng)

采用著名的前向臺(tái)階超聲速繞流算例[24]測(cè)試該算法在復(fù)雜非定常流動(dòng)中的有效性。均勻自由來(lái)流馬赫數(shù)Ma∞=3.0，且初始流場(chǎng)取來(lái)流參數(shù)，通道入口和出口處的高度分別為L(zhǎng)和0.8L，臺(tái)階高度及其上表面長(zhǎng)度分別為0.2L和2.4L。全場(chǎng)采用均勻非結(jié)構(gòu)三角形網(wǎng)格劃分，網(wǎng)格尺度為0.01L。本算例中的時(shí)間均為無(wú)量綱量。

在流場(chǎng)的發(fā)展歷程中，臺(tái)階前首先出現(xiàn)一道向左運(yùn)動(dòng)的弓形激波，因其強(qiáng)度較大，當(dāng)受到上壁面干擾時(shí)會(huì)發(fā)生反射，并從運(yùn)動(dòng)正規(guī)反射逐漸過渡到馬赫反射；而反射激波又從下壁面反射回上壁面，最終在管道壁面形成了4次激波-壁面干擾(Shock-Wall Interaction，SWI)。此外，氣流在臺(tái)階拐角處發(fā)生過度膨脹，過膨脹氣流會(huì)與下壁面相撞形成一道相對(duì)較弱的孤立斜激波，該激波最終與第一道反射激波發(fā)生異側(cè)正規(guī)相交，相交后的弱激波又與第二道反射激波匯聚形成同側(cè)激波-激波干擾(Shock-Shock Interaction，SSI)結(jié)構(gòu)。

圖25依次給出了5個(gè)關(guān)鍵時(shí)刻的辨識(shí)激波帶的聚類結(jié)果，通過二次濾波，在辨識(shí)出運(yùn)動(dòng)強(qiáng)激波的同時(shí)，那道較弱的孤立斜激波也被很好地提取出來(lái)；此外，圖25(d)和25(e)中并沒有錯(cuò)誤地將三波點(diǎn)處的接觸間斷提取出來(lái)。隨后，通過聚類分析和簇合并操作，可以準(zhǔn)確地定位分叉簇與邊界簇的位置，從而清晰地識(shí)別出各個(gè)時(shí)刻流場(chǎng)激波的總個(gè)數(shù)及相互干擾的模式(如SWI或SSI)。

圖25 前向臺(tái)階繞流：不同時(shí)刻激波帶的聚類結(jié)果

以兩時(shí)刻為例，圖26定性比較了最終擬合出的激波線與流場(chǎng)壓強(qiáng)云圖。可以發(fā)現(xiàn)，不僅各條近似直線的反射激波位置準(zhǔn)確，而且臺(tái)階前的大曲率弓形激波也與捕捉激波吻合得很好，反映了該算法對(duì)各種形狀激波的擬合都具有較好的可靠性。在得到不同時(shí)刻的擬合激波線后，可以方便地將其在同一張圖中顯示出來(lái)，進(jìn)而有利于分析激波的演化過程，深化對(duì)運(yùn)動(dòng)激波干擾的認(rèn)識(shí)。從圖27可以看出，相比激波-上壁面干擾點(diǎn)，臺(tái)階后的激波-下壁面干擾點(diǎn)向左移動(dòng)的速度更快；且弓形激波在上壁面反射形成的馬赫干擾的長(zhǎng)度在迅速增大；而在t=1.5時(shí)刻后，弓形激波的位置變化比較緩慢，三波點(diǎn)近似沿弓形激波向下滑動(dòng)。

圖26 前向臺(tái)階繞流：擬合的激波線與壓強(qiáng)云圖對(duì)比

圖27 前向臺(tái)階繞流：4個(gè)時(shí)刻擬合的激波線

3 結(jié) 論

1) 該算法在傳統(tǒng)基于流場(chǎng)當(dāng)?shù)貐?shù)進(jìn)行激波探測(cè)的方法基礎(chǔ)上，對(duì)提取的激波單元進(jìn)一步處理；先后結(jié)合K-means聚類算法和近鄰分析成功實(shí)現(xiàn)了激波干擾點(diǎn)的準(zhǔn)確定位；最后采用Bézier曲線算法擬合得到各條質(zhì)量較高的激波線，數(shù)值試驗(yàn)表明激波線的大部分位置與形狀受流場(chǎng)網(wǎng)格疏密影響不大。

2) 對(duì)于同側(cè)激波干擾，干擾點(diǎn)的位置受流場(chǎng)本身分辨率的影響較大，密網(wǎng)格下一般會(huì)獲得更為精準(zhǔn)的結(jié)果；而對(duì)于異側(cè)激波干擾，其干擾點(diǎn)的識(shí)別結(jié)果通常比較準(zhǔn)確，且網(wǎng)格依賴性不強(qiáng)。

3) 該算法自動(dòng)化程度較高，基本不需要人工干預(yù)，且適用于任意網(wǎng)格類型，對(duì)復(fù)雜非定常流場(chǎng)中的運(yùn)動(dòng)激波辨識(shí)也具有良好的可靠性和準(zhǔn)確性。

總之，該面向全局的激波探測(cè)算法可以有效識(shí)別出更為精確的激波拓?fù)浣Y(jié)構(gòu)，未來(lái)可以用于CFD中某些對(duì)激波位置有特殊需求的技術(shù)中，例如激波裝配[16-19]、網(wǎng)格自適應(yīng)[25]、流場(chǎng)特征可視化[26]等。由于三維激波相交干擾結(jié)構(gòu)比較復(fù)雜，將該方法推廣到三維是今后努力的方向。