王 洋，袁軍婭，王洪興

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

0 引 言

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

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

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

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

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

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

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

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

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

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

其中，

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

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

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

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

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

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

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

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

2 模型以及算例分析

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

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

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

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

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

機翼的氣動熱設(shè)計變量和設(shè)計空間如表1所示。

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

2.2 POD-Kriging代理模型的建立

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

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

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

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

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

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

2.3.1 典型工況算例分析

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

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

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

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

2.3.2 測試樣本點的相對平均誤差

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

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

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

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

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

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

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

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

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

3 熱邊界線性化近似

3.1 熱邊界線性化近似方法

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

熱流計算表達式為

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

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

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

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

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

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

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

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

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

圖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得到的熱流進行比較，圖9為二者計算熱流結(jié)果的相對變化。

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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