代光月，桂業(yè)偉，國義軍，張 勇，童福林

（中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

氣動加熱問題是高超聲速飛行器研制與發(fā)展的最關(guān)鍵問題之一，準確的熱環(huán)境數(shù)據(jù)是高超聲速飛行器防熱設(shè)計的前提。在高超聲速飛行器的概念研究和初步設(shè)計階段，工程設(shè)計需要以低成本和快速度，較為準確地預(yù)測飛行器表面的熱流分布情況。因此，尋求一種快速且精度較高的高超聲速氣動熱環(huán)境計算方法就顯得非常有意義了。

一種較為簡單且有效的方法是由Cooke提出，DeJarnette推廣的軸對稱比擬法［1］。為了采用軸對稱比擬法求解熱流，必須先確定表面流線形狀和尺度因子，這是最困難也是極其重要的部分，直接決定著計算精度［2］。目前，基于軸對稱比擬的經(jīng)典流線法往往先是利用修正牛頓理論或?qū)嶒灥确椒ǐ@得近似的表面壓力，結(jié)合熵值算得其他邊界層外緣參數(shù)，然后在此基礎(chǔ)上利用簡單流線法、精確流線法等方法跟蹤無粘表面流線，計算尺度因子，最后利用理論和半經(jīng)驗公式計算出表面熱流。該方法計算過程相對簡單，計算速度較快，但由于需要用到壓力甚至壓力的二階導數(shù)項，計算精度相對較差；它的另一個缺點是流場處理過于簡單，在處理復雜外形飛行器熱環(huán)境計算時往往就顯得有些力不從心。國內(nèi)，北京航空航天大學的康宏林［3］，西北工業(yè)大學的呂麗麗［4］和南京航空航天大學的方磊［5］等曾在此方面做過一些工作，取得了一些研究成果。

考慮到純數(shù)值模擬雖然在直接計算熱流方面還存在一些問題，但計算的壓力分布是比較準確的，并且能夠比較精細地刻畫流場結(jié)構(gòu)。本文將利用這一特點，基于軸對稱比擬概念，首先利用有限體積法數(shù)值求解Euler方程獲得較為準確的邊界層外緣流場參數(shù)，然后基于有限元的四節(jié)點單元變換方法，直接利用笛卡爾坐標系下的三維速度分量計算無粘表面流線和尺度因子；在獲取無粘表面流線和尺度因子的基礎(chǔ)上，利用Zoby、Moss和Sutton提出的熱流公式計算表面熱流，從而實現(xiàn)數(shù)值算法和工程算法的耦合。由于該方法直接采用流場速度計算無粘表面流線，克服了經(jīng)典流線法需要壓力及其二階導數(shù)且不能保證導數(shù)連續(xù)的不足，而且直接利用笛卡爾坐標系下的速度分量計算尺度因子，提高了對復雜外形的適應(yīng)能力。將上述方法用于求解球錐在不同攻角條件下的表面熱流，并將計算結(jié)果同經(jīng)典流線法結(jié)果及實驗值進行比較，從而對方法進行考核驗證。

1 無粘流場的數(shù)值算法

無粘流場的數(shù)值計算部分，采用已經(jīng)過多次驗證、比較成熟的CFD軟件平臺（高超聲速軟件平臺CHANT［6］）。直角坐標系下，無量綱化后的三維Euler方程可表述如下：

對于空間離散，本文采用一種既能保持解的單調(diào)性，又具有較高精度的NND格式［7］；而對于時間離散，采用隱式LU－SGS［8］時間推進方法。

2 流線和尺度因子的計算

2.1 表面流線的確定

在笛卡爾坐標系下，流線的定義式為：

對式（2），從某一初始位置（x0，y0，z0）積分，則可求得表面流線分布。

2.2 尺度因子的確定

空間一點（x，y，z）在笛卡爾坐標系下可寫為：

若用（s，β，n）代表正交坐標系，其中s為物面上沿流線方向，β為物面上垂直于流線的方向，n為垂直物面方向，那么物面上的點可寫為：

那么，垂直方向的單位矢量eβ可寫為：

因此有

若物面外形可用F（x，y，z）＝x－x（y，z）＝0表示，則有：

（s，β，n）為正交坐標系，即eβ＝en×es，則eβ可寫為：

由式（5）和式（9）可知：

結(jié)合hβ的定義式，可知：

nx，ny，nz的值可利用網(wǎng)格坐標求出位置向量，由向量的矢量積得到?？梢钥吹?，若已知該點處的?x／?β，?y／?β，?z／?β，則可求出該點的尺度因子hβ，又從文獻［9］可知，在物體表面有：

這樣，當給定hβ一個初值，由式（10）即可求得?x／?β，?y／?β，?z／?β的初值，再積分式（12），即可獲得沿流線的?x／?β，?y／?β，?z／?β，進而求得尺度因子的值。

從上面的分析可以看出，在對式（12）進行積分時，需已知（u，v，w）關(guān)于（x，y，z）的偏導數(shù)項，求解方法介紹如下［9］。

圖1為一個計算單元＾Ω到物理單元Ωe的映射。通過有限元的形狀函數(shù)，可以把從（ξ，η）空間到物理空間（u，v，w）的映射表示為：

圖1 從計算單元＾Ω到物理單元Ωe的映射Fig.1 Mapping from master element to a physical element

對照圖1可知，在計算單元每個子區(qū)域上有：

區(qū)域I

區(qū)域II

區(qū)域III

區(qū)域IV

這樣每個區(qū)域中關(guān)于u＝u（ξ，η）的坐標可由式（13）求得，進而可求得每個區(qū)域?u／?ξ，?u／?η，的表達式。以區(qū)域中心點處u的偏導數(shù)?u／?ξ，?u／?η為例，其值應(yīng)為四個區(qū)域的平均：

同理可以計算出v，w關(guān)于ξ，η偏導數(shù)的表達式。

同時，根據(jù)空間的映射變換關(guān)系，在進行一系列矩陣變換的基礎(chǔ)上，可得出ξ，η關(guān)于x，y，z偏導數(shù)的表達式，以?ξ／?x為例：

其中｜J｜是轉(zhuǎn)換的Jacobian行列式，數(shù)值上為物理單元與計算單元的面積比，表達式為：

又因（u，v，w）關(guān)于（x，y，z）的偏導數(shù)具有如下表達式，以u為例：

這樣，（u，v，w）關(guān)于（x，y，z）偏導數(shù)項的值就可以求得了。

2.3 逆向跟蹤流線

上述方法可用于逆向跟蹤流線，只需在式（2）的右邊加一負號即可，這樣，對于任意一個想要考察其熱量值的特定點，都可以通過指定該點為起始點的方式直接求出該點的熱流值，而不需要插值附近點獲得；同時，還大大方便了表面流線的布置，特別是在有攻角情況下，由于流線向背風區(qū)匯集，傳統(tǒng)方法很難使流線均勻覆蓋飛行器表面，而利用逆向跟蹤流線，可使流線分布的十分均勻。

3 高超聲速熱流計算

3.1 駐點熱流計算

本文采用經(jīng)典的Fay－Riddell公式［1］計算駐點熱流，并將其與實驗測得的熱流值進行比較。Fay－Riddell公式的表達式為：

3.2 非駐點熱流計算

對于非駐點區(qū)域繞流，本文選取了由Zoby，Moss和Sutton等提出的應(yīng)用廣泛且具有較好精度的熱流計算公式［10］。層流熱流的計算公式為：

湍流熱流的計算公式為：

其中，Reθ，e是基于邊界層動量厚度的雷諾數(shù)，邊界層動量厚度θ的表達式分別為：

式（25）中，h為尺度因子，s為沿流線的長度。

4 方法驗證

4.1 計算條件

本文選擇了有詳盡實驗數(shù)據(jù)的三維球錐繞流問題進行計算，以達到對方法進行考核驗證的目的。計算所用模型的具體尺寸為：頭部曲率半徑為3.835mm，半錐角為12.84°，錐尾距頭部距離為69.55mm。在單純有攻角情況下，由于流場左右對稱，取半流場進行計算，來流參數(shù)值與實驗數(shù)據(jù)對應(yīng)工況［11］相一致：p∞＝61Pa，T∞＝49.8K，Ma∞＝9.86。網(wǎng)格類型為單塊結(jié)構(gòu)網(wǎng)格。通過網(wǎng)格實驗，兼顧計算精度與計算效率，網(wǎng)格數(shù)為39×41×31。計算時，取攻角分別為0°、8°和16°。

4.2 無粘流場結(jié)果分析

圖2是球錐在攻角分別為0°、8°和16°時球錐表面壓力分布云圖。由圖可以看出，由于來流馬赫數(shù)較高，在球錐的頭部，產(chǎn)生了一道很強的頭激波，波后存在一個較強的高壓區(qū)，隨著計算攻角的增大，迎風面的壓力值逐漸增大，流場結(jié)構(gòu)清晰，符合物理規(guī)律。

4.3 流線分析

圖3是球錐在攻角分別為0°、8°和16°時的壁面流線分布圖?？梢钥闯?，流線形狀隨著攻角增大而扭曲，攻角越大，扭曲程度越厲害，向背風區(qū)匯集，壁面流線分布情況較好，符合實際物理情況。

圖2 攻角為0°、8°和16°時球錐表面壓力分布云圖Fig.2 Pressure distribution of conic，forα＝0°，8°and 16°

圖3 球錐在攻角分別為0°、8°和16°時的壁面流線分布圖Fig.3 Surface streamlines of conic，forα＝0°，8°and 16°

4.4 駐點區(qū)域熱流計算結(jié)果分析

高超聲速飛行器的駐點加熱非常重要，它是最具代表性的點。同時，為了將算得的熱流密度值進行歸一化處理，便于比較，也需首先求出駐點熱流值。本文采取經(jīng)典的Fay－Riddell公式計算駐點熱流，并將其與實驗測得的熱流值進行比較，結(jié)果見下表。

表1 不同攻角下球錐駐點熱流值比較Table 1 Comparison of heating rates at stagnation point

從上表可以看出，計算出的熱流值與實驗值吻合的較好，因此，本文計算得出的熱流密度值是符合物理規(guī)律，可用的。

4.5 非駐點區(qū)域熱流計算結(jié)果分析

圖4是球錐在攻角分別為0°、8°和16°時的迎風子午線上的熱流分布示意圖，圖5是球錐在8°攻角時背風子午線上的熱流分布示意圖。L是參考長度，即：錐尾距頭部距離；圖中熱流密度值為進行歸一化后的熱流值。從圖中可以看出：無論是有攻角還是無攻角條件，迎風和背風子午線上熱流計算結(jié)果與實驗值均吻合的較好，計算精度較經(jīng)典工程算法有較大提高。

圖4 球錐在攻角分別為0°、8°和16°時的迎風子午線上的熱流分布示意圖Fig.4 Windward centerline heating rate distribution，forα＝0°，8°and 16°

圖5 球錐8°時的背風子午線上的熱流分布示意圖Fig.5 leeward centerline heating rate distribution

5 結(jié)束語

基于軸對稱比擬，將數(shù)值算法和工程算法耦合起來，利用有限體積法數(shù)值求解Euler方程獲得較為準確的邊界層外緣無粘流場參數(shù)，然后基于有限元的四節(jié)點單元變換方法，直接利用笛卡爾坐標系下的三維速度分量計算無粘表面流線和尺度因子；在獲取無粘表面流線和尺度因子的基礎(chǔ)上，利用Zoby、Moss和Sutton提出的熱流公式計算表面熱流。該方法直接采用流場速度計算無粘表面流線，克服了經(jīng)典流線法需要壓力及其二階導數(shù)且不能保證導數(shù)連續(xù)的不足，而且直接利用笛卡爾坐標系下的速度分量計算尺度因子，提高了對復雜外形的適應(yīng)能力。為了對方法進行考核驗證，將上述方法用于求解球錐在不同攻角條件下的表面熱流分布情況，將計算結(jié)果同經(jīng)典流線法結(jié)果和實驗值進行對比，結(jié)果表明：本文計算結(jié)果與實驗值吻合的較好，計算精度較經(jīng)典工程算法有較大提高。將程序進一步完善，使其適用于多塊結(jié)構(gòu)網(wǎng)格和對高精度網(wǎng)格技術(shù)的探討將是下一步研究的重點。

［1］中國人民解放軍總裝備部軍事訓練教材編輯工作委員會.高超聲速氣動熱和熱防護［M］.北京：國防工業(yè)出版社，2003.

［2］HAMILTON H H，WEILMUENSTER K J.Application of axisymmetric analogue for calculating heating in threedimensional flows［A］.AIAA 23rd Aerospace Sciences Meeting［C］.January，1985.

［3］康宏琳，閻超，李亭鶴，等.高超聲速再入鈍頭體表面熱流計算［J］.北京航空航天大學學報，2006，32（12）：1395－1398.

［4］呂麗麗.高超聲速氣動熱工程算法研究［D］.［碩士學位論文］.西北工業(yè)大學，2005.

［5］方磊.基于流線跟蹤法的氣動熱工程計算研究［D］.［碩士學位論文］.南京航空航天大學，2008.

［6］中國空氣動力研究與發(fā)展中心計算空氣動力學研究所高超聲速研究室.高超CFD平臺0.0版［M］.2003.

［7］張涵信.無波動、無自由參數(shù)的耗散差分格式［J］.空氣動力學學報，1988，6：143－164.

［8］SHAROV D，NAKAHASHI K.Reordering of 3Dhybrid unstrucrured grids for vectorized LU－SGS navierstokes computations［R］.AIAA－97－2102.

［9］WANG K C.An axisymmetric analog two－layer convective heating procedure with application to the evaluation of space shuttle orbiter wing leading edge and windward surface heating［R］.NASA C R 188343，1994.

［10］ZOBY E V，MOSS J N，SUTTON K.Approximate convective heating equations for hypersonic flows［J］.Journal of Spacecraft and Rockets，1981，18（1）.

［11］MILLER C G.Experimental and predicted heating distributions for biconics at incidence in air at Mach 10［R］.NASA－TP－2334，1984.