橢圓內(nèi)聚激波面間斷化過(guò)程的激波動(dòng)力學(xué)分析

司東現(xiàn)，李祝飛，姬雋澤，張恩來(lái)，楊基明

中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系，合肥 230027

內(nèi)轉(zhuǎn)式進(jìn)氣道憑借壓縮效率高、氣動(dòng)阻力小等性能優(yōu)勢(shì)[1-3]，在新一代高超聲速飛行器研發(fā)中備受關(guān)注。然而，由于采用內(nèi)收縮幾何約束的設(shè)計(jì)特點(diǎn)[4-7]，內(nèi)轉(zhuǎn)式進(jìn)氣道內(nèi)部流動(dòng)往往面臨激波匯聚增強(qiáng)問(wèn)題[8]，因而激波的演化特征從根本上有別于二維平面激波[9]。

近年來(lái)，一系列以軸對(duì)稱內(nèi)錐流動(dòng)為代表的研究工作[8,10-14]相繼展開(kāi)，對(duì)軸對(duì)稱激波匯聚特征的認(rèn)知也不斷加深。然而，在實(shí)際工程應(yīng)用中，理想化的軸對(duì)稱構(gòu)型難以保證，偏離軸對(duì)稱的幾何設(shè)計(jì)方案[15-16]往往更為普遍。為揭示小幅偏離軸對(duì)稱流動(dòng)的典型特征，作者團(tuán)隊(duì)[17]前期提出了一種近軸對(duì)稱橢圓內(nèi)錐構(gòu)型。與軸對(duì)稱激波明顯不同的是，這種橢圓內(nèi)錐前緣激波的初始強(qiáng)度沿周向恒定，但曲率沿周向分布不均勻。在向下游匯聚增強(qiáng)的過(guò)程中，出現(xiàn)沿激波周向的強(qiáng)度非均勻性，并不斷加劇。當(dāng)非均勻性發(fā)展到一定程度，激波面甚至出現(xiàn)間斷，形成拐折結(jié)構(gòu)(Kink)。而一旦激波面出現(xiàn)間斷，即使激波在匯聚過(guò)程中持續(xù)增強(qiáng)，最終也有可能發(fā)生規(guī)則反射[17-18]，與軸對(duì)稱激波的“必然馬赫反射”[11-12,19]相比，出現(xiàn)顛覆性變化。由此可見(jiàn)，激波面從連續(xù)到間斷的轉(zhuǎn)變，對(duì)于橢圓內(nèi)錐流動(dòng)中的激波匯聚問(wèn)題，至關(guān)重要。

隨著三維曲面激波的應(yīng)用不斷深化，相關(guān)理論也不斷充實(shí)。其中具有代表性的是M?lder[20-22]提出的彎曲激波理論，該理論能夠在一定程度上揭示內(nèi)收縮約束下曲面激波的匯聚機(jī)理，尤其在軸對(duì)稱激波匯聚問(wèn)題中展示了可觀的潛力[14]。然而，對(duì)于激波面從連續(xù)到間斷轉(zhuǎn)變的問(wèn)題，現(xiàn)有的理論尚難以定量預(yù)測(cè)，并給出合理解釋。因此，嘗試探索新的研究方法和途徑十分必要。

鑒于三維問(wèn)題的復(fù)雜性，“空間降維”[23-24]不失為一種行之有效的分析方法。經(jīng)典的高超聲速等價(jià)原理[25-26]可以將繞細(xì)長(zhǎng)體的三維定常流動(dòng)等價(jià)轉(zhuǎn)化為二維非定常流動(dòng)。對(duì)于三維內(nèi)錐流動(dòng)中的激波匯聚，則可以等價(jià)為二維平面內(nèi)的激波內(nèi)收縮運(yùn)動(dòng)。此時(shí)，便可望從幾何激波動(dòng)力學(xué)(Geometrical Shock Dynamics，GSD)理論[23-24,27-29]中獲得解決等價(jià)問(wèn)題的思路。GSD理論對(duì)于解決非均勻激波匯聚、激波面從連續(xù)到間斷等問(wèn)題，具有獨(dú)到的優(yōu)勢(shì)。但從目前對(duì)激波動(dòng)力學(xué)的應(yīng)用來(lái)看，已有的研究大都集中于平面激波繞射、平面激波在變截面管道中傳播形狀和強(qiáng)度的演化、聚焦[30-33]或均勻柱激波運(yùn)動(dòng)[34-35]等易于解析的問(wèn)題，對(duì)于本文等價(jià)轉(zhuǎn)化而來(lái)的具有非均勻性的二維內(nèi)收縮運(yùn)動(dòng)激波，少有關(guān)注。此外，由于非均勻性的存在，傳統(tǒng)的GSD方法難以直接應(yīng)用。因此，GSD方法自身，也很有發(fā)展和拓展的必要。

針對(duì)前期[17]采用的橢圓內(nèi)錐構(gòu)型，利用高超聲速等價(jià)原理，將三維定常的橢圓內(nèi)錐流動(dòng)簡(jiǎn)化為二維平面內(nèi)的非定常流動(dòng)，力求對(duì)激波的演化過(guò)程得到另一角度的刻畫(huà)。進(jìn)一步地，發(fā)展出基于GSD理論的“波面-擾動(dòng)追蹤法”，并利用該方法考察了不同長(zhǎng)短軸比下橢圓激波在內(nèi)收縮運(yùn)動(dòng)中的演變過(guò)程，重點(diǎn)分析了非均勻性的發(fā)展和激波面從連續(xù)到間斷的演變特征。

1 模型和方法

1.1 三維橢圓內(nèi)錐模型

本文關(guān)注的三維定常橢圓內(nèi)錐流動(dòng)[17]，來(lái)流馬赫數(shù)Ma∞=6，來(lái)流靜壓p∞=891 Pa，靜溫T∞=101 K。采用如圖1(a)所示的橢圓內(nèi)錐模型，以模型入口橢圓中心為原點(diǎn)O，軸線為x方向，長(zhǎng)軸方向?yàn)閥方向，短軸方向?yàn)閦方向，來(lái)流方向與x方向保持一致。模型軸向長(zhǎng)度L=0.1 m，入口半長(zhǎng)軸長(zhǎng)度a=0.1 m，前緣壓縮角度沿周向保持恒定，為δ0=10°，φ為旋轉(zhuǎn)角。以模型入口橢圓的長(zhǎng)短軸比(Aspect Ratio,AR)為研究參數(shù)，AR=a/b,通過(guò)改變?nèi)肟诎攵梯Sb的長(zhǎng)度得到AR=1.11和AR=1.43這2種典型構(gòu)型，入口型線如圖1(b)所示。

圖1 橢圓內(nèi)錐模型示意圖和前緣型線

1.2 高超聲速等價(jià)原理及非定常數(shù)值模擬方法

針對(duì)小擾動(dòng)情況下的細(xì)長(zhǎng)體構(gòu)型，高超聲速等價(jià)原理[25-26]可以將三維定常流動(dòng)等價(jià)轉(zhuǎn)化為二維平面內(nèi)的非定常流動(dòng)。若自由來(lái)流沿x方向，則三維定常流中的流向位置坐標(biāo)x與二維非定常流中的時(shí)間t具有等價(jià)關(guān)系，即

x～V∞t

(1)

式中：V∞為三維定常流動(dòng)的來(lái)流速度。

對(duì)于本文研究的三維定常橢圓內(nèi)錐流動(dòng)，圖2給出了應(yīng)用高超聲速等價(jià)原理簡(jiǎn)化的示意圖。如圖2所示，在三維定常流場(chǎng)中，橢圓內(nèi)錐模型前緣(x0=0)產(chǎn)生橫向形狀與模型入口重合的三維內(nèi)聚激波。在y-z二維平面內(nèi)，則可以等價(jià)為初始時(shí)刻(t0=0)與橢圓內(nèi)錐模型入口形狀相同的二維橢圓活塞(W0)瞬間開(kāi)始做內(nèi)收縮運(yùn)動(dòng)，產(chǎn)生與W0形狀重合的內(nèi)收縮運(yùn)動(dòng)激波(S0)。在向下游發(fā)展的過(guò)程中，三維橢圓內(nèi)錐模型與二維內(nèi)收縮運(yùn)動(dòng)活塞始終可以建立等價(jià)關(guān)系，且在等價(jià)關(guān)系下，三維內(nèi)聚激波的橫向結(jié)構(gòu)與二維內(nèi)收縮運(yùn)動(dòng)激波結(jié)構(gòu)保持一致。在三維橢圓內(nèi)錐模型的尾緣位置(x=xt)，壁面(Wt)約束消失，激波(St)繼續(xù)向中心運(yùn)動(dòng)。

圖2 橢圓內(nèi)錐模型高超聲速等價(jià)原理示意圖

上述基于高超聲速等價(jià)原理建立的等價(jià)關(guān)系在小擾動(dòng)假設(shè)下(即激波角較小時(shí))成立，滿足βs～sinβs(βs為三維定常激波角)。然而，當(dāng)激波角較大時(shí)，小擾動(dòng)假設(shè)不能完全滿足。因此，在本文的具體案例中，根據(jù)實(shí)際的三維定常激波角對(duì)式(1)中的等價(jià)關(guān)系進(jìn)行了幾何修正。

為驗(yàn)證上述等價(jià)原理的準(zhǔn)確性，使用動(dòng)網(wǎng)格技術(shù)進(jìn)行了二維非定常數(shù)值模擬。1.3節(jié)中將二維數(shù)值結(jié)果與三維定常流動(dòng)結(jié)果[17]進(jìn)行了對(duì)比。二維計(jì)算域的邊界和內(nèi)部離散網(wǎng)格與三維橢圓內(nèi)錐模型入口邊界和離散網(wǎng)格保持一致，以三維定常流場(chǎng)中的來(lái)流靜壓和靜溫條件對(duì)二維計(jì)算域進(jìn)行初始化。

1.3 等價(jià)的二維非定常流場(chǎng)

本節(jié)以AR=1.43構(gòu)型為例，介紹利用高超聲速等價(jià)原理所得到的二維非定常流場(chǎng)基本結(jié)構(gòu)，并與對(duì)應(yīng)的三維定常流場(chǎng)結(jié)構(gòu)[17]進(jìn)行對(duì)比。

圖3(a)～圖3(d)分別給出了三維定常橢圓內(nèi)錐流場(chǎng)中x/L=0.10,1.40,1.60,2.00截面上無(wú)量綱密度(ρ/ρ∞)云圖，詳見(jiàn)文獻(xiàn)[17]，其中，ρ∞為三維定常流場(chǎng)中的來(lái)流密度。圖3(e)～圖3(f)分別展示了修正后對(duì)應(yīng)于圖3(a)～圖3(d)位置的等價(jià)二維非定常流場(chǎng)的無(wú)量綱密度云圖，對(duì)應(yīng)的時(shí)刻分別為t/ΔT=0.10,1.47,1.69,2.12，其中，ΔT滿足式(1)等價(jià)關(guān)系(L～V∞ΔT)。從圖3中的三維定常和二維非定常流場(chǎng)結(jié)構(gòu)可以看出，初始橢圓形激波在匯聚過(guò)程中，非均勻特征逐漸突出。以激波面出現(xiàn)間斷為標(biāo)志，可以將非定常內(nèi)收縮運(yùn)動(dòng)激波(S)的演化過(guò)程劃分為連續(xù)激波非均勻強(qiáng)化和間斷激波平面化發(fā)展兩個(gè)階段。在連續(xù)激波非均勻強(qiáng)化階段，圖3(a)和圖3(b)所示的x/L=0.10,1.40截面上的三維定常流場(chǎng)結(jié)構(gòu)與圖3(e)和圖3(f)所示的t/ΔT=0.10,1.47時(shí)刻的二維非定常流場(chǎng)結(jié)構(gòu)等價(jià)。在內(nèi)收縮運(yùn)動(dòng)過(guò)程中，雖然S的初始強(qiáng)度沿周向相同，但是初始橢圓形激波在長(zhǎng)軸端點(diǎn)處的曲率最大，激波在長(zhǎng)軸端點(diǎn)附近匯聚更快，激波強(qiáng)度和波后密度也會(huì)更快地增大，使得S沿周向逐漸演化出強(qiáng)度的非均勻性，并且隨著S匯聚增強(qiáng)，強(qiáng)度非均勻性也逐漸加劇。尤其在沿長(zhǎng)軸和短軸2個(gè)方向上，S的強(qiáng)度差異最為突出。從t/ΔT=1.47時(shí)刻(圖3(f))可以看出，S在長(zhǎng)軸方向上的強(qiáng)度已經(jīng)明顯大于在短軸方向上的強(qiáng)度，波后呈現(xiàn)為“紅黃色”高密度區(qū)域。總體來(lái)說(shuō)，雖然在這一階段，S的非均勻性逐漸顯著，但激波面始終維持連續(xù)、光滑的形狀。然而，隨著S進(jìn)一步匯聚，激波面難以繼續(xù)維持光滑、連續(xù)的形狀。在三維定常流場(chǎng)中x/L=1.60截面(圖3(c))，等價(jià)的二維非定常流場(chǎng)結(jié)構(gòu)在t/ΔT=1.69時(shí)刻(圖3(g))，強(qiáng)烈的周向不均勻性使S的激波面出現(xiàn)間斷，形成Kink，此后進(jìn)入間斷激波的平面化發(fā)展階段。在t/ΔT=2.12時(shí)刻(圖3(h))，Kink更加明顯，S被中心對(duì)稱的4個(gè)Kink劃分成兩對(duì)強(qiáng)度不同的平面化激波段：沿長(zhǎng)軸運(yùn)動(dòng)的激波段(S1)和沿短軸運(yùn)動(dòng)的激波段(S2)。綜合上述兩個(gè)階段，通過(guò)對(duì)比圖3(a)～圖3(d)和圖3(e)～圖3(h)可以看出，二維運(yùn)動(dòng)激波結(jié)構(gòu)與三維定常激波的橫向結(jié)構(gòu)吻合較好，驗(yàn)證了高超聲速等價(jià)原理在三維橢圓內(nèi)錐流場(chǎng)中的適用性。

圖3 AR=1.43構(gòu)型三維定常及等價(jià)的二維非定常流場(chǎng)結(jié)構(gòu)

GSD對(duì)于解決二維內(nèi)收縮激波運(yùn)動(dòng)問(wèn)題，有著快捷、高效的優(yōu)勢(shì)。利用GSD方法，結(jié)合等價(jià)原理，可望為定量預(yù)測(cè)三維定常橢圓內(nèi)錐激波的非均勻演化過(guò)程，并揭示激波面發(fā)展出間斷的內(nèi)在機(jī)理，提供新的途徑。

2 激波動(dòng)力學(xué)波面-擾動(dòng)追蹤法

由1.3節(jié)分析可知，橢圓內(nèi)收縮運(yùn)動(dòng)過(guò)程中的非均勻演化機(jī)理是值得探究的重點(diǎn)之一。然而，已有的GSD方法還難以解決激波的非均勻擾動(dòng)傳播問(wèn)題。因此，本文從GSD基本方程出發(fā)，發(fā)展出對(duì)激波非均勻擾動(dòng)傳播問(wèn)題具有獨(dú)到適用性的“波面-擾動(dòng)追蹤法”。下面，就其相關(guān)的原理和數(shù)值算法進(jìn)行介紹。

2.1 原 理

最早由Chester[36]、Chisnell[37]以及Whitham[38]建立的CCW關(guān)系是GSD理論的基礎(chǔ)，它描述了激波在一維變截面管道中運(yùn)動(dòng)時(shí)，激波馬赫數(shù)(Mas)隨管道面積(A)變化的關(guān)系：

(2)

式中：K(Mas)為Mas的緩變函數(shù)。進(jìn)一步地，Whitham[39]利用圖4所示的正交曲線坐標(biāo)系(α為常數(shù)，表示不同時(shí)刻的激波面；β為常數(shù)，表示激波面各部分法矢量的積分曲線，即射線)，并結(jié)合微分關(guān)系：

(3)

將GSD推廣到了二維，建立了擾動(dòng)沿激波面?zhèn)鞑サ母拍?，并給出了擾動(dòng)沿激波面?zhèn)鞑サ膬纱靥卣骶€方程。在y-z平面內(nèi)，兩簇特征線方程可表示為

(4)

(5)

式中：A為射線管(相鄰兩條射線之間看作一維變截面管道)面積；θ為射線與水平方向(y軸)的夾角；c為擾動(dòng)傳播速度的系數(shù)；ν為擾動(dòng)軌跡(即特征線)與射線的夾角。若式(4)(式(5))對(duì)左行(右行)特征線在整個(gè)流場(chǎng)成立，則為簡(jiǎn)單波擾動(dòng)，否則為雙向擾動(dòng)[27,34-35]。

圖4 正交曲線坐標(biāo)系(α, β)中的激波位置和射線[27]

在等價(jià)的橢圓激波內(nèi)收縮運(yùn)動(dòng)過(guò)程中，彎曲橢圓激波面上分布著雙向傳播的擾動(dòng)，并且雙向擾動(dòng)始終貫穿激波的非均勻演化過(guò)程，結(jié)合式(4)和式(5)即可計(jì)算出沿特征線的參數(shù)變化。然而，由于特征線式(4)和式(5)不含時(shí)間變量，計(jì)算出的參數(shù)是在空間上分布的，難以將這些參數(shù)對(duì)應(yīng)在任意確定時(shí)刻的激波面上。對(duì)于激波非均勻匯聚問(wèn)題而言，得到不同時(shí)刻的激波面及參數(shù)變化是必要的。為此，本文提出了一種“波面-擾動(dòng)追蹤法”。

圖5以一般的二維曲面運(yùn)動(dòng)激波為例，展示了該方法的基本原理。遵循數(shù)值求解思想，將初始(t0時(shí)刻)激波面用一系列點(diǎn)(圖5中t0時(shí)刻激波面上的綠色點(diǎn))離散。在激波運(yùn)動(dòng)過(guò)程中的任意t時(shí)刻和t+Δt時(shí)刻(Δt為時(shí)間增量)，建立波面-擾動(dòng)追蹤關(guān)系。圖5給出了t時(shí)刻激波面上任意3個(gè)相鄰離散點(diǎn)的位置矢量xt,i-1，xt,i和xt,i+1，它們?cè)趶某跏技げ鎮(zhèn)鞑ザ鴣?lái)的特征線(圖5中藍(lán)色線)上。為區(qū)分?jǐn)_動(dòng)傳播的2個(gè)方向，站在激波面上面向激波的運(yùn)動(dòng)方向看，定義向左傳播的特征線為左行(C+)特征線，向右傳播的特征線為右行(C-)特征線。

圖5 “波面-擾動(dòng)追蹤法”原理示意圖

根據(jù)GSD基本原理，激波從t時(shí)刻運(yùn)動(dòng)到t+Δt時(shí)刻，激波面上的各離散點(diǎn)沿射線運(yùn)動(dòng)。由于射線方向在任意時(shí)刻均與激波面垂直，可以用t時(shí)刻激波面上離散點(diǎn)的法向(圖5中紅色箭頭)代替射線方向，射線的面積為圖5中相鄰黑色虛線(其中，黑色虛線為圖5中相鄰離散點(diǎn)的中垂線)包圍的激波面弧長(zhǎng)。實(shí)際上，圖5中相鄰黑色虛線為一組射線，同樣垂直于t時(shí)刻的激波面，因此可以將激波面看作被射線分割的若干波面微元，每一個(gè)波面微元沿著射線管道的傳播可以看作是激波在準(zhǔn)一維管道[39]中運(yùn)動(dòng)。從t時(shí)刻到t+Δt時(shí)刻，射線管面積改變，波面微元的激波馬赫數(shù)也隨之改變，兩者變化滿足CCW關(guān)系。與此同時(shí)，擾動(dòng)也沿激波面?zhèn)鞑?，因此，可以利用特征線關(guān)系，計(jì)算特征線從t時(shí)刻到t+Δt時(shí)刻的軌跡。再結(jié)合t+Δt時(shí)刻的激波面，可以確定C+和C-兩組特征線傳播到t+Δt時(shí)刻激波面上的位置和馬赫數(shù)。綜上所述，從t0時(shí)刻開(kāi)始，任意時(shí)刻的激波面位置，以及沿激波面?zhèn)鞑サ臄_動(dòng)均可以得到。

2.2 算 法

在均勻靜止介質(zhì)中運(yùn)動(dòng)時(shí)[27]，正交曲線坐標(biāo)系下的激波面由α=a*t描述，其中:a*為波前介質(zhì)聲速。若在y-z平面內(nèi)令激波面位置x=(y,z)，激波面單位法向量n=(cosθ, sinθ)，則可根據(jù)式(3)寫(xiě)出正交坐標(biāo)系下激波面隨時(shí)間推進(jìn)的矢量微分形式，即

(6)

按照2.1節(jié)中介紹的用離散點(diǎn)代替連續(xù)激波面的思想，則式(6)可轉(zhuǎn)化為常微分方程：

(7)

式中：下標(biāo)i表示離散點(diǎn)編號(hào)。對(duì)式(7)應(yīng)用三階Runge-Kutta格式[40]進(jìn)行數(shù)值積分，實(shí)現(xiàn)激波面隨時(shí)間推進(jìn)。在數(shù)值積分中，離散點(diǎn)激波馬赫數(shù)(Masi)變化和射線管面積(Ai)變化滿足

(8)

式中：Ai用離散點(diǎn)鄰近的激波面平均弧長(zhǎng)代替，即

(9)

其中：si為以激波面端點(diǎn)(i=1對(duì)應(yīng)的離散點(diǎn)位置)為起點(diǎn)，沿激波面建立的弧坐標(biāo)下的離散點(diǎn)位置，即

(10)

為確定任意t時(shí)刻激波面離散點(diǎn)法向量，構(gòu)建弧坐標(biāo)下兩組數(shù)據(jù)點(diǎn)(si(t),yi(t))i=1,2,…,N和(si(t),zi(t))i=1,2,…,N，用修正的三次Akima[41]插值得到兩個(gè)函數(shù)Y(s(t))和Z(s(t))。激波面上離散點(diǎn)的單位法向量ni(t)滿足

i=1,2,…,N

(11)

式中：上標(biāo)“′”表示對(duì)弧坐標(biāo)s的微分。對(duì)于激波面端點(diǎn)在實(shí)體固壁上的情況，由于壁面本身可以視為一條射線，端點(diǎn)法向量沿壁面方向。

在數(shù)值計(jì)算中，Δt需要滿足穩(wěn)定性條件(即Courant-Friedrichs-Lew(CFL)條件[33])。同時(shí)，為避免由數(shù)值計(jì)算誤差導(dǎo)致的射線相交(射線管面積收縮為0)，引入射線不相交條件[33]，調(diào)整Δt。

如果已知任意t時(shí)刻的激波面形狀和強(qiáng)度分布(如圖5中t時(shí)刻激波面上的綠色點(diǎn)坐標(biāo)及參數(shù))，按照上述方法，可計(jì)算出t+Δt時(shí)刻的激波面形狀和強(qiáng)度分布(如圖5中t+Δt時(shí)刻激波面上的兩個(gè)邊界點(diǎn)、內(nèi)部黑色點(diǎn)坐標(biāo)及參數(shù))。至此，完成了在一個(gè)時(shí)間步內(nèi)，對(duì)激波面的追蹤。

在激波運(yùn)動(dòng)過(guò)程中，還需要計(jì)算擾動(dòng)沿激波面的傳播，建立波面-擾動(dòng)追蹤關(guān)系。以初始激波面上第m條C-特征線為例，利用式(5)計(jì)算擾動(dòng)由t時(shí)刻傳播至t+Δt時(shí)刻時(shí)，沿特征線的參數(shù)變化。

如圖5所示，若第m條特征線在t時(shí)刻傳播到激波面上x(chóng)t,i點(diǎn)位置，根據(jù)式(5)該特征線在繼續(xù)傳播過(guò)程中的軌跡滿足

z=tan(θt,i-ν(Mast,i))y+zt,i-

tan(θt,i-ν(Mast,i))yt,i

(12)

式中：(yt,i,zt,i)、Mast,i和θt,i分別為t時(shí)刻激波面上x(chóng)t,i點(diǎn)的坐標(biāo)、運(yùn)動(dòng)激波馬赫數(shù)和射線角。其中，θt,i由t時(shí)刻激波面上x(chóng)t,i點(diǎn)的單位法向量nt,i確定，即

(13)

其中：nt,i,y和nt,i,z分別為nt,i在y軸和z軸方向上的分量。

(14)

圖6 擾動(dòng)傳播至邊界求解原理示意圖

(Mast+Δt,i-1-Mast+Δt,i)

(15)

綜合式(12)～式(15)，可以得到初始激波面產(chǎn)生的所有擾動(dòng)從t時(shí)刻到t+Δt時(shí)刻的傳播軌跡和馬赫數(shù)變化(如圖5中t+Δt時(shí)刻所有綠色點(diǎn)的坐標(biāo)和參數(shù))。

除初始激波面外，以任意時(shí)刻激波面上的特征線節(jié)點(diǎn)(如圖5中特征線與激波面的綠色交點(diǎn))和固壁邊界上的端點(diǎn)作為激波面的離散點(diǎn)，可由式(11)求解離散點(diǎn)的法向量；再根據(jù)該時(shí)刻離散點(diǎn)的坐標(biāo)、法向量、馬赫數(shù)分布，以及特征線傳播過(guò)程，求解下一時(shí)刻的激波形狀和參數(shù)。如此進(jìn)行循環(huán)推進(jìn)，最終，不僅可以得到二維非定常激波面隨時(shí)間行進(jìn)過(guò)程中的幾何形狀及參數(shù)變化，同時(shí)還能追蹤擾動(dòng)隨著激波推進(jìn)的傳播過(guò)程。至此，建立了“波面-擾動(dòng)追蹤法”。

另外，當(dāng)同簇特征線趨于相交時(shí)(射線管面積趨于0)，需終止計(jì)算。在數(shù)值計(jì)算中，當(dāng)t時(shí)刻激波面上發(fā)出的兩個(gè)相鄰?fù)靥卣骶€的交點(diǎn)到對(duì)應(yīng)激波面的距離lt，比t時(shí)刻的平均弧長(zhǎng)Δst低3個(gè)數(shù)量級(jí)時(shí)，認(rèn)為同簇特征線趨于相交，終止迭代。

為了便于展示和理解上述“波面-擾動(dòng)追蹤法”的實(shí)施過(guò)程，圖7給出了該算法的流程圖。

圖7 “波面-擾動(dòng)追蹤法”算法流程圖

2.3 算法驗(yàn)證

根據(jù)1.2節(jié)中的等價(jià)關(guān)系，將三維橢圓內(nèi)錐模型前緣(x0=0)處的激波角記為λ0(來(lái)流馬赫數(shù)Ma∞=6，前緣壓縮角度δ0=10°的楔產(chǎn)生的斜激波的激波角)，則初始時(shí)刻(t0=0)二維內(nèi)收縮運(yùn)動(dòng)激波馬赫數(shù)為Mas0=Ma∞tanλ0=1.90。以AR=1.43構(gòu)型為例，取周向[-π/2, π/2]區(qū)間內(nèi)的1/2初始橢圓運(yùn)動(dòng)激波，按等弧長(zhǎng)均勻離散得到N=401個(gè)初始離散點(diǎn)，運(yùn)用發(fā)展的“波面-擾動(dòng)追蹤法”計(jì)算不同時(shí)刻的激波面位置和強(qiáng)度分布。通過(guò)與三維定常數(shù)值模擬[17]結(jié)果對(duì)比，驗(yàn)證“波面-擾動(dòng)追蹤法”。

圖8(a)中虛線展示了從三維定常流場(chǎng)中提取的x/L=0.65,1.30典型截面上的激波面，實(shí)線展示了采用“波面-擾動(dòng)追蹤法”得到的等價(jià)二維非定常流場(chǎng)中對(duì)應(yīng)的激波面(t/ΔT=0.65,1.30)。圖8(b)對(duì)比了沿激波面的激波馬赫數(shù)Mas分布，其中旋轉(zhuǎn)角φ與三維模型保持一致(見(jiàn)圖1(a))。圖8(b)中的實(shí)線為采用“波面-擾動(dòng)追蹤法”得到的Mas分布，圖8(b)中的虛線表示由三維定常流場(chǎng)得到的等價(jià)Mas分布。根據(jù)三維定常流場(chǎng)x/L截面上激波前后的壓比，按斜激波關(guān)系[42]計(jì)算得到激波角λ后，再借助等價(jià)關(guān)系Mas=Ma∞tanλ，換算出等價(jià)的Mas分布。通過(guò)對(duì)比圖8中的實(shí)線和虛線可以看出，盡管GSD方法因僅考慮了激波面，存在一定的誤差[33,43]，但本文基于GSD發(fā)展的“波面-擾動(dòng)追蹤法”的計(jì)算結(jié)果與三維定常結(jié)果吻合良好。

以此為基礎(chǔ)，可以進(jìn)一步地利用“波面-擾動(dòng)追蹤法”定量地分析橢圓內(nèi)收縮運(yùn)動(dòng)激波的非均勻匯聚過(guò)程，以揭示三維定常橢圓內(nèi)聚流場(chǎng)中Kink的形成機(jī)理。

圖8 AR=1.43構(gòu)型不同位置(時(shí)刻)激波面和周向馬赫數(shù)分布對(duì)比

3 二維流場(chǎng)的幾何激波動(dòng)力學(xué)分析

本節(jié)利用“波面-擾動(dòng)追蹤法”計(jì)算了長(zhǎng)短軸比較小的AR=1.11和較大的AR=1.43這2種典型情況，分析激波的形狀及強(qiáng)度的演變過(guò)程和機(jī)理。

圖9(a)給出了AR=1.11構(gòu)型在t/ΔT=0,0.91,1.74,2.44典型時(shí)刻的激波面位置和兩簇特征線，圖9(b)給出了對(duì)應(yīng)于圖9(a)中各時(shí)刻的周向激波馬赫數(shù)(Mas)分布，用以展示激波強(qiáng)度的周向不均勻性。

從圖9(a)和圖9(b)中t/ΔT=0時(shí)刻的激波面和周向Mas分布可以看出，初始時(shí)刻連續(xù)光滑的橢圓激波，雖然沿周向強(qiáng)度均勻分布，但激波自身曲率的不均勻(即幾何上偏離軸對(duì)稱的激波形狀)將在激波傳播過(guò)程中顯現(xiàn)和產(chǎn)生作用，進(jìn)而影響后續(xù)的非均勻演變過(guò)程。由1.3節(jié)可知，初始激波曲率在長(zhǎng)軸端點(diǎn)處最大，激波在長(zhǎng)軸端點(diǎn)附近將更快地匯聚并增強(qiáng)。因此，在后續(xù)的t/ΔT=0.91時(shí)刻(見(jiàn)圖9(b))，激波強(qiáng)度從長(zhǎng)軸到短軸方向呈現(xiàn)逐漸降低的分布趨勢(shì)。隨著激波進(jìn)一步匯聚，強(qiáng)度和幾何非均勻性都不斷強(qiáng)化。對(duì)比圖9(b)中t/ΔT=0.91,1.74,2.44時(shí)刻沿激波面周向的Mas分布可以看出，長(zhǎng)軸和短軸附近的激波強(qiáng)度差異越來(lái)越顯著。在t/ΔT=2.44時(shí)刻，強(qiáng)烈的非均勻性使得激波面難以維持連續(xù)、光滑的形態(tài)，激波面上出現(xiàn)Kink，進(jìn)而從連續(xù)轉(zhuǎn)變至間斷。

圖9 AR=1.11構(gòu)型激波面、特征線以及相應(yīng)的激波馬赫數(shù)分布

增大長(zhǎng)短軸比至AR=1.43時(shí)，結(jié)果如圖10所示。從圖10(a)給出的不同時(shí)刻激波面形狀和圖10(b)給出的激波面周向Mas分布可知，相比于長(zhǎng)短軸比較小的AR=1.11情況而言，激波在匯聚過(guò)程中長(zhǎng)、短軸2個(gè)方向激波強(qiáng)度的差異凸顯得更快，而在激波的匯聚過(guò)程中，Kink也更早地出現(xiàn)。

圖10 AR=1.43構(gòu)型激波面、特征線以及相應(yīng)的激波馬赫數(shù)分布

結(jié)合前期對(duì)三維激波的研究[17]可知，若初始激波曲率均勻，即圓柱激波(AR=1.0)，激波在匯聚過(guò)程中將始終維持連續(xù)、光滑的波面形狀，直至馬赫盤(pán)形成[17]。一旦初始激波偏離圓柱激波形狀(AR增大)，非均勻性在匯聚過(guò)程中的強(qiáng)化不可避免，隨之而來(lái)的便是激波從連續(xù)到間斷化的發(fā)展。激波偏離軸對(duì)稱后的間斷化發(fā)展，改變了原本軸對(duì)稱激波中的“無(wú)限匯聚”模式，最終可能顛覆軸對(duì)稱激波“必然馬赫反射”的規(guī)律[17]?？梢?jiàn)，初始連續(xù)光滑的激波面在發(fā)展過(guò)程中是如何形成Kink的，能否在理論上進(jìn)行描述和預(yù)測(cè)Kink在激波面上的形成位置是探究激波非均勻匯聚內(nèi)在機(jī)理的關(guān)鍵問(wèn)題。下面，借助本文發(fā)展的“波面-擾動(dòng)追蹤法”計(jì)算出的特征線，可以做更進(jìn)一步的GSD分析。

從圖9(a)中特征線分布可以看出，初始橢圓激波面上分布著雙向傳播的擾動(dòng)，激波在內(nèi)收縮運(yùn)動(dòng)過(guò)程中，左行擾動(dòng)沿著激波面以Wd=a*·[(Mas2-1)K(Mas)]1/2(與聲速a*和激波馬赫數(shù)Mas正相關(guān)[27])的速度向φ增大的方向傳播，右行擾動(dòng)沿著激波面以Wd的速度向φ減小的方向傳播，使激波自身不斷增強(qiáng)。雖然初始激波強(qiáng)度均勻分布，但由于幾何非均勻性的影響，以及后續(xù)匯聚過(guò)程中強(qiáng)度和曲率兩方面的強(qiáng)化耦合作用，雙向擾動(dòng)非均勻傳播呈現(xiàn)越來(lái)越明顯的非均勻“Shock-Compression”擾動(dòng)特征[27]。

從圖9(b)中激波面周向Mas分布可以看出，在第四象限內(nèi)(φ∈[-π/2, 0])，激波強(qiáng)度從φ=-π/2～0位置逐漸增大。因此，對(duì)于第四象限內(nèi)某時(shí)刻處于任意φ1和φ2位置的兩同向右行擾動(dòng)(以C-特征線為例)，若φ1<φ2，則φ1位置的右行擾動(dòng)沿著激波面的傳播速度W1小于φ2位置的右行擾動(dòng)沿著激波面的傳播速度W2。雖然激波面在匯聚過(guò)程中持續(xù)收縮，但是由于W1

由上述分析可知，一旦第四象限內(nèi)的同向右行擾動(dòng)跨過(guò)φ=0位置進(jìn)入第一象限，沿著激波面的傳播速度相對(duì)大小關(guān)系便發(fā)生轉(zhuǎn)變。因此，從圖9(a)中可以看出，從第四象限內(nèi)傳播至第一象限內(nèi)的右行擾動(dòng)不斷靠近初始從φ=0位置產(chǎn)生的右行擾動(dòng)，即第四象限內(nèi)C-特征線在圖9(a)中紅色區(qū)域內(nèi)追趕初始從φ=0位置發(fā)出的C-特征線。由于激波的對(duì)稱性，左行擾動(dòng)的傳播過(guò)程與右行擾動(dòng)的傳播過(guò)程呈現(xiàn)對(duì)稱的趨勢(shì)，第一象限內(nèi)C+特征線在圖9(a)中藍(lán)色區(qū)域內(nèi)追趕初始從φ=0位置發(fā)出的C+特征線。在t/ΔT=2.44時(shí)刻附近，同簇C-特征線在第一象限中(0.07, 0.03)位置，同簇C+特征線在第四象限中(0.07, -0.03)位置分別相交(見(jiàn)圖9(a)綠色點(diǎn))。在特征線交點(diǎn)處，激波面出現(xiàn)間斷，并形成Kink。同時(shí)，也意味著激波參數(shù)的間斷。從圖9(b)中t/ΔT=2.44時(shí)刻的激波面周向Mas分布可以看出，激波強(qiáng)度在出現(xiàn)Kink位置附近急劇變化。而激波在長(zhǎng)軸附近強(qiáng)、短軸附近弱的顯著差異，也隨著Kink的出現(xiàn)，被分割成強(qiáng)、弱兩對(duì)激波段。

隨著長(zhǎng)短軸比增大至AR=1.43，特征線的發(fā)展趨勢(shì)與AR=1.11構(gòu)型類似，但由于初始激波的幾何非均勻性增強(qiáng)，擾動(dòng)會(huì)更早地向長(zhǎng)軸附近聚集。如圖10(a)所示，在t/ΔT=1.61時(shí)刻，同簇C-特征線在第一象限中(0.39, 0.05)位置，同簇C+特征線在第四象限中(0.39, -0.05)位置相交(見(jiàn)圖10(a)綠色點(diǎn))。與AR=1.11時(shí)相比，激波面上更早地出現(xiàn)Kink。

由上述分析可知，借助“波面-擾動(dòng)追蹤法”，可以快速地計(jì)算出二維橢圓內(nèi)收縮運(yùn)動(dòng)激波形成Kink的時(shí)間，以及Kink在y-z平面內(nèi)的位置。繼而根據(jù)等價(jià)關(guān)系式(1)，可以確定三維定常橢圓內(nèi)錐流場(chǎng)中激波面出現(xiàn)Kink的位置。因此，“波面-擾動(dòng)追蹤法”能夠揭示三維定常橢圓內(nèi)錐流場(chǎng)中Kink的形成機(jī)理。

結(jié)合1.3節(jié)中三維激波的橫向結(jié)構(gòu)可知，形成Kink之后，激波面被中心對(duì)稱的Kink劃分為兩組趨于平面化發(fā)展的激波段。不過(guò)，一旦激波面上出現(xiàn)Kink，本文發(fā)展的“波面-擾動(dòng)追蹤法”就面臨新的挑戰(zhàn)。激波面上出現(xiàn)Kink后，需要引入“Shock-Shock”關(guān)系[27]來(lái)進(jìn)行有針對(duì)性的描述。以往的“Shock-Shock”作用問(wèn)題，如平面運(yùn)動(dòng)激波在斜楔面上的馬赫反射問(wèn)題[27]，Kink(即馬赫反射的三波點(diǎn))兩側(cè)的激波面(入射激波和馬赫桿)均被認(rèn)為是強(qiáng)度均勻的激波面。然而，對(duì)于本文所面臨的問(wèn)題來(lái)說(shuō)，Kink兩側(cè)的激波面顯然都是非均勻的。換言之，此處所面臨的是更加復(fù)雜的雙向“Shock-Compression”擾動(dòng)與“Shock-Shock”擾動(dòng)相互作用問(wèn)題。如何應(yīng)對(duì)這一新的難題，目前尚處于探索之中。

4 結(jié) 論

利用高超聲速等價(jià)原理，將三維定常橢圓內(nèi)錐流動(dòng)簡(jiǎn)化為二維平面內(nèi)的非定常流動(dòng)。針對(duì)等價(jià)后的二維非定常流動(dòng)，基于GSD原理提出了“波面-擾動(dòng)追蹤法”，定量預(yù)測(cè)和分析了不同長(zhǎng)短軸比下初始強(qiáng)度相同的橢圓激波沿周向非均勻匯聚演變過(guò)程，主要得到以下結(jié)論：

1) 基于激波動(dòng)力學(xué)原理提出了“波面-擾動(dòng)追蹤法”，該方法既能夠得到二維非定常激波在匯聚過(guò)程中的波面演變特征及參數(shù)變化，又能夠追蹤擾動(dòng)沿著激波面的傳播過(guò)程。

2) 在初始橢圓激波自身產(chǎn)生的雙向非均勻“Shock-Compression”擾動(dòng)作用下，非均勻性不斷強(qiáng)化，隨著同向“Shock-Compression”擾動(dòng)的聚集，激波面從連續(xù)發(fā)展出間斷，形成拐折結(jié)構(gòu)(Kink)。

3) 隨著長(zhǎng)短軸比的增加，初始激波面的幾何非均勻性增強(qiáng)，同簇?cái)_動(dòng)會(huì)更快地聚集，激波面上更早形成Kink。