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

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

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

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

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

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

1 計(jì)算方法

1.1 DSMC-IP 方法概述

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

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

雖然DSMC-IP 方法最初保存的信息量僅有速度，而且其信息量演化的理論模型也不明確，但該方法“采用平均量(即信息量)代替隨機(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)動過程中隨分子正常遷移而發(fā)生位置變化；當(dāng)進(jìn)行碰撞時需采用碰撞模型處理信息量的變化過程，Sun 等在文獻(xiàn)[22]中提出了一種唯象論模型，為信息量在碰撞過程中發(fā)生的變化提供了計(jì)算方法。將運(yùn)動與碰撞項(xiàng)剝離后，信息量控制方程僅剩右端關(guān)聯(lián)項(xiàng)，此時，式(6)和式(7)可寫為如下形式：

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

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

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

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

式中：E為總能，即流體微元的內(nèi)能與動能之和，而根據(jù)理想氣體方程，壓力項(xiàng)p/ρ 等于RT，因此如按式(14)將AUSM 格式中的能量對流通量進(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)動項(xiàng)和關(guān)聯(lián)項(xiàng)共同描述了對流輸運(yùn)對宏觀信息量的影響，其中運(yùn)動項(xiàng)主要對應(yīng)于AUSM 通量中的對流部分，如ρucMan和ρEcMan等，表示流體運(yùn)動對物理量的改變，在IP 方法中表現(xiàn)為信息量在分子遷移中的變化；而關(guān)聯(lián)項(xiàng)則對應(yīng)于AUSM 通量中的壓力部分，反映了宏觀壓力對局部流場的影響，在IP 方法中表現(xiàn)為信息量受相鄰單元“類壓力”和“類壓力熱流”影響而發(fā)生的改變。據(jù)此，本文采用AUSM格式對IP 方法中的各類關(guān)聯(lián)項(xiàng)通量計(jì)算格式進(jìn)行重構(gòu)，使其能夠更好的適應(yīng)高超稀薄流場的模擬要求。

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

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

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

2 數(shù)值算例

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

本文首先采用AUSM 分裂式DSMC-IP 方法對超聲速圓柱繞流流場進(jìn)行數(shù)值模擬。該算例來自文獻(xiàn)[35]，圓柱半徑為0.05 m，其表面溫度為500 K，來流氣體為氮?dú)?，來流速?412.5 m/s，溫度300 K，分子轉(zhuǎn)動自由度為2，粘性系數(shù)的溫度指數(shù)為0.74，來流氣體分子數(shù)密度為1.0×1020，其分子自由程約0.0133 m，Kn數(shù)為0.13，屬于典型的稀薄流動。計(jì)算網(wǎng)格單元總數(shù)為3270，模擬分子總數(shù)約116 萬，網(wǎng)格與模型如圖1 所示。本文計(jì)算時間步取1.0×10?8s，每100 時間步統(tǒng)計(jì)一次，計(jì)算至流場基本穩(wěn)定時進(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 方法不同時刻的統(tǒng)計(jì)結(jié)果，并將計(jì)算結(jié)果進(jìn)行了對比。兩種方法的結(jié)果都完整的展示了超聲速繞流流場的發(fā)展過程：高速來流氣體首先在圓柱迎風(fēng)表面受壓縮形成激波區(qū)域，隨后繞經(jīng)圓柱表面在其后方形成膨脹尾流，激波區(qū)與膨脹區(qū)在來流的推動下逐漸發(fā)展，最終相互融合，形成穩(wěn)定的繞流流場。通過對比不難發(fā)現(xiàn)，IP 方法計(jì)算結(jié)果的實(shí)時統(tǒng)計(jì)精度明顯優(yōu)于DSMC 方法，即使在未經(jīng)大量統(tǒng)計(jì)的實(shí)時統(tǒng)計(jì)結(jié)果中，IP 方法結(jié)果中的等值線和區(qū)域均清晰而光滑；DSMC 方法的結(jié)果則呈現(xiàn)出大量的點(diǎn)狀和塊狀結(jié)構(gòu)，其等值線也呈現(xiàn)出較大的擾動，即使經(jīng)過一定時間的統(tǒng)計(jì)，激波邊緣的溫度等值線仍帶有統(tǒng)計(jì)漲落的干擾痕跡，據(jù)此可以證明本文提出的IP 方法可有效降低超聲速流場中的統(tǒng)計(jì)耗散，提高實(shí)時計(jì)算結(jié)果的精度。

圖2 t=1.25×10?4 s 時流場中宏觀量計(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 時流場中宏觀量計(jì)算結(jié)果Fig. 3 The computational results of macro quantities in flow field when t=2.5×10?4 s

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

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

圖5 t=1.25×10?4 s 時駐點(diǎ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 時駐點(diǎn)線上的馬赫數(shù)和壓力分布Fig. 6 The Mach numbers and pressure distribution on the stagnation line when t=2.5×10?4 s

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

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

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

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

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

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

圖9 高超聲速噴管模型結(jié)構(gòu) /mmFig. 9 Structure of the hypersonic nozzle

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

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

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

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

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

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

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

圖15 t=6.25×10?5 s 時噴管中心線上的馬赫數(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ù)和壓力分布Fig. 16 The Mach numbers and pressure distribution on the nozzle center line when t=1.25×10?4 s

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

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

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

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

3 結(jié)論

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

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

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

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