周航，金志光

南京航空航天大學(xué) 能源與動(dòng)力學(xué)院, 南京 210016

在過(guò)去的幾十年間，高超聲速技術(shù)取得了極大進(jìn)步，并不斷推動(dòng)著高超聲速飛行器研究研制工作的變革與發(fā)展[1-2]，高超聲速飛行器前體、進(jìn)氣道及二者的氣動(dòng)一體化設(shè)計(jì)是其中的重要研究方向[3-6]。

乘波體因其顯著的高升阻比和高效預(yù)壓縮特性，長(zhǎng)期以來(lái)被公認(rèn)是高超聲速飛行器設(shè)計(jì)領(lǐng)域的一種理想構(gòu)型[7-8]。乘波理論中，前緣激波與飛行器邊緣在精心設(shè)計(jì)下能夠完美貼合，飛行器猶如“乘坐”在激波面之上，故而得名。在利用乘波概念進(jìn)行高超聲速飛行器前體、進(jìn)氣道設(shè)計(jì)時(shí)，一般首先需要進(jìn)行無(wú)黏基準(zhǔn)流場(chǎng)構(gòu)建，其中關(guān)鍵在于對(duì)前緣激波的設(shè)計(jì)和控制，然后從激波面上選定的前緣線進(jìn)行流線追蹤[9]，即可得到具備乘波特性的氣動(dòng)壓縮面。在乘波設(shè)計(jì)理念的指導(dǎo)下，很自然地逐漸發(fā)展形成了兩種設(shè)計(jì)思路：① 給定激波生成體，通過(guò)數(shù)值計(jì)算得到激波面和波后流場(chǎng)；② 給定激波形狀，通過(guò)求解反問(wèn)題得到壓縮面和波后流場(chǎng)。

第1類(lèi)設(shè)計(jì)思路中，研究者首先基于楔、錐形流的理論解展開(kāi)設(shè)計(jì)，例如分別由Nonweiler[10]和Jones等[11]提出的“Λ”構(gòu)型乘波體和錐導(dǎo)乘波體。Rasmussen[12]應(yīng)用高超聲速小擾動(dòng)理論基于小攻角圓錐繞流和橢圓錐繞流設(shè)計(jì)了乘波體構(gòu)型。隨著CFD技術(shù)的發(fā)展，更多更加復(fù)雜的流場(chǎng)被用作乘波體基準(zhǔn)流場(chǎng)，通過(guò)直接求解三維歐拉方程，Takashima和Lewis[13]研究了楔錐混合乘波體，Cui等[14]基于方錐、花錐、十字錐繞流研究了廣義錐導(dǎo)乘波體。

第2類(lèi)設(shè)計(jì)思路起始于Sobieczky等[15-17]提出的密切軸對(duì)稱(chēng)(Osculating Axisymemtric flow，OA)理論及橫向推進(jìn)求解(Cross-Stream Marching，CSM)方法，其能夠基于廣義激波形狀進(jìn)行乘波體構(gòu)型設(shè)計(jì)。自該理論提出后，直到今天仍引領(lǐng)著密切錐乘波體的研究熱潮，有關(guān)該理論的修正[18-19]、發(fā)展[20-23]和應(yīng)用[24-28]工作層出不窮。其中包括，Chauffour和Lewis[18]與Rodi等[19]先后針對(duì)密切錐乘波體設(shè)計(jì)中的橫向壓力梯度與橫向流動(dòng)進(jìn)行的修正；Qiao等[20]發(fā)展的CSM方法，其給出了一種精確性和魯棒性得到提升的CSMP方法；Zheng等[21]拓展并提出的多密切錐方法(Mutiple-OCs)，通過(guò)在密切面內(nèi)的非共軸特征線(Noncoaxical Method of Characteristics)計(jì)算，使得沿激波面周向不同位置可以有不同的激波強(qiáng)度；Chen等[22]通過(guò)調(diào)整橫截面激波曲率中心位置，給出的高容積率密切錐乘波體設(shè)計(jì)方法；Liu等[23]基于各密切面上不同的來(lái)流馬赫數(shù)設(shè)計(jì)的寬速域密切乘波體。

事實(shí)上，兩種思路均能得到滿(mǎn)足要求的基準(zhǔn)流場(chǎng)，但各有長(zhǎng)短。前者在設(shè)計(jì)時(shí)能夠得到更多復(fù)雜構(gòu)型下的超聲速基準(zhǔn)流場(chǎng)，無(wú)需拘泥于采用二維平面激波、二維軸對(duì)稱(chēng)激波等簡(jiǎn)單幾何形狀，但其生成的激波可控性較差。后者作為一種反向設(shè)計(jì)，能夠?qū)λ思げㄐ螤钸M(jìn)行精細(xì)控制，但實(shí)現(xiàn)難度較大，尤其是三維情況下，現(xiàn)階段僅在某些特殊限制內(nèi)(如激波形狀為三維密切軸對(duì)稱(chēng)曲面，來(lái)流條件為均勻水平氣流等)實(shí)現(xiàn)了對(duì)曲面激波反問(wèn)題的求解，而一般意義上的三維曲面激波反設(shè)計(jì)問(wèn)題目前仍處于探索階段[20-21,29]。工程中，高超聲速飛行器的前體、進(jìn)氣道設(shè)計(jì)往往存在各種限制條件，導(dǎo)致兩者的氣動(dòng)一體化設(shè)計(jì)難度較大，靈活性亟需進(jìn)一步提高，有必要發(fā)展更多針對(duì)該問(wèn)題的求解方法。

本文首先在傳統(tǒng)密切軸對(duì)稱(chēng)方法基礎(chǔ)上，探討了針對(duì)非均勻來(lái)流條件的三維曲面激波設(shè)計(jì)反問(wèn)題，然后提出了一種針對(duì)該問(wèn)題的微元密切軸對(duì)稱(chēng)解法(Micro Osculating Axisymmetric flows，MOA)。為驗(yàn)證求解方法和相關(guān)程序的正確性和可行性，采用該方法重構(gòu)了一個(gè)4°攻角條件下的內(nèi)錐激波生成體，和一個(gè)10°外錐形流中的內(nèi)錐激波生成體，并分別與CFD數(shù)值仿真結(jié)果進(jìn)行了對(duì)比。最后討論了該方法在前體/進(jìn)氣道一體化設(shè)計(jì)中的應(yīng)用前景。

1 傳統(tǒng)密切軸對(duì)稱(chēng)方法

20世紀(jì)90年代，Sobieczky等證明了空間三維流動(dòng)可以用密切平面內(nèi)的軸對(duì)稱(chēng)二維流動(dòng)等效[15]，進(jìn)而由此提出了密切軸對(duì)稱(chēng)理論，解決了三維密切曲面激波的反設(shè)計(jì)問(wèn)題。

圖1以?xún)?nèi)錐形激波為例，闡釋了傳統(tǒng)密切軸對(duì)稱(chēng)理論的基本原理。沿著非軸對(duì)稱(chēng)內(nèi)錐激波周向，三維空間被離散為一系列二維切片，即密切平面，這些周向排列的密切平面與激波曲面垂直，與波前來(lái)流平行。同時(shí)，三維內(nèi)錐激波及波后三維流動(dòng)亦被視作這些密切平面上軸對(duì)稱(chēng)二維激波和二維流動(dòng)的疊加，而根據(jù)軸對(duì)稱(chēng)二維激波曲線反設(shè)計(jì)壓縮面型線可借助二維逆特征線法實(shí)現(xiàn)。

為避免各密切平面之間橫向流動(dòng)的干擾，每個(gè)密切平面內(nèi)軸對(duì)稱(chēng)二維流動(dòng)的激波曲線形狀均相同，即激波強(qiáng)度相同，從而保證了各密切平面內(nèi)的激波后沿壓力相等，這也是上文所述用軸對(duì)稱(chēng)二維流動(dòng)等效空間三維流動(dòng)的前提。

各密切平面上軸對(duì)稱(chēng)流動(dòng)的唯一區(qū)別在于其對(duì)稱(chēng)軸位置各不相同(如OAOA′，OBOB′等)，它們由基準(zhǔn)橫截面上激波曲線的局部曲率中心位置決定。實(shí)際上，這也正是密切軸對(duì)稱(chēng)流動(dòng)與軸對(duì)稱(chēng)流動(dòng)的差別所在：密切軸對(duì)稱(chēng)流動(dòng)是一種更為廣義的軸對(duì)稱(chēng)流動(dòng)，后者只是前者在橫截面上激波曲線形狀為圓形時(shí)的一種特例，此時(shí)，各密切平面即為相同的軸對(duì)稱(chēng)平面，各局部曲率半徑即為同一個(gè)圓弧半徑，各局部曲率中心即為同一個(gè)圓弧中心(見(jiàn)圖2)。

根據(jù)傳統(tǒng)密切軸對(duì)稱(chēng)理論基本原理可知，只有當(dāng)各密切平面之間不存在橫向速度與橫向壓力梯度時(shí)，這種采用周向布置二維密切平面組成三維空間的構(gòu)造方法才是精確的。這說(shuō)明密切軸對(duì)稱(chēng)理論仍然存在一定的局限性，其構(gòu)造的三維激波只能是某種密切曲面，而非真正的任意三維曲面，并且只能基于均勻水平來(lái)流條件。

圖1 傳統(tǒng)密切軸對(duì)稱(chēng)理論示意圖Fig.1 Schematic of traditional osculating axisymmetric flows concept

圖2 軸對(duì)稱(chēng)內(nèi)錐流動(dòng)Fig.2 Axisymmetric internal conical flow

如圖3所示，激波面為滿(mǎn)足密切幾何條件的內(nèi)錐曲面，密切平面OABO′即xy平面，點(diǎn)P位于激波曲面與密切平面OABO′的交線上，P點(diǎn)波前速度矢量為

V=Vx+Vy+Vz

(1)

式中：Vx、Vy、Vz分別為x、y、z方向的速度矢量。

根據(jù)當(dāng)?shù)丶げㄇ嫘螤詈蚏ankine-Hugoniot方程，當(dāng)波前為均勻水平來(lái)流時(shí)，Vy=Vz=0，波前速度V=Vx，波后速度方向仍在密切平面OABO′之內(nèi)；而當(dāng)波前為非均勻來(lái)流時(shí)，如Vz≠0，則波后速度必然存在垂直于密切平面OABO′的速度分量Vcross，若繼續(xù)以密切面OABO′上的軸對(duì)稱(chēng)二維流動(dòng)代替原空間三維流動(dòng)顯然是不合適的。這意味著傳統(tǒng)密切軸對(duì)稱(chēng)理論中構(gòu)造密切平

圖3 傳統(tǒng)密切軸對(duì)稱(chēng)理論適用性分析Fig.3 Applicability analysis of traditional osculating axisymmetric flows theory

面的方法失效，并且不難發(fā)現(xiàn)，在任意給定非均勻波前參數(shù)情況下，波后的三維流場(chǎng)中也不存在其他任何具備此前密切平面特征(例如無(wú)垂直方向上流動(dòng))的整張二維平面。

盡管如此，在三維激波曲面上的任意一點(diǎn)，仍能找到局部的微元密切平面，在該點(diǎn)附近的小范圍內(nèi)，依然能夠延續(xù)用軸對(duì)稱(chēng)二維流動(dòng)代替空間三維流動(dòng)的等效求解思路。將該方法命名為微元密切軸對(duì)稱(chēng)方法，具體原理和實(shí)施方式在第2節(jié)進(jìn)行詳細(xì)分析。

2 微元密切軸對(duì)稱(chēng)方法

如圖4所示，P11為激波面上一點(diǎn)，在非均勻來(lái)流下，其波前速度為某種任意指定大小和方向的Vpre，波后速度Vpost由當(dāng)?shù)丶げㄇ嫘螤詈蚏ankine-Hugoniot方程得到。本文所謂“微元密切面”(Micro Osculating Plane，MOP)之一即點(diǎn)P11附近由矢量Vpre和Vpost確定的微平面，由于Vpre方向的隨機(jī)性，微元密切面并非傳統(tǒng)密切理論中垂直于激波面的軸對(duì)稱(chēng)徑向平面。

根據(jù)微元密切面定義，P點(diǎn)速度方向恰好在微元密切面之內(nèi)，因此可以看作是一種局部的二維流動(dòng)，在該點(diǎn)局部應(yīng)用密切軸對(duì)稱(chēng)原理，即可將復(fù)雜三維流動(dòng)轉(zhuǎn)化為過(guò)該點(diǎn)的微元密切面上的軸對(duì)稱(chēng)二維流動(dòng)進(jìn)行求解。應(yīng)當(dāng)說(shuō)明，在非均勻來(lái)流條件下，按照上述方法生成的激波面上各個(gè)局部微元密切面是以一種相對(duì)無(wú)序的方式排布的，即使是沿流向分布的一系列微元密切面(如圖4中過(guò)點(diǎn)P11和點(diǎn)P12的微元平面)，其上的軸對(duì)稱(chēng)二維流動(dòng)也并非共用一根對(duì)稱(chēng)軸(如軸對(duì)稱(chēng)激波的對(duì)稱(chēng)軸OO′)，而是各自擁有相互獨(dú)立的“微元對(duì)稱(chēng)軸”(如圖4中的O1O2和O2Oi)。微元對(duì)稱(chēng)軸的空間位置和角度各異，是由周向相鄰兩個(gè)微元密切面的交線決定的，這些特征與傳統(tǒng)密切軸對(duì)稱(chēng)理論中整張密切平面對(duì)應(yīng)整條對(duì)稱(chēng)軸不同。微元密切面與微元對(duì)稱(chēng)軸的確定是微元密切軸對(duì)稱(chēng)方法的關(guān)鍵。

圖4 微元密切軸對(duì)稱(chēng)方法示意圖Fig.4 Schematic of micro osculating axisymmetric flows method

由于流場(chǎng)被離散為更多的微元密切面，而非傳統(tǒng)密切理論中數(shù)個(gè)周向排布的密切平面，圖5所示整張密切平面上的軸對(duì)稱(chēng)二維特征線法將難以直接應(yīng)用，但在各個(gè)微元密切平面內(nèi)，均為局部的軸對(duì)稱(chēng)二維流動(dòng)，其控制方程依舊是連續(xù)性方程(式(2))、動(dòng)量方程(式(3)～式(4))和聲速方程(式(5))：

(2)

(3)

(4)

(5)

式中：ρ為密度；p為壓力；Vx、Vy、Vz分別為x、y、z方向的分速度；c為聲速。

特征線方法將這些偏微分方程轉(zhuǎn)化為沿特征線的常微分方程，然后進(jìn)行求解。式(6)～式(7)和式(8)～式(10)分別給出了特征線方程和沿特征線的相容關(guān)系[30]:

(6)

(7)

dp-c2dρ=0

(8)

pVdV+dp=0

(9)

(10)

式中：λ為特征線斜率；θ為流動(dòng)角；μ為馬赫角；V為速度；Ma為馬赫數(shù)；下標(biāo)0和±分別表示流線和左行/右行特征線。

值得注意的是，與在整張軸對(duì)稱(chēng)二維平面上進(jìn)行推進(jìn)求解不同，圖5中前一條特征線的計(jì)算結(jié)果可以直接作為后續(xù)計(jì)算域的初值線(如P1Q1，P2Q2等)，但在圖6所示由一系列微元密切面拼接而成的三維扭轉(zhuǎn)曲面上，則必須對(duì)各微元面內(nèi)的數(shù)據(jù)進(jìn)行處理，即把前一個(gè)微元密切面(Micro Osculating Plane,MOP)內(nèi)結(jié)尾特征線(如MOP1內(nèi)的P1Q1，MOP2內(nèi)的P2Q2等)上的坐標(biāo)與流動(dòng)參數(shù)向后一個(gè)微元密切面當(dāng)?shù)囟S坐標(biāo)系內(nèi)投影變換，從而保證特征線計(jì)算能夠持續(xù)推進(jìn)。這種處理有效避免了三維特征線計(jì)算，三維情況下，沿特征面的相容關(guān)系不再是常微分方程，問(wèn)題將變得十分復(fù)雜，且有研究表明[31]，與有限差分方法相比，三維特征線法在計(jì)算效率上也不再具備特殊優(yōu)勢(shì)。

最終，沿流向排布的一系列微元密切面組成一張三維流面，如圖7所示，三維流動(dòng)曲面相比二

圖6 流向排布微元密切面上的特征線網(wǎng)格Fig.6 Mesh of characteristic lines on curved micro osculating surface in flow direction

圖7 三維密切曲面與傳統(tǒng)二維密切平面對(duì)比Fig.7 Comparison of traditional 2D osculating plane and 3D surface consisting of micro osculating planes

維徑向平面發(fā)生了明顯的偏轉(zhuǎn)，并且隨著流動(dòng)越向下游發(fā)展，偏離程度越大。

3 非均勻來(lái)流條件下三維曲面激波重構(gòu)

基于上述方法編寫(xiě)了非均勻來(lái)流下三維曲面激波構(gòu)建的設(shè)計(jì)程序，并以高超聲速推進(jìn)系統(tǒng)氣動(dòng)設(shè)計(jì)中的兩類(lèi)常見(jiàn)問(wèn)題為例，分析討論微元密切軸對(duì)稱(chēng)方法的設(shè)計(jì)結(jié)果，并與CFD數(shù)值仿真結(jié)果進(jìn)行對(duì)比。三維定常流場(chǎng)的無(wú)黏數(shù)值仿真采用了ANSYS Fluent的密度基求解器，其中歐拉方程的求解采用有限體積法，無(wú)黏通量計(jì)算采用二階AUSM (Advection Upstream Splitting Method)格式，氣體性質(zhì)選擇為理想氣體，比熱比γ為1.4。計(jì)算過(guò)程的收斂依據(jù)為各項(xiàng)殘差下降至少3個(gè)數(shù)量級(jí)并保持穩(wěn)定。計(jì)算所用網(wǎng)格由ANSYS ICEM生成，由于對(duì)稱(chēng)性，僅取模型的一半進(jìn)行網(wǎng)格劃分，并在各壓縮面前緣對(duì)激波根部區(qū)域進(jìn)行局部加密，總網(wǎng)格量約350萬(wàn)。

3.1 帶攻角來(lái)流下的內(nèi)錐曲面激波反設(shè)計(jì)

高超聲速飛行器常伴隨著帶攻角飛行，零攻角條件下設(shè)計(jì)得到的錐導(dǎo)乘波體構(gòu)型在有攻角來(lái)流下，前緣激波將偏離設(shè)計(jì)位置，導(dǎo)致不能完全乘波。本節(jié)針對(duì)這一實(shí)際問(wèn)題，采用微元密切軸對(duì)稱(chēng)方法，在來(lái)流馬赫數(shù)Ma∞=6，攻角α=4°條件下根據(jù)預(yù)設(shè)的軸對(duì)稱(chēng)內(nèi)錐曲面激波(對(duì)稱(chēng)軸水平)反設(shè)計(jì)了三維壓縮面。

設(shè)計(jì)時(shí)指定的內(nèi)錐激波母線形狀如圖8所示，曲線方程為

圖8 內(nèi)錐激波母線形狀Fig.8 Internal conical shock wave geometry

y=A0+A1x+A2x2+A3x3+A4x4+A5x5

(11)

式中:A0、A1、A2、A3、A4、A5為系數(shù)，具體數(shù)值為A0=0.581 2，A1=-0.242 5，A2=-0.172 7，A3= -0.402 2，A4=0.766 4，A5=-0.414 0。

圖9(a)給出了CFD計(jì)算得到的對(duì)稱(chēng)面流場(chǎng)，在4°攻角下，迎風(fēng)側(cè)與背風(fēng)側(cè)激波形狀仍能保持完全相同，并與設(shè)計(jì)之初給定的激波形狀(白色虛線)保持一致。圖9(b)為橫截面流場(chǎng)，圖中V∞為來(lái)流速度，其中的激波形狀亦保持為嚴(yán)格圓形。作為對(duì)照，圖10考察了以均勻水平來(lái)流設(shè)計(jì)的壓縮面在4°攻角下的表現(xiàn)，激波整體形狀出現(xiàn)了明顯的非軸對(duì)稱(chēng)性，橫截面上的波后參數(shù)也呈現(xiàn)出更大的非均勻性，這些現(xiàn)象均對(duì)乘波體設(shè)計(jì)不利。值得注意的是，利用微元密切軸對(duì)稱(chēng)方法，盡管激波形狀在攻角條件下能夠被設(shè)計(jì)為標(biāo)準(zhǔn)內(nèi)錐曲面，波后流場(chǎng)卻并非為軸對(duì)稱(chēng)內(nèi)錐流動(dòng)，而是一種包含顯著橫向速度和參數(shù)梯度的三維流動(dòng)。

為進(jìn)一步衡量計(jì)算準(zhǔn)確度，圖11給出了兩種方法得到的型面經(jīng)CFD計(jì)算得到的激波位置與預(yù)設(shè)位置的誤差曲線，其中Δh為實(shí)際激波位置與預(yù)設(shè)激波位置的縱坐標(biāo)之差，R為進(jìn)口半徑。可以看出，傳統(tǒng)二維方法由于未考慮波前來(lái)流的非均勻性，激波位置誤差達(dá)到了6%，而MOA方法的誤差最大約0.5%，準(zhǔn)確度提升了一個(gè)數(shù)量級(jí)。另外，圖12對(duì)比了采用微元密切軸對(duì)稱(chēng)方法和CFD計(jì)算得到的壁面壓力曲線(圖中p∞為來(lái)流壓力)，兩者基本吻合，迎風(fēng)側(cè)與背風(fēng)側(cè)的最大誤差分別為2.5%和4.3%。事實(shí)上，與其他數(shù)值解法類(lèi)似，微元密切軸對(duì)稱(chēng)方法的計(jì)算精度與在激波曲面上選取微元密切面的數(shù)量有關(guān)，本例中沿流向和周向一共布置了25×19個(gè)微元密切面，隨著離散程度的提升(計(jì)算量也會(huì)相應(yīng)增加)，計(jì)算誤差還會(huì)進(jìn)一步減小。

圖9 利用MOA方法在攻角條件下生成的內(nèi)錐激波的數(shù)值仿真驗(yàn)證Fig.9 Numerical validation of internal conical shock wave designed by MOA in incoming flow with angle of attack

圖10 軸對(duì)稱(chēng)壓縮面在攻角條件下的流場(chǎng)Fig.10 Flowfield of axisymmetric compression wall in incoming flow with angle of attack

圖11 流向激波位置的相對(duì)誤差曲線Fig.11 Relative error curves of shock wave location along flow direction

圖12 MOA與CFD得到的壁面壓力分布對(duì)比Fig.12 Comparison of MOA and CFD results of wall pressure distribution

3.2 外錐形流中的內(nèi)錐激波反設(shè)計(jì)

在高超聲速飛行器前體/進(jìn)氣道一體化設(shè)計(jì)時(shí)，進(jìn)氣道有時(shí)需被放置在經(jīng)歷了前體預(yù)壓縮之后的流場(chǎng)內(nèi)，這意味著進(jìn)氣道入口處氣流不再是常規(guī)設(shè)計(jì)中的遠(yuǎn)前方均勻來(lái)流，因此必須考慮非均勻來(lái)流情況下的進(jìn)氣道設(shè)計(jì)。

基于外錐形流的乘波前體與基于內(nèi)錐形流的內(nèi)轉(zhuǎn)進(jìn)氣道是現(xiàn)階段的研究熱點(diǎn)，圖13給出了二者基準(zhǔn)流場(chǎng)的一種耦合方式，其中，內(nèi)錐激波曲面被包裹在外錐激波以?xún)?nèi)的部分激波面的波前參數(shù)由外錐流動(dòng)決定。本節(jié)針對(duì)這種內(nèi)外錐基準(zhǔn)流場(chǎng)一體化構(gòu)型，采用微元密切軸對(duì)稱(chēng)方法，在來(lái)流馬赫數(shù)Ma∞=6、半錐角ψ=10°的外錐流場(chǎng)中根據(jù)軸對(duì)稱(chēng)內(nèi)錐曲面激波反設(shè)計(jì)了三維壓縮面。內(nèi)錐激波母線形狀仍由式(11)表達(dá)，其中系數(shù)更換為A0=7.600×10-1，A1=-6.653×10-1，A2=2.131× 10-1，A3=-2.281×10-1，A4=7.317×10-2，A5=-9.328×10-4。

圖13 外錐與內(nèi)錐流場(chǎng)的耦合設(shè)計(jì)示意圖Fig.13 Schematic of integrated design of external and internal conical flows

設(shè)計(jì)結(jié)果如圖14所示，其中包括用于產(chǎn)生所需非均勻外錐形流的前體錐面，和采用微元密切軸對(duì)稱(chēng)方法得到的內(nèi)錐激波壓縮面。在前體/進(jìn)氣道一體化設(shè)計(jì)中，進(jìn)氣道的唇口型線通常以橫截面上的激波形狀為基準(zhǔn)，從而實(shí)現(xiàn)進(jìn)氣道的激波封口和全流量捕獲，因此，橫截面上的激波形狀需要額外關(guān)注，圖15在正視圖視角下展示了所設(shè)計(jì)內(nèi)錐激波形狀在橫截面上的準(zhǔn)確度。由于2個(gè)基準(zhǔn)流場(chǎng)不同心，前方外錐形流在內(nèi)錐基準(zhǔn)流場(chǎng)柱坐標(biāo)系下具有明顯的周向速度，即三維流動(dòng)特征，但采用密切軸對(duì)稱(chēng)方法得到的內(nèi)錐激波在橫截面上的形狀依舊保持為嚴(yán)格圓弧形(與標(biāo)準(zhǔn)圓弧虛線重合)。為進(jìn)一步衡量其準(zhǔn)確度，圖16給出了各橫截面上CFD計(jì)算得到的內(nèi)錐激波位置與預(yù)設(shè)位置的誤差曲線，其中橫坐標(biāo)θ為軸對(duì)稱(chēng)內(nèi)錐激波的周向角度(對(duì)稱(chēng)面為0°)，L為內(nèi)錐激波流向總長(zhǎng)度。在前方外錐激波的包裹區(qū)域內(nèi)部，內(nèi)錐激波各位置誤差均不超過(guò)0.5%；當(dāng)θ超出極限角度(內(nèi)、外錐激波在橫截面上的交點(diǎn)所處角度)后，由于內(nèi)錐激波的波前流動(dòng)參數(shù)不再是馬赫數(shù)為6、10°半錐角的錐形流場(chǎng)參數(shù)，而是遠(yuǎn)前方均勻來(lái)流，即與設(shè)計(jì)時(shí)的輸入條件不同，故激波形狀開(kāi)始出現(xiàn)顯著偏差。

綜上所述，采用本文所提出的微元密切軸對(duì)稱(chēng)方法針對(duì)2種非均勻來(lái)流重構(gòu)的內(nèi)錐激波形狀均與數(shù)值仿真結(jié)果完全吻合，從而驗(yàn)證了本文微元密切軸對(duì)稱(chēng)方法的正確性與可行性。還應(yīng)指出，應(yīng)用圖13所示內(nèi)/外錐一體化基準(zhǔn)流場(chǎng)，結(jié)合流線追蹤技術(shù)，可以方便地設(shè)計(jì)出流向布局的三維乘波前體/內(nèi)轉(zhuǎn)進(jìn)氣道氣動(dòng)一體化構(gòu)型：在內(nèi)錐激波面上選取內(nèi)轉(zhuǎn)式進(jìn)氣道前緣線，然后分別在前方外錐流場(chǎng)與后方內(nèi)錐流場(chǎng)中向前/向后流線追蹤，即可得到相應(yīng)的進(jìn)氣道壁面與乘波前體下表面。而在傳統(tǒng)一體化方法中，由于內(nèi)轉(zhuǎn)式進(jìn)氣道(基準(zhǔn)流場(chǎng))基于均勻來(lái)流設(shè)計(jì)，須以平均流動(dòng)參數(shù)代替非均勻的前體波后流場(chǎng)進(jìn)行設(shè)計(jì)[32]，或?qū)?nèi)轉(zhuǎn)進(jìn)氣道直接置于前體前緣捕獲自由來(lái)流(此時(shí)前體位于進(jìn)氣道兩側(cè)，相當(dāng)于周向布局的一體化構(gòu)型)[33-34]。因此，采用基于非均勻來(lái)流下三維激波反設(shè)計(jì)方法得到的內(nèi)外錐耦合流場(chǎng)，能夠顯著增強(qiáng)前體/進(jìn)氣道氣動(dòng)一體化設(shè)計(jì)的靈活性，相關(guān)工作值得進(jìn)一步深入研究。

圖14 內(nèi)外錐耦合流場(chǎng)的數(shù)值仿真結(jié)果Fig.14 Numerical simulation results of integrated outward/inward turning flowfields

圖15 利用MOA方法在外錐流場(chǎng)中生成的內(nèi)錐激波的數(shù)值仿真驗(yàn)證Fig.15 Numerical validation of internal conical shock wave designed by MOA method in external conical incoming flow

圖16 周向激波位置的相對(duì)誤差曲線Fig.16 Relative error curves of shock wave location along circumferential direction

4 結(jié) 論

本文討論了針對(duì)非均勻來(lái)流條件的三維曲面激波反問(wèn)題，提出了針對(duì)該問(wèn)題的微元密切軸對(duì)稱(chēng)方法并進(jìn)行了數(shù)值仿真驗(yàn)證。通過(guò)本文研究，得到以下結(jié)論：

1) 所提出的三維曲面激波反問(wèn)題的微元密切軸對(duì)稱(chēng)方法，解除了傳統(tǒng)密切方法中不能有橫向流動(dòng)的限制條件，可用于非均勻來(lái)流下的三維曲面激波反設(shè)計(jì)。

2) 基于微元密切軸對(duì)稱(chēng)方法分別求解了帶攻角來(lái)流與外錐形流來(lái)流2種來(lái)流條件下的標(biāo)準(zhǔn)內(nèi)錐曲面激波及波后流場(chǎng)，數(shù)值仿真結(jié)果證實(shí)了所設(shè)計(jì)的三維曲面激波形狀與預(yù)設(shè)形狀完全吻合。

3) 該方法為高超聲速飛行器三維乘波前體、內(nèi)轉(zhuǎn)式進(jìn)氣道及二者的氣動(dòng)一體化設(shè)計(jì)提供了新的解決途徑，展示出在該領(lǐng)域內(nèi)的重要應(yīng)用前景。