高超聲速化學(xué)非平衡流動(dòng)Navier-Stokes/DSMC耦合算法

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

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

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

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

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

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

1 數(shù)值模擬方法

1.1 CFD技術(shù)

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

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

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

采用有限體積法離散控制方程，采用無波動(dòng)、無自由參數(shù)的耗散格式(Non-oscillatory, containing No free parameters and Dissipative scheme,NND)離散無黏項(xiàng)、中心差分格式離散黏性項(xiàng)。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)的流場(chǎng)變量，在本文中取控制方程的原始變量；上標(biāo)“L”和“R”分別表示變量是自單元界面的左邊(left)或右邊(right)重構(gòu)得到；minmod為限制器函數(shù)。

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

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

(6)

1.2 DSMC方法

DSMC方法是求解稀薄氣體流動(dòng)常用的數(shù)值模擬方法，采用大量的仿真分子，在微觀水平上模擬分子的運(yùn)動(dòng)和碰撞過程，求解真實(shí)氣體流動(dòng)過程。不論在模擬熱力學(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為單個(gè)反應(yīng)所需的活化能；Λ和η分別與式(4)中的Afr和Bfr意義相同；kB為Boltzmann常數(shù)。

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

(8)

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

(9)

1.3 MPC耦合方法

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

(10)

式中：λ為當(dāng)?shù)亓鲌?chǎng)氣體分子平均自由程；Q為當(dāng)?shù)亓鲌?chǎng)宏觀參數(shù)，可以是壓力、密度或溫度。

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

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

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

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

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

耦合計(jì)算過程中，兩種計(jì)算方法需要進(jìn)行雙向信息交換。CFD向DSMC區(qū)域信息傳遞時(shí)，在CFD與DSMC方法的分界網(wǎng)格面上，根據(jù)CFD計(jì)算的當(dāng)?shù)亓鲌?chǎng)參數(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時(shí)，在耦合邊界上，流動(dòng)已經(jīng)開始出現(xiàn)非平衡現(xiàn)象，此時(shí)速度分布函數(shù)已經(jīng)不再符合Maxwell平衡分布，可以采用Chapman-Enskog分布來描述[41]。

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

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

(11)

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

(12)

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

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

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

(13)

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

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

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

耦合方法的計(jì)算步驟為

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

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

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

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

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

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

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

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

分別采用全DSMC和Navier-Stokes/DSMC耦合方法計(jì)算了空氣5組元(N2、O2、NO、N、O)的化學(xué)非平衡流場(chǎng)。計(jì)算模型為二維圓柱。計(jì)算條件與文獻(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)程并行計(jì)算，全DSMC和Navier-Stokes/DSMC計(jì)算時(shí)間之比約為2.4:1，表明耦合算法比全DSMC計(jì)算能提高計(jì)算效率。

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

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

圖2 CFD殘差曲線對(duì)比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 流場(chǎng)移動(dòng)溫度分布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是圓柱表面壓力和熱流的比較。對(duì)于表面壓力p，全DSMC與耦合算法的結(jié)果基本重合在一起。對(duì)于表面熱流q，在駐點(diǎn)處，由于耦合算法的數(shù)密度稍高(參見圖5(a))，熱流也比全DSMC結(jié)果偏高，其他區(qū)域兩種結(jié)果符合很好。

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

在整體氣動(dòng)力、熱特性方面，對(duì)于圓柱的阻力系數(shù)CD，全DSMC和耦合算法的計(jì)算結(jié)果分別為：CD_DSMC=1.33，CD_Hybrid=1.34，比文獻(xiàn)[17]的結(jié)果CD_Ref=1.32分別高0.75%和1.51%。對(duì)于熱流系數(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方法對(duì)整體氣動(dòng)力、熱特性的模擬比較接近。

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

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

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

圖11是駐點(diǎn)線上N、O的質(zhì)量分?jǐn)?shù)(fmass)分布。隨著高度降低，激波脫體距離減小，N2、O2開始離解的位置向后移動(dòng)。對(duì)于O原子，其質(zhì)量分?jǐn)?shù)最大值差別不大。對(duì)于N原子，隨著高度的降低，質(zhì)量分?jǐn)?shù)增大，表明在計(jì)算范圍內(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 碎片不同高度上的流場(chǎng)壓力分布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é)非平衡流動(dòng)的Navier-Stokes/DSMC耦合算法，在近連續(xù)過渡流區(qū)高超聲速繞流模擬中實(shí)現(xiàn)了化學(xué)非平衡耦合計(jì)算。

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

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

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