高超聲速稀薄流場(chǎng)AUSM 分裂式DSMC-IP 方法研究

許 嘯，馬新建，張 軍，沈 妍

(1. 江蘇科技大學(xué)張家港校區(qū)機(jī)電與動(dòng)力工程學(xué)院，鎮(zhèn)江 212003；2. 上海航天技術(shù)研究院航天動(dòng)力技術(shù)研究院，上海 201109；3. 南京航空航天大學(xué)航空學(xué)院，南京 210001)

高超聲速飛行器是當(dāng)前航空航天領(lǐng)域的研究熱點(diǎn)，該類(lèi)飛行器通常以5Ma～20Ma的速度工作于臨近空間(海拔20 km～100 km)中，執(zhí)行高速巡航、全球快速打擊、空間運(yùn)輸?shù)戎匾蝿?wù)，具有極高的戰(zhàn)略?xún)r(jià)值，其高速運(yùn)動(dòng)狀態(tài)下的稀薄流場(chǎng)結(jié)構(gòu)與特性則是空氣動(dòng)力學(xué)領(lǐng)域的研究重點(diǎn)與難點(diǎn)[1?8]。在臨近空間中，大氣密度隨海拔高度的升高而急劇下降，同時(shí)飛行器表面邊界層的剪切作用會(huì)隨著速度的升高而增加，此時(shí)高超聲速飛行器流場(chǎng)中氣體粒子間斷效應(yīng)將越發(fā)明顯，從而導(dǎo)致流場(chǎng)整體的稀薄特性明顯增強(qiáng)[9]，基于連續(xù)性假設(shè)的Navier-Stokes(N-S)方程逐漸失效，以DSMC(Direct Simulation of Monte Carlo)方法為代表的稀薄氣體動(dòng)力學(xué)方法成為了研究高超稀薄流動(dòng)的主要手段。

DSMC 方法最初由Bird[10]提出，該方法采用大量模擬分子代替真實(shí)分子，將分子的運(yùn)動(dòng)和碰撞進(jìn)行解耦以模擬流動(dòng)變化過(guò)程，并通過(guò)統(tǒng)計(jì)得出流場(chǎng)宏觀信息。該方法基于分子動(dòng)力學(xué)理論，不僅易于實(shí)現(xiàn)真實(shí)的物理運(yùn)動(dòng)模型，而且已證明其收斂于Boltzmann 方程[11]，被認(rèn)為是目前最成功的稀薄流動(dòng)模擬方法之一。在高超聲速稀薄流動(dòng)領(lǐng)域，Reinert 等[12]采用DSMC 方法模擬了三維30°～55°雙楔體的高超聲速流動(dòng)，并與相關(guān)實(shí)驗(yàn)結(jié)果進(jìn)行了比較；Torres 等[13]對(duì)DSMC 方法進(jìn)行了拓展，提出了DMS(Direct Molecular Simulation)方法，并采用該方法模擬了高超稀薄沖擊波加熱氣體混合物的演化過(guò)程，驗(yàn)證了方法的有效性；Jin 等[14]采用DSMC 方法對(duì)長(zhǎng)深比為4 的三維空腔高超聲速稀薄繞流進(jìn)行了數(shù)值模擬，結(jié)果突出反映了空腔內(nèi)部流動(dòng)特性及其表面氣動(dòng)參數(shù)對(duì)Maxwell 調(diào)節(jié)系數(shù)和自由流Knudsen 數(shù)的敏感性；Wang 等[15]采用DSMC 方法模擬了海拔70 km 狀態(tài)下高超稀薄氣體通過(guò)障礙物時(shí)的繞流流場(chǎng)，并分析了不同障礙物結(jié)構(gòu)參數(shù)對(duì)表面氣動(dòng)量變化的影響；李潔等[16]基于DSMC 方法構(gòu)造了相應(yīng)的碰撞、聚合與分離模型，實(shí)現(xiàn)了高超聲速稀薄多相流噴流流場(chǎng)模擬；黃飛等[17]采用DSMC 方法模擬了高超聲速巡航飛行器的尖化前緣流場(chǎng)，并將結(jié)果與N-S 方法的計(jì)算結(jié)果相比較，研究了局部稀薄氣體效應(yīng)對(duì)氣動(dòng)熱特性分布的影響規(guī)律；屈程等[18]將DSMC 方法與結(jié)構(gòu)傳熱法相結(jié)合，對(duì)X-37B軌道飛行器外形長(zhǎng)時(shí)加熱與結(jié)構(gòu)傳熱進(jìn)行了模擬分析；張雋研等[19]采用DSMC 方法對(duì)臨近空間高超聲速飛行器的表面空腔進(jìn)行了模擬等等。這些成果均證明了DSMC 方法在高超聲速稀薄流動(dòng)模擬中的有效性和可靠性。

雖然DSMC 方法獲得了巨大的成功，但該方法采用隨機(jī)模型處理分子生成、碰撞和統(tǒng)計(jì)過(guò)程，導(dǎo)致該方法的計(jì)算結(jié)果不可避免的受到統(tǒng)計(jì)漲落的影響，這不但迫使DSMC 方法必須采用大量增加統(tǒng)計(jì)次數(shù)的方式來(lái)抹平統(tǒng)計(jì)誤差，而且使該方法在模擬流場(chǎng)非定常變化時(shí)始終無(wú)法得到統(tǒng)計(jì)精度較高的結(jié)果。為解決這一問(wèn)題，研究者們提出了兩種不同的思路：樊菁等[20]提出了DSMC-IP(Information Preservation)方法，其核心思想是以宏觀信息量替代隨機(jī)量進(jìn)行統(tǒng)計(jì)，從而降低計(jì)算中的統(tǒng)計(jì)耗散，最初IP 方法僅用于低速平衡流動(dòng)，后經(jīng)多位學(xué)者的研究與發(fā)展，該方法逐漸開(kāi)始適用于多種高速非平衡流動(dòng)的模擬[21?26]；Burt 等[27]提出了DSMC-LD(Low Diffusion)方法，與IP 方法類(lèi)似，該方法同樣將總體速度(bulk velocity)的概念引入分子運(yùn)動(dòng)模型，但不同的是，該方法將整個(gè)流場(chǎng)離散為一系列質(zhì)量恒定而結(jié)構(gòu)可變的封閉式拉格朗日單元，通過(guò)分子與單元界面之間的總碰撞量來(lái)傳遞單元間的動(dòng)量與能量變化，經(jīng)過(guò)多位研究者的開(kāi)發(fā)和優(yōu)化，該方法已成為近平衡流動(dòng)模擬中的重要研究手段[28?31]。

雖然上述兩種方法均取得了一定的成果，但相比于LD 方法，IP 方法無(wú)需處理繁雜的網(wǎng)格單元形變，算法也更易于實(shí)現(xiàn)，同時(shí)文獻(xiàn)[23 ? 26]也證明了IP 方法適用于強(qiáng)非平衡流動(dòng)模擬。據(jù)此，本文在過(guò)往工作[32? 33]的基礎(chǔ)上，對(duì)信息量控制方程中的關(guān)聯(lián)項(xiàng)通量分裂格式進(jìn)行深度優(yōu)化，構(gòu)建AUSM 分裂式DSMC-IP 方法，為高超聲速稀薄流場(chǎng)研究提供新的方案。

1 計(jì)算方法

1.1 DSMC-IP 方法概述

DSMC-IP 方法最初的設(shè)計(jì)目的在于解決DSMC方法模擬微尺度低速流動(dòng)時(shí)，分子熱運(yùn)動(dòng)速度遠(yuǎn)大于宏觀速度，導(dǎo)致統(tǒng)計(jì)噪聲過(guò)大、計(jì)算難以收斂穩(wěn)定的問(wèn)題，主要應(yīng)用于MEMS(Micro Electro Mechanical System)領(lǐng)域的研究。在設(shè)計(jì)之初，該方法基于DSMC 方法中模擬分子代表了大量真實(shí)分子組成的分子團(tuán)這一設(shè)定，將DSMC 方法中模擬分子的實(shí)際運(yùn)動(dòng)速度表述為：

式中：U為模擬分子的實(shí)際速度；V為真實(shí)分子團(tuán)的平均速度；C為分子的熱運(yùn)動(dòng)速度，當(dāng)氣體處于平衡態(tài)時(shí)，其內(nèi)部的分子熱運(yùn)動(dòng)速度之和應(yīng)當(dāng)為0，因此在宏觀狀態(tài)下C=0(C為分子的平均熱運(yùn)動(dòng)速度)，此時(shí)分子團(tuán)實(shí)際速度的平均值即宏觀速度。據(jù)此，DSMC-IP 方法將分子速度分為實(shí)際速度和平均速度兩部分，均由模擬分子攜帶，實(shí)際速度用于執(zhí)行DSMC 方法中模擬分子的位移與碰撞，平均速度則作為宏觀信息保存在每個(gè)模擬分子中，當(dāng)模擬分子運(yùn)動(dòng)時(shí)，采用獨(dú)立的方法對(duì)宏觀速度進(jìn)行更新，最終計(jì)算流場(chǎng)宏觀值時(shí)則采用分子所攜帶的宏觀信息速度進(jìn)行統(tǒng)計(jì)。由于宏觀速度值排除了分子熱運(yùn)動(dòng)的干擾，IP 方法可有效降低分子隨機(jī)信息帶來(lái)的統(tǒng)計(jì)耗散[20]。

雖然DSMC-IP 方法最初保存的信息量?jī)H有速度，而且其信息量演化的理論模型也不明確，但該方法“采用平均量(即信息量)代替隨機(jī)量進(jìn)行演化并統(tǒng)計(jì)從而獲得低統(tǒng)計(jì)耗散解”的思路獲得了成功，并吸引了眾多研究者的關(guān)注。在IP 方法提出之后，Sun 等[24]根據(jù)信息量的特點(diǎn)，將其定義為真實(shí)分子團(tuán)的物理特性的一階原點(diǎn)矩，并由此得出信息量控制方程組的守恒形式：

在DSMC-IP 方法中，分子信息量在運(yùn)動(dòng)過(guò)程中隨分子正常遷移而發(fā)生位置變化；當(dāng)進(jìn)行碰撞時(shí)需采用碰撞模型處理信息量的變化過(guò)程，Sun 等在文獻(xiàn)[22]中提出了一種唯象論模型，為信息量在碰撞過(guò)程中發(fā)生的變化提供了計(jì)算方法。將運(yùn)動(dòng)與碰撞項(xiàng)剝離后，信息量控制方程僅剩右端關(guān)聯(lián)項(xiàng)，此時(shí)，式(6)和式(7)可寫(xiě)為如下形式：

不難發(fā)現(xiàn)，簡(jiǎn)化之后的控制方程組已基本可以通過(guò)離散方法轉(zhuǎn)化為差分方程來(lái)求解，而此時(shí)求解的關(guān)鍵即在于如何處理右端微觀關(guān)聯(lián)項(xiàng)。為此，多位學(xué)者提出了適用于平衡假設(shè)的局部熱力平衡法(local thermodynamic equilibrium, LTE)[24]、適用于非平衡稀薄流動(dòng)模擬的通量分裂法[24]和能夠更加精準(zhǔn)計(jì)算分子間剪切動(dòng)量的八分通量分裂法(octant flux splitting, OFS)[25]。從方法特點(diǎn)上看，LTE 方法將分子分布函數(shù)建立在平衡態(tài)假設(shè)上，而高超稀薄流動(dòng)中存在著大量的非平衡區(qū)域；OFS 方法精度高但計(jì)算量大，適用于研究純熱壓驅(qū)動(dòng)流中的剪切作用。據(jù)此，本文針對(duì)高超聲速流場(chǎng)以激波壓縮效應(yīng)為主的特點(diǎn)，采用FS 方法為基礎(chǔ)并進(jìn)行優(yōu)化，在保證結(jié)果準(zhǔn)確性的條件下提高計(jì)算效率。

1.2 AUSM 分裂式通量計(jì)算方法

根據(jù)文獻(xiàn)[24]中的研究，信息量控制方程可進(jìn)行進(jìn)一步簡(jiǎn)化，其形式為：

將式(10)～式(12)轉(zhuǎn)換為積分形式的離散方程，并采用差分方法即可實(shí)現(xiàn)近似求解，其過(guò)程類(lèi)似于采用基于N-S 方程的有限體積法構(gòu)造通量計(jì)算連續(xù)性流場(chǎng)。IP 方法最初用于求解低速流動(dòng)，因此可采用較為簡(jiǎn)單的中心差分格式構(gòu)造空間離散通量。然而，在高超稀薄流場(chǎng)中，雖然FS 方法能夠較好的體現(xiàn)出非平衡效應(yīng)，但中心差分格式無(wú)法捕捉流動(dòng)中的強(qiáng)激波，必須采用迎風(fēng)格式對(duì)關(guān)聯(lián)項(xiàng)離散通量進(jìn)行重構(gòu)才能獲得較好的計(jì)算結(jié)果，尤其是結(jié)構(gòu)較為復(fù)雜的流場(chǎng)；但困難的是，DSMC-IP 方法的控制方程不同于N-S 方程，即使以宏觀密度ρ 替換方程中的nm，動(dòng)量和能量方程右端關(guān)聯(lián)項(xiàng)中各項(xiàng)變量的宏觀意義仍不十分明確，無(wú)法完全類(lèi)比于N-S 方程中的對(duì)流項(xiàng)或粘性項(xiàng)，因此如何準(zhǔn)確的構(gòu)造關(guān)聯(lián)項(xiàng)離散通量來(lái)滿(mǎn)足DSMC-IP 方法在模擬高超聲速稀薄流動(dòng)時(shí)的計(jì)算要求是該方法面臨的重要問(wèn)題。本文作者曾嘗試采用Van Leer 格式對(duì)IP 關(guān)聯(lián)項(xiàng)控制方程進(jìn)行重構(gòu)[33]，雖然方法初步獲得了成功，但該格式對(duì)動(dòng)量和能量通量的重構(gòu)式與控制方程中對(duì)應(yīng)的關(guān)聯(lián)項(xiàng)契合度不高，因此僅將其用于質(zhì)量通量的重構(gòu)，動(dòng)量和能量在激波處的突變則只能通過(guò)模擬分子在激波處碰撞時(shí)交換其攜帶的信息量來(lái)反映，從而限制了IP 方法對(duì)高超稀薄流動(dòng)的模擬能力。

根據(jù)文獻(xiàn)[34]，AUSM 格式將流動(dòng)中的對(duì)流項(xiàng)拆分為對(duì)流部分和壓力部分，計(jì)算中單元分界面上的通量具體形式如式(13)所示。式中：I為單元編號(hào)；L、R為左右單元標(biāo)識(shí)；Man為局部馬赫數(shù)；ρ 為密度；c為音速；p為壓力；H為總焓；u、v、w分別代表三維方向上的速度；nx、ny、nz分別為單元分界面在三維方向上的單位外法向量。如果將該通量公式與DSMC-IP 方法的關(guān)聯(lián)項(xiàng)控制方程(式(10)～式(12))進(jìn)行對(duì)比則不難發(fā)現(xiàn)：質(zhì)量關(guān)聯(lián)項(xiàng)中的，

nmVc與AUSM 格式中的ρcMan相似；動(dòng)量關(guān)聯(lián)項(xiàng)中的nmRT與AUSM 格式中的P相似；在能量關(guān)聯(lián)項(xiàng)中，雖然nmRTVk與AUSM 格式中的ρHcMan存在一定差異，但總焓H的定義式為：

式中：E為總能，即流體微元的內(nèi)能與動(dòng)能之和，而根據(jù)理想氣體方程，壓力項(xiàng)p/ρ 等于RT，因此如按式(14)將AUSM 格式中的能量對(duì)流通量進(jìn)行拆解，可得其子項(xiàng)ρRTcMan，其與關(guān)聯(lián)項(xiàng)中的nmRTVk是相似的。

綜上所述，DSMC-IP 方法關(guān)聯(lián)項(xiàng)中的各項(xiàng)變量在AUSM 分裂式通量中均有相似的表述。從物理意義上講，IP 方法中的運(yùn)動(dòng)項(xiàng)和關(guān)聯(lián)項(xiàng)共同描述了對(duì)流輸運(yùn)對(duì)宏觀信息量的影響，其中運(yùn)動(dòng)項(xiàng)主要對(duì)應(yīng)于AUSM 通量中的對(duì)流部分，如ρucMan和ρEcMan等，表示流體運(yùn)動(dòng)對(duì)物理量的改變，在IP 方法中表現(xiàn)為信息量在分子遷移中的變化；而關(guān)聯(lián)項(xiàng)則對(duì)應(yīng)于AUSM 通量中的壓力部分，反映了宏觀壓力對(duì)局部流場(chǎng)的影響，在IP 方法中表現(xiàn)為信息量受相鄰單元“類(lèi)壓力”和“類(lèi)壓力熱流”影響而發(fā)生的改變。據(jù)此，本文采用AUSM格式對(duì)IP 方法中的各類(lèi)關(guān)聯(lián)項(xiàng)通量計(jì)算格式進(jìn)行重構(gòu)，使其能夠更好的適應(yīng)高超稀薄流場(chǎng)的模擬要求。

AUSM 分裂式DSMC-IP 方法的關(guān)聯(lián)項(xiàng)通量計(jì)算步驟如下：

4)信息量更新：采用時(shí)間顯式推進(jìn)格式計(jì)算其下一時(shí)間步的信息量：

以上為AUSM 分裂式IP 方法中關(guān)聯(lián)項(xiàng)通量計(jì)算的主要步驟，式(21)中：Q為不同類(lèi)型的信息量；t為時(shí)間；i為各邊界面的編號(hào)；N為單元中的邊界總數(shù)；Δt為時(shí)間步長(zhǎng)度；FQi為單元中各邊界上的總通量；Ai為各邊界的面積(二維情況下為長(zhǎng)度)；Ω為單元體積(二維情況下為面積)，據(jù)此可得關(guān)聯(lián)項(xiàng)更新后的信息量。將該信息量保存于分子中隨其運(yùn)動(dòng)，并在碰撞過(guò)程中采用唯象論碰撞模型進(jìn)行處理，即可得下一時(shí)間步的分子信息量。

2 數(shù)值算例

2.1 超聲速圓柱繞流數(shù)值模擬

本文首先采用AUSM 分裂式DSMC-IP 方法對(duì)超聲速圓柱繞流流場(chǎng)進(jìn)行數(shù)值模擬。該算例來(lái)自文獻(xiàn)[35]，圓柱半徑為0.05 m，其表面溫度為500 K，來(lái)流氣體為氮?dú)?，?lái)流速度1412.5 m/s，溫度300 K，分子轉(zhuǎn)動(dòng)自由度為2，粘性系數(shù)的溫度指數(shù)為0.74，來(lái)流氣體分子數(shù)密度為1.0×1020，其分子自由程約0.0133 m，Kn數(shù)為0.13，屬于典型的稀薄流動(dòng)。計(jì)算網(wǎng)格單元總數(shù)為3270，模擬分子總數(shù)約116 萬(wàn)，網(wǎng)格與模型如圖1 所示。本文計(jì)算時(shí)間步取1.0×10?8s，每100 時(shí)間步統(tǒng)計(jì)一次，計(jì)算至流場(chǎng)基本穩(wěn)定時(shí)進(jìn)行統(tǒng)計(jì)得出收斂結(jié)果。

圖1 超聲速圓柱繞流模型與網(wǎng)格Fig. 1 Model and grid of supersonic flow around the circular cylinder

圖2～圖4 展示了AUSM 分裂式DSMC-IP方法與普通DSMC 方法不同時(shí)刻的統(tǒng)計(jì)結(jié)果，并將計(jì)算結(jié)果進(jìn)行了對(duì)比。兩種方法的結(jié)果都完整的展示了超聲速繞流流場(chǎng)的發(fā)展過(guò)程：高速來(lái)流氣體首先在圓柱迎風(fēng)表面受壓縮形成激波區(qū)域，隨后繞經(jīng)圓柱表面在其后方形成膨脹尾流，激波區(qū)與膨脹區(qū)在來(lái)流的推動(dòng)下逐漸發(fā)展，最終相互融合，形成穩(wěn)定的繞流流場(chǎng)。通過(guò)對(duì)比不難發(fā)現(xiàn)，IP 方法計(jì)算結(jié)果的實(shí)時(shí)統(tǒng)計(jì)精度明顯優(yōu)于DSMC 方法，即使在未經(jīng)大量統(tǒng)計(jì)的實(shí)時(shí)統(tǒng)計(jì)結(jié)果中，IP 方法結(jié)果中的等值線(xiàn)和區(qū)域均清晰而光滑；DSMC 方法的結(jié)果則呈現(xiàn)出大量的點(diǎn)狀和塊狀結(jié)構(gòu)，其等值線(xiàn)也呈現(xiàn)出較大的擾動(dòng)，即使經(jīng)過(guò)一定時(shí)間的統(tǒng)計(jì)，激波邊緣的溫度等值線(xiàn)仍帶有統(tǒng)計(jì)漲落的干擾痕跡，據(jù)此可以證明本文提出的IP 方法可有效降低超聲速流場(chǎng)中的統(tǒng)計(jì)耗散，提高實(shí)時(shí)計(jì)算結(jié)果的精度。

圖2 t=1.25×10?4 s 時(shí)流場(chǎng)中宏觀量計(jì)算結(jié)果Fig. 2 The computational results of macro quantities in flow field when t=1.25×10?4 s

圖3 t=2.5×10?4 s 時(shí)流場(chǎng)中宏觀量計(jì)算結(jié)果Fig. 3 The computational results of macro quantities in flow field when t=2.5×10?4 s

進(jìn)一步對(duì)比兩種方法的等值線(xiàn)結(jié)果，可以發(fā)現(xiàn)，流場(chǎng)發(fā)展初始時(shí)兩者的結(jié)構(gòu)和位置基本一致，但流動(dòng)穩(wěn)定時(shí)激波前緣和尾流區(qū)的部分等值線(xiàn)位置稍有差異，尤其是溫度等值線(xiàn)。這主要是由IP 和DSMC 方法的格式差異造成的：FS 式IP 方法按分子速度大小將單元內(nèi)的分子假定為兩種平衡態(tài)分子群的集合，這一假設(shè)符合了非平衡狀態(tài)下分子速度的分布情況，但在處理分子群內(nèi)的信息溫度時(shí)仍采用了平衡假設(shè)，當(dāng)流場(chǎng)壓縮程度較低且宏觀物理量梯度較大時(shí)，其非平衡特征較為明顯，此時(shí)信息溫度的計(jì)算精度會(huì)低于標(biāo)準(zhǔn)的DSMC 方法。在激波前緣，流場(chǎng)中物理量梯度激增；而在尾流區(qū)中，由于背風(fēng)效應(yīng)，流場(chǎng)密度隨流動(dòng)的發(fā)展不斷下降，這兩處區(qū)域的非平衡特性明顯增強(qiáng)，IP 方法的計(jì)算結(jié)果也因此出現(xiàn)了略低于標(biāo)準(zhǔn)的情況。但在激波中心和近壁面等在飛行器設(shè)計(jì)中至關(guān)重要的區(qū)域里，IP 方法與標(biāo)準(zhǔn)DSMC 方法的計(jì)算結(jié)果符合較好。由此可以看出，AUSM 式DSMC-IP 方法的實(shí)時(shí)計(jì)算精度優(yōu)于DSMC 方法，流場(chǎng)收斂后的計(jì)算結(jié)果與DSMC方法的標(biāo)準(zhǔn)結(jié)果基本符合，僅在流場(chǎng)非平衡效應(yīng)顯著的局部區(qū)域內(nèi)略偏低。從結(jié)果來(lái)看，等值線(xiàn)最大偏差約占流場(chǎng)在X方向上總長(zhǎng)度的3.2%，并不十分明顯，據(jù)此可以初步證明AUSM式DSMCIP 方法的準(zhǔn)確性和有效性。

圖5～圖7 比較了AUSM 分裂式DSMC-IP 方法與DSMC 方法不同時(shí)刻計(jì)算結(jié)果中駐點(diǎn)線(xiàn)上的馬赫數(shù)與壓力值?？梢钥闯?，兩種方法計(jì)算結(jié)果的趨勢(shì)和變化范圍基本相同，兩者之間的差異出現(xiàn)在激波邊緣處，并在激波內(nèi)測(cè)沿來(lái)流方向逐漸減小，在近壁面處接近相等。其原因上一段已經(jīng)分析過(guò)，主要是兩種方法的計(jì)算格式有別，導(dǎo)致結(jié)果受局部流場(chǎng)非平衡效應(yīng)變化的影響。從數(shù)值上看，以統(tǒng)計(jì)結(jié)果為標(biāo)準(zhǔn)，馬赫數(shù)和壓力的最大差距分別約占總數(shù)值變化范圍的7.5%和8%，等值點(diǎn)之間的位置差距最高約占總激波寬度的8%，均在10%以?xún)?nèi)，說(shuō)明兩種方法自身的差異對(duì)結(jié)果造成的影響是比較有限的。另外，在DSMC 方法的實(shí)時(shí)計(jì)算結(jié)果中，激波外層的數(shù)值受統(tǒng)計(jì)漲落的干擾，波動(dòng)性非常明顯，而本文提出的IP 方法則明顯具有更好的穩(wěn)定性，進(jìn)一步證明了該方法的有效性。

圖5 t=1.25×10?4 s 時(shí)駐點(diǎn)線(xiàn)上的馬赫數(shù)和壓力分布Fig. 5 The Mach numbers and pressure distribution on the stagnation line when t=1.25×10?4 s

圖6 t=2.5×10?4 s 時(shí)駐點(diǎn)線(xiàn)上的馬赫數(shù)和壓力分布Fig. 6 The Mach numbers and pressure distribution on the stagnation line when t=2.5×10?4 s

圖7 流動(dòng)穩(wěn)定時(shí)駐點(diǎn)線(xiàn)上的馬赫數(shù)和壓力分布Fig. 7 The Mach numbers and pressure distribution on the stagnation line when the flow is stable

圖8 給出了兩種方法收斂結(jié)果中固壁表面的壓力、熱流和摩阻系數(shù)(Cp、Ch和Cf，L為弧長(zhǎng))，并將結(jié)果與文獻(xiàn)[35]中的對(duì)應(yīng)結(jié)果進(jìn)行了對(duì)比。圖中結(jié)果顯示，壓力與熱流系數(shù)的峰值均出現(xiàn)在圓柱迎風(fēng)面的駐點(diǎn)位置，摩阻峰值則出現(xiàn)在駐點(diǎn)后沿壁面流動(dòng)方向弧長(zhǎng)約0.04 m 的位置，其結(jié)果表明高速來(lái)流在圓柱表面駐點(diǎn)處將動(dòng)能最大程度的轉(zhuǎn)化為流體內(nèi)能，隨后在迎風(fēng)面形成繞流時(shí)切向阻力受圓弧曲率的影響不斷增大，至背風(fēng)區(qū)后在膨脹作用下又不斷減小的特點(diǎn)，該過(guò)程與圓柱繞流的發(fā)展相符。通過(guò)數(shù)值比較可以看出，本文提出的IP 方法與DSMC 以及參考文獻(xiàn)[35]中的數(shù)值吻合較好，各系數(shù)中最大差距值均不超過(guò)總數(shù)值變化幅度的5%，由此可以證明本文方法的可靠性。

圖8 計(jì)算穩(wěn)定時(shí)圓柱表面特征系數(shù)分布Fig. 8 The characteristic coefficients distribution on the cylinder surface when the flow is stable

表1 列出了兩種方法的計(jì)算時(shí)間比較。兩種方法從流動(dòng)初始到流場(chǎng)穩(wěn)定共計(jì)算35000 時(shí)間步，但由于DSMC 方法統(tǒng)計(jì)耗散較大，最終收斂時(shí)共統(tǒng)計(jì)10000 時(shí)間步，而IP 方法則相反，只需統(tǒng)計(jì)3000 時(shí)間步即得到良好的收斂結(jié)果。與DSMC 方法相比，AUSM 分裂式IP 方法增加了分子宏觀信息量的計(jì)算和統(tǒng)計(jì)步驟，因此其算法代碼更多，消耗的內(nèi)存變量更大，計(jì)算時(shí)間自然也更長(zhǎng)。從表中的數(shù)據(jù)來(lái)看，從初始到穩(wěn)態(tài)時(shí)IP 方法的計(jì)算時(shí)間比DSMC 方法多約60%，單位時(shí)間步平均耗時(shí)約為DSMC 方法的1.53 倍；但由于統(tǒng)計(jì)時(shí)間步更少，IP 方法的統(tǒng)計(jì)時(shí)間更低，約為DSMC 方法的41.8%。

表1 圓柱繞流中DSMC 與IP 方法的計(jì)算時(shí)間比較Table 1 Comparison of calculation time between DSMC and IP method in circular cylinder flow

2.2 帶擴(kuò)張角的噴管外高超聲速流動(dòng)數(shù)值模擬

為進(jìn)一步驗(yàn)證AUSM 分裂式DSMC-IP 方法在高超聲速流場(chǎng)模擬中的適用性，本文采用該方法對(duì)帶擴(kuò)張角的噴管外高超聲速流動(dòng)進(jìn)行了模擬。該算例來(lái)自文獻(xiàn)[36]，噴管結(jié)構(gòu)如圖9 所示，模型尺寸單位為mm，其表面溫度為293 K，來(lái)流氣體為氮?dú)?，?lái)流速度1430 m/s，溫度51.9 K，來(lái)流馬赫數(shù)為9.9，來(lái)流氣體密度為4.226×10?4kg/m3，來(lái)流Kn數(shù)為0.0065。本算例中來(lái)流氣體的分子自由程較小(約1.5×10?4m)，來(lái)流馬赫數(shù)極高，且噴管外部構(gòu)型較為復(fù)雜，因此近壁面流場(chǎng)變化比較劇烈，為清晰準(zhǔn)確的描述流動(dòng)特性，本文使用單元數(shù)共計(jì)約4.2 萬(wàn)(如圖10 所示)，分子數(shù)約150 萬(wàn)，計(jì)算時(shí)間步取2.5×10?8s，每100 步統(tǒng)計(jì)一次，計(jì)算至流場(chǎng)基本穩(wěn)定時(shí)進(jìn)行統(tǒng)計(jì)，共統(tǒng)計(jì)5000 步得出最終結(jié)果。

圖9 高超聲速?lài)姽苣Ｐ徒Y(jié)構(gòu) /mmFig. 9 Structure of the hypersonic nozzle

圖10 帶擴(kuò)張角的噴管外高超聲速流場(chǎng)計(jì)算網(wǎng)格Fig. 10 The computational grid of the hypersonic nozzle

圖11～圖13 展示并對(duì)比了AUSM 分裂式DSMCIP 方法與DSMC 方法不同時(shí)刻的計(jì)算結(jié)果?？梢钥闯?，低溫高超來(lái)流首先在噴管外表面形成了繞流區(qū)，受壁面相對(duì)高溫的影響，初始時(shí)該區(qū)域中僅噴管內(nèi)側(cè)的擴(kuò)張段溫度相對(duì)較低；經(jīng)過(guò)發(fā)展后，噴管外側(cè)的來(lái)流在平直段表面形成了類(lèi)似于平板繞流的一次激波，隨后在外側(cè)擴(kuò)張段受其結(jié)構(gòu)影響形成二次斜激波，兩股激波交匯，在外側(cè)擴(kuò)張段末端形成了相互疊加的復(fù)雜激波流動(dòng)；噴管內(nèi)側(cè)來(lái)流繞經(jīng)前端傾斜角后形成小角度斜激波，在對(duì)稱(chēng)面處發(fā)生反射，對(duì)內(nèi)側(cè)擴(kuò)張段表面流場(chǎng)產(chǎn)生壓縮，并在尾流區(qū)與噴管外側(cè)激波后的膨脹流相互作用，形成溫度較高的壓縮流區(qū)域；流動(dòng)最終穩(wěn)定時(shí)，高壓區(qū)集中在管外側(cè)擴(kuò)張段表面與管內(nèi)側(cè)擴(kuò)張段至尾流區(qū)，膨脹流出現(xiàn)在噴管外側(cè)激波后方與噴管內(nèi)側(cè)對(duì)稱(chēng)面上接近下游出口的位置。整體來(lái)看，受?chē)姽軓?fù)雜結(jié)構(gòu)的影響，其內(nèi)外側(cè)流場(chǎng)經(jīng)歷了激波/激波干擾、激波/邊界層干擾、對(duì)稱(chēng)面激波反射等一系列復(fù)雜的高超聲速流動(dòng)過(guò)程；此外，圖14 顯示了IP 方法最終統(tǒng)計(jì)結(jié)果中噴管外側(cè)壁面轉(zhuǎn)角處與噴管尾部的渦旋流線(xiàn)結(jié)構(gòu)，可以看出該流場(chǎng)中還包含了流動(dòng)分離與環(huán)流等現(xiàn)象，而這些復(fù)雜流動(dòng)特征在AUSM 分裂式DSMC-IP 方法的計(jì)算結(jié)果中均有體現(xiàn)，并且與文獻(xiàn)中結(jié)果基本相符，如其中提到噴管外側(cè)的分離泡(separation bubble)，具體可參考文獻(xiàn)[36]。

圖12 t=1.25×10?4 s 時(shí)流場(chǎng)中宏觀量計(jì)算結(jié)果Fig. 12 The computational results of macro quantities in flow field when t=1.25×10?4 s

圖13 流動(dòng)穩(wěn)定后流場(chǎng)中宏觀量統(tǒng)計(jì)結(jié)果Fig. 13 The computational results of macro quantities in flow field when the flow is stable

圖14 流動(dòng)穩(wěn)定時(shí)DSMC-IP 方法計(jì)算結(jié)果中流場(chǎng)不同位置處的渦旋Fig. 14 Vortices at different positions in the flow field computed by the DSMC-IP method when the flow is stable

通過(guò)對(duì)比兩種方法的計(jì)算結(jié)果可以看出，噴管前端區(qū)域的激波等值線(xiàn)位置與數(shù)值幾乎完全吻合，而流動(dòng)膨脹區(qū)的等值線(xiàn)則略有差異，如初始時(shí)噴管壁面尾部、穩(wěn)定時(shí)噴管內(nèi)側(cè)對(duì)稱(chēng)面上靠近流動(dòng)出口的區(qū)域等。該結(jié)果表明，AUSM 式通量分裂型計(jì)算格式使IP 方法具有了良好的激波捕捉能力，這一點(diǎn)在高超聲速流場(chǎng)中體現(xiàn)得尤為明顯。而膨脹區(qū)出現(xiàn)差異的原因則有所不同：初始時(shí)，來(lái)流速度極高，噴管尾部壁面的流體快速向下游流動(dòng)，而上游來(lái)流受壁面阻擋無(wú)法及時(shí)補(bǔ)充，因此在該位置處形成了低壓區(qū)域，DSMC 方法完全依賴(lài)于粒子運(yùn)動(dòng)，在模擬時(shí)該區(qū)域單元極易出現(xiàn)內(nèi)部無(wú)分子的情況，從而出現(xiàn)與宏觀現(xiàn)實(shí)不符的“數(shù)值真空”流場(chǎng)，進(jìn)而導(dǎo)致尾部區(qū)域壓力過(guò)低，周?chē)黧w向“真空區(qū)”內(nèi)卷的情況出現(xiàn)，而本文提出的DSMC-IP 方法繼承了文獻(xiàn)[32]中IP 方法所包含的單元通量計(jì)算法，在計(jì)算無(wú)分子單元時(shí)以單元宏觀量為基礎(chǔ)構(gòu)造數(shù)值通量繼續(xù)計(jì)算，不會(huì)出現(xiàn)“數(shù)值真空”的問(wèn)題；流場(chǎng)穩(wěn)定時(shí)，膨脹區(qū)域密度不斷下降，此時(shí)該區(qū)域中的流場(chǎng)非平衡效應(yīng)增強(qiáng)，IP 方法中信息溫度的局部平衡假設(shè)使其局部計(jì)算精度略有下降。但從結(jié)果來(lái)看，這些差異并不影響整體流場(chǎng)結(jié)構(gòu)的相似，而IP 方法降低統(tǒng)計(jì)耗散、提高計(jì)算精度的優(yōu)勢(shì)同樣明顯，由此可以證明AUSM 分裂式DSMC-IP 方法在高超聲速流動(dòng)模擬中的有效性。

圖15～圖17 顯示了不同時(shí)刻兩種方法計(jì)算結(jié)果中對(duì)稱(chēng)面上的馬赫數(shù)與壓力值。結(jié)果顯示，噴管內(nèi)部的流場(chǎng)變化主要與噴管前傾角所形成的激波發(fā)展過(guò)程相關(guān)：斜激波首先發(fā)展至對(duì)稱(chēng)面，導(dǎo)致該區(qū)域內(nèi)速度驟降，壓力激增；隨后激波受對(duì)稱(chēng)面作用發(fā)生反射，而反射激波在高速來(lái)流和管內(nèi)壁面的共同作用下向出口處偏斜，在對(duì)稱(chēng)面上形成逐漸擴(kuò)張的反射波系；當(dāng)流動(dòng)最終穩(wěn)定時(shí)，來(lái)流激波與反射波發(fā)展至相互平衡的位置。通過(guò)對(duì)比可以發(fā)現(xiàn)，在兩種方法的計(jì)算結(jié)果中，激波前緣位置與數(shù)值吻合較好，進(jìn)一步證明了AUSM分裂格式對(duì)強(qiáng)激波的有效模擬能力；而激波反射點(diǎn)和反射波系的數(shù)值分布則略有差別，這也表明了在模擬壓縮效應(yīng)偏弱的復(fù)雜流場(chǎng)時(shí)，IP 方法和DSMC 方法的結(jié)果之間仍存在少許差異。但總體來(lái)看，兩種方法的結(jié)果吻合度仍比較好，以最終統(tǒng)計(jì)結(jié)果為標(biāo)準(zhǔn)，馬赫數(shù)和壓力值的最大差異均控制在總變化幅值的約5%和7%，而在DSMC 方法的實(shí)時(shí)統(tǒng)計(jì)結(jié)果中，統(tǒng)計(jì)漲落干擾比較明顯，IP 方法具有更高的計(jì)算精度。

圖15 t=6.25×10?5 s 時(shí)噴管中心線(xiàn)上的馬赫數(shù)和壓力分布Fig. 15 The Mach numbers and pressure distribution on the nozzle center line when t=6.25×10?5 s

圖16 t=1.25×10?4 s 時(shí)噴管中心線(xiàn)上的馬赫數(shù)和壓力分布Fig. 16 The Mach numbers and pressure distribution on the nozzle center line when t=1.25×10?4 s

圖17 流動(dòng)穩(wěn)定時(shí)噴管中心線(xiàn)上的馬赫數(shù)和壓力分布Fig. 17 The Mach numbers and pressure distribution on the nozzle center line when the flow is stable

圖18 展示了不同方法計(jì)算所得的噴管外側(cè)表面各特征系數(shù)分布情況，圖中L 為噴管前端平直段的長(zhǎng)度(0.1013 m)，并與文獻(xiàn)[36]中的實(shí)驗(yàn)與參考值進(jìn)行了對(duì)比。圖中清晰的顯示了噴管外側(cè)來(lái)流在其表面形成一次激波和二次激波以及相互干擾的流動(dòng)特征，由于噴管平直段對(duì)高超來(lái)流的影響類(lèi)似于平板，故一次激波前端駐點(diǎn)處壓力系數(shù)與二次激波處的峰值相比偏低，此時(shí)摩擦力的影響占主導(dǎo)；二次激波發(fā)生在噴管轉(zhuǎn)角處，傾斜壁面在法向上對(duì)來(lái)流產(chǎn)生了阻礙，來(lái)流動(dòng)能大量轉(zhuǎn)化為內(nèi)能，因此該處熱流通量和壓力系數(shù)均達(dá)到峰值，而摩擦系數(shù)雖然也增大但未能超過(guò)一次激波駐點(diǎn)處的峰值。從對(duì)比結(jié)果來(lái)看，AUSM 分裂式DSMC-IP 方法與實(shí)驗(yàn)值和參考值均符合得較好，特征量的變化趨勢(shì)和范圍基本一致，僅特征點(diǎn)的位置分布略有差異，其中壓力系數(shù)和熱流通量與實(shí)驗(yàn)值相比，分布差異最大約占總變化范圍的3%和2.2%；摩擦阻力與參考值相比，差異最大約5.3%，均在可接受的誤差范圍內(nèi)，據(jù)此可以證明本文提出的DSMC-IP 方法的準(zhǔn)確性。

圖18 流動(dòng)穩(wěn)定時(shí)噴管外側(cè)表面特征系數(shù)分布Fig. 18 The characteristic coefficients distribution on the nozzle outside surface when the flow is stable

表2 展示了兩種方法的計(jì)算時(shí)間比較。兩種方法從流動(dòng)初始到流場(chǎng)穩(wěn)定共計(jì)算25000 時(shí)間步，DSMC 方法最終收斂時(shí)共統(tǒng)計(jì)10000 時(shí)間步，IP 方法統(tǒng)計(jì)3000 時(shí)間步。從表中可以看出，IP 方法從流動(dòng)初始到流場(chǎng)穩(wěn)定的計(jì)算時(shí)間比DSMC 方法的計(jì)算時(shí)間約高47.9%，總計(jì)算時(shí)間多約16.7%，單位時(shí)間步的平均計(jì)算耗時(shí)約為DSMC 方法的1.45 倍；在收斂統(tǒng)計(jì)時(shí)間方面，IP 方法消耗時(shí)間更少，約為DSMC 方法耗時(shí)的43.8%。該算例中來(lái)流速度極高，流場(chǎng)進(jìn)入穩(wěn)態(tài)的時(shí)間較短，因此IP 方法的總計(jì)算耗時(shí)與圓柱繞流相比略低，但兩種方法計(jì)算時(shí)間的分布趨勢(shì)和大小比較基本一致。由此可見(jiàn)，AUSM 分裂式IP 方法由于自身算法格式比DSMC 方法復(fù)雜，因而總計(jì)算耗時(shí)偏高，但由于IP 方法本身具有良好的統(tǒng)計(jì)性能，該方法不但能夠獲得實(shí)時(shí)計(jì)算精度較高的結(jié)果，而且流場(chǎng)穩(wěn)定后的統(tǒng)計(jì)時(shí)間比DSMC 方法低，在該方面能夠減少部分不必要的時(shí)間消耗。

表2 帶擴(kuò)張角的噴管繞流中DSMC 與IP 方法的計(jì)算時(shí)間比較Table 2 Comparison of calculation time between DSMC and IP method in the flow around nozzle with divergent angel

3 結(jié)論

本文針對(duì)DSMC 方法中統(tǒng)計(jì)耗散嚴(yán)重而傳統(tǒng)DSMC-IP 方法缺乏有效模擬強(qiáng)激波能力的問(wèn)題，采用AUSM 分裂格式對(duì)DSMC-IP 方法進(jìn)行優(yōu)化，提出一種兼具統(tǒng)計(jì)精度與高超聲速稀薄流動(dòng)模擬能力的優(yōu)化型DSMC-IP 方法。通過(guò)對(duì)高速稀薄流場(chǎng)的模擬對(duì)比和分析，得出以下結(jié)論：

(1)在超聲速和高超聲速稀薄流場(chǎng)中，AUSM分裂式DSMC-IP 方法不但能夠降低結(jié)果中的統(tǒng)計(jì)耗散，提高計(jì)算精度，而且能夠?qū)α鲌?chǎng)中的激波進(jìn)行實(shí)時(shí)而有效的模擬，計(jì)算結(jié)果與實(shí)驗(yàn)值和參考值相符較好，主要流動(dòng)特征均能清晰展現(xiàn)，證明了該算法的可靠性和有效性。

(2)由于計(jì)算格式的不同，AUSM 分裂式IP方法在模擬壓縮程度相對(duì)較低且宏觀物理量梯度變化較大的流場(chǎng)時(shí)，計(jì)算精度因信息溫度的局部平衡假設(shè)而略有下降，但其結(jié)果的整體趨勢(shì)仍與DSMC 方法標(biāo)準(zhǔn)結(jié)果符合的較好；而在激波中心、近壁面等重要區(qū)域和壓縮效應(yīng)更為突出的流場(chǎng)中，計(jì)算結(jié)果與標(biāo)準(zhǔn)結(jié)果明顯更加吻合。據(jù)此可以說(shuō)明，AUSM 分裂式DSMC-IP 方法的計(jì)算精度基本滿(mǎn)足高超飛行器的設(shè)計(jì)要求，能夠?qū)﹃P(guān)鍵區(qū)域的流動(dòng)傳熱性能進(jìn)行準(zhǔn)確預(yù)測(cè)。

(3)在計(jì)算耗時(shí)方面，從本文測(cè)試結(jié)果來(lái)看，IP 方法單位時(shí)間步的平均耗時(shí)約為DSMC 方法的1.5 倍～1.6 倍，但由于IP 方法統(tǒng)計(jì)步數(shù)相對(duì)更少，其消耗在收斂統(tǒng)計(jì)過(guò)程中的時(shí)間更低，約為DSMC 方法的40%～43%。綜上所述，DSMC-IP 方法的本質(zhì)是在傳統(tǒng)DSMC 方法的基礎(chǔ)上增加了分子信息量的計(jì)算過(guò)程，因此其在提高實(shí)時(shí)計(jì)算精度的同時(shí)也必然帶來(lái)計(jì)算耗時(shí)的增加，但從測(cè)試結(jié)果來(lái)看，總計(jì)算耗時(shí)的增加幅度在可接受的范圍內(nèi)，且由于IP 方法具有統(tǒng)計(jì)精度更高的優(yōu)勢(shì)，收斂統(tǒng)計(jì)時(shí)反而可以減少部分不必要的時(shí)間開(kāi)銷(xiāo)。在大型的工程計(jì)算中，可將并行計(jì)算方法與DSMC-IP 方法相配合，提高計(jì)算速度。