基于三維特征線理論的曲面激波流場反設計方法

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

(1.非定常空氣動力學與流動控制工業和信息化部重點實驗室，南京 210016；2.南京航空航天大學 航空學院，南京 210016)

0 引 言

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

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

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

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

1 控制方程

依據三維特征線理論[37]，馬赫錐及流線上的特征速度滿足圖1所示的幾何關系，其中馬赫錐由點S2發出，流線 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 速度關系與馬赫錐Fig.1 Velocity relation and Mach cone

2 基本單元的構造及求解方法

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

2.1 基本單元的構型及相關計算幾何算法

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

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

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

假設線段兩端點分別為點1和點2，插值點3位于點1和點2之間，則可用式（19～20）所示的線性插值公式計算點3的物理參數，其中d31為點3到點1的距離，d21為點2到點1的距離。在計算特征線平均物理參數時，插值系數取 ξ =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等黑斜體英文字母均表示坐標系原點到對應點的距離矢量。

圖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物理參數的初值為點 A3和點 A1物理參數的代數平均（式19， ξ =0.5），利用式（26）計算點 S1的坐標，并由式（27）更新插值系數 ξ，最后利用線性插值對點 S1的坐標及物理參數進行更新。通過迭代上述步驟，直到點 S1的坐標變化小于容差 ε ，本文所有算法所用的容差統一取 ε =1.0×10-6。

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

2.1.5 對稱面邊界相關算法

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

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

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

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

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

2.1.6 Cone_X_Curve算法

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

2.1.7 基本單元的求解流程

基于以上6個基本算法，基本單元求解的流程可以表述為：1）用點 S1的物理參數給點 S2賦初值，該兩點為流線的兩端點；2）利用Ray_X_Ray算法計算流線 S1S2和馬赫線 A1S2交點 S2的坐標；3） A3的物理參數的初值采用點 A1和點 S1的平均值，然后利用Ray_X_Line算法計算馬赫線 S2A3和線段 S1A1的交點A3的坐標，并利用線性插值計算點 A3的物理參數；4）利用Cone_X_Curve算法插值求解點 A2和點 A4的坐標及物理參數；5）求解線性方程組，如果該點位于對稱面邊界上，則需要利用對稱面邊界對求解結果進行校正；6）重復步驟b～e直到待求點的空間坐標、壓力和馬赫數在兩次迭代間的變化小于給定的容差 ε，或總迭代數大于預設值imax。最終構成的線性方程組由17個方程構成，如表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轉化為差分項可得式（43～46）；沿流線由任一點1到點2（待求點）對微分方程（14～18）進行時間積分，并將特征算子 ?s轉化為差分項可得式（47～51）。

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

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

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

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

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

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

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

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

當b2≠0時，

當b2=0時：

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

2.2.2 線性方程組的求解

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

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

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

3 三維曲面激波依賴區設計方法

本節包含曲面激波的離散及網格面的推進算法兩部分內容。

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

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

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

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

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

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

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

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

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

3.2 網格面的推進及計算

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

圖14 網格面推進的原理Fig.14 Principle of the mesh plane advancing

圖15 網格面推進流程Fig.15 Flow chart of the mesh plane advancing

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

4 算例驗證及結果討論

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

4.1 圓錐激波依賴區設計驗證

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

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

表2 來流參數Table 2 Inflow parameters

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

圖17 激波面網格及設計所得的流場網格Fig.17 Shock-surface mesh and the designed flow-field mesh

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

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

表3 計算策略Table 3 Calculation methods

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

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

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

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

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

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

4.2 橢圓錐激波依賴區設計驗證

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

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

圖22 激波面網格及設計所得的流場網格Fig.22 Shock-surface mesh and the designed flow-field mesh

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

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

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

圖24 CFD計算網格Fig.24 Computational grid of CFD

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

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

4.3 有攻角圓錐激波依賴區設計驗證

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

圖27 激波面網格及設計所得的流場網格Fig.27 Shock-surface mesh and the designed flow-field mesh

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

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

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

圖29 CFD計算網格Fig.29 Computational grid of CFD

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

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

4.4 一般性三維曲面激波依賴區設計驗證

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

圖32 激波面網格及設計所得的流場網格Fig.32 Shock-surface mesh and the designed flow-field mesh

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

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

圖34 CFD計算網格Fig.34 Computational grid of CFD

圖36 軸向兩切面壓力及馬赫數分布對比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 結 論

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

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

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

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

附錄A