董維中，丁明松，高鐵鎖，江 濤

（中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000）

0 引 言

高超聲速飛行器在大氣層飛行過(guò)程中，不僅存在高熱流的氣動(dòng)熱環(huán)境問(wèn)題，還存在嚴(yán)重的高溫氣體效應(yīng)／非平衡效應(yīng)問(wèn)題，這就給高超聲速飛行器熱防護(hù)設(shè)計(jì)和氣動(dòng)設(shè)計(jì)增加非常大的困難。在這種情況下，準(zhǔn)確預(yù)測(cè)高超聲速飛行器氣動(dòng)力特性和表面氣動(dòng)熱環(huán)境就非常關(guān)鍵。

對(duì)高溫氣體效應(yīng)／非平衡效應(yīng)問(wèn)題研究越來(lái)越得到重視，國(guó)內(nèi)外建立了有關(guān)的計(jì)算軟件和發(fā)表了相當(dāng)多論 文［1－8］。2009 年 前 后，意 大 利 航 天 研 究 中 心（CIRA）的 Votta R等人［7－8］，考慮高溫空氣的化學(xué)非平衡和熱力學(xué)非平衡效應(yīng)、滑移效應(yīng)和表面溫度分布，表面溫度取不同的等溫度值或表面熱輻射平衡條件確定的溫度分布，通過(guò)數(shù)值求解非平衡Navier－Stokes方程的CFD方法（H3NS），對(duì)比分析了熱化學(xué)非平衡效應(yīng)、表面催化效應(yīng)和表面溫度分布對(duì)飛行器氣動(dòng)力／熱特性的影響。

在數(shù)值計(jì)算氣動(dòng)熱中，影響準(zhǔn)確預(yù)測(cè)高超聲速飛行器氣動(dòng)熱環(huán)境的因素非常多，不僅有計(jì)算網(wǎng)格和計(jì)算格式及其邊界處理方法，還有氣體模型、表面溫度和表面催化特性等等。在氣體模型中，不僅有完全氣體模型和化學(xué)反應(yīng)混合氣體模型，還有氣體熱力學(xué)模型。在化學(xué)反應(yīng)混合氣體模型中，常用的有5組分、7組分和11組分空氣化學(xué)反應(yīng)模型。在這些影響準(zhǔn)確預(yù)測(cè)高超聲速飛行器氣動(dòng)熱環(huán)境因素中，目前國(guó)內(nèi)外對(duì)空氣化學(xué)反應(yīng)的組分模型、熱力學(xué)非平衡模型和表面溫度對(duì)氣動(dòng)熱的計(jì)算結(jié)果的影響規(guī)律還不是非常清楚，在高馬赫數(shù)和熱化學(xué)非平衡條件下，氣動(dòng)熱數(shù)值隨表面溫度變化有什么特征，這些問(wèn)題都需要深入分析研究。

本文考慮高溫空氣化學(xué)反應(yīng)效應(yīng)和熱力學(xué)非平衡效應(yīng)，利用自主開(kāi)發(fā)的氣動(dòng)物理計(jì)算軟件系統(tǒng)（AEROPH）中的高超聲速飛行器高溫氣體／非平衡效應(yīng)流場(chǎng)計(jì)算軟件（AEROPH＿Flow），開(kāi)展熱化學(xué)非平衡氣動(dòng)熱數(shù)值計(jì)算，分析高溫空氣化學(xué)反應(yīng)的組分模型、熱力學(xué)非平衡模型和表面催化特性對(duì)氣動(dòng)熱的計(jì)算結(jié)果的影響，考察在高馬赫數(shù)和熱化學(xué)非平衡條件下氣動(dòng)熱數(shù)值隨表面溫度的變化特征，提升準(zhǔn)確預(yù)測(cè)高超聲速飛行器氣動(dòng)熱環(huán)境的能力。

1 非平衡流場(chǎng)控制方程和計(jì)算方法

控制方程是三維熱化學(xué)非平衡Navier－Stokes方程，無(wú)量綱化形式如下［9－10］：

其中，

這里i＝（1，2，…，N），N是所考慮的組分方程數(shù)。ρi是組分i的密度，EV是分子組分的總振動(dòng)能量，Re是Reynolds數(shù)，wi和wV是化學(xué)非平衡源項(xiàng)和振動(dòng)非平衡能量源項(xiàng)。

對(duì)于熱化學(xué)非平衡流動(dòng)情況，狀態(tài)方程和能量關(guān)系式如下［9－10］：

這里ei，tr和ei，e是組分i的平動(dòng)和轉(zhuǎn)動(dòng)能、束縛電子的激發(fā)能。θVj是組分j的振動(dòng)特征溫度，g0i和g1i是組分i的束縛電子第0和第1能級(jí)的簡(jiǎn)并因子，θei是第一能級(jí)的特征溫度。Δ是組分生成能?？諝饨M分的物理化學(xué)數(shù)據(jù)見(jiàn)文獻(xiàn)［9－10］。

對(duì)于化學(xué)非平衡流動(dòng)情況，在控制方程中去掉振動(dòng)非平衡能量方程和相關(guān)項(xiàng)，狀態(tài)方程和能量關(guān)系式如下［9－10］：

完全氣體的控制方程和狀態(tài)方程見(jiàn)文獻(xiàn)［9］。方程（1）采用全隱式的對(duì)稱型TVD格式進(jìn)行差分離散，粘性項(xiàng)用中心差分格式離散。

2 空氣化學(xué)反應(yīng)模型

較常用的空氣化學(xué)反應(yīng)模型有5組分（O2，N2，NO，O，N）、7組分（O2，N2，NO，O，N，NO＋，e）和11組分（O2，N2，NO，O，N，NO＋，e，O＋2，N＋2，O＋，N＋）的模型。在本文研究中，選用了5組分、7組分和11組分的Dunn－Kang空氣化學(xué)反應(yīng)模型。

3 熱力學(xué)非平衡模型

對(duì)于熱力學(xué)非平衡，只考慮振動(dòng)非平衡效應(yīng)，考慮振動(dòng)能量非平衡的組分有O2，N2，NO，NO＋，O＋2和等分子組分。為了減少控制方程數(shù)量，振動(dòng)非平衡效應(yīng)考慮一個(gè)振動(dòng)溫度。振動(dòng)非平衡能量源項(xiàng)有如下形式：

這里、EVj和τVj分別是分子組分j的平衡振動(dòng)能、非平衡振動(dòng)能和振動(dòng)能松弛的特征時(shí)間，τVj具體計(jì)算公式計(jì)算見(jiàn)文獻(xiàn)［9－10］。

4 氣動(dòng)熱計(jì)算方法

對(duì)于完全氣體模型，表面熱流計(jì)算公式為：

對(duì)于化學(xué)非平衡氣體模型，表面熱流由對(duì)流熱流和組分?jǐn)U散熱流兩部分組成，計(jì)算公式為：

對(duì)于熱化學(xué)非平衡氣體模型，表面熱流由平動(dòng)項(xiàng)的對(duì)流熱流、振動(dòng)項(xiàng)對(duì)流熱流和組分?jǐn)U散熱流三部分組成，計(jì)算公式為：

這里k和kV是平動(dòng)項(xiàng)和振動(dòng)項(xiàng)熱傳導(dǎo)系數(shù)。hi和Di是組分i的焓和擴(kuò)散系數(shù)。對(duì)于完全非催化條件，有＝0；對(duì)于完全催化條件，有ciw＝ci∞（來(lái)流條件）。

5 空氣化學(xué)反應(yīng)模型和熱力學(xué)模型對(duì)氣動(dòng)熱計(jì)算影響的分析

針對(duì)高超聲速飛行器，我們建立了氣動(dòng)物理計(jì)算軟件系統(tǒng)，并開(kāi)展了大量的驗(yàn)證與確認(rèn)工作，在高溫氣體／非平衡效應(yīng)流場(chǎng)計(jì)算軟件方面，開(kāi)展了大量的與國(guó)外文獻(xiàn)［11－13］以及國(guó)內(nèi)外的飛行試驗(yàn)、地面試驗(yàn)和理論計(jì)算的結(jié)果對(duì)比分析工作［9－10，14－17］。

計(jì)算模型：球頭模型，RN＝35mm；計(jì)算狀態(tài)：H＝55km，M＝20；熱力學(xué)非平衡模型：一溫度模型（T）和兩溫度模型（T－Tv）；空氣化學(xué)反應(yīng)模型：5組分、7組分和11組分化學(xué)反應(yīng)模型；表面條件：Tw＝300K，完全催化（FCW）和完全非催化（NCW）。

圖1給出了不同組分空氣化學(xué)反應(yīng)模型的熱流分布，從圖中可以看到，7組分模型和11組分模型計(jì)算的熱流分布幾乎沒(méi)有差別，5組分模型計(jì)算的熱流值偏小，在M＝20和完全催化條件下熱流偏小13%左右，原因可能是5組分模型沒(méi)有考慮電離效應(yīng)。圖2給出了不同熱力學(xué)模型的熱流分布，從圖中可以看到，不同熱力學(xué)模型計(jì)算的熱流分布有一定差異，最大差異在10%左右，兩溫度模型計(jì)算的熱流值高于一溫度模型計(jì)算的熱流值。從總的來(lái)看，在同樣催化條件下，不同組分空氣化學(xué)反應(yīng)模型和不同熱力學(xué)模型計(jì)算的熱流分布差異主要集中在球頭模型的駐點(diǎn)區(qū)域，后面差別較小。

6 表面溫度對(duì)氣動(dòng)熱計(jì)算影響的分析

計(jì)算模型：球錐模型，RN＝35mm，θ＝5°，L＝2m；計(jì)算狀態(tài)：H＝60km和55km，M＝20；熱力學(xué)非平衡模型是一溫度模型，空氣化學(xué)反應(yīng)模型是7組分模型；表面條件：Tw＝300K，800K，1500K，2000K，完全催化和完全非催化。

圖1 不同組分空氣化學(xué)反應(yīng)模型的熱流分布（NCW和FCW）Fig.1 Distribution of heat transfer rate on the hemisphere for different species chemical reaction models（NCW and FCW）

圖2 不同熱力學(xué)模型的熱流分布（NCW和FCW）Fig.2 Distribution of heat transfer rate on the hemisphere for different thermodynamic nonequilibrium models（NCW and FCW）

圖3 球錐模型H＝55km的壓力和表面熱流分布云圖Fig.3 Contours of pressure and heat transfer rate for the sphere－cone body at the altitude of 55km

圖3給出Rn＝3 5mm了球錐模型H＝5 5km、M∞＝18、Tw＝300K、Ns＝7、NCW 狀態(tài)熱流分布云圖，從圖可以看出：高熱流區(qū)域和空氣化學(xué)反應(yīng)最強(qiáng)區(qū)域在球錐體的頭部區(qū)域，在后身部區(qū)域熱流低和空氣化學(xué)反應(yīng)弱。

圖4給出了球錐模型駐點(diǎn)熱流和球錐末端熱流及其成分隨表面溫度變化的分布。從圖可以看出：（1）不同表面溫度計(jì)算的熱流分布有一定的差異，在Tw＝300K～800K之間，表面溫度越高，表面熱流越高，在Tw＝800K～2000K之間，隨表面溫度增加，表面熱流先緩慢升高，Tw＝1500K之后，有降低趨勢(shì)；（2）表面溫度增加，組分?jǐn)U散項(xiàng)熱流一直增加，表面溫度從300K到2000K，組分?jǐn)U散項(xiàng)熱流在總熱流中的百分比從33%左右到40%左右；（3）表面溫度變化對(duì)球錐頭部區(qū)域熱流值影響大，但對(duì)球錐身部區(qū)域影響??；（4）表面催化效應(yīng)對(duì)球錐體的頭部區(qū)域熱流影響大，對(duì)后身部區(qū)域熱流影響小。

圖4 球錐模型駐點(diǎn)和末端熱流隨表面溫度變化的分布Fig.4 Variation of the heat transfer rate at the stagnational point and the end of the sphere－cone body for different surface temperatures

從分析可以看出：氣動(dòng)熱隨著表面溫度的變化規(guī)律在高馬赫數(shù)和非平衡條件下就變得非常復(fù)雜，不能再認(rèn)為氣動(dòng)熱遵從隨著表面溫度的升高而降低的規(guī)律。

7 高升阻比外形不同表面溫度氣動(dòng)熱計(jì)算分析

圖5 高升阻比外形駐點(diǎn)熱流隨飛行高度變化分布Fig.5 Variation of the heat transfer rate at the stagnational point of the high lift－to－drag hypersonic vehicle for different flight altitudes and surface temperatures

對(duì)于典型簡(jiǎn)化高升阻比外形，頭部是RN＝35mm的半球，計(jì)算狀態(tài)：H＝45km～70km，M＝11～23，圖5給出了表面溫度Tw＝600K和Tw＝1000K時(shí)高升阻比外形駐點(diǎn)熱流隨飛行高度變化分布。由圖可以看出：（1）從H＝70km到H＝45km，駐點(diǎn)熱流最大值在H＝55km；（2）表面材料催化效應(yīng)對(duì)熱流值影響較大，不同表面溫度熱流值有所不同，當(dāng)表面溫度Tw＝600K時(shí)，對(duì)于完全非催化條件，駐點(diǎn)熱流在1650kW／m2至3000kW／m2之間，而對(duì)于完全催化條件，駐點(diǎn)熱流在1900kW／m2至4600kW／m2之間，當(dāng)表面溫度Tw＝1000K時(shí)，對(duì)于完全非催化條件，駐點(diǎn)熱流在1500kW／m2至3200kW／m2之間，而對(duì)于完全催化條件，駐點(diǎn)熱流在1800kW／m2至5000kW／m2之間；（3）完全氣體（PG）模型計(jì)算的熱流在完全非催化和完全催化熱流之間；（4）表面溫度對(duì)熱流計(jì)算存在一定影響，表面溫度由600K上升到1000K時(shí)，H＝45km時(shí)駐點(diǎn)熱流略微下降，而其他高度熱流都有升高。

圖6給出了H＝55km不同表面溫度和不同飛行馬赫數(shù)高升阻比外形熱流變化分布。由圖可以看到：表面溫度變化對(duì)頭部區(qū)域的熱流影響比較明顯，當(dāng)M＝10時(shí)，熱流值隨著表面溫度的增加而下降；當(dāng)M＝14時(shí)，在Tw＝300K～600K時(shí)，熱流值隨著表面溫度的增加而增加，在Tw＝600K～1500K時(shí)，熱流值隨著表面溫度增加變化緩慢，且有下降趨勢(shì)；當(dāng)M＝18時(shí)，在Tw＝300K～900K時(shí)，熱流值隨著表面溫度的增加而增加，在Tw＝900K～1500K時(shí)，熱流隨著表面溫度增加變化緩慢，且有下降趨勢(shì)。從分析結(jié)果看，高馬赫數(shù)和熱化學(xué)非平衡條件下，在計(jì)算氣動(dòng)熱時(shí)，表面溫度最好取保守上限值，這樣就可以找到最大的駐點(diǎn)熱流值。

圖6 高升阻比外形不同表面溫度和飛行馬赫數(shù)條件下駐點(diǎn)熱流變化分布（FCW）Fig.6 Variation of the heat transfer rate at the stagnational point of the high lift－to－drag hypersonic vehicle for different flight Mach number and surface temperature

8 結(jié)束語(yǔ)

通過(guò)本文的研究工作，可得到如下結(jié)論：

（1）對(duì)高溫氣體效應(yīng)／非平衡效應(yīng)嚴(yán)重的區(qū)域，為了準(zhǔn)確預(yù)測(cè)高超聲速飛行器氣動(dòng)熱環(huán)境，我們要根據(jù)飛行環(huán)境的熱化學(xué)機(jī)制或空氣化學(xué)反應(yīng)和非平衡效應(yīng)的強(qiáng)弱，選擇適當(dāng)?shù)目諝饣瘜W(xué)反應(yīng)模型和熱力學(xué)模型，不能隨意用完全氣體模型、少組分的空氣化學(xué)反應(yīng)模型以及不考慮熱力學(xué)非平衡效應(yīng)，特別是飛行速度非常高的情況。在本文的計(jì)算條件下，數(shù)值分析表明：（a）空氣化學(xué)反應(yīng)模型7組分模型和11組分模型計(jì)算的熱流分布幾乎沒(méi)有差別，而5組分模型計(jì)算的熱流值偏小達(dá)10%以上；（b）不同熱力學(xué)模型計(jì)算的熱流分布差異較小，最大差異在10%左右，兩溫度模型計(jì)算的熱流值高于一溫度模型計(jì)算的熱流值。

（2）高馬赫數(shù)和熱化學(xué)非平衡條件下，氣動(dòng)熱的數(shù)值隨著表面溫度的變化規(guī)律變得非常復(fù)雜，不能再認(rèn)為氣動(dòng)熱的數(shù)值遵從隨著表面溫度的升高而降低的規(guī)律；在計(jì)算氣動(dòng)熱時(shí)，表面溫度最好取接近真實(shí)飛行情況分布和不同的固定值，這樣就可以找到最大的或準(zhǔn)確的熱流值，在高超聲速飛行器的熱防護(hù)設(shè)計(jì)中更有應(yīng)用價(jià)值。

本文只是分析了在再入或飛行速度在7km／s以下的條件，還需要開(kāi)展系統(tǒng)的分析研究，如更高再入或飛行速度、不同模型尺寸和變表面溫度等情況。

［1］GNOFFO P A.Applications of program LAURA to threedimensional AOTV flowfields［R］.AIAA－86－0565.

［2］HARTUNG L C.Development of a nonequilibrium radiative heating prediction method for coupled flowfield solutions［R］.AIAA 91－1406.

［3］KEENAN J A，CANDLER G V.Simulation of graphite ablation and oxidation under re－entry conditions［R］.AIAA－94－2083.

［4］KEENAN J A，CANDLER G V.Simulation of ablation in earth atmospheric entry［R］.AIAA 93－2789，1993.

［5］曾明，杭建，林貞彬，等.不同熱化學(xué)非平衡模型對(duì)高超聲速噴管流場(chǎng)影響的數(shù)值分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2008，24（3）：346－349.

［6］柳軍，劉偉，曾明，等.高超聲速三維熱化學(xué)非平衡流場(chǎng)的數(shù)值模擬［J］.力學(xué)學(xué)報(bào)，2003，35（6）：730－735.

［7］PEZZELLA G，VOTTA R.Finite rate chemistry effects on the high altitude aerodynamics of an apollo－shaped reentry capsule［R］.AIAA 2009－7306.

［8］VOTTA R，et al.Hypersonic low－density aerothermodynamics of orion－like exploration vehicle［J］.Journal ofSpacecraftandRockets，2009，46（4）：781－787.

［9］董維中.熱化學(xué)非平衡效應(yīng)對(duì)高超聲速流動(dòng)影響的數(shù)值計(jì)算與分析［D］.［博士學(xué)位論文］.北京航空航天大學(xué)，1996.

［10］董維中.氣體模型對(duì)高超聲速再入鈍體氣動(dòng)參數(shù)計(jì)算影響的研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19（2）：197－202.

［11］NETTERFIELD M P.Validation of a Navier－Stokes code for thermochemical nonequilibrium flows［R］.AIAA 92－2878.

［12］MUYLAERT J，et al.Standard model testing in the European high facility F4and extrapolation to flight［R］.AIAA－92－3905.

［13］CANDER G V，MACCORMACK R W.The computation of hypersonic ionized flows in chemical and thermal nonequilibium［R］.AIAA－88－0511.

［14］董維中，高鐵鎖，丁明松.高超聲速非平衡流場(chǎng)多個(gè)振動(dòng)溫度模型的數(shù)值研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2007，25（1）：1－6.

［15］董維中，高鐵鎖，張志成.再入體實(shí)驗(yàn)?zāi)Ｐ蜔峄瘜W(xué)非平衡繞流流場(chǎng)的數(shù)值分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2008，26（2）：163－166.

［16］董維中，樂(lè)嘉陵，劉偉雄.駐點(diǎn)壁面催化速率常數(shù)確定的研究［J］.流體力學(xué)實(shí)驗(yàn)與測(cè)量，2000，14（3）：1－6.

［17］董維中，樂(lè)嘉陵，高鐵鎖.鈍體標(biāo)模高焓風(fēng)洞試驗(yàn)和飛行試驗(yàn)相關(guān)性的數(shù)值分析［J］.流體力學(xué)實(shí)驗(yàn)與測(cè)量，2002，16（2）：1－8.