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

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

引 言

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

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

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

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

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

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

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

壓力梯度

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

(1)

流線曲率

D=?δ/?s

(2)

沿流線的動量方程

(3)

沿流線的能量方程

(4)

沿特征線方向的動量方程

=ρV2(Pcosμ?Dsinμ)

(5)

沿特征線方向的流量方程

(6)

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

(7)

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

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

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

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

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

圖1 基于MOCC的兩級壓縮流場反設(shè)計—(已知入射激波和反射激波)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)點的單元過程， 如圖2所示， 此時由1點發(fā)出的流線與由2點發(fā)出的逆右行特征線(C-)交于待求點3。其實內(nèi)外流唯一的區(qū)別在于特征線的選取， 本文研究內(nèi)收縮流場的設(shè)計， 故以內(nèi)流為例， 2階中心差分形式的常微分控制方程如下所示

(8)

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

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

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

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

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

圖3 MOCC步進過程Fig. 3 Step process of MOCC

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

在求得第1道入射波的依賴域ABO后， 取過O′的水平線為中心體，O′即為兩道入射波的共同反射點， 過O′做任意彎曲激波B′O′為第2道入射波。事實上， 即使第2道激波的型線為直線， 由于其波前氣流的不均勻性， 其上各點的激波強度仍是不同的。選取時只需保證B′O′在ABO內(nèi)部且滿足激波角大于Mach角的要求， 由Mach角與激波角的關(guān)系可知， 第2道入射波必在Mach線O′F下游， 這可以為我們的選擇提供一定程度的限制。注意無論是第1道入射波的端點O還是第2道入射波的端點O′， 在選擇時都要考慮到高度(即到軸線的距離)的限制， 太低會導(dǎo)致激波強度過大， 產(chǎn)生Mach反射， 這是我們設(shè)計時不希望看到的。激波B′O′上的波前參數(shù)由插值求得， 顯然其形狀并不受O′O段的影響。通過斜激波關(guān)系式求出B′O′波后內(nèi)點參數(shù)， 進而依照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)格點分為內(nèi)點和反射激波點， 前者仍按照1.2.1節(jié)中內(nèi)點過程進行求解， 激波點的求解則不相同， 相應(yīng)的2階中心差分格式微分控制方程如下所示

(9)

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

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

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

圖4 激波點單元過程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)點單元過程進行求解， 只不過此時選取的是逆左行特征線(C+)， 待求點3由點1發(fā)出的流線與點2發(fā)出的逆左行特征線相交， 如圖5所示。

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

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

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

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

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

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

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

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

(a) Planar internal flowfields

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

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

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

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

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

2.2.1 平面流動

本算例中第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顯示了由激波分布計算得到的平面內(nèi)收縮流場流動特性， 黑色實線代表CFD結(jié)果的等參數(shù)線和激波， MOCC的等參數(shù)線和激波則分別用紅色虛線和白色虛線表示。使用MOCC獲得的平面流場與CFD結(jié)果具有十分相似的流動特征， 兩道入射激波按照預(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ù)進行對比， 如圖9(a)所示， 其中散點和線條分別代表CFD和MOCC的結(jié)果， 壁面和中心體處壓力分布則用黑色和紅色進行區(qū)分?？梢园l(fā)現(xiàn)無論是數(shù)值還是壓升規(guī)律， 無論是壁面還是中心體， MOCC的結(jié)果均和CFD吻合得很好。此外， 如圖9(b)所示， 分別提取了CFD結(jié)果流場在反射波后和喉道處的流向角分布， 紅色實線代表喉道處的流向角分布， 反射波后的流向角分布則用黑色虛線表示。可以發(fā)現(xiàn)反射波后的流向角與水平偏差不超過0.015， 即使繼續(xù)延伸至喉道， 這個差值也不超過 0.04， 很好地保持了給定的流向角分布規(guī)律。

2.2.2 軸對稱流動

本算例中第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è)計參數(shù)對比Fig. 9 Comparison of design parameters of planar internal flowfield

圖10顯示了由激波分布計算得到的軸對稱內(nèi)收縮流場流動特性， 相關(guān)線條定義與平面情況一致。使用MOCC獲得的軸對稱流場與CFD結(jié)果具有十分相似的流動特征， 兩道入射激波按照預(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ù)進行對比， 如圖11(a)所示， 相關(guān)線條定義與平面情況一致。MOCC的數(shù)值和壓升規(guī)律結(jié)果在壁面和中心體上的分布均與CFD一致。分別提取了CFD結(jié)果流場在反射波后和喉道處的流向角分布情況， 如圖11(b)， 發(fā)現(xiàn)反射波后的流向角與水平偏差不超過0.006， 即使繼續(xù)延伸至喉道， 這個差值也不超過 0.035， 流向角保持設(shè)計的水平分布。

(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è)計參數(shù)對比Fig. 11 Comparison of design parameters of axisymmetric internal flowfield

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

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

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

雙彎曲入射兩級壓縮基準(zhǔn)流場是最普適的情況， 本文提出的設(shè)計方法可以用于對任意波系結(jié)構(gòu)的兩級壓縮流場的反設(shè)計。本節(jié)對理論求解的雙彎曲入射激波流場進行數(shù)值驗證， 分別以平面和軸對稱流動為例。平面算例中第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可以看出， 此時兩道入射激波都是曲線， 盡管第2道激波的彎曲程度并不明顯， 這是因為要滿足激波角大于Mach角的要求， 因此若激波彎曲程度過大， 會導(dǎo)致在靠近壁面處激波角過小。Mach數(shù)云圖和等值線吻合良好， CFD結(jié)果的3道激波及反射點、 喉道點、 壁面轉(zhuǎn)折點等對應(yīng)特征點均與理論值十分接近。以反射點為例， 由于本算例中平面和軸對稱流場的兩道入射激波在設(shè)計時具有相同的交點， 故圖13(a)中平面流場和圖13(b)中軸對稱流場的反射點均為(2.226， 0.406)， 與設(shè)計一致。綜合上述數(shù)值驗證的結(jié)果， 可以證明本文的反設(shè)計方法對任意滿足物理解的激波型線組合都是適用的。

(a) Planar

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

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

本章利用上述程序設(shè)計了一系列軸對稱兩級壓縮基準(zhǔn)流場， 并研究了第2道入射激波激波角分布(以彎曲系數(shù)為代表)和反射激波激波角分布(以激波角沿流向變化規(guī)律為代表)對基準(zhǔn)流場的氣動和幾何性能的影響。在對比基準(zhǔn)流場的性能時， 以總收縮比、 內(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è)計了一系列不同第2道激波激波角分布的基準(zhǔn)流場用來做對比試驗。在選擇算例時為控制變量唯一， 本節(jié)固定第2道入射波的起始點(圖1中B′)和反射點(圖1中O′)位置， 這里選取B′位置為(1.081 9， 0.875 5)，O′位置為(1.637 9， 0.591 9)。通過改變第2道入射波彎曲程度來探究激波角分布對流場氣動性能的影響， 這里彎曲程度用系數(shù)coc(coefficient of curve)表示， 定義如下

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

由圖14(a)可以看出， 隨著第2道入射波彎曲程度的增大， 增壓比呈現(xiàn)單調(diào)上升趨勢， 而總壓恢復(fù)系數(shù)則隨之下降， 且兩者隨coc的變化接近線性。這是由于隨第2道入射波愈發(fā)彎曲， 激波強度隨之增大， 對氣流的壓縮增強， 壓比增加； 與此同時激波壓縮區(qū)占比增加， 流場的總壓損失增大， 總壓恢復(fù)系數(shù)減小。另一方面， 如圖14(b)所示， 隨著coc的增加， 無論是總收縮比還是內(nèi)收縮比都顯著增加， 這是因為隨著第2道入射波在靠近中心體處的強度增加， 氣流在經(jīng)過第2道入射波時的偏轉(zhuǎn)角隨之增大， 第2段壁面轉(zhuǎn)折增加， 導(dǎo)致等熵壓縮區(qū)減小， 反射激波所打肩點隨之降低， 所以盡管入口型面保持不變， 總收縮比和內(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è)計了一系列不同反射激波激波角分布的基準(zhǔ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)可以看出， 隨著反射激波的增強， 增壓比隨之增加， 總壓恢復(fù)系數(shù)隨之減小且變化速度有所下降， 基準(zhǔn)流場壓縮效率略微降低。這是由于兩道入射波保持不變， 因此前兩組激波壓縮區(qū)保持不變， 隨著反射激波激波角分布系數(shù)的增加， 反射激波在靠近肩點時的激波強度大幅增加。肩點位置隨之前移， 使得等熵壓縮區(qū)略微減小， 因此總壓比和總壓恢復(fù)系數(shù)的變化越來越慢。另一方面， 如圖15(b)所示， 隨著cor的增加， 總收縮比隨之增加， 但到極限值時基本不變， 而內(nèi)收縮比則先增加后減小， 這也是由于反射激波激波角分布的變化影響了等熵壓縮區(qū)和反射激波波后依賴域的分配比例， 從而導(dǎo)致肩點位置和喉道高度的變化。

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

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