黃 飛，陳 智，程曉麗，沈 清

（中國航天空氣動力技術研究院，北京 100074）

0 引 言

DSMC方法自出現(xiàn)以來已成功應用于航天、微機電、材料加工等領域［1－5］。在DSMC模擬方法中，網(wǎng)格的主要作用有兩個方面，一是為了采樣統(tǒng)計出宏觀流場變量，另一方面主要是為了在模擬過程中選取碰撞對，即兩個相互碰撞的分子必須在同一網(wǎng)格中。研究表明，DSMC模擬中必須保證網(wǎng)格尺寸在當?shù)胤肿幼杂沙痰娜种涣考壊拍苡行ПＷC模擬結果的正確性，即在DSMC模擬過程中碰撞對的選取在空間上必須滿足一定的限制，這一要求給DSMC模擬相對較高的密度區(qū)域（自由程相對較?。砹撕艽箅y度，并且在全流場模擬過程中很難保證全流場網(wǎng)格尺寸滿足這一要求。

為了降低模擬過程中對網(wǎng)格尺寸的限制，同時在不損失精度的前提下提高計算效率，Bird［6］引入了子網(wǎng)格的概念，該方法將流場中的計算網(wǎng)格在不同方向分為若干子網(wǎng)格，流場中的碰撞對在子網(wǎng)格中選取，而背景網(wǎng)格用來進行流場的采樣統(tǒng)計，這一方法的成功應用被許多DSMC研究者所采納。然而，這一靜態(tài)子網(wǎng)格的方法很難保證子網(wǎng)格中的粒子至少在兩個以上，尤其是在對高密度區(qū)域進行較多子網(wǎng)格劃分時，為了解決大量靜態(tài)子網(wǎng)格中可能少于兩個模擬粒子的問題（尤其在三維情況下），Bird［7］引入了動態(tài)子網(wǎng)格方法，該方法根據(jù)背景網(wǎng)格中的平均分子數(shù)動態(tài)分配子網(wǎng)格數(shù)，這種方法大大改善了靜態(tài)子網(wǎng)格的缺陷，但仍然難以消除此類問題出現(xiàn)的可能性；另一方面，無論是靜態(tài)子網(wǎng)格還是動態(tài)子網(wǎng)格都一定程度上增加了碰撞對選取時此類特殊情況下的邏輯判斷。鑒于此，LeBeau［8］等人引入了虛擬子網(wǎng)格方法，該方法從網(wǎng)格中首先選取碰撞對中的第一個分子，然后遍歷其它所有分子，選擇離第一個分子最近的分子作為碰撞對中的另一個分子，此種選取方式有效避免了Bird子網(wǎng)格方法中的缺陷，然而該方法需要對網(wǎng)格中的分子數(shù)N進行有效控制，不能太大，否則其在網(wǎng)格中進行的遍歷搜索，引入了N2次方量級的搜索次數(shù)及距離計算，大大增加了計算量。國內(nèi)關于DSMC模擬技術方面的工作主要集中于IP方法［9］、并行算法［10］、位置元方法［11］、網(wǎng)格處理［12］及 LD－DSMC［13］等方面，尚未發(fā)現(xiàn)關于權因子方面的研究報道。

綜上分析可以發(fā)現(xiàn)，無論是子網(wǎng)格還是虛擬子網(wǎng)格，其本質(zhì)都是盡可能的降低選取碰撞對時兩個分子之間的空間距離，盡可能保證相對較近的分子得到有效碰撞?；诖?，本文在虛擬子網(wǎng)格的方法上引入自適應的碰撞距離思想，將相對較粗的網(wǎng)格中的碰撞對限制在合理的距離范圍內(nèi)以滿足近距離的碰撞要求，與LeBeau虛擬子網(wǎng)格方法不同的是，本文在選取碰撞對時不需要遍歷網(wǎng)格中所有的分子，因而也無需對網(wǎng)格中的模擬分子數(shù)進行多余的控制，從而能夠在滿足計算精度的要求下放寬對網(wǎng)格的限制。

1 計算方法

1.1 DSMC 方法介紹［6］

DSMC方法的主要思想是采用有限個模擬分子代替真實分子，在計算機中進行分子與分子之間、分子與物面之間的相互碰撞模擬，整個模擬中分子之間及其與物面之間不斷交換動能與內(nèi)能，待模擬足夠的時間步以后，采樣統(tǒng)計給出每個網(wǎng)格中的宏觀流場結果。計算中分子之間采用VHS碰撞模型，能量交換采用L－B模型［14］，按照Bird的能量按自由度分配原則采用取舍法進行抽樣分配物面邊界條件采用完全漫反射條件。

1.2 碰撞對選取的主要關鍵技術點

碰撞距離的求解方法：

由于本文計算代碼為非結構網(wǎng)格，故而主要以非結構四面體網(wǎng)格為例進行介紹。從前面分析可以發(fā)現(xiàn)，本文的方法關鍵在于自適應碰撞距離d的確定，d的合理選取直接關系到模擬結果的準確性。計算求解中d由下式確定：

其中，L＝Lmax／（／n），λi為網(wǎng)格中的局部自由程，由流場參數(shù)確定，Lmax為四面體中的最大邊長，n為每個等效虛擬子網(wǎng)格中的最小分子數(shù)，為可調(diào)參數(shù)，為采樣一定步數(shù)后網(wǎng)格中平均分子數(shù)，即為：

從以上可以發(fā)現(xiàn)，碰撞距離d由兩個參數(shù)L、λi確定，λi／3的限制保證了碰撞距離限制在三分之一局部自由程范圍之內(nèi)，并且λi隨流場宏觀參數(shù)的變化而改變，本文在DSMC采樣前流場每隔500步更新一次，同時對相應的自由程λi、碰撞距離d都進行更新，從而能夠有效做到碰撞距離隨流場的變化進行自適應更新；L主要是為了解決在模擬過程中局部高密度網(wǎng)格中模擬分子數(shù)相對較少時碰撞對難以選取的問題。

DSMC模擬中碰撞對選取的主要步驟：

1）給定碰撞對中第二個分子選取的循環(huán)次數(shù)限制數(shù)m及第一個分子重新選取的循環(huán)次數(shù)限制數(shù)p。

2）首先在網(wǎng)格中隨機選取碰撞對中的第一個分子A。

3）在網(wǎng)格中隨機選取第二個分子B，計算兩個分子之間的距離di。

今天，我們主要來了解一下中國古代最具代表性的十大樂器，分別是 ：琵琶、二胡、編鐘、簫、笛、瑟、琴、塤（）、笙、鼓。

4）判斷0＜di≤d是否滿足，滿足則停止搜索進行分子之間的碰撞；不滿足接著判斷循環(huán)次mi數(shù)是否滿足mi＜m，滿足返回至3），不滿足則繼續(xù)判斷pi是否滿足pi＜p，滿足則返回至2），不滿足則停止搜索，進行分子之間的碰撞。

1.3 氣動力參數(shù)定義

2 算例驗證一：平板外形

2.1 計算模型及計算條件

本文采用不同網(wǎng)格尺寸下的平板外形對計算方法進行了系統(tǒng)的計算分析。其中計算模型及尺寸如圖1所示，板長L＝1.2043m，圖中邊界面A′ABB′、C′CD′D及 A′C′FE為對稱邊界，ABDC及 ACC′A′為遠場來流條件，EFD′B′為壁面邊界，BDD′B′為真空出口邊界；來流氣體為氮氣，速度1412.5m／s，溫度300K，密度4.65×10－6，壁溫500K。

圖1 計算模型示意圖Fig.1 Model in simulation

2.2 結果分析

為了對本文所發(fā)展的虛擬子網(wǎng)格方法進行驗證計算，本文在三維計算域中共選取了三種不同尺寸的網(wǎng)格進行了計算分析，其三組網(wǎng)格數(shù)分別為10850、74897及334275個單元，計算中的全場模擬分子數(shù)約分別為25萬、120萬及350萬。其中第一組網(wǎng)格為最粗的網(wǎng)格，對其分別采用虛擬子網(wǎng)格方法與不采用該方法兩種工況進行了計算分析。

圖2為計算給出的流場壓力等值線云圖，從中可以看出超聲速來流在平板前形成了斜激波，激波邊界層的相互干擾使得波后壓力升高，激波后在靠近前緣點附近壓力升至最高，之后沿壁面逐漸降低，在出口處受真空邊界的影響，壓力減小變快。圖3為不同網(wǎng)格下的物面壓力系數(shù)分布結果。從中可以看出，物面壓力系數(shù)的整體分布規(guī)律與前面所述規(guī)律類似，在前緣附近先升至最高，之后沿板面逐漸降低，在板尾部降低最快。從中還可發(fā)現(xiàn)，相對較密的后兩組網(wǎng)格所得結果與Bird結果吻合較好，而10850個計算單元下的粗網(wǎng)格在不采用虛擬子網(wǎng)格方法時的計算結果明顯偏離Bird所得計算結果。當加入子網(wǎng)格方法后此粗網(wǎng)格亦能得到與Bird一致的計算結果。此外，在板前緣壓力系數(shù)結果吻合較差，這主要是由于在板前緣流場壓力梯度相對較大，相對較粗的網(wǎng)格統(tǒng)計結果難以描述這種梯度變化。

圖2 流場壓力等值線云圖（334275cell）Fig.2 Pressure contours

圖3 不同網(wǎng)格下的物面壓力系數(shù)分布Fig.3 Pressure coefficient distribution along the wall

圖4、圖5分別為不同網(wǎng)格下的物面摩擦系數(shù)及物面熱流系數(shù)的分布結果。從中可以看出，在前緣處摩阻系數(shù)、熱流系數(shù)達最大值，隨后沿壁面降低，在尾緣處由于氣體真空膨脹加速，使摩擦系數(shù)、熱流系數(shù)在尾緣處有所增加。從中亦可發(fā)現(xiàn)，相對較密的后兩組74897、334275個單元的網(wǎng)格計算結果都與Bird計算結果吻合很好，與壓力系數(shù)的吻合規(guī)律類似，10850個計算單元下的粗網(wǎng)格的計算結果在采用本文所發(fā)展的虛擬子網(wǎng)格方法之后亦能得到滿意的結果。

圖4 不同網(wǎng)格下的物面摩擦系數(shù)分布Fig.4 Friction coefficient on the surface

圖5 不同網(wǎng)格下物面熱流系數(shù)分布Fig.5 Heat flux coefficient on the surface

3 算例驗證二：鈍錐外形

3.1 計算模型及計算條件

本文算例二主要針對鈍錐外形進行了計算分析，其中計算網(wǎng)格共采用了三種不同的尺寸。計算模型尺寸如圖6所示，鈍錐頭部半徑為R＝35mm，半錐角5°，長1m。來流氣體為空氣，來流速度為5000m／s，飛行高度90km，壁溫300K。

圖6 計算模型示意圖Fig.6 Model in simulation

3.2 結果分析

計算分析中采用的三種不同尺寸的網(wǎng)格數(shù)分別為654114、128454及34173個單元，其中最后一組最粗的網(wǎng)格采用本文所發(fā)展的虛擬子網(wǎng)格與不采用該方法兩種工況進行了計算分析。圖7、圖8分別給出了流場壓力云圖及不同網(wǎng)格尺寸下的表面壓力分布結果，從中可以看出在鈍錐頭部形成了明顯的激波結構，在駐點區(qū)域壓力較高，隨后沿錐面降低，在駐點下游附近由于膨脹作用，壓力有所下降，而在錐尾部區(qū)域由于出口邊界為真空邊界條件，流動在此區(qū)域加速膨脹使得壓力在尾部有所下翹；從中還可發(fā)現(xiàn)，在粗網(wǎng)格下采用本文所發(fā)展的虛擬子網(wǎng)格方法能夠得到與密網(wǎng)格一致的結果，而不采用該方法時，較粗網(wǎng)格下的結果明顯偏離其它密網(wǎng)格的結果。

圖7 流場壓力云圖（654114cell）Fig.7 Pressure contours

圖8 不同網(wǎng)格下的物面壓力分布Fig.8 Pressure on the surface at different grid sizes

圖9、圖10分別給出了不同網(wǎng)格尺寸下表面熱流及摩阻分布，從中可以看出，相對較密的前兩組網(wǎng)格654114、128454個單元的網(wǎng)格結果一致，而采用虛擬子網(wǎng)格方法34173個單元的粗網(wǎng)格亦能得到與密網(wǎng)格一致的熱流、摩阻分布結果，同壓力類似，不采用虛擬子網(wǎng)格方法的粗網(wǎng)格單元很難與密網(wǎng)格結果吻合。通過鈍錐外形的算例分析進一步證明了本文所發(fā)展的方法能夠在相對較粗的網(wǎng)格下取得相對較為滿意的結果。

圖9 不同網(wǎng)格下的物面熱流分布Fig.9 Heat flux on the surface at different grid sizes

圖10 不同網(wǎng)格下的摩阻分布Fig.10 Friction on the surface at different grid sizes

4 結 論

本文針對DSMC模擬過程中的碰撞特征，經(jīng)過計算分析給出了一種可進一步提高DSMC方法模擬效率的自適應虛擬子網(wǎng)格方法。通過在粗網(wǎng)格下的計算分析可以看出，本文所發(fā)展的方法能夠在相對較大的粗網(wǎng)格下得到與密網(wǎng)格一致的結果。本文的方法一定程度上放寬了DSMC模擬中的網(wǎng)格尺度，有效降低了全流場的模擬分子數(shù)，從而節(jié)省了內(nèi)存，提高了計算效率。

［1］SHEN Q.DSMC method and the calculation of rarefied gas flow［J］.AdvanceinMechanics，1996，26（1）：1－13.（in Chinese）沈青.DSMC方法與稀薄氣體計算的發(fā)展［J］.力學進展，1996，26（1）：1－13.

［2］WU Q F，CHEN W F.DSMC method of nonequilibrium thermodynamics and chemical reaction in the high temperature rarefied gas flow［M］.Changsha：National University of Defense Technology Press，1999.（in Chinese）吳其芬，陳偉芳.高溫稀薄氣體熱化學非平衡流動的DSMC方法［M］.長沙：國防科技大學出版社，1999.

［3］RAULT D F G.Study of shuttle body flap effectiveness at high altitudes using DSMC［R］.AIAA－94－2021，1994

［4］HUANG F，ZHAO B，CHENG X L，et al.Real gas effects of reentry vehicle aerodynamics under hypersonic［J］.Journalof Astronautics，2014，35（3）：283－290.（in Chinese）黃飛，趙波，程曉麗，等.返回器高空稀薄氣動特性的真實氣體效應研究［J］.宇航學報，2014，35（3）：283－290.

［5］LVANOV M，GIMELSHEIN S.Influence of real gas effects on control surfaces efficiency at high flight altitudes［R］.AIAA 93－5116.

［6］BIRD G A.Molecular gas dynamics and direct simulation of gas flow［M］.London：OxfordUniversity Press，1994，438－451.

［7］BIRD G A.Forty years of DSMC，and Now？［A］.22ndInternational Symposium of Rarefied Gas Dynamics［C］.Sydney Australia，2000.

［8］LEBEAU G J，BOYLES K A，LUMPKIN F E.Virtual subcells for the direct simulation Monte Carlo method［R］.AIAA 2003－1031.

［9］FAN J，SHEN C.Micro－Scale gas flows［J］.AdvanceinMechanics，2002，32（3）：321－336.（in Chinese）樊菁，沈青.微尺度氣體流動［J］.力學進展，2002，32（3）：321－336.

［10］WANG X D，WU Y Z，XIA J，et al.A paprallel algorithm of 3Dunstructured DSMC method and its application［J］.Journal ofAstronautics，2007，28（6）：1500－1505.（in Chinese）王學德，伍貽兆，夏健，等.三維非結構網(wǎng)格DSMC并行算法及應用研究［J］.宇航學報，2007，28（6）：1500－1505.

［11］FAN J，LIU H L，SHEN C，et al.A molecular reflection deterministic criterion used in the position element algorithm of direct statistical simulation［J］.ACTAAerodynamicaSinica，2000，18（2）：180－187.（in Chinese）樊菁，劉宏立，沈青等.直接統(tǒng)計模擬位置元方法中的分子表面反射確定論判據(jù)［J］.空氣動力學學報，2000，18（2）：180－187.

［12］LIANG J，YAN C，DU B Q.An algorithm study of three－dimensional DSMC simulation based on two－level cartesian coordinates grid structure［J］.ACTAAerodynamicaSinica，2010，28（4）：466－471.（in Chinese）梁杰，閻超，杜波強.基于兩級直角網(wǎng)格結構的三維DSMC算法研究［J］.空氣動力學學報，2010，28（4）：466－471.

［13］SU W，HE X Y，CAI G B.Extension of the low diffusion particle method for near－continuum two－phase flow simulations［J］.ChineseJournalofAeronautics，2013，26（1）：37－46.

［14］BORGANOFF C，LARSEN P S.Statistical collision model for Monte Carlo simulation of polyatomic gas mixture［J］.Journal ofComputationalPhysics，1975，18：405－420.