基于代理模型的高超聲速氣動(dòng)熱模型降階研究

陳鑫， 劉莉， 岳振江

(北京理工大學(xué) 宇航學(xué)院，飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081)

?

基于代理模型的高超聲速氣動(dòng)熱模型降階研究

陳鑫， 劉莉， 岳振江

(北京理工大學(xué) 宇航學(xué)院，飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081)

高超聲速飛行器氣動(dòng)熱的快速準(zhǔn)確預(yù)測(cè)是當(dāng)前高超聲速氣動(dòng)熱彈性分析的重要前提. 針對(duì)當(dāng)前高超聲速氣動(dòng)熱工程計(jì)算、高精度數(shù)值計(jì)算和實(shí)驗(yàn)研究均不能很好適應(yīng)工程應(yīng)用的問(wèn)題，結(jié)合代理模型的基本思想，提出了基于代理模型的高超聲速氣動(dòng)熱模型降階方法，建立了一種高超聲速氣動(dòng)熱模型降階框架. 以典型高超聲速三維翼面為例，對(duì)比拉丁超立方采樣方法lhsdesign函數(shù)和改進(jìn)的逐次枚舉的拉丁超立方方法SLE，利用相同的設(shè)計(jì)樣本點(diǎn)和代理模型構(gòu)造方法，SLE方法構(gòu)造的降階模型預(yù)測(cè)翼面溫度平均誤差、L∞和eNRSME均小于lhsdesign方法，SLE采樣方法有助于提高降階模型的精度；對(duì)比Kriging和RBF兩種代理模型構(gòu)造方法，Kriging方法構(gòu)造降階模型優(yōu)于RBF方法. 針對(duì)典型的高超聲速三維翼面氣動(dòng)熱預(yù)測(cè)表明，本文高超聲速氣動(dòng)熱降階方法具有較高的精度和效率.

高超聲速；氣動(dòng)熱；代理模型；模型降階

氣動(dòng)熱彈性分析涉及氣動(dòng)、結(jié)構(gòu)、控制和推進(jìn)等多個(gè)系統(tǒng)，在高超聲速飛行器的分析和優(yōu)化設(shè)計(jì)過(guò)程中越來(lái)越受到重視. 高超聲速飛行器氣動(dòng)熱的快速準(zhǔn)確預(yù)測(cè)是氣動(dòng)熱彈性分析的重要前提. 目前，主要利用氣動(dòng)加熱工程計(jì)算、計(jì)算流體力學(xué)(computational fluid dynamics， CFD)數(shù)值計(jì)算和實(shí)驗(yàn)研究3種方法來(lái)預(yù)測(cè)高超聲速氣動(dòng)熱環(huán)境. 基于簡(jiǎn)單幾何假設(shè)的氣動(dòng)加熱工程計(jì)算忽略了真實(shí)氣體效應(yīng)、氣流黏性等，僅適用于平板或圓柱等簡(jiǎn)單飛行器構(gòu)型的氣動(dòng)熱預(yù)測(cè)[1-2]. 已有文獻(xiàn)研究表明忽略真實(shí)氣體黏性和真實(shí)氣體效應(yīng)甚至?xí)苯佑绊懙椒治鼍萚3]. CFD數(shù)值計(jì)算能夠充分考慮氣流黏性、真實(shí)氣體效應(yīng)等，捕捉激波及邊界層轉(zhuǎn)捩等現(xiàn)象，能較好地求解耦合分析中廣泛存在的非線(xiàn)性方程，王衛(wèi)星等[4]研究指出邊界層轉(zhuǎn)捩對(duì)壁面熱流密度分布影響較大. 由于巨大的分析自由度以及由于不確定性和設(shè)計(jì)優(yōu)化需要的重復(fù)計(jì)算導(dǎo)致CFD計(jì)算量巨大[3]. 高超聲速氣動(dòng)熱的實(shí)驗(yàn)研究成本較高，同時(shí)實(shí)驗(yàn)研究大部分僅關(guān)注凸起等局部特性對(duì)流場(chǎng)邊界層的影響[5]，當(dāng)前應(yīng)用并不廣泛. 因此，開(kāi)展相應(yīng)地氣動(dòng)熱降階模型(reduced order model， ROM)研究已經(jīng)成為一個(gè)較為活躍的領(lǐng)域.

代理模型是通過(guò)數(shù)學(xué)方法，構(gòu)造出一個(gè)計(jì)算量較小，但計(jì)算結(jié)果與數(shù)值分析結(jié)果或真實(shí)物理試驗(yàn)結(jié)果相近的近似數(shù)學(xué)模型，以代替原數(shù)值分析模型或真實(shí)物理試驗(yàn). 代理模型以其較好的近似精度和魯棒性，近10年來(lái)被廣泛應(yīng)用于飛行器分析設(shè)計(jì)領(lǐng)域. 楊華等[6]利用徑向基函數(shù)，解決了復(fù)雜形狀機(jī)翼的二維氣動(dòng)力代理模型的構(gòu)造問(wèn)題. 梁煜等[7]基于計(jì)算流體力學(xué)分析結(jié)果構(gòu)建Kriging氣動(dòng)力代理模型，用于氣動(dòng)布局參數(shù)匹配優(yōu)化設(shè)計(jì)，提高設(shè)計(jì)效率并保證了可信度. 夏露等[8]提出了一種基于Kriging自適應(yīng)代理模型的氣動(dòng)優(yōu)化方法，成功應(yīng)用于翼型氣動(dòng)性能優(yōu)化設(shè)計(jì).

針對(duì)高超聲速氣動(dòng)熱的模型降階研究目前較少，本文結(jié)合代理模型基本思想，提出了一種基于代理模型的高超聲速氣動(dòng)熱的降階模型方法，建立了高超聲速氣動(dòng)熱模型降階框架. 典型三維翼面算例驗(yàn)證了降階模型的可行性，該降階模型有助于推動(dòng)CFD技術(shù)推廣到實(shí)際工程應(yīng)用.

1 代理模型

代理模型的基本思想是利用原始高精度模型或試驗(yàn)獲得一組樣本點(diǎn)，利用代理模型構(gòu)造方法從這組樣本點(diǎn)中提取輸入?yún)?shù)和輸出參數(shù)的近似關(guān)系進(jìn)而建立一個(gè)模擬原始高精度模型的模型[9]. 試驗(yàn)設(shè)計(jì)方法和代理模型構(gòu)造方法是代理模型的兩個(gè)重要組成部分.

1.1 試驗(yàn)設(shè)計(jì)方法

試驗(yàn)設(shè)計(jì)DoE(design of experiment)方法是一種科學(xué)合理的數(shù)學(xué)安排，在設(shè)計(jì)空間內(nèi)獲取能反映精確模型特征樣本點(diǎn)的方法. 試驗(yàn)設(shè)計(jì)的目標(biāo)是通過(guò)盡可能少的樣本點(diǎn)盡可能地反映空間內(nèi)精確模型的特征. 空間均勻性和投影均勻性是衡量試驗(yàn)設(shè)計(jì)方法的主要標(biāo)準(zhǔn). 常用的試驗(yàn)設(shè)計(jì)方法有均勻試驗(yàn)設(shè)計(jì)、中心復(fù)合設(shè)計(jì)和拉丁超立方設(shè)計(jì)等. 本文選取常用的Matlab平臺(tái)的拉丁超立方試驗(yàn)設(shè)計(jì)方法lhsdesign函數(shù)和基于Maxmin準(zhǔn)則逐次局部枚舉的改進(jìn)拉丁超立方試驗(yàn)設(shè)計(jì)方法[10](successive local enumeration， SLE).

1.1.1 拉丁超立方設(shè)計(jì)方法

拉丁超立方試驗(yàn)設(shè)計(jì)是一種分層抽樣方法，其試驗(yàn)設(shè)計(jì)點(diǎn)在設(shè)計(jì)空間內(nèi)分布較均勻，且試驗(yàn)次數(shù)等于水平數(shù). 拉丁超立方試驗(yàn)設(shè)計(jì)點(diǎn)的生成方法為

(1)

式中:1≤j≤n，1≤i≤m，n為因素個(gè)數(shù)或設(shè)計(jì)變量個(gè)數(shù)，m為樣本點(diǎn)個(gè)數(shù)；U為[0，1]區(qū)間內(nèi)的隨機(jī)數(shù)；i為第i次試驗(yàn)；j為第j個(gè)設(shè)計(jì)變量；π為0，1，…，m-1的獨(dú)立隨機(jī)排列.

標(biāo)準(zhǔn)的拉丁超立方設(shè)計(jì)方法僅具有投影均勻性的特點(diǎn)，為了獲得好的空間均勻性，學(xué)者開(kāi)始研究最優(yōu)拉丁超立方采樣方法(optimal latin hypercube design， OLHD)， 本文采用Matlab中l(wèi)hsdesign函數(shù)優(yōu)化采樣方法，lhsdesign函數(shù)采樣標(biāo)準(zhǔn)參數(shù)設(shè)置為“maximin”，迭代次數(shù)為100.

1.1.2 基于Maxmin準(zhǔn)則逐次局部枚舉的改進(jìn)拉丁超立方設(shè)計(jì)方法(SLE)

基于Maxmin準(zhǔn)則的逐次局部枚舉拉丁超立方試驗(yàn)設(shè)計(jì)方法SLE的采樣過(guò)程是一個(gè)簡(jiǎn)單的局部?jī)?yōu)化過(guò)程，目標(biāo)為待選定的樣本點(diǎn)與已生成的樣本點(diǎn)之間的最小距離最大化. 與其它現(xiàn)有拉丁超立方采樣算法最大的不同在于，其它算法采樣過(guò)程中有全局目標(biāo)函數(shù)，而SLE方法采樣過(guò)程中的目標(biāo)函數(shù)為局部的最小距離最大化.

對(duì)比研究表明SLE采樣方法具有很高的采樣效率和精度，均優(yōu)于lhsdesign，且不受采樣點(diǎn)個(gè)數(shù)和采樣空間維數(shù)的限制，具有很好的通用性[10].

1.2 代理模型構(gòu)造方法

常用的代理模型方法有多項(xiàng)式相應(yīng)面、移動(dòng)最小二乘法、徑向基函數(shù)、Kriging模型、神經(jīng)網(wǎng)絡(luò)等. 本文中分別采用徑向基函數(shù)(radial basis function， RBF)和Kriging模型.

1.2.1 徑向基函數(shù)

徑向基函數(shù)(RBF)是一類(lèi)以未知點(diǎn)與已知數(shù)據(jù)點(diǎn)之間的歐式距離作為自變量的函數(shù). 基本形式為

(2)

(3)

式中c=1/σ2，σ為高斯函數(shù)的寬度，本文取為2.0.

1.2.2Kriging函數(shù)

Kriging模型是對(duì)空間中分布的數(shù)據(jù)點(diǎn)進(jìn)行求線(xiàn)性最優(yōu)、無(wú)偏內(nèi)插估計(jì)的一種方法. 表達(dá)式為

(4)

式中:g(X)為設(shè)計(jì)空間范圍內(nèi)關(guān)于X的全局近似模型;z(X)為均值為0，方差為σ2，協(xié)方差不為0的隨機(jī)過(guò)程. 協(xié)方差矩陣為

(5)

(6)

引入相關(guān)向量，

(7)

Kriging模型可表示為

(8)

式中:β為未知參數(shù)；σ2和R都為θ的函數(shù)；y為由采樣點(diǎn)響應(yīng)值組成的ns維列向量.β和σ2由最小二乘估計(jì)得到

(9)

相關(guān)參數(shù)θ可通過(guò)優(yōu)化得到

(10)

Kriging代理模型需要復(fù)雜的全局性?xún)?yōu)化方法求解帶約束的優(yōu)化問(wèn)題. 本文中選遺傳算法求解θ.

2 基于代理模型的模型降階方法

代理模型用來(lái)建立輸入?yún)?shù)與單輸出參數(shù)或少量輸出參數(shù)的近似關(guān)系，例如利用代理模型建立機(jī)翼氣動(dòng)外形的幾何參數(shù)蒙皮厚度、翼梁厚度等與升力系數(shù)或阻力系數(shù)等氣動(dòng)參數(shù)的近似擬合關(guān)系. 通常僅知道簡(jiǎn)單的氣動(dòng)參數(shù)不能滿(mǎn)足實(shí)際工程設(shè)計(jì)要求，如高超聲速氣動(dòng)熱彈性分析和氣動(dòng)熱防護(hù)設(shè)計(jì)時(shí)需要得到整個(gè)翼面各處溫度分布情況. 因此如何準(zhǔn)確高效得到整個(gè)流場(chǎng)參數(shù)具有較大研究意義.

對(duì)于確定構(gòu)型的繞流流場(chǎng)，描述構(gòu)型流場(chǎng)的數(shù)量為構(gòu)型表面流場(chǎng)網(wǎng)格單元數(shù)與流場(chǎng)變量個(gè)數(shù)的乘積. 針對(duì)高超聲速飛行器氣動(dòng)熱預(yù)測(cè)，溫度為標(biāo)量，描述高超聲速飛行器氣動(dòng)熱流場(chǎng)的個(gè)數(shù)為飛行器表面流場(chǎng)網(wǎng)格單元數(shù).

如圖1所示，利用代理模型預(yù)測(cè)氣動(dòng)熱步驟如下.

步驟1 確定設(shè)計(jì)空間和設(shè)計(jì)變量. 針對(duì)高超聲速機(jī)翼氣動(dòng)熱預(yù)測(cè)，設(shè)計(jì)變量選為飛行馬赫數(shù)、飛行高度和飛行攻角，設(shè)計(jì)空間為飛行馬赫數(shù)、飛行高度和飛行攻角的上下限范圍.

步驟2 運(yùn)用試驗(yàn)設(shè)計(jì)方法獲得設(shè)計(jì)空間內(nèi)的樣本點(diǎn)I(i)，i=1，2，…，n，n為樣本點(diǎn)個(gè)數(shù)；樣本點(diǎn)在設(shè)計(jì)空間內(nèi)保證空間均勻性和空間正交性.

步驟4 利用代理模型方法構(gòu)造設(shè)計(jì)樣本點(diǎn)與翼面上的每一個(gè)節(jié)點(diǎn)溫度之間的近似擬合關(guān)系. 翼面共有p個(gè)節(jié)點(diǎn)，即建立p個(gè)近似擬合關(guān)系.

3 高超聲速氣動(dòng)熱模型降階框架

針對(duì)高超聲速飛行器氣動(dòng)熱計(jì)算遇到的復(fù)雜耗時(shí)的非線(xiàn)性數(shù)值求解難題，基于代理模型提出一種快速高效的高超聲速氣動(dòng)熱模型降階(ROM)框架，如圖2所示.

降階模型框架主要包括:物理模型的選擇；設(shè)計(jì)變量及設(shè)計(jì)空間的確定；高精度數(shù)值計(jì)算；surrogate方法構(gòu)造降階模型；降階模型近似精度評(píng)估.

4 算例驗(yàn)證及分析

4.1 模型描述及設(shè)計(jì)空間的確定

本文中選取氣動(dòng)熱結(jié)構(gòu)分析中常用的典型的F104戰(zhàn)斗機(jī)的機(jī)翼，該機(jī)翼被多次運(yùn)用來(lái)做典型的高超聲速分析源于其有良好的工程應(yīng)用背景[11]. 機(jī)翼相關(guān)參數(shù)如圖3所示.

F104機(jī)翼的氣動(dòng)熱設(shè)計(jì)變量和設(shè)計(jì)空間如下. 馬赫數(shù)：5.0≤Ma≤10.0;攻角:-8.0°≤α≤8.0°;飛行高度:20 km≤H≤40 km.

初始選取100個(gè)樣本點(diǎn)，采用拉丁超立方試驗(yàn)設(shè)計(jì)方法Matlab的lhsdesign函數(shù)以及基于Maxmin準(zhǔn)則的逐次局部枚舉改進(jìn)的拉丁超立方(SLE)方法分別產(chǎn)生.

4.2 高超聲速氣動(dòng)熱數(shù)值計(jì)算

本文高超聲速氣動(dòng)熱計(jì)算采用CFD-Fastran求解器計(jì)算F104翼面的溫度分布和熱流分布.

Mcnamra J J[11]指出在計(jì)算高超聲速氣動(dòng)熱建議采用B-L湍流模型，B-L湍流模型是一種最基本的湍流模型，能更好地反應(yīng)高超聲速流中的非線(xiàn)性特征. 董素君等[12]利用CFD-Fastran求解器針對(duì)典型雙錐面得到不同湍流模型計(jì)算高超聲速氣動(dòng)熱流，指出B-L湍流模型計(jì)算快，能達(dá)到較好的準(zhǔn)確性，建議氣動(dòng)熱計(jì)算采用B-L模型，同時(shí)指出近壁面第一層網(wǎng)格間距y-plus需保證在1以下. 閻超等[13]通過(guò)探討了CFD熱流計(jì)算的格式效應(yīng)及網(wǎng)格效應(yīng)，研究指出AUSM格式在熱流計(jì)算精確性方面具有優(yōu)勢(shì). 潘沙等[14]和覃文潔等[15]指出氣動(dòng)熱的數(shù)值模擬中，網(wǎng)格因素非常關(guān)鍵，計(jì)算結(jié)果對(duì)網(wǎng)格，特別是網(wǎng)格壁面附近的法向網(wǎng)格間距十分敏感，同時(shí)指出氣動(dòng)熱計(jì)算收斂的判斷比壓力收斂緩慢的多，建議采用直接觀察氣動(dòng)熱數(shù)據(jù)來(lái)判斷收斂.

因此，本文計(jì)算采用B-L湍流模型，空間格式采用高階AUSM，并且嚴(yán)格保證近壁面網(wǎng)格的大小，計(jì)算時(shí)保證氣動(dòng)熱收據(jù)的收斂性.

董素君等[12]通過(guò)典型雙錐面得到不同湍流模型計(jì)算高超聲速氣動(dòng)熱流，并與試驗(yàn)數(shù)據(jù)對(duì)比指出了Fastran求解器計(jì)算高超聲速飛行器氣動(dòng)熱的可行性. 如圖4所示，進(jìn)一步通過(guò)與文獻(xiàn)[11]中CFL3D求解器對(duì)比得到，文獻(xiàn)[11]中CFL3D求解器翼面前緣駐點(diǎn)溫度為938.0 K，F(xiàn)astran計(jì)算駐點(diǎn)溫度為988.3 K，翼面上溫度梯度分布趨勢(shì)相近，翼面溫度分布情況總體相差并不大，因此本文均采用CFD-Fastran求解器計(jì)算F104翼面的溫度分布和熱流分布.

本文采用計(jì)算流體網(wǎng)格如圖5所示，分別針對(duì)42.0萬(wàn)、81.9萬(wàn)和158.0萬(wàn)的計(jì)算域進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，如表1所示，當(dāng)網(wǎng)格計(jì)算量為42.0萬(wàn)時(shí)，平均y-plus為3.82，駐點(diǎn)溫度為1 078.6 K. 隨著網(wǎng)格計(jì)算量增大到81.9萬(wàn)，近壁面網(wǎng)格平均y-plus減小為0.76，駐點(diǎn)溫度為988.3 K. 進(jìn)一步增加網(wǎng)格數(shù)量，平均y-plus減小到0.45，駐點(diǎn)溫度相比于81.9萬(wàn)網(wǎng)格并未明顯減小，為987.9 K. 網(wǎng)格無(wú)關(guān)性驗(yàn)證計(jì)算結(jié)果得到，計(jì)算樣本點(diǎn)均可采用選取81.9萬(wàn)網(wǎng)格作為計(jì)算網(wǎng)格.

Tab.1 Leading edge temperature with different number of grids

網(wǎng)格計(jì)算量/萬(wàn)駐點(diǎn)溫度/K平均y?plus42010786382819988307615809879045

4.3 算例計(jì)算結(jié)果及分析

Culler A J等[16]指出氣動(dòng)熱彈性分析時(shí)氣動(dòng)彈性系統(tǒng)響應(yīng)與熱傳導(dǎo)之間巨大的時(shí)間尺度差別，故在氣動(dòng)熱彈性分析中計(jì)算氣動(dòng)熱時(shí)可以只考慮穩(wěn)態(tài)氣動(dòng)熱分布. 因此，本文在利用CFD計(jì)算氣動(dòng)熱時(shí)僅僅考慮機(jī)翼的穩(wěn)態(tài)氣動(dòng)熱流分布情況. 所有翼型的網(wǎng)格拓?fù)浣Y(jié)構(gòu)和網(wǎng)格節(jié)點(diǎn)數(shù)均相同，為81.9萬(wàn)結(jié)構(gòu)化網(wǎng)格，物面邊界條件設(shè)置為輻射熱平衡，輻射系數(shù)0.85.

為了定量衡量預(yù)測(cè)結(jié)果的好壞，采用均方根誤差eNRMSE和最大值誤差L∞，其公式為

(11)

(12)

式中：i為預(yù)測(cè)工況下的第i個(gè)節(jié)點(diǎn)的溫度值；TROM為在Surrogate方法得到的降階模型的預(yù)測(cè)值；TFull為CFD計(jì)算得到的溫度值.

4.3.1 測(cè)試樣本點(diǎn)相對(duì)平均誤差分析

利用拉丁超立方設(shè)計(jì)方法lhsdesign函數(shù)在設(shè)計(jì)空間內(nèi)隨機(jī)選取10個(gè)測(cè)試樣本點(diǎn). 利用高超聲速氣動(dòng)熱模型降階框架，針對(duì)10個(gè)測(cè)試樣本點(diǎn)，翼面上表面溫度平均相對(duì)誤差如圖6所示.

如圖6所示，圖6(a)、6(b)為利用lhsdesign函數(shù)選取樣本點(diǎn)，分別選取Kriging和RBF代理模型構(gòu)造方法的到的10個(gè)預(yù)測(cè)點(diǎn)翼面上表面平均相對(duì)誤差圖，對(duì)比圖6(a)和6(b)可以看到：利用相同的采樣方法，Kriging方法得到翼面上表面溫度平均相對(duì)誤差最大為5%，而RBF方法翼面上表面溫度平均相對(duì)誤差最大達(dá)到14%，對(duì)比圖6(c)及6(d)同樣發(fā)現(xiàn)利用相同的SLE采樣方法，Kriging方法得到翼面上表面溫度平均相對(duì)誤差最大為1%，而RBF方法翼面上表面溫度平均相對(duì)誤差最大達(dá)到20%，在處理高超聲速氣動(dòng)熱預(yù)測(cè)時(shí)kriging方法精度要優(yōu)于RBF方法. 同時(shí)對(duì)比圖6(a)6(c)可知，取相同設(shè)計(jì)樣本點(diǎn)，具有較好空間均勻性和投影均勻性的SLE方法得到翼面上表面溫度平均相對(duì)誤差均要小于lhsdesign方法，圖6(c)的SLE采樣方法翼面最大平均相對(duì)誤差為1%，而圖6(a)中l(wèi)hsdesign方法最大平均相對(duì)誤差達(dá)到5%，圖6(c)中翼面各處翼面溫度平均相對(duì)誤差均小于圖6(a)中翼面各處翼面溫度平均相對(duì)誤差.

4.3.2 測(cè)試樣本點(diǎn)eNRSME和L∞誤差分析

表2給出了不同采樣方法(lhsdesign/SLE)和不同代理模型構(gòu)造方法(Kriging/RBF)的10個(gè)測(cè)試點(diǎn)上、下翼面預(yù)測(cè)溫度的eNRMSE和L∞誤差分布情況，Avg為測(cè)試樣本點(diǎn)的平均誤差.

由表2對(duì)比可知，利用相同lhsdesign試驗(yàn)設(shè)計(jì)方法，Kriging方法的上下翼面預(yù)測(cè)溫度eNRSME誤差均小于RBF方法，上表面eNRSME平均誤差Kriging方法為2.6%，小于RBF方法的7.6%，下表面eNRSME平均誤差Kriging方法為3.1%，而RBF方法的eNRSME平均誤差為9.6%；Kriging方法上表面翼面溫度預(yù)測(cè)平均誤差L∞為5.0%，而RBF方法上表面翼面溫度預(yù)測(cè)平均誤差L∞為10.2%，是Kriging方法的2倍. Kriging方法下表面翼面溫度預(yù)測(cè)平均誤差L∞為5.6%，同樣小于RBF方法的溫度平均誤差7.7%. 同樣對(duì)比分析得，相同的SLE試驗(yàn)設(shè)計(jì)方法，得到Kriging方法的上、下表面翼面溫度預(yù)測(cè)eNRSME平均誤差分別為0.85%,0.95%，遠(yuǎn)小于RBF方法的11.1%,12.6%，L∞平均誤差1.9%,2.6%，小于RBF方法的15.6%,8.2%. 分析可知，Kriging方法在處理高超聲速氣動(dòng)熱預(yù)測(cè)時(shí)比RBF方法更具有優(yōu)勢(shì).

表2 不同采樣方法和不同代理模型構(gòu)造方法測(cè)試點(diǎn)eNRMSE和L∞

Tab.2eNRSMEandL∞of test cases by different sampling methods and surrogate methods

采樣方法lhsde?signSLE替代方法KrigingRBFKrigingRBF翼面eNRMSE/%L∞/%上029～960(Avg．26)046～1870(Avg．50)下034～1540(Avg．31)057～2440(Avg．56)上11～143(Avg．76)17～232(Avg．102)下08～216(Avg．96)07～271(Avg．77)上017～420(Avg．085)063～580(Avg．19)下017～390(Avg．095)08～62(Avg．26)上35～195(Avg．111)37～305(Avg．156)下31～238(Avg．126)066～2960(Avg．82)

對(duì)比表2中的lhsdesign試驗(yàn)設(shè)計(jì)方法和SLE試驗(yàn)設(shè)計(jì)方法，利用相同的100個(gè)試驗(yàn)設(shè)計(jì)樣本點(diǎn)，lhsdesign試驗(yàn)設(shè)計(jì)方法選取樣本點(diǎn)，利用Kriging方法上下表面溫度預(yù)測(cè)的eNRSME平均誤差分別為2.6%和3.1%，L∞平均誤差分別為5.0%和5.6%；SLE試驗(yàn)設(shè)計(jì)方法選取樣本點(diǎn)，同樣利用Kriging方法上下表面溫度預(yù)測(cè)的eNRSME平均誤差分別為0.85%和0.95%，L∞平均誤差分別為1.9%和2.6%. 顯然利用相同的試驗(yàn)設(shè)計(jì)方法，SLE試驗(yàn)設(shè)計(jì)方法方法能夠明顯提高降階模型的預(yù)測(cè)精度.

4.3.3 典型算例結(jié)果與分析

選取典型工況下的上、下翼面溫度CFD值與降階模型(ROM)預(yù)測(cè)結(jié)果如圖7所示，典型工況的飛行速度為1 798.1 m/s，高度為36 612.0 m，攻角為-1.0°. 如圖7所示，選取SLE采樣樣本點(diǎn)分別利用Kriging方法和RBF方法構(gòu)造翼面上下表面溫度降階模型. 通過(guò)ROM與對(duì)應(yīng)真實(shí)值對(duì)比可知，kriging方法構(gòu)造翼面上下表面如圖7(a)和7(c)降階模型預(yù)測(cè)結(jié)果均好于RBF方法，如圖7(b)和7(d).

典型算例中Kriging方法構(gòu)造高超聲速氣動(dòng)熱降階模型預(yù)測(cè)結(jié)果與CFD計(jì)算溫度值符合，降階模型較好，具有較高的精度，驗(yàn)證了本文高超聲速氣動(dòng)熱降階模型的可行性.

4.3.4 降階模型效率分析

如表3所示，翼面上下表面各有4 800個(gè)節(jié)點(diǎn)，即利用代理模型方法構(gòu)造4 800個(gè)近似擬合關(guān)系才能預(yù)測(cè)完整翼面溫度分布. 利用數(shù)值計(jì)算(CFD)完成一個(gè)樣本點(diǎn)的計(jì)算需約19.1 h. 選取100個(gè)設(shè)計(jì)樣本點(diǎn)，利用Kriging方法構(gòu)造高超聲速氣動(dòng)熱降階模型耗時(shí)約45.94 s，而RBF耗時(shí)僅僅為0.009 s，這是由于Kriging方法構(gòu)造近似關(guān)系是需要全局優(yōu)化過(guò)程，這樣比徑向基函數(shù)(RBF)更耗時(shí)，構(gòu)造4 800個(gè)近似擬合關(guān)系耗時(shí)遠(yuǎn)遠(yuǎn)大于RBF.

表3 降階模型與數(shù)值計(jì)算對(duì)比

注：在i5-2320上計(jì)算，3.0 GHz， 4.00 GB RAM.

針對(duì)單個(gè)預(yù)測(cè)工況，利用Kriging方法的高超聲速氣動(dòng)熱降階模型(ROM)預(yù)測(cè)耗時(shí)1.17 s，RBF方法構(gòu)造的降階模型(ROM)耗時(shí)1.03 s，兩種方法構(gòu)造的ROM預(yù)測(cè)耗時(shí)相近，與數(shù)值計(jì)算耗時(shí)約69 760 s(19.1 h)相比，本文建立的高超聲速氣動(dòng)熱降階模型(ROM)具有很高的效率.

5 結(jié) 論

① 提出一種基于代理模型的運(yùn)用于高超聲速氣動(dòng)熱預(yù)測(cè)的模型降階方法，建立了一種快速高效的高超聲速氣動(dòng)熱降階模型框架. 典型三維算例表明，該模型降階方法和降階模型框架能夠成功應(yīng)用于高超聲速氣動(dòng)熱計(jì)算，具有較高的精度和效率；

② 對(duì)比lhsdesign計(jì)算試驗(yàn)方法和SLE計(jì)算試驗(yàn)方法構(gòu)造高超聲速氣動(dòng)熱ROM，指出相同的設(shè)計(jì)樣本點(diǎn)，SLE計(jì)算試驗(yàn)方法的氣動(dòng)熱ROM具有更高的精度，空間均勻性和投影均勻性更好的SLE試驗(yàn)設(shè)計(jì)方法有助于提高氣動(dòng)熱ROM的預(yù)測(cè)精度；對(duì)比Kriging方法和RBF方法精度和效率，指出Kriging方法在構(gòu)造高超聲速氣動(dòng)熱降階模型時(shí)更具有優(yōu)勢(shì).

[1] 陳鑫，劉莉，李昱霖，等.高超聲速飛行器翼面氣動(dòng)加熱的工程計(jì)算方法[J].彈箭與制導(dǎo)學(xué)報(bào)，2013，33(3):133-137.

Chen Xin， Liu Li， Li Yunlin， et al. Engineering calculation of aerodynamic heating for airfoils of hypersonic vehicles[J]. Journal of Projectiles， Rockets， Missiles and Guidance， 2013，33(3)：133-137. (in Chinese)

[2] 車(chē)競(jìng)，唐碩，何開(kāi)鋒.類(lèi)乘波體飛行器氣動(dòng)加熱的工程計(jì)算方法[J]. 彈道學(xué)報(bào)，2006，18(4):58-65.

Che Jing， Tang Shuo， He Kaifeng. Engineering calculation of aerodynamic heating for quasi-waverider vehicle[J]. Journal of Ballistics， 2006， 18(4):58-65. (in Chinese)

[3] Mcnamara J J， Crowell A R. Approximate modeling of unsteady aerodynamics for hypersonic aeroelastic[J]. Journal of Aircraft， 2010，47(6):1932-1945.

[4] 王衛(wèi)星，郭榮偉.基于邊界層轉(zhuǎn)捩的高超聲速進(jìn)氣道特性研究[J].航空學(xué)報(bào)，2012，33(10):1772-1780.

Wang Weixing， Guo Rongwei. Study of flow characteristicsof hypersonic inlet based on boundary layer transition[J]. Acta Aeronautica et Astronautica Sinca, 2012，33(10):1772-1780. (in Chinese)

[5] Estruch S D，卜雪琴.高超聲速下表面凸起干擾氣動(dòng)熱實(shí)驗(yàn)研究[J].航空學(xué)報(bào)，2012，33(9):1578-1586.

Estruch S D， Bu Xueqin. Experimental investigation on hypersonic interface heating around surface protuberance[J]. Acta Aeronautica et Astronautica Sinca, 2012，33(9):1578-1586. (in Chinese)

[6] 楊華，姚衛(wèi)星.基于徑向基函數(shù)的機(jī)翼二維氣動(dòng)代理模型設(shè)計(jì)[J].計(jì)算力學(xué)學(xué)報(bào)，2008，25(6):797-802.

Yang Hua， Yao Weixing. 2D surrogate model of wing lift distribution based on radial function[J]. Chinese Journal of Compuational Mechanics， 2008，25(6):797-802. (in Chinese)

[7] 梁煜，程小全，酈正能，等.基于代理模型的氣動(dòng)外形平面參數(shù)多目標(biāo)匹配設(shè)計(jì)[J].航空學(xué)報(bào)，2010，31(6):1141-1148.

Liang Yu， Cheng Xiaoquan， Li Zhengneng， et al. Multi-object aerodynamic configuration parameter design using kriging approximation[J]. Acta Aeronautica et Astronautica Sinca， 2010，31(6):1141-1148. (in Chinese)

[8] 夏露，王丹.基于Kriging自適應(yīng)代理模型的氣動(dòng)優(yōu)化方法[J].航空計(jì)算技術(shù)，2013，43(1):13-17.

Xia Lu， Wang Dan. Aerodynamic optimization method based on Kriging adaptive surrogate model[J]. Aeronautical Computing Technique， 2013，43(1):13-17. (in Chinese)

[9] Simpson T W， Booker A J， Ghosh D， et al. Approximation methods in multidisciplinary analysis and optimization: a panel discussion[J]. Structural and Multidisplinary Optimization， 2004，27(5):303-313.

[10] Zhu H G， Liu L， Long T， et al. A novel algorithm of maximin Latin Hypercube design using successive local emumeration[J]. Engineering Optimization， 2012，44(5):551-564.

[11] Mcnamara J J. Aeroelastic and aerothermoelastic behavior of two and three dimensional lifting surfaces in hypersonic flow[D]. Ann Arbor: University of Michigan， 2005.

[12] 董素君，居世超，齊玢，等.CFD-Fastran氣動(dòng)熱計(jì)算模型及網(wǎng)格效應(yīng)分析[J].航空計(jì)算技術(shù)，2011，41(2):40-42.

Dong Sujun， Ju Shichao， Qi Fen. et al. Model and grid dependency accuracy by CFD-FASTRAN software[J]. Aeronautical Computing Technique， 2011，41(2):40-42. (in Chinese)

[13] 閻超，禹建軍，李君哲.熱流CFD計(jì)算中格式和網(wǎng)格效應(yīng)若干問(wèn)題研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2006，24(1):125-130.

Yan Chao， Yu Jianjun， Li Junzhe. Scheme effect and grid dependency in CFD computations of heat transfer[J]. Acta Aerodynamica Sinica, 2006，24(1):125-130. (in Chinese)

[14] 潘沙，馮定華，丁國(guó)昊，等.氣動(dòng)熱數(shù)值模擬中的網(wǎng)格相關(guān)性及收斂[J].航空學(xué)報(bào)，2010，31(3):493-499.

Pan Sha， Feng Dinghua， Ding Guohao， et al. Grid dependency and convergence of hypersonic arothermal simulation[J]. Acta Aeronautica Et Astronautica Sinca， 2010，31(3)：493-499. (in Chinese)

[15] 覃文潔，胡春光，郭良平，等.近壁面網(wǎng)格尺寸對(duì)湍流計(jì)算的影響[J].北京理工大學(xué)學(xué)報(bào)，2006，26(5):388-392.

Qin Wenjie， Hu Chunguang， Guo Liangping， et al. Effect of near-wall grid size on turbulent flow solutions[J]. Transactions of Beijing Institude of Technology， 2006，26(5):388-392. (in Chinese)

[16] Culler A J， Mc. Namara J J. Studies on fluid-thermal-structural coupling for aerothermoelasticity in hypersonic flow[J]. AIAA Journal， 2010，48(5):1721-1738.

(責(zé)任編輯：劉雨)

A Reduced Order Modeling for Aerothermodynamic of Hypersonic Vehicles Based on Surrogate Method

CHEN Xin， LIU Li， YUE Zhen-jiang

(Key Laboratory of Dynamics and Control of Flight Vehicle， Ministry of Education，School of Aerospace Engineering， Beijing Institute of Technology， Beijing 100081， China)

Accurate and efficient estimation of aerodynamic heating is the basic prerequisite of theaerothermo-elasticity analysis. In order to solve the engineering adaptability problems of hypersonic aerothermodynamics calculation, numerical superprecision calculation and experimental investigation, a novel estimation method based on surrogate method for aerothermdynamic was proposed. Furthermore, a reduced order modeling (ROM) framework for aerothermodynamic was developed. Test results for the three-dimensional aerothermodynamics over a typical hypersonic control surface indicate that the average absolute errors,L∞and NRSME by SLE are all smaller than those by lhsdesign using the same samlping points and the same sorrogate methods. So the SLE sampling method can help to improve the precision of ROM. Comparing the results of ROMs for the three-dimensional aerothermodynamics over a hypersonic control surface, it is indicated that the precision of Kriging is better than that of RBF. The developed ROMs for the three-dimensional surface show higher precision and efficiency for hypersonic aerothermodynamic.

hypersonic; aerothermodynamic; surrogate model; reduced oder model

2014-06-11

國(guó)家自然科學(xué)基金資助項(xiàng)目(11372036)

陳鑫(1988—)，男，博士生，E-mail:blingkx@hotmail.com.

劉莉(1964—)，女，教授，博士生導(dǎo)師，E-mail:liuli@bit.edu.cn.

V 215.3； V 215.4

A

1001-0645(2016)04-0340-08

10.15918/j.tbit1001-0645.2016.04.002