王 澤，王梓伊，王 旭，宋述芳，張偉偉，*

(1. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000)

0 引 言

高超聲速飛行器是現(xiàn)今各國(guó)研究的熱點(diǎn)，無(wú)論是飛行器的熱防護(hù)設(shè)計(jì)還是熱氣動(dòng)彈性分析，氣動(dòng)熱的高效精準(zhǔn)預(yù)示都是極為關(guān)鍵的問(wèn)題[1-2]。傳統(tǒng)的氣動(dòng)熱預(yù)示手段主要包括飛行試驗(yàn)、地面風(fēng)洞試驗(yàn)、工程算法和數(shù)值模擬等方法。飛行試驗(yàn)和地面風(fēng)洞試驗(yàn)雖然可以準(zhǔn)確預(yù)示飛行器的氣動(dòng)熱，但是受成本限制，難以作為常規(guī)預(yù)示手段[3-4]。具有試驗(yàn)和理論基礎(chǔ)的工程方法可以高效地進(jìn)行氣動(dòng)熱預(yù)示，但是因?yàn)槠浒T多假設(shè)，對(duì)復(fù)雜外形、流動(dòng)的氣動(dòng)熱預(yù)示精度低。數(shù)值模擬方法雖然可以實(shí)現(xiàn)對(duì)復(fù)雜外形飛行器較高精度的氣動(dòng)熱預(yù)示，但是對(duì)網(wǎng)格質(zhì)量要求高、收斂性差、計(jì)算量大，依舊無(wú)法滿足工程設(shè)計(jì)中氣動(dòng)熱的快速精準(zhǔn)預(yù)示需求[5-6]。

為解決高超聲速飛行器氣動(dòng)熱高效精準(zhǔn)預(yù)示問(wèn)題，研究者們進(jìn)行了大量探索。呂麗麗等[7]通過(guò)在邊界層外求解Euler 方程獲得無(wú)黏參數(shù)，邊界層內(nèi)使用工程算法預(yù)示熱流，提高了工程算法對(duì)復(fù)雜外形的適用性。Mcnamara 等[8]通過(guò)兩次數(shù)值求解獲得物面絕熱溫度分布Taw和對(duì)流換熱系數(shù)h，再利用Newton 冷卻公式可以快速獲得任意壁面溫度Tw下的熱流分布。Falkiewicz[9]通過(guò)本征正交分解（proper orthogonal decomposition, POD）對(duì) F104 機(jī)翼表面瞬態(tài)溫度場(chǎng)和熱流場(chǎng)進(jìn)行降階，構(gòu)建輸入為來(lái)流參數(shù)和溫度場(chǎng)模態(tài)系數(shù)、輸出為氣動(dòng)熱降階分布的Kriging 代理模型，實(shí)現(xiàn)氣動(dòng)熱的快速預(yù)示。陳鑫等[10-12]構(gòu)建了輸入為來(lái)流參數(shù)、輸出為F104 機(jī)翼物面溫度場(chǎng)POD 基系數(shù)的代理模型，并對(duì)該模型的構(gòu)建方法和采樣方法進(jìn)行了深入研究。聶春生等[13]結(jié)合POD 方法和RBF 插值，建立了一種適宜復(fù)雜外形表面熱流預(yù)示的代理模型，實(shí)現(xiàn)了沿給定彈道三維熱環(huán)境快速預(yù)示。張智超等[14]提出了一種基于徑向基神經(jīng)網(wǎng)絡(luò)的氣動(dòng)熱快速預(yù)示代理模型方法，輸入為來(lái)流參數(shù)，輸出為每個(gè)網(wǎng)格節(jié)點(diǎn)熱流預(yù)示結(jié)果。以上幾種思路雖然都實(shí)現(xiàn)了氣動(dòng)熱的快速精準(zhǔn)預(yù)示，但是仍存在一些弊端：CFD 結(jié)合工程算法，雖然提升了對(duì)復(fù)雜外形的勝任性，但因工程模型的固有假設(shè)，精度仍有限；POD 降階方法需要大量的訓(xùn)練樣本，而且很難實(shí)現(xiàn)對(duì)幾何外形的泛化；使用來(lái)流參數(shù)映射壁面熱流，無(wú)法體現(xiàn)氣動(dòng)熱的局部性特點(diǎn)。

近些年，機(jī)器學(xué)習(xí)技術(shù)得到了飛速發(fā)展，特別是可以依據(jù)有限樣本信息實(shí)現(xiàn)最大泛化性能的支持向量機(jī)（support vector machine, SVM）已經(jīng)被用于如偏微分方程求解[15]、圖像識(shí)別[16]和氣動(dòng)力建模[17-19]等諸多領(lǐng)域。

本文提出了一種基于支持向量機(jī)和邊界層理論的當(dāng)?shù)鼗瘹鈩?dòng)熱預(yù)示建模方法（boundary layer theory and support vector machine localized aeroheating modeling method, BSLM）。由邊界層理論得知，壁面熱流具有很強(qiáng)的當(dāng)?shù)匦?，是一個(gè)局部量，而且其分布很大程度上取決于邊界層外緣信息[20]?；谶吔鐚油饩壭畔⒑捅诿鏌崃鏖g的密切關(guān)聯(lián)，本文采用了數(shù)據(jù)驅(qū)動(dòng)的當(dāng)?shù)鼗Ｋ悸?，利用SVM 構(gòu)建了從當(dāng)?shù)剡吔鐚油饩壧卣鞯疆?dāng)?shù)乇诿鏌崃鞯臍鈩?dòng)熱預(yù)示模型。論文首先探究了邊界層外緣特征對(duì)模型性能的影響，驗(yàn)證了特征選擇方法的有效性；然后針對(duì)二級(jí)壓縮面及雙橢球開展氣動(dòng)熱預(yù)示，分析了模型預(yù)示精度、泛化能力以及對(duì)非均勻分布壁面溫度等當(dāng)?shù)鼗吔鐥l件的適用性，并與傳統(tǒng)POD 降階模型方法進(jìn)行了簡(jiǎn)要對(duì)比；最后，總結(jié)了研究結(jié)果，給出了相關(guān)結(jié)論。

1 氣動(dòng)熱預(yù)示模型

1.1 支持向量機(jī)

SVM 是一種建立在VC 維（Vapnik-Chervonenkis dimension）理論和結(jié)構(gòu)風(fēng)險(xiǎn)最小原理基礎(chǔ)上的機(jī)器學(xué)習(xí)算法[21]。其核心思想是將回歸問(wèn)題經(jīng)過(guò)嚴(yán)格數(shù)學(xué)推導(dǎo)轉(zhuǎn)化為二次規(guī)劃問(wèn)題，并通過(guò)非線性映射將樣本從原始空間映射到高維特征空間實(shí)現(xiàn)樣本的線性可分。設(shè)樣本訓(xùn)練集D={(x1,y1), (x2,y2), ···, (xm,ym)},xi∈Rn,yi∈R,i=1, 2, ···,m，xi為真實(shí)輸入，yi為真實(shí)輸出?？紤]近似回歸模型f(x)=ωT?(x)+b， ω和b是模型參數(shù)，?(x)是樣本x的非線性映射。給定真實(shí)輸出y和模型輸出f(x)的容忍偏差?，支持向量回歸問(wèn)題可寫為：

其中，C為正則化常數(shù)，??為損失函數(shù)：

式（1）引入拉格朗日乘子αi≥0 、≥0，由拉格朗日乘子法得到式（1）的對(duì)偶問(wèn)題：

對(duì)應(yīng)KKT（Karush-Kuhn-Tucker）條件為：

其中ξi代表松弛變量，是支持向量機(jī)的優(yōu)化參數(shù)。求解該對(duì)偶問(wèn)題，得到拉格朗日乘子αi、及參數(shù)b，則近似回歸模型f(x)可表示為：

1.2 特征選擇

輸入特征對(duì)機(jī)器學(xué)習(xí)模型性能的影響十分關(guān)鍵。特征不足難以全面展現(xiàn)數(shù)據(jù)的特性，不利于模型精度的提高；相反，特征冗余可能使模型過(guò)擬合，降低模型泛化性能[22]。為了挑選出有效特征，需要進(jìn)行特征選擇。常用的特征選擇方法主要有經(jīng)驗(yàn)選擇法、逐一遍歷法、過(guò)濾法、打包法和嵌入法等。本文使用過(guò)濾法和嵌入法結(jié)合的特征選擇方法：首先對(duì)候選特征進(jìn)行重要性排序，為避免某些包含關(guān)鍵信息的特征被丟棄，需將重要性低但包含關(guān)鍵信息的特征排序人為提前；然后從最重要的特征開始按順序逐個(gè)地增加特征直到模型精度不再提高或遍歷所有特征；最后選擇模型預(yù)示性能最好的特征組合。其中，特征重要性排序使用隨機(jī)森林（random forest, RF）。隨機(jī)森林是一種并行式集成機(jī)器學(xué)習(xí)方法，內(nèi)嵌平均不純度減少算法用于計(jì)算特征的重要性。本文特征選擇方法具體流程如圖1 所示。

圖1 特征選擇流程圖Fig. 1 Flow diagram of the feature selection

1.3 模型構(gòu)建

氣動(dòng)熱預(yù)示模型的構(gòu)建過(guò)程如圖2 所示。首先，求解Euler 方程獲得邊界層外緣信息。之后，根據(jù)邊界層外緣信息構(gòu)建出候選特征并進(jìn)行特征選擇，得到輸入特征組合。最后，使用訓(xùn)練數(shù)據(jù)對(duì)支持向量機(jī)進(jìn)行訓(xùn)練，并測(cè)試模型性能。至此，氣動(dòng)熱預(yù)示模型構(gòu)建完畢。

圖2 氣動(dòng)熱預(yù)示模型構(gòu)建流程圖Fig. 2 Flow diagram of the aeroheating prediction model construction

2 算例驗(yàn)證及結(jié)果分析

2.1 氣動(dòng)熱數(shù)值計(jì)算方法驗(yàn)證

本文使用RANS 計(jì)算氣動(dòng)熱樣本數(shù)據(jù)，為保證數(shù)據(jù)精度需對(duì)RANS 計(jì)算方法進(jìn)行驗(yàn)證。二維外形氣動(dòng)熱計(jì)算采用AUSM+空間格式、 LU-SGS 偽時(shí)間推進(jìn)和k?ωSST湍流模型，CFL 數(shù)取0.5，以高超聲速圓柱作為驗(yàn)證算例[23]；三維外形采用Steger 空間格式、Minmod 限制器和k?ωSST湍流模型，CFL 數(shù)取0.5，以雙橢球作為驗(yàn)證算例。

高超聲速圓柱半徑0.038 m，來(lái)流馬赫數(shù)為6.5，雷諾數(shù)1.6×106/m，總溫1 900 K，壁面溫度294 K[24]。雙橢球標(biāo)模算例與文獻(xiàn)[25]一致，來(lái)流馬赫數(shù)為10.02，雷諾數(shù)2.2×106/m，總溫1 457 K，壁面溫度294 K，迎角?5°。二維圓柱和三維雙橢球計(jì)算網(wǎng)格如圖3 所示，第一層網(wǎng)格高度都為1×10?6m。圖4 和圖5 分別給出了圓柱和雙橢球的壁面熱流RANS 計(jì)算結(jié)果，兩者都與試驗(yàn)結(jié)果高度吻合，本文采用的數(shù)值計(jì)算方法可以準(zhǔn)確計(jì)算二維及三維外形的壁面熱流。

圖3 高超聲速圓柱和雙橢球計(jì)算網(wǎng)格Fig. 3 Computational grids of the hypersonic cylinder and the double ellipsoid

圖4 圓柱壁面熱流分布計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比Fig. 4 Heat flux distribution on the cylinder surface compared between the computational and experimental results

圖5 雙橢球表面中心線熱流分布計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比Fig. 5 Heat flux distribution along the centerline of the double ellipsoid surface compared between the computational and experimental results

2.2 雙橢球上表面中心線熱流預(yù)示

本節(jié)利用雙橢球上表面中心線熱流預(yù)示算例對(duì)特征選擇方法進(jìn)行驗(yàn)證，對(duì)本文氣動(dòng)熱建模方法（BSLM）進(jìn)行初步測(cè)試，并與傳統(tǒng)POD 降階模型方法進(jìn)行簡(jiǎn)單對(duì)比。

在每個(gè)迎角的雙橢球上表面中心線抽取65 個(gè)點(diǎn)的邊界層外緣信息和熱流數(shù)據(jù)，取點(diǎn)位置及數(shù)量需保證表面熱流的準(zhǔn)確刻畫。以迎角0°、10°、20°和25°的數(shù)據(jù)作為訓(xùn)練集，迎角15°的數(shù)據(jù)作為驗(yàn)證集，迎角5°和30°的數(shù)據(jù)作為測(cè)試集。根據(jù)邊界層外緣信息和泛化對(duì)象構(gòu)造候選特征，如表1 所示。s為流線長(zhǎng)度，即雙橢球上表面中心線任意一點(diǎn)沿壁面到前緣點(diǎn)的距離；αs=tan?1(ve/ue)為當(dāng)?shù)赜?；下?biāo)e 代表邊界層外緣，下標(biāo) ∞代表來(lái)流，下標(biāo)w 代表壁面。

表1 候選特征Table 1 Feature candidates

使用隨機(jī)森林對(duì)候選特征進(jìn)行重要性排序，結(jié)果如圖6 所示。特征ρe/ρ∞和特征Pe/P∞的重要性大小遠(yuǎn)超其他特征，這兩個(gè)特征作為首要選擇。特征s和αs分別表示流線長(zhǎng)度和當(dāng)?shù)赜?，雖然包含信息簡(jiǎn)單、特征重要性低，但是它們與其他特征組合后能刻畫非常豐富的流場(chǎng)特性，是非常關(guān)鍵的特征。因此，根據(jù)經(jīng)驗(yàn)將特征s和αs的排序提前到第3 位和第4 位，最終特征排序?yàn)閒2＞f3＞f1＞f6＞f5＞f8＞f9＞f4＞f7。以f2 和f3 兩個(gè)特征為初始特征組合，按順序逐個(gè)增加特征，進(jìn)行建模得到訓(xùn)練集和驗(yàn)證集的均方根誤差（RMSE），結(jié)果如圖7 所示。可以看出，隨著特征的逐個(gè)增加，訓(xùn)練集和驗(yàn)證集的RMSE 值呈現(xiàn)先減小后增大的趨勢(shì)，特征數(shù)量為5 的特征組合對(duì)應(yīng)著最小的RMSE 值。因此，本算例選擇特征組合[f2，f3，f1，f6，f5]進(jìn)行熱流預(yù)示。

圖6 特征重要性柱狀圖Fig. 6 Bar chart for the feature importance

圖7 不同特征組合下模型訓(xùn)練和驗(yàn)證精度Fig. 7 Model training and validation accuracy under different feature combinations

POD 方法使用的熱流數(shù)據(jù)集與BSLM 方法相同，代理模型采用Kriging 模型。不同迎角雙橢球上表面中心線BSLM 熱流預(yù)示結(jié)果、POD 熱流預(yù)示結(jié)果和RANS 結(jié)果對(duì)比如圖8 所示，可以看出，BSLM方法預(yù)示的熱流曲線與RANS 計(jì)算的熱流曲線吻合很好，而且在前緣點(diǎn)處BSLM 方法預(yù)示結(jié)果明顯好于POD 方法預(yù)示結(jié)果。如表2 所示，三組迎角下前緣點(diǎn)處BSLM 方法預(yù)示的熱流值與RANS 計(jì)算值的相對(duì)誤差（RE）都小于POD 方法；當(dāng)α=30°時(shí)，BSLM方法的RE 值僅有POD 方法的8.71%。三組迎角下BSLM 方法預(yù)示的熱流值與RANS 計(jì)算值的均方根誤差（RMSE）也均小于POD 方法；當(dāng)α=30°時(shí)，BSLM方法的RMSE 值僅為POD 方法的18.80%。BSLM 方法預(yù)示結(jié)果精度優(yōu)于POD 方法的原因主要是：BSLM 方法是一種當(dāng)?shù)鼗７椒ǎ蠠崃鞯漠?dāng)?shù)匦蕴攸c(diǎn)，而POD 方法構(gòu)建的是全局模型；BSLM 方法不僅使用了熱流信息，還使用了邊界層外緣信息，相比于POD 方法信息更豐富。

表2 BSLM 和POD 方法的RMSE 值和前緣點(diǎn)RE 值對(duì)比Table 2 RE at the leading-edge point and RMSE compared between BSLM and POD

圖8 BSLM、POD 熱流預(yù)示值和RANS 計(jì)算值對(duì)比Fig. 8 Heat flux comparison between the calculation of RANS and the prediction of BSLM and POD

2.3 二級(jí)壓縮面熱流預(yù)示

在實(shí)際氣動(dòng)加熱過(guò)程中，結(jié)構(gòu)表面溫度往往是非均勻分布的，因此，有必要進(jìn)一步考察非均布壁面溫度條件下BSLM 方法的適用性。

采用二級(jí)壓縮面作為測(cè)試算例，該算例改動(dòng)自某熱氣動(dòng)彈性實(shí)驗(yàn)?zāi)Ｐ?，其幾何尺寸和?jì)算網(wǎng)格如圖9 所示，第一層網(wǎng)格高度1×10?3mm。來(lái)流參數(shù)選自30 km 高度的大氣，以馬赫數(shù)Ma和壁面溫度Tw作為變量，數(shù)值計(jì)算格式與二維圓柱相同。訓(xùn)練狀態(tài)和測(cè)試狀態(tài)如表3 所示，標(biāo)注“√”的狀態(tài)為訓(xùn)練狀態(tài)，標(biāo)注“△”的狀態(tài)為測(cè)試狀態(tài)。此外，測(cè)試狀態(tài)還包括Ma= 9、Tw非均勻分布的狀態(tài)，Tw在二級(jí)壓縮面的分布如圖10 所示。

表3 訓(xùn)練和測(cè)試狀態(tài)劃分Table 3 Training and test states setting

圖9 二級(jí)壓縮面幾何尺寸和計(jì)算網(wǎng)格（單位：mm）Fig. 9 Geometry dimension and computational grid of the twostage compression surface (unit：mm)

圖10 二級(jí)壓縮面壁面溫度分布Fig. 10 Temperature distribution on the two-stage compression surface

在每個(gè)狀態(tài)的二級(jí)壓縮面上選擇78 個(gè)點(diǎn)的熱流數(shù)據(jù)和邊界層外緣信息。構(gòu)造候選特征，進(jìn)行特征選擇，得到特征組合完成模型構(gòu)建，進(jìn)行模型測(cè)試。二級(jí)壓縮面熱流預(yù)示結(jié)果和RANS 計(jì)算結(jié)果對(duì)比如圖11 所示，可以看出，熱流預(yù)示值和RANS 計(jì)算值基本吻合。二級(jí)壓縮面上一些監(jiān)測(cè)點(diǎn)的熱流預(yù)示值與RANS 計(jì)算值的對(duì)比結(jié)果如表4 所示，相對(duì)誤差均小于5%。由以上結(jié)果可知，BSLM 方法可以較為精確地進(jìn)行氣動(dòng)熱預(yù)示，適用于非均布壁面溫度邊界條件，且具備對(duì)來(lái)流馬赫數(shù)和壁面溫度的泛化性。BSLM 氣動(dòng)熱預(yù)示模型對(duì)于非均布壁面溫度邊界條件的適用能力主要來(lái)自該模型的當(dāng)?shù)匦浴．?dāng)?shù)匦阅Ｐ湍茌^為輕松地考慮非均布壁面溫度等當(dāng)?shù)鼗吔鐥l件，而POD 等方法建立的全局性模型對(duì)于當(dāng)?shù)鼗吔鐥l件的處理則比較困難。

表4 典型位置熱流預(yù)示值和RANS 計(jì)算值相對(duì)誤差Table 4 Relative error of heat flux at typical locations of the two-stage compression surface between the BSLM prediction and the RANS calculation

圖11 熱流預(yù)示結(jié)果和RANS 計(jì)算結(jié)果對(duì)比Fig. 11 Comparison of heat flux between the RANS calculation and the BSLM prediction

2.4 三維雙橢球熱流預(yù)示

本節(jié)使用三維雙橢球算例測(cè)試BSLM 方法對(duì)三維構(gòu)形的適用性。

以三維雙橢球迎角0°、10°、20°和25°作為訓(xùn)練狀態(tài)，迎角15°作為測(cè)試狀態(tài)，選擇雙橢球表面的20441個(gè)數(shù)值模擬網(wǎng)格節(jié)點(diǎn)的熱流數(shù)據(jù)和邊界層外緣數(shù)據(jù)進(jìn)行訓(xùn)練和預(yù)示。選擇的特征組合為特征θ=tan?1(y/z)，x、y和z為空間坐標(biāo)。預(yù)示結(jié)果如圖12 所示，圖12（a）為RANS 計(jì)算的三維雙橢球熱流分布云圖，圖12（b）為BSLM 氣動(dòng)熱模型預(yù)示的三維雙橢球熱流分布云圖。可以看出，氣動(dòng)熱模型準(zhǔn)確辨識(shí)了大橢球頭部的高熱流區(qū)和小橢球頭部的次高熱流區(qū)，并且熱流分布與RANS 計(jì)算結(jié)果基本一致，而且RANS 計(jì)算和模型預(yù)示得到的駐點(diǎn)熱流分別為396.95 kW/m2和401.50 kW/m2，相對(duì)誤差僅為1.15%。由上述結(jié)果可以推斷，本文提出的BSLM 氣動(dòng)熱預(yù)示方法可以較為精確地實(shí)現(xiàn)三維構(gòu)形的氣動(dòng)熱預(yù)示。

圖12 RANS 計(jì)算和BSLM 預(yù)示的雙橢球表面熱流分布云圖Fig. 12 Heat flux contours on the double ellipsoid surface compared between the RANS calculation and the BSLM prediction

3 結(jié) 論

本文針對(duì)氣動(dòng)熱的高效、高精度預(yù)示問(wèn)題，在Euler 方程計(jì)算的邊界層外緣信息的基礎(chǔ)上，結(jié)合少量RANS 結(jié)果，使用支持向量機(jī)算法發(fā)展了一種數(shù)據(jù)驅(qū)動(dòng)的當(dāng)?shù)鼗瘹鈩?dòng)熱預(yù)示建模方法，實(shí)現(xiàn)了不同來(lái)流狀態(tài)、壁面溫度下二維及三維構(gòu)型的氣動(dòng)熱高效精準(zhǔn)預(yù)示。通過(guò)雙橢球上表面中心線熱流預(yù)示、二級(jí)壓縮面熱流預(yù)示以及三維雙橢球熱流預(yù)示結(jié)果驗(yàn)證，可以得到以下結(jié)論：

1）氣動(dòng)熱預(yù)示模型輸入特征選擇對(duì)氣動(dòng)熱預(yù)示精度有很大的影響，需要事先進(jìn)行合理的特征選擇以保證模型性能。本文所提出的特征選擇方法能夠篩選出合適的特征組合，使得模型具備較高精度和較好的泛化能力。

2）論文所提出的BSLM 氣動(dòng)熱預(yù)示方法具有較高的預(yù)示精度和良好的泛化性能，熱流峰值及典型位置的熱流預(yù)示結(jié)果和數(shù)值仿真結(jié)果的相對(duì)誤差均小于5%。

3）BSLM 方法與傳統(tǒng)POD 方法相比，可以充分利用邊界層外緣信息，具備更高的預(yù)示精度，特別是針對(duì)熱流峰值的預(yù)示精度可以提升4 倍以上。另外，由于采用當(dāng)?shù)鼗２呗?，?shí)現(xiàn)非均勻壁面溫度泛化方面，難度遠(yuǎn)低于POD 類全局性建模方法。