稀薄流中基于碰撞分子特定次序的DSMC改進(jìn)方法

張 利,史金光,韓 宇

(1.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094;2.海軍駐沈陽(yáng)彈藥專業(yè)軍事代表室,遼寧 沈陽(yáng) 110045)

現(xiàn)代戰(zhàn)爭(zhēng)對(duì)彈箭突防能力提出了更高的要求,彈箭的射程越來(lái)越遠(yuǎn),飛行高度越來(lái)越高。隨著飛行高度的增加,大氣密度隨之降低,當(dāng)空氣密度降低到一定程度時(shí),氣體的間斷粒子效應(yīng)使得連續(xù)介質(zhì)假設(shè)不再成立,無(wú)法采用傳統(tǒng)的N-S方程求解彈箭的氣動(dòng)力特性,而需要采取稀薄氣體動(dòng)力學(xué)[1]的方法進(jìn)行研究。Bird于1963年提出了采用直接模擬蒙特卡洛(direct simulation Monte Carlo,DSMC)方法解決高空飛行器過(guò)渡流區(qū)的氣動(dòng)特性問(wèn)題[2-4]。DSMC方法的主要思想是:使用大量模擬分子代替真實(shí)氣體分子,在計(jì)算機(jī)中存儲(chǔ)模擬分子的位置坐標(biāo)、速度和內(nèi)能等信息,這些運(yùn)動(dòng)信息隨著模擬分子的運(yùn)動(dòng)、與邊界的作用以及分子間碰撞而改變,通過(guò)空間網(wǎng)格對(duì)模擬分子的諸量求和以得到流場(chǎng)的宏觀量。DSMC方法已經(jīng)得到了實(shí)驗(yàn)驗(yàn)證。超聲速鈍頭體彈箭前方的繞流流場(chǎng)中,氣體分子碰撞將導(dǎo)致各種化學(xué)反應(yīng),如氮?dú)狻⒀鯕?、一氧化氮、氮原子、氧原子等組分的離解、復(fù)合、置換等復(fù)雜的化學(xué)反應(yīng)[5-10],這些氣體分子間的化學(xué)反應(yīng)會(huì)對(duì)彈箭的受力、受熱分析產(chǎn)生不可忽視的影響。氣體分子碰撞時(shí)可能發(fā)生的化學(xué)反應(yīng)并不唯一,比如:雙原子分子AB與雙原子分子CD進(jìn)行碰撞后,既可能是AB發(fā)生離解生成2個(gè)新的單原子分子,也可能是CD發(fā)生離解形成2個(gè)新的單原子分子,等等。但事實(shí)上,任何碰撞分子可能發(fā)生的化學(xué)反應(yīng)都是唯一的,化學(xué)反應(yīng)類型與碰撞對(duì)分子的次序是一一對(duì)應(yīng)的關(guān)系。因此,本文在Bird的氣體分子化學(xué)反應(yīng)DSMC模擬程序的基礎(chǔ)上,提出了一種化學(xué)反應(yīng)類型取決于碰撞對(duì)中分子特定次序的DSMC方法,并利用FORTRAN90[11]編制了DSMC仿真程序,對(duì)過(guò)渡區(qū)域高超聲速圓柱繞流進(jìn)行模擬分析,驗(yàn)證了本文方法的可行性和合理性。

1 空氣中的化學(xué)反應(yīng)類型

高超聲速稀薄氣體流動(dòng)中存在著多種不同形式的熱化學(xué)非平衡現(xiàn)象。本文研究的是五組元(N2,O2,N,O,NO)氣體分子間的離解、置換和復(fù)合3種化學(xué)反應(yīng)類型34種化學(xué)反應(yīng)過(guò)程:

N2+X2N+X
O2+X2O+X
NO+XN+O+X
N2+ONO+N
NO+OO2+N

(1)

式中:催化物X可以是N2,O2,N,O,NO。對(duì)于化學(xué)反應(yīng)模擬,Bird提出了化學(xué)反應(yīng)概率模型,當(dāng)組分A的分子與組分B的分子發(fā)生碰撞時(shí),以一定的概率發(fā)生化學(xué)反應(yīng)。通常引入反應(yīng)截面σR,而σR與總碰撞截面σT的比值σR/σT代表彈性碰撞導(dǎo)致化學(xué)反應(yīng)的概率,并稱為位阻因子。同時(shí),碰撞中的總能量應(yīng)超過(guò)分子的活化能。

①離解反應(yīng)。空氣中的離解反應(yīng)有:

O2+X→O+O+X
NO+X→N+O+X
N2+X→N+N+X

(2)

式中:催化物X可以是N2,O2,N,O,NO,所以共有15種離解反應(yīng)。

離解反應(yīng)在DSMC模擬中的實(shí)現(xiàn)通常是利用Larsen-Bergnakke的振動(dòng)松弛理論,當(dāng)雙原子分子振動(dòng)級(jí)數(shù)激發(fā)到大于該分子離解能所對(duì)應(yīng)的振動(dòng)級(jí)數(shù),則認(rèn)為該雙原子分子發(fā)生了離解,并標(biāo)記該雙原子分子以進(jìn)行離解反應(yīng)的計(jì)算。

②復(fù)合反應(yīng)?？諝庵械膹?fù)合反應(yīng)主要有:

N+N+X→ N2+X
O+O+X→ O2+X
N+O+X→NO+X

(3)

復(fù)合反應(yīng)在DSMC模擬中的實(shí)現(xiàn)通常是利用配分函數(shù)與細(xì)致平衡原理確定抽樣概率,并根據(jù)積累概率法判斷復(fù)合反應(yīng)是否發(fā)生,當(dāng)網(wǎng)格中的累計(jì)概率大于1時(shí),即認(rèn)為發(fā)生一次復(fù)合反應(yīng),標(biāo)記該網(wǎng)格中的2個(gè)原子以進(jìn)行復(fù)合反應(yīng)的計(jì)算,并將該網(wǎng)格中的累計(jì)概率減去1。計(jì)算公式為

(4)

式中:a1,b1為復(fù)合反應(yīng)常數(shù),T為溫度,n為數(shù)密度。

③置換反應(yīng)?？諝庵械慕粨Q原子化學(xué)反應(yīng)都是雙原子分子和單原子分子間進(jìn)行碰撞交換原子?？諝庵械闹脫Q反應(yīng)有以下4種類型:

N2+O→NO+N
NO+N→N2+O
NO+O→O2+N
O2+N→NO+O

(5)

置換反應(yīng)在DSMC模擬中的實(shí)現(xiàn)通常是利用配分函數(shù)與細(xì)致平衡原理確定置換反應(yīng)抽樣概率,并根據(jù)取舍法判斷是否發(fā)生離解反應(yīng)。若抽樣概率大于隨機(jī)數(shù)Rf,則發(fā)生置換反應(yīng),并標(biāo)記該雙原子分子和自由原子,若抽樣概率小于隨機(jī)數(shù)Rf,則進(jìn)行下一碰撞對(duì)的判斷。計(jì)算公式為

(6)

式中:a2,b2為置換反應(yīng)常數(shù)。

2 基于碰撞分子特定次序的DSMC改進(jìn)方法

在DSMC方法中,化學(xué)反應(yīng)發(fā)生在分子對(duì)的碰撞過(guò)程中,Bird給出的方法是逐一判斷碰撞對(duì)分子所有可能發(fā)生的化學(xué)反應(yīng)。值得注意的是,分子碰撞時(shí)可能發(fā)生的化學(xué)反應(yīng)并不唯一,沒(méi)有必要逐一判斷碰撞對(duì)所有可能的化學(xué)反應(yīng),且該過(guò)程效率低下。為此,本文提出限定碰撞對(duì)中分子的特定次序的DSMC改進(jìn)方法,提高計(jì)算效率。

①判斷碰撞對(duì)中第1個(gè)分子是否為雙原子分子,若不是雙原子分子,則轉(zhuǎn)到②;若是雙原子分子,根據(jù)抽樣反應(yīng)概率判斷該分子是否發(fā)生離解反應(yīng),若發(fā)生了離解,則標(biāo)記該碰撞對(duì),若未發(fā)生離解反應(yīng),則進(jìn)行下一碰撞對(duì)的判斷。

②第1個(gè)分子為單原子分子,判斷第2個(gè)分子是否為雙原子分子,若不是雙原子分子,則轉(zhuǎn)到③;若是雙原子分子,則根據(jù)抽樣反應(yīng)概率判斷第2個(gè)分子是否發(fā)生置換反應(yīng),如果發(fā)生了置換,則標(biāo)記該碰撞對(duì),如果未發(fā)生置換,則進(jìn)行下一碰撞對(duì)的判斷。

③在上述步驟中判斷出碰撞對(duì)均為單原子分子,則利用抽樣概率判斷該碰撞對(duì)是否發(fā)生復(fù)合反應(yīng),如果發(fā)生了單原子復(fù)合,則標(biāo)記該碰撞分子對(duì),如果未發(fā)生置換反應(yīng),則進(jìn)行下一碰撞對(duì)的判斷。

仿真過(guò)程如圖1所示。圖中,Rf為隨機(jī)數(shù),Pcol,Prec,Pdis,Pexch分別表示碰撞概率、復(fù)合反應(yīng)概率、離解反應(yīng)概率和置換反應(yīng)概率。

圖1 化學(xué)反應(yīng)的DSMC流程圖

3 數(shù)值算例

3.1 計(jì)算模型與計(jì)算條件

現(xiàn)代戰(zhàn)爭(zhēng)對(duì)彈箭突防能力提出了更高的要求,彈箭的射程越來(lái)越遠(yuǎn),飛行高度越來(lái)越高。隨著彈箭飛行高度的增加,流場(chǎng)的稀薄效應(yīng)隨之增強(qiáng),彈箭的受力、受熱分析無(wú)法采用連續(xù)流理論,需要采取稀薄氣體動(dòng)力學(xué)的方法。本節(jié)在直角坐標(biāo)網(wǎng)格框架下,在間斷能級(jí)的Larsen-Borgnakke碰撞模型和Bird的化學(xué)反應(yīng)概率模型基礎(chǔ)上,利用DSMC方法模擬了彈箭飛行速度為7.5 km/s,飛行高度h分別為80 km,85 km,90 km時(shí)的圓柱繞流,為后續(xù)彈箭在稀薄流區(qū)的受力、受熱分析提供理論基礎(chǔ)。本文來(lái)流條件參考了文獻(xiàn)[12],在80 km和85 km處,來(lái)流組分摩爾分?jǐn)?shù)為x(N2)=0.763,x(O2)=0.237,在90 km處,x(N2)=0.788,x(O2)=0.209,x(O)=0.003,物面溫度均為1 000 K。

3.2 計(jì)算結(jié)果及分析

將壁面處的熱流密度模擬結(jié)果無(wú)量綱化,得到熱流系數(shù)Ch的表達(dá)式:

(7)

式中:q為熱流密度,ρ為來(lái)流密度,u∞為來(lái)流速度。

圖2是本文計(jì)算的沿圓柱表面熱流系數(shù)Ch和沿駐點(diǎn)線流場(chǎng)溫度T的分布曲線。圖2(a)～2(c)中原點(diǎn)處為駐點(diǎn),橫坐標(biāo)是圓柱表面不同位置θ,縱坐標(biāo)是熱流系數(shù)Ch;圖2(d)中橫坐標(biāo)是距壁面的距離d。表1為本文數(shù)值計(jì)算結(jié)果和文獻(xiàn)[12]的比較,結(jié)果表明,在80 km,85 km和90 km處,熱流系數(shù)峰值的誤差分別1.1%,2.0%和1.9%,結(jié)果吻合得較好,并且熱流系數(shù)峰值隨高度增加而增大,其值從80 km處的0.380增大到90 km處的0.632。從熱流系數(shù)曲線還可以看出,在駐點(diǎn)附近熱流密度有最大值,過(guò)了駐點(diǎn)后熱流密度逐漸減小,在尾部熱流密度非常小,接近于0。圖3是壁面熱流系數(shù)Ch和摩阻系數(shù)Cf隨高度變化曲線,可以看出熱流密度和摩阻系數(shù)隨高度變化比較敏感。

表1 不同高度圓柱表面熱流密度和熱流系數(shù)峰值

圖4(a)是流場(chǎng)平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度沿駐點(diǎn)線的分布,可以看出壁面附近存在明顯的熱力學(xué)非平衡現(xiàn)象。由此不難得出,這一區(qū)域存在脫體的弓形激波,激波內(nèi)流場(chǎng)振動(dòng)溫度明顯低于平動(dòng)溫度和轉(zhuǎn)動(dòng)溫度,且振動(dòng)溫度在此兩者之后開始上升,這是由于振動(dòng)松弛速率要遠(yuǎn)小于平動(dòng)與轉(zhuǎn)動(dòng)松弛的速率,氣體內(nèi)分子的振動(dòng)平衡需要更多次的分子碰撞,振動(dòng)達(dá)到平衡就會(huì)比平動(dòng)和轉(zhuǎn)動(dòng)要遲一些。

圖2 計(jì)算結(jié)果比較

圖3 稀薄程度對(duì)壁面熱流系數(shù)和摩阻系數(shù)的影響

圖4 稀薄程度對(duì)駐點(diǎn)線物理量的影響

圖4(b)給出了80 km高空處沿駐點(diǎn)線各氣體組分的摩爾分?jǐn)?shù)xB(B分別代表N2,O2,NO,N)分布曲線。從圖中可以看出,靠近壁面過(guò)程中,由于氣體分子密度逐漸增加,分子碰撞頻率逐漸增大,離解反應(yīng)、復(fù)合反應(yīng)和置換反應(yīng)發(fā)生的概率增大,使得氮?dú)夂脱鯕饽柗謹(jǐn)?shù)逐漸降低,氮原子和氧原子逐漸增加,靠近壁面時(shí)冷壁效應(yīng)導(dǎo)致氮?dú)夂脱鯕饽柗謹(jǐn)?shù)略有升高。值得注意的是,來(lái)流氣體分子在流場(chǎng)溫度達(dá)到峰值前就發(fā)生了離解反應(yīng),表明在高空稀薄流條件下,高超聲速?gòu)椉诿媲胺郊げ▽痈浇嬖诿黠@的化學(xué)反應(yīng),這對(duì)彈頭部的受力、受熱有著很大影響。圖4(c)是氮原子摩爾分?jǐn)?shù)隨高度的分布情況,在90 km高空處,基本上沒(méi)有氮原子的生成,流場(chǎng)化學(xué)反應(yīng)變得很微弱,可以得出彈箭表面的化學(xué)反應(yīng)隨高度增加而減弱,飛行高度達(dá)到自由分子流區(qū)時(shí)可以不再考慮化學(xué)反應(yīng)。

4 結(jié)論

本文詳細(xì)討論了包括離解反應(yīng)、復(fù)合反應(yīng)和置換反應(yīng)在內(nèi)的3種化學(xué)反應(yīng)類型,基于Bird建議的化學(xué)反應(yīng)抽樣概率模型,給出了這3種化學(xué)反應(yīng)在DSMC方法中的實(shí)現(xiàn)過(guò)程。在Bird給出的模擬二維圓柱繞流的化學(xué)反應(yīng)DSMC程序的基礎(chǔ)上,提出了化學(xué)反應(yīng)類型取決于分子碰撞對(duì)特定次序的DSMC改進(jìn)方法,并以高超聲速圓柱繞流為例,在80 km,85 km,90 km高空時(shí)模擬了五組元混合氣體的32種化學(xué)反應(yīng),在80 km,85 km和90 km時(shí)的熱流系數(shù)峰值的誤差分別為1.1%,2.0%和1.9%,與理論數(shù)據(jù)吻合得較好,表明本文提出的改進(jìn)方法是合理可行的,且效率較高,為后續(xù)高超聲速?gòu)椉诟呖障”×鲄^(qū)的受力、受熱分析提供了理論基礎(chǔ)。