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

(1.非定?？諝鈩?dòng)力學(xué)與流動(dòng)控制工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，南京 210016；2.南京航空航天大學(xué) 航空學(xué)院，南京 210016)

0 引 言

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

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

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

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

1 控制方程

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

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

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

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

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

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

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

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

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

2.1.2 特征線速度矢量及法矢的計(jì)算方法

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

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

2.1.3 Ray_X_Ray算法

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

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

2.1.4 Ray_X_Line算法

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

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

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

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

2.1.5.1 對(duì)稱面鏡像點(diǎn)坐標(biāo)及物理參數(shù)的計(jì)算

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

2.1.5.2 對(duì)稱面上點(diǎn)的坐標(biāo)及物理參數(shù)的校正

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

2.1.6 Cone_X_Curve算法

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

圖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 初值線鏡像點(diǎn)的構(gòu)造Fig.8 Construction of the mirror points of the initial value curve

2.1.7 基本單元的求解流程

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

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

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

2.2.1 控制方程的離散

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

式（43～51）中： δt表示變量t的一個(gè)增長(zhǎng)量；變量上方的單橫線表示該段特征線上的參數(shù)平均值（式19）。將式（43～49，51）中的時(shí)間積分項(xiàng)處理為（52～59）。其中，馬赫線和流線的時(shí)間項(xiàng)分別按式（60～61）計(jì)算；δl為特征線的長(zhǎng)度，且定義特征線上的速度矢量方向?yàn)橛沙踔迭c(diǎn)指向待求點(diǎn)時(shí) δl為正值，否則為負(fù)值。離散后總的控制方程見(jiàn)附錄。

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

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

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

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

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

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

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

當(dāng)b2≠0時(shí)，

當(dāng)b2=0時(shí)：

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

2.2.2 線性方程組的求解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 圓錐激波依賴區(qū)設(shè)計(jì)驗(yàn)證

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

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

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

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

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

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

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

表3 計(jì)算策略Table 3 Calculation methods

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

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

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

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

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

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

4.2 橢圓錐激波依賴區(qū)設(shè)計(jì)驗(yàn)證

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

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

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

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

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

表5 計(jì)算誤差對(duì)比Table 5 Comparison of the calculation errors

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

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

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

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

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

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

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

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

表6 計(jì)算誤差對(duì)比Table 6 Comparison of the calculation errors

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

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

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

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

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

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

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

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

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

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

表7 計(jì)算誤差對(duì)比Table 7 Comparison of the calculation errors

5 結(jié) 論

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

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

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

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

附錄A