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

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

1)(中國空氣動力研究與發(fā)展中心,超高速空氣動力研究所,綿陽 621000)2)(中國空氣動力研究與發(fā)展中心,空氣動力學(xué)國家重點實驗室,綿陽 621000)3)(國家計算流體力學(xué)實驗室,北京 100191)

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

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

1)(中國空氣動力研究與發(fā)展中心,超高速空氣動力研究所,綿陽 621000)2)(中國空氣動力研究與發(fā)展中心,空氣動力學(xué)國家重點實驗室,綿陽 621000)3)(國家計算流體力學(xué)實驗室,北京 100191)

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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中,

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

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

代入方程(6),有

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

則任意時刻三個約化速度分布函數(shù)對宏觀流動量的表示式為:

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

則控制方程可寫為

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

2.2 守恒性及H定理的證明

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

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

由于

則可知(17)式前兩項自動滿足,同時可得

代入(17)式第三項得

則

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

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

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

由于

而

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

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

2.3 物面邊界條件的可計算建模

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

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

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

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

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

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

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

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

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

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

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

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

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

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

最終控制方程可以寫成

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

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

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

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

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

圖1 跨流域圓柱繞流流場與物面壓力分布GKUA計算與DSMC、連續(xù)流經(jīng)典數(shù)據(jù)比較驗證 (a)Kn∞ =0.3時駐點線密度分布;(b)Kn∞=0.0001時沿物面迎風(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.

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

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

圖4 (網(wǎng)刊彩色)不同非平衡效應(yīng)不同方法模型計算得到的物面熱流和壓力分布 (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ù)對比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的振動特征溫度3371 K低很多,不足以使與物面發(fā)生碰撞的氣體分子振動能大量激發(fā),因此盡管動理學(xué)模型方程中考慮了振動能的激發(fā),但對與物面相關(guān)的氣動力系數(shù)、物面熱流分布等影響很小,在流場內(nèi)部尤其是激波后的核心區(qū)域,溫度很高,足以引起振動能的激發(fā),導(dǎo)致流場參數(shù)與流動結(jié)構(gòu)發(fā)生明顯變化.

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

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

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

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

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

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

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

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

[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稀薄氣體動力學(xué) (北京:國防工業(yè)出版社)第38,83—88頁]

[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é)報 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é)報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é)位論文 (綿陽:中國空氣動力研究與發(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é)報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)[李志輝,吳俊林,蔣新宇,馬強2015航空學(xué)報36 201]

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

[30]Li H Y,Li Z H,Luo W Q,Li M 2014Sci.Sin.:Phys.Mech.Astron.44 194(in Chinese)[李海燕,李志輝,羅萬清,李明2014中國科學(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氣體分子運動論 (北京:原子出版社)第71—75頁]

[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é)報46 336]

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

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

[49]Xu A G,Zhang G C,Li Y J,Li H 2014Prog.Phys.Sin.34 136(in Chinese)[許愛國,張廣財,李英駿,李華 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īng)用 (北京:清華大學(xué)出版社)第62頁]

[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

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

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

?2017中國物理學(xué)會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