非均勻來流下三維激波反問題的微元密切軸對稱解法

周航，金志光

南京航空航天大學 能源與動力學院, 南京 210016

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

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

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

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

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

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

1 傳統(tǒng)密切軸對稱方法

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

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

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

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

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

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

圖2 軸對稱內(nèi)錐流動Fig.2 Axisymmetric internal conical flow

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

V=Vx+Vy+Vz

(1)

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

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

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

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

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

2 微元密切軸對稱方法

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

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

圖4 微元密切軸對稱方法示意圖Fig.4 Schematic of micro osculating axisymmetric flows method

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

(2)

(3)

(4)

(5)

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

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

(6)

(7)

dp-c2dρ=0

(8)

pVdV+dp=0

(9)

(10)

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

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

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

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

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

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

3 非均勻來流條件下三維曲面激波重構

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

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

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

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

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

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

圖10 軸對稱壓縮面在攻角條件下的流場Fig.10 Flowfield of axisymmetric compression wall in incoming flow with angle of attack

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

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

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

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

基于外錐形流的乘波前體與基于內(nèi)錐形流的內(nèi)轉(zhuǎn)進氣道是現(xiàn)階段的研究熱點，圖13給出了二者基準流場的一種耦合方式，其中，內(nèi)錐激波曲面被包裹在外錐激波以內(nèi)的部分激波面的波前參數(shù)由外錐流動決定。本節(jié)針對這種內(nèi)外錐基準流場一體化構型，采用微元密切軸對稱方法，在來流馬赫數(shù)Ma∞=6、半錐角ψ=10°的外錐流場中根據(jù)軸對稱內(nèi)錐曲面激波反設計了三維壓縮面。內(nèi)錐激波母線形狀仍由式(11)表達，其中系數(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)錐流場的耦合設計示意圖Fig.13 Schematic of integrated design of external and internal conical flows

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

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

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

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

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

4 結(jié) 論

本文討論了針對非均勻來流條件的三維曲面激波反問題，提出了針對該問題的微元密切軸對稱方法并進行了數(shù)值仿真驗證。通過本文研究，得到以下結(jié)論：

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

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

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