Maxwell氣固相互作用模型對(duì)稀薄高超聲速凹腔繞流流動(dòng)特征和熱環(huán)境的影響

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

1. 中國(guó)航天空氣動(dòng)力技術(shù)研究院,北京 100074 2. 清華大學(xué) 航天航空學(xué)院,北京 100084 3. 中國(guó)航天科技集團(tuán)有限公司 航天飛行器氣動(dòng)熱防護(hù)實(shí)驗(yàn)室,北京 100074

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

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

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

對(duì)于連續(xù)流區(qū)的高速凹腔繞流問(wèn)題,前人在實(shí)驗(yàn)和計(jì)算方面都做過(guò)大量的研究[3-4,7-11],Lawson和Barakos[12]在其綜述文章中做了比較全面的介紹。

早期的研究發(fā)現(xiàn)凹腔流動(dòng)特征與凹腔入口處的速度和凹腔幾何形狀息息相關(guān),幾何形狀主要通過(guò)凹腔長(zhǎng)深比(長(zhǎng)度與深度之比,L/D)表征[12]。一般而言,凹腔繞流分成3種模式：開(kāi)式、閉式和過(guò)渡式[9]。對(duì)于高速凹腔繞流,長(zhǎng)度較小而深度較大的凹腔通常為開(kāi)式凹腔,長(zhǎng)度較大而深度較小的凹腔為閉式凹腔。一般情況下,長(zhǎng)深比小于10的凹腔為開(kāi)式,長(zhǎng)深比大于13的凹腔為閉式,上述兩類流動(dòng)之間的中間狀態(tài)稱為過(guò)渡式凹腔流動(dòng),其長(zhǎng)深比為10～13[10]。需要注意的是,上述凹腔流動(dòng)模式的界限只是名義上的,僅僅具有指導(dǎo)意義,因?yàn)椴煌芯空邷y(cè)量得到的界限實(shí)際上是一個(gè)范圍,開(kāi)式凹腔長(zhǎng)深比的上邊界為9～11,閉式凹腔長(zhǎng)深比的下邊界為12～15[11]。

對(duì)于開(kāi)式凹腔繞流,外部主流在凹腔前緣分離,外部主流和凹腔之間形成的剪切層橫跨整個(gè)凹腔,并附著到凹腔后緣,在凹腔內(nèi)部形成一個(gè)單渦結(jié)構(gòu)[7]。對(duì)于閉式凹腔,外部主流在凹腔前緣分離后,沒(méi)有足夠的能量跨越整個(gè)凹腔,最終附著在凹腔底面；流動(dòng)再次在凹腔底面分離,并再次附著到凹腔后緣,在凹腔內(nèi)部形成了雙渦結(jié)構(gòu)[8]。過(guò)渡式凹腔兼具開(kāi)式和閉式凹腔的部分特點(diǎn),一般不太穩(wěn)定,在一定條件下容易向開(kāi)式或閉式凹腔轉(zhuǎn)化[4]。

關(guān)于凹腔熱環(huán)境方面的研究,邱波等[13]通過(guò)求解可壓縮的Navier-Stokes(N-S)方程,對(duì)橫縫流場(chǎng)結(jié)構(gòu)和熱流分布進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)縫隙壁面上形成了多個(gè)局部高熱流區(qū)；同時(shí),他們還研究了來(lái)流馬赫數(shù)、雷諾數(shù)和攻角等對(duì)防熱瓦橫縫旋渦結(jié)構(gòu)及熱環(huán)境的影響。

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

在現(xiàn)有的稀薄高超聲速凹腔繞流研究中,氣固相互作用(Gas-Surface Interaction,GSI)模型均為完全漫反射模型,國(guó)內(nèi)外的研究人員尚未考慮過(guò)GSI模型對(duì)凹腔繞流的影響。然而,大多數(shù)稀薄流問(wèn)題都涉及氣固相互作用,因?yàn)闅夤滔嗷プ饔媚Ｐ蜑橄”×鞯睦碚摲治龊蛿?shù)值模擬提供適當(dāng)?shù)倪吔鐥l件[20]。稀薄流中速度滑移和溫度跳躍大小就與氣固相互作用息息相關(guān)[21]。而且,由于凹腔的存在,流場(chǎng)的表面積與體積之比變大,氣固相互作用的影響更為明顯[22]。

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

1 氣固相互作用模型

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

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

cr=ci-2(n·ci)n

(1)

當(dāng)氣體分子發(fā)生漫反射時(shí),反射后的速度服從基于固體表面溫度Tw的Maxwell分布。在徑向?yàn)楸砻嫱夥ㄏ虻木植壳蜃鴺?biāo)系中,反射分子的速度為

(2)

式中：cr為反射速度cr的大??；cw為基于表面溫度Tw的最概然熱運(yùn)動(dòng)速率；Ri(1 ≤i≤ 4)為均勻分布于區(qū)間(0, 1)的隨機(jī)數(shù)；θ和φ分別為該球坐標(biāo)系中的仰角和方位角。

Maxwell氣固相互作用模型為完全漫反射和鏡面反射的加權(quán)平均,即σ部分的入射分子發(fā)生漫反射,剩下的1-σ部分反生鏡面反射,其中σ稱為Maxwell適應(yīng)系數(shù)[23]。

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

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

2 計(jì)算方法

2.1 DSMC方法

目前,稀薄氣體流動(dòng)最為廣泛使用的計(jì)算方法是由Bird[27]提出的DSMC方法,其前提是稀疏氣體假設(shè),即氣體分子之間的距離遠(yuǎn)大于其直徑?；谶@個(gè)假設(shè),氣體分子運(yùn)動(dòng)和碰撞實(shí)現(xiàn)解耦。已經(jīng)證明,當(dāng)仿真分子數(shù)趨于無(wú)窮大時(shí),DSMC收斂到稀薄氣體動(dòng)力學(xué)的Boltzmann方程[28]。

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

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

2.2 亞松弛采樣技術(shù)

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

首先,基于初始條件,實(shí)現(xiàn)宏觀量的初始化估計(jì),表示為A0。

(3)

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

(4)

3 物理模型

3.1 研究對(duì)象

考察帶有凹腔的三維平板模型,如圖1所示。凹腔長(zhǎng)度為L(zhǎng),深度為D,寬度為W,凹腔上游的平板長(zhǎng)度為L(zhǎng)u,下游長(zhǎng)度為L(zhǎng)d,平板寬度為Wp。坐標(biāo)原點(diǎn)位于平板凹腔幾何中心在平板表面的投影點(diǎn),x軸沿流向?yàn)檎?y軸垂直向上為正,z軸由右手法則確定。

需要注意的是,本文考慮的凹腔尺寸代表了典型的防熱瓦縫隙尺寸和航天器表面的缺陷尺寸,為了進(jìn)行驗(yàn)證,這些尺寸和文獻(xiàn)[14]中的一致。此外,為了深入研究GSI模型的影響,排除凹腔幾何因素的干擾,僅僅考慮長(zhǎng)深比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 來(lái)流及計(jì)算參數(shù)

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

表1 自由來(lái)流條件Table 1 Free-stream parameters

4 方法驗(yàn)證和確認(rèn)

準(zhǔn)確DSMC需要網(wǎng)格尺寸不大于當(dāng)?shù)胤肿悠骄杂沙?也要求時(shí)間步長(zhǎng)不大于當(dāng)?shù)胤肿悠骄鲎矔r(shí)間。此外,為了降低統(tǒng)計(jì)誤差,單個(gè)網(wǎng)格單元內(nèi)部的仿真分子數(shù)也應(yīng)足夠多。因此,設(shè)計(jì)一系列算例對(duì)研究方法進(jìn)行驗(yàn)證和確認(rèn),涉及所有算例的GSI模型均為完全漫反射模型,即Maxwell適應(yīng)系數(shù)σ= 1.00。

4.1 方法驗(yàn)證和網(wǎng)格尺寸

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

(5)

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

圖2 方法驗(yàn)證和網(wǎng)格無(wú)關(guān)性確認(rèn)Fig.2 Method validation and grid-independence verification

4.2 時(shí)間步長(zhǎng)

為進(jìn)行時(shí)間步長(zhǎng)無(wú)關(guān)性研究,采用3個(gè)時(shí)間步長(zhǎng)(Δt= 1.5×10-8s, 1.0×10-8s, 0.5×10-8s)計(jì)算,其中Δt= 1.0×10-8s為標(biāo)準(zhǔn)時(shí)間步長(zhǎng)。需要注意的是,上述3個(gè)時(shí)間步長(zhǎng)均遠(yuǎn)遠(yuǎn)小于基于來(lái)流條件的分子平均碰撞時(shí)間,也小于當(dāng)?shù)鼐W(wǎng)格尺寸與仿真分子速率的比值。

采用上述3個(gè)時(shí)間步長(zhǎng)計(jì)算出的凹腔下游壁面中心線上的表面壓力和傳熱系數(shù)分布如圖3所示。顯然,無(wú)論是表面壓力系數(shù)還是傳熱系數(shù),上述3個(gè)時(shí)間步長(zhǎng)給出的結(jié)果幾乎沒(méi)有變化,確認(rèn)了計(jì)算結(jié)果的時(shí)間步長(zhǎng)無(wú)關(guān)性。

圖3 時(shí)間步長(zhǎng)無(wú)關(guān)性確認(rèn)Fig.3 Verification for independence of time step

4.3 仿真分子數(shù)

為進(jìn)行仿真分子數(shù)無(wú)關(guān)性研究,采用3個(gè)仿真分子數(shù)(Particle Per Cell(PPC) = 40, 80, 120,PPC表示每個(gè)網(wǎng)格單元的平均仿真分子數(shù))計(jì)算,其中PPC = 80為標(biāo)準(zhǔn)仿真分子數(shù)。

圖4為采用上述3套網(wǎng)格計(jì)算出的凹腔下游壁面中心線上的表面壓力和傳熱系數(shù)分布。顯然,無(wú)論是表面壓力系數(shù)還是傳熱系數(shù),上述3個(gè)仿真分子數(shù)給出的結(jié)果幾乎都沒(méi)有變化,確認(rèn)了計(jì)算結(jié)果的仿真分子數(shù)無(wú)關(guān)性。

圖4 網(wǎng)格單元仿真分子數(shù)無(wú)關(guān)性確認(rèn)Fig.4 Verification for independence of average number of simulation particles in each cell

5 流場(chǎng)特征

圖5為不同Maxwell適應(yīng)系數(shù)給出的對(duì)稱面上凹腔內(nèi)部的無(wú)量綱密度(當(dāng)?shù)孛芏圈雅c來(lái)流密度之比,ρ*=ρ/ρ∞)等值線和流線圖。由圖5可見(jiàn),當(dāng)σ=1.00時(shí),即GSI為完全漫反射時(shí),外部流動(dòng)直接繞過(guò)凹腔,在凹腔前緣分離的剪切層再次附著在凹腔后緣,在凹腔內(nèi)部形成一個(gè)充滿腔體的單渦結(jié)構(gòu),分離點(diǎn)附近出現(xiàn)局部低密度區(qū),再附點(diǎn)附近則出現(xiàn)局部高密度區(qū)域。

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

圖5 對(duì)稱面凹腔內(nèi)部無(wú)量綱密度ρ/ρ∞等值線及流場(chǎng)結(jié)構(gòu)Fig.5 Density ratio contours ρ/ρ∞ with streamline traces inside cavity on symmetric plane

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

6 表面壓力和熱流

6.1 表面壓力

圖6是不同Maxwell適應(yīng)系數(shù)給出的3個(gè)凹腔表面中心線壓力的分布,圖中的pfp表示GSI為完全漫反射時(shí)對(duì)應(yīng)平板附著流動(dòng)(沒(méi)有凹腔的情形)的表面壓力,取值為12.48 Pa。

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

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

相對(duì)于平板壓力值pfp,凹腔底面(CF)的壓力都較低。隨著GSI從完全漫反射變化為純鏡面反射,由于外部流動(dòng)被卷入凹腔的程度減弱,渦結(jié)構(gòu)逐漸消失,底部演變?yōu)椤八浪畢^(qū)”,故壓力沿著CF的分布逐漸變得均勻。

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

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

6.2 表面熱流

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

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

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

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

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

7 結(jié) 論

采用DSMC方法模擬稀薄流區(qū)高超聲速凹腔繞流問(wèn)題,考慮GSI模型的影響,探索Maxwell適應(yīng)系數(shù)對(duì)凹腔流場(chǎng)特征和表面壓力、熱流的變化規(guī)律,得到結(jié)論如下：

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

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

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

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