劉超宇，屈峰，*，孫迪，劉傳振，錢戰(zhàn)森，白俊強(qiáng)

1．西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

2．中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074

3．中國(guó)航空工業(yè)空氣動(dòng)力研究院，沈陽(yáng) 110034

乘波體作為一種典型的高超聲速飛行器構(gòu)型，能夠通過(guò)附著于前緣的激波將高壓氣流限制在飛行器下表面，在保證該區(qū)域流動(dòng)均勻的同時(shí)可有效阻止溢流，從而實(shí)現(xiàn)較高的升阻比并便于實(shí)現(xiàn)機(jī)體與發(fā)動(dòng)機(jī)的一體化設(shè)計(jì)［1-2］。因此，乘波體的研究受到了日益廣泛的關(guān)注。傳統(tǒng)的乘波體，如楔形流乘波體［3］、錐導(dǎo)乘波體［4］及密切錐乘波體［5］等，均是從已有的激波形狀或流場(chǎng)出發(fā)，通過(guò)追蹤流線生成乘波面。劉傳振等［6-9］從密切錐方法出發(fā)，發(fā)展了定平面乘波體設(shè)計(jì)方法，提高了傳統(tǒng)乘波體的設(shè)計(jì)靈活性和整體升阻特性。但是該類乘波體的設(shè)計(jì)方法是基于二維無(wú)黏理論推導(dǎo)出來(lái)的，因?yàn)楹雎匀S以及黏性效應(yīng)，其在設(shè)計(jì)工況下仍會(huì)出現(xiàn)溢流，升阻比難以達(dá)到最優(yōu)；另外，該類乘波體仍具有傳統(tǒng)乘波體在偏離設(shè)計(jì)條件下，氣動(dòng)特性會(huì)出現(xiàn)惡化的不足［10］；同時(shí)，工程化的外形需要在理論模型的基礎(chǔ)上實(shí)現(xiàn)頭部/前緣鈍化以及側(cè)緣設(shè)計(jì)，而這些局部變化會(huì)對(duì)其氣動(dòng)特性產(chǎn)生較大的不利影響。因此，為使定平面形狀的密切錐乘波體更具工程應(yīng)用價(jià)值，有必要在考慮黏性的情況下，對(duì)其開(kāi)展精細(xì)的氣動(dòng)優(yōu)化設(shè)計(jì)［11］。

目前，國(guó)內(nèi)外學(xué)者針對(duì)高超聲速乘波體的優(yōu)化 設(shè)計(jì)已開(kāi)展了大量研究。例如：Rodi［12-13］采 用遺傳算法和粒子群算法對(duì)相切流場(chǎng)乘波體的貝塞爾曲線前緣進(jìn)行了優(yōu)化設(shè)計(jì)研究；國(guó)內(nèi)的張鋒濤等［14］采用神經(jīng)網(wǎng)絡(luò)對(duì)乘波體進(jìn)行優(yōu)化，該方法在保證計(jì)算穩(wěn)定的前提下有效提高了優(yōu)化效率；徐佳勝［15］和吳功名［16］基于Kriging 代理模型方法分別對(duì)密切錐乘波體和類HTV-2 飛行器進(jìn)行優(yōu)化設(shè)計(jì)。但是，上述工作都是基于智能優(yōu)化算法或代理模型方法進(jìn)行的。這些優(yōu)化方法在針對(duì)設(shè)計(jì)變量眾多的乘波體三維整機(jī)開(kāi)展氣動(dòng)優(yōu)化設(shè)計(jì)時(shí)，會(huì)出現(xiàn)計(jì)算量過(guò)大的不足，一些情況下甚至可能得不到收斂結(jié)果［17-19］。與上述優(yōu)化方法相比，基于伴隨方程的梯度類優(yōu)化設(shè)計(jì)方法可實(shí)現(xiàn)計(jì)算量與設(shè)計(jì)變量之間的基本解耦，并且精度較高，在大規(guī)模設(shè)計(jì)變量問(wèn)題中優(yōu)勢(shì)明顯［20］，比較適合于設(shè)計(jì)變量眾多的高超聲速乘波體整機(jī)氣動(dòng)外形優(yōu)化設(shè)計(jì)。

一般來(lái)說(shuō)，伴隨類優(yōu)化方法可分為連續(xù)伴隨優(yōu)化方法和離散伴隨優(yōu)化方法。而這2 類方法在高超聲速氣動(dòng)優(yōu)化設(shè)計(jì)中都得到了廣泛的應(yīng)用。連續(xù)伴隨優(yōu)化方法的研究方面，Kline 等［21］采用連續(xù)伴隨方法針對(duì)高超聲速進(jìn)氣道開(kāi)展了優(yōu)化研究；高昌等［22-23］通過(guò)建立基于連續(xù)伴隨的優(yōu)化方法，針對(duì)二維高超聲速進(jìn)氣道和Sanger 飛行器機(jī)翼開(kāi)展了優(yōu)化設(shè)計(jì)研究。相較于連續(xù)伴隨優(yōu)化方法，離散伴隨方程的優(yōu)勢(shì)為它與Navier-Stokes 方程的導(dǎo)數(shù)關(guān)系更清晰，實(shí)現(xiàn)起來(lái)比較方便且獲取的梯度信息更為精確［24］。因此，諸多學(xué)者針對(duì)離散伴隨方法在高超聲速優(yōu)化設(shè)計(jì)問(wèn)題中的應(yīng)用開(kāi)展了研究。例如，Damm 等［25］基于離散伴隨方法針對(duì)NASA P2 高超聲速進(jìn)氣道開(kāi)展了氣動(dòng)優(yōu)化設(shè)計(jì)。高正紅等［26］基于離散伴隨方程法對(duì)某高超聲速導(dǎo)彈的前體進(jìn)行了氣動(dòng)外形優(yōu)化。宋紅超等［27］發(fā)展了基于離散伴隨方法的高超聲速單邊膨脹噴管優(yōu)化設(shè)計(jì)方法。但是，上述將伴隨方法應(yīng)用于高超聲速優(yōu)化問(wèn)題中的工作多是針對(duì)飛行器局部簡(jiǎn)單部件開(kāi)展優(yōu)化設(shè)計(jì)的，缺乏在復(fù)雜飛行器整機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)當(dāng)中的應(yīng)用。另一方面，離散伴隨優(yōu)化多點(diǎn)設(shè)計(jì)的方式一般為加權(quán)平均，設(shè)計(jì)結(jié)果依賴于選取的權(quán)重系數(shù)，且需要進(jìn)行多次嘗試，而結(jié)構(gòu)化網(wǎng)格求解器在相同資源條件下計(jì)算效率及精度更高［24］。因此，有必要實(shí)現(xiàn)基于結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)的高超聲速飛行器離散伴隨氣動(dòng)優(yōu)化設(shè)計(jì)方法。

針對(duì)上述問(wèn)題，本文結(jié)合基于全速域通量求解方法和RANS（Reynolds Averaged Navier-Stokes equation）湍流模型的結(jié)構(gòu)化網(wǎng)格求解器、魯棒的結(jié)構(gòu)網(wǎng)格變形方法、適用于高維設(shè)計(jì)變量的自由變形 （Free Form Deformation， FFD） 參數(shù)化方法、離散伴隨方法與序列二次規(guī)劃 （Sequential Quadratic Programming， SQP） 算法，實(shí)現(xiàn)了基于離散伴隨的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)方法?；谠摲椒?，針對(duì)定平面形狀的密切錐乘波體構(gòu)型分別開(kāi)展了單一設(shè)計(jì)工況和多工況的三維整機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)。

1 基于離散伴隨的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)方法

1.1 高超聲速數(shù)值模擬方法

本文采用自主研發(fā)的高分辨率高超聲速RANS 求解器針對(duì)高超聲速流場(chǎng)開(kāi)展數(shù)值模擬。該求解方法基于三維可壓縮Navier-Stokes方程：

式中：q為流場(chǎng)守恒變量；f(q)表征無(wú)黏通量；fv(q)為黏性通量；Ω代表控制體；?Ω為控制體邊界；ΔSi為相應(yīng)單元界面的面積。關(guān)于該控制方程的詳細(xì)介紹，可參考文獻(xiàn)［28］。

針對(duì)上述控制方程，采用基于格心的有限體積法構(gòu)建了多塊結(jié)構(gòu)網(wǎng)格的流場(chǎng)求解器［29］。時(shí)間推進(jìn)采用隱式LUSGS-GMRES（Lower Upper Symmetric Gauss Seidel-Generalized Minimum Residual）混合方法以保證定常流動(dòng)求解的穩(wěn)定性和收斂效率［30-31］。空間離散方法采用二階MUSCL（Monotonic Upstream-centered Scheme for Conservation Laws）重構(gòu)方法以及minmod 限制器［32］，通量求解方法采用適用于全速域流動(dòng)計(jì)算的HLLEMS（modified HLLE scheme with a switch function）格式［33］。湍流模型采用一方程SA 模型［34］。為了提高求解效率，該求解器采用了MPI（Message Passing Interface）并行和多重網(wǎng)格技術(shù)。

1.2 離散伴隨梯度求解方法

采用離散伴隨方程法求解梯度是本文所實(shí)現(xiàn)的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)平臺(tái)的核心。離散伴隨方程法首先對(duì)非線性化的偏微分流動(dòng)控制方程進(jìn)行離散，然后對(duì)其進(jìn)行線性化，最后通過(guò)離散形式的線化控制方程，推導(dǎo)出離散伴隨方程及其邊界條件。如果伴隨方程的形式和原始流場(chǎng)控制方程的離散形式完全對(duì)應(yīng)，則可以求解出目標(biāo)函數(shù)和控制方程離散形式的精確導(dǎo)數(shù)［35］。相比于連續(xù)伴隨方程法，離散伴隨方程法得到的梯度往往更接近真實(shí)值，可以同時(shí)用于自適應(yīng)網(wǎng)格生成和誤差分析［36］。因此，本文采用離散伴隨方法對(duì)氣動(dòng)目標(biāo)函數(shù)的梯度進(jìn)行求解。

假設(shè)氣動(dòng)優(yōu)化目標(biāo)函數(shù)為F(G(x)，Q(x))，其中：x為設(shè)計(jì)變量，表征FFD 方法中各個(gè)控制點(diǎn)的位移；G(x)為由設(shè)計(jì)變量x所確定的 CFD計(jì)算網(wǎng)格，包括表面網(wǎng)格Gs和空間網(wǎng)格Gv；Q為流場(chǎng)解向量。通過(guò)FFD 外形參數(shù)化，可由設(shè)計(jì)變量x獲得表面網(wǎng)格Gs(x)，再由網(wǎng)格變形法獲得空間網(wǎng)格Gv(x)。在給定一組設(shè)計(jì)變量x的情況下，通過(guò)求解式（2）便可得到相應(yīng)的一組空間流場(chǎng)解Q。

式中：R為流場(chǎng)控制方程的殘差，對(duì)于收斂的定常流場(chǎng)，可認(rèn)為控制方程的殘差近似為0。

將F、R分別對(duì)設(shè)計(jì)變量x求全導(dǎo)數(shù)，則有

式中：dQdx為空間流場(chǎng)解對(duì)設(shè)計(jì)變量的全導(dǎo)數(shù)。若采用傳統(tǒng)的梯度求解方法求解，需對(duì)x的每一個(gè)分量進(jìn)行擾動(dòng)，然后分別求解流場(chǎng)控制方程。在設(shè)計(jì)變量較多時(shí)，計(jì)算量毫無(wú)疑問(wèn)是巨大的。為了提高求解效率，將式（4）作恒等變換為

將式（5）代入式（3）可得

則可以將矩陣求逆轉(zhuǎn)化為求解線性方程組，即轉(zhuǎn)化為求解伴隨變量：

此時(shí)，式（6）可以寫為

式中：?F?G只需要在收斂的流場(chǎng)上對(duì)網(wǎng)格變量求偏導(dǎo)數(shù)即可；dGdx為空間網(wǎng)格對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)，通過(guò)網(wǎng)格變形算法確定；?R?G為流場(chǎng)殘差對(duì)于空間網(wǎng)格的偏導(dǎo)數(shù)，在流場(chǎng)計(jì)算中可求得此項(xiàng)，計(jì)算量也相對(duì)較小。因此準(zhǔn)確高效求解式（9）主要體現(xiàn)在偏導(dǎo)數(shù)矩陣的獲取以及伴隨方程式的求解。偏導(dǎo)數(shù)矩陣的獲取通過(guò)自動(dòng)微分技術(shù)來(lái)實(shí)現(xiàn)，具體可參考文獻(xiàn)［37］；伴隨方程式（8）本質(zhì)上是一個(gè)大規(guī)模、高度稀疏的線性方程組，采用廣義最小殘量（Generalized Minimum Residual， GMRES）算法來(lái)對(duì)其進(jìn)行求解［30］。

1.3 氣動(dòng)優(yōu)化設(shè)計(jì)平臺(tái)

基于1.1 節(jié)所述的高精度RANS 流場(chǎng)求解器，通過(guò)結(jié)合適用于高維的FFD 參數(shù)化方法［38］、魯棒的結(jié)構(gòu)網(wǎng)格變形方法［39-40］、離散伴隨方法和SQP 算法［41］，構(gòu)建了基于離散伴隨的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)平臺(tái)，其流程圖如圖1所示。

圖1 基于離散伴隨的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)方法流程圖Fig. 1 Flow chart of aerodynamic shape optimization design software for hypersonic vehicle based on discretized adjoint method

首先，基于高超聲速飛行器初始外形生成用于CFD 計(jì)算的空間網(wǎng)格，并從空間網(wǎng)格中提取出目標(biāo)表面網(wǎng)格。其次，根據(jù)優(yōu)化的問(wèn)題建立包裹飛行器外形的FFD 控制體，將控制體上的坐標(biāo)點(diǎn)位移作為外形設(shè)計(jì)變量。對(duì)控制體坐標(biāo)點(diǎn)進(jìn)行擾動(dòng)之后，輸出新外形的表面網(wǎng)格。表面網(wǎng)格進(jìn)一步傳遞至網(wǎng)格變形模塊，根據(jù)逆距離權(quán)重插值算法實(shí)現(xiàn)了空間網(wǎng)格的變形更新，生成了新的空間計(jì)算網(wǎng)格。隨后，對(duì)于新的空間網(wǎng)格進(jìn)行高精度CFD 流場(chǎng)計(jì)算，得到收斂的流場(chǎng)解向量以及氣動(dòng)目標(biāo)函數(shù)的值。根據(jù)收斂的流場(chǎng)解向量構(gòu)造伴隨方程，求解伴隨方程獲得目標(biāo)函數(shù)相對(duì)于設(shè)計(jì)變量的梯度。對(duì)于多目標(biāo)優(yōu)化問(wèn)題，采用加權(quán)平均方法將多目標(biāo)優(yōu)化問(wèn)題轉(zhuǎn)化為單目標(biāo)優(yōu)化問(wèn)題。不同權(quán)重系數(shù)的設(shè)置能夠改變優(yōu)化的側(cè)重點(diǎn)，以達(dá)到期望的優(yōu)化結(jié)果。最后，目標(biāo)函數(shù)以及它們對(duì)應(yīng)的梯度一起反饋給優(yōu)化算法，判斷是否滿足優(yōu)化收斂準(zhǔn)則。若不滿足，優(yōu)化算法SQP 將計(jì)算搜索方向和步長(zhǎng)，獲得新的設(shè)計(jì)變量，進(jìn)行下一步迭代，如此循環(huán)直至優(yōu)化迭代收斂。

2 高超聲速定平面密切錐乘波體

選取基于定平面形狀方法設(shè)計(jì)得到的密切錐乘波體作為初始構(gòu)型［6-8］。圖2 為該乘波體全模的幾何外形，圖3 為乘波體半模的幾何尺寸，其中：X軸為流向，Y軸為升力 向，Z軸為展向。該構(gòu)型的平面形狀為小彎頭單后掠，后掠角為75°，并通過(guò)前緣鈍化、上表面容積擴(kuò)充、增加后體等方法將其工程化。通常，前緣鈍化處理有2 種方式：① 根據(jù)鈍化半徑，對(duì)前緣上下表面進(jìn)行修形；② 下表面前緣不變，整體上移上表面直到滿足鈍化半徑［42-43］。本文采用第1 種方式對(duì)乘波體前緣進(jìn)行鈍化處理，鈍化半徑為 2 cm。由于定平面密切錐乘波體在設(shè)計(jì)時(shí)未考慮底阻的影響，本文為研究三維黏性效應(yīng)以及工程化處理對(duì)乘波體升阻特性的影響，將在優(yōu)化設(shè)計(jì)時(shí)基于分部件積分并扣除掉底阻。本文所選取的乘波體初始構(gòu)型的設(shè)計(jì)工況為：馬赫數(shù)Ma=5，迎角α=6°。

圖2 乘波體初始構(gòu)型Fig. 2 Initial model for waverider

圖3 乘波體平面尺寸Fig. 3 Plane size for waverider

采用粗、中、細(xì)3 套網(wǎng)格對(duì)本文所選用的定平面密切錐乘波體進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，網(wǎng)格量分別為212 萬(wàn)、403 萬(wàn)和801 萬(wàn)。表1 給出了在密切錐乘波體設(shè)計(jì)狀態(tài)下的網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果，表中CL、CD分別為升力系數(shù)和阻力系數(shù)；圖4 給出了3 種網(wǎng)格在對(duì)稱面處壓力系數(shù)（Cp）對(duì)比?？梢钥吹?，3 種網(wǎng)格的阻力系數(shù)CD的計(jì)算結(jié)果相差較小，因此，為保證優(yōu)化設(shè)計(jì)過(guò)程中的計(jì)算精度和計(jì)算效率，本文選用中等網(wǎng)格進(jìn)行后面的氣動(dòng)優(yōu)化設(shè)計(jì)。

表1 乘波體網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果Table 1 Grid independent compute results of waverider

圖4 對(duì)稱面處表面壓力系數(shù)對(duì)比Fig. 4 Comparison of surface pressure coefficient on plane of symmetry

3 離散伴隨方法精度驗(yàn)證

在開(kāi)展氣動(dòng)優(yōu)化設(shè)計(jì)之前，采用離散伴隨方法針對(duì)密切錐乘波體進(jìn)行梯度求解，對(duì)伴隨方法求解結(jié)果進(jìn)行計(jì)算精度的驗(yàn)證。圖5 給出了乘波體的FFD 控制框。該控制框采用20 個(gè)控制剖面，每個(gè)控制剖面沿弦向上下表面各有17 個(gè)控制點(diǎn)，共34 個(gè)控制點(diǎn)。每個(gè)控制點(diǎn)只能在Y向移動(dòng)，用來(lái)對(duì)乘波體型面進(jìn)行擾動(dòng)。另外，為了保持前后緣線為直線，每個(gè)控制剖面上下表面沿弦向第一個(gè)和最后一個(gè)點(diǎn)保持不動(dòng)。因此共有600 個(gè)FFD 控制點(diǎn)位移設(shè)計(jì)變量。

圖5 密切錐乘波體FFD 控制框Fig. 5 FFD box of osculating-cone waverider

對(duì)FFD 控制點(diǎn)進(jìn)行擾動(dòng)，并分別用有限差分法和離散伴隨方程法計(jì)算各氣動(dòng)力系數(shù)相對(duì)于設(shè)計(jì)變量的偏導(dǎo)數(shù)，計(jì)算過(guò)程中差分步長(zhǎng)選取1×10-4、1×10-5和1×10-6。采用一組隨機(jī)數(shù)選取15 個(gè)FFD 控制點(diǎn)的計(jì)算結(jié)果，圖6 給出了2 種方法的計(jì)算結(jié)果對(duì)比，其中，Xi為隨機(jī)選取的第i個(gè)設(shè)計(jì)變量，CMZ為俯仰力矩系數(shù)?？梢钥闯鲇邢薏罘址ǎ‵inite Difference method， FD）和伴隨方程法的計(jì)算結(jié)果吻合較好，說(shuō)明伴隨方程法在求解梯度時(shí)具有較高的計(jì)算精度。

圖6 伴隨方程法和有限差分法計(jì)算結(jié)果對(duì)比Fig. 6 Comparison of computed results between adjoint method and finite difference method

4 高超聲速乘波體氣動(dòng)優(yōu)化設(shè)計(jì)

基于第1 節(jié)實(shí)現(xiàn)的高超聲速飛行器離散伴隨氣動(dòng)優(yōu)化設(shè)計(jì)方法，本節(jié)針對(duì)定平面形狀的密切錐乘波體構(gòu)型分別開(kāi)展了設(shè)計(jì)工況下的三維整機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)以及包含非設(shè)計(jì)點(diǎn)的多工況氣動(dòng)優(yōu)化設(shè)計(jì)。

4.1 單點(diǎn)優(yōu)化

單點(diǎn)優(yōu)化設(shè)計(jì)選取初始構(gòu)型的設(shè)計(jì)工況，即Ma=5，高度H=34 km，α=6°。目標(biāo)函數(shù)及設(shè)計(jì)約束為

式中：CL為升力系數(shù)；CD為阻力系數(shù)；CD0為初始阻力系數(shù)；ti0為初始構(gòu)型第i個(gè)站位處的厚度；ti為優(yōu)化后構(gòu)型第i個(gè)站位處的厚度；V為優(yōu)化構(gòu)型體積；V0為初始構(gòu)型體積。該算例的設(shè)計(jì)目標(biāo)為升阻比最大，而設(shè)計(jì)約束為：飛行器阻力值不超過(guò)初始值、厚度不小于初始的95%、容積不降低。

計(jì)算效率方面，本文基于MPI 并行方法在56核CPU 下開(kāi)展整機(jī)單點(diǎn)優(yōu)化設(shè)計(jì)，并僅耗時(shí)約27 h 即完成13 次梯度優(yōu)化迭代以及59 次流場(chǎng)求解，進(jìn)而得到了收斂的優(yōu)化構(gòu)型。相較于初始構(gòu)型，單點(diǎn)優(yōu)化構(gòu)型滿足厚度約束且其容積保持基本不變。

表2 對(duì)比了優(yōu)化前后不同構(gòu)型的升阻力特性，其中：CDp表征壓差阻力系數(shù)；CDf為摩擦阻力系數(shù)；CL CD為升阻比?？梢钥闯?，本文優(yōu)化得到的構(gòu)型可以顯著改善飛行器的升阻特性，優(yōu)化后飛行器升阻比提升了4.9%。

表2 優(yōu)化前后升阻力特性對(duì)比（單點(diǎn)優(yōu)化）Table 2 Comparison of lift and drag characteristics before and after optimization （single point optimization）

表3 給出了單點(diǎn)優(yōu)化前后壓心xcp位置的變化情況。其中：xcp為縱向壓心位置，即壓心距頭部的距離占全機(jī)長(zhǎng)度的比例?？梢钥闯?，單點(diǎn)優(yōu)化構(gòu)型的壓心位置無(wú)明顯變化。

表3 優(yōu)化前后壓心位置對(duì)比（單點(diǎn)優(yōu)化）Table 3 Comparison of pressure center positions before and after optimization （single point optimization）

圖7 和圖8 給出了初始外形和優(yōu)化后外形的空間流場(chǎng)及壓力系數(shù)云圖對(duì)比。圖9 為選取的乘波體5 個(gè)站位示意圖。圖10 給出了乘波體在不同站位處的外形變化和表面壓力系數(shù)分布對(duì)比。可以看出，初始外形在頭部隆起，使得氣流在背風(fēng)面的流動(dòng)會(huì)受到壓縮進(jìn)而產(chǎn)生激波，激波后的頭部背風(fēng)面處會(huì)有局部高壓區(qū)。而優(yōu)化后的外形在頭部隆起的部位變得更加平緩，使得背風(fēng)面激波強(qiáng)度減弱，激波相對(duì)上移，波后壓力減小。初始外形迎風(fēng)面激波依附在前緣，流場(chǎng)較均勻。優(yōu)化后外形在頭部下表面（靠近對(duì)稱面處）下凸，物面角增大，激波強(qiáng)度更強(qiáng)，在波后出現(xiàn)局部高壓區(qū)，使得頭部處升力增大。因初始外形前緣進(jìn)行鈍化處理，在進(jìn)行黏性計(jì)算后前緣處會(huì)出現(xiàn)溢流，迎風(fēng)面高壓氣流轉(zhuǎn)移到背風(fēng)面，使得背風(fēng)面前緣處的低壓區(qū)減小，導(dǎo)致氣動(dòng)性能的損失。優(yōu)化后的外形則對(duì)背風(fēng)面前緣附近的溢流有所抑制。

圖7 優(yōu)化前后空間流場(chǎng)對(duì)比（單點(diǎn)優(yōu)化）Fig. 7 Comparison of flow field of waverider before and after optimization（single point optimization）

圖8 優(yōu)化前后壓力系數(shù)云圖對(duì)比（單點(diǎn)優(yōu)化）Fig. 8 Comparison of pressure coefficient contours before and after optimization （single point optimization）

圖9 乘波體5 個(gè)站位Fig. 9 Five stations of waverider configuration

圖10 優(yōu)化前后各截面外形和表面壓力系數(shù)分布對(duì)比（單點(diǎn)優(yōu)化）Fig. 10 Comparison of sectional shapes and surface pressure coefficient distributions before and after optimization（single point optimization）

在乘波體的中后部，優(yōu)化前后外形在背風(fēng)面壓力變化較小，優(yōu)化后外形在類似于翼身過(guò)渡處略微上凸，相應(yīng)處的壓力也會(huì)較大。而優(yōu)化后外形的下表面相對(duì)上移，氣流流經(jīng)物面的物面角減小，激波相對(duì)下移，且激波強(qiáng)度減弱，導(dǎo)致波后的壓力也會(huì)更小。由于乘波體中后部的面積較大，這也使得優(yōu)化后外形整體的升力有所下降。但是，由于激波強(qiáng)度的減弱，波阻也會(huì)更小，從而使得優(yōu)化后外形阻力明顯下降。另外，優(yōu)化后外形在前緣處相對(duì)下偏，表面壓力波動(dòng)較大，在一定程度上限制了迎風(fēng)面附近的高壓氣流。

4.2 多點(diǎn)優(yōu)化

本節(jié)選取如下?tīng)顟B(tài)開(kāi)展多點(diǎn)氣動(dòng)優(yōu)化設(shè)計(jì)：Ma=5，H=34 km，αn(n=1，2，3)。其中，非設(shè)計(jì)狀態(tài)為α1=0°和α2=2°，設(shè)計(jì)狀態(tài)為α3=6°。目標(biāo)函數(shù)及設(shè)計(jì)約束為

式中：CL0為初始升力系數(shù)。

多點(diǎn)優(yōu)化問(wèn)題的目標(biāo)函數(shù)為各迎角狀態(tài)下升阻比的線性加權(quán)和。其中，權(quán)重系數(shù)按照各迎角狀態(tài)下初始升阻比的比例，分別取為ε1=0.1、ε2=0.4、ε3=0.5。多點(diǎn)優(yōu)化設(shè)計(jì)的約束為：各迎角狀態(tài)升力值不低于初始值、阻力值不超過(guò)初始值、幾何約束為厚度不小于初始的95%、容積不降低。

計(jì)算效率方面，基于MPI 并行方法在56 核CPU 下開(kāi)展整機(jī)多點(diǎn)優(yōu)化設(shè)計(jì)，僅耗時(shí)約60 h 即完成9 次梯度優(yōu)化迭代以及93 次流場(chǎng)求解，進(jìn)而得到了收斂的優(yōu)化構(gòu)型。相較于初始構(gòu)型，多點(diǎn)優(yōu)化構(gòu)型也滿足厚度約束且其容積保持基本不變。圖11 給出了多點(diǎn)優(yōu)化迭代收斂歷史，圖中L/D為升阻比。在優(yōu)化的前20 步主迭代內(nèi)，各迎角狀態(tài)的升力系數(shù)和升阻比提升明顯，非設(shè)計(jì)狀態(tài)阻力系數(shù)有不同程度的下降，在設(shè)計(jì)狀態(tài)6°迎角阻力系數(shù)有小幅提升。20 步之后，各氣動(dòng)參數(shù)變化趨于平緩。

圖11 多點(diǎn)優(yōu)化迭代收斂歷史Fig. 11 Convergence history of multi-point optimization iterations

表4 對(duì)比了多點(diǎn)優(yōu)化前后不同狀態(tài)升阻力特性。可以看出，非設(shè)計(jì)迎角狀態(tài)的升阻比明顯增大。優(yōu)化效果的來(lái)源主要為升力的大幅增加，同時(shí)阻力有所降低；在設(shè)計(jì)狀態(tài)升阻比小幅提升，阻力略有升高。相對(duì)于初始外形，氣動(dòng)優(yōu)化得到的外形可以保證設(shè)計(jì)狀態(tài)升阻比不下降的同時(shí)，顯著提升非設(shè)計(jì)狀態(tài)的整機(jī)升阻特性。

表4 優(yōu)化前后升阻力特性對(duì)比（多點(diǎn)優(yōu)化）Table 4 Comparison of lift and drag characteristics before and after optimization （multi-point optimization）

表5 給出了多點(diǎn)優(yōu)化前后壓心xcp位置的變化情況。其中，CM為俯仰力矩系數(shù)的絕對(duì)值，力矩參考點(diǎn)為機(jī)頭對(duì)應(yīng)的坐標(biāo)?？梢钥闯?，相較于初始構(gòu)型，多點(diǎn)優(yōu)化構(gòu)型在各迎角下的壓心位置均有所前移。在非設(shè)計(jì)0°迎角狀態(tài)，壓心變化幅度較大，原因可能是在0°迎角下，俯仰力矩系數(shù)絕對(duì)值較小，壓心位置對(duì)氣動(dòng)外形變化的敏感度較大。另外，本文選取研究對(duì)象是定平面形狀密切錐乘波體，并未添加平垂尾、V 尾等氣動(dòng)控制部件，因此壓心變化幅度較大。

表5 優(yōu)化前后壓心位置對(duì)比（多點(diǎn)優(yōu)化）Table 5 Comparison of pressure center positions before and after optimization （multi-point optimization）

圖12 展示了優(yōu)化前后外形變化情況，其中，灰白色為初始外形，藍(lán)色為優(yōu)化外形，5 個(gè)典型站位中黑色實(shí)線為初始外形的截面曲線，紅色實(shí)線為優(yōu)化外形的截面曲線。圖13～圖16 給出了初始外形和優(yōu)化后外形在2 個(gè)非設(shè)計(jì)迎角狀態(tài)下的空間流場(chǎng)及壓力系數(shù)云圖對(duì)比。圖17 和圖18 給出了乘波體在2 個(gè)非設(shè)計(jì)迎角狀態(tài)下各站位處的外形變化和表面壓力系數(shù)分布對(duì)比。與單點(diǎn)優(yōu)化類似，優(yōu)化后的外形在頭部隆起的部位變得更加平緩，使得背風(fēng)面激波強(qiáng)度減弱，激波相對(duì)上移，波后壓力減小，這在非設(shè)計(jì)狀態(tài)優(yōu)化前后的流場(chǎng)對(duì)比中可以明顯看出。初始外形下表面激波依附在前緣，流場(chǎng)較均勻。優(yōu)化后外形在頭部下表面（靠近對(duì)稱面處）因型面的變化產(chǎn)生二次激波，在波后出現(xiàn)局部高壓區(qū)，使得頭部處升力增大。但相較于下壓縮面，上表面的型面變化導(dǎo)致優(yōu)化外形對(duì)來(lái)流的壓縮減弱，背風(fēng)面壓力明顯下降，因此其對(duì)非設(shè)計(jì)迎角狀態(tài)下的氣動(dòng)性能產(chǎn)生主要的影響。因初始外形前緣進(jìn)行鈍化處理，在進(jìn)行黏性計(jì)算后前緣處會(huì)出現(xiàn)不同程度的溢流，在非設(shè)計(jì)狀態(tài)處溢流較為明顯，迎風(fēng)面高壓氣流轉(zhuǎn)移到背風(fēng)面，使得背風(fēng)面前緣處的低壓區(qū)減小，導(dǎo)致氣動(dòng)性能的損失。而優(yōu)化后的外形顯著抑制了背風(fēng)面前緣附近的溢流量，特別是在2°迎角時(shí)，高壓氣流能更好地被限制在前緣下表面。

圖12 優(yōu)化前后外形比較Fig. 12 Comparison of waverider shape before and after optimization

圖13 非設(shè)計(jì)狀態(tài)α1=0°優(yōu)化前后空間流場(chǎng)對(duì)比Fig. 13 Comparison of flow field of waverider before and after optimization at non-design state with α1=0°

圖14 非設(shè)計(jì)狀態(tài)α1=0°優(yōu)化前后壓力系數(shù)云圖對(duì)比Fig. 14 Comparison of pressure coefficient contours before and after optimization at non-design state with α1=0°

圖15 非設(shè)計(jì)狀態(tài)α2=2°優(yōu)化前后空間流場(chǎng)對(duì)比Fig. 15 Comparison of flow field of waverider before and after optimization at non-design state with α2=2°

圖16 非設(shè)計(jì)狀態(tài)α2=2°優(yōu)化前后壓力系數(shù)云圖對(duì)比Fig. 16 Comparison of pressure coefficient contours before and after optimization at non-design state with α2=2°

圖17 優(yōu)化前后各截面外形和表面壓力系數(shù)分布對(duì)比（α1=0°）Fig. 17 Comparison of sectional shapes and surface pressure coefficient distributions before and after optimization （α1=0°）

圖18 優(yōu)化前后各截面外形和表面壓力系數(shù)分布對(duì)比（α2=2°）Fig. 18 Comparison of sectional shapes and surface pressure coefficient distributions before and after optimization （α2=2°）

在乘波體的中后部，優(yōu)化后外形在中部前緣明顯上凸，激波后的氣流流經(jīng)背風(fēng)面又產(chǎn)生一道膨脹波，且相對(duì)于初始外形，膨脹波前移，波后氣流加速，因此中后部背風(fēng)面的壓力更小。相應(yīng)的，優(yōu)化后外形上表面（靠近前緣處）也略微上移，激波強(qiáng)度相對(duì)更強(qiáng)，波后速度更小，升力更大，這也是升力增量的主要來(lái)源。此外，優(yōu)化后外形下表面（靠近對(duì)稱面處）相對(duì)上移，氣流流經(jīng)物面的物面角減小，激波相對(duì)下移，且激波強(qiáng)度減弱，導(dǎo)致波后的壓力略有減小，這使得波阻有所下降。

綜上所述，在非設(shè)計(jì)狀態(tài)下優(yōu)化后外形背風(fēng)面壓力下降，迎風(fēng)面壓力增加，同時(shí)優(yōu)化后的前緣外形能更好地將高壓氣流限制在迎風(fēng)面。因此，優(yōu)化后外形非設(shè)計(jì)狀態(tài)下的升阻比提升明顯。

圖19 和圖20 給出了初始外形和優(yōu)化后外形在設(shè)計(jì)狀態(tài)6°迎角的空間流場(chǎng)及壓力系數(shù)云圖對(duì)比。圖21 給出了乘波體在設(shè)計(jì)迎角狀態(tài)下各站位處的外形變化和表面壓力系數(shù)分布對(duì)比。與單點(diǎn)優(yōu)化不同的是，優(yōu)化后外形在頭部下表面（靠近對(duì)稱面處）先相對(duì)上移，而后又下凸，產(chǎn)生二次激波，在波后出現(xiàn)局部高壓區(qū)，頭部處壓力相對(duì)單點(diǎn)優(yōu)化外形更大；另外，通過(guò)空間流場(chǎng)對(duì)比和表面壓力系數(shù)分布對(duì)比也可以看出，多點(diǎn)優(yōu)化后的外形能夠更好地阻止前緣迎風(fēng)面的高壓氣流泄漏到背風(fēng)面。在乘波體的中后部，與單點(diǎn)優(yōu)化不同的是，優(yōu)化后外形在中后部前緣更加上翹，激波強(qiáng)度相對(duì)更強(qiáng)，波后速度更小，因此前緣處的背風(fēng)面壓力更小，迎風(fēng)面壓力更大；但相應(yīng)的波阻也會(huì)增加，因此阻力也相對(duì)更大。

圖19 設(shè)計(jì)狀態(tài)α3=6°優(yōu)化前后空間流場(chǎng)對(duì)比Fig. 19 Comparison of flow field of waverider before and after optimization at design state with α3=6°

圖20 設(shè)計(jì)狀態(tài)α3=6°優(yōu)化前后壓力系數(shù)云圖對(duì)比Fig. 20 Comparison of pressure coefficient contours before and after optimization at design state with α3=6°

圖21 優(yōu)化前后各截面外形和表面壓力系數(shù)分布對(duì)比（α3=6°）Fig. 21 Comparison of sectional shapes and surface pressure coefficient distributions before and after optimization （α3=6°）

綜上所述，相對(duì)于單點(diǎn)優(yōu)化外形，多點(diǎn)優(yōu)化后的外形在設(shè)計(jì)狀態(tài)下頭部和前緣處的壓力分布明顯改善，同時(shí)優(yōu)化后的前緣外形能更好地將高壓氣流限制在迎風(fēng)面。但是因?yàn)椴ㄗ璧脑黾?，阻力相?duì)更大。因此，多點(diǎn)優(yōu)化時(shí)飛行器在設(shè)計(jì)狀態(tài)下的升阻比提升幅度較小。

5 結(jié) 論

1）本文通過(guò)結(jié)合RANS 求解方法、結(jié)構(gòu)網(wǎng)格變形方法、FFD 外形參數(shù)化方法、離散伴隨方法與序列二次規(guī)劃算法，實(shí)現(xiàn)了基于離散伴隨的高超聲速飛行器氣動(dòng)優(yōu)化設(shè)計(jì)方法。

2）定平面密切錐乘波體的梯度計(jì)算結(jié)果表明，本文所采用的離散伴隨方程法具有較高的梯度求解精度。針對(duì)定平面形狀的密切錐乘波體的氣動(dòng)優(yōu)化設(shè)計(jì)表明，本文所采用的基于離散伴隨的優(yōu)化設(shè)計(jì)方法在針對(duì)大規(guī)模設(shè)計(jì)變量的高超聲速飛行器開(kāi)展氣動(dòng)優(yōu)化時(shí)具有較高的優(yōu)化效率，具體為：在400 萬(wàn)多塊結(jié)構(gòu)網(wǎng)格、600 個(gè)設(shè)計(jì)變量以及303 個(gè)設(shè)計(jì)約束條件下，該方法僅分別花費(fèi)2 240CPU 小時(shí)和3 360CPU 小時(shí)即完成單一設(shè)計(jì)狀態(tài)和多工況的整機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)。

3）定平面密切錐乘波體在單設(shè)計(jì)點(diǎn)優(yōu)化后，阻力下降明顯，升阻比提升了近5%；在多點(diǎn)優(yōu)化設(shè)計(jì)后，可以保證各個(gè)飛行狀態(tài)下的氣動(dòng)性能均有所提升，具體為：多點(diǎn)優(yōu)化得到的構(gòu)型在保證設(shè)計(jì)點(diǎn)狀態(tài)升阻特性不下降的同時(shí)，將非設(shè)計(jì)點(diǎn)的升阻比提升10%以上。

4）本文所實(shí)現(xiàn)的氣動(dòng)優(yōu)化設(shè)計(jì)方法在一定程度上改善了定平面密切錐乘波體的理論局限性，并在工程化外形整機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)中展現(xiàn)出良好的優(yōu)化性能。