采用態－態模型的熱化學非平衡噴管流數值研究

徐 丹， 曾 明， 張 威， 柳 軍

（國防科學技術大學航天科學與工程學院，湖南長沙 410073）

采用態－態模型的熱化學非平衡噴管流數值研究

徐 丹， 曾 明， 張 威， 柳 軍

（國防科學技術大學航天科學與工程學院，湖南長沙 410073）

數值求解耦合態－態模型的N2／N準一維噴管非平衡流方程，獲得駐室溫度3000 K～10000 K，壓力1 atm～10 atm條件下的流場參數分布和詳細的振動能級分布，分析噴管流動中的宏觀和微觀特性，并對這類復合占優的流動考察雙溫度模型中假設的合理性.

態－態模型；熱化學非平衡；振動能級分布；高超聲速噴管流動；數值模擬

0 引言

高超聲速高焓流動中常同時存在熱力和化學非平衡現象.對熱力非平衡流動，氣體各種模式的內能無法用統一的溫度表征.一般而言，轉動與平動內能模式較易達到平衡，可以用同一溫度描述，記作平動／轉動溫度或平動溫度T，且平動和轉動能級分布為該溫度下的Boltzmann分布.在電子能可以忽略的情況下，熱力非平衡主要表現為振動非平衡.有兩類描述流動中熱力非平衡現象和熱力－化學耦合的方法.一類是多溫度模型法，以雙溫度模型法［1］為代表，通過振動能變化率方程描述振動非平衡過程，以不同于平動溫度的振動溫度Tv表征非平衡振動能，將化學反應控制溫度設為平動溫度與振動溫度的某種平均值以反映熱力－化學耦合；另一類是態－態模型法［2－3］，直接描述分子在各振動能級間的非平衡躍遷過程，給出各能級粒子數變化率方程，熱力－化學耦合則通過分別描述各能級上分子的化學反應過程自然體現.

雙溫度模型法引入振動能變化率方程和流動控制方程耦合求解，得到流場各點的非平衡振動能后借用平衡振動能和溫度的關系式，折算出振動溫度Tv，它與平動溫度的差別能直觀地反映流動中的振動非平衡現象.振動溫度進入化學反應控制溫度也是一個簡便地體現熱力－化學耦合的方法.雙溫度模型及其后來拓展的多溫度模型一定程度上反映了流動中的熱力非平衡現象和熱力－化學耦合，在高溫非平衡流的數值模擬中獲得了廣泛的應用，在許多條件下也獲得了良好的數值結果，但是也存在局限［4］，例如在復合區的效果不如離解區，無法合理解釋AVCO實驗數據，等.其實化學反應控制溫度（即化學反應速率系數Arrhenius公式中的溫度）中包含振動溫度的內涵是值得進一步分析的.Arrhenius速率系數公式中的溫度參數一方面是反映粒子碰撞頻率，與隨機熱運動（平動）速度有關；另一方面則是反映碰撞能夠引起化學反應的概率，與參與碰撞粒子的內能（包括平動、轉動、振動、電子能）有關，這里以溫度反映粒子內能的物理內涵其實是基于參與反應的粒子能級分布滿足該溫度下的平衡分布（Boltzmann分布）的.因此從這個角度看，在化學反應控制溫度中引入振動溫度，可以說隱含著對振動非平衡流采用了振動能級分布滿足振動溫度下Boltzmann分布的假設，這個假設是否合理，需要進一步研究.另外，采用何種方式對平動溫度和振動溫度取平均以確定化學反應控制溫度主要依據經驗，也具有相當的不確定性.

態－態模型法描述振動非平衡及其與化學反應耦合的方法則更為精細和直接.它將粒子在能級間的躍遷與化學反應視為同一性質的過程研究，能夠得到流體組元詳細能級分布而不只是總的非平衡振動能的變化規律，掌握了更多的細節；它直接研究各能級上粒子的化學反應，避免了多溫度模型中經驗性的熱力－化學耦合方法.不過，盡管態－態模型在理論上展現出諸多優點，但各種內能模式能級間的躍遷過程復雜，使得計算中處理的信息量巨大，并且計算結果強烈依賴于所使用的微觀模型和動力學給出的相關數據，所以目前還無法實現態－態模型與多維流場程序的耦合［5］，國外采用態－態模型的研究一般集中在不涉及流動的零維問題［5］，或是正激波后、邊界層［6－7］和準一維噴管非平衡流動［8－9］問題.

國內在態－態模型的研究方面剛剛起步，一些學者對模型基本原理和實驗方法進行了總結和討論［10］，筆者曾采用態－態模型研究零維問題［11］，對等溫等壓或等溫等容的N2／N混合物系統，詳細研究其熱化學非平衡過程，揭示了非平衡過程中的一些微觀分布信息，加深了對態－態模型和非平衡過程的理解.

進一步將態－態模型與流動耦合，研究N2／N混合物的準一維熱化學非平衡噴管流動.國外學者求解耦合態－態模型的準一維非平衡噴管流控制方程時是采用空間推進方法［8－9］，考慮到非平衡噴管流中喉部非聲速等原因造成的入口條件試湊的困難，本文采用時間推進法.在駐室溫度3 000 K～10 000 K，壓力1 atm～10 atm條件下數值求解N2／N非平衡噴管流動，獲得流場參數分布特別是詳細的振動能級分布，分析噴管流動中的宏觀和微觀特性變化特點.對這類復合占優的流動，考察雙溫度模型中假設的合理性.

1 基本方程與計算方法

態－態模型將位于不同振動能級上的分子視為不同組元，位于振動能級v上的N2記作N2（v）.本文采用有68（0～67）個振動能級的N2能級模型，因此N2（v）／N混合物中含68個分子組元和1個N原子組元，共69個組元.態－態模型直接描述分子在不同振動能級間的躍遷和各能級上分子N2（v）與原子N的化學反應，共5類過程：分子間的振動－平動能量交換過程（vTm）、分子與原子間的振動－平動能量交換過程（vTa）、振動－振動能量交換過程（vv）；原子作為碰撞參與者的離解－復合反應（dra）和分子作為碰撞參與者的離解－復合反應（drm）.具體表示如下［12］

其中v和w代表N2的不同振動能級.根據上述各類躍遷或化學反應的速率方程就可給出N2（v）和N共計69個組元的質量生成率.各類過程的速率系數具體計算公式和相關數據見文獻［6，12－13］.

對非平衡流動，需要69個組元連續方程.加上動量、能量方程后，守恒形式的準一維熱化學非平衡噴管流控制方程為

其中ρi為N2（v）和N各組元的密度，A為噴管截面面積，E為單位體積的總能，˙wi為的各組元的質量生成率，由上述（1）～（5）式給出的各類過程的速率方程表達，其中包含各組元密度ρi．

對控制方程（6）式進行差分離散時，無粘通量項E和源項W采用全隱格式，右端項F采用顯式離散，隱式求解方法為LU-SGS方法，具體參見文獻［14］.源項隱式求解時需要的（?W／?U）矩陣需根據態－態模型各躍遷過程的速率系數公式推導計算［15］.

定解條件中，噴管入口給出總溫總壓條件，流速由內點外插得到，假設駐室至入口第一點流動接近平衡因而等熵.入口點的N原子和各振動能級N2（v）的質量分數采用駐室溫度壓力下的態－態模型平衡結果給定.噴管超聲速出口邊界上所有物理量由內點外插得到.初場的確定首先根據噴管各點膨脹比分布馬赫數，以熱化學平衡條件給出對應馬赫數的宏觀流動參數，N和N2（v）組元的質量分數則取為駐室條件下的態－態模型的平衡結果.

為驗證筆者編制的“態－態模型”數值計算程序，曾對文獻［12］的算例，進行了N2／N系統熱化學非平衡過程的數值模擬，程序結果和文獻結果吻合很好［11］.

2 數值算例與結果分析

噴管流動計算針對一軸對稱的噴管進行.噴管半徑沿軸線的分布曲線為拋物線

噴管長度為1 m，喉部位置在x＝0.5 m處.

噴管駐室溫度T0在3 000 K～10 000 K范圍，駐室壓力p0在1 atm～10 atm范圍，本文選取了多個算例條件求解N2／N非平衡噴管流場，下面以T0＝8 000 K（N2離解程度較高）和T0＝3 000 K（離解程度很低）兩組算例為代表給出數值計算結果和分析.

2.1 T0＝8 000 K算例

采用了態－態模型進行噴管非平衡流計算得到的宏觀流動參數如壓力、溫度、密度、速度在氣流膨脹過程中的變化特點與雙溫度模型得到的結果一致，這里不再給出，下面主要討論組元質量分數或N2（v）與N2的粒子數之比沿軸線的分布，噴管中不同位置的振動能級分布及其偏離Boltzmann分布的特點.

2.1.1 參數沿軸線分布

圖1給出了N質量分數CN沿噴管軸線的分布.p0＝1 atm條件下CN約在0.7～0.8之間，5 atm時約為0.4～0.5，10 atm時約為0.3～0.4.這是總壓升高抑制駐室N2離解，而膨脹過程中提高壓力又促進N復合生成N2的緣故.三個總壓條件下噴管喉道下游都表現出化學凍結的現象，但總壓越大，凍結出現得越晚.

圖2以第0（基態）和第30振動能級為代表給出了N2（v）質量分數CN2（v）沿軸線的分布.各振動能級N2（v）的質量分數在入口到x／L＝0.4這段范圍內變化都較小，直到喉道附近才發生明顯變化，最后在喉道下游的不同位置趨于凍結.喉道前溫度變化較小，各能級的CN2（v）也相對穩定，到喉道附近后溫度快速降低，CN2（v）發生明顯變化.喉道下游氣流膨脹過程中溫度繼續下降，但由于壓力和密度的快速下降，各類躍遷和化學反應速率都急劇下降，而流速又快速增大，因此喉道下游不遠處就表現出流動熱化學凍結的特點.不同振動能級出現凍結的位置略有差別，但大致都在x／L＝0.6～0.8之間.總壓越高流動凍結出現得越晚，這是由于壓力增加（密度增大）增大了粒子碰撞頻率，使流動的非平衡程度降低.振動基態的CN2（0）在膨脹過程中隨溫度降低而增加（見圖2（a）），而其它較高振動能級的CN2（v）變化趨勢相反（見圖2（b））.從CN2（v）凍結值來看，總壓越大，基態及低能級的凍結值越高，而高能級的凍結值越低，這與總壓升高推遲流動凍結的特點是一致的.

圖1 T0＝8 000 K條件下CN沿噴管軸線分布Fig.1 Distribution of CNalong the axis of nozzle at T0＝8 000 K

圖2 T0＝8 000 K條件下CN2（v）沿噴管軸線分布Fig.2 Distribution of CN2（v）along the axis of nozzle at T0＝8 000 K

圖2反映了總壓條件對CN2（v）在膨脹過程中變化規律的影響，但由于不同總壓條件下N2／N混合物的組成不同，質量分數還不能反映不同振動能級上粒子分布.圖3分別給出了基態和第30振動能級的N2（v）粒子數在總的N2粒子數中所占比例沿軸線的分布.由圖3（a）可見，在流動凍結之前，p0＝1 atm條件下N2（0）所占比例最高，p0＝10 atm條件下最低，這對應于總壓較低時氣流溫度較低，更多的粒子通過躍遷進入低振動能級；從N2（0）比例的凍結值來看，則是p0＝10 atm條件下最高，p0＝1 atm條件下最低.這是由于壓力較低條件下能級分布較早地進入凍結狀態.較高振動能級的粒子比例分布的變化特點及總壓條件的影響則與基態能級情況相反，這可由圖3（b）看出.

圖3 T0＝8 000 K條件下N2（v）占N2比例沿噴管軸線分布Fig.3 Distribution of ratio of N2（v）in N2along the axis of nozzle at T0＝8 000 K

2.1.2 噴管不同位置的振動能級分布

圖4 T0＝8 000 K條件下CN2（v）分布與Boltzmann分布Fig.4 Distributions of CN2（v）and Boltzmann distribution at T0＝8 000 K

圖4分別給出了p0＝1 atm、5 atm條件下，噴管軸線位置為x／L＝0.2、0.5、0.8處的各振動能級上的質量分數.圖中實線代表態－態模型計算得到的實際振動能級分布，虛線代表該位置處流場平動溫度下的Boltzmann分布.可見在x／L＝0.2處，兩分布基本重合，說明此處流動仍接近振動平衡，這是由于流動速度還很低，溫度壓力和駐室相比也沒有明顯下降，流動過程中分子間可進行足夠多的碰撞以使振動能級分布達到當地溫度下的平衡態.在x／L＝0.5（喉道）處，p0＝1 atm條件下的兩種分布出現差別，高能級上N2（v）的質量分數要高于當地溫度下的Boltzmann分布值，這表明氣流在膨脹降溫和加速過程中，粒子并沒有及時完成由高能級向低能級的躍遷，流動出現振動非平衡；而在p0＝5 atm的條件下，兩種分布仍基本重合，這是由于壓力增加提高了碰撞頻率，推遲了噴管流動中非平衡的出現.在x／L＝0.8處，實際振動能級分布與當地溫度下Boltzmann分布的差別在p0＝1 atm（圖4（a））、5 atm（圖4（b））和10 atm（圖不再給出）條件下都非常明顯，除基態外，各能級N2（v）質量分數都要高于Boltzmann分布值，并且總壓越低，差別越大.這是因為喉道下游氣流溫度和壓力（密度）的降低與速度的增加都非常迅速，粒子間碰撞頻率明顯降低，非平衡程度進一步加劇，流動甚至趨近于化學和熱力學凍結.

圖4明顯地反映出噴管流動中的振動非平衡現象.雙溫度或多溫度模型通過引入表征非平衡振動能的振動溫度描述振動非平衡，隱含假設振動能級分布滿足振動溫度下的Boltzmann分布.這里對該假設的合理性做進一步分析.首先根據本文采用了態－態模型計算得到的細致振動能級分布和各能級能值求得總的非平衡振動能，再采用和雙溫度模型中類似的辦法折算出振動溫度，然后比較本文得到的振動能級分布和該振動溫度下的Boltzmann分布.

圖5給出了T0＝8 000 K，p0＝1 atm、5 atm和10 atm條件下，平動溫度和振動溫度沿噴管軸線的分布.可見喉道前和喉道后部分區域振動溫度與平動溫度都是基本重和的（高總壓條件下重合區域更大），在氣流繼續流向下游的過程中二者分離，平動溫度仍然有明顯下降，而振動溫度下降緩慢并趨于凍結不變.由此可知引入振動溫度確實能在一定程度上方便地反映出振動非平衡現象，但是否能夠反映出細致的非平衡能級分布還需繼續考察.

圖6對比了T0＝8 000 K，p0＝1 atm、5 atm條件下采用了態－態模型算得的實際振動能級分布和振動溫度下的Boltzmann分布.可見二者有明顯偏離，與圖4非常類似.當然，在x／L＝0.8處，實際振動能級分布與振動溫度下Boltzmann分布的差別要比與平動溫度下Boltzmann分布的差別稍小.這和喉道前及喉道后近下游平動溫度與振動溫度重合而在遠下游處振動溫度凍結在高于平動溫度值的特點（見圖5）是一致的.該現象表明，對實際的非平衡流動，即使引入振動溫度表征非平衡振動能，也不能簡單地用振動溫度下的Boltzmann分布描述細致的能級分布.

圖5 T0＝8 000 K條件下平動溫度和振動溫度沿噴管軸線分布Fig.5 Distributions of translational and vibrational temperature along the axis of nozzle at T0＝8 000 K

圖6 T0＝8 000 K條件下CN2（v）分布與振動溫度下Boltzmann分布Fig.6 Distributions of CN2（v）and Boltzmann distribution at vibrational temperature，T0＝8 000 K

2.2 T0＝3 000 K算例

該算例條件下N2的離解非常微弱.計算中發現N質量分數在喉道后出現振蕩，并略有上升（見圖7）.這可能是因為算例條件下N為微量組元，而68個振動能級上的N2（v）的離解或其逆反應均可能影響到它的質量分數，容易引起數值振蕩.

根據這個情況，對該算例條件另外進行了只計及N2的vTm和vv過程的噴管非平衡流場計算，分析化學反應和vTa過程對宏觀流動特性和微觀粒子能級分布的影響.

圖8對比了計及所有振動躍遷和化學反應與只計及N2的vTm和vv躍遷過程算得的溫度沿噴管軸線分布，其中振動溫度是根據非平衡振動能折算得出.可見化學反應和vTa過程對溫度分布的影響非常有限，除入口附近無化學反應情況下算得的值略高外，其它位置是否計及化學反應得到的結果基本重合.流速、壓力、密度特性也基本不受化學反應和vTa過程影響.

圖7 T0＝3 000 K、p0＝1 atm條件下CN沿軸線分布Fig.7 Distribution of CNalong the axis of nozzle at T0＝3 000 K，p0＝1 atm

圖8 T0＝3 000 K、p0＝1 atm條件下溫度沿軸線分布（包含或不含化學反應）Fig.8 Distribution of temperature along the axis with or without chemical reaction at T0＝3 000 K，p0＝1 atm

但是，是否計及化學反應和vTa過程獲得的粒子能級分布結果有明顯差異.圖9給出了噴管3個位置上的振動能級分布情況.不考慮化學反應時，振動能級分布與當地平動溫度下的Boltzmann分布在喉道前及喉道附近完全重合，在下游偏離Boltzmann分布的程度也比考慮化學反應的情況下弱得多.這主要是因為N與N2高振動能級間的躍遷速率較高，而高振動能級的質量分數本身就很小，存在化學反應時N的復合會對其產生較明顯的影響；而沒有化學反應時，vTm和vv過程會使能級分布較快趨近于平衡態.同時可以看到此時高振動能級的質量分數很小，在10－10以下，所以雖然考慮化學反應與否得到的高振動能級分布有較明顯差異，但基本不影響宏觀流動特性.

圖9 T0＝3 000 K、p0＝1 atm條件下CN2（v）分布與Boltzmann分布Fig.9 Distributions of CN2（v）and Boltzmann distribution at T0＝3 000 K，p0＝1 atm

3 結論

用耦合態－態模型與流動控制方程數值研究準一維熱化學非平衡噴管流，在駐室溫度3 000 K～10 000 K，壓力1 atm～10 atm條件下數值求解了N2／N非平衡噴管流動，分析流動宏觀和微觀特性，得到如下結論：

1）采用時間推進法求解耦合態－態模型的準一維非平衡噴管流，可有效解決非平衡噴管流求解時由于喉道非聲速等原因造成的入口條件試湊困難問題.

2）噴管流動中喉道后振動非平衡現象顯著.雖然引入振動溫度能夠表征非平衡振動能，卻無法細致描述各振動能級上的粒子分布，真實的粒子分布與振動溫度下的Boltzmann分布仍有明顯差別，復合占優的膨脹流中粒子偏離平衡分布的程度比離解占優的激波后流動中更強.膨脹過程中不同振動能級分布的變化特點不同，細致的微觀能級分布規律必須通過采用態－態模型獲得.

3）在較低的總溫條件（如T0＜3 000 K）下N2的離解非常微弱，化學反應對流場宏觀特性的影響可以忽略，因此流動計算時是否考慮化學反應均可得到合理的流場宏觀參數分布.但是否計及化學反應得到的高振動能級分布偏離平衡分布的程度有明顯不同，本文初步分析認為不計及化學反應得到的結果更加合理，而計及微量化學反應的計算結果具有較大的數值不確定性.該問題在未來的研究中可以進一步深入.

4）態－態模型較雙溫度或多溫度模型更具理論優勢，但由于其計算量巨大，目前用于多維非平衡流動模擬還有困難.耦合態－態模型進行準一維噴管流研究，一方面更深入了解了重要非平衡區域微觀過程對宏觀流動特性的影響，可為改進現有的宏觀熱化學模型提供思路，另一方面也為未來將態－態模型用于二維或三維非平衡流動研究做了一定技術儲備.

［1］Park C.Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles［J］.Progress in Astronautics and Aeronautics，1985，96：511－537.

［2］Capitelli M，Armenise I，Gorse C.State-to-state approach in the kinetics of air components under re-entry conditions［J］. Journal of Thermophysics and Heat Transfer，1997，11（4）：570－578.

［3］Silva M L，Guerra V，Loureiro J.State-resolved dissociation rate for extremely nonequilibrium atmospheric entries［J］.Journal of Thermophysics and Heat Transfer，2007，21（1）：40－49.

［4］Park C.The limits of two-temperature model［R］.AIAA 2010－911，2010.

［5］Colonna G，Armenise I，Bruno D，et al.Reduction of state-to-state kinetics to macroscopic models in hypersonic flows［J］．Journal of Thermophysics and Heat Transfer，2006，20（3）：477－486.

［6］Capitelli M，Armenise I，Colonna G，Koudriavtsev N，et al.Nonequilibrium vibrational kinetics during hypersonic flow of a solid body in nitrogen and its influence on the surface heat flux［J］.Plasma Chemistry and Plasma Processing，1995，15（3）：501－527.

［7］Armenise I，Capitelli M，Gorse C.Nitrogen nonequilibrium vibrational distribution and non-Arrhenius dissociation constants in hypersonic boundary layers［J］.Journal of Thermophysics and Heat Transfer，1998，12（1）：45－51.

［8］Colonna G，Tuttafesta M，Capitelle M，et al.NO formation in one-dimensional nozzle air flow with state-to-state nonequilibrium vibrational kinetics［R］.AIAA 98－2951，1998.

［9］Colonna G，Tuttafesta M，Capitelle M，et al.Non-Arrhenius no formation rate in one-dimensional nozzle air flow［J］.Journal of Thermophysics and Heat Transfer，1999，13（3）：372－375.

［10］戴東旭，楊學明.基元化學反應態態動力學研究［J］.化學進展，2007，19（11）：1633－1645.

［11］徐丹，張威，曾明，等.態－態模型下N2／N系統的熱化學非平衡過程研究［J］.空氣動力學學報，2014，32（3）：280－288.

［12］Colonna G，Pietanza L D，Capitelle M.Recombination-assisted nitrogen dissociation rates under nonequilibrium conditions ［J］.Journal of Thermophysics and Heat Transfer，2008，22（3）：399－406.

［13］Esposito F，Armenise I，Capitelli M.N-N2State to state vibrational-relaxation and dissociation rates based on quasiclassicalcalculations［J］.Chemical Physics，2006，331：1－8.

［14］柳軍.熱化學非平衡流及其輻射現象的實驗和數值計算研究［D］.長沙：國防科學技術大學，2004.

［15］徐丹.高溫熱化學非平衡模型研究［D］.長沙：國防科學技術大學，2012.

Numerical Study of Thermochemical Nonequilibrium Nozzle Flow in State-to-State Model

XU Dan，ZENG Ming，ZHANG Wei，LIU Jun
（College of Aerospace Science and Engineering，National University of Defense Technology，Changsha 410073，China）

A state-to-state model is used to simulate thermochemical nonequilibrium nozzle flow of N2／N mixture.Quasi-onedimensional nozzle flow at reservoir temperature from 3 000 K to 8 000 K and pressure from 1 atm to 10 atm is simulated numerically to obtain flow properties，especially detail vibrational population distribution.Evolution of macro-and micro-properties of nonequilbrium flow during expanding is analyzed.Rationality of two-temperature or multi-temperature approach is investigated for recombination dominating flows.

state-to-state model；thermochemical nonequilibrium；vibrational population distribution；hypersonic nozzle flow；numerical simulation

date：2013－10－09；Revised date：2014－02－21

V211.22

A

2013－10－09；

2014－02－21

國家自然科學基金（11102231）資助項目

徐丹（1987－），男，遼寧鳳城，博士生，主要從事計算流體力學、高溫氣體動力學研究，E-mail：449761494＠qq.com

1001-246X（2014）05-0531-08