屈程　鄒鑫

摘 要：建立了基于非結(jié)構(gòu)網(wǎng)格的高超聲速稀薄流DSMC計算方法。發(fā)展了5組分有限速率化學反應模型、反應發(fā)生狀態(tài)判定方法與求解技術(shù)，發(fā)展了DSMC熱流密度高效求解方法。以高超聲速圓柱外形為研究對象，針對不同飛行高度下2組分混合氣體模型（不含化學反應）和5組分混合氣體模型（含化學反應）的流動開展了數(shù)值模擬，給出了兩種流態(tài)下的繞流流場以及熱流密度分布，比較并分析了化學反應效應對流場特性，尤其是對熱流密度的影響。

關(guān)鍵詞：化學反應 熱流密度 非結(jié)構(gòu)網(wǎng)格 DSMC

中圖分類號：V411.3；O356 文獻標識碼：A 文章編號：1672-3791（2018）11（b）-0092-05

氣動加熱問題是高超聲速稀薄流域飛行器面臨的重要問題，準確獲得飛行器的氣動熱環(huán)境能有助于研究者對此類飛行器開展熱防護方案設(shè)計。在高超聲速稀薄流域飛行器物面附近和前緣弓形激波后的流動區(qū)域，來流通過激波被急劇加熱，在此狀態(tài)下，空氣分子會發(fā)生離解、復合、交換等一系列化學反應[1-4]，Boyd研究表明，化學反應效應會影響高超聲速稀薄流的流場特征[5]。目前，類似參數(shù)下的地面試驗尚不完善，而Bird[2]提出的DSMC（Direct Simulation Monte Carlo）方法是目前公認的能夠很好地模擬稀薄流流動問題的數(shù)值方法，因此，開展高空高超聲速稀薄流化學反應效應DSMC氣動熱特性研究，對于高空高超聲速飛行器熱防護設(shè)計具有重大意義。

基于位置元方案的直角坐標網(wǎng)格利用與物面相交的子網(wǎng)格表面近似表征物面，而適體網(wǎng)格對任意復雜外形有高度貼體性，能夠提高物面氣動力、氣動熱的數(shù)值模擬精度。本文基于二維非結(jié)構(gòu)適體網(wǎng)格，編寫了考慮稀薄流化學反應效應的DSMC計算程序。針對稀薄流區(qū)高超聲速流動，建立高超聲速稀薄流熱化學反應效應影響下的DSMC熱流計算方法，同時，討論了5組分有限速率化學反應模型與求解技術(shù)，編寫了基于非結(jié)構(gòu)網(wǎng)格二維DSMC程序，針對圓柱外形繞流流場進行了數(shù)值模擬與分析。

1 稀薄流DSMC熱流計算方法

工程實際中通常采用耦合的方法通過求解物面熱流密度得到飛行器的結(jié)構(gòu)傳熱情況，進而設(shè)計出合理的熱防護方案，這就對熱流密度求解的精確度提出了很高的要求。隨著連續(xù)介質(zhì)模型的失效，高超聲速稀薄流問題中的熱流已經(jīng)不再能夠由低階的宏觀溫度來表征。此時需要對流場中的微觀粒子進行宏觀統(tǒng)計，從而得到與分子能量相關(guān)聯(lián)的通過氣體某一位置上某一小面積元的熱流通量表達式。

如圖1所示，某一面積微元的面積為dS，其單位法向矢量為，分子速度分布函數(shù)為，觀察速度在附近的中的分子，時間t到t+dt內(nèi)通過dS的分子在dt開始瞬間位于以dS為基底邊長為cdt的柱形里面，于是，dt內(nèi)穿過dS的速度在附近的分子數(shù)目為：

每個分子攜帶的能量為：，εin是與一個分子相聯(lián)系的內(nèi)自由度能量，包括轉(zhuǎn)動能和振動能。這時，dt時間內(nèi)穿過dS的熱流通量為：

通過對入射分子和反射分子的能量通量進行統(tǒng)計計算就能夠得到單位時間傳遞到飛行器單位表面積上的氣動加熱熱流，在DSMC算法中，熱流密度qw可以表示為以下形式：

其中，i代表入射分子，r代表反射分子，Δt代表加熱時間，S表示加熱面積。

2 分子搜索技術(shù)

DSMC方法需要在計算中不斷更新模擬分子的位置信息，快速而準確地模擬分子跟蹤算法能夠保證該方法的計算精確度，本文采用直線搜索技術(shù)來跟蹤流場中的模擬分子，假設(shè)模擬分子從單元ABC中的P處運動到Q處，如圖2所示。

（1）按邊循環(huán)，計算邊的始末點和Q點構(gòu)成三角形的有向面積，如果有向面積全部為正，說明Q點仍在初始單元內(nèi)，停止搜索。如果有向面積出現(xiàn)負值，說明P點在對應邊外側(cè)，判斷PQ與對應邊是否相交，如果PQ與對應邊不相交，則繼續(xù)判定初始單元的下一條邊，如果相交，判定該邊是否具有相鄰單元，如果沒有相鄰單元轉(zhuǎn)到第4步，如果該邊有相鄰單元，則向該邊相鄰的單元內(nèi)搜索，轉(zhuǎn)到第2步。

（2）利用面積元方法[4]尋找該單元除該邊以外的對應邊，如果找不到對應邊，說明Q點在該單元內(nèi)，停止搜索。如果找到對應邊，判斷PQ與對應邊是否相交，如果PQ與對應邊不相交，則繼續(xù)判定該單元的下一條邊，如果相交，判定該邊是否具有相鄰單元，如果沒有相鄰單元轉(zhuǎn)到第4步，如果該邊有相鄰單元，則向該邊相鄰的單元內(nèi)搜索，轉(zhuǎn)到第3步。

（3）重復類似第2步的方法，直至找到Q所在單元，停止搜索。

（4）根據(jù)邊的性質(zhì)判定該邊是物面還是遠場邊界。如果是遠場邊界，那么模擬分子直接越出邊界，將其刪除，停止搜索；如果是物面邊界，那么根據(jù)PQ與邊的交點得到模擬分子撞擊物面的位置，轉(zhuǎn)到第5步。

（5）調(diào)用模擬分子與物面相互作用的子程序，返回模擬分子與物面作用的準確信息，得到模擬分子所在網(wǎng)格單元，停止搜索。

3 DSMC算法介紹

DSMC算法的本質(zhì)是對稀薄流場中的分氣體子運動和氣體分子間碰撞進行解耦運算，該方法采用大量的模擬分子對真實氣體進行模擬，每個模擬分子代表特定數(shù)目的真實氣體分子。在流場模擬過程中，模擬分子與模擬分子以及物面不斷通過碰撞的方式進行能量交換，經(jīng)過一定模擬時間，流場中的分子數(shù)量趨于穩(wěn)定，此時，采用統(tǒng)計采樣的方式得到流場宏觀計算結(jié)果。網(wǎng)格在DSMC算法中能夠?qū)崿F(xiàn)對流場宏觀流動參數(shù)進行空間離散以及分子碰撞對的選擇。本文采用非結(jié)構(gòu)網(wǎng)格進行DSMC數(shù)值模擬，模擬氣體分子采用可變硬球（VHS）模型，分子碰撞對采用非時間計數(shù)器（NTC）法選取，分子對碰撞過程中的能量交換采用Larsen-Borgnakke唯象論模型處理，物面采用基于完全漫反射模型的恒溫邊界條件。

4 化學反應模型及實現(xiàn)

本文在研究化學反應效應對飛行器表面熱流的影響時考慮的是5組分無電離模型，在模型中考慮如下離解、置換和復合3種化學反應。

（1）N2+M2N+M

（2）O2+M2O+M

（3）NO+MN+0+M （4）

（4）N2+0NO+N

（5）NO+OO2+N

上述反應式中，M為催化劑，可以為5組分中任意一個，化學反應不會改變催化劑的性質(zhì)。本文認為空氣中的化學反應與非彈性碰撞中的內(nèi)能松弛過程是耦合在一起的，在DSMC模擬中采用與不同的反應類型相對應的方式判定其化學反應的發(fā)生。

（1）離解反應：本文利用L-B模型的振動松弛理論，當分子的振動級數(shù)激發(fā)到高于分子分裂能所對應的振動級數(shù)時，認為分子發(fā)生離解反應。

（2）復合反應：本文引入Bird的唯象化學反應模型[1]，利用配分函數(shù)與平衡碰撞理論確定復合反應的抽樣幾率[6]，計算公式為，其中A、B是常數(shù)，文獻[6]中給出了對于空氣中不同復合反應A、B的常用值。

（3）置換反應：本文利用文獻[5]中引入的二元碰撞理論計算碰撞發(fā)生置換反應的概率，在每次碰撞中積累每個單元格內(nèi)各置換反應的反應概率，當某種置換反應的概率累加大于1時，認為發(fā)生一次與之相對應的置換反應。

5 算例和結(jié)果分析

5.1 計算條件

采用本文發(fā)展的計算方法，針對二維高超聲速圓柱外形繞流開展了數(shù)值模擬分析。來流速度設(shè)為7500m/s，攻角為0°，壁面溫度Tw=300K，分別采用考慮化學反應效應的5組分氣體模型和不考慮空氣中化學反應的2組分氣體模型，表1給出了其他具體計算參數(shù)。圓柱的半徑為0.08m，外形幾何及網(wǎng)格邊界如圖3所示。

5.2 結(jié)果比較與分析

在80km計算高度，采用不同組分氣體模型得到的流場密度云圖如圖4所示，觀察該圖可以發(fā)現(xiàn)，受稀薄氣體效應的影響，激波的過渡區(qū)相對較大，這和連續(xù)流場有較為明顯的差異，對比兩種流態(tài)結(jié)果可以看到，采用5組分混合氣體模型得到的激波位置更貼近物面。

在80km計算高度，采用不同組分氣體模型得到的流場駐點線平動溫度如圖5所示，經(jīng)過觀察可以發(fā)現(xiàn)，采用5組分氣體模型計算的到的流場激波結(jié)果更靠近物面，本算例中，化學反應效應的使駐點線平動溫度峰值降低了約9000K，對駐點線平動溫度影響顯著。這主要是因為駐點處溫度很高，流場中的大量熱量被劇烈化學反應吸收，因此平動溫度降低，此時，受脹冷縮效應的影響，空氣分子的振動幅度隨之減小，增加了氣體的壓縮性，因此，激波離物面更近。

80km和90km計算高度下，駐點線的U向速度分布如圖5所示，經(jīng)過觀察可以發(fā)現(xiàn)，在駐點線上，速度沿來流方向不斷降低，在駐點處降低至0，同時，采用考慮流場化學反應的5組分混合氣體模型計算得到的流場激波位置更為靠近物面。由圖6（a）（b）明顯可以發(fā)現(xiàn)，在稀薄氣體效應的作用下，激波層厚度隨計算高度的增加不斷增厚。

80km和90km計算高度下，采用不同組分氣體模型得到的物面熱流密度沿軸線方向分布如圖7所示，可以看出在駐點處熱流密度均出現(xiàn)了最高值，在迎風面，熱流密度整體相對較高，沿X軸線方向熱流密度急劇下降，在背風面，熱流密度整體相對很低，這是因為相對迎風面來說，圓柱體尾部背風區(qū)氣體分子密度很低，這樣通過統(tǒng)計方法得到的壁面熱流密度也很小。對比圖7（a）和圖7（b）可以發(fā)現(xiàn)，采用兩種組分混合氣體模型得到的熱流密度變化趨勢大體一致，但是采用5組分氣體模型計算得到的物面熱流密度比采用兩組分氣體模型的結(jié)果低，同時，在90km計算高度，化學反應效應對熱流密度的影響比80km高度更弱，這表明計算高度越低，空氣越稀薄，化學反應效應對熱流密度影響有減弱的趨勢。

80km和90km計算高度下，采用不同組分氣體模型得到的圓柱外形駐點熱流密度如表2所示，觀察該表可以發(fā)現(xiàn)稀薄流化學反應效應降低了駐點處的熱流密度，同時隨著高度增加，化學反應效應對駐點熱流密度的影響有降低的趨勢，在80km計算高度，化學反應效應使駐點熱流密度降低了35.748%，在90km計算高度，化學反應效應使駐點熱流密度降低了7.358%，可以預測，隨著高度進一步增加，化學反應效應對駐點熱流密度的影響將更加微弱。

6 結(jié)語

本文開展了基于非結(jié)構(gòu)網(wǎng)格的高超聲速稀薄流化學反應效應影響下的DSMC氣動熱數(shù)值分析研究，建立了基于非結(jié)構(gòu)網(wǎng)格的高超聲速稀薄流的DSMC熱流計算方法與程序。采用本文發(fā)展的方法對圓柱外形數(shù)值算例進行了模擬，對比和分析了化學反應效應對流場特性以及熱流特性的影響，算例結(jié)果表明高溫化學反應效應會對氣動熱特性產(chǎn)生較大的影響，是高空高超聲速飛行器熱防護設(shè)計必須考慮的重要因素。

參考文獻

[1] 沈青.稀薄氣體動力學[M].北京：國防工業(yè)出版社，2003.

[2] BIRD GA. Molecular gas dynamics and the direct simulation of gas flow[M].Oxford： Clarendon Press， 1994.

[3] 吳其芬.高溫稀薄氣體熱化學反應流動的DSMC方法[M].長沙：國防科技大學出版社，1999.

[4] 王學德.高超聲速稀薄氣流非結(jié)構(gòu)網(wǎng)格DSMC及并行算法研究[D].南京航空航天大學，2006.

[5] Lain D.Boyd. Rotational and vibrational nonequilibrium effects in rarefied， hypersonic flows[J].Journal of Thermophysics and Heat Transfer，1990，4（4）：478-484.

[6] Ann B.Carlson，G.A.Bird.Implementation of a Vibrationally Linked Chemical Reaction Model for DSMC.NASA technical memorandum-109109[Z].