Maxwell氣固相互作用模型對稀薄高超聲速凹腔繞流流動特征和熱環境的影響

靳旭紅,黃飛,程曉麗,蘇鵬輝,*

1. 中國航天空氣動力技術研究院,北京 100074 2. 清華大學 航天航空學院,北京 100084 3. 中國航天科技集團有限公司 航天飛行器氣動熱防護實驗室,北京 100074

航天器再入過程之初,速度通常高達數十千米每秒。因此,在再入過程中,航天器與地球大氣產生劇烈的摩擦,一方面使其速度降低到著落速度,另一方面卻導致其外表面溫度高達1 700 K,這個溫度遠遠超過不銹鋼的熔點[1]。因此,為了保護航天器外部表面和內部設備,熱防護系統至關重要。比如,航天飛機表面一般都覆蓋了超過32 000片防熱瓦,這些防熱瓦通常以對齊或者交錯的方式相連,不可避免地產生了防熱瓦縫隙[2]。

此外,制造公差、傳感器安裝、非相似材料的不均勻膨脹或燒蝕等諸多原因也會導致航天器表面存在縫隙和缺陷等凹腔結構[3]。航天器表面凹腔結構的存在會影響其流動狀態和傳熱特性：首先,凹腔入口處產生邊界層分離和再附,導致局部熱流升高；其次,凹腔干擾可增加湍流度,促進邊界層轉捩,導致整個表面的熱流增加；最后,由于凹腔狹小,輻射散熱效應被阻塞,即使凹腔內熱流很低,也可能產生較高的凹腔表面溫度[4]。因此,表面凹腔的存在對局部熱防護系統會產生重要的影響,甚至導致局部防熱結構的破壞[5]。

然而,在熱防護系統熱載荷的計算和分析中,出于簡化問題的考慮,通常假設航天器具有光滑連續的表面。根據哥倫比亞航天飛機事故調查委員會(Columbia Accident Investigation Board,CAIB)的調查結果,造成哥倫比亞航天飛機墜毀事故最可能的原因是左部機翼前緣熱防護結構的破裂,致使機翼直接暴露于高速高能氣流中[6]。該調查結果不但得到了美國國家航空航天局事故調查組(NASA Accident Investigation Team,NAIT)的支持,同時還表明了再入條件下熱防護表面凹腔結構繞流的復雜性和對其展開研究的必要性。

對于連續流區的高速凹腔繞流問題,前人在實驗和計算方面都做過大量的研究[3-4,7-11],Lawson和Barakos[12]在其綜述文章中做了比較全面的介紹。

早期的研究發現凹腔流動特征與凹腔入口處的速度和凹腔幾何形狀息息相關,幾何形狀主要通過凹腔長深比(長度與深度之比,L/D)表征[12]。一般而言,凹腔繞流分成3種模式：開式、閉式和過渡式[9]。對于高速凹腔繞流,長度較小而深度較大的凹腔通常為開式凹腔,長度較大而深度較小的凹腔為閉式凹腔。一般情況下,長深比小于10的凹腔為開式,長深比大于13的凹腔為閉式,上述兩類流動之間的中間狀態稱為過渡式凹腔流動,其長深比為10～13[10]。需要注意的是,上述凹腔流動模式的界限只是名義上的,僅僅具有指導意義,因為不同研究者測量得到的界限實際上是一個范圍,開式凹腔長深比的上邊界為9～11,閉式凹腔長深比的下邊界為12～15[11]。

對于開式凹腔繞流,外部主流在凹腔前緣分離,外部主流和凹腔之間形成的剪切層橫跨整個凹腔,并附著到凹腔后緣,在凹腔內部形成一個單渦結構[7]。對于閉式凹腔,外部主流在凹腔前緣分離后,沒有足夠的能量跨越整個凹腔,最終附著在凹腔底面；流動再次在凹腔底面分離,并再次附著到凹腔后緣,在凹腔內部形成了雙渦結構[8]。過渡式凹腔兼具開式和閉式凹腔的部分特點,一般不太穩定,在一定條件下容易向開式或閉式凹腔轉化[4]。

關于凹腔熱環境方面的研究,邱波等[13]通過求解可壓縮的Navier-Stokes(N-S)方程,對橫縫流場結構和熱流分布進行了數值模擬,發現縫隙壁面上形成了多個局部高熱流區；同時,他們還研究了來流馬赫數、雷諾數和攻角等對防熱瓦橫縫旋渦結構及熱環境的影響。

然而,幾乎所有高超聲速凹腔繞流的研究都只考慮連續流區的層流和湍流,而涉及臨近空間稀薄流區凹腔繞流現象的研究較少[14]。考慮到航天飛機等可重復飛行器再入時在高空仍然面臨嚴峻的熱防護挑戰,Palharini等[15-16]采用直接模擬Monte Carlo(Direct Simulation Monte Carlo,DSMC)方法研究了80 km高度稀薄流區高超聲速凹腔繞流問題,探索凹腔長深比對凹腔內部流場結構、壁面壓力和熱流的影響規律；他們發現稀薄過渡流區幾何尺寸對高超聲速凹腔繞流的影響與連續流區存在不同,主要表現為當長深比等于4時凹腔內部就出現了雙渦結構和剪切層在凹腔底面再附,而連續流區產生類似現象對應的長深比為14。Paolicchi和Santos[17]也進行了類似的研究,只不過他們考慮的凹腔更深。靳旭紅等[18]采用DSMC方法對稀薄流區高超聲速凹腔繞流問題進行了數值模擬,研究了稀薄氣體效應和三維效應對凹腔內部流場結構和壁面熱流的影響規律。Guo和Luo[19]同樣采用DSMC方法研究了一系列凹腔繞流流場特征,他們的來流條件覆蓋了從近連續流到稀薄流的寬流域。

在現有的稀薄高超聲速凹腔繞流研究中,氣固相互作用(Gas-Surface Interaction,GSI)模型均為完全漫反射模型,國內外的研究人員尚未考慮過GSI模型對凹腔繞流的影響。然而,大多數稀薄流問題都涉及氣固相互作用,因為氣固相互作用模型為稀薄流的理論分析和數值模擬提供適當的邊界條件[20]。稀薄流中速度滑移和溫度跳躍大小就與氣固相互作用息息相關[21]。而且,由于凹腔的存在,流場的表面積與體積之比變大,氣固相互作用的影響更為明顯[22]。

本文在前人研究的基礎上,采用DSMC方法對稀薄高超聲速三維凹腔繞流問題進行模擬,考慮氣固相互作用模型對凹腔流動特征和氣動表面量(表面壓力和熱流)的影響,探索流動特征、表面壓力和熱流隨著GSI模型的變化規律,并從氣體動理論的角度闡明其作用機理,為航天飛機等可重復飛行器的熱防護系統精細化設計提供數據支持。

1 氣固相互作用模型

氣固相互作用的研究介于氣體動理論與固體物理之間,歷史上,第一個GSI模型由Maxwell提出,這就是著名的Maxwell氣固相互作用模型,也稱為Maxwell邊界條件[23]。Maxwell模型考慮了兩類相互作用：鏡面反射和漫反射。

當氣體分子發生鏡面反射時,其切向速度分量保持不變,而法向速度分量反向。因此,入射角等于反射角,反射前后氣體分子的動能不變。令入射分子的速度為ci,固體表面單位外法向量為n,則反射后的速度cr為

cr=ci-2(n·ci)n

(1)

當氣體分子發生漫反射時,反射后的速度服從基于固體表面溫度Tw的Maxwell分布。在徑向為表面外法向的局部球坐標系中,反射分子的速度為

(2)

式中：cr為反射速度cr的大小；cw為基于表面溫度Tw的最概然熱運動速率；Ri(1 ≤i≤ 4)為均勻分布于區間(0, 1)的隨機數；θ和φ分別為該球坐標系中的仰角和方位角。

Maxwell氣固相互作用模型為完全漫反射和鏡面反射的加權平均,即σ部分的入射分子發生漫反射,剩下的1-σ部分反生鏡面反射,其中σ稱為Maxwell適應系數[23]。

1971年,Cercignani和Lampis[24]提出一種新的GSI模型,并由Lord[25]在DSMC中實現,這就是所謂的Cercignani-Lampis-Lord(CLL)模型。CLL模型滿足互易性原理,通過引入對應的適應系數,考慮了氣體分子切向動量和法向能量在固體表面的適應程度,并且能復現高速稀薄流條件下的實驗結果。Padilla和Boyd[26]采用平板模型,通過對比速度分布、邊界層形狀和氣動表面量,對Maxwell模型和CLL模型做了評估,發現當氣體分子對表面的適應程度高于50%時,Maxwell模型和CLL模型給出的邊界層形狀和氣動表面量類似。

在Maxwell模型和CLL模型之后,雖然又有多種GSI模型相繼被提出,但稀薄高超聲速來流條件下氣體與固體表面相互作用中的基礎物理和化學過程仍然幾乎未知。在這種情況下,研究GSI模型對流動特征的影響,并評估氣動表面量對GSI模型參數的敏感性尤其具有指導意義。由于Maxwell模型在稀薄氣體動力學研究中應用最為廣泛,加之其具有簡潔的形式且滿足互易性原理,本文通過改變Maxwell適應系數σ,采用參數化研究的手段評估氣固相互作用對稀薄高超聲速凹腔繞流流動特征和氣動表面量的影響。

2 計算方法

2.1 DSMC方法

目前,稀薄氣體流動最為廣泛使用的計算方法是由Bird[27]提出的DSMC方法,其前提是稀疏氣體假設,即氣體分子之間的距離遠大于其直徑。基于這個假設,氣體分子運動和碰撞實現解耦。已經證明,當仿真分子數趨于無窮大時,DSMC收斂到稀薄氣體動力學的Boltzmann方程[28]。

DSMC方法通過模擬大量仿真分子實現對氣體流動的計算。每一個仿真分子都代表大量真實氣體分子,并具有位置、速度和內能等屬性。隨著模擬過程的進行,通過跟蹤仿真分子的軌跡、計算仿真分子之間的碰撞以及仿真分子和固體表面之間的相互作用,仿真分子的信息不斷更新[29]。跟蹤足夠多的時間步以保證流場收斂且統計誤差足夠小之后,統計出流場宏觀量(速度、溫度及壓力場等)和氣動表面量(表面壓力及熱流等)。

基于DSMC方法,課題組自主開發了一套通用的三維并行稀薄氣體流動分析軟件,其核心計算模塊為一套非結構網格DSMC程序,并已成功應用于大量工程問題。在該DSMC軟件中,氣體分子模型為變徑硬球(Variable Hard Sphere,VHS)模型[23],氣體分子之間的碰撞采用非時間計數(No Time Counter,NTC)采樣技術[23],氣體分子之間的內能交換采用Larsen-Borgnakke模型[30]。

2.2 亞松弛采樣技術

由于凹腔內部局部為低速流動,采用了孫泉華和Boyd[31]提出的宏觀量亞松弛采樣技術,通過引入松弛因子,等效地增加了采樣數量,有效地降低了DSMC方法模擬低速流動宏觀量計算的統計誤差。亞松弛采樣技術已被證明對局部為低速流動的問題(比如凹腔繞流和管道流動等問題)特別有用。定常流動的亞松弛宏觀量采樣技術通過如下過程實現[31]：

首先,基于初始條件,實現宏觀量的初始化估計,表示為A0。

(3)

最后,增加校正步以消除初始條件和較早時間步瞬時值引入的誤差,即

(4)

3 物理模型

3.1 研究對象

考察帶有凹腔的三維平板模型,如圖1所示。凹腔長度為L,深度為D,寬度為W,凹腔上游的平板長度為Lu,下游長度為Ld,平板寬度為Wp。坐標原點位于平板凹腔幾何中心在平板表面的投影點,x軸沿流向為正,y軸垂直向上為正,z軸由右手法則確定。

需要注意的是,本文考慮的凹腔尺寸代表了典型的防熱瓦縫隙尺寸和航天器表面的缺陷尺寸,為了進行驗證,這些尺寸和文獻[14]中的一致。此外,為了深入研究GSI模型的影響,排除凹腔幾何因素的干擾,僅僅考慮長深比L/D= 1的情形。為便于下文分析,把凹腔的上游壁面表示為FW,其位置xFW=-1.5 mm,凹腔底面表示為CF,其位置yCF=-3.0 mm,下游壁面AW,其位置為xAW= 1.5 mm。

圖1 帶有凹腔的平板示意圖Fig.1 Schematic drawing of flat plate with cavity

3.2 來流及計算參數

考慮80 km高度的大氣環境,來流大氣為體積比例為76.3%的N2和23.7%的O2,來流馬赫數Ma∞=26.77,攻角α= 0°,來流溫度T∞、壓力p∞和密度ρ∞采用美國標準大氣(1976)確定,來流速度u∞由來流馬赫數Ma∞確定,具體的來流參數列于表1中。平板及凹腔表面為完全漫反射表面,其溫度均保持為1 000 K,接近再入飛行器駐點附近的溫度。取凹腔長度為特征長度,則對應的來流Reynolds數和Knudsen數分別為Re∞=126.875和Kn∞=0.367,表明該凹腔繞流問題是一個稀薄過渡流區的高超聲速層流問題。

表1 自由來流條件Table 1 Free-stream parameters

4 方法驗證和確認

準確DSMC需要網格尺寸不大于當地分子平均自由程,也要求時間步長不大于當地分子平均碰撞時間。此外,為了降低統計誤差,單個網格單元內部的仿真分子數也應足夠多。因此,設計一系列算例對研究方法進行驗證和確認,涉及所有算例的GSI模型均為完全漫反射模型,即Maxwell適應系數σ= 1.00。

4.1 方法驗證和網格尺寸

為進行網格無關性研究,設計3套網格(稀疏、標準和稠密網格),網格單元數分別為544 040、4 330 480和14 640 000。每套網格均為非均勻網格,在平板和凹腔附近實施了局部加密。令表面壓力和熱流為p和q,則表面壓力和傳熱系數的定義為

(5)

圖2為采用上述3套網格計算出的凹腔下游壁面(AW)中心線上的表面壓力和傳熱系數分布,同時還給出了Palharini等[14]的計算結果。顯然,無論是表面壓力系數還是傳熱系數,本文結果與Palharini等[14]的結果都比較一致,表明課題組發展的DSMC軟件能準確預測表面壓力和熱流等宏觀量,具備可靠性。此外,稀疏、標準和稠密網格給出的壓力和熱流幾乎都沒有變化,確認了計算結果的網格無關性。

圖2 方法驗證和網格無關性確認Fig.2 Method validation and grid-independence verification

4.2 時間步長

為進行時間步長無關性研究,采用3個時間步長(Δt= 1.5×10-8s, 1.0×10-8s, 0.5×10-8s)計算,其中Δt= 1.0×10-8s為標準時間步長。需要注意的是,上述3個時間步長均遠遠小于基于來流條件的分子平均碰撞時間,也小于當地網格尺寸與仿真分子速率的比值。

采用上述3個時間步長計算出的凹腔下游壁面中心線上的表面壓力和傳熱系數分布如圖3所示。顯然,無論是表面壓力系數還是傳熱系數,上述3個時間步長給出的結果幾乎沒有變化,確認了計算結果的時間步長無關性。

圖3 時間步長無關性確認Fig.3 Verification for independence of time step

4.3 仿真分子數

為進行仿真分子數無關性研究,采用3個仿真分子數(Particle Per Cell(PPC) = 40, 80, 120,PPC表示每個網格單元的平均仿真分子數)計算,其中PPC = 80為標準仿真分子數。

圖4為采用上述3套網格計算出的凹腔下游壁面中心線上的表面壓力和傳熱系數分布。顯然,無論是表面壓力系數還是傳熱系數,上述3個仿真分子數給出的結果幾乎都沒有變化,確認了計算結果的仿真分子數無關性。

圖4 網格單元仿真分子數無關性確認Fig.4 Verification for independence of average number of simulation particles in each cell

5 流場特征

圖5為不同Maxwell適應系數給出的對稱面上凹腔內部的無量綱密度(當地密度ρ與來流密度之比,ρ*=ρ/ρ∞)等值線和流線圖。由圖5可見,當σ=1.00時,即GSI為完全漫反射時,外部流動直接繞過凹腔,在凹腔前緣分離的剪切層再次附著在凹腔后緣,在凹腔內部形成一個充滿腔體的單渦結構,分離點附近出現局部低密度區,再附點附近則出現局部高密度區域。

但是,隨著σ的減少(漫反射的比例逐漸降低),凹腔內部流場特征和密度分布出現了一些變化。比如,渦結構不再充滿整個凹腔,而是向凹腔左上角移動,渦形狀也隨之變化；當GSI為純鏡面反射時,渦結構完全消失,凹腔底部基本不存在流動,成為所謂的“死水區”。同時,凹腔后緣的再附現象逐漸消失,表現為后緣附近的高密度區逐漸縮小甚至消失。此外,凹腔內部的密度也逐漸減小：GSI模型為完全漫反射時,凹腔再附點附近區域的密度均大于來流密度,其他區域的密度則與來流密度相當；當σ減小到0.75時,再附點高密度區域面積減小的同時,峰值密度亦降低；當σ進一步減小到0.50(即一半漫反射一般鏡面反射)時,在再附點高密度區面積進一步縮小和峰值密度繼續減小的同時,凹腔內部大部分區域的密度均低于來流密度；σ=0.25時,再附現象及其引起的局部高密度區完全消失；σ=0,即GSI為純鏡面反射時,凹腔內部的密度呈現分層分布,底部出現了“死水區”。

圖5 對稱面凹腔內部無量綱密度ρ/ρ∞等值線及流場結構Fig.5 Density ratio contours ρ/ρ∞ with streamline traces inside cavity on symmetric plane

上述流場特征和密度分布隨著Maxwell適應系數σ的變化與氣體分子和凹腔表面之間的剪切作用息息相關。隨著σ的減小,氣體分子與凹腔表面之間的切向動量交換減弱,即黏性作用減弱,導致外部氣流被卷入凹腔內部的程度減弱。GSI為完全漫反射(σ=1.00)時,由于黏性剪切作用較強,在凹腔前緣分離的剪切層在后緣發生再附之后,一部分氣流沿著凹腔下游壁面被卷入凹腔內部,使得凹腔內部的密度較大的同時,還形成一個飽滿的渦結構。隨著σ的減小,GSI從完全漫反射向純鏡面反射變化,氣體分子與平板之間的黏性剪切作用減弱,凹腔頂部的剪切層逐漸消失,凹腔后緣的再附現象也隨之減弱甚至消失,外流流動與凹腔之間的相互干擾僅僅作用在凹腔頂部區域,導致凹腔內部的密度逐漸減小,渦結構也逐漸消失。

6 表面壓力和熱流

6.1 表面壓力

圖6是不同Maxwell適應系數給出的3個凹腔表面中心線壓力的分布,圖中的pfp表示GSI為完全漫反射時對應平板附著流動(沒有凹腔的情形)的表面壓力,取值為12.48 Pa。

圖6 不同GSI模型對應的凹腔表面中心線表面壓力分布Fig.6 Distributions of surface pressure along centerlines of cavity surfaces for different GSI models

沿著凹腔上游壁面(FW),不同Maxwell適應系數對應的壓力分布表現出相似性：壓力在FW頂部較低,沿著FW向下逐漸增大,并在y/D=-0.05 處達到峰值,此后一致單調降低直至FW底部。GSI模型對FW的壓力分布影響頗為明顯：GSI從完全漫反射變化為純鏡面反射時,FW表面壓力逐漸增大。具體而言,σ=1.00, 0.75, 0.50, 0.25, 0對應的FW壓力峰值與平板壓力的比值分別為0.48、1.69、5.37、13.27和34.29,非完全漫反射給出的FW壓力峰值均大于對應平板壓力,而且鏡面反射的峰值壓力是漫反射的70余倍,明顯強調了FW峰值壓力分布對于GSI模型的敏感性。

相對于平板壓力值pfp,凹腔底面(CF)的壓力都較低。隨著GSI從完全漫反射變化為純鏡面反射,由于外部流動被卷入凹腔的程度減弱,渦結構逐漸消失,底部演變為“死水區”,故壓力沿著CF的分布逐漸變得均勻。

不同Maxwell適應系數對應的凹腔下游壁面(AW)的壓力分布也呈現類似性：壓力在AW底部較低,此后沿著AW向上直至凹腔后緣一直單調增大。同樣,AW的壓力分布對GSI模型也特別敏感：隨著GSI從完全漫反射變化為純鏡面反射,AW表面壓力逐漸增大。具體而言,σ=1.00, 0.75, 0.50, 0.25, 0對應的AW壓力峰值與平板壓力的比值分別為12.25、19.49、28.15、39.00和56.37,均大于對應的平板壓力,鏡面反射的峰值壓力是漫反射的4倍多。

在飛行器氣動設計中,如果航天器或者防熱瓦表面比較光滑或飛行速度較快導致一部分氣體分子在其表面發生鏡面反射時,需要特別注意表面缺陷或縫隙上游壁面和下游壁面的壓力載荷。

6.2 表面熱流

圖7為不同Maxwell適應系數給出的3個凹腔表面中心線壓力的分布,圖中的qfp表示GSI為完全漫反射時對應平板附著流動(沒有凹腔的情形)的表面熱流,取值為79.54 kW/m2。由于GSI為純鏡面反射(σ= 0)時氣體分子與凹腔表面之間沒有能量交換,故熱流均為0。

圖7 不同GSI模型對應的凹腔表面中心線熱流分布Fig.7 Distributions of surface heat flux along centerlines of cavity surfaces for different GSI models

不同Maxwell適應系數對應的凹腔上游壁面(FW)熱流分布與壓力分布相似。除了σ=0.25之外,其他GSI模型對應的峰值熱流均小于平板熱流。GSI模型對FW的熱流分布影響頗為明顯：隨著σ從1.00降低到0.25,FW的表面熱流逐漸增大,鏡面反射的峰值熱流是漫反射的250余倍。因此,在飛行器氣動設計中,如果一部分氣體分子在其表面發生鏡面反射時,需要特別留意表面缺陷或縫隙上游壁面的熱載荷。

凹腔底面(CF)的熱流分布與壓力類似,隨著σ從1.00降低到0.25,由于黏性剪切作用減弱導致外流被卷入的程度減弱,熱流沿著CF的分布逐漸變得均勻。

不同Maxwell適應系數對應的凹腔下游壁面(AW)的熱流分布同樣呈現類似性：熱流在AW底部較低,此后沿著AW向上直至凹腔后緣一直單調增大。與GSI模型對AW表面壓力的影響不同,隨著σ從1.00降低到0.25,AW表面熱流逐漸降低。這是因為AW的表面熱流主要來源于外流與該表面之間的黏性剪切,隨著σ的降低,黏性剪切作用減弱,熱流自然降低。

7 結 論

采用DSMC方法模擬稀薄流區高超聲速凹腔繞流問題,考慮GSI模型的影響,探索Maxwell適應系數對凹腔流場特征和表面壓力、熱流的變化規律,得到結論如下：

1) 稀薄流條件(80 km)下,GSI為完全漫反射時,外部流動直接繞過凹腔,在凹腔前緣分離的剪切層再次附著在后緣,在凹腔內部形成一個充滿腔體的單渦結構。

2) 隨著GSI從完全漫反射向純鏡面反射變化,渦結構不斷減小直至消失,凹腔底部逐漸出現所謂的“死水區”；凹腔內部的密度也逐漸減小,凹腔后緣的再附現象逐漸減弱甚至消失。

3) 隨著GSI從完全漫反射向純鏡面反射變化,氣體與凹腔表面之間的切向動量交換減弱,即黏性剪切作用減弱,外部氣流被卷入凹腔的程度減弱,導致凹腔底面的壓力和熱流分布逐漸變得均勻。

4) 與GSI為完全漫反射相比,鏡面反射導致凹腔上游側面的峰值壓力增大70余倍、下游側面的峰值壓力增大4倍多,近鏡面反射(σ= 0.25)導致上游壁面的峰值熱流增大250余倍,表明在飛行器設計中,如果一部分氣體分子在其表面發生鏡面反射時,應特別留意表面缺陷或縫隙上游壁面壓力載荷、熱載荷和下游壁面的壓力載荷。