李志輝，方 明，唐少強

（1．中國空氣動力研究與發(fā)展中心超高速研究所，四川 綿陽 621000；

2．北京大學工學院，北京 100871；3．國家計算流體力學實驗室，北京 100191）

0 引 言

基于分子運動與碰撞隨機統(tǒng)計模擬的直接模擬Monte Carlo（DSMC）方法［1］自1963年 G．A．Bird將其發(fā)展用于稀薄氣體流動模擬以來，經過四十多年的發(fā)展，已在求解稀薄過渡領域飛行器繞流中取得巨大成功［2－9］，從復雜高超聲速氣動力／熱流動模擬到細觀速度分布函數的非平衡表征均得到了實驗的驗證與支持?；?對 DSMC 方 法 應 用 研 究 實 踐［7，10－11］體會到，該方法因其概率統(tǒng)計模擬實質，其誤差來自兩個方面，一是計算中的網格尺寸和時間步長，二是統(tǒng)計抽樣中的漲落。如何控制稀薄氣體動力學DSMC方法統(tǒng)計噪聲，提高計算效率與精度，一直是人們研究感興趣的問題［12－16］。文獻［13］中分析了輸運系數和網格尺寸之間的關系，其結論是網格尺度不應大于氣體分子平均自由程。文獻［14］分析了輸運系數和時間步長之間的關系，其結論是時間步長不應大于氣體分子平均碰撞時間尺度。這些與實際模擬計算的經驗相吻合。關于DSMC方法模擬中抽樣的統(tǒng)計漲落，文獻［12］給出了統(tǒng)計量收斂速度的一般規(guī)律，即收斂速度為樣本容量 N 的O（N－1／2）。文獻［15］中給出了傅里葉熱流計算中溫度的統(tǒng)計漲落，文獻［16］從統(tǒng)計力學角度分析了統(tǒng)計量的漲落現象。在DSMC方法中，最重要的兩個統(tǒng)計宏觀量是代表動量的速度和代表能量的溫度，本文嘗試從數理統(tǒng)計角度對這兩個量的統(tǒng)計噪聲進行初步分析，給出控制和減輕DSMC模擬結果統(tǒng)計漲落的方法和途徑。

1 DSMC方法

氣體分子運動論的基本控制方程是Boltzmann方程［17］：

其考察的是微觀粒子的速度分布函數f基于位置空間、速度空間在任一時刻由非平衡態(tài)向平衡態(tài)的演化［17－18，20］，這 是 一 個 高 度 非 線 性 高 維 積 分 － 微 分 方程。該方程的理論解只在極其簡單的情況下存在，一般問題的數值求解一直是研究的熱點問題［18－23］。

DSMC方法巧妙地利用上述方程的碰撞松弛演化特征，從微觀分子運動論概率統(tǒng)計原理出發(fā)，將分子的運動和它們之間的碰撞解耦，用有限個仿真分子代替大量的真實氣體分子，仿真分子的位置坐標、速度分量以及內能被記錄下來，其值因仿真分子的運動、與邊界的相互作用以及仿真分子間碰撞而隨時間改變，最后通過統(tǒng)計網格內仿真分子運動狀態(tài)實現對真實氣體流動的模擬［1－11］。該方法盡管不是對Boltzmann方程的離散與直接求解，但其與Boltzmann方程關于分子混沌與二體碰撞的假設是一致的，文獻［24］研究結果表明，當模擬分子數趨于無窮時，DSMC方法得到的粒子分布函數收斂到Boltzmann方程的一種修正形式。DSMC方法是對Boltzmann方程或其模型方程的一種統(tǒng)計模擬，為了保證隨機統(tǒng)計模擬結果的真實性，該方法要求用于分子碰撞取樣的物理空間網格尺度小于氣體分子平均自由程，且用于將仿真分子運動與碰撞解耦的時間間隔Δt小于平均碰撞時間。該方法得到的是模擬分子的微觀信息，人們通常關心的宏觀流動量可以從微觀信息的統(tǒng)計平均得到。如宏觀速度可定義為

這里的括號〈〉表示系綜平均；記分子熱運動速度c＝v－V0。宏觀溫度可定義為

這里的c表示熱運動速度的大小。

上述宏觀量均定義為系綜平均，實際的DSMC模擬過程中，如果每個網格里的仿真分子數很有限，不可避免地導致統(tǒng)計上的噪聲波動，以致掩蓋有用信息。從數學角度分析統(tǒng)計噪聲對統(tǒng)計量的影響程度并提出控制統(tǒng)計噪聲的理論規(guī)則與途徑，正是本文希望探索的工作。

2 統(tǒng)計噪聲的數學分析

在流場的DSMC方法模擬中，總是先讓流場演化足夠長的時間，使其充分發(fā)展，然后開始做統(tǒng)計抽樣。為了減小統(tǒng)計漲落的影響，一般采用多次統(tǒng)計求平均值的方法。在統(tǒng)計的過程中，各網格中模擬分子的數目變化不大。為了數學分析的方便，假設每個網格中模擬分子的數目保持不變，對應第j個網格中模擬分子的個數為Nj。根據平衡氣體統(tǒng)計理論［12，25］，第j個網格中模擬子的速度分量ui、vi及wi應該滿足正態(tài)分布體的宏觀流動速度，為該網格處的氣體流動溫度，k為玻爾茲曼常數，m為分子質量。

2．1 噪聲、統(tǒng)計相對漲落及置信度的定義

上面所定義的相對統(tǒng)計漲落ε是針對流場的每一個網格處。對于整個流場，我們定義統(tǒng)計量兩種形式的全局相對統(tǒng)計漲落。第一種形式是全局最大相對統(tǒng)計漲落Emax

其反映的是整個流場中噪聲影響的極大值。第二種形式是全局平均相對漲落Eave

這里F表示整個流場區(qū)域，V（F）為區(qū)域F的體積，其反映的是噪聲對于整個流場之影響的總體平均效應。可以根據流場特征，選擇適當的相對統(tǒng)計漲落來表示。

由于統(tǒng)計量存在隨機性，統(tǒng)計量的相對統(tǒng)計漲落是基于概率意義下表述的。為此，我們定義統(tǒng)計量在相對統(tǒng)計漲落ε下的置信度［25］α滿足

其中α為0＜α＜1，稱1－α為置信水平，∑為在此置信水平下相對統(tǒng)計漲落的上界。

2．2 速度的相對統(tǒng)計漲落

第n次抽樣時，第j個網格處某方向上（比如y方向）的宏觀速度Vj，n由式（2）給定為

其中，vi是網格中第i個模擬分子在該方向上的速度。模擬共進行N次抽樣，得到的第j個網格處該方向上的宏觀流動速度Vj定義為

此處U為參考速度。若取參考速度U為該網格處氣體流動的當地速度j，且鑒于第j個網格處的聲速ac滿足

即有

其中，Ma為第j個網格處氣體流動當地馬赫數，γ為氣體的絕熱指數。

結合式（4）、（7）及（12），可得

上式表明，∑越大越好。相對統(tǒng)計漲落允許波動的范圍越大，上限∑越大，其置信水平便越高。另一方面，從概率論角度來說，∑是在一定置信水平下相對統(tǒng)計漲落的上限，則越小越好。平衡二者，上式取等號。這樣，在置信水平不低于1－α的條件下，網格j處速度Vj的相對統(tǒng)計漲落為

該公式意味著，當取參考速度U為網格j處當地氣流速度ˉVj時，速度的相對統(tǒng)計漲落與流動馬赫數成反比，與抽樣次數成平方根反比例關系。

若（10）式中U取為一般參考速度，則有

即相對統(tǒng)計漲落與參考速度成反比，與流場溫度的平方根成正比，且與抽樣次數的平方根成反比。

2．3 溫度的相對統(tǒng)計漲落

同上，第n次抽樣時，第j個網格處的宏觀溫度Tj，n由式（3）給定為

通過N次抽樣，第j個網格處的宏觀溫度Tj定義為雪夫不等式，同樣可以得到，若取參考溫度為流場溫度，此時，在置信水平為1－α下，溫度Tj的相對統(tǒng)計漲落為

即相對統(tǒng)計漲落除與抽樣次數有關外，還與流場溫度成正比、與參考溫度成反比。

需要指出的是，公式（14）和（18）的結論較為簡單，但其參考值為流場參數的準確值，實際計算不便操作；式（15）和（19）結論形式復雜些，但參考值的選取根據流場情況易于實現。鑒于流場中某些區(qū)域速度很小，尤其是相當低馬赫數流動問題，需要控制DSMC模擬過程中速度的相對統(tǒng)計漲落量，可采用（15）式為模擬準則。對于速度較均勻變化的高溫非平衡流場，需要考察流場溫度的相對統(tǒng)計漲落，可采用（19）式。一般來說，相對統(tǒng)計漲落變化較大，實際應用時，通常考察最大相對統(tǒng)計漲落Emax或全局平均相對漲落Eave為宜。計算中，主要通過選擇適當的參考值，控制Emax、Eave的大小，可以采用兩種方式，一是增大網格內模擬分子的個數，然而這種方式又會大大增加計算內存為代價，二是增大抽樣統(tǒng)計的次數（增大N），通過增加計算時間來實現，這種方式是降低DSMC統(tǒng)計漲落更適合選取的途徑。

該公式意味著，溫度的相對統(tǒng)計漲落只與統(tǒng)計抽樣次數的平方根成反比，而與流場的宏觀參數無關。

若取參考溫度為Tref，則溫度的相對統(tǒng)計漲落為

3 算例與結果分析

3．1 Couette流動計算分析

Couette流動是一種較為簡單而經典的邊界值問題，兩個相對運動的無限大平行平板相距H、置于相同溫度Tw，向相反方向沿各自平面分別以Uw速度運動。取兩平板距離的中央為坐標原點，x軸為運動方向，y軸與之垂直。以具有溫度Tref＝273K、壓力Pref＝0．01atm的氬氣作為模擬氣體，假設Tw＝Tref，變徑軟球（VSS）分子模型ω＝0．81、α＝1．4被用來描述分子間相互作用。為檢驗上述結論，選擇式（14）、（15）作為控制DSMC方法模擬結果相對統(tǒng)計漲落的準則關系式，其中U取為平板速度。需要說明的是，在模擬過程中，控制每個網格中的模擬分子個數Nj是比較困難的，通常的做法是控制模擬的流場分子總數Nt，可以認為這兩者是等效的。根據上述公式，我們通過設置板速Uw、模擬分子總數Nt和抽樣次數N的變化來觀察DSMC方法模擬得到的Couette流動速度的統(tǒng)計噪聲變化特點。圖1～圖4分別繪出Kn＝0．1128情況下，不同板速、不同模擬分子數與抽樣次數得到的Couette流動上半槽道內無量綱流動速度沿板間距離分布比較。由式（15）知板速Uw對DSMC模擬結果相對統(tǒng)計漲落的影響是反比例關系，圖1中的速度分布可以很好地說明這一點。圖1（a）繪出固定模擬分子總數Nt＝1000和抽樣次數 N＝20000，不同板速Uw＝0．025 Ma、0．05 Ma、0．1 Ma、0．2 Ma引起的 Couette流動速度分布模擬結果，相互比較看出，在固定Nt和N 的情況下，Couette流動速度相對統(tǒng)計漲落即模擬精度高低與板速Uw的大小關系密切，Uw越低，DSMC結果統(tǒng)計漲落越大，如Uw＝0．025 Ma情況，Couette流動速度的相對統(tǒng)計漲落甚為嚴重，幾乎掩蓋真實的氣體流動速度分布；而Uw增大，同樣情況下得到的DSMC結果統(tǒng)計漲落明顯減小，如Uw＝0．2 Ma對應的Couette流動速度呈現相當光滑分布，見圖中虛線表示；由于不同計算條件下相對統(tǒng)計漲落ε形式表現為雜亂無章，為了定量刻畫其變化幅度，圖1（b）繪出給定上述Nt和N值條件下，不同板速引起的Couette流動速度的全局最大相對統(tǒng)計漲落Emax和全局平均相對統(tǒng)計漲落Eave，該圖是對數圖像，其中虛線由不同情況下Emax擬合得到、點線由Eave擬合得到，結果表明Emax和Eave均在斜率為－1的擬合直線附近散布。該圖還說明，Emax與Eave之間差一個常數量，控制Emax相對于控制Eave更為困難。如若把相對統(tǒng)計漲落控制在1%以下，在Nt＝1000和N＝20000的條件下，采用控制Eave的辦法，對于板速大于0．1 Ma的Couette流動就能滿足要求，而式（15）表明若采用控制Emax的辦法，只有板速大于0．4 Ma的Couette流動才能達到要求。圖2（a）繪出固定板速Uw＝0．025 Ma和抽樣次數N＝20000，使用不同模擬分子總數Nt＝1000、4000、16000、64000得到的結果，圖中看出，在固定板速Uw和N 的情況下，模擬分子總數Nt越少，如Nt＝1000，DSMC結果統(tǒng)計漲落越大，對此較低板速Uw＝0．025 Ma引起的Couette流動，只有在模擬分子總數增加到一定程度Nt＝64000時，DSMC模擬得到的流動速度才較為光滑可行；圖2（b）給出相應的全局最大相對統(tǒng)計漲落Emax和全局平均相對統(tǒng)計漲落Eave，該圖也是對數圖像，結果表明Emax和

圖1 不同Uw引起Couette流動速度DSMC模擬值與相對統(tǒng)計漲落變化關系Fig．1 Velocity distribution of Couette flow and relative statistical fluctuations under different Uwfrom DSMC simulation

圖2 不同模擬分子數得到的Couette流動速度DSMC模擬值與相對統(tǒng)計漲落Fig．2 Velocity distribution of Couette flow and relative statistical fluctuations under different Ntfrom DSMC simulation

圖3 不同抽樣次數得到的Couette流動速度DSMC模擬值與相對統(tǒng)計漲落Fig．3 Velocity distribution of Couette flow and relative statistical fluctuations under different Nfrom DSMC simulation

圖4 不同模擬條件下得到的Couette流動速度DSMC模擬值與統(tǒng)一算法結果比較Fig．4 DSMC simulated velocity distribution under different conditions with the comparison of results from the unified algorithm

Eave均在斜率為－1／2的直線附近散布，兩者之間也是差一個常數量。從圖中可以看出，對此Uw＝0．025 Ma的Couette流動模擬，若固定N＝20000，則Nt增加到Nt＝32000時在95%的置信水平下，Eave就可控制在1%以下；如若控制Emax在1%以下，根據公式（15），此時模擬分子總數Nt須在105以上，這需要比控制Eave多得多的計算內存，顯示出本文開展DSMC方法統(tǒng)計噪聲控制規(guī)則研究的現實意義。關于速度的相對統(tǒng)計漲落公式（14）、（15）也說明統(tǒng)計抽樣次數對流動速度統(tǒng)計噪聲的影響呈平方根反比例關系，這一點可以由圖3的計算結果得到很好地說明。圖3（a）、（b）分別繪出固定板速Uw＝0．025 Ma和模擬分子總數Nt＝1000，取不同抽樣次數N＝20000、80000、320000、1280000得到的 DSMC 模擬結果及其全局最大相對統(tǒng)計漲落Emax和全局平均相對統(tǒng)計漲落Eave，顯示出DSMC方法對用于統(tǒng)計抽樣的次數具有較強的依賴性，如果N太少，如圖中N＝20000所得到的速度分布統(tǒng)計漲落就相當大，只有當抽樣次數達到N＝1280000這樣的量級，所得到的DSMC統(tǒng)計量就表現出較為光滑分布。圖3（b）繪出Emax、Eave的對數圖像分布，結果表明對此Uw＝0．025 Ma的Couette流動模擬，若將全場模擬分子總數固定在Nt＝1000的條件下，Emax和Eave在斜率為－1／2的擬合線附近散布，兩者之間依然差一個常數量。如希望把統(tǒng)計漲落控制在1%以下，采用Eave則抽樣次數達到640000就可滿足，而若采用控制Emax的辦法則需統(tǒng)計抽樣至少6×106次以上，這將以巨大的計算時間開支為代價。為了考察圖1～圖3描述的三種情況下DSMC模擬值的準確性，圖4繪出相應于（a）Uw＝0．2 Ma（設置 N 為20000及 Nt為1000）、（b）Nt＝64000（設置Uw為0．025 Ma及N 為20000）、（c）N＝1280000（設置Uw為0．025 Ma及Nt為1000）三種情況得到的Couette平板間速度分布與直接求解Boltzmann模型方程的統(tǒng)一算法［26］結果定量化比較，看出只要按照本文提出的相對統(tǒng)計漲落理論規(guī)則（14）、（15）式確定模擬參數，上述三種情況所得到的DSMC模擬值就彼此吻合且與文獻［26］統(tǒng)一算法計算分布在誤差范圍內相當一致。另一方面，由式（15）可計算出上述三種情況下Couette流動速度DSMC模擬值的最大相對統(tǒng)計漲落Emax在95%的置信水平下為大約2．3%，全局平均相對統(tǒng)計漲落Eave則低于1%，三種情況下得到的DSMC模擬值均趨于統(tǒng)一算法［26］結果。圖1～圖4證實本文關于DSMC模擬結果統(tǒng)計噪聲的理論推導公式（14）、（15）的可靠實用性。鑒于一般情況下計算流場的來流速度是給定的，為了減小流場速度統(tǒng)計噪聲的影響，應該選擇控制模擬分子總數Nt和抽樣次數N 的大小。根據圖2、圖3所反映的計算結果證實，控制這兩個模擬參數的取值，對于減小統(tǒng)計噪聲的效果是一樣的。在實際計算中，為了減小計算機內存負荷，一般可采取如圖3所示通過增大抽樣次數N，延長DSMC模擬計算時間，來獲得較小統(tǒng)計噪聲的DSMC模擬結果。通過上述過渡流區(qū)Couette剪切流動算例的計算比較與分析也看出，來流速度大小對DSMC統(tǒng)計噪聲的影響很大；來流速度越小，獲得小噪聲的流場速度分布越困難，可根據式（15）所反映的計算規(guī)則來選取DSMC模擬參數，并將統(tǒng)計噪聲控制在所希望的計算精度范圍內。

3．2 激波內流動計算研究

正激波內流動作為從均勻超聲速上游（狀態(tài)1）到均勻亞聲速下游（狀態(tài)2）流動的過渡，由于其一維屬性以及不考慮氣面邊界相互作用，常用作一種檢驗數值方法與計算規(guī)則的偏離熱力學平衡態(tài)流動算例。擬定來自文獻［20，27］相同的計算條件，使用變軟球（VSS）分子模型，對馬赫數Ms＝1．55的激波內流動問題進行DSMC模擬研究。為了檢驗2．3節(jié)推導出的流場溫度相對統(tǒng)計漲落表達式，使用式（18）作為控制DSMC方法模擬結果相對統(tǒng)計漲落計算規(guī)則。在處理溫度的DSMC模擬統(tǒng)計噪聲時，若取參考溫度Tref為流場溫度ˉTj，流場宏觀流動參數對溫度的相對統(tǒng)計漲落影響不明顯。根據3．1節(jié)模擬體會，由于模擬分子總數Nt和抽樣次數N 對DSMC模擬結果相對統(tǒng)計漲落的影響效果一樣，此處考慮相同的計算流場模擬分子總數Nt＝20841（計算過程中總的模擬分子數存在一定波動），選取不同的統(tǒng)計抽樣次數N＝500、1000、2000、4000、8000、16000、32000和64000，考察統(tǒng)計噪聲對DSMC模擬結果的影響。圖5繪出三種不同抽樣次數N＝500、2000和8000得到的溫度分布，看出對于這種較低馬赫數激波內流動模擬，統(tǒng)計抽樣次數對DSMC模擬結果噪聲影響很大；對于選取太少的N＝500，所得到的溫度分布統(tǒng)計漲落表現得很嚴重；當N值增大，激波輪廓的統(tǒng)計波動明顯減?。蝗舨捎幂^大的統(tǒng)計抽樣次數如N＝8000時，所得到的溫度分布輪廓線較光滑，見圖中虛線表示。圖6給出了分別在三種更高抽樣次數N＝8000、16000和32000下得到的溫度分布與來自文獻［20，27］的統(tǒng)一算法結果比較，圖中顯示出當抽樣次數大于8000時統(tǒng)計噪聲對模擬結果的影響已經不大。綜合圖5和圖6不同抽樣次數下激波結構內流動溫度分布DSMC模擬結果變化特點，可看出N值由小到大得到的DSMC模擬值均收斂并完全吻合于直接求解Boltzmann模型方程的統(tǒng)一算法結果，反映出統(tǒng)計量的相對統(tǒng)計漲落隨著抽樣次數增大而減弱，當統(tǒng)計抽樣次數增大到一定程度，DSMC模擬將收斂于共同的真值。為了將DSMC模擬值與實驗數據比較，圖7繪出相應于N＝32000模擬得到的密度分布與統(tǒng)一算法結果［20，27］及實驗數據［28］相比較，可看出三種結果變化趨勢相當一致。圖8繪出不同抽樣次數下DSMC模擬得到的激波內流動溫度分布的全局最大相對統(tǒng)計漲落Emax與全局平均相對統(tǒng)計漲落Eave隨抽樣次數N的變化關系對數圖像，可看出Emax、Eave分別散布在兩條平行直線附近，該直線分別由不同情況下Emax、Eave值擬合而來，其斜率均為－1／2，這直接證實關于溫度分布DSMC模擬結果相對統(tǒng)計漲落的表達式（18）的正確性，該式給出了DSMC模擬熱力學非平衡溫度場問題所需要統(tǒng)計抽樣次數的選取準則。統(tǒng)計漲落給出一個波動范圍，最終真實計算的漲落值都落在這個范圍之內。圖中數據分析表明，在抽樣次數N上升到104量級時，最大相對統(tǒng)計漲落可以控制在1%以內，而此時平均相對統(tǒng)計漲落已降至千分之一以內。該圖還反映了全局最大相對統(tǒng)計漲落Emax和全局平均相對統(tǒng)計漲落Eave具有完全相似的變化關系，Eave較之于Emax更為苛刻，一般計算只要將Eave控制在所要求計算精度就可以了，在實際應用中可根據客觀情況選擇判斷。

圖5 不同抽樣次數N＝500、2000、8000得到的Ms＝1．55激波內流動溫度分布Fig．5 Temperature distribution for Ms＝1．55shock wave under different N＝500，2000，8000

圖6 N＝8000、16000、32000得到的激波結構溫度分布DSMC模擬值與統(tǒng)一算法結果比較Fig．6 DSMC simulated temperature distribution for Ms＝1．55under different N＝8000，16000，32000 with the comparison of results from the unified algorithm

圖7 N＝32000得到的激波結構密度分布DSMC模擬值與統(tǒng)一算法、實驗數據比較Fig．7 DSMC simulated density distribution under N＝32000with comparison of the results from the unified algorithm and the experiment data

圖8 DSMC模擬激波內流動溫度分布的Emax、Eave隨統(tǒng)計抽樣次數N變化關系Fig．8 Relative statistical fluctuations Emaxand Eaveof temperature distribution under different Nfrom DSMC simulation

4 結束語

本文在定義DSMC模擬結果相對統(tǒng)計漲落的基礎上，引入概率統(tǒng)計置信度概念，通過數學分析推導出流場速度與溫度的相對統(tǒng)計漲落表達式。在此基礎上，以過渡流區(qū)Couette剪切流和激波內流動為研究對象，進行DSMC模擬研究，分別計算分析了Couette流動中速度的相對統(tǒng)計漲落、激波結構內流場溫度的相對統(tǒng)計漲落及其變化規(guī)律與影響因素，證實本文關于DSMC方法中統(tǒng)計噪聲的理論推導與控制規(guī)則的正確性。研究表明，在一定的置信度下，DSMC模擬結果的相對統(tǒng)計漲落除與網格內模擬分子數及統(tǒng)計抽樣次數的平方根成反比例，流場速度的相對統(tǒng)計漲落還與流場馬赫數成反比，而溫度的相對統(tǒng)計漲落與流場宏觀參數無關。得到了有關速度與溫度相對統(tǒng)計漲落更具一般性的結論，為DSMC方法模擬中控制并降低統(tǒng)計量的相對統(tǒng)計漲落提供了理論依據。而且計算體會到，為了控制流場參數的相對統(tǒng)計漲落，需要增大網格內模擬分子數或增大DSMC模擬過程中統(tǒng)計抽樣的次數，且二者是等效的，可根據實際計算需要，進行選取。一般為了控制計算內存，通常以增大DSMC統(tǒng)計抽樣次數，增加計算時間來降低統(tǒng)計量的相對統(tǒng)計漲落。

［1］ BIRD G A．Approach to translational equilibrium in a rigid sphere gas［J］．Phys．Fluids，1963，6（10）：1518－1519．

［2］ BORGNAKKE C，LARSEN P S．Statistical collision model for Monte Carlo simulation of polyatomic gas mixture［J］．J．of Comput．Phys．，1975，18（4）：405－420．

［3］ PHAM－Van DIEP G，ERWIN D，MUNTZ E P．Nonequilibrium molecular motion in a hypersonic shock wave［J］．Science，1989，245（4918）：624－626．

［4］ CARLSON A B，HASSAN H A．Direct simulation of reentry flows with ionization［R］．AIAA Paper 90－144，1990．

［5］ BIRD G A．Molecular gas dynamics and the direct simulation of gas flows［M］．Clarendon Press，Oxford，1994．

［6］ 樊菁，沈青．過渡領域高超聲速圓柱繞流的直接模擬［J］．空氣動力學學報，1995，13（4）：405－413．

［7］ 李志輝，吳振宇．阿波羅指令艙稀薄氣體動力學特征的蒙特卡羅數值模擬［J］．空氣動力學學報，1996，14（2）：230－233．

［8］ 吳其芬，陳偉芳．高溫稀薄氣體熱化學非平衡流動的DSMC方法［M］．長沙：國防科技大學出版社，1999．

［9］ IVANOV MS，VASHCHENKOV P，KASHKOVSKY A，Numerical investigation of the EXPERT reentry vehicle aerothermodynamics along the descent trajectory［R］．AIAA 2007－4145，2007．

［10］LI Z H，XIE Y R．Technique of molecular indexing applied to the direct simulation Monte Carlo method［A］．Proc．of 20th International Symposium on Rarefied Gas Dynamics［C］．ed．by Shen C，Peking Univ．Press，pp205－209，1997．

［11］LIANG J，LI Z H，LI Z H．DSMC computation of rarefied flow around a three－dimensional re－entry capsule［A］．3th Sino－Russian Hypersonic Flow Conference［C］．SRHFC－3，Wulumuqi，China，Sep．，2002．

［12］裴鹿成，張孝澤．蒙特卡羅方法及其在粒子輸運問題中的應用［M］．北京：科學出版社，1980．

［13］FRANCIS J．ALEXANDER．Cell size dependence of transport coefficients in stochastic particle algorithms［J］．Phys．Fluids，1998，10（6）：1540－1542．

［14］ALEJANDRO L．GARCIA．Time step truncation error in direct simulation Monte Carlo［J］．Phys．Fluids，2000，12（10）：2621－2633．

［15］RADER D J，GALLIS M A，TORCZYNSKI J R．Direct simulation Monte Carlo convergence behavior of the hard－sphere－gas thermal conductivity for Fourier heat flow［J］．Phys．Fluids，2000，18（7）：1－16．

［16］NICOLAS G．HADJICONSTANTINOU，ALEJANDRO L．GARCIA，MARTIN Z．BAZANT，GANG HE．Statistical error in particle simulations of hydrodynamic phenomena［J］．J．of Comput．Phys．，2003，187（1）：274－297．

［17］CHAPMAN S，COWLING T G．The mathematical theory of non－uniform gases［M］．3rd ed．，Cambridge Univ．Press，1990．

［18］YEN S M．Numerical solution of the nonlinear Boltzmann equation for nonequilibrium gas flow problems［J］．Ann．Rev．Fluid．Mech，1984，16（1）：67－97．

［19］YANG J Y，HUANG J C．Rarefied flow computations using nonlinear model boltzmann equations［J］．J．of Comput．Phys．，1995，120（2）：323－339．

［20］LI Z H，ZHANG H X．Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］．J．of Comput．Phys．，2004，193（2）：708－738．

［21］ALEXEI H，PIOTR K．Fast numerical method for the Boltzmann equation on non－uniform grids［J］．J．of Comput．Phys．，2008，227（13）：6681－6695．

［22］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，Jason Hall（Ed．），InTech Publisher，Chapter 14，2011：273－332．

［23］MORRIS A B，VARGHESE P L，GOLDSTEIN D B．Monte Carlo solution of the Boltzmann equation via a discrete velocity model［J］．J．of Comput．Phys．，2011，230（4）：1265－1280．

［24］WAGNER W．A convergence proof for Bird＇s direct simulation Monte Carlo method for the Boltzmann equation［J］．J．Sat．Phys．，1992，66（3／4）：1011－1044．

［25］汪仁官．概率論引論［M］．北京大學出版社，1994．

［26］李志輝．稀薄流到連續(xù)氣體流動問題統(tǒng)一算法應用研究［D］．［博士后研究報告］．北京：清華大學，2003．

［27］李志輝，張涵信．激波結構內流動問題的氣體運動論描述［J］．空氣動力學學報，2007，25（4）：411－418．

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