高超聲速化學非平衡流動Navier-Stokes/DSMC耦合算法

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

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

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

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

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

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

1 數(shù)值模擬方法

1.1 CFD技術

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

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

(1)

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

具有ns個組元的化學反應方程可以表示為

(2)

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

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

(3)

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

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

(4)

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

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

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

(5)

表1 5組元空氣反應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為單元界面重構的流場變量，在本文中取控制方程的原始變量；上標“L”和“R”分別表示變量是自單元界面的左邊(left)或右邊(right)重構得到；minmod為限制器函數(shù)。

高超聲速化學非平衡流動的特征時間尺度和化學反應特征時間尺度相差較大，各個化學反應的特征時間尺度也有較大差異，由此帶來所謂的“剛性問題”，使得CFD求解收斂困難。目前主要通過3類方法解決“剛性”問題：全隱算法[9]、點隱算法[10]和解耦算法[11]，其中全隱算法將控制方程中的對流項、黏性項和化學反應源項都進行隱式處理，因此時間步長不受穩(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等對角化方法中的近似矩陣，避免了對角化近似對實際時間步長的影響；采用精確的源項Jacobian矩陣代替對角化方法中的近似矩陣，能夠反映不同化學反應的快慢；采用內迭代松弛過程代替LU-SGS中的掃描過程，消除了相鄰網(wǎng)格單元隱式求解中的數(shù)據(jù)依賴。作為一種收斂速率快、適合于并行計算的隱式方法，全矩陣DPLUR是美國高溫非平衡流動CFD軟件US3D主要的隱式方法，是美國另一高溫非平衡CFD軟件DPLR作為補充的隱式方法。

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

(6)

1.2 DSMC方法

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

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

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

(7)

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

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

(8)

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

(9)

1.3 MPC耦合方法

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

(10)

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

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

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

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

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

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

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

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

根據(jù)連續(xù)流失效判斷方法，當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為無量綱化的分子熱運動速度矢量，C=c/(2kBT/m′)1/2，c為 分子熱運動速度，m′為氣體分子質量。Γ(C)表達式為

(12)

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

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

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

(13)

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

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

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

耦合方法的計算步驟為

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

步驟2進行CFD和DSMC獨立計算。為了盡量保持CFD和DSMC計算的同步，本文設置DSMC每計算mD步進行1步CFD計算(mD可調)。

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

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

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

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

2 算例結果與討論

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

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

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

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

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

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

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

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

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

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

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

圖6 O摩爾分數(shù)分布Fig.6 Distribution of O mole fraction

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

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

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

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

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

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

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

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

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

圖12是軸向力系數(shù)Ca和駐點熱流q0隨高度的變化。隨著高度降低，Ca減小。在85～75 km 范圍內，Ca降幅較大，在75～70 km范圍內，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 駐點線N、O質量分數(shù)分布Fig.11 Distribution of mass fraction of N and O on stagnation streamline

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

3 結 論

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

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

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

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