基于軸對稱比擬的高超聲速飛行器表面熱環境數值與工程耦合算法研究

代光月，桂業偉，國義軍，張 勇，童福林

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

0 引 言

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

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

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

1 無粘流場的數值算法

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

對于空間離散，本文采用一種既能保持解的單調性，又具有較高精度的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）可知：

結合hβ的定義式，可知：

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

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

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

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

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

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

區域I

區域II

區域III

區域IV

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

同理可以計算出v，w關于ξ，η偏導數的表達式。

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

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

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

這樣，（u，v，w）關于（x，y，z）偏導數項的值就可以求得了。

2.3 逆向跟蹤流線

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

3 高超聲速熱流計算

3.1 駐點熱流計算

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

3.2 非駐點熱流計算

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

湍流熱流的計算公式為：

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

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

4 方法驗證

4.1 計算條件

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

4.2 無粘流場結果分析

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

4.3 流線分析

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

圖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 駐點區域熱流計算結果分析

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

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

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

4.5 非駐點區域熱流計算結果分析

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

圖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 結束語

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

［1］中國人民解放軍總裝備部軍事訓練教材編輯工作委員會.高超聲速氣動熱和熱防護［M］.北京：國防工業出版社，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］.［碩士學位論文］.西北工業大學，2005.

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

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

［7］張涵信.無波動、無自由參數的耗散差分格式［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.