連續(xù)伴隨方法在二維高超聲速進(jìn)氣道優(yōu)化中的應(yīng)用

高 昌, 張小慶, 賀元元, 吳穎川, 樂(lè)嘉陵

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 高超聲速?zèng)_壓發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000)

0 引 言

以超燃沖壓發(fā)動(dòng)機(jī)為動(dòng)力的吸氣式高超聲速飛行器是未來(lái)實(shí)現(xiàn)空天飛行和全球快速到達(dá)的重要途徑。進(jìn)氣道作為壓縮部件，其工作效率是影響沖壓發(fā)動(dòng)機(jī)工作狀態(tài)和性能的重要因素。研究表明，對(duì)于Ma=5～7的碳?xì)淙剂习l(fā)動(dòng)機(jī)，進(jìn)氣道壓縮效率每提高1%，發(fā)動(dòng)機(jī)比沖可提高3%～5%[1]。高超聲速進(jìn)氣道是機(jī)體和發(fā)動(dòng)機(jī)強(qiáng)耦合部件之一，已經(jīng)試飛成功的美國(guó)X-43A和X-51A飛行器分別采用了二維曲面壓縮和內(nèi)轉(zhuǎn)進(jìn)氣道。這類進(jìn)氣道是基于無(wú)黏特征線反設(shè)計(jì)得到，在真實(shí)飛行條件下，受黏性和激波邊界層干擾的影響，進(jìn)氣道的性能惡化。對(duì)于二維簡(jiǎn)單構(gòu)型進(jìn)氣道可以采用黏性修正方法提高進(jìn)氣道性能，但對(duì)于三維復(fù)雜構(gòu)型進(jìn)氣道，黏性修正方法很難得到理想的結(jié)果。在結(jié)構(gòu)設(shè)計(jì)上進(jìn)氣道構(gòu)型還應(yīng)滿足吸氣式高超聲速飛行器一體化設(shè)計(jì)的幾何約束[2]。從提高進(jìn)氣道性能和一體化設(shè)計(jì)的角度考慮，進(jìn)氣道優(yōu)化設(shè)計(jì)是吸氣式高超聲速飛行器技術(shù)研究的重要內(nèi)容。

以遺傳算法為代表的隨機(jī)類算法是目前高超聲速飛行器優(yōu)化設(shè)計(jì)所采用的主要方法。遺傳算法魯棒性強(qiáng)，可以得到目標(biāo)函數(shù)的全局最優(yōu)解，但遺傳算法在優(yōu)化循環(huán)中搜索面廣，需要評(píng)價(jià)每個(gè)個(gè)體的優(yōu)劣，計(jì)算開(kāi)銷較大，對(duì)于多設(shè)計(jì)變量的優(yōu)化問(wèn)題存在困難。為了提高優(yōu)化效率，在工程應(yīng)用中通常采用簡(jiǎn)化分析模型或通過(guò)少量的樣本點(diǎn)構(gòu)建代理模型，但針對(duì)多設(shè)計(jì)變量問(wèn)題構(gòu)建精確的代理模型仍需要較大的計(jì)算開(kāi)銷。

另一類優(yōu)化方法是基于目標(biāo)函數(shù)梯度的局部尋優(yōu)算法，如擬牛頓方法和逐步二次規(guī)劃方法，這一類算法是對(duì)目標(biāo)函數(shù)采用局部二次多項(xiàng)式近似，根據(jù)目標(biāo)函數(shù)的梯度信息確定尋優(yōu)方向，逐步向目標(biāo)函數(shù)的局部最優(yōu)逼近。梯度信息的獲取可以通過(guò)有限差分法，有限差分法需要分別計(jì)算目標(biāo)函數(shù)對(duì)每個(gè)設(shè)計(jì)變量的梯度，計(jì)算開(kāi)銷隨設(shè)計(jì)變量的數(shù)目增長(zhǎng)。另一種方法是伴隨方法，伴隨方法通過(guò)構(gòu)建流動(dòng)控制方程的伴隨方程，通過(guò)一次伴隨方程的求解可以獲得目標(biāo)函數(shù)對(duì)所有設(shè)計(jì)變量的梯度值，計(jì)算開(kāi)銷與設(shè)計(jì)變量數(shù)目無(wú)關(guān)，伴隨方程的求解計(jì)算量與流動(dòng)控制方程大致相當(dāng)。

伴隨方程最早由Pironneau在控制理論中引入[3]，Jameson從應(yīng)用的角度出發(fā)，推導(dǎo)了勢(shì)流方程和Euler方程的伴隨方程，首次給出了基于伴隨方法的二維翼型反設(shè)計(jì)算例[4]，隨后Jameson[5-7]、Anderson[8]和Giles[9-10]等對(duì)連續(xù)伴隨方法和離散伴隨方法的推廣以及在非結(jié)構(gòu)網(wǎng)格上的應(yīng)用都做出了大量的工作，目前伴隨方法在CFD誤差評(píng)估[11]、不確定度分析[12]和自適應(yīng)網(wǎng)格[13]等問(wèn)題中都有重要應(yīng)用。

在氣動(dòng)優(yōu)化研究中，伴隨方法主要應(yīng)用于翼型、葉輪機(jī)械的優(yōu)化[14]。目前，國(guó)外研究人員已經(jīng)開(kāi)展了伴隨方法在高超聲速優(yōu)化問(wèn)題中的應(yīng)用探索，Copeland給出了考慮化學(xué)非平衡NS方程的連續(xù)伴隨方程的推導(dǎo)和求解，并針對(duì)再入航天器熱流優(yōu)化問(wèn)題對(duì)比了伴隨方法和有限差分法計(jì)算得到的梯度值[15]；Kline采用連續(xù)伴隨方法針對(duì)單楔壓縮的二維和準(zhǔn)三維高超聲速進(jìn)氣道優(yōu)化問(wèn)題進(jìn)行了研究[16-17]。

伴隨方法在高超聲速進(jìn)氣道優(yōu)化中的應(yīng)用處于初步嘗試和驗(yàn)證階段，為了進(jìn)一步探索伴隨方法的應(yīng)用潛力，本文選取二維曲面壓縮無(wú)黏進(jìn)氣道構(gòu)型，采用FFD方法對(duì)進(jìn)氣道外壓縮面實(shí)現(xiàn)參數(shù)化控制，在考慮黏性的情況下以進(jìn)氣道出口流量為目標(biāo)函數(shù)求解連續(xù)伴隨方程，采用三套網(wǎng)格驗(yàn)證網(wǎng)格密度對(duì)外壓縮面壁面靈敏度的影響，通過(guò)與有限差分法對(duì)比驗(yàn)證連續(xù)伴隨方法計(jì)算結(jié)果的可靠性。選取外壓縮面為優(yōu)化對(duì)象，采用BFGS算法優(yōu)化進(jìn)氣道出口流量，驗(yàn)證伴隨方法在高超聲速進(jìn)氣道優(yōu)化中的適用性。

1 二維曲面壓縮進(jìn)氣道

二維曲面壓縮進(jìn)氣道設(shè)計(jì)采用特征線理論[18]，根據(jù)求解問(wèn)題類型的不同，進(jìn)氣道分為4個(gè)區(qū)進(jìn)行反設(shè)計(jì)：前緣激波依賴區(qū)(A區(qū))、主壓縮區(qū)(B區(qū))、末端激波依賴區(qū)(C區(qū))和穩(wěn)定區(qū)(D區(qū))(圖1)。A區(qū)的求解依賴于來(lái)流條件、初始?jí)嚎s角和壁面壓升規(guī)律，B區(qū)的求解依賴于A區(qū)出口參數(shù)分布和壁面壓升規(guī)律，C區(qū)的求解依賴于B區(qū)出口參數(shù)分布并在肩點(diǎn)c處實(shí)現(xiàn)反射激波消波，D區(qū)對(duì)氣流方向進(jìn)一步梳理。設(shè)計(jì)時(shí)，給定的設(shè)計(jì)馬赫數(shù)為6.3，設(shè)計(jì)高度為27 km，初始?jí)嚎s角為6°，壁面壓升規(guī)律選擇雙曲正切函數(shù)。設(shè)計(jì)得到的進(jìn)氣道總長(zhǎng)為L(zhǎng)=849.8 mm，捕獲高度為H=139.1 mm。圖2給出了特征線方法和Euler方程計(jì)算得到的進(jìn)氣道流場(chǎng)馬赫數(shù)分布，Euler方程計(jì)算得到的流場(chǎng)結(jié)構(gòu)與特征線方法符合較好，在Euler方程計(jì)算結(jié)果中，內(nèi)壓縮段肩點(diǎn)c處出現(xiàn)了微弱的反射波。

圖1 進(jìn)氣道設(shè)計(jì)方法Fig.1 Inlet design method

圖2 特征線法和Euler方程計(jì)算結(jié)果Fig.2 Flowfield of characteristics method and Euler equation

2 伴隨優(yōu)化方法

2.1 伴隨方程

連續(xù)伴隨方法是基于流動(dòng)控制方程推導(dǎo)而來(lái)，以進(jìn)氣道出口流量為目標(biāo)函數(shù)，其可表示為進(jìn)氣道出口的積分:

(1)

定常流動(dòng)控制方程的差分形式表示為:

(2)

其中R為殘差，U為流場(chǎng)守恒變量，A為Jacobi矩陣。優(yōu)化問(wèn)題可以等價(jià)為以流動(dòng)控制方程為約束，求目標(biāo)函數(shù)的極值。通過(guò)Lagrange乘子將流動(dòng)控制方程引入目標(biāo)函數(shù):

(3)

目標(biāo)函數(shù)關(guān)于壁面的靈敏度?T/?S表示了壁面法向變形對(duì)目標(biāo)函數(shù)的影響量，可以用目標(biāo)函數(shù)的變分δT/δS進(jìn)行近似:

(4)

運(yùn)用變分原理將流動(dòng)控制方程代入公式(4)可得:

(5)

為了計(jì)算δT/δS，可以通過(guò)構(gòu)建伴隨方程及相應(yīng)的邊界條件來(lái)消除公式(5)中包含δU的項(xiàng):

(6)

上式即為流動(dòng)控制方程(2)的伴隨方程，ψ稱為伴隨變量。通過(guò)構(gòu)建伴隨方程和相應(yīng)的邊界條件，目標(biāo)函數(shù)的靈敏度就可以表示為:

(7)

在應(yīng)用伴隨方法時(shí)，首先對(duì)流動(dòng)控制方程進(jìn)行求解得到的流場(chǎng)變量，代入伴隨方程，迭代求解得到伴隨變量，將伴隨變量代入公式(7)即得到進(jìn)氣道流量關(guān)于壁面的靈敏度信息。

2.2 曲面參數(shù)化和網(wǎng)格變形

在氣動(dòng)外形優(yōu)化問(wèn)題中，需要通過(guò)正交、完備的基函數(shù)對(duì)優(yōu)化型面的設(shè)計(jì)空間進(jìn)行參數(shù)化表示，可以減少優(yōu)化過(guò)程中的設(shè)計(jì)變量，同時(shí)在優(yōu)化循環(huán)中保證型面的光滑和連續(xù)。采用FFD方法對(duì)進(jìn)氣道外壓縮面進(jìn)行參數(shù)化控制，F(xiàn)FD方法通過(guò)構(gòu)建包裹變形曲面的控制體，利用基函數(shù)建立變形曲面與控制點(diǎn)的映射關(guān)系:

(8)

圖3 FFD控制體(黑色為變形前，紅色為變形后)Fig.3 FFD control box(black：original one, red：deformed one)

網(wǎng)格變形采用彈性比擬法，將計(jì)算網(wǎng)格等效為彈性體，以外壓縮面的位移為邊界條件來(lái)求解網(wǎng)格點(diǎn)的變形量，通過(guò)求解線彈性Poisson方程得到各個(gè)網(wǎng)格點(diǎn)的位移量:

(9)

其中Tr(ε)為ε的跡，λ=νE(1-2ν)/(1+ν)，μ=2E/(1+ν)，E為彈性模型，ν為泊松比。在求解網(wǎng)格變形過(guò)程中，邊界層和局部加密區(qū)域網(wǎng)格由于網(wǎng)格縱橫比大、體積小，這些網(wǎng)格容易過(guò)渡扭曲產(chǎn)生負(fù)體積，在求解過(guò)程中采用隨網(wǎng)格體積成反比變化的彈性模量，增強(qiáng)這些網(wǎng)格的剛度，保持原有的網(wǎng)格質(zhì)量。

3 外壓縮面靈敏度分析

進(jìn)氣道的初始構(gòu)型為采用特征線法設(shè)計(jì)的無(wú)黏構(gòu)型，在黏性情況下由于受邊界層擠壓的影響，唇口處會(huì)產(chǎn)生溢流，外壓縮面形狀影響到前體流場(chǎng)波系結(jié)構(gòu)進(jìn)而影響進(jìn)氣道的捕獲流量。計(jì)算構(gòu)型前緣和唇口進(jìn)行倒圓(倒圓直徑1 mm)，采用三套網(wǎng)格(如圖4所示)對(duì)目標(biāo)函數(shù)進(jìn)氣道流量關(guān)于外壓縮面靈敏度的影響因素進(jìn)行分析。第一套網(wǎng)格(Grid 1)壁面網(wǎng)格第一層高度為0.1 mm，網(wǎng)格總數(shù)為3.9萬(wàn)；第二套網(wǎng)格(Grid 2)壁面第一層網(wǎng)格高度為0.01 mm，流場(chǎng)中間網(wǎng)格分布與第一套網(wǎng)格基本保持一致，網(wǎng)格總數(shù)為4.6萬(wàn)；第三套網(wǎng)格(Grid 3)在第二套網(wǎng)格的基礎(chǔ)上在對(duì)流場(chǎng)中間網(wǎng)格進(jìn)行了加密，壁面網(wǎng)格基本保持不變，網(wǎng)格總數(shù)為6.5萬(wàn)。

(a) Grid 1

(b) Grid 2

(c) Grid 3

在計(jì)算外壓縮面靈敏度時(shí)，伴隨方程考慮了湍流黏性項(xiàng)的影響。流動(dòng)控制方程和伴隨方程求解過(guò)程中，無(wú)黏通量計(jì)算采用Roe格式，湍流模型采用SA模型。

圖5給出了在三套網(wǎng)格上流動(dòng)控制方程(Flow)求解過(guò)程中進(jìn)氣道流量的收斂歷程和伴隨方程(Adjoint)求解過(guò)程中外壓縮面總靈敏度的收斂歷程。對(duì)于同一套網(wǎng)格，流動(dòng)控制方程和伴隨方程的收斂歷程保持一致；三套網(wǎng)格得到的進(jìn)氣道出口流量相差很小，Grid 1計(jì)算得到的外壓縮面總靈敏度相對(duì)偏?。蝗拙W(wǎng)格迭代60 000步后流動(dòng)控制方程和伴隨方程都基本收斂。

圖5 流動(dòng)控制方程和伴隨方程收斂歷程Fig.5 Convergence history of governing equations and adjoint equation

圖6給出了三套網(wǎng)格計(jì)算得到的外壓縮面靈敏度的分布，在進(jìn)氣道前緣處受前緣曲率的影響，靈敏度出現(xiàn)異常峰值。為了避免對(duì)優(yōu)化過(guò)程產(chǎn)生誤導(dǎo)，在伴隨方程求解過(guò)程中設(shè)置光順因子(Smooth Ratio，SR)，將前緣附近的靈敏度進(jìn)行光順處理。對(duì)比Grid 2得到的兩條靈敏度曲線(Grid 2SR=0.01和Grid 2SR=0.03)，當(dāng)光順因子較小時(shí)(SR=0.01)，前緣附近出現(xiàn)了明顯的異常峰值，當(dāng)設(shè)置光順因子SR=0.03，前緣峰值消失，同時(shí)正常峰值(x=0.47 m處)也有所降低。

圖6 外壓縮面靈敏度分布Fig.6 External compression surface sensitivity

各條靈敏度曲線在x=0.47m處均出現(xiàn)峰值，分析進(jìn)氣道流場(chǎng)，該處產(chǎn)生的壓縮波與入射激波在唇口處交匯，直接影響了進(jìn)氣道的捕獲流量。在峰值之后產(chǎn)生的壓縮波入射進(jìn)入進(jìn)氣道內(nèi)部，對(duì)捕獲流量不再產(chǎn)生影響，所以當(dāng)x>0.47 m后靈敏度迅速下降到0。Grid 2和Grid 3計(jì)算得到的靈敏度分布基本一致，Grid 1計(jì)算得到峰值相對(duì)較小，說(shuō)明壁面靈敏度結(jié)果對(duì)流場(chǎng)中間網(wǎng)格分布不敏感，對(duì)邊界層網(wǎng)格分布依賴性較強(qiáng)。

為了對(duì)靈敏度的合理性進(jìn)行驗(yàn)證，根據(jù)公式(11)將靈敏度分布映射到9個(gè)設(shè)計(jì)變量，得到目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的梯度。同時(shí)根據(jù)有限差分法，分別在每個(gè)設(shè)計(jì)變量上施加差分量，計(jì)算目標(biāo)函數(shù)對(duì)該設(shè)計(jì)變量的梯度，差分步長(zhǎng)分別取1×10-3m，5×10-4m和2.5×10-4m，計(jì)算得到的結(jié)果如圖7所示。在不同步長(zhǎng)下差分方法得到的結(jié)果相互吻合，說(shuō)明有限差分結(jié)果具備收斂性。對(duì)于伴隨方法，Grid 2和Grid 3得到的結(jié)果略高于差分方法，變化規(guī)律一致，第一套網(wǎng)格得到的結(jié)果明顯偏小，說(shuō)明伴隨方法可以得到有效的靈敏度信息。

圖7 目標(biāo)函數(shù)梯度Fig.7 Objective function gradient

4 優(yōu)化結(jié)果

根據(jù)以上建立的參數(shù)化幾何模型和伴隨方程求解方法，以進(jìn)氣道流量為目標(biāo)，對(duì)進(jìn)氣道外壓縮面進(jìn)行優(yōu)化，計(jì)算網(wǎng)格選擇第二套網(wǎng)格，優(yōu)化算法采用BFGS方法，BFGS方法是一種改進(jìn)的擬牛頓法，是利用目標(biāo)函數(shù)梯度求解無(wú)約束優(yōu)化問(wèn)題的最有效的一類方法。圖8給出了進(jìn)氣道流量的優(yōu)化歷程,在第3次優(yōu)化循環(huán)后流量達(dá)到穩(wěn)定，流量由12.37 kg/s增加到13.15 kg/s，提高了6.3%。

圖9給出了優(yōu)化前后的流場(chǎng)馬赫數(shù)分布。在優(yōu)化前，對(duì)于無(wú)黏進(jìn)氣道構(gòu)型，在黏性情況下氣流受到邊界層排擠，唇口處產(chǎn)生了溢流。優(yōu)化后，入射激波位置上移實(shí)現(xiàn)封口，消除了溢流，基本實(shí)現(xiàn)了全流量捕獲。圖10給出了優(yōu)化前后外壓縮面型線和壓力系數(shù)的分布變化，優(yōu)化后外壓縮面前部出現(xiàn)了略微的下凸，在中后部出現(xiàn)明顯的內(nèi)凹，與靈敏度的分布規(guī)律一致，最大法向變形為5.6 mm。優(yōu)化后外壓縮面壓力梯度整體增大，外壓縮面壓縮能力提升。表1給出了優(yōu)化前后進(jìn)氣道出口流量系數(shù)φ，質(zhì)量平均馬赫數(shù)Maexit，壓升比pexit/p∞和總壓恢復(fù)系數(shù)σ的變化，優(yōu)化后進(jìn)氣道出口馬赫數(shù)提高約2.4%，由于外壓縮面壓縮能力的提升，進(jìn)氣道出口壓升比提高約5.3%，外壓縮面以等熵壓縮為主，整體上減少了進(jìn)氣道的總壓損失，總壓恢復(fù)系數(shù)提高了約10.6%，進(jìn)氣道整體性能得到提高。

圖8 進(jìn)氣道出口流量變化Fig.8 Mass flow rate in optimization evaluations

(a) 優(yōu)化前

(b) 優(yōu)化后

圖10 優(yōu)化前后外壓縮型面和壓力系數(shù)分布Fig.10 Original and optimized external compression surfaces and Cp distribution

表1 進(jìn)氣道出口參數(shù)Table 1 Parameters of the inlet exit

5 結(jié) 論

本文針對(duì)連續(xù)伴隨方法在二維高超聲速進(jìn)氣道優(yōu)化設(shè)計(jì)中的應(yīng)用開(kāi)展了研究。進(jìn)氣道構(gòu)型采用基于特征線法設(shè)計(jì)的曲面壓縮進(jìn)氣道。為了分析壁面靈敏度的影響因素，采用了三套疏密不同網(wǎng)格進(jìn)行計(jì)算對(duì)比分析。對(duì)于目標(biāo)函數(shù)靈敏度，采用有限差分結(jié)果進(jìn)行對(duì)比驗(yàn)證?；谒捎玫倪B續(xù)伴隨方法和FFD參數(shù)化方法，以進(jìn)氣道外壓縮面為優(yōu)化對(duì)象，采用BFGS算法對(duì)進(jìn)氣道流量進(jìn)行優(yōu)化，得到了以下結(jié)論：

1) 以進(jìn)氣道流量為目標(biāo)函數(shù)，伴隨方法可以有效得到外壓縮面靈敏度信息，靈敏度對(duì)邊界層網(wǎng)格分布有較強(qiáng)的依賴性。

2) 采用BFGS算法對(duì)進(jìn)氣道流量?jī)?yōu)化，優(yōu)化后唇口處實(shí)現(xiàn)激波封口，流量提高6.3%，外壓縮面最大法向位移為5.6mm，外壓縮面壓縮能力增強(qiáng)，進(jìn)氣道整體性能得到提高。

3) 從優(yōu)化結(jié)果來(lái)看，伴隨方法可有效應(yīng)用于高超聲速進(jìn)氣道優(yōu)化。