李中華，黨雷寧，李志輝,2

1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所，綿陽 621000 2. 國家計算流體力學(xué)實(shí)驗(yàn)室，北京 100083

高超聲速飛行器在大氣中飛行，經(jīng)過激波壓縮的高溫氣體，激波后的溫度可達(dá)上萬度，此時飛行器前方的繞流流場，特別是弓形激波附近流場，氣體分子狀態(tài)遠(yuǎn)離平衡態(tài)，不僅氣體分子的移動溫度、轉(zhuǎn)動溫度和振動溫度存在著顯著差別，而且會伴隨著劇烈的化學(xué)反應(yīng)[1-4]。特別是在稀薄過渡流區(qū)，化學(xué)反應(yīng)具有強(qiáng)烈的非平衡特征，會對飛行器的力、熱特性產(chǎn)生顯著影響。隨著載人航天和臨近空間飛行器的發(fā)展，各類航天器在稀薄過渡流區(qū)飛行的時間越來越長，化學(xué)非平衡的影響尤為越突出[5-8]。

對于化學(xué)非平衡的處理，在數(shù)值模擬方法上，傳統(tǒng)的計算流體力學(xué)(Computational Fluid Dynamics，CFD)技術(shù)一般通過求解連續(xù)流(克努森數(shù)Kn≤0.01)的三維Navier-Stokes流體動力學(xué)方程,耦合化學(xué)反應(yīng)動力學(xué)模型對化學(xué)反應(yīng)過程求解，發(fā)展出了多種高效的數(shù)值計算方法[9-12]。隨著氣體稀薄程度的增大，連續(xù)介質(zhì)假設(shè)失效，基于Navier-Stokes方程的數(shù)值模擬方法會出現(xiàn)較大誤差。在稀薄流區(qū)(Kn≥1)，基于粒子隨機(jī)統(tǒng)計模擬的直接仿真蒙特卡羅(Direct Simulation Monte Carlo，DSMC)方法是處理各類非平衡流動的常用方法，得到了廣泛的發(fā)展和應(yīng)用[13-17]。隨著Kn的減小，DSMC計算效率會越來越低，直至計算難以進(jìn)行[18-20]。

處于中等Kn的過渡流區(qū)流動，無論在試驗(yàn)技術(shù)還是數(shù)值計算方面均是難于處理的一種流動[21-22]。為了研究過渡流區(qū)的數(shù)值模擬方法，眾多研究者進(jìn)行了各種嘗試?；贐oltzmann模型方程的統(tǒng)一算法，取得了較好的效果，但目前還沒有辦法處理化學(xué)非平衡的流動[23-27]?；贑FD技術(shù)和DSMC方法的Navier-Stokes/DSMC耦合計算方法，在理想氣體的耦合計算中，取得了很好的結(jié)果[28-30]，能夠處理包含移動、轉(zhuǎn)動、振動等非平衡過程，計算結(jié)果與全DSMC方法的結(jié)果很接近。在過渡流區(qū)，采用Navier-Stokes/DSMC耦合算法，能夠代替全DSMC方法模擬高超聲速非平衡流動[31-35]。

在理想氣體的耦合算法基礎(chǔ)上[36-40]，本文嘗試開展高超聲速化學(xué)非平衡流動的Navier-Stokes/DSMC耦合算法研究。基于包含化學(xué)非平衡流動的CFD和DSMC方法和軟件，建立了過渡流區(qū)化學(xué)非平衡流動的Navier-Stokes/DSMC耦合算法。利用耦合計算方法，對高超聲速二維圓柱繞流進(jìn)行了模擬，與相關(guān)文獻(xiàn)結(jié)果進(jìn)行比較，驗(yàn)證化學(xué)非平衡Navier-Stokes/DSMC耦合計算的有效性。

1 數(shù)值模擬方法

1.1 CFD技術(shù)

對高超聲速飛行器繞流問題，在連續(xù)流區(qū)域一般采用CFD算法求解Navier-Stokes方程組，可以得到準(zhǔn)確的飛行器氣動力特性，而且計算效率較高。

對于ns種組元的混合氣體，三維化學(xué)非平衡流動控制方程包含ns個組元質(zhì)量守恒方程，3個方向的動量守恒方程和1個能量守恒方程，具體形式為

(1)

式中：Q為守恒變量矢量；E、F、G為無黏通量矢量；Ev、Fv、Gv為黏性通量矢量；S為化學(xué)反應(yīng)源項；t為時間坐標(biāo)；x、y、z為3個方向坐標(biāo)。

具有ns個組元的化學(xué)反應(yīng)方程可以表示為

(2)

式中：αri、βri分別表示反應(yīng)物和生成物的當(dāng)量系數(shù)，下標(biāo)“r”和“i”分別代表第r個反應(yīng)以及第i個組元；Xi為第i個組元的分子式。

根據(jù)質(zhì)量反應(yīng)定律，式(2)中組元Xi的凈生成速率為

(3)

式中：[Xi]=ρi/Mi為組元Xi的摩爾濃度，ρi為組元密度，Mi為組元分子相對質(zhì)量；Rfr、Rbr為正向反應(yīng)速率和逆向反應(yīng)速率；kfr、kbr為正向反應(yīng)速率系數(shù)和逆向反應(yīng)速率系數(shù)。

本文使用Dunn&Kang的化學(xué)動力學(xué)模型。在Dunn&Kang的化學(xué)動力學(xué)模型中，正向、逆向反應(yīng)速率系數(shù)根據(jù)修正的Arrhenius經(jīng)驗(yàn)公式得到，即

(4)

式中：T為流場溫度；A為置前因子；B為溫度指數(shù)；C為反應(yīng)活化能與普適氣體常數(shù)的比值；下標(biāo)“f”和“b”分別表示正向、逆向反應(yīng)。

表1給出了5組元空氣反應(yīng)Dunn&Kang模型的正向和逆向反應(yīng)速率系數(shù)計算中的各常數(shù)值[5]。

采用有限體積法離散控制方程，采用無波動、無自由參數(shù)的耗散格式(Non-oscillatory, containing No free parameters and Dissipative scheme,NND)離散無黏項、中心差分格式離散黏性項。NND格式為

(5)

表1 5組元空氣反應(yīng)Dunn&Kang模型[5]Table 1 Dunn&Kang model for 5 species air reaction[5]

注：M1=N,NO;M2=O,O2,NO;M3=O2,N2;M4=O,N,NO

式中：W為單元界面重構(gòu)的流場變量，在本文中取控制方程的原始變量；上標(biāo)“L”和“R”分別表示變量是自單元界面的左邊(left)或右邊(right)重構(gòu)得到；minmod為限制器函數(shù)。

高超聲速化學(xué)非平衡流動的特征時間尺度和化學(xué)反應(yīng)特征時間尺度相差較大，各個化學(xué)反應(yīng)的特征時間尺度也有較大差異，由此帶來所謂的“剛性問題”，使得CFD求解收斂困難。目前主要通過3類方法解決“剛性”問題：全隱算法[9]、點(diǎn)隱算法[10]和解耦算法[11]，其中全隱算法將控制方程中的對流項、黏性項和化學(xué)反應(yīng)源項都進(jìn)行隱式處理，因此時間步長不受穩(wěn)定性條件限制。從并行計算的角度看，以LU-SGS(Lower-Upper Symmetric Gauss-Seidel)為代表的傳統(tǒng)隱式方法具有數(shù)據(jù)依賴，不適合并行計算，在并行分區(qū)數(shù)較大時收斂速率下降。全矩陣DPLUR(Full Matrix Data-Parallel Lower-Upper Relaxation)[9]隱格式是一種全隱算法，它采用精確的分裂后的無黏Jacobian矩陣代替LU-SGS等對角化方法中的近似矩陣，避免了對角化近似對實(shí)際時間步長的影響；采用精確的源項Jacobian矩陣代替對角化方法中的近似矩陣，能夠反映不同化學(xué)反應(yīng)的快慢；采用內(nèi)迭代松弛過程代替LU-SGS中的掃描過程，消除了相鄰網(wǎng)格單元隱式求解中的數(shù)據(jù)依賴。作為一種收斂速率快、適合于并行計算的隱式方法，全矩陣DPLUR是美國高溫非平衡流動CFD軟件US3D主要的隱式方法，是美國另一高溫非平衡CFD軟件DPLR作為補(bǔ)充的隱式方法。

本文采用全矩陣DPLUR方法進(jìn)行時間推進(jìn)，計算過程為：在m=1,mmax中開展循環(huán)得到

(6)

1.2 DSMC方法

DSMC方法是求解稀薄氣體流動常用的數(shù)值模擬方法，采用大量的仿真分子，在微觀水平上模擬分子的運(yùn)動和碰撞過程，求解真實(shí)氣體流動過程。不論在模擬熱力學(xué)非平衡還是化學(xué)非平衡方面，DSMC方法均具有很強(qiáng)的能力。本文采用的是Bird的DSMC方法[13]，分子的模擬采用VHS(Variable Hard Sphere)分子模型，分子內(nèi)部能量交換模擬采用L-B(Larsen-Borgnakke)模型。

在處理化學(xué)反應(yīng)過程中，化學(xué)反應(yīng)模型采用1.1節(jié)CFD中相同的模型和數(shù)據(jù)?；瘜W(xué)反應(yīng)速率系數(shù)式(4)改寫為

kD(T)=ΛTηexp(-Ea/kBT)

(7)

式中：Ea為單個反應(yīng)所需的活化能；Λ和η分別與式(4)中的Afr和Bfr意義相同；kB為Boltzmann常數(shù)。

高溫空氣中，一個組元A分子與一個組元B分子發(fā)生碰撞，其發(fā)生離解或置換反應(yīng)的概率為

(8)

復(fù)合反應(yīng)概率與第三體的數(shù)密度nT成正比。通常情況下，如果活化能為零，復(fù)合反應(yīng)的概率為

(9)

1.3 MPC耦合方法

在高超聲速流動中，雖然采用CFD技術(shù)，能夠達(dá)到較高的計算效率。但是非平衡流動中，往往會出現(xiàn)Navier-Stokes方程失效的情況。一般采用當(dāng)?shù)豄n判斷Navier-Stokes方程的失效問題，其表達(dá)式為

(10)

式中：λ為當(dāng)?shù)亓鲌鰵怏w分子平均自由程；Q為當(dāng)?shù)亓鲌龊暧^參數(shù)，可以是壓力、密度或溫度。

該參數(shù)同時考慮了當(dāng)?shù)孛芏?平均自由程)和流場參數(shù)梯度的影響。取梯度最大的流動參數(shù)來計算Knl。當(dāng)Knl≥0.02，認(rèn)為連續(xù)流方程失效。在跨激波區(qū)域(特別是強(qiáng)激波)和邊界層內(nèi)，雖然密度較高，但流場梯度較大，容易出現(xiàn)連續(xù)流方程失效的情況。尾跡流區(qū)域內(nèi)，由于氣流膨脹迅速，而且密度較低，也容易出現(xiàn)稀薄氣體效應(yīng)(參見圖1)。

在連續(xù)流方程失效的區(qū)域，采用DSMC方法能夠很好地模擬稀薄氣體效應(yīng)。但是DSMC方法計算效率較低，限制了其應(yīng)用范圍。如果在一個流場中同時包含連續(xù)流區(qū)域和稀薄氣體效應(yīng)區(qū)域，采用兩種算法耦合計算，是一個有效的方法。

在一個流場中同時使用CFD和DSMC方法，需要采用一定的耦合方法，把兩套軟件耦合在一起。本文采用一種稱為MPC(Modular Particle-Continuum)的耦合處理技術(shù)[31-33]。

MPC的基本思想是把整個流場分為兩種區(qū)域，CFD區(qū)域和DSMC區(qū)域。兩種方法各自獨(dú)立計算。在分界面上，兩種方法進(jìn)行信息交換。

MPC耦合技術(shù)中，對計算區(qū)域劃分統(tǒng)一的網(wǎng)格。DSMC方法的計算區(qū)域需要根據(jù)計算過程中流場參數(shù)自動選取，選取那些連續(xù)流方程失效的區(qū)域。采用式(10)來判斷連續(xù)流方程是否失效。連續(xù)流方程失效的網(wǎng)格，組合成DSMC計算區(qū)域，在這些網(wǎng)格內(nèi)，采用DSMC方法進(jìn)行模擬。在計算過程中，根據(jù)CFD結(jié)果不斷調(diào)整DSMC方法計算區(qū)域，直到流場達(dá)到穩(wěn)定。圖1繪出按上述方法進(jìn)行自動調(diào)整的流場分區(qū)示意圖，其中藍(lán)色為Navier-Stokes計算區(qū)域，綠色為DSMC計算區(qū)域。

耦合計算過程中，兩種計算方法需要進(jìn)行雙向信息交換。CFD向DSMC區(qū)域信息傳遞時，在CFD與DSMC方法的分界網(wǎng)格面上，根據(jù)CFD計算的當(dāng)?shù)亓鲌鰠?shù)，向DSMC區(qū)域隨機(jī)補(bǔ)充仿真分子，仿真分子攜帶耦合邊界的宏觀信息。

圖1 Navier-Stokes/DSMC耦合算法分區(qū)示意圖Fig.1 Sketch of regions in Navier-Stokes/DSMC hybrid algorithm

根據(jù)連續(xù)流失效判斷方法，當(dāng)Knl≥0.02時，在耦合邊界上，流動已經(jīng)開始出現(xiàn)非平衡現(xiàn)象，此時速度分布函數(shù)已經(jīng)不再符合Maxwell平衡分布，可以采用Chapman-Enskog分布來描述[41]。

對于偏離Maxwell分布的非平衡分布，其一階展開形式為

f(C)=f0(C)Γ(C)

(11)

式中：f0為平衡態(tài)Maxwell分布函數(shù)；C為無量綱化的分子熱運(yùn)動速度矢量，C=c/(2kBT/m′)1/2，c為 分子熱運(yùn)動速度，m′為氣體分子質(zhì)量。Γ(C)表達(dá)式為

(12)

其中：C為分子速度，下標(biāo)x、y、z為3個坐標(biāo)方向；q為熱流；τ為剪切應(yīng)力項；xi為i方向的坐標(biāo)；κ為熱傳導(dǎo)系數(shù)；μ為黏性系數(shù)；δij為Kronecker符號(當(dāng)i=j時，δij=1，當(dāng)i≠j，δij=0)；u、v、w為3個方向的流動速度；p為壓力。

本文采用文獻(xiàn)[41]給出的方法，在邊界面上給出符合Chapman-Enskog分布的分子速度。該方法基于Maxwell分布的速度采樣方法，根據(jù)當(dāng)?shù)氐姆瞧胶舛龋捎谩熬芙^-接受”的隨機(jī)步驟，產(chǎn)生符合Chapman-Enskog分布的分子速度。

DSMC向CFD區(qū)域進(jìn)行信息傳遞時，由于DSMC方法的結(jié)果有一定的統(tǒng)計波動，在化學(xué)非平衡的流場中，在某些區(qū)域，部分組元的含量很少，統(tǒng)計波動很大。DSMC統(tǒng)計波動如果傳遞到CFD求解區(qū)域，可能會導(dǎo)致CFD求解不穩(wěn)定。為了克服這一困難，本文采用了一種稱為“亞松弛”的技術(shù)來抑制DSMC方法的統(tǒng)計波動對CFD計算的影響。

(13)

式中：θ為網(wǎng)格中新參數(shù)比較小的權(quán)重；相應(yīng)地，1-θ為比較大的權(quán)重。

在本文的CFD計算中，采用的是單溫度模型，在DSMC結(jié)果向CFD傳遞信息時，把移動溫度、轉(zhuǎn)動溫度和振動溫度按自由度加權(quán)平均，作為一個溫度傳遞給CFD計算。

仿真實(shí)踐表明，這種“亞松弛”技術(shù)，能夠比較好地抑制DSMC方法的統(tǒng)計波動對CFD計算的影響。特別是化學(xué)非平衡流場中，含量較低的組元的統(tǒng)計波動，也能夠得到很好的抑制。

耦合方法的計算步驟為

步驟1計算區(qū)域劃分統(tǒng)一的網(wǎng)格。為了加快計算進(jìn)度，本文在邊界層內(nèi)預(yù)設(shè)nD層網(wǎng)格為DSMC計算區(qū)域(nD可調(diào))，其余為CFD計算區(qū)域。為各計算區(qū)域設(shè)置邊界條件和初始條件。

步驟2進(jìn)行CFD和DSMC獨(dú)立計算。為了盡量保持CFD和DSMC計算的同步，本文設(shè)置DSMC每計算mD步進(jìn)行1步CFD計算(mD可調(diào))。

步驟3信息交換。每進(jìn)行1步CFD計算，開展一次信息交換。耦合邊界上，DSMC根據(jù)CFD結(jié)果改變邊界信息，CFD根據(jù)DSMC結(jié)果采用亞松弛改變相應(yīng)網(wǎng)格點(diǎn)上的宏觀參數(shù)。

步驟4每隔kD個時間步長(kD可調(diào))調(diào)整一次DSMC計算區(qū)域。根據(jù)CFD結(jié)果，采用當(dāng)?shù)乜伺瓟?shù)判斷連續(xù)流方程的失效網(wǎng)格，標(biāo)注為DSMC計算的網(wǎng)格。

步驟5重復(fù)步驟2～步驟4，直到流場穩(wěn)定。

步驟6流場穩(wěn)定后，不再進(jìn)行分區(qū)調(diào)整，繼續(xù)重復(fù)步驟2～步驟3，直到DSMC統(tǒng)計到足夠的樣本。

2 算例結(jié)果與討論

2.1 二維圓柱高超聲速化學(xué)非平衡繞流

分別采用全DSMC和Navier-Stokes/DSMC耦合方法計算了空氣5組元(N2、O2、NO、N、O)的化學(xué)非平衡流場。計算模型為二維圓柱。計算條件與文獻(xiàn)[17]一致。來流參數(shù)：速度U∞=7 810 m/s， 溫度T∞=178.2 K，Kn∞=0.01，來流組元(N2、O2、O)的摩爾分?jǐn)?shù)分別為0.775、0.204、0.021，壁面溫度Tw=1 000 K。

本算例中，在采用相同的仿真粒子密度的條件下，均用24進(jìn)程并行計算，全DSMC和Navier-Stokes/DSMC計算時間之比約為2.4:1，表明耦合算法比全DSMC計算能提高計算效率。

圖2是采用全CFD和耦合算法時兩種殘差曲線的對比。采用全CFD計算時，殘差降低迅速，很快能下降4個量級，說明本文采用的方法能夠解決化學(xué)反應(yīng)的剛性問題。在耦合算法中，CFD的殘差降低較慢，有兩個原因：① 為了能夠盡量和DSMC計算保持同步，計算中采用了較小的CFL(Courant-Friedrichs-Levy)數(shù)(DSMC計算中采用的是真實(shí)的時間)，而且DSMC計算兩步，CFD計算一步；② 受DSMC統(tǒng)計波動影響，CFD殘差也有一定波動，而且在流場穩(wěn)定后，殘差與初始計算相比也下降3個量級后不再變化，說明耦合計算中CFD部分已經(jīng)收斂。

圖3是耦合算法和全DSMC計算得到的移動溫度(Ttr)場云圖。對比兩種結(jié)果，可以看出，在空間分布上，兩種結(jié)果基本一致。圖4是駐點(diǎn)線上N2移動溫度和振動溫度(Tvib)的比較。在駐點(diǎn)線上，耦合算法的結(jié)果中，由于引入了CFD的計算，移動溫度激波的起始位置稍微落后于全DSMC的結(jié)果，激波厚度更薄。振動溫度與全DSMC的結(jié)果基本一致，而且與文獻(xiàn)[17]結(jié)果也比較吻合。雖然在本文的CFD算法中，所采用的熱力學(xué)模型為熱力學(xué)平衡模型，但在熱力學(xué)非平衡的計算中，在誤差許可范圍內(nèi)，可以得到與全DSMC一致的結(jié)果。

圖2 CFD殘差曲線對比Fig.2 Comparison of CFD residual curves

圖5是駐點(diǎn)線上數(shù)密度(n)和壓力(p)分布。比較表明，兩種結(jié)果的數(shù)密度和壓力在駐點(diǎn)線上分布吻合較好，激波起始位置和波后參數(shù)基本一致。在駐點(diǎn)處，耦合算法的數(shù)密度比全DSMC結(jié)果稍高。

圖6是兩種算法結(jié)果中組元O摩爾分?jǐn)?shù)的分布云圖。在激波后，由于氣體溫度很高，O2的離解度較高。從圖7中給出的駐點(diǎn)線上各組元的摩爾分?jǐn)?shù)(fmole)分布圖中可以看出，波后O2的摩爾分?jǐn)?shù)很低，接近為0，說明O2基本上全部離解。N2的離解度也很高，3種結(jié)果中(文獻(xiàn)[17]、全DSMC、耦合算法)波后N2摩爾分?jǐn)?shù)的最小值分別為0.238、0.212、0.198，耦合算法的結(jié)果中，N2的離解度更高一些。相應(yīng)地，波后N的摩爾分?jǐn)?shù)耦合算法結(jié)果更高。比較表明，兩種算法給出了基本一致的結(jié)果，與文獻(xiàn)的結(jié)果吻合很好。

圖3 流場移動溫度分布Fig.3 Distribution of translational temperature in flow field

圖4 駐點(diǎn)線溫度分布比較(N2)Fig.4 Comparison of temperatures distribution on stagnation streamline (N2)

圖5 駐點(diǎn)線參數(shù)分布比較Fig.5 Comparison of parameter distribution on stagnation streamline

圖6 O摩爾分?jǐn)?shù)分布Fig.6 Distribution of O mole fraction

圖7 駐點(diǎn)線組元摩爾分?jǐn)?shù)分布比較Fig.7 Comparison of species mole fraction distribution on stagnation streamline

以上比較表明，在化學(xué)非平衡的模擬方面，耦合算法能夠得到與全DSMC模擬較為一致的結(jié)果。

圖8是圓柱表面壓力和熱流的比較。對于表面壓力p，全DSMC與耦合算法的結(jié)果基本重合在一起。對于表面熱流q，在駐點(diǎn)處，由于耦合算法的數(shù)密度稍高(參見圖5(a))，熱流也比全DSMC結(jié)果偏高，其他區(qū)域兩種結(jié)果符合很好。

圖8 壁面參數(shù)分布比較Fig.8 Comparison of surface parameters distribution

在整體氣動力、熱特性方面，對于圓柱的阻力系數(shù)CD，全DSMC和耦合算法的計算結(jié)果分別為：CD_DSMC=1.33，CD_Hybrid=1.34，比文獻(xiàn)[17]的結(jié)果CD_Ref=1.32分別高0.75%和1.51%。對于熱流系數(shù)Ch，全DSMC和耦合算法的結(jié)果分別為：Ch_DSMC=0.059，Ch_Hybrid=0.056，比文獻(xiàn)[17]的結(jié)果Ch_Ref=0.062分別低4.8%和9.6%。結(jié)果表明，耦合算法與全DSMC方法對整體氣動力、熱特性的模擬比較接近。

2.2 航天器碎片過渡區(qū)繞流耦合計算

采用所建立的耦合算法，對某航天器解體碎片在過渡區(qū)的化學(xué)非平衡繞流開展了計算。解體碎片簡化為球柱模型，模型長度為1.3 m，直徑為0.246 m。模擬參數(shù)見表2，迎角為0°。

圖9是不同高度上(H=80、70 km)流場壓力分布?？梢钥闯?，隨著高度降低，激波脫體距離逐漸減小。圖10是典型的表面壓力和熱流分布，在這種外形簡化條件下，最大壓力和熱流均出現(xiàn)在駐點(diǎn)區(qū)域。

圖11是駐點(diǎn)線上N、O的質(zhì)量分?jǐn)?shù)(fmass)分布。隨著高度降低，激波脫體距離減小，N2、O2開始離解的位置向后移動。對于O原子，其質(zhì)量分?jǐn)?shù)最大值差別不大。對于N原子，隨著高度的降低，質(zhì)量分?jǐn)?shù)增大，表明在計算范圍內(nèi)，高度降低，化學(xué)反應(yīng)越來越強(qiáng)烈。

圖12是軸向力系數(shù)Ca和駐點(diǎn)熱流q0隨高度的變化。隨著高度降低，Ca減小。在85～75 km 范圍內(nèi)，Ca降幅較大，在75～70 km范圍內(nèi)，Ca降幅較緩。隨著高度降低，q0大幅升高。

表2 碎片模擬高度和速度Table 2 Altitude and velocity of fragment simulation

圖9 碎片不同高度上的流場壓力分布Fig.9 Pressure distribution of flow field of fragment at different altitudes

圖10 碎片典型表面壓力和熱流分布Fig.10 Typical distribution of pressure and heating flux on fragment surface

圖11 駐點(diǎn)線N、O質(zhì)量分?jǐn)?shù)分布Fig.11 Distribution of mass fraction of N and O on stagnation streamline

圖12 軸向力系數(shù)和駐點(diǎn)熱流隨高度的變化Fig.12 Variation of Ca and q0 with altitude

3 結(jié) 論

1) 基于相同的化學(xué)反應(yīng)模型和數(shù)據(jù)，采用MPC耦合技術(shù)，發(fā)展了化學(xué)非平衡流動的Navier-Stokes/DSMC耦合算法，在近連續(xù)過渡流區(qū)高超聲速繞流模擬中實(shí)現(xiàn)了化學(xué)非平衡耦合計算。

2) 通過高超聲速圓柱繞流的算例與其他結(jié)果的比較研究，表明耦合算法不論在流場結(jié)構(gòu)、流場非平衡現(xiàn)象、還是飛行器表面參數(shù)、飛行器整體氣動力/熱特性方面，都能夠得到全DSMC計算吻合的結(jié)果，證實(shí)本文設(shè)計的耦合方法是可行的。

3) 通過對碎片在過渡流區(qū)的仿真，得到碎片在過渡流區(qū)的力/熱特性。在過渡區(qū)，軸向力系數(shù)隨高度下降而減小，駐點(diǎn)熱流隨高度降低而大幅增加。

本文的CFD算法中，所采用的熱力學(xué)模型為單溫度模型，該模型假設(shè)氣體處于熱力學(xué)平衡狀態(tài)。而DSMC方法是熱非平衡的模型，在耦合的過程中會出現(xiàn)一定的偏差。對于更精確的耦合算法，需要開展熱力學(xué)非平衡的CFD技術(shù)和相應(yīng)的耦合技術(shù)研究。