金 浩，方 明，李埌全，劉 沙，鐘誠文

(1. 西北工業(yè)大學 航空學院，西安 710072；2. 中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引言

2020年8月SpaceX使用載人龍飛船完成NASA首次商業(yè)載人任務，民用航天邁入新的階段。商用航天的快速發(fā)展，以及降低飛行成本、保證飛行安全的迫切需求，對高超聲速飛行器氣動和防熱設計提出了更高要求。高超聲速飛行的典型特征是激波，導致飛行器周圍空氣強烈壓縮，產生大量熱量，轉化為內能后使得空氣溫度急劇升高，引起飛行器表面的嚴重氣動加熱現(xiàn)象[1-3]。為保證飛行安全，降低飛行器防熱系統(tǒng)體積和重量，需要對飛行器氣動熱環(huán)境進行準確預測。

高超聲速氣動熱研究方法主要有工程算法、數(shù)值模擬、風洞試驗以及飛行試驗[4-5]。工程算法通常只適用于簡單構型，且計算誤差較大；飛行試驗時間和經費成本高昂，通常作為地面氣動熱預測的驗證手段。在當前和可預見的未來，風洞試驗和數(shù)值模擬是氣動熱地面預測的主要手段。當前，氣動熱試驗[4]研究一般在激波風洞等高焓設備中進行，其能夠獲取較為復雜外形飛行器的大面積和局部氣動熱數(shù)據(jù)，但地面風洞條件下無法完全復現(xiàn)真實飛行環(huán)境。而且，地面氣動熱試驗時間較短，模型并未充分加熱，與飛行器的長時間飛行存在顯著差異，因而試驗獲得的熱流數(shù)據(jù)也與真實飛行情況存在較大差異。

隨著對高超聲速流動機理的深入認知和計算技術與資源的快速發(fā)展，數(shù)值模擬成為高超聲速氣動力/熱研究的重要手段。對于連續(xù)流區(qū)，基于N-S方程的計算體系已經較為成熟[6-7]。在本文重點關注的稀薄過渡流區(qū)，直接模擬蒙特卡羅[8-10]（Direct Simulation Monte Carlo，DSMC）方法是研究非平衡流動問題最為有效的手段之一，常用于高超聲速飛行器的氣動力/熱計算。在DSMC方法中，氣體分子與壁面碰撞后，以一定的方式[11]（如漫反射、CLL(Cercignani，Lampis and Lord)模型）發(fā)生反射并進入流場，而確定分子反射行為的一個重要因素就是壁面溫度。在工程實踐中，通常采用人為設定的恒溫壁面（冷壁或熱壁），而反射后的分子速度等重要信息則根據(jù)初始設定的恒溫壁面抽樣給出。然而，飛行器在臨近空間的長時間飛行，必然導致壁面溫度的顯著變化，任何初始設定的壁面恒溫都不會準確，甚至飛行器不同部位的溫度也會存在明顯差異，因而壁面溫度的恒溫處理必將導致氣動熱環(huán)境預測的較大偏差。

從物理上講，飛行器在受熱后，壁面溫度將逐漸升高，并將熱流傳導到飛行器內部，同時壁面將以黑體輻射的形式將部分能量輻射出去?？梢约僭O：當飛行器溫度達到平衡后，來自繞流氣體的氣動加熱，將與輻射釋放出去的能量達到平衡，即壁面輻射平衡。正是基于這一假設，國內外在壁面輻射平衡方面開展了大量研究工作。在連續(xù)流N-S方程方面，Edquist[12]基于LAURA軟件，發(fā)展了輻射平衡壁面溫度條件，研究了熱力學非平衡下的火星飛行器后體熱流分布；周欣欣[13]對壁面熱傳導下輻射平衡的溫度邊界條件進行了研究，發(fā)現(xiàn)壁面假設對數(shù)值計算熱流結果影響較大，采用輻射平衡壁面溫度邊界可以正確計算鈍頭體沿軌道駐點的溫度和熱流變化；董維中等[14]基于表面熱輻射效應、催化效應和燒蝕效應，使用AEROPH_Flow軟件，對球錐模型進行了計算，分析了壁面溫度分布對氣動熱環(huán)境的影響；原志超[15]采用等溫和輻射平衡壁面，計算了高超聲速圓柱完全催化下的表面熱流，驗證了其提出的“松弛”方法的有效性。目前，與DSMC方法相關的壁面輻射平衡研究工作處于起步階段，屈程[16-17]的工作具有重要意義，他發(fā)展了一種基于輻射平衡的壁面溫度邊界條件，并以簡單外形進行了驗證分析。然而，該方法初始分布時需要分布較多模擬分子，無法應用于真實飛行器氣動熱計算中。

受上述工作啟發(fā)，本文將在對傳統(tǒng)DSMC修正極小、計算量成本增量可控的約束條件下，發(fā)展DSMC方法壁面輻射平衡邊界條件，在此基礎上開發(fā)具有一定工程實用性的軸對稱構型DSMC求解器，以鈍錐構型驗證其正確性。重點考察雙錐構型試驗條件（冷壁）和輻射平衡條件下的分離區(qū)、壁面壓力、熱流和溫度分布，找尋對飛行器氣動和防熱設計有指導意義的普遍規(guī)律。

1 基于壁面輻射平衡的DSMC方法

1.1 DSMC氣動熱計算方法

DSMC方法計算氣動熱，需要先確定溫度邊界條件以及氣體分子-表面相互作用模型。目前，DSMC方法常采用恒溫壁面溫度邊界條件，在此邊界條件下，壁面溫度為恒定壁面溫度，由初始設定給出，氣體分子與物面作用反射進入流場的能量信息由該恒溫確定。

氣體分子-表面相互作用有多種模型，目前使用廣泛的包括鏡面反射模型、完全漫反射模型和Maxwell反射模型等[8-9]。對于鏡面反射模型，分子與物面碰撞后平行于物面的切向速度分量保持不變，而法向速度分量的大小不變，方向與碰撞前相反；對于漫反射模型，氣體分子與壁面碰撞后的速度與碰撞前的速度沒有關系，反射分子作為一個整體在速度上滿足Maxwell速度分布，即:

式中，c′為 分子的熱運動速度，k為Boltzmann常數(shù)，m為分子質量，Tr為反射溫度，基于完全漫反射邊界以及恒溫物面條件，Tr=Tw。根據(jù)分布函數(shù)取樣，反射后分子的速度如下:

對于Maxwell反射模型，物面反射由漫反射和鏡面反射組成，假設氣體分子發(fā)生漫反射和鏡面反射比例為 σ:(1?σ)，本文計算中σ取為0.85。

根據(jù)DSMC計算的基本原理，表面熱流密度qw定 義為:

式中，下標i、r分別代表入射分子和反射分子；E代表分子總能量，包括平動、轉動、振動三種能量；Δt、A分別為作用時間與作用面積。

1.2 輻射平衡壁面邊界

考慮表面熱輻射效應、催化效應、燒蝕效應和熱傳導效應，利用準定常假設，由交界面一直延伸到物體深處的能量守恒方程[14,18]為:

式中，qw為 表面熱流，qc為通過熱傳導從表面?zhèn)鞯奖诿嫔钐師崃鳎攀潜砻姘l(fā)射系數(shù)，若為絕對黑體，則 ε=1， σ為Stefan-Boltzmann常數(shù)，其值為5.6704×10?8Js?1m?2K?4，Tw為 壁面溫度，為表面燒蝕率，Hw和Ha分別為表面總焓和材料生成焓。

Hirschel[19]研究得出在來流速度小于8 km/s，飛行高度小于100 km時，對流換熱占主導地位，來自激波層的輻射熱流可忽略不計，本研究中qw為非催化條件下的對流熱流。加熱一段時間后，壁面處在完全“熱透”狀態(tài)，不再往壁面深處傳熱，表面熱傳導熱流qc= 0。因此對于非燒蝕表面，不考慮表面催化效應條件下，方程（4）可化簡為:

式（5）所表示的即為輻射平衡邊界條件。本文主要研究DSMC方法在該邊界條件下的熱流、壓力和溫度以及與恒溫物面邊界的計算差異。

1.3 輻射平衡壁面建模及實現(xiàn)

根據(jù)DSMC方法熱流計算流程，考慮完全將壁面“熱透”最終狀態(tài)，采用穩(wěn)定取樣的熱流值作為式（5）中qw值，經過幾次迭代后，得出最終平衡時表面溫度和熱流值，該邊界條件結合恒溫冷壁邊界條件，可以確定出實際飛行器表面熱流的范圍。具體實現(xiàn)步驟如下：

（1）根據(jù)來流條件進行初始化，將壁面初始溫度Tw,0設為恒溫冷壁溫度；

（2）由輻射平衡溫度邊界式（5），在壁面邊界編號為m的單元，表面的氣動加熱熱流與該壁面輻射出的熱流相等，流動經充分長的n1時間步進入穩(wěn)態(tài)，通過合理的n2時間步后抽樣得到流場及壁面熱流、壓力等信息。壁面溫度由:

計算得到，并將該溫度作為下一次迭代的初值。

在本文研究中，假設物面為絕對黑體，即 ε=1；根據(jù)來流條件和計算模型，確定流動穩(wěn)定時間步長n1和抽樣時間步長n2；若式（6）得到的壁面溫度低于Tmin， 將壁面溫度設為Tmin， 其中Tmin為迭代初始溫度和壁面溫度的最小值；收斂條件為迭代的壁面溫度兩次之間的相對誤差在1%以內。

結合上述步驟更新壁面溫度的邊界條件，基于輻射平衡的壁面熱流DSMC方法計算流程如圖1。

圖1 DSMC壁面輻射平衡計算流程Fig. 1 Flow chart of wall radiation equilibrium in DSMC

在此基礎上，課題組開發(fā)了一套基于流場笛卡爾網格及自適應技術的MPI并行軸對稱構型DSMC求解器（簡稱為DSMC2A）。

2 基于鈍錐構型的算例驗證

稀薄氣體壁面輻射平衡研究文獻較少，可以參考的是高超聲速鈍錐繞流[17]。鈍錐長1 m，頭部半徑35 mm，半錐角為5°，如圖2所示。來流條件對應飛行高度為100 km，速度7000 m/s，其他自由來流條件如表1所示。計算時，空氣取為N2和O2雙組分氣體，其摩爾百分比為0.763∶0.237，初始壁面溫度為350 K，采用輻射平衡溫度邊界和Maxwell反射模型。

表1 自由來流條件Table 1 Freestream conditions

圖2 鈍錐構型 ( 單位：mm)Fig. 2 Blunted cone mode (unit: mm)

圖3給出了輻射平衡條件下的物面熱流系數(shù)Ch，并和文獻[17]中結果進行了比較，可以看出兩者吻合較好，驗證了軸對稱程序和輻射平衡壁面建模和實現(xiàn)的正確性。

圖3 物面熱流系數(shù)驗證Fig. 3 Validation of Surface heat flux coefficient distribution

同時，還計算了恒溫冷壁（壁面溫度為350 K）完全漫反射條件下的駐點熱流值，為93.9 kW/m2，與輻射平衡下Maxwell模型的熱流值81.32 kW/m2進行了比較，兩者偏差達13.397%，表明壁面輻射平衡下的熱流值較冷壁情況顯著降低。

3 在雙錐構型中的應用

3.1 計算條件

采用CUBRC 48英寸激波風洞中進行試驗的25°/55°雙錐構型[20-22]，如圖4所示，試驗氣體為N2，來流馬赫數(shù)為15.6。

圖4 CUBRC尖雙錐模型 ( 單位：mm)Fig. 4 CUBRC sharp double-cone model (unit: mm)

計算迎角0°的試驗狀態(tài)，試驗來流條件如表2所示。分子碰撞采用VHS模型，考慮平動能、轉動能和振動能之間的能量交換，動能和內能能量交換采用Larsen-Borgnakke模型，使用無化學反應氣體模型，壁面采用完全漫反射模型，轉動弛豫碰撞次數(shù)取為5，振動弛豫碰撞次數(shù)[8-9]取為溫度T的表達式：

表2 雙錐構型計算的自由來流和表面條件[21]Table 2 Freestream and surface conditions for sharp double-cone model computation[21]

式中， ω為黏性指數(shù)，C1和C2為常數(shù)。對于氮氣，C1取為9.1，C2取為220.0。

3.2 試驗條件下的結果分析

激波風洞試驗時間短，模型壁面加熱并不嚴重，可以將初始壁面溫度下（即圖1中輻射平衡計算尚未開始）的計算視為試驗工況。圖5給出了流場分離區(qū)壓力云圖和流線，從中可以看到左、右分離點以及分離區(qū)范圍；與軟件SMILE和DS2V計算結果[21]進行比較，見表3，可以看出三種手段計算出的分離區(qū)范圍和大小，并無明顯差異。

圖5 壓力云圖和分離區(qū)流線(試驗條件)Fig. 5 Pressure contours and streamlines in separation region

表3 分離區(qū)計算結果比較Table 3 Comparison of calculated separation regions

圖6給出了模型壁面沿中軸線的壓力分布DSMC計算和CUBRC Run 7風洞試驗結果對比，可以看出壓力在錐前體、分離區(qū)、由于激波/激波相互作用產生的高壓區(qū)以及第二段錐上均與試驗值吻合得很好。

圖6 計算和試驗測量的壁面壓力驗證Fig. 6 Validation between calculated and measured surface pressure

圖7給出了DSMC方法預測的壁面中軸線熱流和風洞試驗結果，兩者在絕大部分測量點處表現(xiàn)出良好的一致性，DSMC計算的干擾區(qū)熱流峰值稍高于試驗測量值，但在可以接受的合理范圍內。

圖7 計算和試驗測量的壁面熱流驗證Fig. 7 Validation between calculated and measured surface heat flux

3.3 輻射平衡下的結果分析

圖8給出了輻射平衡（對應圖1中的迭代結束）條件下分離區(qū)與冷壁計算對比，可以看出輻射平衡條件下的左分離點為75.8 mm，相對于冷壁條件下前移；右分離點為103.9 mm，相對于冷壁條件下后移；分離區(qū)大小為28.1 mm，較冷壁條件下顯著增大。該圖還表明，相對于冷壁，輻射平衡條件下的分離區(qū)附近壓力下降。

圖8 冷壁和輻射平衡分離區(qū)的壓力云圖比較Fig. 8 Comparison of pressure contour in separation region between cold wall and radiative equilibrium condition

圖9給出了冷壁和輻射平衡表面壓力沿中軸線的計算結果，可以看出輻射平衡條件下的壓力峰值明顯低于冷壁條件壓力峰值，兩者差異約為15.4%，同時壓力峰值后移。在兩種邊界條件下，壓力在分離區(qū)和再附點處差異較大，但在錐前段和錐后凸臺段基本相同。

圖9 冷壁和輻射平衡壁面壓力比較Fig. 9 Comparison of surface pressure between cold wall and radiative equilibrium boundary condition

表4給出了冷壁和輻射平衡條件下的阻力系數(shù)（凸臺表面面積為參考面積），兩種壁面條件下的阻力系數(shù)差異約為0.33%，表明壁面溫度變化不會導致飛行器整體氣動力特性的顯著差異。

表4 阻力系數(shù)計算結果比較Table 4 Comparison of drag coefficient for different conditions

圖10給出了輻射平衡條件下的中軸線壁面溫度，是來流對壁面充分加熱，達到完全“熱透”的結果。可以看到，輻射平衡下的壁面溫度顯著高于冷壁溫度（297.2 K），但前緣和再附點處的溫度峰值均顯著低于來流總溫（2116 K）。對于整個模型而言，前緣位置溫度最高，再附點次之，這與直觀理解相吻合。這意味著，對于高超聲速飛行器，防熱系統(tǒng)設計的重點應該放在上述兩個區(qū)域。

圖10 輻射平衡壁面溫度Fig. 10 Wall temperatures for radiative equilibrium boundary condition

圖11給出了從冷壁迭代到輻射平衡的表面熱流結果，可以看出，受壁溫的影響，冷壁熱流與輻射平衡壁面熱流的差別非常大，后者遠小于前者，前緣處的熱流峰值比試驗狀態(tài)低約50%，再附點處的熱流峰值比試驗狀態(tài)低2/3左右。即使真實飛行中的壁面沒有達到輻射平衡，至少長時間的飛行也使得壁面充分加熱，即飛行條件下的熱流也應該顯著低于激波風洞給出的預測值。在實際應用中，可以認為真實熱流值在冷壁和輻射平衡壁面的計算結果之間，因此兩種條件相結合，可以給出真實飛行的熱流預測范圍估計。

圖11 冷壁和輻射平衡壁面熱流比較Fig. 11 Comparison of surface heat flux between cold wall and radiative equilibrium boundary condition

4 結論

本文在傳統(tǒng)CFD方法輻射平衡壁面模型基礎上，提出了適用于DSMC方法的輻射平衡壁面溫度條件，開發(fā)了軸對稱構型DSMC求解器，以鈍錐構型驗證了其正確性?；陔p錐構型，重點考察了恒溫冷壁（試驗條件）和輻射平衡條件下的壁面壓力分布和熱流。數(shù)值計算結果表明：恒溫冷壁條件得到的壁面壓力分布和熱流與風洞試驗結果吻合，兩種溫度邊界下的壓力峰值差異約為15.4%，但是整體氣動力特性差異僅約為0.33%。相對于冷壁，輻射平衡計算得到的前緣處熱流峰值降低約50%，再附點處的熱流峰值降低約三分之二；兩種條件相結合，可以給出壁面熱流的預測范圍。該項研究可應用于真實飛行器外形和飛行條件下的表面熱流和溫度分布預測，為臨近空間長航時飛行器設計提供可靠依據(jù)。本文研究并未考慮氣體輻射、壁面催化/氧化/燒蝕等其他復雜過程對壁面熱流大小的影響，后續(xù)將深入研究。