李志輝，吳俊林，蔣新宇，張涵信

（1.中國空氣動力研究與發(fā)展中心，綿陽 621000；2.國家計算流體力學實驗室，北京 100191）

0 引 言

飛行器往返大氣層先后經歷連續(xù)流、近連續(xù)流、滑移過渡流以至高稀薄自由分子流［1］等，不同流域氣體熱力學性質及繞流狀態(tài)互不相同，特別是位于連續(xù)流與自由分子流之間的過渡區(qū)流動，無論在實驗技術還是數(shù)值方法、理論描述方面均是相當難于處理的流動。在稀薄過渡流區(qū)，由于連續(xù)介質假設失效，N－S方程已不能有效描述飛行器繞流流動現(xiàn)象；而基于微觀分子運動與碰撞隨機統(tǒng)計模擬的DSMC方法因受網格劃分、時間步長限制，難以對低克努森Kn∞數(shù)高空近連續(xù)過渡流進行數(shù)值仿真［2－4］。如何準確模擬近連續(xù)滑移過渡流區(qū)高超聲速氣動力／熱環(huán)境，一直是人們所關心的問題［4－5］，急需發(fā)展相關基礎理論與數(shù)值計算方法。

氣體分子運動論（氣體動理學理論）的基本方程——Boltzmann（玻爾茲曼）方程［6］本身可描述各個流域的氣體分子輸運現(xiàn)象，通過求解Boltzmann方程可以研究各流域氣體動力學問題。但該方程是一個高度非線性積分－微分方程，除Maxwell分布等少數(shù)幾個解析解外，幾乎不可能求出精確解。為避免求解復雜碰撞項所帶來的困難，眾多學者基于守恒不變量，由數(shù)學上較簡單的統(tǒng)計和碰撞松弛模型代替Boltzmann方程碰撞項，提出許多類似Boltzmann方程各階矩的氣體運動論模型方程。最簡單的BGK模型方程于1954年提出［7］，隨后，人們發(fā)展了多種形式的Boltzmann模型方程，諸如 Holway的橢球統(tǒng)計模型［8］、Shakhov基于BGK方程的修正而得到的高階推廣模型方程［9］等。因此，求解從Boltzmann方程簡化而得的氣體運動論模型方程，是一個較為經濟有效的解決稀薄氣體流動問題的途徑?；谖墨I［7－10］基礎，我們［11－15］從研究描述各流域均適用的氣體分子速度分布函數(shù)方程出發(fā)，吸收計算數(shù)學指數(shù)型積分求解原理，提出發(fā)展氣體分子運動論離散速度坐標法，研制經改進的Gauss－Hermite積分方法等系列離散速度數(shù)值積分技術，消除原分布函數(shù)對速度空間的連續(xù)依賴性，將Boltzmann模型方程化為在各個離散速度坐標點處基于時間、位置空間具有非線性源項的雙曲型守恒方程，利用二階Runge－Kutta方法數(shù)值求解碰撞松弛源項方程，引入NND格式數(shù)值求解對流運動項。發(fā)展可用于氣體分子速度分布函數(shù)有限差分求解的氣體運動論邊界條件數(shù)學模型并進行數(shù)值處理。由此建立求解一維、二維、三維Boltzmann模型方程的氣體運動論統(tǒng)一算法（Gas－Kinetic Unified Algorithm，GKUA），并開展跨流域氣體繞流問題模擬應用［16－18］。為了開展近空間飛行器跨流域非平衡繞流問題研究，近年來，基于對轉動自由度松弛變化［19］、氣體分子運動論 Rykov模型［20－21］研究，在氣體分子速度分布函數(shù)演化更新求解中考慮轉動自由度影響，把氣體分子轉動能作為分布函數(shù)的自變量，在統(tǒng)一算法理論計算框架下，提出考慮轉動非平衡效應Boltzmann模型方程數(shù)值算法［22］。本文在此基礎上，針對高、低不同馬赫數(shù)（1.5≤Ms≤25）氮氣一維激波結構和全局 Knudsen數(shù)（9×10－4≤Kn∞≤10）豎直平板、Ramp制動器跨流域繞流、尖雙錐再入體繞流算例，進一步研究含轉動非平衡效應Boltzmann模型方程在全飛行流域繞流問題中的模擬研究，計算剖析飛行器從高稀薄自由分子流到連續(xù)流跨越中間過渡帶非平衡流動變化規(guī)律與復雜繞流現(xiàn)象。

1 控制方程及計算方法

基于轉動松弛特性 Rykov模型［20－21］，采用轉動慣量來描述氣體分子的自旋運動，利用分子總角動量守恒作為一個新的碰撞不變量?；跉怏w分子速度分布函數(shù)f＝f（r，ξ，t，e），這里r、ξ分別是位置空間和速度空間的坐標，e為內能，在求解Boltzmann模型方程統(tǒng)一算法框架［11－18］下，基于權因子1和e對速度分布函數(shù)進行無窮積分，引入約化速度分布函數(shù)f0與f1，發(fā)展含轉動非平衡效應的Boltzmann模型方程，其無量綱形式為：

宏觀流動參數(shù)可以由f0和f1對速度空間積分求得：

粘性系數(shù)μ和碰撞次數(shù)Z可定義為：

這里，T∞和T＊均為有量綱的值，但B＝T∞／T＊是無量綱量［23］。對于氮氣N2而言，T＊＝91.5K，1／δ＝1.55。

速度分布函數(shù)f0、f1是依賴于時間t、位置空間r和速度空間ξ的函數(shù)，需將速度空間離散降維，去掉分布函數(shù)對速度空間的連續(xù)依賴性。氣體運動論離散速度坐標法［10－11，13－14，24］通過一套與積分規(guī)則相一致的離散速度坐標點分布，用具有N個元素的函數(shù)簇代替原分布函數(shù)對積分變量的連續(xù)依賴性，把模型方程轉換為在各個離散速度坐標點處彼此獨立的關于位置空間和時間的雙曲型守恒方程，而宏觀流動參數(shù)則可以通過相應于離散速度坐標法的積分求解規(guī)則動態(tài)確定。如經過速度空間離散化處理的速度分布函數(shù)方程在計算平面內矩陣形式為：

其中，U＝JF，E＝·U，G＝·U，H＝·U，

采用文獻［11，13－14，24］所建立二階有限差分格式，將時間分裂數(shù)值計算方法應用于方程（4）中位置空間的三個對流運動方程和 碰撞松弛源項方程數(shù)值求解，得到在每個離散速度坐標點處直接求解分子速度分布函數(shù)的氣體運動論數(shù)值格式：

2 邊界條件數(shù)值處理方法

對于流場邊界條件數(shù)值處理，按照統(tǒng)一算法基于分子速度分布函數(shù)邊界處理數(shù)學模型［14，17，24］，對于特征速度小于零的無窮遠處來流條件，我們以無量綱的來流宏觀參數(shù)n∞＝1、U∞、V∞、W∞和T∞＝1表征的平衡態(tài)分布函數(shù)確定：

而在流場出口邊界上，認為分布函數(shù)趨于均勻變化而不存在變化梯度，即：

物面邊界條件由氣體分子與物體表面碰撞時，遵循無穿透通量守恒條件［13－15，17］和平動、轉動能量平衡關系，推導得到：

其中，nw為無穿透通量守恒邊界條件確定的物面反射氣體分子數(shù)密度，

物體表面氣體的平動溫度和轉動溫度可表示為：

公式（10）、（11）中的適應系數(shù)αt、、αr和可以有兩種賦值方法。第一種，對應于完全熱適應邊界條件：αt＝1，＝0，αr＝1，＝0。第二種，適應系數(shù)通過分析氮氣中熱傳導實驗數(shù)據(jù)得到［25］：αt＝0.8，＝0.1，αr＝0.5＝0.8。

3 數(shù)值計算結果與討論

3.1 高、低不同馬赫數(shù)定常正激波內流動

正激波相當于一個從均勻超聲速上游到一個均勻亞聲速下游流動過渡［26］。引入無量綱參數(shù)：

為檢驗考慮轉動非平衡影響的Boltzmann模型方程統(tǒng)一算法對高、低不同馬赫數(shù)激波內流動問題模擬能力，計算氮氣中Ms＝1.53、3.2、10、25激波內流動。圖1分別繪出馬赫數(shù)Ms＝1.53、3.2、10、25激波結構無量綱密度ρ＊分布，圖中實線“Cal.”表示統(tǒng)一算法結果，符號“○”為來自文獻［27］實驗數(shù)據(jù)，圖1（a）、（b）中符號“▽”為來自文獻［21］模擬值，圖1（c）、（d）中符號“△”為來自文獻［19］使用廣義退化Boltzmann方程（GBE）計算結果?？煽闯?，對于高、低不同馬赫數(shù)1.53～25激波內流動密度分布統(tǒng)一算法結果均與實驗數(shù)據(jù)、典型文獻模擬值吻合很好。圖1關于弱激波Ms＝1.53到強激波Ms＝25的計算，證實本文方法在計算分析雙原子氮氣高、低不同馬赫數(shù)激波內流動問題方面的準確可靠性。圖2給出了Ms＝7激波內流動平動溫度和轉動溫度分布，其中符號“Ref.”為來自文獻［21］模擬值，圖2（b）標示出激波轉動溫度實驗測試數(shù)據(jù)［28］?？煽闯?，從整個激波結構的變化情況來看，不僅平動溫度而且轉動溫度分布本文計算值均與考慮轉動能自由度影響的Rykov模型結果相一致，且與轉動溫度分布實驗測試數(shù)據(jù)相吻合，驗證了統(tǒng)一算法程序求解考慮轉動非平衡影響的一維流動問題可行性。

圖1 高低不同馬赫數(shù)下激波內流動密度分布計算比較Fig.1 Density distribution in shock wave of various Mach numbers

圖2 Ms＝7激波結構內流動溫度分布計算比較Fig.2 Temperature profiles of shock wave at Ms＝7

3.2 二維再入制動器全飛行流域繞流計算分析

為了檢驗與應用本文發(fā)展的統(tǒng)一算法計算再入不同Knudsen數(shù)二維高超聲速繞流問題可靠性，擬定豎直平板繞流試驗狀態(tài)［29］：來流馬赫數(shù)Ma∞＝8.3666，壁溫與總溫比Tw／T0＝0.6667，γ＝1.4。

表1繪出利用壁面完全熱適應邊界條件計算Kn∞＝0.1、0.01豎直平板繞流阻力與熱傳導系數(shù)與文獻［29］數(shù)據(jù)比較，可看出，在高稀薄流Kn∞＝0.1與近連續(xù)流Kn∞＝0.01跨越不同流區(qū)上述豎直平板繞流物面力／熱系數(shù)GKUA結果均與典型文獻吻合很好，二者最大偏差不超過4.9%。

圖3、圖4分別繪出Kn∞＝0.01平板繞流流場平動溫度Tt、轉動溫度Tr等值線分布與文獻［29］數(shù)據(jù)比較，可看出兩種結果變化趨勢相當吻合。對此近連續(xù)流場，平板繞流前端的減速區(qū)內出現(xiàn)明晰的脫體激波，由于Knudsen數(shù)不是小量，流場稀薄效應明顯存在，平板前端脫體激波較厚。圖5描繪了該平板繞流Kn∞＝0.01、Ma∞＝8.3666流線分布，給出了氣體繞流方向與流動變化特點，并與圖5（b）顯示的文獻［29］中的流線圖進行了對比。結果表明，GKUA流線分布（圖5a）與文獻［29］計算結果之間具有較好的一致性。對此近連續(xù)繞流Kn∞＝0.01狀態(tài)，氣體流經平板，在物體前部一定區(qū)域因跨越脫體激波而流向改變，繞過物體最高點，進入背風區(qū)，出現(xiàn)流動分離、回流渦區(qū)。從圖3～圖5統(tǒng)一算法計算得到的流場結構與文獻［29］結果定性化比較看出，GKUA結果流場分辨率較好，加之表1物面力／熱系數(shù)定量化比較，進一步證實本文考慮轉動非平衡效應的Boltzmann模型方程統(tǒng)一算法計算再入不同流區(qū)高超聲速流動問題準確可行性。

表1 跨流區(qū)豎直平板繞流阻力和熱傳導系數(shù)比較Table 1 Comparison of drag and heat transfer coefficients past erectplate covering various flow regimes

圖3 平板繞流Kn∞＝0.01平動溫度Tt等值線分布計算比較Fig.3 Translational temperature contours in erect plate flow with Kn∞＝0.01

圖4 平板繞流Kn∞＝0.01轉動溫度Tr等值線分布計算比較Fig.4 Rotational temperature contours in erect plate flow with Kn∞＝0.01

圖5 平板繞流Kn∞＝0.01、Ma∞＝8.3666流場流線分布Fig.5 Streamline structures past erect plate with Kn∞ ＝0.01

針對如圖6所示的一個平板與25°傾角斜板構成二維Ramp再入制動器外形，計算Ma∞＝10的高馬赫數(shù)繞流流場。設置：克努森數(shù)Kn∞＝6×10－4、來流氣體溫度T∞＝52K、密度ρ∞＝3.91×10－4kg／m3、速度V∞＝1477m／s、壁面溫度Tw＝290K。

圖6 Ramp制動器外形及尺寸（單位：cm）Fig.6 Ramp body geometries（unit：cm）

本文計算得到該Ramp制動器繞流流場馬赫數(shù)等值線與來自文獻［30］的Boltzmann－ES橢球統(tǒng)計模型結果比較，見圖7所示，可看出二者吻合較好，流場變化規(guī)律較為一致，計算再現(xiàn)了上表面水平板前端尖前緣產生斜激波與斜板上部壓縮斜激波相互干擾繞流現(xiàn)象，得到了較高分辨率的流場計算結果，驗證了本文方法程序計算的正確性。

該狀態(tài)下Ramp外形繞流上部流場x＝0.1m和x＝0.2m截線處氣流密度、溫度與馬赫數(shù)隨z軸橫向分布本文計算與文獻［30］模擬結果比較情況，見圖8、圖9所示?？煽闯鲈诹鲌霾煌恢帽疚挠嬎阒蹬c文獻［30］橢球統(tǒng)計模型結果變化趨勢吻合較好，揭示了Ramp凹面體外形前、后部流場出現(xiàn)不同強度壓縮斜激波流動變化規(guī)律，驗證了含轉動非平衡效應Boltzmann模型方程數(shù)值算法對二維高超聲速連續(xù)流區(qū)Kn∞＝6×10－4、Ma∞＝10繞流問題計算適應性。

圖8 Ramp上部流場x＝0.1m截線處密度、溫度與馬赫數(shù)橫向分布Fig.8 Density（ρ），temperature（T）and Mach（M）distribution in the upside field of Ramp at x＝0.1m

圖9 Ramp上部流場x＝0.2m截線處密度、溫度與馬赫數(shù)橫向分布Fig.9 Density（ρ），temperature（T）and Mach（M）distribution in the upside field of Ramp at x＝0.2m

圖10繪出該二維Ramp制動器上表面壓力系數(shù)Cp和摩阻系數(shù)Cf沿x軸線分布本文計算值（Cal.）與文獻［30］中Boltzmann－ES模型結果（Ref.）比較情況，可看出兩者相當吻合，證實所發(fā)展含轉動非平衡效應Boltzmann模型方程統(tǒng)一算法對飛行器物面氣動特性計算的可靠性與精度。

圖10 Ramp上表面壓力、摩阻系數(shù)分布計算比較Fig.10 Distribution of pressure and frictional drag coefficients on the upper surface of Ramp

為了檢驗求解轉動非平衡效應Boltzmann模型方程統(tǒng)一算法對二維Ramp制動器外形再入各流域繞流問題計算能力，擬定高稀薄自由分子流到連續(xù)介質繞流狀態(tài)Kn∞＝10～0.0009、Ma∞＝1.8進行計算分析。圖11給出了跨流區(qū)Kn∞＝10、0.8、0.01、0.0009四個不同克努森數(shù)下計算得到Ramp制動器再入繞流流場尾部矢量流線與流動局部結構?？煽闯觯瑢τ贙n∞＝10的高稀薄自由分子流，氣流完全附著物面流動；而當流動進入稀薄過渡流，如Kn∞＝0.8，在飛行器尾部上端尖角后部繞流區(qū)開始出現(xiàn)較小的回流渦結構；對于Kn∞＝0.01的近連續(xù)流區(qū)，飛行器尾部繞流出現(xiàn)雙渦結構，Kn∞＝0.0009的連續(xù)流區(qū)，尾部流場又褪變?yōu)榇嬖谝粋€巨大的回流渦流動結構，該圖直觀展示了二維制動器再入過程跨越飛行各流域所出現(xiàn)的不同繞流變化規(guī)律。

圖11 Ramp制動器再入過程Kn∞＝10、0.8、0.01、0.0009繞流結構（Ma＝1.8，CFL＝2）Fig.11 Flow structures of Ramp re－entering for all flow regimes with Kn∞ ＝10，0.8，0.01and 0.0009（Ma＝1.8，CFL＝2）

3.3 雙維外形體再入稀薄流區(qū)繞流計算分析

為了檢驗考慮轉動非平衡影響的Boltzmann模型方程統(tǒng)一算法對三維飛行器再入繞流問題的模擬能力，擬定雙錐再入體外形及尺寸如圖12所示，為了與適于高稀薄流模擬的GBE解算器［31］計算比較，擬定來流馬赫數(shù)Ma∞＝3、Kn∞＝1尖雙錐繞流狀態(tài)，圖13繪出統(tǒng)一算法計算該狀態(tài)駐點線密度分布與文獻［31］結果（符號“Δ”）對比情況，可看出兩種方法計算結果很好一致，證實考慮轉動非平衡效應的統(tǒng)一算法用于計算三維再入飛行器繞流問題可行性。

圖14繪出Ma∞＝3繞流Kn∞＝1和Kn∞＝0.1尖雙錐流場密度等值線分布，圖中看出，這兩種狀態(tài)克努森數(shù)較高，流動均呈現(xiàn)較強的稀薄流特征，不存在明顯的激波，流動參數(shù)分布沒有強間斷現(xiàn)象，而且圖14（a）對應Kn∞＝1的高稀薄繞流物體前方擾動區(qū)域較圖14（b）對應更低克努森數(shù)Kn∞＝0.1的稀薄過渡區(qū)繞流要大得多；Kn∞＝0.1對應的氣體流動趨近物體前端出現(xiàn)了激波強擾動現(xiàn)象。

圖12 雙錐再入體外形Fig.12 Bicone geometries

圖13 駐點線密度分布Fig.13 Density distribution along the stagnation line for Ma∞＝3，Kn∞＝1

圖14 不同克努森數(shù)尖雙錐Ma＝3繞流場密度分布Fig.14 Density distribution around bicone with different Kn

圖15直觀顯示Kn∞＝1和Kn∞＝0.1兩種狀態(tài)下尖雙錐繞流流線，反映出兩個狀態(tài)的繞流面貌呈現(xiàn)明顯不同的流動結構，Kn∞＝1對應的高稀薄近自由分子流，氣流完全附著物面流動過錐體后，平滑過渡到背風下游流動狀態(tài)，呈現(xiàn)明顯的高稀薄流特征；而Kn∞＝0.1對應的更低克努森數(shù)稀薄過渡區(qū)繞流，則在錐體后端面背風區(qū)開始出現(xiàn)流動分離，形成上下兩個回流渦，鞍結點拓撲結構，開始呈現(xiàn)出近連續(xù)流繞流面貌。

圖15 不同克努森數(shù)尖雙錐繞流（Ma∞＝3）流線結構Fig.15 Streamline distribution around bicone at different Kn with Ma∞＝3

為了剖析跨流區(qū)流動的非平衡效應，圖16、圖17分別繪出了Kn∞＝1、0.1兩種狀態(tài)對應的繞流流場溫度分布。在Kn∞＝1對應的高稀薄流區(qū)，平動溫度與轉動溫度分布彼此差別很大，計算得到的物體前部繞流區(qū)平動溫度最大值約為3.2，而轉動溫度最大為2.4左右；對應Kn∞＝0.1的稀薄過渡區(qū)繞流，平動溫度與轉動溫度之間的差異明顯減小，物體前部高溫區(qū)平動溫度最大值下降為2.8，轉動溫度升高約為2.65；而從圖3、圖4所繪豎直平板繞流看出，對于Kn∞＝0.01的近連續(xù)繞流，流場平動溫度與轉動溫度分布趨于一致，流場等值線結構差別不大，流動屬于近平衡態(tài)流動。這體現(xiàn)了高Kn∞數(shù)稀薄過渡流區(qū)存在嚴重的流動非平衡，且克努森數(shù)越大非平衡效應越明顯，說明了在稀薄流區(qū)考慮轉動非平衡效應的必要性。

圖16 Ma∞＝3，Kn∞＝1尖雙錐繞流溫度分布Fig.16 Temperature contours around bicone at Ma∞＝3，Kn∞＝1

圖17 Ma∞＝3，Kn∞＝0.1尖雙錐繞流溫度分布Fig.17 Temperature contours around bicone at Ma∞＝3，Kn∞＝0.1

4 結 論

基于對轉動自由度松弛變化、氣體分子運動論Rykov模型研究，應用稀薄流到連續(xù)流的氣體運動論統(tǒng)一算法數(shù)值求解考慮轉動非平衡影響的Boltzmann模型方程。針對高、低不同馬赫數(shù)定常正激波結構、從自由分子流到連續(xù)流二維豎直平板、Ramp制動器及三維尖雙錐外形體再入稀薄流區(qū)繞流問題計算分析，得到如下結論：

（1）通過對弱激波到強激波1.53≤Ms≤25考慮轉動非平衡影響的氮氣激波結構內流動計算結果與實驗數(shù)據(jù)、典型文獻模擬值比較分析，驗證了含轉動非平衡效應Boltzmann模型方程統(tǒng)一算法對求解強、弱激波一維流動問題的準確可靠性。

（2）通過對二維豎直平板、Ramp再入制動器從自由分子流到連續(xù)流跨流域繞流流場轉動非平衡效應模擬研究，驗證了統(tǒng)一算法對全飛行流域模擬計算的一致適用性。

（3）通過對三維尖雙錐外形再入稀薄氣體繞流狀態(tài)數(shù)值計算，對比分析繞流流場平動溫度與轉動溫度分布，揭示了來流Kn∞數(shù)越高，過渡區(qū)非平衡效應越嚴重，飛行器繞流流動非平衡影響越明顯。證實考慮轉動非平衡影響的Boltzmann模型方程統(tǒng)一算法對三維復雜飛行器再入繞流問題模擬能力。

本文工作僅是含轉動非平衡效應Boltzmann模型方程統(tǒng)一算法在跨流域繞流問題方面的初步應用與計算檢驗，將其應用于跨流區(qū)復雜高超聲速氣動熱力學問題以及考慮轉動能與振動能、混合氣體非平衡效應適應性研究，均有待進一步研究。

致謝：本文部分三維繞流計算在國家超級計算天津中心進行，特此感謝。

［1］TSIEN H S.Superaerodynamics，mechanics of rarefied gases［J］.JournalofAeronauticsScience，1946，13（12）：653－664.

［2］KOPPENWALLNER G，LEGGE H.Drag of bodies in rarefied hypersonic flow［R］.AIAA Paper 85－0998.Progress in Astronautics and Aeronautics：Thermophysical Aspects of Reentry Flows.New York，1985，103：44－59.

［3］LI Z H，WU Z Y.DSMC simulation on rarefied aerodynamics of Apollo－CM［J］.ACTAAerodynamicaSinica，1996，14（2）：230－233.（in Chinese）李志輝，吳振宇.阿波羅指令艙稀薄氣體動力學特征的蒙特卡羅數(shù)值模擬［J］.空氣動力學學報，1996，14（2）：230－233.

［4］YOON S，GNOFFO P A，WHITE J A，et al.Computational challenges in hypersonic flow simulations［R］.AIAA 2007－4265.

［5］ZHUANG F G，CUI E J，ZHANG H X.Some development of future spacecrafts and aerodynamics tasks［C］.Proc.of First Aerodynamics and Aerothermodynamics，2006，8：1－12.（in Chinese）莊逢甘，崔爾杰，張涵信.未來空間飛行器的某些發(fā)展和空氣動力學的任務［C］.中國第一屆近代空氣動力學與氣動熱力學會議論文集，2006，8：1－12.

［6］CHAPMANN S，COWLING T G.The mathematical theory of non－uniform hases［M］.3rded.Cambridge，1990.

［7］BHATNAGAR P L，GROSS E P，KROOK M.A model collision processes in gases［J］.Phys.Rev.，1954，94：511－525.

［8］HOLWAY Jr L H.New statistical models for kinetic theory，methods of construction［J］.Phys.Fluids，1966，9：1658－1673.

［9］SHAKHOV E M.Kinetic model equations and numerical results［M］.In：Oguchi H Ed.Proceedings of 14thInternational Symposium on Rarefied Gas Dynamics，Tokyo：University of Tokyo Press，1984，pp137－148.

［10］YANG J Y，HUANG J C.Rarefied flow computations using nonlinear model Boltzmann equations［J］.J.ofComput.Phys.，1995，120：323－339.

［11］LI Z H，ZHANG H X.Pilot study on gas－kinetic unified numerical algorithm for flows from rarefied transition to continuum［J］.ACTAAerodynamicaSinica，2000，18（3）：251－259.（in Chinese）李志輝，張涵信.稀薄流到連續(xù)流的氣體運動論統(tǒng)一數(shù)值算法初步研究［J］.空氣動力學學報，2000，18（3）：251－259.

［12］LI Z H，ZHANG H X.Study on gas kinetic unified algorithm for flows from rarefied to continuum［J］.ACTAAerodynamica Sinica，2003，21（3）：255－266.（in Chinese）李志輝，張涵信.稀薄流到連續(xù)流的氣體運動論統(tǒng)一算法研究［J］.空氣動力學學報，2003，21（3）：255－266.

［13］LI Z H，ZHANG H X.Study on gas kinetic algorithm for flows from rarefied transition to continuum using boltzmann model equation［J］.ActaMechanicaSinica，2002，34（2）：145－155.（in Chinese）李志輝，張涵信.稀薄流到連續(xù)流的氣體運動論模型方程算法研究［J］.力學學報，2002，34（2）：145－155.

［14］LI Z H，ZHANG H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］.J.ofComput.Phys.，2004，193：708－738.

［15］LI Z H.Gas－kinetic unified algorithm for re－entering complex problems covering various flow regimes by solving Boltzmann model equation［M］.Advances in Spacecraft Technologies.In－Tech Publisher，2011：273－332.

［16］LI Z H，ZHANG H X，F(xiàn)U S，et al.A gas kinetic algorithm for flows in microchannel［J］.InternationalJournalofNonlinear SciencesandNumericalSimulation，2005，6（3）：261－270.

［17］LI Z H，ZHANG H X.Gas－kinetic numerical studies of threedimensional complex flows on spacecraft re－entry［J］.J.Comput.Phys.，2009，228：1116－1138.

［18］LI Z H，ZHANG H X.Study on the unified algorithm for threedimensional complex problems covering various flow regimes using Boltzmann model equation［J］.Sci.China－Phys.Mech.＆Astron.，2009，52（1）：124－138.（in Chinese）李志輝，張涵信.基于Boltzmann模型方程各流域三維復雜繞流問題統(tǒng)一算法研究［J］.中國科學，2009，39（3）：414－427.

［19］AGRWAL R K，CHEN R，CHEREMISIN F G.Computation of hypersonic shock wave flows of a diatomic using generalized Boltzmann equation［R］.AIAA 2007－4541.39th AIAA Thermophysics Conference.Miami，2007.

［20］RYKOV V A.Model kinetic equation of a gas with rotational degrees of freedom［J］.FluidDynamics，1975，10：959－966.

［21］RYKOV V A，TITAREV V A，SHAKHOV E M.Shock wave structure in a diatomic gas based on a kinetic model［J］.FluidDynamics，2008，43（2）：316－326.

［22］LI Z H，WU J L，PENG A P，ZHANG H X.Study of the unified algorithm solving Boltzmann model equation with the effect of rotational non－equilibrium［C］.The 8th Cross－straits Conference on Aeronautics and Astronautics.Beijing，2012.（in Chinese）李志輝，吳俊林，彭傲平，等.考慮轉動非平衡效應Boltzmann模型方程統(tǒng)一算法應用研究［C］.第八屆海峽兩岸航空航天學術研討會論文集，北京，2012.

［23］LARINA I N，RYKOV V A.Diatomic gas flow past a sphere based on kinetic equations［J］.Dokl.Akad.NaukSSSR，1976，227：60－62.

［24］LI Z H.Study on the unified algorithm from rarefied flow to continuum［D］.［PhD thesis］.Mianyang：China Aerodynamics Research and Development Center，2001.（in Chinese）李志輝.從稀薄流到連續(xù)流的氣體運動論統(tǒng)一數(shù)值算法研究［D］.［博士學位論文］.中國空氣動力研究與發(fā)展中心，2001.

［25］LARINA I N，RYKOV V A.The boundary condition on the body surface for a diatomic gas［C］.Proceedings of 15thInternational Symposium on Rarefied Gas Dynamics，Stuttgart，1986，1：635－643.

［26］LI Z H.Applications of gas kinetic unified algorithm for flows from rarefied transition to continuum［D］.［Post－doctor disser－tation］.Beijing：Tsinghua University，2003.（in Chinese）李志輝.稀薄流到連續(xù)流氣體流動問題統(tǒng)一算法應用研究［D］.［博士后研究報告］.清華大學，2003，9.

［27］ALSMEYER H.Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam［J］.J.Fluid.Mech，1976，74：497－513.

［28］ROBBEN F，TALBOT L.Experimental study of the rotational distribution function of nitrogen in a shock wave［J］.Physicsof Fluids，1966，9：653－662.

［29］RYKOV V A，TITAREV V A，SHAKHOV E M.Numerical study of the transverse supersonic flow of a diatomic rarefied gas past a plate［J］.ComputationalMathematicsandMathematical Physics，2007，47（1）：136－150.

［30］ANDRIES P，BOURGAT J F.Numerical comparison between the Boltzmann and ES－BGK models for rarefied gases［J］.ComputerMethodsinAppliedMechanicsandEngineering，2002，191（31）：3369－3390.

［31］WILSON C D，AGARWAL R K，TCHEREMISSINE F G.Computation of hyersonic flow of a diatomic gas in rotationa nonequilibrium past 3Dblunt bodies using the generalized Boltzmann equation［R］.AIAA 2009－3836.