DSMC方法中的統(tǒng)計噪聲分析

李志輝，方 明，唐少強(qiáng)

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

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

0 引 言

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

1 DSMC方法

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

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

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

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

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

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

2 統(tǒng)計噪聲的數(shù)學(xué)分析

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

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

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

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

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

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

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

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

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

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

此處U為參考速度。若取參考速度U為該網(wǎng)格處氣體流動的當(dāng)?shù)厮俣萰，且鑒于第j個網(wǎng)格處的聲速ac滿足

即有

其中，Ma為第j個網(wǎng)格處氣體流動當(dāng)?shù)伛R赫數(shù)，γ為氣體的絕熱指數(shù)。

結(jié)合式（4）、（7）及（12），可得

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

該公式意味著，當(dāng)取參考速度U為網(wǎng)格j處當(dāng)?shù)貧饬魉俣取j時，速度的相對統(tǒng)計漲落與流動馬赫數(shù)成反比，與抽樣次數(shù)成平方根反比例關(guān)系。

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

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

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

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

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

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

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

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

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

3 算例與結(jié)果分析

3．1 Couette流動計算分析

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

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

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

圖3 不同抽樣次數(shù)得到的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)一算法結(jié)果比較Fig．4 DSMC simulated velocity distribution under different conditions with the comparison of results from the unified algorithm

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

3．2 激波內(nèi)流動計算研究

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

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

圖6 N＝8000、16000、32000得到的激波結(jié)構(gòu)溫度分布DSMC模擬值與統(tǒng)一算法結(jié)果比較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得到的激波結(jié)構(gòu)密度分布DSMC模擬值與統(tǒng)一算法、實驗數(shù)據(jù)比較Fig．7 DSMC simulated density distribution under N＝32000with comparison of the results from the unified algorithm and the experiment data

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

4 結(jié)束語

本文在定義DSMC模擬結(jié)果相對統(tǒng)計漲落的基礎(chǔ)上，引入概率統(tǒng)計置信度概念，通過數(shù)學(xué)分析推導(dǎo)出流場速度與溫度的相對統(tǒng)計漲落表達(dá)式。在此基礎(chǔ)上，以過渡流區(qū)Couette剪切流和激波內(nèi)流動為研究對象，進(jìn)行DSMC模擬研究，分別計算分析了Couette流動中速度的相對統(tǒng)計漲落、激波結(jié)構(gòu)內(nèi)流場溫度的相對統(tǒng)計漲落及其變化規(guī)律與影響因素，證實本文關(guān)于DSMC方法中統(tǒng)計噪聲的理論推導(dǎo)與控制規(guī)則的正確性。研究表明，在一定的置信度下，DSMC模擬結(jié)果的相對統(tǒng)計漲落除與網(wǎng)格內(nèi)模擬分子數(shù)及統(tǒng)計抽樣次數(shù)的平方根成反比例，流場速度的相對統(tǒng)計漲落還與流場馬赫數(shù)成反比，而溫度的相對統(tǒng)計漲落與流場宏觀參數(shù)無關(guān)。得到了有關(guān)速度與溫度相對統(tǒng)計漲落更具一般性的結(jié)論，為DSMC方法模擬中控制并降低統(tǒng)計量的相對統(tǒng)計漲落提供了理論依據(jù)。而且計算體會到，為了控制流場參數(shù)的相對統(tǒng)計漲落，需要增大網(wǎng)格內(nèi)模擬分子數(shù)或增大DSMC模擬過程中統(tǒng)計抽樣的次數(shù)，且二者是等效的，可根據(jù)實際計算需要，進(jìn)行選取。一般為了控制計算內(nèi)存，通常以增大DSMC統(tǒng)計抽樣次數(shù)，增加計算時間來降低統(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］ 樊菁，沈青．過渡領(lǐng)域高超聲速圓柱繞流的直接模擬［J］．空氣動力學(xué)學(xué)報，1995，13（4）：405－413．

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

［8］ 吳其芬，陳偉芳．高溫稀薄氣體熱化學(xué)非平衡流動的DSMC方法［M］．長沙：國防科技大學(xué)出版社，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］裴鹿成，張孝澤．蒙特卡羅方法及其在粒子輸運問題中的應(yīng)用［M］．北京：科學(xué)出版社，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］．北京大學(xué)出版社，1994．

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

［27］李志輝，張涵信．激波結(jié)構(gòu)內(nèi)流動問題的氣體運動論描述［J］．空氣動力學(xué)學(xué)報，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．