彭傲平李志輝 吳俊林蔣新宇

1)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,超高速空氣動(dòng)力研究所,綿陽(yáng) 621000)2)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,綿陽(yáng) 621000)3)(國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室,北京 100191)

含振動(dòng)能激發(fā)Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法驗(yàn)證與分析?

彭傲平1)李志輝1)2)3)?吳俊林1)蔣新宇1)

1)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,超高速空氣動(dòng)力研究所,綿陽(yáng) 621000)2)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,綿陽(yáng) 621000)3)(國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室,北京 100191)

振動(dòng)能激發(fā),Boltzmann方程,氣體動(dòng)理論統(tǒng)一算法,熱力學(xué)非平衡效應(yīng)

為模擬研究高溫高馬赫數(shù)下多原子氣體內(nèi)能激發(fā)對(duì)跨流域非平衡流動(dòng)的影響,將轉(zhuǎn)動(dòng)能、振動(dòng)能分別作為氣體分子速度分布函數(shù)的自變量,把轉(zhuǎn)動(dòng)能和振動(dòng)能處理為連續(xù)分布的能量模式,將Boltzmann方程的碰撞項(xiàng)分解成彈性碰撞項(xiàng)和非彈性碰撞項(xiàng),同時(shí)將非彈性碰撞按一定松弛速率分解為平動(dòng)-轉(zhuǎn)動(dòng)能松弛過程和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能松弛過程,構(gòu)造了一類考慮振動(dòng)能激發(fā)的Boltzmann模型方程,并證明了其守恒性和H定理.基于內(nèi)部能量變量對(duì)分布函數(shù)無(wú)窮積分,引入三個(gè)約化速度分布函數(shù),得到一組考慮振動(dòng)能激發(fā)的約化速度分布函數(shù)控制方程組,使用離散速度坐標(biāo)法,基于LU-SGS隱式格式和有限體積法求解離散速度分布函數(shù),建立含振動(dòng)能激發(fā)的氣體動(dòng)理論統(tǒng)一算法.通過開展高稀薄流到連續(xù)流圓柱繞流問題統(tǒng)一算法與直接模擬蒙特卡羅法模擬結(jié)果對(duì)比分析,特別是過渡流區(qū)平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)非平衡效應(yīng)對(duì)繞流流場(chǎng)與物面力熱特性的影響機(jī)制,證實(shí)了所建立的含振動(dòng)能激發(fā)的Boltzmann模型方程及氣體動(dòng)理論統(tǒng)一算法的準(zhǔn)確可靠性.

1 引 言

航天器從外層空間再入大氣層跨流域高超聲速、高溫繞流流場(chǎng)中,氣體分子的微觀自由度(平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)和電子態(tài))因受到一定程度的激勵(lì),會(huì)出現(xiàn)能量彼此傳遞,使分子和原子間發(fā)生化學(xué)和電離反應(yīng)[1].氣體的宏觀運(yùn)動(dòng)和狀態(tài)變化同相應(yīng)的微觀物理化學(xué)過程相互影響呈現(xiàn)復(fù)雜的非平衡現(xiàn)象.根據(jù)表征分子微觀自由度之間能量傳遞或組元之間進(jìn)行化學(xué)反應(yīng)的特征弛豫時(shí)間與流動(dòng)特征時(shí)間大小尺度的不同,可將非平衡流動(dòng)分為平動(dòng)和轉(zhuǎn)動(dòng)非平衡流、振動(dòng)和化學(xué)非平衡流以及電離輻射非平衡流[2].如果流動(dòng)特征時(shí)間極小或流場(chǎng)的物理量變化梯度極大,平動(dòng)、轉(zhuǎn)動(dòng)與振動(dòng)非平衡效應(yīng)表現(xiàn)為與分子性質(zhì)相關(guān)聯(lián)的氣體介質(zhì)特性,如比熱、比熱比等不再保持常數(shù),出現(xiàn)黏性、熱傳導(dǎo)和擴(kuò)散的非平衡變化,這是一種最基本且接近高超聲速再入多尺度非平衡流動(dòng)的現(xiàn)實(shí)環(huán)境.在統(tǒng)計(jì)熱力學(xué)研究中,一個(gè)氣體分子的總能量是平動(dòng)能、轉(zhuǎn)動(dòng)能、振動(dòng)能、電子能等幾種能量模式的總和,由于電子能模式過于復(fù)雜,已不能采用幾何和熱自由度來描述,且電子能激發(fā)需要大量能量,對(duì)于大部分氣體分子而言,在電子激發(fā)前分子結(jié)構(gòu)就已經(jīng)被破壞,通常在不太高的溫度下電子能一般不會(huì)激發(fā),可忽略不計(jì).所以在后面的討論中,我們假設(shè)分子只具有三種能量模態(tài):平動(dòng)能、轉(zhuǎn)動(dòng)能和振動(dòng)能.

對(duì)于不考慮內(nèi)能影響的簡(jiǎn)單氣體,僅有平動(dòng)能,氣體分子運(yùn)動(dòng)論(氣體動(dòng)理學(xué)理論)的基本方程——Boltzmann方程[3?8]能很好地描述從連續(xù)流到自由分子流各個(gè)流動(dòng)區(qū)域的氣體分子輸運(yùn)現(xiàn)象,該方程是一個(gè)高度非線性積分-微分方程,除Maxwell分布等少數(shù)幾個(gè)解析解外,幾乎不可能求出精確解.為此簡(jiǎn)化而來的Bhatnagar-Gross-Krook(BGK)模型方程[9]、橢球統(tǒng)計(jì)模型方程[10],Shakhov模型方程[11]等能較好地用于簡(jiǎn)單氣體流動(dòng)現(xiàn)象的模擬,這些模型方程結(jié)構(gòu)簡(jiǎn)單、便于數(shù)值模擬[12?19],特別是近十余年,文獻(xiàn)[13,16,19—27]從分析研究氣體分子速度分布函數(shù)與宏觀流動(dòng)量變化關(guān)系出發(fā),使用氣體分子碰撞松弛參數(shù)和當(dāng)?shù)仄胶鈶B(tài)分布函數(shù),對(duì)Boltzmann方程碰撞積分進(jìn)行物理分析與可計(jì)算建模,確立描述各流域復(fù)雜流動(dòng)輸運(yùn)現(xiàn)象統(tǒng)一的Boltzmann模型方程,提出并發(fā)展了可用于高、低不同馬赫數(shù)繞流模擬的離散速度坐標(biāo)法與動(dòng)態(tài)確定物理空間宏觀流動(dòng)量的離散速度數(shù)值積分技術(shù),先后建立起求解跨流域繞流問題的氣體動(dòng)理論統(tǒng)一算法(GKUA)與航天器再入跨流域多尺度繞流氣動(dòng)力/熱問題并行計(jì)算應(yīng)用研究平臺(tái)[19,28].在針對(duì)航天器試驗(yàn)與飛行狀態(tài)的氣動(dòng)特性研究中,氣體介質(zhì)均為雙原子氣體或者多原子氣體,隨著溫度的增加氣體分子內(nèi)能的激發(fā)將會(huì)對(duì)其輸運(yùn)系數(shù)、飛行器繞流流場(chǎng)參數(shù)的分布產(chǎn)生影響,尤其是航天器再入過程中,因高溫高速引起的氣體振動(dòng)能激發(fā)等現(xiàn)象顯著,為數(shù)值模擬研究這種非平衡流動(dòng)問題,在連續(xù)流區(qū)域通常采用基于多溫度模型的N-S方程[29,30]、稀薄流區(qū)采用基于Larsen-Bergnakke唯象論模型的直接模擬蒙特卡羅法(DSMC)方法[31?33],而過渡流區(qū)尚未有準(zhǔn)確可靠的非平衡流動(dòng)模擬方法.

根據(jù)Boltzmann方程描述氣體流動(dòng)過程中微觀分子速度分布函數(shù)基于位置空間、速度空間在任一時(shí)刻由非平衡態(tài)向平衡態(tài)的演化屬性,可以考慮從其出發(fā)進(jìn)行模型化設(shè)計(jì),研究全流域尤其是過渡流區(qū)的非平衡流動(dòng)問題.為此,Wang等[34,35]提出了一種處理具有內(nèi)能影響的多原子氣體的半經(jīng)典方法,即平動(dòng)能根據(jù)自由度按經(jīng)典方法處理,而內(nèi)自由度按量子力學(xué)觀點(diǎn)處理,由此得到了Boltzmann方程的推廣形式——王承書-烏倫貝克(WCU)方程,該方程在分子內(nèi)能為非簡(jiǎn)并態(tài)時(shí)成立.Morse[36]在Wang等研究成果的基礎(chǔ)上,將彈性碰撞與非彈性碰撞解耦用于處理能量松弛過程,依此構(gòu)造了一種類似BGK模型的考慮分子內(nèi)能間斷能級(jí)的模型方程.在這些研究中,內(nèi)能是作為單一模式處理的,即沒有區(qū)分轉(zhuǎn)動(dòng)能和振動(dòng)能.如果按量子力學(xué)觀點(diǎn)處理分子能級(jí)的躍遷,每個(gè)能級(jí)均為一個(gè)獨(dú)立的分布函數(shù),勢(shì)必會(huì)加大數(shù)值處理的難度.盡管在高溫下分子結(jié)構(gòu)具有量子效應(yīng),但在實(shí)際的氣體流動(dòng)中其量子效應(yīng)是可以忽略的[4],因此可以將分子的內(nèi)能態(tài)看成是連續(xù)的,可采用經(jīng)典熱力學(xué)的方法進(jìn)行處理.根據(jù)這種思想建立的雙原子氣體ES-BGK模型[10,37,38]、Rykov模型[39?42]等僅考慮了轉(zhuǎn)動(dòng)能的影響.但隨著溫度的進(jìn)一步增加,振動(dòng)能開始激發(fā),其對(duì)氣體熱力學(xué)特性的影響將越來越明顯,例如對(duì)空氣而言當(dāng)溫度超過600 K時(shí)振動(dòng)能的影響將不可忽略.Holway[10]在Morse研究的基礎(chǔ)上,將內(nèi)能中的轉(zhuǎn)動(dòng)能和振動(dòng)能分別按經(jīng)典方法和間斷能級(jí)處理,建立了適于多原子氣體的ES模型.Tantos等[43,44]根據(jù)量子效應(yīng)是可以忽略的思想[4,45]對(duì)多原子氣體ES模型進(jìn)行改造后用于平板間熱傳導(dǎo)研究,考察了氣體轉(zhuǎn)動(dòng)能、振動(dòng)能激發(fā)對(duì)熱傳導(dǎo)的影響.

本文根據(jù)連續(xù)能級(jí)的處理方式,在文獻(xiàn)[46—48]開展雙原子氣體含轉(zhuǎn)動(dòng)非平衡效應(yīng)Boltzmann模型方程統(tǒng)一算法研究的基礎(chǔ)上,將轉(zhuǎn)動(dòng)能、振動(dòng)能分別作為氣體分子速度分布函數(shù)中的自變量,引入轉(zhuǎn)動(dòng)能、振動(dòng)能能模,將Boltzmann方程的碰撞項(xiàng)分解成彈性碰撞項(xiàng)和非彈性碰撞項(xiàng)兩部分,同時(shí)根據(jù)Holway[10]處理碰撞過程中能量松弛的思想,將非彈性碰撞項(xiàng)按一定的松弛速率分解為平動(dòng)-轉(zhuǎn)動(dòng)能松弛和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能松弛,構(gòu)造一類考慮振動(dòng)能激發(fā)的Boltzmann模型方程.為數(shù)值求解該模型方程,引入三個(gè)約化速度分布函數(shù),得到一組考慮振動(dòng)能激發(fā)的約化速度分布函數(shù)控制方程組,在氣體動(dòng)理論統(tǒng)一算法[19?28]體系下,使用離散速度坐標(biāo)法,結(jié)合LU-SGS隱式格式和有限體積法對(duì)約化速度分布函數(shù)控制方程組進(jìn)行求解[26,27],捕捉離散速度分布函數(shù)演化更新,擬研究建立含振動(dòng)能激發(fā)的氣體動(dòng)理論統(tǒng)一算法.以圓柱繞流問題為研究對(duì)象通過開展稀薄流到連續(xù)流區(qū)統(tǒng)一算法與DSMC方法結(jié)果對(duì)比分析,特別是過渡流區(qū)平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)非平衡效應(yīng)對(duì)繞流流場(chǎng)與物面力熱特性的影響機(jī)制,對(duì)所建立的含振動(dòng)能激發(fā)的Boltzmann模型方程進(jìn)行驗(yàn)證分析.

2 建立含振動(dòng)能激發(fā)的Boltzmann模型方程

2.1 氣體分子速度分布函數(shù)對(duì)宏觀流動(dòng)量的動(dòng)態(tài)表征關(guān)系

在不考慮氣體離解和化學(xué)反應(yīng)的情況下,引入轉(zhuǎn)動(dòng)能和振動(dòng)能的氣體分子速度分布函數(shù)可定義為f(r,ξ,t;erot,evib),其中t為時(shí)間,r為空間位置坐標(biāo),ξ為分子運(yùn)動(dòng)速度;erot,evib分別為轉(zhuǎn)動(dòng)能和振動(dòng)能能模,在量子態(tài)下為對(duì)應(yīng)的各個(gè)能級(jí)能量值,這里將其視為連續(xù)能級(jí)中的點(diǎn),取值范圍均為[0,+∞),并假定二者為相互獨(dú)立的變量.宏觀流動(dòng)量可以通過對(duì)分布函數(shù)取矩的方法得到[20,21,24,49,50].設(shè)nrot,vib為能模erot,evib下的數(shù)密度,nvib為能模evib下的數(shù)密度,n為氣體流動(dòng)當(dāng)?shù)財(cái)?shù)密度,U為宏觀氣體流動(dòng)速度;Ttr,Trot,Tvib分別為平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)溫度;qtr,qrot,qvib分別為分子平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)能輸運(yùn)產(chǎn)生的熱流矢量,而總的熱流矢量qtot為三者之和;τij為脅強(qiáng)張量在i,j方向的分量;p為流場(chǎng)壓力.則各宏觀流動(dòng)量的表達(dá)式為

i,j對(duì)應(yīng)x,y,z,p=nkBTtr,qtot=qtr+qrot+qvib,其中,m為分子質(zhì)量;δrot,δvib分別為氣體分子的轉(zhuǎn)動(dòng)和振動(dòng)自由度;kB為Boltzmann常數(shù);C=ξ?U為分子隨機(jī)熱運(yùn)動(dòng)速度.

由于氣體的轉(zhuǎn)動(dòng)特征溫度Θrot很低(如O2為2.07 K,N2為2.88 K),因而在實(shí)際的氣體動(dòng)力學(xué)研究領(lǐng)域可以認(rèn)為分子的轉(zhuǎn)動(dòng)能是完全激發(fā)的,轉(zhuǎn)動(dòng)自由度δrot可以確定,對(duì)于線性分子如O2,N2,CO2等,δrot=2,對(duì)于非線性分子如H2O等,δrot=3.而振動(dòng)特征溫度Θvib較高(如O2為2256 K,N2為3371 K,CO2有三個(gè),分別為960,1930,3390 K),一般情況下振動(dòng)能難以完全激發(fā),其自由度δvib為溫度的函數(shù),假設(shè)分子具有N個(gè)振動(dòng)特征溫度即N個(gè)振動(dòng)模式,對(duì)于線性分子N=3M?5,對(duì)于非線性分子N=3M?6,這里M為分子中的原子個(gè)數(shù),則溫度T時(shí)總的振動(dòng)自由度為

其中,Θv,i為第i振動(dòng)能模式下的振動(dòng)特征溫度.

在分布函數(shù)中引入轉(zhuǎn)動(dòng)能和振動(dòng)能能模后,需要首先確定平衡態(tài)時(shí)的分布函數(shù)fM.文獻(xiàn)[3,4]中給出了轉(zhuǎn)動(dòng)能和振動(dòng)能能模的分布函數(shù)為

結(jié)合麥克斯韋速度分布函數(shù),則有[51]

同時(shí)定義平動(dòng)-轉(zhuǎn)動(dòng)實(shí)效溫度Tt,r和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)實(shí)效溫度Tov(即流場(chǎng)等效溫度)分別為

(5)式中,Tov可通過迭代的方法進(jìn)行求解.

在實(shí)效溫度的定義中引入了平動(dòng)溫度Ttr、轉(zhuǎn)動(dòng)溫度Trot和振動(dòng)溫度Tvib這三種溫度,因此可以建立一種考慮振動(dòng)能激發(fā)的三溫度Boltzmann模型方程.參照現(xiàn)有多原子氣體動(dòng)力學(xué)模型方程如Rykov等[39?41]模型方程,Holway[10]的間斷能級(jí)模型方程等的構(gòu)造原理,將分子運(yùn)動(dòng)的松弛過程簡(jiǎn)化為平動(dòng)能松弛、平動(dòng)-轉(zhuǎn)動(dòng)能松弛和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能松弛,而忽略轉(zhuǎn)動(dòng)能間、振動(dòng)能間以及轉(zhuǎn)動(dòng)-振動(dòng)能間等的松弛.由此可以將Boltzmann模型方程寫成如下形式:

其中,

這里,υtot=Pr·nkBTtr/μ,Zrot和Zvib分別為轉(zhuǎn)動(dòng)和振動(dòng)松弛碰撞數(shù),Pr為氣體Prandtl數(shù),μ為黏性系數(shù).在DSMC方法中,通常取但是在動(dòng)理學(xué)模型方程中,需要通過溫度的松弛過程來確定其中的松弛碰撞數(shù).在實(shí)際計(jì)算中通常取

由于分布函數(shù)f(r,ξ,t;erot,evib)中包含了ξ,erot,evib,如果對(duì)這幾個(gè)變量均采用離散坐標(biāo)法,計(jì)算量將很大,為此引入如下三個(gè)約化速度分布函數(shù):

代入方程(6),有

其中,i=1,2,3,且有

則任意時(shí)刻三個(gè)約化速度分布函數(shù)對(duì)宏觀流動(dòng)量的表示式為:

為數(shù)值求解方便,引入如下無(wú)量綱參數(shù)將氣體動(dòng)力學(xué)模型方程進(jìn)行無(wú)量綱化:

則控制方程可寫為

宏觀流動(dòng)量的表達(dá)式可寫為:

2.2 守恒性及H定理的證明

在確定氣體動(dòng)理學(xué)模型方程后需要對(duì)其基本屬性進(jìn)行分析,這里考慮無(wú)量綱后的方程,并省略頂標(biāo) “?”. 首先設(shè)

這里,υrot=υtot/Zrot,υvib=υtot/Zvib. 則由質(zhì)量、動(dòng)量和動(dòng)能守恒關(guān)系,有

由于

則可知(17)式前兩項(xiàng)自動(dòng)滿足,同時(shí)可得

代入(17)式第三項(xiàng)得

則

與(5)式一致,說明(5)式的定義是合理的,滿足能量守恒.

由Boltzmann方程的H定理可知,氣體動(dòng)理學(xué)模型方程碰撞項(xiàng)Hi應(yīng)滿足

對(duì)于i=1的情形,假設(shè)υtr=υtot? υrot?υvib,(22)式變?yōu)?/p>

由于

而

由于υtr,υrot,υvib均大于零,當(dāng)f1,不同號(hào)時(shí)(26)式右邊第一項(xiàng)的被積函數(shù)小于零,其積分值也應(yīng)小于零,同樣當(dāng)f1與分別不同號(hào)時(shí)(26)式右邊第二項(xiàng)、第三項(xiàng)均小于零,當(dāng)且僅當(dāng)時(shí)(26)式等于零,即(23)式成立.

對(duì)于i=2和3的情形可類似得到.于是該動(dòng)理學(xué)模型方程的H定理得以證明.

2.3 物面邊界條件的可計(jì)算建模

在物面邊界上,氣體分子從物面反射進(jìn)入氣流,假定物面對(duì)分子沒有吸附作用,且反射是瞬時(shí)完成的,反射分子按完全漫反射處理,即以完全適應(yīng)于物面溫度Tw和速度Uw的平衡態(tài)分布散射,其分布函數(shù)為

約化后的速度分布函數(shù)為

其中,Cw=ξ?Uw,對(duì)于固定物面,Uw=0,nw為物面數(shù)密度.沿物面法線方向,流體質(zhì)量流量為零,有

可得出經(jīng)物面反射的氣體分子數(shù)密度nw,其中,Cnw=Cw·vwall,vwall為物面單位外法向矢量,由物面指向流體內(nèi)部.可由內(nèi)場(chǎng)分布函數(shù)插值得到.

3 構(gòu)造含振動(dòng)能激發(fā)的Boltzmann模型方程隱式格式

在確定氣體動(dòng)理學(xué)模型方程和物面邊界條件后,可采用離散速度坐標(biāo)法、LU-SGS隱式格式和有限體積法進(jìn)行求解[25?27].

以二維氣體流動(dòng)為例,經(jīng)離散速度坐標(biāo)法在離散速度點(diǎn)(ξxσ,ξyδ)處對(duì)速度空間離散降維的Boltzmann模型方程一般式可寫為

其中,σ,δ為離散速度網(wǎng)格點(diǎn),Sυσ,δ為控制方程右手項(xiàng).

在采用有限體積法求解上述方程時(shí),首先在網(wǎng)格中心型單元控制體積?IJ內(nèi)積分,有

1．PERT指征：臨床確診或疑診PEI，即可行PERT[1]?？筛鶕?jù)患者基礎(chǔ)疾病、PEI臨床癥狀、胰腺外分泌功能檢測(cè)、營(yíng)養(yǎng)不良的客觀證據(jù)等進(jìn)行綜合評(píng)估。

其中,F=(ξxσfσ,δ)i+(ξyδfσ,δ)j,v為單元控制體邊界沿位置空間網(wǎng)格I,J增大方向的法向量,符號(hào)表示量XIJ在控制體積?IJ內(nèi)的平均值.

根據(jù)積分變換原則,(31)式左邊第二項(xiàng)可以改寫為:

這里,A=?F?f=(ξx)i+(ξy)j.

利用Steger-Warming流通矢量分裂,將控制體積界面上的通量分解為正通量和負(fù)通量,其中正通量由計(jì)算得出,負(fù)通量由計(jì)算得出,即

最終控制方程可以寫成

若采用顯式格式,上式可由二階或三階具有TVD性質(zhì)的Runge-Kutta方法計(jì)算[23?27,30].為提高計(jì)算效率本文基于LU-SGS計(jì)算原理,構(gòu)造求解Boltzmann模型方程的隱式方法.根據(jù)隱式格式構(gòu)造原理,推導(dǎo)后有

由于網(wǎng)格中心型有限體積法所得到的分布函數(shù)位于控制體中心,因此物面分布函數(shù)應(yīng)首先通過靠近物面若干排控制體中心分布函數(shù)插值得到,再通過物面質(zhì)量流量為零和物面分子反射后按Maxwell分布的原則對(duì)物面分布函數(shù)進(jìn)行修正,進(jìn)而得到物面上各個(gè)宏觀繞流參數(shù).

4 算法驗(yàn)證與結(jié)果分析

為了考察直接求解Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法對(duì)稀薄流到連續(xù)流跨流域氣體流動(dòng)問題模擬能力與準(zhǔn)確收斂性,擬定來流馬赫數(shù)Ma∞=1.8、克努森數(shù)分別為Kn∞=0.3與Kn∞=0.0001的兩組圓柱繞流狀態(tài),圖1(a)繪出了本文統(tǒng)一算法對(duì)高稀薄圓柱繞流Kn∞=0.3流場(chǎng)沿駐點(diǎn)線密度分布計(jì)算結(jié)果與DSMC模擬值的比較情況.由圖可見,兩種結(jié)果符合較好,GKUA計(jì)算結(jié)果光滑穩(wěn)定,與DSMC模擬變化規(guī)律一致.圖1(b)展示了Kn∞=0.0001連續(xù)介質(zhì)流動(dòng)狀態(tài)下圓柱繞流流場(chǎng)物面壓力分布GKUA計(jì)算結(jié)果與連續(xù)介質(zhì)流理論的經(jīng)典分析數(shù)據(jù)比較情況,其中符號(hào)“□”表示GKUA計(jì)算結(jié)果,符號(hào)“Δ”為取自文獻(xiàn)[20]引用的連續(xù)流理論解經(jīng)典數(shù)據(jù),可看出二者整體上符合較好.這檢驗(yàn)了直接求解Boltzmann模型方程的GKUA在計(jì)算高稀薄流到連續(xù)流區(qū)繞流狀態(tài)正確性.只是在偏離駐點(diǎn)較遠(yuǎn)位置,存在一定的系統(tǒng)偏差,這可能是由于文獻(xiàn)[20]數(shù)據(jù)是由完全連續(xù)介質(zhì)流理論得到,而本文GKUA計(jì)算是以Kn∞=0.0001作為連續(xù)介質(zhì)流動(dòng)狀態(tài)參與比較,而真實(shí)連續(xù)流是克努森數(shù)趨于零所對(duì)應(yīng)的流動(dòng)狀態(tài).

為考察本文所建立的含振動(dòng)能激發(fā)Boltzmann模型方程及其在氣體動(dòng)理論統(tǒng)一算法中實(shí)現(xiàn)的正確可靠性,擬定非平衡效應(yīng)嚴(yán)重的稀薄過渡流區(qū)圓柱繞流問題進(jìn)行統(tǒng)一算法計(jì)算檢驗(yàn),并與DSMC方法結(jié)果對(duì)比分析.所有DSMC方法的結(jié)果均由Bird[52]的DS2V軟件模擬而得到.設(shè)置來流氣體為N2,Ma∞=5,T∞=Tw=500 K,Kn∞=0.01.圖2繪出了本文含振動(dòng)能激發(fā)的Boltzmann模型方程統(tǒng)一算法計(jì)算結(jié)果(圖中用“GKUA_vibration”表示)與DSMC方法結(jié)果(圖中用“DSMC”表示)的對(duì)比情況.從圖2中可以看出,兩者激波位置略有差別,這是由于DSMC方法采用的是自適應(yīng)的非結(jié)構(gòu)網(wǎng)格,在宏觀流動(dòng)量梯度大的位置進(jìn)行了自適應(yīng)加密,可以更好地捕捉激波,而本文GKUA采用的是近物面加密的結(jié)構(gòu)網(wǎng)格,在激波位置網(wǎng)格稍寬,對(duì)激波的捕捉能力稍弱,進(jìn)一步研究可通過加密網(wǎng)格得到更好的結(jié)果,但會(huì)增加計(jì)算開銷.從宏觀流動(dòng)量的分布云圖來看,兩種結(jié)果符合較好,壓力和等效溫度分布幾乎完全一致,表明本文建立的考慮振動(dòng)能激發(fā)的Boltzmann模型方程是合理可行的.從平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)溫度的分布來看,在波后駐點(diǎn)區(qū)域,平動(dòng)溫度最高,轉(zhuǎn)動(dòng)溫度次之,振動(dòng)溫度最低,本文計(jì)算所得的振動(dòng)溫度略高于DSMC方法的結(jié)果.

圖1 跨流域圓柱繞流流場(chǎng)與物面壓力分布GKUA計(jì)算與DSMC、連續(xù)流經(jīng)典數(shù)據(jù)比較驗(yàn)證 (a)Kn∞ =0.3時(shí)駐點(diǎn)線密度分布;(b)Kn∞=0.0001時(shí)沿物面迎風(fēng)面壓力分布Fig.1.Validation of fl ow fi eld and surface pressure distribution of Ma∞=1.8 around cylinder covering fl ow regimes with comparison of GKUA,DSMC and classical data of continuum fl ow:(a)Stagnation line density distribution for Kn∞=0.3;(b)wind-surface pressure distribution for Kn∞=0.0001.

對(duì)于多原子氣體,在不考慮離解等化學(xué)反應(yīng)的情況下,氣體分子除具有平動(dòng)能外,還具有內(nèi)能,即轉(zhuǎn)動(dòng)能和振動(dòng)能.氣體的轉(zhuǎn)動(dòng)特征溫度很低(一般小于10 K),一般情況下轉(zhuǎn)動(dòng)能完全激發(fā),而振動(dòng)特征溫度較高,振動(dòng)能隨溫度升高逐漸激發(fā).對(duì)于N2而言,其振動(dòng)特征溫度為3371 K,當(dāng)溫度大于750 K時(shí)其振動(dòng)自由度大于0.1,振動(dòng)能激發(fā)的影響開始顯現(xiàn).從圖2中各種溫度分布來看,波后絕大部分區(qū)域內(nèi)溫度都高于750 K,內(nèi)能激發(fā)將會(huì)對(duì)流場(chǎng)結(jié)構(gòu)產(chǎn)生影響,圖3(a)—圖3(c)給出了GKUA基于簡(jiǎn)單氣體的Shakhov模型計(jì)算結(jié)果(圖中用“GKUA_Shakhov”表示)與本文考慮振動(dòng)能激發(fā)的Boltzmann模型方程統(tǒng)一算法結(jié)果(圖中用“GKUA_vibration”表示)的對(duì)比情況.在考慮內(nèi)能激發(fā)后,激波更貼近物面,氣體分子一部分平動(dòng)能轉(zhuǎn)化成內(nèi)能,使得平動(dòng)溫度明顯降低,同時(shí)可以看出基于簡(jiǎn)單氣體Shakhov模型的統(tǒng)一算法結(jié)果的激波角略大于考慮內(nèi)能激發(fā)后的結(jié)果.圖3(d)—圖3(f)給出了僅考慮轉(zhuǎn)動(dòng)能激發(fā)的多原子氣體ES模型統(tǒng)一算法(圖中用“GKUA_ES”表示)與本文考慮振動(dòng)能激發(fā)的計(jì)算結(jié)果(圖中用“GKUA_vibration”表示)的對(duì)比情況.可以看出,兩種途徑得到的宏觀流動(dòng)量等值線云圖基本相同,但基于ES模型的GKUA計(jì)算所得的平動(dòng)溫度和轉(zhuǎn)動(dòng)溫度在波后要高于考慮振動(dòng)能激發(fā)的統(tǒng)一算法結(jié)果,說明氣體分子的部分能量轉(zhuǎn)化成振動(dòng)能使得溫度降低.圖3直觀清晰地展示了振動(dòng)激發(fā)對(duì)流場(chǎng)結(jié)果的影響.

圖2 (網(wǎng)刊彩色)本文含振動(dòng)能激發(fā)的Boltzmann模型方程GKUA計(jì)算所得宏觀流動(dòng)量與DSMC方法結(jié)果對(duì)比(a)馬赫數(shù);(b)壓力;(c)等效溫度;(d)平動(dòng)溫度;(e)轉(zhuǎn)動(dòng)溫度;(f)振動(dòng)溫度Fig.2.(color online)Macro-parameters simulated in present method comparing with DSMC method:(a)Mach number;(b)pressure;(c)overall temperature;(d)translational temperature;(e)rotational temperature;(f)vibrational temperature.

圖3 (網(wǎng)刊彩色)本文含振動(dòng)能激發(fā)的Boltzmann模型方程GKUA計(jì)算所得宏觀流動(dòng)量分別與簡(jiǎn)單氣體Shakhov模型、多原子氣體含轉(zhuǎn)動(dòng)能激發(fā)ES模型的GKUA結(jié)果對(duì)比 (a)馬赫數(shù);(b)壓力;(c)平動(dòng)溫度;(d)馬赫數(shù);(e)平動(dòng)溫度;(f)轉(zhuǎn)動(dòng)溫度Fig.3.(color online)Macroscopic fl ow variables simulated by the present GKUA with comparison of Shakhov model and ES model:(a)Mach number;(b)pressure;(c)translational temperature;(d)Mach number;(e)translational temperature;(f)rotational temperature.

圖4給出了分別考慮平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)不同非平衡效應(yīng),使用本文GKUA與DSMC計(jì)算所得的物面熱流與壓力分布曲線對(duì)比情況,其中“GKUA-Shakhov”表示將N2作為簡(jiǎn)單氣體采用基于Shakhov模型的GKUA計(jì)算結(jié)果,“GKUA-ES”表示僅考慮轉(zhuǎn)動(dòng)能激發(fā)而基于多原子氣體ES模型的GKUA計(jì)算結(jié)果,“GKUA-vibration”表示考慮N2振動(dòng)能激發(fā)的本文Boltzmann模型方程GKUA所得結(jié)果,“DSMC-translation”表示將N2作為簡(jiǎn)單氣體不考慮內(nèi)能激發(fā)的DSMC結(jié)果,“DSMC-rotation”表示僅考慮N2轉(zhuǎn)動(dòng)能影響的DSMC結(jié)果,“DSMC-vibration”表示考慮N2振動(dòng)能激發(fā)的DSMC結(jié)果.從圖4(a)所示物面熱流分布曲線可以看出,在駐點(diǎn)附近,“DSMC-translation”結(jié)果明顯高于其他幾種結(jié)果,且存在嚴(yán)重的統(tǒng)計(jì)波動(dòng),在繞流圓柱物面中部,“GKUA-vibration”偏離其他幾種結(jié)果,最大偏差不超過15%.而圖4(b)所示物面壓力分布曲線顯示在駐點(diǎn)附近“DSMC-translation”結(jié)果最小,“GKUA-Shakhov”次之,其他結(jié)果幾乎重合.表1給出了不同非平衡效應(yīng)內(nèi)能激發(fā)的Boltzmann模型方程統(tǒng)一算法與DSMC方法模擬得到的圓柱氣動(dòng)力系數(shù)對(duì)比情況,從中可以看出,相較于簡(jiǎn)單氣體模型,分別考慮轉(zhuǎn)動(dòng)、振動(dòng)能影響的Boltzmann模型方程GKUA所得軸向力系數(shù)與法向力系數(shù)均有所增大,但幅度在2.5%以內(nèi),說明內(nèi)能激發(fā)對(duì)氣動(dòng)力系數(shù)影響很小.

圖4 (網(wǎng)刊彩色)不同非平衡效應(yīng)不同方法模型計(jì)算得到的物面熱流和壓力分布 (a)熱流;(b)壓力Fig.4.(color online)Heat fl ux and pressure on the wall with different non-equilibrium e ff ects and method models:(a)Heat fl ux;(b)pressure.

表1 不同氣體模型GKUA算法與DSMC方法所得的軸向力系數(shù)和法向力系數(shù)對(duì)比Table 1.Axial force coefficients and normal force coefficients of cylinder by different gas models comparing GKUA method with DSMC method.

根據(jù)上述結(jié)果分析,由于物面溫度設(shè)置為500 K,較N2的振動(dòng)特征溫度3371 K低很多,不足以使與物面發(fā)生碰撞的氣體分子振動(dòng)能大量激發(fā),因此盡管動(dòng)理學(xué)模型方程中考慮了振動(dòng)能的激發(fā),但對(duì)與物面相關(guān)的氣動(dòng)力系數(shù)、物面熱流分布等影響很小,在流場(chǎng)內(nèi)部尤其是激波后的核心區(qū)域,溫度很高,足以引起振動(dòng)能的激發(fā),導(dǎo)致流場(chǎng)參數(shù)與流動(dòng)結(jié)構(gòu)發(fā)生明顯變化.

為了進(jìn)一步分析內(nèi)能激發(fā)、振動(dòng)特征溫度高低不同對(duì)繞流的影響,這里假設(shè)N2的振動(dòng)特征溫度為500 K,而其他條件不變,使用GKUA計(jì)算上述同一狀態(tài),得到的圓柱繞流軸向力系數(shù)為0.7043、法向力系數(shù)為0.4061.對(duì)比表1中振動(dòng)能一行,可以看出,即使把振動(dòng)特征溫度設(shè)置更低為500 K,氣體分子振動(dòng)能激發(fā)更為明顯時(shí),所計(jì)算的軸向、法向氣動(dòng)力系數(shù)僅由0.7044,0.4063降低為0.7043,0.4061,幾乎沒有變化,可以認(rèn)為在計(jì)算誤差范圍內(nèi),其變化幅度為0.049%.進(jìn)一步分析表1表明,分別使用GKUA與DSMC計(jì)算得到的結(jié)果偏差范圍0.06%—1.67%,證實(shí)兩種方法各自計(jì)算實(shí)現(xiàn)的正確性,同時(shí)可以看出從簡(jiǎn)單氣體到轉(zhuǎn)動(dòng)能再到振動(dòng)能影響,所帶來的軸向力、法向力系數(shù)相對(duì)變化量在1.97%—2.97%.圖5繪出了設(shè)置不同振動(dòng)特征溫度Θvib=500 K與Θvib=3371 K條件下,GKUA計(jì)算的流場(chǎng)中平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度、振動(dòng)溫度和等效溫度分布等值線云圖對(duì)比情況.由圖5可看出,各圖上半部分若設(shè)置振動(dòng)特征溫度Θvib=500 K較小時(shí),氣體分子的振動(dòng)能激發(fā)在較低溫度時(shí)發(fā)生,振動(dòng)非平衡效應(yīng)過于明顯,需要消耗更多氣體流動(dòng)能量,使得流場(chǎng)內(nèi)部的溫度偏低.尤其是在激波后的核心區(qū)域,溫度接近2500 K,以Θvib=500 K計(jì)則振動(dòng)自由度達(dá)到1.81,而實(shí)際Θvib=3371 K對(duì)應(yīng)的振動(dòng)自由度為0.946,前者接近完全激發(fā),從圖5(c)和圖5(d)可以看出,內(nèi)能完全激發(fā),分別在上述兩個(gè)振動(dòng)特征溫度設(shè)置下GKUA計(jì)算振動(dòng)能顯著激發(fā)后的繞流物體周圍的溫度分布呈現(xiàn)明顯變化.

圖5 (網(wǎng)刊彩色)不同振動(dòng)特征溫度設(shè)置下圓柱繞流流場(chǎng)溫度等值線云圖 (a)平動(dòng)溫度;(b)轉(zhuǎn)動(dòng)溫度;(c)振動(dòng)溫度;(d)等效溫度Fig.5.(color online)Temperature distribution in the fl ow fi eld with different sets of characteristic temperature of vibration:(a)Translational temperature;(b)rotational temperature;(c)vibrational temperature;(d)overall temperature.

圖6給出了上述兩種振動(dòng)特征溫度設(shè)置下圓柱繞流物面熱流和壓力分布.由于物面溫度保持為500 K,若設(shè)置Θvib=500 K所對(duì)應(yīng)的振動(dòng)自由度為1.164,而實(shí)際振動(dòng)特征溫度Θvib=3371 K時(shí)振動(dòng)自由度僅為0.016.從圖中可以看出,隨著振動(dòng)特征溫度的降低,振動(dòng)能的激發(fā)對(duì)駐點(diǎn)區(qū)域的熱流和壓力影響非常明顯,其中駐點(diǎn)熱流降低了約15%,壓力降低約8%.結(jié)合圖4對(duì)分別考慮平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)非平衡效應(yīng)所得到物面熱流與壓力分布的對(duì)比分析,可以推測(cè),當(dāng)物面溫度增加時(shí),即使振動(dòng)特征溫度較高,也會(huì)對(duì)物面的熱流和壓力分布產(chǎn)生影響.

圖6 (網(wǎng)刊彩色)不同振動(dòng)特征溫度設(shè)置下圓柱繞流物面熱流和壓力分布 (a)熱流;(b)壓力Fig.6.(color online)Heat fl ux and pressure along the wall surface with different sets of characteristic temperature of vibration:(a)Heat fl ux;(b)pressure.

5 結(jié) 論

根據(jù)表征分子微觀自由度之間能量傳遞的特征弛豫時(shí)間與流動(dòng)特征時(shí)間大小尺度的不同與氣體分子碰撞過程能量松弛演化特點(diǎn),將Boltzmann方程的碰撞項(xiàng)分解成彈性碰撞項(xiàng)和非彈性碰撞項(xiàng)兩部分,并將非彈性碰撞項(xiàng)按一定松弛速率分解為平動(dòng)-轉(zhuǎn)動(dòng)能松弛和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能松弛過程,將轉(zhuǎn)動(dòng)能、振動(dòng)能分別作為氣體分子速度分布函數(shù)的自變量,把轉(zhuǎn)動(dòng)能和振動(dòng)能處理為連續(xù)分布的能量模式,構(gòu)造了一類考慮振動(dòng)能激發(fā)的Boltzmann模型方程,并證明了其守恒性和H定理.通過對(duì)稀薄流到連續(xù)流跨流域圓柱繞流問題Boltzmann模型方程統(tǒng)一算法模擬驗(yàn)證與過渡流區(qū)考慮平動(dòng)、轉(zhuǎn)動(dòng)與振動(dòng)能激發(fā)的圓柱繞流模擬,得出以下結(jié)論.

1)通過將相同狀態(tài)下本文統(tǒng)一算法結(jié)果與DSMC方法結(jié)果對(duì)比,兩者符合較好,驗(yàn)證了本文建立的含振動(dòng)能激發(fā)的非平衡Boltzmann模型方程是準(zhǔn)確可靠的.

2)通過對(duì)簡(jiǎn)單氣體、僅考慮轉(zhuǎn)動(dòng)能激發(fā)的氣體以及考慮振動(dòng)能激發(fā)的氣體圓柱繞流進(jìn)行模擬,研究發(fā)現(xiàn)內(nèi)能激發(fā)對(duì)流場(chǎng)參數(shù)分布影響較大,激波位置更貼近物面,而溫度的分布變化最為顯著,激波后核心區(qū)的溫度因能量轉(zhuǎn)化使振動(dòng)能激發(fā)而有所下降.分子內(nèi)能激發(fā)熱力學(xué)非平衡效應(yīng)的實(shí)質(zhì)是平動(dòng)能向轉(zhuǎn)動(dòng)能、振動(dòng)能的傳遞以及各個(gè)能量間的交換,而根據(jù)經(jīng)典熱力學(xué)的原理,能量的宏觀表現(xiàn)為溫度,氣體分子各個(gè)能量間的傳遞達(dá)到平衡時(shí),必然會(huì)導(dǎo)致宏觀溫度的降低.

3)考慮分子內(nèi)能激發(fā)后,相比簡(jiǎn)單氣體,氣動(dòng)力系數(shù)會(huì)小幅增加,但即使氣體分子振動(dòng)能幾乎完全激發(fā),影響也較小.氣動(dòng)力的微觀本質(zhì)為氣體分子與物面碰撞引起的動(dòng)量變化,宏觀上為應(yīng)力張量在法向、切向作用的體現(xiàn),而在能量松弛過程中平動(dòng)能間的松弛時(shí)間比其他幾種要短得多,導(dǎo)致內(nèi)能的激發(fā)對(duì)應(yīng)力張量的變化影響較小,正應(yīng)力和切應(yīng)力小幅增加.

4)當(dāng)物面上氣體分子的振動(dòng)能激發(fā)顯著增加時(shí),駐點(diǎn)區(qū)域物面壓力和熱流會(huì)出現(xiàn)較為明顯的下降.壓力是通過平動(dòng)溫度計(jì)算,而熱流為溫度梯度的變化,對(duì)于等溫物面,在駐點(diǎn)附近溫度最高,內(nèi)能的激發(fā)更為顯著,平動(dòng)能降低導(dǎo)致平動(dòng)溫度降低,進(jìn)而壓力較小,內(nèi)能激發(fā)導(dǎo)致駐點(diǎn)區(qū)域等效溫度降低,溫度梯度也隨之降低,熱流也減小.

本文僅是含內(nèi)能激發(fā)熱力學(xué)非平衡效應(yīng)Boltzmann模型方程數(shù)值算法的初步工作與階段成果,下一步將開展更多狀態(tài)尤其是高馬赫數(shù)的模擬研究,并拓展到三維問題的計(jì)算分析,深入探索高超聲速條件下振動(dòng)能激發(fā)對(duì)飛行器氣動(dòng)力/熱特性的影響.

[1]Votta R,Schettino A,Bon fi glioli A 2013Aerosp.Sci.Technol.25 253

[2]Shevyrin A A,Vashchenkov P V,Bondar Y A,Ivanov M S 2014Proceedings of the 29th International Symposium on Rare fi ed Gas DynamicsXi’an,China,July 13–18,2014 p155

[3]Bird G A 1994Molecular Gas Dynamics and The Direct Simulation of Gas Flows(Oxford:Oxford University Press)pp50–54

[4]Shen Q 2003Rare fi ed Gas Dynamics(Beijing:National Defense Industry Press)pp38,83–88(in Chinese)[沈青 2003稀薄氣體動(dòng)力學(xué) (北京:國(guó)防工業(yè)出版社)第38,83—88頁(yè)]

[5]Struchtrup H 2005Macroscopic Transport Equations for Rare fi ed Gas Flows(Berlin:Springer)p27

[6]Chapman S,Cowling T G 1970The Mathematical Theory of Non-uniform Gases(Cambridge:Cambridge University Press)pp46–48

[7]Cercignani C 1988The Boltzmann Equation and Its Applications(New York:Springer Science Business Media)pp64–66

[8]Kremer G M 2010An Introduction to the Boltzmann Equation and Transport Processes in Gases(Berlin:Springer)p37

[9]Bhatnagar P L,Gross E P,Krook M 1954Phys.Rev.94 511

[10]Holway L H 1966Phys.Fluids9 1658

[11]Shakhov E M 1968Fluid Dynam.3 95

[12]Yang J Y,Huang J C 1995J.Comput.Phys.120 232

[13]Li Z H,Zhang H X 2002Acta Mech.Sin.34 145(in Chinese)[李志輝,張涵信 2002力學(xué)學(xué)報(bào) 34 145]

[14]Olga I R,Alexey P P,Irina A G 2013Comput.Fluids80 71

[15]Titarev V,Dumbser M,Utyuzhnikov S 2014J.Comput.Phys.256 17

[16]Li Z H,Peng A P,Fang F,Li S X,Zhang S Y 2015Acta Phys.Sin.64 224703(in Chinese)[李志輝,彭傲平,方方,李四新,張順玉2015物理學(xué)報(bào)64 224703]

[17]Xu K,Huang J C 2010J.Comput.Phys.229 7747

[18]Cai Z N,Li R 2014J.Comput.Phys.267 63

[19]Li Z H,Peng A P,Zhang H X,Yang J Y 2015Prog.Aerosp.Sci.74 81

[20]Li Z H 2001Ph.D.Dissertation(Mianyang:China Aerodynamics Research and Development Center)(in Chinese)[李志輝 2001博士學(xué)位論文 (綿陽(yáng):中國(guó)空氣動(dòng)力研究與發(fā)展中心)]

[21]Li Z H,Zhang H X 2004J.Comput.Phys.193 708

[22]Li Z H,Zhang H X 2005Adv.Mech.Sin.35 559(in Chinese)[李志輝,張涵信 2005力學(xué)進(jìn)展 35 559]

[23]Li Z H,Zhang H X 2007Acta Mech.Sin.:PRC23 121

[24]Li Z H,Zhang H X 2009J.Comput.Phys.228 1116

[25]Li Z H,Peng A P,Zhang H X,Deng X G 2011Sci.Sin.:Phys.Mech.Astron.54 1687

[26]Peng A P,Li Z H,Wu J L,Jiang X Y 2016Chin.J.Theor.Appl.Mech.48 95(in Chinese)[彭傲平,李志輝,吳俊林,蔣新宇2016力學(xué)學(xué)報(bào)48 95]

[27]Peng A P,Li Z H,Wu J L,Jiang X Y 2016J.Comput.Phys.327 919

[28]Li Z H,Wu J L,Jiang X Y,Ma Q 2015Acta Aeronaut.Astron.Sin.36 201(in Chinese)[李志輝,吳俊林,蔣新宇,馬強(qiáng)2015航空學(xué)報(bào)36 201]

[29]Li H Y 2007Ph.D.Dissertation(Mianyang:China Aerodynamics Research and Development Center)(in Chinese)[李海燕 2007博士學(xué)位論文 (綿陽(yáng):中國(guó)空氣動(dòng)力研究與發(fā)展中心)]

[30]Li H Y,Li Z H,Luo W Q,Li M 2014Sci.Sin.:Phys.Mech.Astron.44 194(in Chinese)[李海燕,李志輝,羅萬(wàn)清,李明2014中國(guó)科學(xué):物理學(xué)力學(xué)天文學(xué)44 194]

[31]Boyd I D,Josyula E 2011Phys.Fluids23 057101

[32]Yang H S 2013M.S.Thesis(Shanghai:Shanghai Jiaotong University)(in Chinese)[楊浩森2013碩士學(xué)位論文(上海:上海交通大學(xué))]

[33]Li Z,Zhu T,Levin D A 2013AIAA Paper AIAA2013-1201

[34]Wang C S,Uhlenbeck G E,Boer J D 1964Studies in Statistical Mechanics(Amsterdam:North-Holland Publishing Company)p2

[35]Wang C S(translated by Ying C T,Zhang C Z)1994The Kinetic Theory of a Gas(Beijing:Atom Energy Press)pp71–75(in Chinese)[王承書著 (應(yīng)純同,張存鎮(zhèn)譯)1994氣體分子運(yùn)動(dòng)論 (北京:原子出版社)第71—75頁(yè)]

[36]Morse T F 1964Phys.Fluids7 2012

[37]Andries P,Le Tallec P,Perlat J P,Perthame B 2000Eur.J.Mech.B:Fluid19 813

[38]Brull S,Schneider J 2009Continuum Mech.Thermodyn.20 489

[39]Rykov V A 1975Fluid Dynam.+10 959

[40]Rykov V A,Titarev V A,Shakhov E M 2008Fluid Dynam.+43 316

[41]Rykov V A,Titarev V A,Shakhov E M 2007Comp.Math.Math.Phys.+47 136

[42]Wu L,White C,Thomas J S,Reese J M,Zhang Y H 2015J.Fluid Mech.763 24

[43]Tantos C,Valougeorgis D,Frezzotti A 2015Int.J.Heat Mass Trans.88 636

[44]Tantos C,Ghiroldi G P,Valougeorgis D,Frezzotti A 2016Int.J.Heat Mass Trans.102 162

[45]Allu P,Mazumder S 2016Int.J.Heat Mass Tran.100 165

[46]Li Z H,Jiang X Y,Wu J L,Peng A P 2014Chin.J.Theor.Appl.Mech.46 336(in Chinese)[李志輝,蔣新宇,吳俊林,彭傲平2014力學(xué)學(xué)報(bào)46 336]

[47]Wu J L,Peng A P,Li Z H,Fang M 2015Acta Aerodynam.Sin.33 5(in Chinese)[吳俊林,彭傲平,李志輝,方明2015空氣動(dòng)力學(xué)學(xué)報(bào)33 5]

[48]Jiang X Y,Li Z H,Wu J L 2014Chin.J.Comput.Phys.31 403(in Chinese)[蔣新宇,李志輝,吳俊林2014計(jì)算物理31 403]

[49]Xu A G,Zhang G C,Li Y J,Li H 2014Prog.Phys.Sin.34 136(in Chinese)[許愛國(guó),張廣財(cái),李英駿,李華 2014物理學(xué)進(jìn)展34 136]

[50]Ying C T 1990Theory and Application of Gases Transport(Beijing:Tsinghua University Press)p62(in Chinese)[應(yīng)純同 1990氣體輸運(yùn)理論及應(yīng)用 (北京:清華大學(xué)出版社)第62頁(yè)]

[51]Laurent B,Bérénice G,Milana P C,Francesco S 2013Proc.Appl.Math.Mech.13 353

[52]Bird G A 2005Proceedings of the 24th International Symposium on Rare fi ed Gas DynamicsMelville,Canada 2005 p541

Validation and analysis of gas-kinetic uni fi ed algorithm for solving Boltzmann model equation with vibrational energy excitation?

Peng Ao-Ping1)Li Zhi-Hui1)2)3)?Wu Jun-Lin1)Jiang Xin-Yu1)

1)(Hypervelocity Aerodynamics Institute,China Aerodynamics Research and Development Center,Mianyang 621000,China)
2)(State Key Laboratory of Aerodynamics,China Aerodynamics Research and Development Center,Mianyang 621000,China)
3)(National Laboratory for Computational Fluid Dynamics,Beijing 100191,China)

2 May 2017;revised manuscript

19 May 2017)

With the increase of temperature in fl ow fi eld,gas molecules possess not only rotational degree of freedom,but also vibrational energy excitation.In order to simulate and study the in fl uence of internal energy excitation on polyatomic gas fl ow with high temperature and high Mach number,according to the general Boltzmann equation,we consider the rotational and vibrational energy modes as the independent variables of gas molecular velocity distribution function.It is assumed that the rotational and vibrational energy modes are described by continuous distribution with degree of freedom and temperature.Based on the Borgnakke-Larsen collision model used in direct simulation Monte Carlo(DSMC)method,the collision term of Boltzmann equation with internal energy excitation is divided into elastic and inelastic collision terms.The inelastic collision is decomposed into translational-rotational energy relaxation and translational-rotationalvibrational energy relaxation according to a certain relaxation rate obtained from the reciprocalities of rotational and vibrational collisions numbers per one elastic collision.Then a kind of Boltzmann model equation considering the excitation of vibrational energy is constructed.For showing the consistency between the present model equation and Boltzmann equation,the conservation of summational invariants and the H-theorem of this model are proved.When solving the present model equation with numerical methods,because of the continuous energy modes,it is difficult to simulate this model equation directly.In this paper,three control equations are derived and solved by the LU-SGS(lower-upper symmetric Gauss-Seidel)method,and the cell-centered fi nite volume method with multi-block patched grid technique in physical space.As a result,these gas-kinetic uni fi ed algorithm(GKUA)with vibrational energy excitation has been developed.Results are presented for N2with different Knudsen numbers around cylinder from continuum to rare fi ed gas fl ow by using the present Boltzmann model equation,GKUA with simple gas model,and DSMC method.Very good agreement between the present model and DSMC results is obtained,which shows that the accuracy and reliability of the present model.Comparing the translational,rotational,vibrational,and total temperatures computed by different methods,the e ff ects of the rotational and vibrational degrees of freedom are demonstrated.For the simple gas model,the translational temperature is much higher than those for the other two models with internal energy excitation.At the same time,the distance from shock wave to wall for the simple gas model is about twice those for the other two models.On the other hand,the obtained aerodynamic force coefficients of the cylinder are increasing according to the sequence from the simple gas model to the rotational energy excitation model to the vibrational energy excitation model,but the variation range is very small.By reducing the gas characteristic vibrational temperature,the temperature after the shock wave is much lower,and the heat fl ux declines evidently at the stagnation point with the same temperature as the wall temperature.This implies that with the wall temperature increasing the heat fl ux declines.

Boltzmann equation,vibrational energy excitation,gas-kinetic uni fi ed algorithm,thermodynamics non-equilibrium e ff ect

(2017年5月2日收到;2017年5月19日收到修改稿)

10.7498/aps.66.204703

?國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(批準(zhǔn)號(hào):2014CB744100)和國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):11325212,91016027)資助的課題.

?通信作者.E-mail:zhli0097@x263.net

?2017中國(guó)物理學(xué)會(huì)Chinese Physical Society

http://wulixb.iphy.ac.cn

PACS:47.45.Ab,47.45.–n,47.85.Gj,34.50.EzDOI:10.7498/aps.66.204703

*Project supported by the National Basic Research Program of China(Grant No.2014CB744100)and the National Natural Science Foundation of China(Grant Nos.11325212,91016027).

?Corresponding author.E-mail:zhli0097@x263.net