張 利,史金光,韓 宇

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

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

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

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

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

(1)

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

①離解反應(yīng)?？諝庵械碾x解反應(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模擬中的實現(xiàn)通常是利用Larsen-Bergnakke的振動松弛理論,當(dāng)雙原子分子振動級數(shù)激發(fā)到大于該分子離解能所對應(yīng)的振動級數(shù),則認(rèn)為該雙原子分子發(fā)生了離解,并標(biāo)記該雙原子分子以進行離解反應(yīng)的計算。

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

(4)

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

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

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

(5)

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

(6)

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

2 基于碰撞分子特定次序的DSMC改進方法

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

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

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

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

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

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

3 數(shù)值算例

3.1 計算模型與計算條件

現(xiàn)代戰(zhàn)爭對彈箭突防能力提出了更高的要求,彈箭的射程越來越遠,飛行高度越來越高。隨著彈箭飛行高度的增加,流場的稀薄效應(yīng)隨之增強,彈箭的受力、受熱分析無法采用連續(xù)流理論,需要采取稀薄氣體動力學(xué)的方法。本節(jié)在直角坐標(biāo)網(wǎng)格框架下,在間斷能級的Larsen-Borgnakke碰撞模型和Bird的化學(xué)反應(yīng)概率模型基礎(chǔ)上,利用DSMC方法模擬了彈箭飛行速度為7.5 km/s,飛行高度h分別為80 km,85 km,90 km時的圓柱繞流,為后續(xù)彈箭在稀薄流區(qū)的受力、受熱分析提供理論基礎(chǔ)。本文來流條件參考了文獻[12],在80 km和85 km處,來流組分摩爾分?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 計算結(jié)果及分析

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

(7)

式中:q為熱流密度,ρ為來流密度,u∞為來流速度。

圖2是本文計算的沿圓柱表面熱流系數(shù)Ch和沿駐點線流場溫度T的分布曲線。圖2(a)～2(c)中原點處為駐點,橫坐標(biāo)是圓柱表面不同位置θ,縱坐標(biāo)是熱流系數(shù)Ch;圖2(d)中橫坐標(biāo)是距壁面的距離d。表1為本文數(shù)值計算結(jié)果和文獻[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ù)曲線還可以看出,在駐點附近熱流密度有最大值,過了駐點后熱流密度逐漸減小,在尾部熱流密度非常小,接近于0。圖3是壁面熱流系數(shù)Ch和摩阻系數(shù)Cf隨高度變化曲線,可以看出熱流密度和摩阻系數(shù)隨高度變化比較敏感。

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

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

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

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

圖4 稀薄程度對駐點線物理量的影響

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

4 結(jié)論

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