王 丁，王江峰，*，李龍飛

(1.非定?？諝鈩恿W(xué)與流動控制工業(yè)和信息化部重點實驗室，南京 210016；2.南京航空航天大學(xué) 航空學(xué)院，南京 210016)

0 引 言

乘波氣動布局由于其高升阻比、方便進行內(nèi)外流一體化設(shè)計等優(yōu)勢在高超聲速滑翔飛行器[1-2]及巡航飛行器[3-4]氣動設(shè)計領(lǐng)域得到了廣泛關(guān)注。由于基準流場前緣激波依賴區(qū)的設(shè)計決定了飛行器腹部下方的流場特性及乘波特性，因此前緣激波依賴區(qū)流場的設(shè)計是乘波布局設(shè)計過程中的重點之一。盡管基于傳統(tǒng)二維特征線法及吻切理論的激波流場反問題的求解在近年來已得到充分的研究，但該方法無法有效地求解含有顯著三維流動效應(yīng)的流場，因此針對三維曲面激波流場反問題的求解仍有待進一步的研究。

根據(jù)基準流場設(shè)計原理的不同，將以往設(shè)計方法分為兩類進行介紹。第一類方法將三維流場簡化為二維流場切片的疊加，屬于準三維方法。早期研究所用的基準流場主要包括斜劈流場[5]及圓錐激波流場[6]，近年來吻切理論在此類方法中占據(jù)主導(dǎo)地位。Sobieczky等[7]提出了吻切錐設(shè)計方法，該方法將三維流場等效為吻切面上二維錐形流場切片的疊加，隨后又將吻切面流場進一步拓展為任意軸對稱流場切片，即吻切軸對稱法[8]，并在內(nèi)外流一體化設(shè)計領(lǐng)域得到了廣泛應(yīng)用，例如全乘波內(nèi)外流一體化布局[9-10]、吻切內(nèi)錐前體進氣道一體化布局[11]等。以上方法中，每個吻切面共用同一種流場構(gòu)型，而吻切流場法[12-13]在每個吻切面內(nèi)獨立設(shè)計各自的二維流場，進一步提高了設(shè)計的靈活性，該方法近年來在寬速域滑翔飛行器和內(nèi)外流一體化飛行器設(shè)計領(lǐng)域得到了大量應(yīng)用，例如展向變馬赫數(shù)乘波體[14-16]、展向變激波角乘波體[17]、多級壓縮內(nèi)外流一體化乘波布局[18-19]、雙乘波內(nèi)外流一體化布局[20-21]。前述方法在設(shè)計過程中忽略了流場中的三維效應(yīng)，為了增強對三維流動效應(yīng)的模擬能力，實現(xiàn)基于復(fù)雜三維激波流場的乘波布局設(shè)計，Zheng等[22]提出了多吻切錐設(shè)計方法，該方法將吻切面內(nèi)流場處理為多個不共軸的圓錐流場的組合，從而顯著提高了對帶有較大橫向壓力梯度流場的模擬能力。隨后，使用流面替代吻切面實現(xiàn)對三維流場的逼近，進一步提出了局部偏轉(zhuǎn)吻切錐方法[23]，該方法對帶側(cè)滑橢圓錐激波流場、變截面橢圓錐流場的三維流動效應(yīng)均表現(xiàn)出了良好的模擬能力，但該方法與真正意義上的三維模擬仍有一定區(qū)別。文中僅給出了其設(shè)計得到的壁面壓力分布與數(shù)值模擬的對比（最大誤差在2%左右），但未披露壁面馬赫數(shù)的計算誤差。

第二類方法是全三維方法，對于已知三維壁面形狀的三維基準流場正問題[24-25]來說，求解的關(guān)鍵在于對三維激波的準確模擬，相關(guān)數(shù)值模擬方法主要分為激波裝配法和激波捕捉法兩類。激波裝配法[26]通過對激波面進行追蹤并利用Rankine-Hugoniot關(guān)系式求解激波前后信息，保證了激波模擬的準確度，但激波面的追蹤過程涉及復(fù)雜的網(wǎng)格操作，因此與非結(jié)構(gòu)網(wǎng)格相結(jié)合以增強激波裝配算法的通用性和魯棒性是該方法發(fā)展的重點之一[27-28]。對于激波捕捉法而言，激波面由計算格式在給定的網(wǎng)格中進行模擬捕捉，因此激波形狀無法被精確描述。近似求解一維Riemann問題的通量格式（例如Roe格式及AUSM系列格式等）在低階激波捕捉格式中應(yīng)用最為廣泛，此類格式計算效率高且物理意義明確，但計算結(jié)果中出現(xiàn)的“紅寶石”現(xiàn)象、過熱現(xiàn)象、全速域計算問題[29-30]及基于維度分裂向多維格式推廣時忽略了單元切向物理信息[31]等問題仍有待進一步的研究。此外，基于疊波法的Godunov型大時間步長格式[32]也值得進一步關(guān)注，該方法充分考慮了精確解的依賴域，通過波束的線性疊加假設(shè)實現(xiàn)了CFL（Courant-Friedrichs-Lewy）數(shù)的顯著提高，且具有對間斷模擬的精度隨著CFL數(shù)提高而增加的優(yōu)點。提高激波分辨率的有效途徑之一是使用低數(shù)值耗散的高階格式，目前主要包括高階有限差分格式（例如WENO[33]格式）、基于有限體積的高階格式（如有限體積高階WENO格式[34]）和有限元類高階格式（例如DG[35]格式），文獻[36]對不同高階格式的特點及發(fā)展給出了詳細討論?；谌S特征線理論的數(shù)值算法是三維曲面激波反問題的主要求解方法之一。目前相對較為實用的三維特征線方法主要包括兩種：一種是校正六面體雙特征線法[37]，該方法基于三個初值點，通過求解三個馬赫錐的交點得到待求點，Huang等[38]給出了該方法在非均勻來流曲激波流場設(shè)計上的應(yīng)用，但未詳細評估該方法的計算誤差，且該方法目前僅限于求解激波依賴區(qū)流場，尚未見到應(yīng)用于基準流場其他模塊（如壓縮區(qū)等）的設(shè)計；另一種是五面體雙特征線方法[39-40]，該方法通過沿流線進行時間積分得到待求點，由于該方法要求初值面和求解面相互平行，壁面及激波形狀均已知，且需要構(gòu)造在局部坐標系中分別相差90°的四根雙特征線，因此主要應(yīng)用在內(nèi)流道流場[41-42]及激波流場[43]正問題的求解中，并不適用于激波流場反問題的求解。

為了克服傳統(tǒng)吻切理論在設(shè)計全三維曲面激波流場時的缺陷，實現(xiàn)對全三維復(fù)雜激波流場反問題的求解，本文基于三維特征線理論重構(gòu)了描述三維流動的控制方程，并在此基礎(chǔ)上發(fā)展了用于求解三維激波依賴區(qū)流場的全三維設(shè)計方法。采用圓錐激波流場、橢圓錐激波流場、小攻角圓錐激波流場和由Bezier曲面描述的一般性曲面激波流場等四個典型算例進行了方法驗證。通過與二維特征線方法對比，分析了本文設(shè)計方法的計算誤差及并行效率。與數(shù)值模擬結(jié)果對比，驗證了本文設(shè)計方法對由橫向壓力梯度及攻角引起的三維流動效應(yīng)的模擬能力。

1 控制方程

依據(jù)三維特征線理論[37]，馬赫錐及流線上的特征速度滿足圖1所示的幾何關(guān)系，其中馬赫錐由點S2發(fā)出，流線 S1S2為馬赫錐的軸線，若線段 S1S2的長度為流動速度V的大小，則線段 S1B1的長度為聲速a，其方向與特征法矢n一 致，點 S1和點 A1位于垂直于馬赫錐軸線的平面上，線段 S1A1的長度為特征聲速c，線段 A1S2的長度為特征速度Vc（式1）的大小。若定義 ?s為沿流線的變化率算子（式2）， ?c為沿馬赫線的變化率算子（式3），則對描述超聲速有旋無黏流動的控制方程（式4～8）中的微分算子進行特征分解（式2和式9），可得沿馬赫線的控制方程（式10～13）以及沿流線的控制方程（式14～18）。

圖1 速度關(guān)系與馬赫錐Fig.1 Velocity relation and Mach cone

2 基本單元的構(gòu)造及求解方法

本節(jié)通過構(gòu)造基本單元求解特定點的壓力P、密度 ρ 、速度矢量V、壓力梯度 ?P、速度張量 ?V等17個物理參數(shù)。

2.1 基本單元的構(gòu)型及相關(guān)計算幾何算法

圖2給出了基本單元構(gòu)型及其求解流程。圖中，初值線 D1D2和 D3D4為已知的三維網(wǎng)格線，綠點的物理參數(shù)均已知，藍點 S2為待求點，基本單元包含4條馬赫線和1條流線。在單元求解過程中共涉及6個基礎(chǔ)算法，下面分別進行介紹。

圖2 基本單元構(gòu)型及其求解流程Fig.2 Basic cell and the solving procedure

2.1.1 線性插值及特征線平均物理參數(shù)計算方法

假設(shè)線段兩端點分別為點1和點2，插值點3位于點1和點2之間，則可用式（19～20）所示的線性插值公式計算點3的物理參數(shù)，其中d31為點3到點1的距離，d21為點2到點1的距離。在計算特征線平均物理參數(shù)時，插值系數(shù)取 ξ =0.5。

2.1.2 特征線速度矢量及法矢的計算方法

參考圖3，若約定變量上方加雙橫線代表對矢量進行歸一化處理，則流線的方向矢量l可以表示為式（21），馬赫線（以線段 S2A3為例）的特征速度、方向矢量l和法矢n可 分別表示為式（22～24），其中dA3S1為點 A3到點 S1的矢量，但對于圖2中的馬赫線 S2A1來說則是點 S1到點 A1。

圖3 特征速度計算方法Fig.3 Calculation method of the characteristic velocity

2.1.3 Ray_X_Ray算法

該算法用于求解空間兩射線的交點坐標。參考圖4，以計算兩射線 S2A1和 S2A3的交點 S2為例，首先計算兩射線 S2A1和 S2A3的方向矢量l，則交點 S2的坐標為式（25），如果點 S2位于對稱面邊界上，則還需要對計算得到的坐標進行對稱面邊界校正。這里S2、A1等黑斜體英文字母均表示坐標系原點到對應(yīng)點的距離矢量。

圖4 Ray_X_Ray算法原理Fig.4 Principle of the Ray_X_Ray algorithm

2.1.4 Ray_X_Line算法

該算法用于求解空間射線和直線段的交點。參考圖5，以計算射線 S2S1和線段 A1A3的交點 S1為例，首先計算射線 S2S1的方向矢量和線段 A3A1的矢量，點S1物理參數(shù)的初值為點 A3和點 A1物理參數(shù)的代數(shù)平均（式19， ξ =0.5），利用式（26）計算點 S1的坐標，并由式（27）更新插值系數(shù) ξ，最后利用線性插值對點 S1的坐標及物理參數(shù)進行更新。通過迭代上述步驟，直到點 S1的坐標變化小于容差 ε ，本文所有算法所用的容差統(tǒng)一取 ε =1.0×10-6。

圖5 Ray_X_Line算法原理Fig.5 Principle of the Ray_X_Line algorithm

2.1.5 對稱面邊界相關(guān)算法

對稱面邊界由一個對稱面上任意選取的參考點Qs及 法矢n確 定，n的方向定義為由流場計算域內(nèi)部指向外部。對稱邊界條件包括對稱面處點參數(shù)的校正及鏡像點的計算，下面分別進行介紹。

2.1.5.1 對稱面鏡像點坐標及物理參數(shù)的計算

假設(shè)已知流場內(nèi)側(cè)一點 QL，若dLs為 點 QL到點Qs的 距離矢量，則該點關(guān)于對稱面的鏡像點 QR的坐標QR、速度VR、速度張量 ?VR（若z平面為對稱面）、速度張量 ?VR（若y平面為對稱面）、壓力梯度 ?PR，分別按式（28～32）計算，其中矩陣右下角的“L”表示矩陣內(nèi)的變量取點 QL的值。

2.1.5.2 對稱面上點的坐標及物理參數(shù)的校正

對于對稱面上任一點 Qr，若drs為 點 Qr到點 Qs的距離矢量，則該點校正后的坐標Qrc、 速度Vrc、速度張量?Vrc（ 若z平面為對稱面）、速度張量 ?Vrc（若y平面為對稱面）、壓力梯度 ?Prc，分別按式（3 3～37）計算，其中矩陣右下角的“r”表示矩陣內(nèi)的變量取點 Qr的值。

2.1.6 Cone_X_Curve算法

該算法用于求解已知點發(fā)出的馬赫錐和初值線的交點，圖6為交點的求解流程。參考圖7，以交點A4的計算為例，初值線為特征網(wǎng)格線 D1D2（點數(shù)為Nspan）。如果存在對稱面邊界，則初值線在對稱面處需要構(gòu)造Nextend（式38）個鏡像點，來統(tǒng)一對稱面附近基本單元的計算（如圖8所示）。式（39）為馬赫線A4S2的特征矢量lA4S2在平面 A4S2S1上的投影。式（40）為矢量dA4S2和投影方向矢量lproj分別與矢量dS1S2的夾角的差值，用來評估當前點與真實交點的誤差。假設(shè)點 S1在初值線上的指標為istart，沿初值線由點 S1向點 D1方向的搜索步長為istep=1，利用istart和istart+istep點 按插值系數(shù) ξ =1.0×10-4進行線性插值，得到istart點在初值線上極小鄰域內(nèi)的一點，并用該點計算 d θstart，且有 d θstart＜0。交點的求解過程分為兩步：1）由點 S1沿初值線向點 D1方向搜索，在每點計算式（39～40），直到相鄰兩點 d θi·dθi-istep＜0，則交點位于該區(qū)間內(nèi)；2）按照式（41～42）進行加權(quán)插值，其中初始的 ξL=0對應(yīng)于初值線上指標為i-istep的點，ξR=1對應(yīng) 于 指標為i的 網(wǎng) 格點 ， 在迭 代 過 程中 ， ξL、ξR、 d θL和 d θR按照圖6所示流程進行更新，直到收斂。

圖6 Cone_X_Curve算法求解流程Fig.6 Flow chart of the Cone_X_Curve algorithm

圖7 Cone_X_Curve算法原理Fig.7 Principle of the Cone_X_Curve algorithm

圖8 初值線鏡像點的構(gòu)造Fig.8 Construction of the mirror points of the initial value curve

2.1.7 基本單元的求解流程

基于以上6個基本算法，基本單元求解的流程可以表述為：1）用點 S1的物理參數(shù)給點 S2賦初值，該兩點為流線的兩端點；2）利用Ray_X_Ray算法計算流線 S1S2和馬赫線 A1S2交點 S2的坐標；3） A3的物理參數(shù)的初值采用點 A1和點 S1的平均值，然后利用Ray_X_Line算法計算馬赫線 S2A3和線段 S1A1的交點A3的坐標，并利用線性插值計算點 A3的物理參數(shù)；4）利用Cone_X_Curve算法插值求解點 A2和點 A4的坐標及物理參數(shù)；5）求解線性方程組，如果該點位于對稱面邊界上，則需要利用對稱面邊界對求解結(jié)果進行校正；6）重復(fù)步驟b～e直到待求點的空間坐標、壓力和馬赫數(shù)在兩次迭代間的變化小于給定的容差 ε，或總迭代數(shù)大于預(yù)設(shè)值imax。最終構(gòu)成的線性方程組由17個方程構(gòu)成，如表1所示，其中馬赫線 S2A1和S2A3按式（10～13）分別提供4個方程，沿馬赫線 S2A2和S2A4的式（10～13）相減提供4個方程，流線 S1S2按式（14～18）提供最后的5個方程。

表1 基本單元的線性方程組Table 1 Linear equations of the basic cell

2.2 控制方程的離散及線性方程組的求解

2.2.1 控制方程的離散

首先介紹控制方程（10～13）和（14～18）的離散：沿馬赫線由任一點1到點2（待求點）對微分方程（10～13）進行時間積分，并將特征算子 ?c轉(zhuǎn)化為差分項可得式（43～46）；沿流線由任一點1到點2（待求點）對微分方程（14～18）進行時間積分，并將特征算子 ?s轉(zhuǎn)化為差分項可得式（47～51）。

式（43～51）中： δt表示變量t的一個增長量；變量上方的單橫線表示該段特征線上的參數(shù)平均值（式19）。將式（43～49，51）中的時間積分項處理為（52～59）。其中，馬赫線和流線的時間項分別按式（60～61）計算；δl為特征線的長度，且定義特征線上的速度矢量方向為由初值點指向待求點時 δl為正值，否則為負值。離散后總的控制方程見附錄。

下面介紹三種不同的 ω 取值方案。

2.2.1.1 “ ω 01”方案（ ω1=0， ω2=1）

該方案假設(shè)速度張量 ?V和壓力梯度 ?P沿特征線為常數(shù)，并用待求點的速度張量 ?V和壓力梯度 ?P作為整個特征線上的值來進行時間積分。

2.2.1.2 “ ω 12”方案（ ω1=1/2， ω2=1/2）

該方案假設(shè)速度張量 ?V和壓力梯度 ?P沿特征線隨運動時間為線性變化。

2.2.1.3 “ ω -V”方案（ ω1和 ω2為特征速度和距離的函數(shù)）

將沿馬赫線和流線的特征速度的大小分別表示為式（62～63），那么速度張量 ?V和壓力梯度 ?P的各分量（統(tǒng)一用 Λ表示）沿特征線上點1到點2的時間積分可以表示為式（64～68），該式的解析解為式（69～70）。為了加強數(shù)值穩(wěn)定性，將式（69～70）統(tǒng)一用式（69）的四階泰勒展開（式71～73）表示。

當b2≠0時，

當b2=0時：

“ ω 12”和“ ω -V”方案中均用到了初值點1的速度張量 ?V和壓力梯度 ?P，在流場的計算過程中，采用最小二乘法進行計算。假設(shè)在流場中與點1由網(wǎng)格線直接相連的點有m個，若記點1的待求變量為?Λ0，則由式（74～78）所示的最小二乘法計算 ? Λ0。其中，di為點1到第i個相鄰點的距離矢量， Φi為最小二乘法的殘差函數(shù)，f為最小二乘法的目標函數(shù)。

2.2.2 線性方程組的求解

將表1中的線性方程組表示為Ax=b，其具有如下特點：1）系數(shù)矩陣為病態(tài)矩陣，使得A和b的小誤差會導(dǎo)致x產(chǎn)生較大的誤差；2）系數(shù)矩陣的最小奇異值（第17個奇異值，線性方程組的秩為17）與次小奇異值（第16個）相差較大。以上兩個特點使得該線性方程組的求解存在較大的穩(wěn)定性問題，需要進行特殊處理。

本文基于Tikhonov正則化和Lagrange插值方法發(fā)展了適用于該線性方程組求解的Tikhonov-Lagrange擬合法。若線性方程組系數(shù)矩陣的奇異值分解為式（79），則參考Tikhonov正則化方法[44]，本文構(gòu)造的方程組的解為式（80），其中對角矩陣 Σ*的第i個對角項為式（81），r為調(diào)節(jié)指數(shù)，當α為0時，方程組的解恢復(fù)為未經(jīng)任何修正的原始解。以第4.1節(jié)圓錐激波算例中的第一個基本單元的求解為例，時間積分 選 “ ω 01” 方 案 ， 則 求 解 的 變 量P、 ρ和V隨α在間的變化曲線如圖9左側(cè)所示。為了防止坐標軸過多產(chǎn)生干擾，圖中隱去了5個變量的y軸。當r=1時，P、ρ和w曲線是近似重合的；當r=2時，u和v曲線也近似重合。r=1 時的速度分量u和w在α=-Σ1r7附近出現(xiàn)震蕩，則系數(shù)矩陣的奇異值分解的小誤差會引起所求物理量的較大誤差，這也是使用高斯消去法直接求解此類線性方程組會有嚴重的穩(wěn)定性問題的原因，但這一震蕩情況在r=2得到顯著改善。對r=2的變化曲線的光滑部分進行采樣，建立P、ρ和V的Lagrange插值曲線，則擬合曲線在α=0處的值即為Tikhonov-Lagrange擬合法得到的線性方程組的解。本文在α=0的左右兩側(cè)按照式（82～85）各取對稱的7個點，當采樣點α確定時，對應(yīng)的解由式（80～81）解出，從而完成5個變量（P、ρ和V）的Lagrange插值曲線的擬合。由于壓力梯度和速度張量不是流場設(shè)計中的關(guān)心變量，因此不進行擬合求解。

圖9 典型變量隨系數(shù)的變化曲線Fig.9 Variation curves of typical variables with the coefficient

3 三維曲面激波依賴區(qū)設(shè)計方法

本節(jié)包含曲面激波的離散及網(wǎng)格面的推進算法兩部分內(nèi)容。

3.1 曲面激波的離散及數(shù)值解存在的必要條件

本文關(guān)注的重點在于激波流場求解算法本身，因此不對存在反問題解的激波曲面的數(shù)學(xué)表達式及最優(yōu)網(wǎng)格單元劃分方法做深入探討。簡單起見，對于圓錐及橢圓錐等曲面激波直接采用三維建模軟件生成模型，并在網(wǎng)格劃分軟件中離散成均勻的四邊形網(wǎng)格。對網(wǎng)格點周圍相鄰的四個四邊形網(wǎng)格單元的法矢進行反距插值可以得到網(wǎng)格點處的法矢。對于復(fù)雜的激波曲面算例采用3 × 1階Bezier曲面進行準確描述（見式86～88），其中 P0,0等為Beizer曲面的控制點。將Bezier曲面均勻離散為Nu×Nv個點（如圖10所示），則第（i,j）個網(wǎng)格點的參數(shù)坐標為（iδu,jδv），參數(shù)步長的定義見式（89～90），網(wǎng)格點處單位法矢的計算見式（91）。對于離散后的激波網(wǎng)格，網(wǎng)格點的波后物理參數(shù)可由Rankine-Hugoniot關(guān)系式（式92～98）計算得到，其中下標1代表波前參數(shù)，下標2代表波后參數(shù)，下標n代表激波曲面的法向，下標t代表激波曲面的切向。

對于得到的離散激波網(wǎng)格，激波反問題有解需要保證兩個條件，從而保證待求點可解：1）激波流場本身是存在的且流場內(nèi)馬赫波不會匯聚成激波；2）流場內(nèi)任一點的基本單元均能夠成功建立。下面針對這兩個條件進行展開討論。

要使激波流場是存在的，則需要使激波曲面上任一點的波前法向馬赫數(shù)大于1，且該點的激波角 β 小于脫體激波角。參考圖10，由波前法向馬赫數(shù)大于1可得式（99），其中激波角由式（100）計算。該點的脫體激波角對應(yīng)于波后氣流偏折角 θ的極值，即式（101），因此激波角的范圍可表示為式（102～105）。此外，激波曲面波后任一點向下游發(fā)射的馬赫錐必須與激波曲面相交，該條件保證了圖2中逆馬赫線A1S2存在，該條件可表示為式（106）。代入波后馬赫數(shù)及對應(yīng)的馬赫角α2的計算公式（107～108），可化簡為式（99），即只要波前法向馬赫數(shù)大于1，則該條件成立。要使流場內(nèi)馬赫波不會匯聚成激波，需要約束激波曲面在兩個相互正交的典型方向上的曲率，參考圖11，第一個方向是來流速度V1及激波曲面法矢n組成的平面與激波曲面的交線方向，由正交關(guān)系可得第二個方向。以第一個方向為例討論激波曲率對馬赫線匯聚點位置的影響。點 A2為點 A1在激波曲面上極小鄰域內(nèi)的一點，兩點間距為 δs， 激波曲面在點 A1處的曲率中心為點 C1，曲率半徑為R，兩點發(fā)出的逆馬赫線交于點 D1（若兩條逆馬赫線平行，則交點在無窮遠處），線段 A1D1的長度為dA1D1，則由圖11中的幾何關(guān)系可得式（109～110），以上兩式按照式（111～112）的推導(dǎo)可得式（113），其中波后馬赫角α2-A1及氣流偏折角 θA1的計算公式見式（108，101）。顯然兩馬赫線交點的距離與當?shù)丶げㄇ拾霃匠烧?。若該交點落在流場內(nèi)，則流場內(nèi)存在馬赫波匯聚為激波的現(xiàn)象，因此該式對激波曲率半徑的最小值進行了約束。對于反問題而言，由于壁面位置在求解前是未知的，因此如何在該條件的基礎(chǔ)上發(fā)展激波形狀的顯式約束表達式是后續(xù)研究的重點之一。

圖10 激波曲面的離散Fig.10 Discretization of the shock surface

圖11 激波曲率的估計Fig.11 Estimation of the shock curvature

離散激波流場有解等效為流場中所有離散單元有解。下面介紹流場內(nèi)任一點的基本單元能夠成功建立的必要條件。圖12給出了展向相鄰兩網(wǎng)格單元及數(shù)值基本單元的結(jié)構(gòu)與理想情況下典型角度間的關(guān)系。式（114）和式（115）保證了交點 S2的存在性，顯式約束了線段 S1A1的方向，同時隱式約束了其長度，其中式（114）表示點 S1位于點 A1發(fā)出的逆馬赫錐內(nèi)，式（115）表示流線 S1S2與馬赫錐的交點 S2位于當前網(wǎng)格層的下游。為了完成基本單元控制方程的構(gòu)造，還需要保證馬赫線 S2A4及 S2A2與初值線 D1D3的交點存在（參考圖12），該條件難以顯式給出。但當四邊形D1S1A1D4和S1D3D6A1為激波網(wǎng)格單元時，可以根據(jù)圖12中建立的理想化角度關(guān)系對激波初值面上網(wǎng)格長寬比的理想值給出估計。假設(shè)點 S1和點 A1的物理參數(shù)一致，則流線 S1S2的長度 δsl、線段 S1A1的長度 δl和線段 S1A4的長度 δw間滿足式（116）和（117）的關(guān)系，合并得式（118），其中 ta n(β-θ)由式（119）計算。為了方便起見將稱為理想長寬比。對于激波單元S1D3D6A1來說，當激波網(wǎng)格長寬比小于理想長寬比時，交點A4位于激波單元的橫邊 S1D3之內(nèi)，即基本單元的求解被限制在展向相鄰兩激波單元內(nèi)，因此該參數(shù)為激波面網(wǎng)格的劃分方案提供了直觀的參考。圖13給出了典型馬赫數(shù)下激波網(wǎng)格長寬比隨激波角的變化，其中虛線為激波角取馬赫角時的極限值。隨著激波角的增大，理想長寬比逐漸降低，即當激波單元橫向?qū)挾炔蛔儠r，縱向長度需要不斷降低來保證數(shù)值解存在。

圖12 縱向步長的預(yù)估方法Fig.12 Prediction method of the lengthwise step

圖13 理想長寬比與激波角的關(guān)系Fig.13 Relationship between the ideal length-to-width ratio and the shock angle

3.2 網(wǎng)格面的推進及計算

圖14所示為網(wǎng)格面的推進原理，具體流程如圖15所示。圖14中的坐標系標注了圖15中i（縱向）、j（展向）和k（層推進方向）循環(huán)指標對應(yīng)的方向。以3 × 4個網(wǎng)格點的激波面網(wǎng)格為例，此時Nlengthwise=4、Nspanwise=3、mmax為每一層求解時的最

圖14 網(wǎng)格面推進的原理Fig.14 Principle of the mesh plane advancing

圖15 網(wǎng)格面推進流程Fig.15 Flow chart of the mesh plane advancing

大迭代數(shù)。對于“ ω 01”方案，由于不涉及前一層網(wǎng)格的速度張量及壓力梯度的值，因此mmax=1。對于另外兩種格式，則需要通過當前層的物理參數(shù)校正前一層的速度張量及壓力梯度，本文給定mmax=30。網(wǎng)格推進的流程為按照i（縱向）、j（展向）和k（層推進方向）三個循環(huán)指標反復(fù)求解基本單元。計算流程中的殘值定義為式（120），其中i、j和k指標遍歷當前網(wǎng)格面的所有網(wǎng)格點。由于沿展向求解基本單元時，單元間相互獨立，因此該部分循環(huán)可以進行OpenMP并行加速。

4 算例驗證及結(jié)果討論

本文選取圓錐激波依賴區(qū)設(shè)計問題、橢圓錐激波依賴區(qū)設(shè)計問題、小攻角圓錐激波依賴區(qū)設(shè)計問題和基于Bezier曲面的三維曲面激波依賴區(qū)設(shè)計問題4個算例來分別討論和驗證以下內(nèi)容：三維設(shè)計方法的計算誤差和并行效率、對由橫向壓力梯度和攻角引起的三維流動效應(yīng)的模擬能力、當前算法對一般性曲面激波反問題的求解能力。

4.1 圓錐激波依賴區(qū)設(shè)計驗證

算例為半錐角15°的圓錐激波流場的設(shè)計（如圖16所示）。坐標系原點位于左側(cè)半徑為100 mm的圓截面圓心處，激波錐長度為1 000 mm，來流條件見表2。該算例用于驗證當前設(shè)計方法對三維激波依賴區(qū)的設(shè)計能力，討論不同求解策略設(shè)計結(jié)果的計算誤差，并在此基礎(chǔ)上討論當前并行策略的計算效率。

圖16 圓錐激波計算模型Fig.16 Computational model of the conical shock

表2 來流參數(shù)Table 2 Inflow parameters

圖17所示為三維特征線所用的激波面網(wǎng)格（紅色網(wǎng)格，縱向 × 展向 = 100 × 50）和計算得到的流場網(wǎng)格。其中，激波面為位于yz坐標系第四象限的1/4圓錐，藍色部分為反設(shè)計得到的壁面網(wǎng)格，流場網(wǎng)格量約為24.3萬，兩個對稱面邊界分別為y= 0 mm和z=0 mm平面。圖18所示為“ ω 01-TL”方案計算得到的流場，實體部分為計算得到的流場壓力分布，切面由鏡像后的完整流場中截取而來。由圖可知，當前方法能夠得到光滑的流場設(shè)計結(jié)果。

圖17 激波面網(wǎng)格及設(shè)計所得的流場網(wǎng)格Fig.17 Shock-surface mesh and the designed flow-field mesh

圖18 設(shè)計所得流場的壓力分布Fig.18 Pressure distribution of the designed flow field

圖19所示為四種計算策略計算得到的壁面壓力及馬赫數(shù)分布。四種計算策略由不同時間積分、調(diào)節(jié)指數(shù)和線性方程組求解方法組合得到（如表3所示，其中“FP”和“TL”為該表右側(cè)列出的線性方程組求解方法的簡寫）。由于列主元消去法無法穩(wěn)定求解當前流場，因此選擇全主元消去法（采用開源矩陣庫Eigen中的fullPivLu函數(shù)進行求解）作為對比項。

表3 計算策略Table 3 Calculation methods

圖19中實線部分的激波面網(wǎng)格量為100 × 50，虛線部分的激波面網(wǎng)格量為50 × 50。黑色實線為激波線上取100個點時的二維特征線方法的計算結(jié)果。如果認為與二維特征線的計算結(jié)果越近則計算精度越高，則可以得到以下結(jié)論：1）四種方法的計算精度排序由高到低依次為：“ ω 12-TL” ≈“ ω -V-TL”＞“ ω 01-TL” ＞“ ω 01-FP” ；2）時間積分方案“ ω 12”與“ ω 01”對當前算例的計算結(jié)果幾乎沒有區(qū)別；3）使用Tikhonov-Lagrange擬合法求解線性方程組相比fullPivLu方法能夠有效地提高計算精度，并且改變網(wǎng)格量對計算結(jié)果的影響較小，而fullPivLu方法則對網(wǎng)格量極為敏感，使用50 × 50的粗網(wǎng)格和100 × 50的精細網(wǎng)格的計算結(jié)果完全不同，且粗網(wǎng)格計算時反倒得到較好的結(jié)果；4）除了“ ω 01-FP”方案，其余三種方案計算得到的壁面尾部馬赫數(shù)誤差低于0.16%、壓力誤差低于0.17%，因此能夠滿足一般流場的設(shè)計需求。

“ ω 12-TL”和“ ω -V-TL”方案考慮了壓力梯度和速度張量沿特征線的變化，因而獲得了更為精確的結(jié)果，但當求解前一層網(wǎng)格面（圖20中第i-1層）尾部格點的壓力梯度和速度張量時，由于下游點6和點7的缺失，導(dǎo)致該點處實際所用的模板與周圍點存在較大差異，難以得到光滑的最小二乘解，導(dǎo)致了圖19中計算結(jié)果在尾部的振蕩，從而引發(fā)穩(wěn)定性問題，因此改良前一層網(wǎng)格的壓力梯度及速度張量的計算方法是后續(xù)研究的重點之一?；谝陨涎芯拷Y(jié)果，選擇“ ω 01-TL”作為后續(xù)激波依賴區(qū)的設(shè)計方法。

圖19 不同計算方法設(shè)計結(jié)果的對比Fig.19 Comparison of the results from different calculation methods

圖20 網(wǎng)格面尾部點的最小二乘計算模板Fig.20 Calculation stencil of the least square method for the ending points on the mesh surface

選擇網(wǎng)格量為100 × 50的激波面網(wǎng)格來研究本文并行策略的計算效率，計算格式為“ ω 01-TL”，測試所用CPU為Intel i7-9750H，對比結(jié)果見表4。計算結(jié)果表明，本文提出的設(shè)計方法具有較高的并行效率。

表4 并行計算加速比和計算效率Table 4 Speedup and efficiency of the parallel computation

4.2 橢圓錐激波依賴區(qū)設(shè)計驗證

算例為橢圓錐激波依賴區(qū)的設(shè)計（圖21所示）。橢圓半錐角分別為15°和17°，坐標系原點位于左側(cè)橢圓截面的中心，且距橢圓錐頂點200 mm，激波錐的兩截面間距為1 000 mm，來流條件見表2。圖22為所用的激波面網(wǎng)格（紅色網(wǎng)格，縱向×展向 = 100 × 50），兩個對稱面邊界分別為y= 0 mm和z= 0 mm平面，右側(cè)為計算得到的流場網(wǎng)格，藍色部分為反設(shè)計得到的壁面形狀。流場的壓力分布如圖23所示。本算例用于驗證當前設(shè)計方法對由橫向壓力梯度引起的流場三維效應(yīng)的模擬能力。

圖21 橢圓錐激波計算模型Fig.21 Computational model of the elliptical-cone shock

圖22 激波面網(wǎng)格及設(shè)計所得的流場網(wǎng)格Fig.22 Shock-surface mesh and the designed flow-field mesh

圖23 設(shè)計所得流場的壓力分布Fig.23 Pressure distribution of the designed flow field

為了驗證本文設(shè)計方法的準確性，采用商業(yè)軟件Fluent對反設(shè)計得到的壁面在表2來流條件下進行無黏數(shù)值模擬，計算格式為AUSM+，網(wǎng)格如圖24所示，網(wǎng)格量約為136萬。圖25所示為沿x軸x= 240 mm和x= 500 mm切面的無量綱壓力云圖對比，其中CFD表示Fluent給出的數(shù)值模擬結(jié)果，3D-MOC為本文的設(shè)計方法，可以看到兩種方法計算得到的截面壓力分布基本一致，激波輪廓在相接處吻合良好。圖26給出了兩切面對應(yīng)的壁面壓力分布和馬赫數(shù)分布，兩種方法計算得到的壓力分布曲線幾乎重合，最大誤差位于z= 0 mm（橢圓截面短軸）位置處。相比而言，兩種方法計算得到的馬赫數(shù)有較大的誤差，且最大誤差同樣位于z= 0 mm處。表5給出了z=0 mm處具體計算數(shù)值的對比，可以看到兩截面處的最大壓力誤差不超過0.3%，最大馬赫數(shù)誤差不超過1%，能夠滿足一般的飛行器設(shè)計需求。

表5 計算誤差對比Table 5 Comparison of the calculation errors

圖24 CFD計算網(wǎng)格Fig.24 Computational grid of CFD

圖25 軸向兩切面壓力分布對比Fig.25 Comparison of pressure contours in two slices along the axial direction

圖26 軸向兩切面壓力及馬赫數(shù)分布對比Fig.26 Comparison of wall pressure and Mach number distributions in two slices along the axial direction

4.3 有攻角圓錐激波依賴區(qū)設(shè)計驗證

算例為攻角α= 3°、半錐角15°的圓錐激波依賴區(qū)的設(shè)計（如圖16所示），幾何參數(shù)與4.1節(jié)一致，其余來流條件與表2一致。圖27為所用的激波面網(wǎng)格（紅色網(wǎng)格，縱向 × 展向 = 100 × 100），對稱面邊界為z= 0 mm平面，右側(cè)為計算得到的流場網(wǎng)格，網(wǎng)格量約48.5萬，藍色部分為反設(shè)計得到的壁面形狀。流場的壓力分布如圖28所示。本算例用于驗證當前設(shè)計方法對由攻角引起的流場三維效應(yīng)的模擬能力。

圖27 激波面網(wǎng)格及設(shè)計所得的流場網(wǎng)格Fig.27 Shock-surface mesh and the designed flow-field mesh

圖28 設(shè)計所得流場的壓力分布Fig.28 Pressure distribution of the designed flow field

采用如圖29所示的網(wǎng)格對設(shè)計得到的壁面流場進行數(shù)值模擬，網(wǎng)格量約為145萬。圖30所示為沿x軸x= 120 mm和x= 380 mm切面的無量綱壓力云圖的對比?？梢钥吹?，兩種方法計算得到的截面壓力分布基本一致，激波位置在相接處吻合良好。圖31給出了兩切面對應(yīng)的壁面壓力分布和馬赫數(shù)分布，兩種方法計算得到的壓力分布曲線幾乎重合，但馬赫數(shù)分布有一定差異。壁面壓力和馬赫數(shù)的最大誤差均位于y= 0 mm位置處。表6給出了y= 0 mm處具體計算數(shù)值的對比，可以看到，兩截面處的最大壓力誤差約0.3%，最大馬赫數(shù)誤差不超過1.7%，能夠滿足一般的飛行器設(shè)計需求。

表6 計算誤差對比Table 6 Comparison of the calculation errors

圖29 CFD計算網(wǎng)格Fig.29 Computational grid of CFD

圖30 軸向兩切面壓力分布對比Fig.30 Comparison of pressure contours in two slices along the axial direction

圖31 軸向兩切面壓力及馬赫數(shù)分布對比Fig.31 Comparison of wall pressure and Mach number distributions in two slices along the axial direction

4.4 一般性三維曲面激波依賴區(qū)設(shè)計驗證

算例為一般性曲面激波依賴區(qū)的設(shè)計驗證，采用3 × 1次Bezier曲面構(gòu)造激波形狀，式（121）為Bezier曲面的控制點，來流條件見表2。圖32為所用的激波面網(wǎng)格（紅色網(wǎng)格，縱向×展向 = 100 × 50），兩個對稱面邊界分別為y= 0 mm和z= 0 mm平面，右側(cè)為計算得到的流場網(wǎng)格，藍色部分為反設(shè)計得到的壁面形狀。流場的壓力分布如圖33所示。

圖32 激波面網(wǎng)格及設(shè)計所得的流場網(wǎng)格Fig.32 Shock-surface mesh and the designed flow-field mesh

圖33 設(shè)計所得流場的壓力分布Fig.33 Pressure distribution of the designed flow field

采用圖34所示的網(wǎng)格對設(shè)計得到的壁面流場進行數(shù)值模擬，網(wǎng)格量約為229萬。圖35所示為沿x軸x= 200 mm和x= 450 mm切面的無量綱壓力云圖對比，可以看到兩種方法計算得到的截面壓力分布基本一致，激波輪廓在相接處吻合良好。圖36給出了兩切面對應(yīng)的壁面壓力分布和馬赫數(shù)分布。兩種方法計算得到的壓力分布曲線幾乎重合，x= 200 mm和x= 450 mm切面的最大壓力和馬赫數(shù)誤差分別近似位于z= 127 mm和z= 147 mm位置處。表7給出了最大誤差處具體計算數(shù)值的對比?？梢钥吹?，兩截面處的最大壓力誤差不超過0.2%，最大馬赫數(shù)誤差不超過1.3%，能夠滿足一般的飛行器設(shè)計需求。

圖34 CFD計算網(wǎng)格Fig.34 Computational grid of CFD

圖36 軸向兩切面壓力及馬赫數(shù)分布對比Fig.36 Comparison of wall pressure and Mach number distributions in two slices along the axial direction

表7 計算誤差對比Table 7 Comparison of the calculation errors

5 結(jié) 論

本文基于三維特征線理論，提出了一種針對三維曲面激波流場的反設(shè)計方法。首先基于特征分解重構(gòu)了三維流場的控制方程，針對控制方程的離散形式發(fā)展了三種不同的時間積分策略，并提出了能夠穩(wěn)定求解基本單元線性方程組的Tikhonov-Lagrange擬合法。通過與二維特征線方法及數(shù)值模擬結(jié)果對比，分析了當前設(shè)計方法對典型三維激波流場的計算誤差及并行效率。結(jié)論如下：

1）利用圓錐激波流場算例對比了不同計算方案對設(shè)計結(jié)果的影響。一方面，“ ω 12”和“ ω -V”時間積分方案給出了基本一致的壁面壓力和馬赫數(shù)分布，且其誤差約為“ ω 01”方案的一半，但兩種格式在每個推進層的末尾由于最小二乘計算模板的變化容易產(chǎn)生數(shù)值振蕩。另一方面，針對線性方程組的求解，本文提出的Tikhonov-Lagrange擬合法在網(wǎng)格加密前后能給出規(guī)律基本一致的壁面參數(shù)分布，而標準的列主元消去法無法穩(wěn)定計算流場，全主元消去法對于激波面的網(wǎng)格量較為敏感，無法給出一致的計算結(jié)果，因此本文提出的Tikhonov-Lagrange擬合法在穩(wěn)定求解控制方程組方面具有較大的優(yōu)勢。

2）采用本文提出的“ ω 01-TL”方案研究了當前設(shè)計方法的并行計算效率。針對激波面網(wǎng)格為100×50的圓錐激波流場算例，當并行核數(shù)為6時計算時間為 9 min 53 s，加速比為4.24，并行效率為70.7%，因此本文提出的設(shè)計方法具有較高的并行計算效率。

3）針對橢圓錐激波流場算例、小攻角圓錐激波流場算例及由Bezier曲面描述的一般性曲面激波流場算例，驗證了本文設(shè)計方法對三維流動效應(yīng)的模擬能力。與無黏數(shù)值模擬結(jié)果相比，橢圓錐激波流場在x= 240 mm和500 mm切面處的最大壓力誤差不超過0.3%，最大馬赫數(shù)誤差不超過1%；小攻角圓錐激波流場在x= 120 mm和380 mm切面處的最大壓力誤差約0.3%，最大馬赫數(shù)誤差不超過1.7%；Bezier曲面激波流場在x= 200 mm和450 mm切面處的最大壓力誤差約0.2%，最大馬赫數(shù)誤差不超過1.3%。因此，本文提出的設(shè)計方法能夠針對典型三維曲面激波流場的反問題給出合理的設(shè)計結(jié)果。

附錄A