湯祎麒， 施崇廣， 鄭曉剛， 朱呈祥， 尤延鋮

(廈門大學(xué)航空航天學(xué)院， 福建廈門 361101)

引 言

21世紀(jì)以來， 高超聲速飛行器逐漸成為全球航空航天工業(yè)的前沿陣地。作為核心結(jié)構(gòu)的進(jìn)氣道， 在捕獲足夠多空氣的同時(shí)要實(shí)現(xiàn)高效的壓縮， 保證整個(gè)推進(jìn)系統(tǒng)產(chǎn)生足夠的推力來滿足寬Mach數(shù)的工作范圍。有研究表明， 對Ma=5～7的碳?xì)淙剂巷w行器， 每提高1%的進(jìn)氣道壓縮效率， 就可以增加3%～5%的推進(jìn)系統(tǒng)比沖[1]。20世紀(jì)40年代， 國內(nèi)外研究者提出了高流量捕獲、 低外部阻力、 寬Mach數(shù)工作范圍的三維內(nèi)收縮式進(jìn)氣道， 相比于傳統(tǒng)的二元式[2]、 軸對稱式[3]和側(cè)壓式進(jìn)氣道[4]， 壓縮效率顯著提高。這類進(jìn)氣道往往通過在基準(zhǔn)流場上流線追蹤得到設(shè)計(jì)型面， 基準(zhǔn)流場直接決定了進(jìn)氣道的氣動(dòng)性能。1942年， Busemann[5]提出了一種錐形內(nèi)收縮流場， 包含一簇等熵壓縮波和一道末端錐形激波， 總壓恢復(fù)較高。但過長壁面和完全內(nèi)收縮的特性導(dǎo)致其無法在低Mach數(shù)下啟動(dòng)， M?lder等[6-8]遂采用截短的型面降低黏性損失， 并通過流線追蹤提升起動(dòng)特性。岳連捷等[9]利用NURBS樣條曲線和ICFA(internal conical flow A)流場組合出總壓恢復(fù)最大的基準(zhǔn)流場。為提高流場出口均勻度， 郭軍亮等[10]基于Busemann流場與ICFA流場組合設(shè)計(jì)了ICFC(internal conical flow C)和ICFD(internal conical flow D)流場。Matthews等[11]基于特征線法(method of characteristics， MOC)設(shè)計(jì)了等內(nèi)收縮角、 等壓力的軸對稱基準(zhǔn)流場。南向軍等[12]結(jié)合特征線法反設(shè)計(jì)了壁面壓力可控的基準(zhǔn)流場， 對得到的進(jìn)氣道開展風(fēng)洞試驗(yàn)[13-14]。朱偉等[15]基于此提出了壁面Mach數(shù)分布可控的基準(zhǔn)流場設(shè)計(jì)方法。李永洲等[16]設(shè)計(jì)了一種中心體可控的雙波入射內(nèi)收縮基準(zhǔn)流場。Malo-Molina等[17-19]采用能夠在偏航和俯仰平面對氣流雙重壓縮的三維無黏流場設(shè)計(jì)了一種新型的Jaws進(jìn)氣道。劉燚[20]和方興軍[21]設(shè)計(jì)了由出口Mach數(shù)或速度分布計(jì)算流場的特征線法， 韓偉強(qiáng)等[22]據(jù)此提出了由給定下游邊界流場參數(shù)逆向求解基準(zhǔn)流場的特征線法， 為基準(zhǔn)流場的求解提供了新思路。衛(wèi)鋒等[23]發(fā)展了可實(shí)現(xiàn)消波的基準(zhǔn)流場設(shè)計(jì)方法。美國Astrox公司的Kothari等[24]基于三維內(nèi)收縮軸對稱基準(zhǔn)流場設(shè)計(jì)得到了Funnel型進(jìn)氣道。隨著高超聲速飛行器逐漸向一體化發(fā)展， 幾何過渡變截面內(nèi)收縮進(jìn)氣道逐漸興起。Smart等[25-26]基于截短角度的倒置等熵噴管流場設(shè)計(jì)了矩形進(jìn)口、 橢圓出口的內(nèi)收縮進(jìn)氣道。譚慧俊等[27]發(fā)展了基于樣條曲面的新型內(nèi)通道設(shè)計(jì)方法， 不會(huì)對入射波依賴域內(nèi)的型面產(chǎn)生影響。Gollan等[28]綜合典型的設(shè)計(jì)方法設(shè)計(jì)了非常規(guī)進(jìn)口、 橢圓出口的模塊化進(jìn)氣道， 可用于錐形前體的一體化設(shè)計(jì)。上述方法只能達(dá)到幾何上的變截面設(shè)計(jì)， 為實(shí)現(xiàn)設(shè)計(jì)狀態(tài)的全流量捕獲， 研究者逐漸把目光轉(zhuǎn)向氣動(dòng)過渡， 為復(fù)雜乘波外形與內(nèi)收縮進(jìn)氣道的一體化提供了可能。洛克西德·馬丁公司設(shè)計(jì)的FALCON飛行器[29]采用了一種變截面的翼身融合進(jìn)氣道方案。尤延鋮等[30]采用吻切軸對稱理論基于非對稱內(nèi)收縮基準(zhǔn)流場設(shè)計(jì)了變截面內(nèi)乘波進(jìn)氣道， 可同時(shí)控制進(jìn)、 出口截面形狀。尤延鋮等[31]還基于由部分Busemann流場和ICFA流場組成的基準(zhǔn)流場提出了激波形狀可以任意選取的變截面內(nèi)乘波式進(jìn)氣道的設(shè)計(jì)方法。Xiao等[32]采用等收縮比的方法設(shè)計(jì)了一種新型變截面內(nèi)乘波進(jìn)氣道， 通過較均勻的喉道截面壓力分布實(shí)現(xiàn)了隔離段內(nèi)的消波。

內(nèi)收縮基準(zhǔn)流場正朝著精細(xì)化設(shè)計(jì)的趨勢發(fā)展， 但前文提到的基準(zhǔn)流場及其設(shè)計(jì)方法仍存在一定缺陷： 一方面反設(shè)計(jì)用到的特征線法無法獲得高階的氣動(dòng)參數(shù)， 計(jì)算耗時(shí)較長， 不利于對流場的優(yōu)化求解； 另一方面絕大部分流場是由一道入射激波和一道反射激波組成的雙波結(jié)構(gòu)， 激波壓縮比例小， 且設(shè)計(jì)時(shí)需要給定壁面型線或沿程氣動(dòng)參數(shù)， 無法直接控制波系分布。Shi等[33]提出的彎曲激波特征線法(method of curved-shock characteristics， MOCC)能很好地解決前一個(gè)問題： 對于同樣數(shù)量的網(wǎng)格節(jié)點(diǎn)， 相同精度下MOCC相比于MOC可以節(jié)省約90%的計(jì)算時(shí)間， 而優(yōu)化問題往往需要求解上萬個(gè)流場， 這時(shí)MOCC可以節(jié)省大約90%的計(jì)算資源。針對后一個(gè)問題， 在保持其他參數(shù)不變的情況下通過多級壓縮可以在滿足增壓比和總壓恢復(fù)的平衡下提高流場在高型面設(shè)計(jì)Mach數(shù)時(shí)的壓縮效率。考慮到不同形狀的激波曲線會(huì)對進(jìn)氣道的幾何形狀和氣動(dòng)性能以至發(fā)動(dòng)機(jī)流量產(chǎn)生很大影響， 本文基于MOCC提出了一種激波型線和出口參數(shù)分布可控的兩級壓縮內(nèi)收縮基準(zhǔn)流場的反設(shè)計(jì)方法， 借助數(shù)值模擬驗(yàn)證了設(shè)計(jì)方法對任意彎曲波系結(jié)構(gòu)的準(zhǔn)確性。隨后改變設(shè)計(jì)參數(shù)， 探索流場氣動(dòng)特征與性能參數(shù)的變化規(guī)律， 以及不同設(shè)計(jì)參數(shù)間的相互影響， 對涉及的關(guān)鍵問題進(jìn)行研究。該方法的提出有望對后續(xù)任意三維激波的內(nèi)乘波流場及進(jìn)氣道的設(shè)計(jì)與優(yōu)化提供理論支撐。

1 基準(zhǔn)流場設(shè)計(jì)方法

1.1 彎曲激波理論及彎曲激波特征線法

彎曲激波特征線法(MOCC)是基于彎曲激波理論建立的。該理論由M?lder提出， 1階形式如下

其中， 系數(shù)A，B，C，E，G是來流條件和激波角的表達(dá)式， 具體可參考文獻(xiàn)[34]。下標(biāo)1, 2分別代表波前和波后的流動(dòng)條件,Sa和Sb分別為激波在流動(dòng)面和流動(dòng)垂直面的曲率, Shi等將方程兩邊對σ求導(dǎo)， 得到2階彎曲激波方程， 方程及系數(shù)的推導(dǎo)詳見文獻(xiàn)[35]。

A″1P′1+B″1D′1+E″1?！?=
A″2P′2+B″2D′2+C″S′a+G″S′b+const″

A?1P′1+B?1D′1+E?1?！?=
A?2P′2+B?2D′2+C?S′a+G?S′b+const?

當(dāng)來流條件已知， 通過彎曲激波理論可以求解彎曲激波波后氣動(dòng)參數(shù)及梯度， 但其作為高階R-H理論僅適于激波附近區(qū)域。Shi等提出了MOCC來擴(kuò)展其應(yīng)用范圍， MOCC采用流線-特征線坐標(biāo)系， 方程(1)～(6)組成了此法的偏微分方程

壓力梯度

P=(?p/?s)/(ρV2)

(1)

流線曲率

D=?δ/?s

(2)

沿流線的動(dòng)量方程

(3)

沿流線的能量方程

(4)

沿特征線方向的動(dòng)量方程

=ρV2(Pcosμ?Dsinμ)

(5)

沿特征線方向的流量方程

(6)

由于氣動(dòng)參數(shù)沿流線/特征線偏微分等價(jià)于全微分， 故可將上述控制方程轉(zhuǎn)換為全微分方程

(7)

對1階常微分方程組(7)兩邊沿特征線或流線方向求導(dǎo)可得2階常微分方程組。

利用MOCC， 彎曲激波理論可應(yīng)用于彎曲激波波后全流場氣動(dòng)參數(shù)及導(dǎo)數(shù)的求解。不同于特征線法需要1條流線和2條特征線迭代求解相容方程， MOCC僅需要1條流線和1條特征線迭代求解控制方程， 消除了計(jì)算誤差的累積。通過建立壓力梯度、 流線曲率與氣動(dòng)參數(shù)的直接聯(lián)系， 合理利用梯度信息等高階變量， 不僅可以減少網(wǎng)格數(shù)量、 提高計(jì)算效率， 還可以用來研究相關(guān)物理現(xiàn)象， 具有高精度的優(yōu)勢。

1.2 兩級壓縮內(nèi)收縮基準(zhǔn)流場反設(shè)計(jì)

鑒于MOCC在求解平面和軸對稱流場中的高精度， 可用于超聲速流場的反設(shè)計(jì)。進(jìn)氣道設(shè)計(jì)分為兩個(gè)部分： 1)確定無黏超聲速流場； 2)采用流線追蹤技術(shù)設(shè)計(jì)無黏進(jìn)氣道。第一步是設(shè)計(jì)的核心， 因?yàn)檫M(jìn)氣道性能與超聲速流場直接相關(guān)， 而入射和反射波又直接決定了流場的性能。所以在已知激波或出口參數(shù)的情況下， 可在生成進(jìn)氣道之前算出無黏進(jìn)氣道的總壓損失?？紤]到實(shí)際應(yīng)用價(jià)值， 超聲速流場的反設(shè)計(jì)至關(guān)重要。由于激波與能量損失有著直接關(guān)系， 因此在理論設(shè)計(jì)中可以通過給定入射激波直接控制氣動(dòng)性能。內(nèi)部流動(dòng)時(shí)入射激波撞擊壁面會(huì)發(fā)生激波反射， 為了控制總能量損失， 發(fā)展了已知入射和反射激波的反設(shè)計(jì)方法。本文的設(shè)計(jì)方法給定了兩道入射激波的型線方程、 反射激波的激波角及其波后氣流角分布。

基于MOCC的兩級壓縮內(nèi)收縮流場反設(shè)計(jì)流場可分為“三波四區(qū)”， 如圖1所示。流場主要包含3組激波壓縮區(qū)(AO′B′,B′O′C和O′DE)和1組等熵壓縮區(qū)(CO′D), 相比于李永洲等[16]的雙彎曲入射波基準(zhǔn)流場少了一個(gè)等熵壓縮區(qū)， 所以流場長度進(jìn)一步縮短， 壓比也會(huì)顯著提升。圖1中紅色實(shí)線分別代表2道入射激波(AO和B′O′)和1道反射激波(O′D), 綠色虛線代表各部分結(jié)尾特征線， 紫色實(shí)線分別代表壓縮壁面和中心體型面， 黑色點(diǎn)劃線為中心軸(與中心體型線之間的距離定義為中心體高度)。超聲速流場的計(jì)算只是單元過程的不斷重復(fù)， 反設(shè)計(jì)同樣基于各單元過程， 下面對各區(qū)域的設(shè)計(jì)過程進(jìn)行詳細(xì)說明。

圖1 基于MOCC的兩級壓縮流場反設(shè)計(jì)—(已知入射激波和反射激波)Fig. 1 Inverse design of two-stage compression flowfield based on MOCC(known incident and reflected shocks)

1.2.1 第1道入射激波波后依賴域ABO

內(nèi)流中激波波后依賴域的求解主要涉及內(nèi)點(diǎn)的單元過程， 如圖2所示， 此時(shí)由1點(diǎn)發(fā)出的流線與由2點(diǎn)發(fā)出的逆右行特征線(C-)交于待求點(diǎn)3。其實(shí)內(nèi)外流唯一的區(qū)別在于特征線的選取， 本文研究內(nèi)收縮流場的設(shè)計(jì)， 故以內(nèi)流為例， 2階中心差分形式的常微分控制方程如下所示

(8)

其中， 上標(biāo)“-”代表離散點(diǎn)氣動(dòng)參數(shù)的平均值。圖2給出了MOCC平面和軸對稱流場內(nèi)點(diǎn)單元過程，紅色、 紫色和綠色曲線分別代表入射波、 流線和特征線。給定來流條件和激波形狀， 利用R-H關(guān)系式可以求解激波點(diǎn)1和2波后諸如壓力和氣流角等氣動(dòng)參數(shù)?；谇蟮玫臍鈩?dòng)參數(shù)和激波曲率， 結(jié)合彎曲激波方程組， 可以獲得點(diǎn)1和2波后氣動(dòng)參數(shù)的導(dǎo)數(shù)。當(dāng)點(diǎn)1和2的物理特征確定之后， 待求點(diǎn)3的位置可根據(jù)如下幾何方程求解

其中, 下標(biāo)p1,p2,p3代表單元點(diǎn)1， 2和3。一旦單元點(diǎn)3的位置確定， 相關(guān)氣動(dòng)參數(shù)初值便可通過點(diǎn)1的Taylor展開求解

圖2 內(nèi)點(diǎn)單元過程Fig. 2 Interior unit process of MOCC

最后利用點(diǎn)1， 2， 3的氣動(dòng)參數(shù)求解方程組(8)， 更新點(diǎn)3的氣動(dòng)參數(shù)。如果點(diǎn)3的新舊氣動(dòng)參數(shù)差值大于設(shè)定殘差， 就根據(jù)離散點(diǎn)之間的平均值迭代求解點(diǎn)3的位置及氣動(dòng)參數(shù)， 直到誤差小于設(shè)定殘差。經(jīng)過幾輪迭代， 點(diǎn)3的氣動(dòng)參數(shù)及其導(dǎo)數(shù)將逐漸收斂至真實(shí)值附近。

在給定來流參數(shù)和激波型線AO后， 通過斜激波關(guān)系式定義了一條非特征線的邊值線AO。利用MOCC即可求出區(qū)域ABO的流場參數(shù)并得到壁面型線AB。 步進(jìn)過程如圖3所示： 已知激波點(diǎn)B1， …，Bn， 由B1和B2依照上述的單元過程確定A2， 重復(fù)上述過程可以依次確定A3， …，An。過Bn和An的流線即為產(chǎn)生給定激波的內(nèi)收縮壁面。

圖3 MOCC步進(jìn)過程Fig. 3 Step process of MOCC

1.2.2 次級激波壓縮區(qū)域B′O′C

在求得第1道入射波的依賴域ABO后， 取過O′的水平線為中心體，O′即為兩道入射波的共同反射點(diǎn)， 過O′做任意彎曲激波B′O′為第2道入射波。事實(shí)上， 即使第2道激波的型線為直線， 由于其波前氣流的不均勻性， 其上各點(diǎn)的激波強(qiáng)度仍是不同的。選取時(shí)只需保證B′O′在ABO內(nèi)部且滿足激波角大于Mach角的要求， 由Mach角與激波角的關(guān)系可知， 第2道入射波必在Mach線O′F下游， 這可以為我們的選擇提供一定程度的限制。注意無論是第1道入射波的端點(diǎn)O還是第2道入射波的端點(diǎn)O′， 在選擇時(shí)都要考慮到高度(即到軸線的距離)的限制， 太低會(huì)導(dǎo)致激波強(qiáng)度過大， 產(chǎn)生Mach反射， 這是我們設(shè)計(jì)時(shí)不希望看到的。激波B′O′上的波前參數(shù)由插值求得， 顯然其形狀并不受O′O段的影響。通過斜激波關(guān)系式求出B′O′波后內(nèi)點(diǎn)參數(shù)， 進(jìn)而依照1.2.1節(jié)的方法求出激波B′O′的依賴域B′O′C及壓縮型面B′C， 注意B′C與第1段壁面的B′B部分并不重合。

1.2.3 等熵壓縮區(qū)域CO′D

同樣利用MOCC可求得區(qū)域CO′D的流場參數(shù)及壁面CD。首先將離散網(wǎng)格點(diǎn)分為內(nèi)點(diǎn)和反射激波點(diǎn)， 前者仍按照1.2.1節(jié)中內(nèi)點(diǎn)過程進(jìn)行求解， 激波點(diǎn)的求解則不相同， 相應(yīng)的2階中心差分格式微分控制方程如下所示

(9)

其中， 上標(biāo)“-”表示離散點(diǎn)氣動(dòng)參數(shù)的平均值。

圖4展示了MOCC激波點(diǎn)求解的單元過程。反射激波， 流線和特征線分別用紅色、 紫色和綠色曲線表示。盡管單元點(diǎn)1和2的氣動(dòng)參數(shù)及其導(dǎo)數(shù)為已知邊界條件， 激波點(diǎn)單元過程仍然需要一個(gè)額外的條件來完成求解。該條件需要與反射激波相關(guān)， 例如反射激波形狀或者激波波后參數(shù)等?？紤]到實(shí)際應(yīng)用， 大部分超聲速流場都要求出口氣流水平。因此， 本節(jié)以反射激波出口氣流水平作為額外邊界條件介紹單元過程。需要注意的是， 其他與反射激波相關(guān)的邊界條件也可以作為相應(yīng)已知條件。單元過程的具體細(xì)節(jié)也會(huì)因?yàn)檫吔鐥l件的不同有些許差異, 但整體流程是類似的。由于單元點(diǎn)1和2處的物理信息已知， 因此， 激波點(diǎn)單元過程首先需要確定單元點(diǎn)3的位置和曲率。通過弦截法可以獲得單元點(diǎn)3的激波角， 進(jìn)而利用如下所示的幾何方程可以確定單元點(diǎn)3的位置

確定位置后， 利用R-H關(guān)系式結(jié)合激波已知的物理信息可以得到單元點(diǎn)3處的初始?xì)鈩?dòng)參數(shù)。下一步， 利用單元點(diǎn)1， 2， 3的氣動(dòng)參數(shù)求解2階中心差分形式的控制方程組(9)， 更新單元點(diǎn)3處的氣動(dòng)參數(shù)。單元點(diǎn)3的氣動(dòng)參數(shù)將會(huì)不斷迭代更新直到殘差小于設(shè)定值。幾輪迭代之后， 相應(yīng)的氣動(dòng)參數(shù)及其導(dǎo)數(shù)將逐漸收斂至真實(shí)值附近。

圖4 激波點(diǎn)單元過程Fig. 4 Shock unit process of MOCC

1.2.4 反射激波波后依賴域O′DE

確定反射激波O′D后， 利用斜激波關(guān)系式求出波后流場參數(shù)， 區(qū)域O′DE以及中心體型線O′E同樣只需采用類似于1.2.1節(jié)中的MOCC內(nèi)點(diǎn)單元過程進(jìn)行求解， 只不過此時(shí)選取的是逆左行特征線(C+)， 待求點(diǎn)3由點(diǎn)1發(fā)出的流線與點(diǎn)2發(fā)出的逆左行特征線相交， 如圖5所示。

圖5 內(nèi)點(diǎn)求解單元過程(外流)Fig. 5 Interior unit process of MOCC(outflow)

MOCC通過構(gòu)建高階控制方程， 搭配彎曲激波方程組， 可用于高精度求解平面和軸對稱超聲速流場。同時(shí)， 本文的兩級壓縮設(shè)計(jì)使得在第2道激波和反射激波前的來流都是非均勻的， 這也證明了MOCC的普適性。注意到基于MOCC得到的兩級壓縮內(nèi)收縮流場與最左側(cè)來流的均勻性無關(guān)， 但設(shè)計(jì)流場求解時(shí)有可能出現(xiàn)網(wǎng)格線相交， 也即同簇特征線相交的情況， 這說明給定的激波和流場參數(shù)分布不存在對應(yīng)的物理解， 需要調(diào)整輸入重新設(shè)計(jì)。

綜上所述， 本文基于MOCC的兩級壓縮內(nèi)收縮流場設(shè)計(jì)方法不僅能夠確保兩道激波符合給定型線方程， 實(shí)現(xiàn)流場波系結(jié)構(gòu)的主動(dòng)控制； 還能保證反射激波波后流向角的水平分布， 便于參數(shù)均勻的反設(shè)計(jì)； 與通過給定壁面沿程參數(shù)分布或直接給定等熵壓縮段壁面型線來求解流場的方法存在很大不同。由于進(jìn)氣道的流量特性取決于入射激波的形狀和位置， 因此只要基準(zhǔn)流場的入射激波完全將進(jìn)口封閉， 就可以實(shí)現(xiàn)全流量捕獲。反射激波的強(qiáng)度和位置僅影響下游的氣動(dòng)特征， 并不會(huì)改變流場的內(nèi)乘波特性。

2 設(shè)計(jì)方法數(shù)值校驗(yàn)暨超聲速流場反設(shè)計(jì)

隨著高超聲速飛行器的工作范圍不斷擴(kuò)大， 為了兼顧寬來流Mach數(shù)下的總體性能(尤其是巡航點(diǎn))， 其型面設(shè)計(jì)Mach數(shù)將不斷增加[16]。對于內(nèi)收縮進(jìn)氣道， 此時(shí)氣流僅通過兩道波主壓縮(等熵壓縮區(qū)域的壓縮能力并不強(qiáng))， 造成的激波損失將顯著增大。而且其出口壓比受進(jìn)氣道長度的限制并不會(huì)很高。在二維和軸對稱進(jìn)氣道設(shè)計(jì)中， 兩級壓縮內(nèi)收縮流場便能較好地解決上述問題。

2.1 超聲速流場反設(shè)計(jì)

基于彎曲激波理論的MOCC可用于高精度求解平面和軸對稱超聲速流場， 來流經(jīng)過壁面偏轉(zhuǎn)時(shí)會(huì)在前緣產(chǎn)生入射激波， 隨后入射波撞擊中心體產(chǎn)生反射激波， 反設(shè)計(jì)分為入射波和反射波兩部分。首先給定兩道入射波的激波型線方程， 利用R-H方程和彎曲激波方程根據(jù)來流條件可以求解波后氣動(dòng)參數(shù)及梯度， 本文中來流條件均設(shè)置為 26 km 海拔理想大氣， Mach數(shù)6。對于每段離散激波， 通過內(nèi)點(diǎn)過程不斷迭代獲得收斂的氣動(dòng)參數(shù)及梯度， 然后沿流向推進(jìn)對第1道入射波部分進(jìn)行求解， 被第2道入射波截取后的流場特征如圖6上半部分所示。在得到第2道入射波部分后， 結(jié)尾特征線的參數(shù)將作為求解反射激波部分的邊界條件。由于中心體高度不同， 反射波幾何方程也不相同， 這里通過給定反射激波角變化來定義反射波的形狀， 而反射波后的流動(dòng)方向保持水平。使用激波過程可求解反射波附近的氣動(dòng)參數(shù)， 進(jìn)而利用內(nèi)點(diǎn)過程計(jì)算壁面、 反射波和特征線包圍的區(qū)域， 計(jì)算過程終止于壁面和反射波的交點(diǎn)處。由于入射波的不確定性， 壁面長度無法在求解之前限定而只能在求解結(jié)束后獲得。Mach桿后的亞聲速區(qū)域會(huì)降低流場的氣動(dòng)性能， 所以內(nèi)收縮流動(dòng)中通常設(shè)計(jì)中心體來避免Mach反射， 中心體高度的不同也會(huì)導(dǎo)致求解域范圍產(chǎn)生差異。最后借助反射波附近參數(shù)通過內(nèi)點(diǎn)過程求解反射波波后流場， 圖6下半部分展示了完整的內(nèi)收縮流場流動(dòng)特征。

由于流動(dòng)形式的不同， 即使在相同的入射波下平面和軸對稱反設(shè)計(jì)流場也具有不同的流動(dòng)特征。梯度是造成這種現(xiàn)象的主要原因， 相關(guān)分析見文獻(xiàn)[33]， 這里不再贅述。

(a) Planar internal flowfields

(b) Axially symmetric internal flowfields圖6 均勻來流已知入射和反射激波超聲速流場反設(shè)計(jì)Fig. 6 Inverse design of planar and axially symmetric internal flowfields with shocks in uniform flows

本節(jié)將分別針對第1道入射波彎曲、 第2道入射波平直和雙彎曲入射激波的情況， 以均勻來流已知激波分布的平面和軸對稱內(nèi)收縮流動(dòng)為例， 根據(jù)1.2 節(jié)的設(shè)計(jì)方法， 將MOCC求解的流場與CFD結(jié)果進(jìn)行比較。進(jìn)行數(shù)值驗(yàn)證的同時(shí)完成超聲速流場的反設(shè)計(jì)。

CFD在正設(shè)計(jì)中具有廣泛應(yīng)用， 已被證實(shí)是一種精確的方法， 本文利用CFD模擬了給定激波分布的內(nèi)收縮流動(dòng)， 以驗(yàn)證設(shè)計(jì)方法的精度。由于基本流場的求解中忽略了黏性， 因此CFD計(jì)算中選擇無黏模型， 無黏通量利用Roe-FDS算法計(jì)算， 控制方程采用完全隱式格式的2階迎風(fēng)離散法求解。結(jié)構(gòu)化網(wǎng)格及邊界條件如圖7所示， 為了清晰起見， 實(shí)際網(wǎng)格密度已降低。對近壁面網(wǎng)格采用局部加密， 并在出口水平延伸一段距離用來觀察結(jié)尾的反射激波串。經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證后， 全計(jì)算域?qū)嶋H網(wǎng)格數(shù)量由各自的尺寸決定， 分辨率足以獲得收斂結(jié)果。入口和出口分別設(shè)置為壓力遠(yuǎn)場和壓力出口邊界條件， 來流設(shè)置為理想氣體， 壁面采用滑移邊界條件。收斂標(biāo)準(zhǔn)是監(jiān)控變量殘差下降3個(gè)量級且進(jìn)出口流量守恒。

圖7 結(jié)構(gòu)網(wǎng)格和邊界條件Fig. 7 Grid structure and boundary conditions

2.2 彎直入射兩級壓縮基準(zhǔn)流場

2.2.1 平面流動(dòng)

本算例中第1道入射波的型線方程為y=-0.1x2-0.2x+1， 第2道入射激波的型線方程為y=-0.892 1x+1.865 3， 反射激波角的變化率為dθre/dx=-0.2， 中心體高度為0.404 1， 質(zhì)量平均的出口壓比和總壓恢復(fù)系數(shù)分別為 128.48 和 0.350 7。

圖8顯示了由激波分布計(jì)算得到的平面內(nèi)收縮流場流動(dòng)特性， 黑色實(shí)線代表CFD結(jié)果的等參數(shù)線和激波， MOCC的等參數(shù)線和激波則分別用紅色虛線和白色虛線表示。使用MOCC獲得的平面流場與CFD結(jié)果具有十分相似的流動(dòng)特征， 兩道入射激波按照預(yù)設(shè)匯聚于(1.638， 0.404)， 證明MOCC具有足夠高的精度求解已知參數(shù)分布邊界的平面超聲速流場。

圖8 彎直入射波平面內(nèi)收縮流場Fig. 8 Planar internal flowfield with one direct incident shock and one curved shock

分別提取壁面和中心體處的壓比參數(shù)進(jìn)行對比， 如圖9(a)所示， 其中散點(diǎn)和線條分別代表CFD和MOCC的結(jié)果， 壁面和中心體處壓力分布則用黑色和紅色進(jìn)行區(qū)分。可以發(fā)現(xiàn)無論是數(shù)值還是壓升規(guī)律， 無論是壁面還是中心體， MOCC的結(jié)果均和CFD吻合得很好。此外， 如圖9(b)所示， 分別提取了CFD結(jié)果流場在反射波后和喉道處的流向角分布， 紅色實(shí)線代表喉道處的流向角分布， 反射波后的流向角分布則用黑色虛線表示?？梢园l(fā)現(xiàn)反射波后的流向角與水平偏差不超過0.015， 即使繼續(xù)延伸至喉道， 這個(gè)差值也不超過 0.04， 很好地保持了給定的流向角分布規(guī)律。

2.2.2 軸對稱流動(dòng)

本算例中第1道入射波的型線方程為y=-0.03x2-0.2x+1， 第2道入射激波的型線方程為y=-0.510 0x+1.427 2， 反射激波角的變化率為dθre/dx=0.2， 中心體高度為0.591 9， 質(zhì)量平均的出口壓比和總壓恢復(fù)系數(shù)分別為 50.36 和 0.628 7。

(a) Pressure distribution along the wall and centerbody

(b) Deflection angle distribution along the throat and just behind the reflected shock圖9 平面內(nèi)收縮流場設(shè)計(jì)參數(shù)對比Fig. 9 Comparison of design parameters of planar internal flowfield

圖10顯示了由激波分布計(jì)算得到的軸對稱內(nèi)收縮流場流動(dòng)特性， 相關(guān)線條定義與平面情況一致。使用MOCC獲得的軸對稱流場與CFD結(jié)果具有十分相似的流動(dòng)特征， 兩道入射激波按照預(yù)設(shè)匯聚于(1.638， 0.592)， 證明MOCC具有足夠高的精度求解已知參數(shù)分布邊界的軸對稱超聲速流場。

圖10 彎直入射波軸對稱內(nèi)收縮流場Fig. 10 Axisymmetric internal flowfield with one direct incident shock and one curved shock

分別提取壁面和中心體處的壓比參數(shù)進(jìn)行對比， 如圖11(a)所示， 相關(guān)線條定義與平面情況一致。MOCC的數(shù)值和壓升規(guī)律結(jié)果在壁面和中心體上的分布均與CFD一致。分別提取了CFD結(jié)果流場在反射波后和喉道處的流向角分布情況， 如圖11(b)， 發(fā)現(xiàn)反射波后的流向角與水平偏差不超過0.006， 即使繼續(xù)延伸至喉道， 這個(gè)差值也不超過 0.035， 流向角保持設(shè)計(jì)的水平分布。

(a) Pressure distribution along the wall and centerbody

(b) Deflection angle distribution along the throat and just behind the reflected shock圖11 軸對稱內(nèi)收縮流場設(shè)計(jì)參數(shù)對比Fig. 11 Comparison of design parameters of axisymmetric internal flowfield

此外， 觀察CFD結(jié)果中延伸段內(nèi)的流場， 如圖12所示， 可以發(fā)現(xiàn)反射激波準(zhǔn)確打在肩點(diǎn)(1.862， 0.642)處， 同時(shí)水平段內(nèi)存在不斷反射的激波串， 但由于激波強(qiáng)度較弱， 所以暫未做消波處理。

圖12 水平延伸段反射激波串Fig. 12 Reflected shock train in the horizontal extension

2.3 雙彎曲入射兩級壓縮基準(zhǔn)流場

雙彎曲入射兩級壓縮基準(zhǔn)流場是最普適的情況， 本文提出的設(shè)計(jì)方法可以用于對任意波系結(jié)構(gòu)的兩級壓縮流場的反設(shè)計(jì)。本節(jié)對理論求解的雙彎曲入射激波流場進(jìn)行數(shù)值驗(yàn)證， 分別以平面和軸對稱流動(dòng)為例。平面算例中第1道入射波的型線方程為y=-0.03x2-0.2x+1， 第2道入射激波的型線方程為y=-0.023 1x2-0.376 3x+1.358 1， 反射激波角的變化率為dθre/dx=-0.15， 中心體高度為 0.406 3， 質(zhì)量平均的出口壓比和總壓恢復(fù)系數(shù)分別為13.37和0.747 2。軸對稱算例中第1道入射波的型線方程為y=-0.03x2-0.2x+1， 第2道入射激波的型線方程為y=-0.060 1x2-0.225 1x+1.204 9， 反射激波角的變化率為 dθre/dx=0.2， 中心體高度為0.406 3， 質(zhì)量平均的出口壓比和總壓恢復(fù)系數(shù)分別為53.10和0.695 8。

由圖13可以看出， 此時(shí)兩道入射激波都是曲線， 盡管第2道激波的彎曲程度并不明顯， 這是因?yàn)橐獫M足激波角大于Mach角的要求， 因此若激波彎曲程度過大， 會(huì)導(dǎo)致在靠近壁面處激波角過小。Mach數(shù)云圖和等值線吻合良好， CFD結(jié)果的3道激波及反射點(diǎn)、 喉道點(diǎn)、 壁面轉(zhuǎn)折點(diǎn)等對應(yīng)特征點(diǎn)均與理論值十分接近。以反射點(diǎn)為例， 由于本算例中平面和軸對稱流場的兩道入射激波在設(shè)計(jì)時(shí)具有相同的交點(diǎn)， 故圖13(a)中平面流場和圖13(b)中軸對稱流場的反射點(diǎn)均為(2.226， 0.406)， 與設(shè)計(jì)一致。綜合上述數(shù)值驗(yàn)證的結(jié)果， 可以證明本文的反設(shè)計(jì)方法對任意滿足物理解的激波型線組合都是適用的。

(a) Planar

(b) Axisymmetric圖13 雙彎曲入射波內(nèi)收縮流場Fig. 13 Internal flowfield with double curved incident shocks

3 兩級壓縮基準(zhǔn)流場性能研究

本章利用上述程序設(shè)計(jì)了一系列軸對稱兩級壓縮基準(zhǔn)流場， 并研究了第2道入射激波激波角分布(以彎曲系數(shù)為代表)和反射激波激波角分布(以激波角沿流向變化規(guī)律為代表)對基準(zhǔn)流場的氣動(dòng)和幾何性能的影響。在對比基準(zhǔn)流場的性能時(shí)， 以總收縮比、 內(nèi)收縮比、 質(zhì)量平均的喉道總壓恢復(fù)系數(shù)和出口壓比為主要標(biāo)準(zhǔn)。

3.1 第2道入射激波激波角分布的影響

本節(jié)根據(jù)先前研究結(jié)果選取第1道彎曲激波型線方程為y=-0.03x2-0.2x+1， 反射激波激波角系數(shù)取為 dθre/dx=-0.379， 隨后分別設(shè)計(jì)了一系列不同第2道激波激波角分布的基準(zhǔn)流場用來做對比試驗(yàn)。在選擇算例時(shí)為控制變量唯一， 本節(jié)固定第2道入射波的起始點(diǎn)(圖1中B′)和反射點(diǎn)(圖1中O′)位置， 這里選取B′位置為(1.081 9， 0.875 5)，O′位置為(1.637 9， 0.591 9)。通過改變第2道入射波彎曲程度來探究激波角分布對流場氣動(dòng)性能的影響， 這里彎曲程度用系數(shù)coc(coefficient of curve)表示， 定義如下

其中，xmid為彎曲激波B′O′中點(diǎn)的橫坐標(biāo)， 由定義可以看出， 隨著彎曲程度的增加， 激波角在靠近中心體的部分增加得更快， 而在靠近壁面附近則變化平緩。由于受到Mach角的限制， 顯然coc的選取不是任意的， 對于給定的來流，coc存在對應(yīng)的上下限。針對本組算例，coc的取值范圍為0～0.038。如圖14(a)所示， 黑色虛線和紅色實(shí)線分別代表流場增壓比和總壓恢復(fù)系數(shù)隨coc的變化規(guī)律。流場幾何參數(shù)隨coc的變化規(guī)律見圖14(b)， 流場總收縮比和內(nèi)收縮比分別用黑色虛線和紅色實(shí)線表示。

由圖14(a)可以看出， 隨著第2道入射波彎曲程度的增大， 增壓比呈現(xiàn)單調(diào)上升趨勢， 而總壓恢復(fù)系數(shù)則隨之下降， 且兩者隨coc的變化接近線性。這是由于隨第2道入射波愈發(fā)彎曲， 激波強(qiáng)度隨之增大， 對氣流的壓縮增強(qiáng)， 壓比增加； 與此同時(shí)激波壓縮區(qū)占比增加， 流場的總壓損失增大， 總壓恢復(fù)系數(shù)減小。另一方面， 如圖14(b)所示， 隨著coc的增加， 無論是總收縮比還是內(nèi)收縮比都顯著增加， 這是因?yàn)殡S著第2道入射波在靠近中心體處的強(qiáng)度增加， 氣流在經(jīng)過第2道入射波時(shí)的偏轉(zhuǎn)角隨之增大， 第2段壁面轉(zhuǎn)折增加， 導(dǎo)致等熵壓縮區(qū)減小， 反射激波所打肩點(diǎn)隨之降低， 所以盡管入口型面保持不變， 總收縮比和內(nèi)收縮比依然不斷升高。需要注意的是， 總收縮比隨coc增加得越來越快， 而內(nèi)收縮比則越來越慢。

(a) Boost ratio and total pressure recovery coefficient

(b) Total shrinkage ratio and internal shrinkage ratio圖14 流場性能隨第2道入射激波激波角分布的變化趨勢Fig. 14 Trend of flowfield performance with the shock angle distribution of the second incident shock

3.2 反射激波激波角分布

本節(jié)根據(jù)先前研究結(jié)果選取第1道彎曲激波型線方程為y=-0.03x2-0.2x+1， 第2道彎曲激波型線方程為y=-0.060 1x2-0.225 1x+1.204 9(相應(yīng)的coc為0.04)， 隨后分別設(shè)計(jì)了一系列不同反射激波激波角分布的基準(zhǔn)流場用來做對比試驗(yàn)。反射激波激波角分布用系數(shù)cor(coefficient of reflect)表示， 定義如下

其中，θre為反射激波O′D的激波角， 由定義可以看出， 隨著cor的增加， 激波角沿x方向增加得更快。同樣地， 由于受到物理解的限制，cor的選取不是任意的， 對于給定的來流條件，cor存在對應(yīng)的上下限。針對本組算例，cor的取值范圍為 -0.165～3.82， 相應(yīng)結(jié)果見圖15， 相關(guān)線條定義與圖14一致。

(a) Boost ratio and total pressure recovery coefficient

(b) Total shrinkage ratio and internal shrinkage ratio圖15 流場性能隨反射激波激波角分布的變化趨勢Fig. 15 Trend of flowfield performance with the shock angle distribution of the reflected shock

由圖15(a)可以看出， 隨著反射激波的增強(qiáng)， 增壓比隨之增加， 總壓恢復(fù)系數(shù)隨之減小且變化速度有所下降， 基準(zhǔn)流場壓縮效率略微降低。這是由于兩道入射波保持不變， 因此前兩組激波壓縮區(qū)保持不變， 隨著反射激波激波角分布系數(shù)的增加， 反射激波在靠近肩點(diǎn)時(shí)的激波強(qiáng)度大幅增加。肩點(diǎn)位置隨之前移， 使得等熵壓縮區(qū)略微減小， 因此總壓比和總壓恢復(fù)系數(shù)的變化越來越慢。另一方面， 如圖15(b)所示， 隨著cor的增加， 總收縮比隨之增加， 但到極限值時(shí)基本不變， 而內(nèi)收縮比則先增加后減小， 這也是由于反射激波激波角分布的變化影響了等熵壓縮區(qū)和反射激波波后依賴域的分配比例， 從而導(dǎo)致肩點(diǎn)位置和喉道高度的變化。

3.3 不同設(shè)計(jì)參數(shù)對彼此取值范圍的影響

通過上述研究也可以看出， 在給定的設(shè)計(jì)參數(shù)下， 無論是coc還是cor都存在各自的取值范圍， 超出這個(gè)范圍則不存在對應(yīng)的物理解。顯然設(shè)計(jì)參數(shù)A的變化對設(shè)計(jì)參數(shù)B的取值范圍是會(huì)有影響的， 進(jìn)而會(huì)影響到流場的氣動(dòng)和幾何性能， 這顯然是流場優(yōu)化設(shè)計(jì)需要考慮的， 本節(jié)就以coc和cor的變化及其對彼此取值范圍的影響進(jìn)行分析討論。

首先保證3.1節(jié)中的其余設(shè)計(jì)參數(shù)不變， 分別取不同的cor， 如表1所示， 其中case 5即為3.1節(jié)中的算例。注意到在cor>-0.379時(shí)，coc的取值可以為負(fù)， 也即在凹向來流的第1道激波后， 第2道激波凸向來流。這雖然和平時(shí)的印象相悖， 但要注意到此時(shí)第2道激波波前來流非均勻， 因此不能僅根據(jù)第2道激波幾何形狀的凹凸來判斷激波的強(qiáng)度變化， 而要結(jié)合激波角的變化。觀察到隨著cor的減小，coc的取值下限不斷增加， 而上限則不變。我們知道激波角應(yīng)大于波前Mach數(shù)對應(yīng)的Mach角， 而若coc較大則會(huì)使得在靠近壓縮面處第2道激波的激波角過小。因而在確定的第1道激波后條件下， 第2道激波coc的上限不變， 也即coc的上限主要受制于來流條件。但coc的下限則并不僅僅依賴于來流條件， 由于等熵區(qū)CO′D的存在而同時(shí)受到反射激波角分布的影響。隨著反射激波沿流向強(qiáng)度遞減速度逐漸提升， 靠近壓縮面處的反射激波強(qiáng)度逐漸減弱， 順流管向前追蹤至第2道激波， 第2道激波靠近壁面處強(qiáng)度必須隨之減小， 因而要求第2道激波的彎曲程度增加， 即coc的下限提升。

表1 不同cor下coc的取值范圍

保證3.2節(jié)中的其余設(shè)計(jì)參數(shù)不變， 分別取不同的coc， 如表2所示， 其中工況2即為3.2中算例。受制于來流條件和出口不能跨聲速的要求，coc的取值僅限于0.039～0.043， 觀察到隨著coc的減小，cor的取值上下限均有略微增加。與上述原因一致， 隨著第2道激波彎曲程度的增加， 反射激波沿流向強(qiáng)度遞減速度逐漸提升， 即cor上下限隨之遞減。

表2 不同coc下cor的取值范圍

由本節(jié)的研究可以看出， 流場的氣動(dòng)參數(shù)隨設(shè)計(jì)參數(shù)的改變而變化， 不同設(shè)計(jì)參數(shù)的取值范圍也會(huì)受到彼此的影響， 進(jìn)而影響流場性能。在無法直接由設(shè)計(jì)要求構(gòu)造流場時(shí)， 基于MOCC給定激波的反設(shè)計(jì)在不需要進(jìn)行復(fù)雜而耗時(shí)的CFD驗(yàn)證前就可以得到流場性能的預(yù)先估計(jì)， 從而在選擇設(shè)計(jì)參數(shù)時(shí)不再盲目， 而是有的放矢。在給定進(jìn)氣道基準(zhǔn)流場的氣動(dòng)和幾何參數(shù)的需求后， 可以改變設(shè)計(jì)參數(shù)進(jìn)行一系列的流場設(shè)計(jì)， 得到參數(shù)的變化規(guī)律， 在結(jié)果區(qū)間快速選取設(shè)計(jì)參數(shù)， 完成流場的優(yōu)化設(shè)計(jì)。

4 結(jié)論

本文基于MOCC提出了一種兩級壓縮內(nèi)收縮基準(zhǔn)流場的結(jié)構(gòu)與設(shè)計(jì)方法， 并通過數(shù)值模擬檢驗(yàn)了設(shè)計(jì)方法對于任意彎曲激波波系結(jié)構(gòu)的平面和軸對稱流場的準(zhǔn)確性， 主要得出以下結(jié)論：

1)基于MOCC的設(shè)計(jì)方法靈活且高效， 在具有高精度的同時(shí)可節(jié)省計(jì)算資源， 適用于任何波系組成的平面和軸對稱流場設(shè)計(jì)。

2)給定入反射激波參數(shù)和出口流向角分布反設(shè)計(jì)流場壓縮型面可以控制基準(zhǔn)流場的波系結(jié)構(gòu)和壓縮效率， 保證雙彎曲入射激波封口且出口氣流均勻。

3)分別改變第2道入射激波和反射激波的激波角分布， 求解了一系列具有中心體的軸對稱流場， 對比發(fā)現(xiàn)設(shè)計(jì)參數(shù)改變對流場氣動(dòng)和幾何性能的影響規(guī)律， 以及不同設(shè)計(jì)參數(shù)的變化對彼此取值范圍的影響規(guī)律， 便于設(shè)計(jì)前對流場參數(shù)的選取。

4)高求解精度和快設(shè)計(jì)速度使得MOCC成為超聲速流場分析和反設(shè)計(jì)領(lǐng)域的良好候選者； 多波入射的結(jié)構(gòu)使得流場長度縮短、 壓縮效率提升， 對于拓寬基準(zhǔn)流場的設(shè)計(jì)思路具有重要意義。未來將進(jìn)一步開展流場設(shè)計(jì)的優(yōu)化選取和基于此的內(nèi)乘波進(jìn)氣道的相關(guān)研究。

致謝感謝國家自然科學(xué)基金(U20A3069, U21B6003)的資助, 也感謝1912項(xiàng)目組對本文工作的支持。