姬雋澤，李祝飛，*，司東現(xiàn)，張 濤，施崇廣，尤延鋮，楊基明

(1. 中國科學(xué)技術(shù)大學(xué)，合肥 230027；2. 廈門大學(xué)，廈門 361102)

0 引 言

近年來，內(nèi)轉(zhuǎn)式進(jìn)氣道[1-4]憑借諸多性能優(yōu)勢，在吸氣式高超聲速飛行器的研發(fā)中展示了誘人的潛力。然而，由于內(nèi)轉(zhuǎn)式進(jìn)氣道普遍是基于軸對稱內(nèi)收縮基準(zhǔn)流場進(jìn)行設(shè)計(jì)的，其內(nèi)部流動(dòng)存在匯聚效應(yīng)[5-8]，與二元進(jìn)氣道相比更為復(fù)雜，甚至存在本質(zhì)性差異。一般而言，二元進(jìn)氣道內(nèi)部的激波相互作用通常表現(xiàn)為平面激波的入射和反射[9-10]，并存在規(guī)則反射和馬赫反射兩種反射類型[11]。然而，在軸對稱內(nèi)收縮流場中[12-14]，內(nèi)錐形的入射激波會(huì)持續(xù)向中心匯聚增強(qiáng)，最終必然在軸線上發(fā)生馬赫反射[15-16]。在工程實(shí)際中，馬赫反射區(qū)域具有強(qiáng)壓縮損失和弱抗反壓能力等特性[12,17]，會(huì)給進(jìn)氣道性能帶來嚴(yán)重的負(fù)面影響，因此應(yīng)極力避免。

長久以來，通過在軸對稱內(nèi)收縮基準(zhǔn)流場中設(shè)置中心體將軸線區(qū)域掩蓋，進(jìn)而規(guī)避馬赫反射的措施[12,18]，在內(nèi)轉(zhuǎn)式進(jìn)氣道設(shè)計(jì)中被廣泛采納。實(shí)際上，這種措施在給進(jìn)氣道設(shè)計(jì)帶來便利的同時(shí)，并未促進(jìn)流動(dòng)機(jī)理認(rèn)識的提升，甚至制約了設(shè)計(jì)理論及方法的持續(xù)創(chuàng)新。因此，隨著對內(nèi)轉(zhuǎn)式進(jìn)氣道流場品質(zhì)的要求日益嚴(yán)苛，越來越多的研究者開始重視軸對稱內(nèi)錐流動(dòng)中的激波反射問題，并嘗試探索理論求解方法。譚廉華等[19-20]從軸對稱內(nèi)錐形激波和馬赫盤入手，較早地給出了確定它們形狀的理論關(guān)系式。近年來，彎曲激波理論[21-25]（Curved shock theory，簡稱CST）逐步得到廣泛應(yīng)用。該理論可以在已知內(nèi)錐流動(dòng)中入射激波和馬赫盤形狀的前提下，計(jì)算出反射激波的局部形狀以及波后流動(dòng)的高階參數(shù)。將上述兩種理論相結(jié)合，似乎可以逐次確定軸對稱內(nèi)錐流動(dòng)中原本難以解析的激波形狀，給理論求解軸對稱內(nèi)錐形激波反射問題帶來了新的曙光。然而，應(yīng)用這兩種理論均涉及到馬赫盤位置這一關(guān)鍵參數(shù)。由于目前缺乏對軸對稱內(nèi)錐流場中馬赫盤形成機(jī)理及馬赫盤下游流動(dòng)的深入認(rèn)識，從理論上確定馬赫盤位置的研究仍鮮見報(bào)道。

根據(jù)平面激波反射理論可知[11]，當(dāng)入射激波的角度β大于von Neumann準(zhǔn)則規(guī)定的激波角βN時(shí)，有可能發(fā)生馬赫反射。對于平面激波的馬赫反射結(jié)構(gòu)，入射激波、反射激波和馬赫桿的波后氣流偏轉(zhuǎn)角關(guān)系（即三激波理論）是唯一確定的。由于平面激波的β保持恒定，三激波理論可以在平面激波的任意位置成立。于是，馬赫桿的形成位置僅由其下游流動(dòng)決定[11,26]。

與平面激波相比，軸對稱內(nèi)錐形激波隨著向中心匯聚，自身的激波角β具有持續(xù)增大的特點(diǎn)，使得影響馬赫盤位置的因素變得更為復(fù)雜。在軸對稱內(nèi)錐形激波的入射過程中，當(dāng)β尚未增大至βN時(shí)，不能形成馬赫反射結(jié)構(gòu)。換言之，軸對稱內(nèi)錐形入射激波自身的匯聚增強(qiáng)過程會(huì)限制馬赫盤可能存在的空間范圍。進(jìn)一步地，雖然三激波理論在β超出βN之后成立，但是三激波理論中的反射激波后的氣流偏轉(zhuǎn)角不是唯一確定的，而是與當(dāng)?shù)厝肷浼げǖ摩戮o密相關(guān)。由此可見，在確定軸對稱內(nèi)錐形激波反射流場中馬赫盤位置時(shí)，入射激波自身的強(qiáng)度變化需要予以考慮。然而，關(guān)于馬赫盤位置與激波匯聚增強(qiáng)過程之間的關(guān)系，目前鮮見報(bào)道。

此外，與平面激波馬赫反射類似的是，在軸對稱內(nèi)錐形激波反射流場中馬赫盤下游流動(dòng)也會(huì)影響馬赫盤位置。Gounko[27]通過改變軸對稱內(nèi)錐模型的壁面長度，獲得了不同的馬赫盤下游流動(dòng)條件及馬赫盤位置，但并沒有闡明馬赫盤下游流動(dòng)對其形成位置的影響機(jī)理。在下游流動(dòng)影響下，馬赫盤的形成及演化規(guī)律仍有待系統(tǒng)地研究。

綜上所述，本文綜合考慮軸對稱內(nèi)錐形激波匯聚過程及馬赫盤下游流動(dòng)兩方面的因素，采用無黏數(shù)值模擬，考察內(nèi)收縮直錐的壁面前緣角度和壁面長度對馬赫盤位置的影響，揭示形成馬赫盤的內(nèi)在機(jī)制及其演化規(guī)律，以期豐富對軸對稱基準(zhǔn)流場的認(rèn)識并支撐內(nèi)轉(zhuǎn)式進(jìn)氣道設(shè)計(jì)潛力的挖掘。

1 模型和數(shù)值方法

如圖1所示，利用軸對稱內(nèi)收縮直錐生成內(nèi)錐形激波，采用原點(diǎn)位于模型入口圓心處的軸對稱坐標(biāo)系，其中x為軸線方向，r為半徑方向。該模型的壁面前緣半徑R= 100 mm，前緣壓縮角度為θw，無量綱的壁面長度為w/R。參考內(nèi)轉(zhuǎn)式進(jìn)氣道基準(zhǔn)流場設(shè)計(jì)參數(shù)[28]，選取的壁面前緣角度范圍為θw= 8°～15°，通過改變?chǔ)葁得到具有不同初始強(qiáng)度的激波，繼而調(diào)節(jié)內(nèi)錐形激波的匯聚增強(qiáng)過程。在每個(gè)給定的θw工況下，通過選取不同的壁面長度w/R來改變壁面尾緣膨脹波的位置，進(jìn)而產(chǎn)生不同的馬赫盤下游流動(dòng)條件。當(dāng)壁面太短時(shí)，尾緣膨脹波靠近上游，會(huì)對入射的內(nèi)錐形激波產(chǎn)生干擾；而壁面太長時(shí)，內(nèi)收縮比過大，超聲速流場難以建立。因此，w/R取值約束在一定范圍內(nèi)[11,26]。表1列出了本文所有工況中θw和w/R的取值。表1中，w/R最大值和最小值為w/R可取的極限值，在本文中不作為主要討論的工況。

圖1 模型及計(jì)算域Fig. 1 Model and computational domain

表1 數(shù)值計(jì)算工況Table 1 Calculation conditions

鑒于本文主要關(guān)注內(nèi)收縮直錐軸線附近的波系結(jié)構(gòu)，黏性和邊界層效應(yīng)對遠(yuǎn)離壁面的主要激波及其反射結(jié)構(gòu)的影響可以忽略[5,24]，因而通過求解理想氣體（比熱比γ= 1.4）的歐拉方程，對軸對稱內(nèi)收縮直錐流場進(jìn)行CFD （Computational Fluid Dynamics）定常模擬。空間差分采用總變差減?。═VD）的迎風(fēng)格式[29]，數(shù)值通量采用HLLC（Harten-Lax-van Leer-contact）近似黎曼求解器[30-31]計(jì)算。在所有計(jì)算工況下，均采用結(jié)構(gòu)化網(wǎng)格將計(jì)算域離散。計(jì)算域如圖1所示，除壁面和軸線邊界外，其他邊界均采用無反射的遠(yuǎn)場條件[19,24]。遠(yuǎn)場來流條件與文獻(xiàn)[6]中保持一致，馬赫數(shù)Ma∞= 6，靜壓p∞= 0.9 kPa。在計(jì)算過程中監(jiān)測各方程殘差，待殘差下降4～5個(gè)數(shù)量級且波系結(jié)構(gòu)穩(wěn)定后，認(rèn)為流場計(jì)算收斂。

為考核網(wǎng)格無關(guān)性，分別采用單元尺度為Δx/R≈Δr/R≈ 0.002和0.001的兩套網(wǎng)格進(jìn)行計(jì)算，其中Δx和Δr分別為單個(gè)網(wǎng)格沿x和r方向的長度。將獲得的流場進(jìn)行對比，發(fā)現(xiàn)馬赫盤位置沿x方向的變化量僅約為0.001R，兩套網(wǎng)格的主要波系結(jié)構(gòu)幾乎重合。因此，本文所有工況的數(shù)值模擬都采用單元尺度為Δx/R≈ Δr/R≈ 0.001的網(wǎng)格。

為驗(yàn)證本文數(shù)值模擬方法的可靠性，以θw= 10°、w/R= 1的工況為例，與文獻(xiàn)[6]中的實(shí)驗(yàn)結(jié)果進(jìn)行了對比。如圖2所示，將數(shù)值模擬得到的馬赫數(shù)等值線（紅色實(shí)線）疊加在風(fēng)洞實(shí)驗(yàn)[6]紋影照片上可以看出，數(shù)值模擬給出的馬赫盤位置稍靠下游，與實(shí)驗(yàn)結(jié)果相距僅約為0.0008R，兩者的主要激波反射結(jié)構(gòu)基本重合。這也表明本文的數(shù)值模擬結(jié)果是可靠的。

2 結(jié)果與討論

2.1 軸對稱內(nèi)錐流動(dòng)基本結(jié)構(gòu)

鑒于w/R不影響軸線附近的激波反射類型，本節(jié)統(tǒng)一選取w/R= 1.5，并以θw= 10°工況為例，簡要介紹軸對稱內(nèi)錐流動(dòng)中的基本激波結(jié)構(gòu)。

圖3（a）展示了θw= 10°、w/R= 1.5工況的流場馬赫數(shù)（Ma）云圖及等值線，圖3（b）將激波反射結(jié)構(gòu)所在區(qū)域（虛線框Ⅰ）進(jìn)行了局部放大顯示。結(jié)合圖3（a）和圖3（b）可以看出，壁面前緣產(chǎn)生的軸對稱內(nèi)錐形入射激波IS在向下游發(fā)展的過程中持續(xù)地匯聚增強(qiáng)，尤其在接近軸線（r/R= 0）時(shí)出現(xiàn)較為劇烈的彎曲。最終，IS在軸線附近發(fā)生馬赫反射，匯聚過程被馬赫盤MD終結(jié)，同時(shí)在三波點(diǎn)T處形成反射激波RS和剪切層SL（在本文討論中不考慮SL的厚度，將其視為滑移線[11]）。壁面尾緣產(chǎn)生的膨脹波EW會(huì)透射過RS并入射在SL上，對SL產(chǎn)生干擾。圖3（c）將MD下游的局部流場（圖3（a）中虛線框Ⅱ）進(jìn)行了放大顯示，可以看出，在緊鄰MD下游的區(qū)域，SL下方與軸線包裹的流管內(nèi)是亞聲速氣流；隨著SL向下游延伸，流管內(nèi)的氣流逐漸加速至超聲速。換言之，在SL下方與軸線包裹的流管內(nèi)，會(huì)形成聲速喉道。在上、下游流動(dòng)的影響下，SL的形態(tài)在聲速喉道的形成過程中，扮演著重要的角色。

在來流Ma∞= 6條件下，由平面激波反射理論可知[11]，對于氣流偏轉(zhuǎn)角同樣為θw（見表1）的平面激波而言，理論上均不存在馬赫反射。對于軸對稱內(nèi)錐形激波而言，雖然其初始強(qiáng)度等同于氣流偏轉(zhuǎn)角為θw的平面激波（壁面前緣處的波后氣流偏轉(zhuǎn)角完全由壁面約束），但是在匯聚增強(qiáng)的過程中，內(nèi)錐形入射激波IS會(huì)超出von Neumann準(zhǔn)則對應(yīng)的激波強(qiáng)度[11]（稱為von Neumann強(qiáng)度[32]），使得馬赫反射對應(yīng)的入射激波強(qiáng)度條件得以滿足（見圖3）。鑒于馬赫盤位置應(yīng)當(dāng)在入射激波強(qiáng)度達(dá)到von Neumann強(qiáng)度之后，允許馬赫盤存在的空間范圍取決于軸對稱內(nèi)錐形激波的匯聚增強(qiáng)過程。此外，在該空間范圍內(nèi)，三波點(diǎn)需與馬赫盤下游的流動(dòng)條件相匹配，因此馬赫盤的位置受其下游流動(dòng)的影響。

為深入探究影響馬赫盤位置的機(jī)制，需要從兩個(gè)方面入手：首先，通過改變壁面前緣角度，考察激波匯聚增強(qiáng)過程對馬赫盤位置的影響；然后，通過改變壁面長度，進(jìn)一步考察馬赫盤下游流動(dòng)對馬赫盤位置的影響。

2.2 壁面前緣角度對馬赫盤位置的影響

由2.1節(jié)可知，壁面前緣角度θw決定了軸對稱內(nèi)錐形入射激波的初始強(qiáng)度，并通過影響激波匯聚過程來限制允許馬赫盤存在的空間范圍。因此，闡明θw對內(nèi)錐形激波匯聚過程的影響是十分必要的。

根據(jù)譚廉華等[19-20]提出的理論，軸對稱內(nèi)錐形入射激波形狀r=r(x)可以表達(dá)為二階非線性常微分方程式（1）和式（2）：

其中，rx和rxx分別表示r關(guān)于x的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，入射激波后的氣流偏轉(zhuǎn)角θi、壓力pi以及馬赫數(shù)Mai均是入射激波角βi的函數(shù)。對于本文而言，式（1）存在形如式（3）和式（4）的兩個(gè)邊界條件，

其中，壁面前緣處的激波角βw由Ma∞和θw確定?？梢姡葁是影響激波匯聚增強(qiáng)過程的重要參數(shù)。

以w/R= 1.5工況為例，采用Runge-Kutta法[33]對式（1）～式（4）進(jìn)行積分，得到了不同θw下的入射激波形狀，如圖4所示。同時(shí)，在圖4中還疊加了本文數(shù)值模擬獲得的入射激波及馬赫盤形狀。對比圖4中的結(jié)果可知，在不考慮發(fā)生馬赫反射的情況下，譚廉華等[19-20]的理論能夠較為準(zhǔn)確地給出入射激波的形狀。由于式（1）～式（4）是基于入射激波波后流線為直線這一假設(shè)[19-20]推導(dǎo)而來的，而實(shí)際流場中入射激波波后流線存在曲率[21]，因此，圖4中的理論和數(shù)值模擬結(jié)果存在一定的偏差。然而，這種偏差很小，式（1）～式（4）仍不失為一種描述入射激波形狀的有效理論。需要特別注意的是，式（1）～式（4）在軸線處（r/R= 0）存在奇性，在入射激波到達(dá)軸線之前，必然會(huì)發(fā)生馬赫反射并形成馬赫盤[15-16]。由于本節(jié)主要關(guān)注遠(yuǎn)離軸線處的入射激波，不影響對內(nèi)錐形激波匯聚過程的討論。

圖4 不同θw下軸對稱內(nèi)錐形入射激波形狀(w/R = 1.5)Fig. 4 Shapes of axisymmetric internal conical shocks at different θw (w/R = 1.5)

鑒于式（1）～式（4）對內(nèi)錐形激波匯聚增強(qiáng)機(jī)理的反映不夠直觀，因而結(jié)合CST對圖4中的結(jié)果進(jìn)一步分析。針對任意曲面激波，CST方程[21]給出了波后當(dāng)?shù)匮亓骶€的無量綱壓力梯度的表達(dá)式。對于均勻來流下的軸對稱內(nèi)錐形入射激波而言，其波后沿流線的無量綱壓力梯度P2i可表示為式（5）[22]：

式（5）中a1和a2均為與Ma∞和βi有關(guān)的系數(shù)，Sai和Sbi分別為流平面（激波當(dāng)?shù)貋砹骱筒ê髿饬髌D(zhuǎn)方向共同決定的平面）和法平面（與流平面內(nèi)激波當(dāng)?shù)厍邢虼怪钡钠矫妫﹥?nèi)的激波曲率[22]。在壁面前緣處（x/R= 0），可以將式（5）化簡為式（6）[8]：

式（6）中系數(shù)c與初始激波強(qiáng)度正相關(guān)，Sci為激波的初始橫向曲率。由于在壁面前緣處R是確定的（Sci確定），根據(jù)式（6）可知，c越大（θw越大），初始時(shí)激波匯聚增強(qiáng)越快（P2i越大）。

在向下游的匯聚過程中，由于激波波后流動(dòng)方向不完全受壁面約束，式（6）將不再精確適用。因此，根據(jù)式（1）～式（4）計(jì)算得到的激波形狀r=r(x)（見圖4），再結(jié)合式（5）可以理論給出P2i的沿程分布。如圖5所示，在壁面前緣處（x/R= 0），P2i隨著θw的增加而增大，這與式（6）的描述相一致。前文圖4已經(jīng)表明，在下游同一x/R位置處，當(dāng)θw越大時(shí)，入射激波角增大，并且入射激波到軸線的距離更近。這也意味著激波面的曲率Sai和Sbi均增大。相應(yīng)地，根據(jù)式（5）可知P2i也隨著增大（見圖5）。進(jìn)一步表明了隨著θw的增大，內(nèi)錐形激波匯聚增強(qiáng)得更快。換言之，在激波強(qiáng)度和激波面曲率的共同作用下，初始較強(qiáng)的激波（θw大）匯聚增強(qiáng)始終快于初始較弱的激波（θw?。?。

圖5 不同θw下軸對稱內(nèi)錐形激波波后沿流線的無量綱壓力梯度分布(w/R = 1.5)Fig. 5 Distributions of the normalized pressure gradient behind the axisymmetric internal conical shocks at different θw (w/R = 1.5)

為了認(rèn)識內(nèi)錐形激波匯聚過程中激波自身強(qiáng)度的變化，如圖6所示，在數(shù)值模擬流場中沿著入射激波前緣到三波點(diǎn)提取了激波前后的壓比pi/p∞，并利用譚廉華等[19-20]的理論，由式（1）～式（4）給出了pi/p∞沿x/R的分布。圖6中理論和數(shù)值模擬獲得的pi/p∞分布整體符合較好，在接近三波點(diǎn)時(shí)，由于激波波后流線曲率較大[21]，理論結(jié)果存在一定的偏差。隨著下游激波強(qiáng)度的持續(xù)增大，在內(nèi)錐形激波的匯聚過程中，必然會(huì)達(dá)到von Neumann強(qiáng)度，圖6中用虛線給出了本文Ma∞= 6條件下von Neumann準(zhǔn)則對應(yīng)的激波壓比pN/p∞≈ 9.67。在不同的θw下，內(nèi)錐形激波匯聚增強(qiáng)的快慢不同（見圖5），導(dǎo)致入射激波上對應(yīng)于von Neumann強(qiáng)度（pi/p∞=pN/p∞）的位置存在差異。在給定的θw下，根據(jù)式（1）～式（4）描述的入射激波形狀，可以得到入射激波上von Neumann強(qiáng)度的對應(yīng)點(diǎn)（xN/R，rN/R）。

圖6 不同θw下軸對稱內(nèi)錐形激波壓比的沿程分布(w/R = 1.5)Fig. 6 Pressure ratio distributions of the axisymmetric internal conical shocks at different θw (w/R = 1.5)

如圖7（a）所示，在θw= 7°～16°范圍內(nèi)，利用式（1）～式（4）分別計(jì)算出了xN/R和rN/R隨θw的變化曲線，并且標(biāo)記出了本文數(shù)值模擬流場中各θw工況對應(yīng)的xN/R和rN/R位置。從圖7（a）可以看出，隨著θw的增大，xN/R減小，而rN/R迅速增大。進(jìn)一步地，將這些由式（1）～式（4）計(jì)算出的（xN/R，rN/R），在（x/R，r/R）平面內(nèi)依次展示，可以得到如圖7（b）所示的曲線。鑒于當(dāng)入射激波強(qiáng)度達(dá)到von Neumann強(qiáng)度時(shí)，馬赫反射才有可能發(fā)生，因而圖7（b）中的（xN/R，rN/R）曲線和軸線（r/R= 0）共同限定了不同θw下允許馬赫盤存在的空間范圍。換言之，式（1）～式（4）為從理論上確定馬赫盤的位置提供了一種可能的途徑。

圖7 不同θw下von Neumann強(qiáng)度對應(yīng)點(diǎn)和三波點(diǎn)(w/R = 1.5)Fig. 7 Point of von Neumann strength and triple point at different θw (w/R = 1.5)

為對比三波點(diǎn)與（xN/R，rN/R）之間的位置關(guān)系，從數(shù)值模擬流場中提取了不同θw下的三波點(diǎn)位置和von Neumann強(qiáng)度對應(yīng)點(diǎn)的坐標(biāo)，疊加在圖7（b）中。通過數(shù)值模擬結(jié)果可以看出，在不同的θw工況下，三波點(diǎn)均被約束在von Neumann強(qiáng)度對應(yīng)點(diǎn)（xN/R，rN/R）的下游，兩者相距很近，但并不重合。隨著θw增大，軸對稱內(nèi)錐形激波更快地達(dá)到von Neumann強(qiáng)度，但馬赫盤的位置并非完全依賴于激波匯聚增強(qiáng)過程。由此可見，僅通過式（1）～式（4）結(jié)合CST方程難以完全確定馬赫盤的位置，還需要考慮下游流動(dòng)條件。

2.3 馬赫盤下游流動(dòng)對其位置的影響

針對馬赫盤下游的流動(dòng)問題，一般可以將滑移線SL包裹的流管（見圖3）視為準(zhǔn)一維流動(dòng)[11,34-36]。根據(jù)馬赫盤面積Am和流管內(nèi)聲速喉道面積As，可以建立形如式（7）的等熵關(guān)系[11,26]：

式（7）中為緊鄰馬赫盤下游的氣流平均馬赫數(shù)。馬赫盤位置與其下游流動(dòng)之間的匹配，往往可以通過聲速喉道來實(shí)現(xiàn)。

對于平面激波反射問題，已經(jīng)發(fā)展出相應(yīng)的幾何模型[11,26]來描述馬赫桿下游的流動(dòng)。如圖8所示，起初從三波點(diǎn)T′發(fā)出的滑移線SL′具有恒定的角度（即三波點(diǎn)后氣流偏轉(zhuǎn)角）。緊接著，壁面尾緣產(chǎn)生的膨脹波穿過反射激波RS′后，作用在SL′上，使得SL′的角度連續(xù)減小（D′A′段），進(jìn)而形成了聲速喉道。將這一幾何模型及相應(yīng)的二維激波/膨脹波關(guān)系與式（7）相結(jié)合，即可確定馬赫桿MS′的位置[11,26]。這種方法為軸對稱內(nèi)錐形激波反射問題中馬赫盤位置的確定，提供了極大的啟示。然而，在軸對稱內(nèi)錐形激波反射流場中，馬赫盤下游的流動(dòng)在受到壁面尾緣膨脹波干擾之前，先需要匹配反射激波后的非均勻壓力，導(dǎo)致滑移線經(jīng)歷的變化比二維情況要復(fù)雜得多。因此，需要對軸對稱內(nèi)錐形激波反射流場中馬赫盤下游流動(dòng)進(jìn)行深入考察。

圖8 平面激波反射示意圖Fig. 8 Schematic of planar shock reflection

2.3.1 馬赫盤下游流動(dòng)特征

首先關(guān)注壁面尾緣膨脹波干擾之前，馬赫盤下游的流動(dòng)特征。對于不同的θw，當(dāng)w/R充分大時(shí)（見表1），下游壁面尾緣產(chǎn)生的膨脹波離馬赫盤較遠(yuǎn)，使得滑移線SL包裹的流管充分發(fā)展。根據(jù)流管自身能否演化形成聲速喉道，將馬赫盤下游流動(dòng)劃分為兩類。

第一類工況：在壁面尾緣膨脹波干擾之前，馬赫盤下游滑移線SL包裹的流管能夠形成聲速喉道。θw= 8°或10°且w/R充分大時(shí)，均屬于此類工況。以θw= 10°、w/R= 3.5工況為例進(jìn)行詳細(xì)分析。圖9展示了該工況下的流場局部無量綱壓力（p/p∞）云圖及其等值線，可以看到，入射激波IS和反射激波RS的波后壓力都是不均勻的。以穿過三波點(diǎn)T1的流線表征馬赫盤MD下游的滑移線SL，并以MD下游軸線（r/R= 0）上的馬赫數(shù)表征流管內(nèi)的馬赫數(shù)分布，量化分析馬赫盤下游的流動(dòng)特征。

圖9 θw = 10°、w/R = 3.5工況數(shù)值模擬流場的局部無量綱壓力云圖Fig. 9 Local normalized pressure contours of the numerical flow field at θw = 10° and w/R = 3.5

圖10給出了SL包裹的流管面積A= πr2與馬赫盤面積Am之比的沿程變化情況，其中，xm為馬赫盤MD在軸線處的x坐標(biāo)。同時(shí)，圖10中還疊加了流管內(nèi)的馬赫數(shù)（Ma）分布，并標(biāo)明了流管內(nèi)首個(gè)聲速喉道的位置。結(jié)合圖9和圖10可知，SL包裹的流管從T1開始快速收縮，使得緊鄰馬赫盤下游的亞聲速氣流加速。在A1點(diǎn)處，流管內(nèi)形成了聲速喉道。值得注意的是，在該工況下，壁面尾緣膨脹波與SL干擾的起始點(diǎn)D1不僅遠(yuǎn)離聲速喉道A1，而且兩者之間還有超聲速區(qū)相阻隔?？梢?，聲速喉道A1的形成原因，是獨(dú)立于壁面尾緣膨脹波的。這與圖8中的平面激波反射流場存在本質(zhì)差異。

圖10 θw = 10°、w/R = 3.5工況流管面積變化及其內(nèi)馬赫數(shù)分布Fig. 10 Area ratio of the stream tube and distribution of Mach number in the stream tube at θw = 10° and w/R = 3.5

進(jìn)一步地，圖11給出了反射激波RS（從T1點(diǎn)到F1點(diǎn)）波前、波后的壓力以及滑移線SL上（從 T1點(diǎn)到D1點(diǎn)）的壓力分布，以便于通過分析馬赫盤下游流動(dòng)如何與反射激波后的非均勻壓力場相匹配，進(jìn)而揭示馬赫盤下游滑移線SL包裹的流管內(nèi)聲速喉道的形成原因。對比圖11和圖6可以看出，RS的波前壓力與入射激波IS波后壓力的變化趨勢并不一致。這是因?yàn)樵谳S對稱內(nèi)收縮流場中，IS波后不同位置的氣流在到達(dá)RS之前，經(jīng)歷了不同程度的等熵壓縮過程（見圖9）。這些氣流經(jīng)過RS再次壓縮后，壓力呈先降低（見圖11的T1E1段）再升高（見圖11的E1F1段）的變化趨勢，而RS波后非均勻分布的壓力將直接影響馬赫盤下游的滑移線SL及流管形態(tài)。

圖11 θw = 10°、w/R = 3.5工況反射激波前后及滑移線上無量綱壓力分布Fig. 11 Pre-shock and post-shock normalized pressure distribution of the reflected shock, and normalized pressure distribution along the slip line at θw = 10° and w/R = 3.5

從圖11可以看出，滑移線SL上T1B1段的壓力連續(xù)下降，以匹配RS在T1E1段的波后壓力；而滑移線SL上B1D1段的壓力連續(xù)上升，以匹配RS在E1F1段的波后壓力。結(jié)合圖10和圖11可知，在滑移線SL上T1B1段的降壓過程中，SL包裹的流管內(nèi)原本亞聲速的氣流連續(xù)加速至超聲速。換言之，流管內(nèi)聲速喉道的形成是匹配RS波后非均勻壓力的結(jié)果。RS和SL之間（見圖9中T1D1F1區(qū)域）通過等熵膨脹或壓縮來實(shí)現(xiàn)壓力變化，盡管這種變化關(guān)系不如二維平面激波反射流場那樣清晰，但從理論上建立分析方法仍然是將來值得嘗試的研究課題。

第二類工況：在壁面尾緣膨脹波干擾之前，馬赫盤下游滑移線SL包裹的流管無法形成聲速喉道。當(dāng)θw增加至12°或更大（θw= 15°）且w/R充分大時(shí)，均屬于此類工況。圖12以θw= 12°、w/R= 3.0工況為例，展示了流場的無量綱壓力云圖。圖13給出了滑移線SL包裹的流管面積的沿程變化情況(A/Am)以及流管內(nèi)的馬赫數(shù)分布。圖14展示了反射激波RS（從T2點(diǎn)到F2點(diǎn)）波前、波后的壓力以及滑移線SL上（從 T2點(diǎn)到A2點(diǎn)）的壓力分布。

圖12 θw = 12°、w/R = 3.0工況數(shù)值模擬流場的局部無量綱壓力云圖Fig. 12 Local normalized pressure contours of the numerical flow field at θw = 12° and w/R = 3.0

圖13 θw = 12°、w/R = 3.0工況流管面積變化及其內(nèi)馬赫數(shù)分布Fig. 13 Area ratio of the stream tube and distribution of Mach number in the stream tube at θw = 12° and w/R = 3.0

與前文θw= 10°工況類似，馬赫盤下游流管內(nèi)的壓力也需要匹配RS的波后壓力。有所不同的是，θw= 12°時(shí)RS在T2E2段的波后壓力（見圖14）不像θw= 10°工況那樣迅速降低（見圖11）。因而θw= 12°時(shí)馬赫盤下游流管內(nèi)氣流的加速更為緩慢，以至于在壓力谷值點(diǎn)B2處，氣流仍未達(dá)到聲速（見圖13）。在B2點(diǎn)下游，由于要匹配RS在E2F2段波后壓力的升高，流管內(nèi)的亞聲速氣流連續(xù)減速增壓（B2G2段），直至壁面尾緣膨脹波在D2點(diǎn)處對滑移線SL產(chǎn)生干擾。在馬赫盤下游流管內(nèi)的亞聲速流動(dòng)中，由于下游擾動(dòng)能夠傳到上游，在D2點(diǎn)上游不遠(yuǎn)處的G2點(diǎn)已經(jīng)能夠感受到壁面尾緣膨脹波對SL的影響。由膨脹波帶來的壓力降低（見圖14），使得流管內(nèi)的氣流急劇加速，直至在A2點(diǎn)處形成聲速喉道（見圖13）?？梢姡葁= 12°時(shí)馬赫盤下游流管內(nèi)的氣流無法通過匹配RS的波后壓力達(dá)到聲速，而聲速喉道的形成依賴于壁面尾緣膨脹波對SL的干擾。

圖14 θw = 12°、w/R = 3.0工況反射激波前后及滑移線上無量綱壓力分布Fig. 14 Pre-shock and post-shock normalized pressure distribution of the reflected shock, and normalized pressure distribution along the slip line at θw = 12° and w/R = 3.0

綜上所述，在兩類馬赫盤下游流動(dòng)中，反射激波RS波后壓力下降的快慢程度不同，導(dǎo)致馬赫盤下游流管內(nèi)聲速喉道的形成過程存在是否依賴于壁面尾緣膨脹波的差異性。由于馬赫盤與聲速喉道之間具有如式（7）所示的平衡關(guān)系，馬赫盤位置對壁面尾緣膨脹波的依賴性也將不同。因此，需進(jìn)一步考察壁面尾緣膨脹波對馬赫盤位置的影響。

2.3.2 膨脹波對馬赫盤位置的影響

在給定壁面前緣角度θw和入口半徑R的條件下，w/R決定了流場的內(nèi)收縮比以及壁面尾緣膨脹波的存在區(qū)域。在本文關(guān)注的w/R范圍內(nèi)，壁面尾緣膨脹波對馬赫盤下游流動(dòng)產(chǎn)生干擾的位置隨著w/R的增加而向下游移動(dòng)。

以θw= 10°為例，不同w/R工況下的入射激波及馬赫盤形狀如圖15所示。從圖15中馬赫盤局部區(qū)域（虛線框）的放大圖可以看出，不同w/R工況下，入射激波的形狀均重合，但馬赫盤位置存在差異。這種差異與2.3.1節(jié)討論的聲速喉道的形成機(jī)制緊密相關(guān)。

圖15 不同w/R下軸對稱內(nèi)錐形入射激波及馬赫盤形狀(θw = 10°)Fig. 15 Shapes of axisymmetric internal conical shocks and Mach disks at different w/R (θw = 10°)

圖16提取了θw= 10°不同w/R工況中馬赫盤下游流管內(nèi)的馬赫數(shù)分布，其中T1標(biāo)明了三波點(diǎn)的位置。以壁面長度充分大的w/R= 3.5工況為參照，A1點(diǎn)標(biāo)明了其聲速喉道的位置，而其他w/R工況下流管內(nèi)的馬赫數(shù)分布從T1下游不同位置開始與w/R= 3.5工況產(chǎn)生差異。如2.3.1節(jié)所述，這種差異是由壁面尾緣膨脹波對滑移線產(chǎn)生干擾而導(dǎo)致的。在w/R= 1.1～2.0范圍內(nèi)，壁面尾緣膨脹波對滑移線產(chǎn)生干擾的位置（圖16中三角符號）在A1點(diǎn)上游，并輔助馬赫盤下游流管迅速形成聲速喉道（圖16中矩形符號）。在w/R= 2.5～3.7范圍內(nèi)，壁面尾緣膨脹波對滑移線產(chǎn)生干擾的位置在A1點(diǎn)下游，此時(shí)馬赫盤下游流管中聲速喉道形成于A1點(diǎn)，不依賴于壁面尾緣膨脹波。

圖16 θw = 10°工況不同w/R下的流管內(nèi)馬赫數(shù)分布Fig. 16 Distributions of Mach number in the stream tubes with different w/R at θw = 10°

進(jìn)一步地，在θw= 10°工況下圖17（a）分別以w/R= 3.5時(shí)的馬赫盤位置xm*/R和聲速喉道位置xs*/R為參考值，給出了不同w/R工況的馬赫盤位置和聲速喉道位置的變化情況?？梢钥闯觯趙/R=1.1～2.0范圍內(nèi)，隨著w/R增大，聲速喉道明顯向下游移動(dòng)；同時(shí)，馬赫盤向上游移動(dòng)。需要注意的是，馬赫盤位置的變化幅度遠(yuǎn)小于聲速喉道位置的變化幅度。馬赫盤的位置僅在w/R＜ 1.5時(shí)變化明顯，在w/R＞ 1.5之后變化微小。結(jié)合前文2.2節(jié)可知，由于在軸線附近入射激波的形狀將越來越陡峭，馬赫盤位置的明顯移動(dòng)，必然引起馬赫盤面積的大幅變化。

圖17 θw = 10°工況不同w/R下的馬赫盤、聲速喉道的位置和面積Fig. 17 Positions and areas of Mach disks and sonic throats with different w/R at θw = 10°

圖17（b）以w/R= 3.5時(shí)的馬赫盤面積Am*和聲速喉道面積As*為參考值，給出了不同w/R工況下的馬赫盤面積和聲速喉道面積的變化情況?？梢钥闯?，由于兩者遵循式（7）的關(guān)系，兩者的變化曲線幾乎重合。在w/R＜ 1.5時(shí)，馬赫盤面積和聲速喉道面積顯著減?。辉趙/R= 1.5～2.5范圍內(nèi)，盡管聲速喉道的位置隨著w/R增大而向下游大幅移動(dòng)（見圖17（a）），其面積變化卻很?。ㄒ妶D17（b））。這主要是由于在w/R= 1.5～2.0范圍內(nèi)，當(dāng)壁面尾緣膨脹波入射滑移線時(shí)，馬赫盤下游流管內(nèi)的氣流已經(jīng)接近聲速，流管的微小面積變化也會(huì)使得流管內(nèi)快速形成聲速喉道（見圖16）。在w/R= 2.5～3.7范圍內(nèi)，聲速喉道的形成不依賴壁面尾緣膨脹波，因此，從圖17（a）和圖17（b）中可以看出，聲速喉道及馬赫盤的面積和位置，均基本保持恒定。類似地，在θw= 8°、w/R= 1.5～4.3工況下，壁面尾緣膨脹波不影響聲速喉道的形成。

對于θw= 12°和15°工況，在本文研究的w/R范圍內(nèi)（見表1），馬赫盤下游流管內(nèi)聲速喉道的形成均依賴于壁面尾緣膨脹波。圖18以θw= 12°工況為例，給出了不同w/R工況馬赫盤下游流管內(nèi)的馬赫數(shù)分布。圖19（a）和圖19（b）分別展示了不同w/R工況馬赫盤、聲速喉道的位置和面積，其中，參考值xm*、Am*和As*均選自w/R= 3.0工況。結(jié)合圖18和圖19（a）可以看出，隨著w/R增加，壁面尾緣膨脹波入射位置向下游移動(dòng)，聲速喉道的位置也隨之向下游移動(dòng)，幾乎呈線性變化。然而，馬赫盤位置的變化幅度遠(yuǎn)小于聲速喉道位置的變化幅度，馬赫盤向上游移動(dòng)的趨勢逐漸減緩，在w/R＞ 1.5之后變化微小。從圖19（b）可以看出，馬赫盤和聲速喉道的面積同樣在w/R＞ 1.5之后變化很小。其原因與θw= 10°、w/R= 1.5～2.0工況下馬赫盤、聲速喉道面積趨于穩(wěn)定的原因類似。

圖18 θw = 12°工況不同w/R下的流管內(nèi)馬赫數(shù)分布Fig. 18 Distributions of Mach number in the stream tubes with different w/R at θw = 12°

圖19 θw = 12°工況不同w/R下的馬赫盤、聲速喉道的位置和面積Fig. 19 Positions and areas of Mach disks and sonic throats with different w/R at θw = 12°

在本文研究的軸對稱內(nèi)錐流場中，馬赫盤下游聲速喉道的形成機(jī)制存在是否依賴于壁面尾緣膨脹波的兩類情況。對于不依賴壁面尾緣膨脹波的情況（如θw= 8°、w/R= 1.5～4.3工況以及θw= 10°、w/R=2.5～3.7工況），馬赫盤下游流管通過與反射激波后的非均勻壓力相匹配，進(jìn)而形成聲速喉道；對于依賴壁面尾緣膨脹波的情況（如θw= 8°、w/R= 1.2工況，θw= 10°、w/R= 1.1～2.0工 況 以及θw= 12°和15°工況），馬赫盤下游流管還需要進(jìn)一步匹配壁面尾緣膨脹波引起的壓力降低，才能夠形成聲速喉道。這種匹配關(guān)系，為理論求解馬赫盤位置提供了啟示。

2.4 馬赫盤位置的理論求解

結(jié)合2.2節(jié)和2.3節(jié)的分析可知，馬赫盤位于入射激波上von Neumann強(qiáng)度對應(yīng)位置所限定的范圍內(nèi)，而馬赫盤下游的流動(dòng)還需要同時(shí)滿足反射激波下游的壓力變化以及準(zhǔn)一維等熵關(guān)系式（7）。換言之，確定馬赫盤的位置需要求解反射激波下游的超聲速流場。最近，Shi等[25]發(fā)展的彎曲激波特征線方法（Method of curved-shock characteristics，MOCC）具有高效且精確求解超聲速流場的優(yōu)勢。本節(jié)在給定來流和壁面前緣入口條件下，借助MOCC求解馬赫盤位置。

首先，假定馬赫盤處于入射激波上的von Neumann強(qiáng)度對應(yīng)點(diǎn)（見2.2節(jié)）。利用MOCC，可以準(zhǔn)確求解入射激波從壁面前緣到三波點(diǎn)的形狀，而且可以進(jìn)一步獲得入射激波下游的超聲速流場，包括入射激波與反射激波之間的流場、反射激波形狀及其波后流場?；凭€作為反射激波波后邊界上的流線，其形狀在利用MOCC的求解過程中可以確定[37]。根據(jù)滑移線形狀，計(jì)算出滑移線包裹的流管面積，并找到流管面積的第一個(gè)最小值A(chǔ)sc。

上述滑移線形狀僅是基于反射激波下游的流場得到的，而反射激波下游壓力變化還需要與馬赫盤下游流管內(nèi)的壓力相匹配。因此，利用準(zhǔn)一維流動(dòng)的等熵關(guān)系式（7），驗(yàn)證馬赫盤下游流動(dòng)是否符合2.3.2節(jié)中的描述。由于事先假定了馬赫盤位置，式（7）中的馬赫盤面積Am及流管內(nèi)初始馬赫數(shù)均為已知量[19]，利用式（7）可以得到聲速喉道面積的理論值A(chǔ)st。將Asc與Ast進(jìn)行對比，若兩者之差在設(shè)定的誤差范圍內(nèi)，則認(rèn)為Asc是準(zhǔn)確的，表明馬赫盤下游流動(dòng)既能匹配反射激波波后的非均勻壓力變化，又滿足準(zhǔn)一維流動(dòng)的等熵關(guān)系，即假定的馬赫盤位置正確；若在設(shè)定的誤差范圍之外，則需要將假定的馬赫盤位置沿x方向調(diào)整后，重新計(jì)算一遍滑移線包裹的流管面積變化，并再次對比新獲得的Asc與Ast的值。如此迭代求解，直至最終得到馬赫盤的位置。

圖20以θw= 12°工況為例，利用上述理論方法求解了馬赫盤位置隨w/R的變化情況。作為對比，圖20中還給出了CFD流場中的馬赫盤位置?？梢钥闯?，理論與CFD結(jié)果吻合良好，兩者之間的最大差異僅為0.002。這表明，本節(jié)的理論方法，能夠準(zhǔn)確地預(yù)測出馬赫盤隨w/R的變化。鑒于MOCC在理論求解效率和精度方面的獨(dú)特優(yōu)勢，可望在更寬的來流條件和內(nèi)錐幾何參數(shù)范圍內(nèi)，快速預(yù)測馬赫盤的位置。

圖20 θw = 12°工況不同w/R下馬赫盤位置的理論和數(shù)值對比Fig. 20 Theoretical and numerical positions of Mach disks with different w/R at θw = 12°

3 結(jié) 論

在來流Ma∞= 6條件下，采用無黏數(shù)值模擬研究了軸對稱直內(nèi)錐流動(dòng)中影響馬赫盤位置的因素，主要得到了以下結(jié)論：

1）軸對稱內(nèi)錐形入射激波自身的匯聚增強(qiáng)過程影響馬赫盤的位置。結(jié)合入射激波形狀理論與彎曲激波理論，闡明了隨著壁面前緣角度θw增大，入射激波的匯聚增強(qiáng)越快，并且能夠更快地達(dá)到von Neumann強(qiáng)度，該強(qiáng)度的對應(yīng)位置限定了允許馬赫盤存在的空間范圍。

2）馬赫盤下游聲速喉道的形成機(jī)制存在是否依賴于壁面尾緣膨脹波的差異性，進(jìn)而對馬赫盤位置產(chǎn)生不同的影響。當(dāng)聲速喉道的形成不依賴于壁面尾緣膨脹波時(shí)，聲速喉道及馬赫盤的位置和面積均不隨w/R增大而變化。當(dāng)聲速喉道的形成依賴于壁面尾緣膨脹波時(shí)，隨著w/R的增加，聲速喉道的位置向下游大幅移動(dòng)，但聲速喉道及馬赫盤的面積逐漸趨于穩(wěn)定。

3）將彎曲激波特征線方法與馬赫盤下游流管的準(zhǔn)一維流動(dòng)模型相結(jié)合，可以快速求解馬赫盤的位置。