化學(xué)反應(yīng)對(duì)高超聲速稀薄流氣動(dòng)熱影響研究

屈程　鄒鑫

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

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

中圖分類(lèi)號(hào)：V411.3；O356 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1672-3791（2018）11（b）-0092-05

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

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

1 稀薄流DSMC熱流計(jì)算方法

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

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

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

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

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

2 分子搜索技術(shù)

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

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

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

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

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

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

3 DSMC算法介紹

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

4 化學(xué)反應(yīng)模型及實(shí)現(xiàn)

本文在研究化學(xué)反應(yīng)效應(yīng)對(duì)飛行器表面熱流的影響時(shí)考慮的是5組分無(wú)電離模型，在模型中考慮如下離解、置換和復(fù)合3種化學(xué)反應(yīng)。

（1）N2+M2N+M

（2）O2+M2O+M

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

（4）N2+0NO+N

（5）NO+OO2+N

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

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

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

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

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

5.1 計(jì)算條件

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

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

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

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

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

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

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

6 結(jié)語(yǔ)

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

參考文獻(xiàn)

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

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

[3] 吳其芬.高溫稀薄氣體熱化學(xué)反應(yīng)流動(dòng)的DSMC方法[M].長(zhǎng)沙：國(guó)防科技大學(xué)出版社，1999.

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