王 洋，袁軍婭，王洪興

（1. 北京航空航天大學(xué)，北京，100191；2. 航空工業(yè)第一飛機(jī)設(shè)計(jì)研究院，西安，710089；3. 北京衛(wèi)星環(huán)境工程研究所，北京，100029）

0 引 言

高超聲速飛行器一般是指飛行馬赫數(shù)大于 5的飛行器，飛行器與大氣劇烈摩擦，產(chǎn)生強(qiáng)烈的氣動(dòng)加熱。同時(shí)與結(jié)構(gòu)之間發(fā)生復(fù)雜的耦合作用。在流固耦合分析中，表面熱流是流固邊界傳遞的關(guān)鍵參數(shù)。高超聲速氣動(dòng)熱與來流條件（飛行高度、馬赫數(shù)、攻角等）、表面溫度和結(jié)構(gòu)形狀密切相關(guān)。目前，在高超聲速流固耦合研究中，氣動(dòng)加熱的計(jì)算方法主要分為工程計(jì)算和計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）數(shù)值計(jì)算。工程算法基于邊界層理論或者實(shí)驗(yàn)經(jīng)驗(yàn)關(guān)系式，主要用于求解簡單外形飛行器的氣動(dòng)熱問題，對于復(fù)雜外形，由于其流場復(fù)雜，存在流動(dòng)分離、激波-附面層干擾等現(xiàn)象，單純的工程方法難以滿足要求[1,2]。隨著計(jì)算機(jī)和CFD數(shù)值計(jì)算方法的不斷發(fā)展，已經(jīng)能夠針對復(fù)雜外形開展全尺寸三維流場仿真計(jì)算，但是效率仍舊較低，不能滿足快速工程研制的需求。CFD降階模型的研究成為解決這類高超聲速氣動(dòng)熱問題的一種有效途徑，可以克服工程方法不夠精確和CFD計(jì)算耗時(shí)的問題。

本文從兩個(gè)方面著手，減少流固耦合分析中熱邊界的計(jì)算耗時(shí)。一方面應(yīng)用本征正交分解（Proper Orthogonal Decomposition，POD）降階代理模型，建立一種快速計(jì)算高超聲速氣動(dòng)熱的方法。另一方面引入基于換熱系數(shù)的線性近似方法進(jìn)行熱邊界傳遞，通過減少因壁溫引起的流場計(jì)算迭代次數(shù)來縮短流固耦合計(jì)算時(shí)間。以三維翼面為例，驗(yàn)證了該方法有效性。

1 本征正交分解與代理模型

降階就是利用一個(gè)簡單的模型來代替CFD的來流條件和計(jì)算結(jié)果之間復(fù)雜的非線性關(guān)系，其中POD方法已經(jīng)被很多研究人員采納，并得到一定的成果。Burkardt等人[3]用該方法得到了N-S方程的降階模型；Carlson和Roveda[4]將弱化形式的動(dòng)量守恒方程投影到由速度相關(guān)的POD基函數(shù)形成的降階空間，建立了一組低維的POD系數(shù)為狀態(tài)變量的微分方程組來預(yù)測和控制翼面上的氣動(dòng)力，對 CFD和降階模型（Reduced Order Model，ROM）的數(shù)據(jù)結(jié)果進(jìn)行對比研究指出ROM利用少量的POD基模態(tài)能夠準(zhǔn)確地反映系統(tǒng)的動(dòng)態(tài)特征。除此之外，POD方法還被應(yīng)用在最優(yōu)控制[5～7]、流體力學(xué)[8～10]、熱傳導(dǎo)[11,12]分析中，極大地縮短了計(jì)算時(shí)間，提高計(jì)算效率。

POD方法是一種將復(fù)雜系統(tǒng)分解為數(shù)個(gè)基本模態(tài)的方法，在基本模態(tài)中包含著原始系統(tǒng)的主要特征，利用少量的基本模態(tài)就可以表達(dá)復(fù)雜的原始系統(tǒng)的特征，得到原始系統(tǒng)的ROM。利用實(shí)驗(yàn)或者CFD數(shù)值計(jì)算得到一系列不同輸入條件下的響應(yīng)輸出采用列向量的形式構(gòu)成向量組，通過矩陣運(yùn)算得到一組規(guī)范正交基，使得數(shù)據(jù)集合在這組基上的投影最大。通過對矩陣進(jìn)行奇異值分解就可以得到正交基Φ以及特征值，由于矩陣S并非方陣，所以得到的特征值應(yīng)該是方陣的特征值。特征值的大小表征了該特征值對應(yīng)的特征向量，即POD基所包含的響應(yīng)輸出集合的特征的多少。通常通過廣義“能”E來表征[13]前n個(gè)POD基所包含原系統(tǒng)特征的權(quán)重，廣義能定義為

式中 λ為特征值；n為POD基的個(gè)數(shù)；k為系統(tǒng)樣本點(diǎn)個(gè)數(shù)。

通常特征值從大到小衰減很快，這樣廣義能 E會(huì)趨近于1，前面幾個(gè)POD基就能包含絕大部分的集合特征。通過對特征向量的個(gè)數(shù)進(jìn)行截?cái)?，取前面部分基來表達(dá)整個(gè)系統(tǒng)，從而將全階系統(tǒng)投影到低維空間。

代理模型的基本思想是通過數(shù)學(xué)方法，利用已有的輸入和輸出參數(shù)，構(gòu)造出一個(gè)和原始計(jì)算結(jié)果相近的數(shù)學(xué)模型[14]。目前常用代理模型有多項(xiàng)式響應(yīng)面模型、Kriging模型和徑向基函數(shù)等，本文采用 Kriging模型。

式中 f(X)為回歸模型，代表確定的函數(shù)； z(X)為一個(gè)均值為0，方差為σ2的隨機(jī)過程，其協(xié)方差為

式中 R為兩個(gè)輸出參量 X(i)和 X(j)的相關(guān)函數(shù)，與兩點(diǎn)在空間的位置相關(guān)，關(guān)系式為

其中，

式中 g為所有元素為1的向量。這樣Kriging模型可以表示為

給定一個(gè)新的輸入?yún)?shù)向量X，即可通過該模型預(yù)測輸出參數(shù)。因此只需要將帶約束的優(yōu)化問題解決，得到θ，就可以構(gòu)建出此模型來。一般通過遺傳算法來求解最優(yōu)的θ向量。

POD方法通過利用全階模型計(jì)算結(jié)果建立一組最佳充分描述全階系統(tǒng)動(dòng)力學(xué)特性的正交基來描述系統(tǒng)的主要特性[16]。代理模型的目的是得到原系統(tǒng)樣本點(diǎn)與截?cái)嗪驪OD基投影系數(shù)的一一對應(yīng)關(guān)系。本文以高超聲速飛行器機(jī)翼氣動(dòng)熱計(jì)算為例，建立POD-Kriging代理模型，過程包括：

a）選取馬赫數(shù)、飛行高度和飛行攻角為設(shè)計(jì)變量，根據(jù)飛行軌道設(shè)計(jì)變量變化范圍，確定設(shè)計(jì)變量空間，采用拉丁超立方試驗(yàn)設(shè)計(jì)方法在給定的變量范圍內(nèi)獲得設(shè)計(jì)空間樣本I=[Ma,α,H]。

b）采用數(shù)值計(jì)算方法得到各樣本點(diǎn)對應(yīng)響應(yīng)輸出，以表面熱流為例，輸出即為飛行器表面所有網(wǎng)格點(diǎn)給定壁溫條件下的熱流q，構(gòu)建系統(tǒng)的特征矩陣

式中 n為樣本點(diǎn)個(gè)數(shù)；m為網(wǎng)格點(diǎn)個(gè)數(shù)。

c）將系統(tǒng)矩陣Q進(jìn)行奇異值分解，利用公式[15]：

對于給定樣本空間范圍內(nèi)的任意計(jì)算狀態(tài)，利用這個(gè)近似擬合關(guān)系得到POD基系數(shù)，從而得到q的近似響應(yīng)值。

2 模型以及算例分析

2.1 模型和設(shè)計(jì)空間的選取

本文以戰(zhàn)斗機(jī) F-104機(jī)翼模型為例，分析基于POD-Kriging代理模型獲得高超聲速氣動(dòng)熱的過程、精度和計(jì)算效率。該模型常被用來做高超聲速分析，機(jī)翼相關(guān)參數(shù)如圖1所示。

圖1 機(jī)翼外形Fig.1 Shape of the Wing c—弦長

機(jī)翼外流場模型采用ICEM網(wǎng)格繪制軟件劃分，分為8個(gè)block，對機(jī)翼附近網(wǎng)格進(jìn)行附面層加密，網(wǎng)格如圖2所示，網(wǎng)格數(shù)為197.6萬。

圖2 機(jī)翼外流場網(wǎng)格Fig.2 Wing Outflow Field Grid

機(jī)翼的氣動(dòng)熱設(shè)計(jì)變量和設(shè)計(jì)空間如表1所示。

表1 機(jī)翼氣動(dòng)熱設(shè)計(jì)變量和設(shè)計(jì)空間Tab.1 Aerodynamic Heat Design Variables and Space of Wing

2.2 POD-Kriging代理模型的建立

在給定的設(shè)計(jì)空間內(nèi)，經(jīng)過拉丁超立方抽樣，得到的100個(gè)樣本點(diǎn)的初始來流條件，如表2所示。

表2 樣本點(diǎn)數(shù)據(jù)Tab.2 Sample Point Data

對于每一個(gè)樣本點(diǎn)，本文采用CFD-Fastran軟件計(jì)算翼面熱流分布。CFD-Fastran常用于高超聲速飛行器外流場計(jì)算，文獻(xiàn)[17]通過數(shù)值計(jì)算與試驗(yàn)對比，得出結(jié)論CFD-FASTRAN軟件的k-ε及B-L兩種模型求解高超聲速氣動(dòng)熱問題均能達(dá)到較好的準(zhǔn)確度，且 B-L模型計(jì)算速度快、使用方便。本文采用 B-L湍流模型，空間采用高階AUSM格式，時(shí)間采用隱式迭代格式。以計(jì)算得到的不同來流條件下恒壁溫的熱流分布Q為基礎(chǔ)，建立POD-Kringing代理模型。

2.3 代理模型計(jì)算結(jié)果分析

以熱流Q為例，說明代理模型的計(jì)算精度。經(jīng)過POD方法分解得到100個(gè)基模態(tài)的熱流分布，圖3為每一個(gè)基模態(tài)對應(yīng)的特征值的對數(shù)曲線，特征值的大小表征了基模態(tài)對于熱流分布的貢獻(xiàn)，可以看出前10階包含了大部分熱流的特性。

圖3 基模態(tài)對應(yīng)的特征值Fig.3 Eigenvalues of the Corresponding base Mode

2.3.1 典型工況算例分析

在給定的樣本空間里隨機(jī)取某一個(gè)測試點(diǎn)，進(jìn)行CFD計(jì)算，然后與POD-Kriging代理模型的預(yù)測結(jié)果進(jìn)行比較，在每個(gè)網(wǎng)格點(diǎn)進(jìn)行相對誤差計(jì)算，誤差計(jì)算公式為

式中 ROM為用代理模型方法得到的熱流預(yù)測值；CFD為用CFD-Fastran計(jì)算得到的熱流值。

圖 4為在流條件為飛行高度 H=20 km，攻角α=0°，馬赫數(shù)Ma=7.5的情況下CFD計(jì)算的熱流值與 POD-Kriging模型預(yù)測的熱流值的相對誤差分布。由圖4可以看出，在機(jī)翼前緣處，二者的結(jié)果比較接近，在后緣部分，誤差有所增加，誤差最高的地方集中在前緣后部靠近翼根處，達(dá)到7%左右。結(jié)果表明，對于單個(gè)測試點(diǎn)，POD-Kriging代理模型的預(yù)測結(jié)果比較準(zhǔn)確。

圖4 馬赫數(shù)為7.5來流翼面平均誤差分布Fig.4 Average Error Distribution of the Wing Surface with the Maher Number of 7.5

2.3.2 測試樣本點(diǎn)的相對平均誤差

在輸入?yún)?shù)的空間內(nèi)隨機(jī)選取10個(gè)測試樣本點(diǎn)，測試樣本點(diǎn)的輸入?yún)?shù)如表 3所示。分別采用CFD-Fastran計(jì)算和代理模型計(jì)算熱流分布，利用式（11）計(jì)算每個(gè)測試樣本點(diǎn)的誤差值Pi，10個(gè)樣本點(diǎn)的相對平均誤差為

表3 隨機(jī)選取的10個(gè)樣本點(diǎn)輸入?yún)?shù)Tab.3 Input Parameters of Randomly Selected 10 Sample Points

用POD-Kriging方法得到的10個(gè)測試樣本點(diǎn)的翼面熱流平均相對誤差云圖，如圖5所示。

圖5 10個(gè)樣本點(diǎn)翼面熱流密度平均誤差分布Fig.5 Average Error Distribution of Heat Flow Density of Wing Surface at 10 Sample Points

從圖5中可以看出，翼面誤差基本小于6%，前緣處誤差較小，誤差較大的地方集中在上翼面凸起處和下翼面局部。

2.3.3 POD基模態(tài)個(gè)數(shù)影響

預(yù)測的準(zhǔn)確程度與POD基模態(tài)的個(gè)數(shù)有關(guān)，圖6顯示的是基模態(tài)為15，20，25，30時(shí)，代理模型的熱流預(yù)測值和CFD計(jì)算值的相對誤差。

圖6 不同基模態(tài)個(gè)數(shù)下平均誤差分布的變化云圖Fig.6 Variation of Mean Error Distribution under Different Fundamental Mode Numbers

由圖 6可知，隨著基模態(tài)個(gè)數(shù)的提高，平均誤差逐漸變小，而且當(dāng)基模態(tài)個(gè)數(shù)在20以后，平均誤差的變化較小，誤差最明顯的地方集中在翼面凸起處，且最高誤差在6%左右。

3 熱邊界線性化近似

3.1 熱邊界線性化近似方法

在一定樣本基礎(chǔ)上，通過建立 POD-Kriging代理模型，可以快速得到高超聲速飛行器不同飛行條件下給定壁面溫度的表面熱流分布。而在一般的高超聲速飛行器流固耦合研究中，通常會(huì)采用松耦合的方式，即將流場、固體結(jié)構(gòu)場和溫度場分開計(jì)算。以假定表面溫度得到熱流，再進(jìn)行結(jié)構(gòu)計(jì)算，不斷迭代到誤差在要求范圍內(nèi)。每一次CFD計(jì)算都要耗費(fèi)大量的計(jì)算時(shí)間，給耦合仿真的工程應(yīng)用帶來困難。引入線性化近似的方法[18]可以改善CFD迭代耗時(shí)問題，即：傳遞熱邊界時(shí)，采用線性近似來計(jì)算新壁溫條件下高超聲速飛行器表面氣動(dòng)熱，以代替CFD迭代計(jì)算熱流，必要時(shí)再采用CFD或者代理模型迭代計(jì)算，某些情況下，這種近似在不進(jìn)行CFD迭代的情況下也可以獲得滿足工程需要的計(jì)算結(jié)果。

熱流計(jì)算表達(dá)式為

式中 h為換熱系數(shù)；wT為壁面溫度；awT為絕熱壁面溫度。

如果壁面溫度對換熱系數(shù)的影響較小，就可以以飛行器的換熱系數(shù)和絕熱壁溫作為對象，建立POD-Kriging模型，那么熱流就是壁面溫度的線性函數(shù)。在耦合計(jì)算過程中，當(dāng)表面溫度變化時(shí)，可以通過換熱系數(shù)和絕熱表面溫度直接得到變化后的表面熱流，而不需要通過新的CFD計(jì)算得到表面熱流。

建立絕熱壁溫的 POD-Kriging模型過程與建立熱流的 POD-Kriging模型過程相同，只是樣本點(diǎn)計(jì)算的響應(yīng)輸出是絕熱壁溫。建立換熱系數(shù)的 POD-Kriging模型，需要同時(shí)計(jì)算定壁溫表面熱流和絕熱溫度來獲得。

假設(shè)壁面溫度與來流溫度T∞相同時(shí)，表面熱流為

忽略壁面溫度對換熱系數(shù)的影響，則換熱系數(shù)可以通過下式計(jì)算：

對于每個(gè)樣本點(diǎn)，需要計(jì)算壁溫為來流溫度的熱流分布和絕熱壁溫分布，以此數(shù)據(jù)為基礎(chǔ)，建立換熱系數(shù)和絕熱壁溫的POD-Kriging模型。

3.2 換熱系數(shù)誤差分析

取來流條件為飛行高度 H=29.9 km，攻角α=0.9°，馬赫數(shù) Ma=9.08，采用CFD分別計(jì)算此條件下絕熱壁面的溫度分布awT和給定溫度為300 K條件下的熱流分布由式（15）求出換熱系數(shù)1h。改變壁溫為500 K，求出此時(shí)的熱流同樣求出換熱系數(shù)對比熱流和換熱系數(shù)的變化。

圖7為改變壁溫前后，各網(wǎng)格點(diǎn)的熱流變化。

圖7 改變壁溫引起的熱流的變化Fig.7 Change of Heat Flux Caused by Wall Temperature

圖8為換熱系數(shù)在溫度改變前后的變化。

圖8 改變壁溫引起換熱系數(shù)的變化Fig.8 Change of Heat Transfer Coefficient Caused by Changing Wall Temperature

由圖7、圖8可以看出，與對熱流的影響相比，溫度變化對于換熱系數(shù)的影響較小，變化率在6.5%以內(nèi)。

溫度改變后，由式（13）利用線性化方法求解新的熱流，與CFD得到的熱流進(jìn)行比較，圖9為二者計(jì)算熱流結(jié)果的相對變化。

圖9 線性化處理和CFD結(jié)果的相對變化Fig.9 Relative Change of Linearization and CFD Results

由圖9可以看出，線性化處理得到的壁溫變化后熱流與CFD計(jì)算結(jié)果相對誤差基本在7%以下。在計(jì)算過程中，CFD計(jì)算一次需要10 h以上，而采用線性化處理只需不到1 s，極大提高了計(jì)算效率。

3.3 代理模型結(jié)合線性化近似結(jié)果分析

取飛行高度 H=29.9 km，攻角α=0.9°，馬赫數(shù)Ma=9.08作為來流條件，利用POD-Kriging代理模型快速計(jì)算該狀態(tài)下冷壁熱流和絕熱壁的溫度，再利用式（13）可以求出壁溫改變?yōu)?00 K后的表面熱流，同時(shí)采用CFD-Fastran計(jì)算壁溫改變?yōu)?00 K后該狀態(tài)的表面熱流，圖10顯示的是二者的相對誤差。從圖10中可以看出，前緣誤差較小，在2.5%以內(nèi)，后緣誤差較大，整個(gè)翼面的誤差值在7%以內(nèi)。和圖9相比，誤差較大的網(wǎng)格點(diǎn)有所增多，這是因?yàn)榻^熱壁溫是經(jīng)過POD代理模型的預(yù)測值，增加了誤差。

圖10 POD-Kriging模型和CFD結(jié)果的相對誤差Fig.10 Relative Error of the POD-Kriging Model and the CFD Result

4 結(jié) 論

本文以減少流固耦合過程中熱邊界的求解時(shí)間為目標(biāo)，一方面基于 POD-Kriging代理模型建立快速計(jì)算不同狀態(tài)的氣動(dòng)熱參數(shù)分布的方法；另一方面提出基于換熱系數(shù)的線性近似方法確定熱邊界，以三維機(jī)翼為例，分析了方法的可行性，得到以下結(jié)論：

a）對高超聲速飛行器機(jī)翼氣動(dòng)熱提出的POD-Kriging代理模型結(jié)合的方法，成功實(shí)現(xiàn)模型的降階，代理模型的熱流計(jì)算誤差在7%之內(nèi)；

b）在POD-Kriging建模過程中選取的基模態(tài)個(gè)數(shù)對代理模型預(yù)測結(jié)果會(huì)產(chǎn)生影響，結(jié)果表明增加基模態(tài)個(gè)數(shù)在一定程度上能夠提高預(yù)測精度，當(dāng)其達(dá)到20后，預(yù)測誤差變化很微??；

c）引入換熱系數(shù)線性近似的方法，能夠節(jié)省耦合計(jì)算過程中因壁溫變化引起的流場迭代步驟，在不降低熱流計(jì)算精度的前提下，可以縮短計(jì)算時(shí)間。

本文建立的方法在滿足一定計(jì)算精度的條件下，極大地提高計(jì)算效率，促進(jìn)高超聲速飛行器再入過程中的流固耦合分析或多物理場優(yōu)化設(shè)計(jì)的工程應(yīng)用。