一種基于聚類分析的二維激波模式識別算法

常思源，白曉征，劉君

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

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

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

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

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

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

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

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

1 二維激波模式識別算法

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

1.1 辨識激波單元

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

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

(1)

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

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

δp<ξ1δpmax

(2)

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

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

(3)

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

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

1.2 激波單元的聚類分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 激波的模式識別與擬合

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

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

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

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

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

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

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

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

(4)

如4個點對應(yīng)的三階Bézier曲線表達式為

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

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

(5)

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

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

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

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

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

1.4 算法小結(jié)

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

圖10 算法流程圖

從激波捕捉流場出發(fā)，通過流向壓強梯度閾值法進行二級濾波辨識提取出具有一定厚度的激波帶，到最終準(zhǔn)確識別出激波-激波干擾點和激波-邊界干擾點等關(guān)鍵位置，且擬合出的各條激波線的形狀和位置都具有較高的精度。總之，該算法的最大優(yōu)勢是基本不需要人工參與，且算法簡單高效，可以準(zhǔn)確實現(xiàn)從“激波帶”到“激波線”的擬合識別。

2 算例驗證

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

2.1 定常斜激波

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

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

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

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

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

2.2 定常激波反射

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

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

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

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

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

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

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

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

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

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

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

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

圖19 激波干擾：幾何模型與壓強云圖

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

圖20 激波干擾：辨識出的激波帶

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

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

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

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

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

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

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

2.4 前向臺階非定常流動

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

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

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

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

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

圖26 前向臺階繞流：擬合的激波線與壓強云圖對比

圖27 前向臺階繞流：4個時刻擬合的激波線

3 結(jié) 論

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

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

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

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