李志輝, 梁 杰, 李中華, 李海燕, 吳俊林, 戴金雯, 唐志共

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

0 引 言

可回收類(lèi)航天器如飛船返回艙，從空間軌道返回地球表面著陸場(chǎng)，先后經(jīng)歷在軌自由分子流，到進(jìn)入、下降段稀薄過(guò)渡流，以至高空近連續(xù)滑移流，最后到達(dá)近地面連續(xù)流，是一個(gè)跨流域多尺度非平衡變化過(guò)程[1-8]。如何準(zhǔn)確模擬返回艙以極高速度再入跨流域復(fù)雜多物理場(chǎng)非平衡繞流問(wèn)題？特別是準(zhǔn)確預(yù)測(cè)近連續(xù)滑移過(guò)渡流區(qū)高超聲速氣動(dòng)力/熱環(huán)境、姿態(tài)配平、升阻比、熱流變化等氣動(dòng)特征，對(duì)航天器精細(xì)化氣動(dòng)設(shè)計(jì)與成功回收著陸具有至關(guān)重要作用與指導(dǎo)意義[8-17]。載人飛船返回著陸過(guò)程，一旦失效可能造成船毀人亡，1967年前蘇聯(lián)“聯(lián)盟一號(hào)”載人飛船開(kāi)傘失敗，航天員以身殉職；2008年美國(guó)“獵戶號(hào)”飛船回收系統(tǒng)空投測(cè)試，飛船模型發(fā)生翻滾墜毀。以美國(guó)航天飛機(jī)為例，由地面風(fēng)洞模擬得到的機(jī)身體襟翼配平角為7.5°，而實(shí)際飛行試驗(yàn)結(jié)果為16°，導(dǎo)致這種極端不一致源于跨流域復(fù)雜多尺度流動(dòng)稀薄非平衡效應(yīng)，同樣2003年2月1日美國(guó)“哥倫比亞”號(hào)航天飛機(jī)返航再入到離地面60 km高空因左翼受損解體失事，直接與高空高溫稀薄非平衡流動(dòng)效應(yīng)相關(guān)聯(lián)[11-12,17-19]。對(duì)于具有大尺寸復(fù)雜結(jié)構(gòu)的跨大氣層航天飛行器，再入飛行走廊中大部分區(qū)域都為稀薄、過(guò)渡流區(qū)所覆蓋(200～40 km)，當(dāng)航天器返回地面時(shí)，通常會(huì)在太空邊緣的亞軌道狀態(tài)進(jìn)行推返分離、調(diào)姿配平，如飛船返回艙再入到130 km左右進(jìn)行推返分離，需要在飛行高度約125 km到105 km期間根據(jù)預(yù)裝的配平迎角進(jìn)行調(diào)姿配平，圖1所示飛船返回艙再入過(guò)程跨流域飛行示意圖；飛船返回艙一般采用橫偏質(zhì)心位置的方法來(lái)提供再入配平迎角和實(shí)現(xiàn)飛行軌跡機(jī)動(dòng)控制所需的配平升阻比，因此準(zhǔn)確預(yù)測(cè)配平迎角隨再入高度變化對(duì)控制系統(tǒng)以及返回艙落點(diǎn)精度至關(guān)重要[3-8,10-16]。

圖1 飛船返回艙再入跨流域高超聲速繞流過(guò)程示意圖Fig.1 Sketch on hypersonic flows of reentry spacecraft covering various flow regimes

正是由于航天器以極高超聲速再入跨流域高稀薄過(guò)渡到連續(xù)流區(qū)復(fù)雜多尺度真實(shí)氣體非平衡流動(dòng)特點(diǎn)，使用空氣動(dòng)力學(xué)研究賴(lài)以依靠的三種手段：理論計(jì)算、風(fēng)洞實(shí)驗(yàn)與飛行試驗(yàn)氣動(dòng)辨識(shí)，對(duì)飛船返回艙跨流域稀薄氣動(dòng)特性模擬揭示，特別是中間過(guò)渡流區(qū)高空配平迎角與實(shí)際飛行遙測(cè)結(jié)果差異較大。工程上通過(guò)控制系統(tǒng)余量設(shè)計(jì)來(lái)解決上述問(wèn)題，再入返回使用姿控發(fā)動(dòng)機(jī)進(jìn)行軌控保守設(shè)計(jì)，造成RCS系統(tǒng)燃料消耗與防熱結(jié)構(gòu)質(zhì)量大大增加，有效載荷受到限制。深空探測(cè)返回器與飛船返回艙有相似的大鈍頭體外形(如我國(guó)月地高速再入返回器)，以近第二宇宙速度、半彈道跳躍式再入大氣層[6,12-13,19-21]，涉及稀薄、過(guò)渡、連續(xù)等多個(gè)流態(tài)，經(jīng)歷化學(xué)平衡、熱化學(xué)電離輻射非平衡等跨尺度多物理場(chǎng)復(fù)雜流動(dòng)狀態(tài)。文獻(xiàn)[6,13,19-20]綜述了國(guó)內(nèi)外載人登月艙的研究發(fā)展概況，對(duì)我國(guó)發(fā)展載人登月艙涉及的關(guān)鍵技術(shù)進(jìn)行了分析。文獻(xiàn)[12, 21]對(duì)載人登月返回的再入方式選擇和再入彈道設(shè)計(jì)進(jìn)行討論，認(rèn)為必須考慮升阻比等問(wèn)題。文獻(xiàn)[3-4,7-8]分析了返回艙再入稀薄流域的配平特性和主要影響因素。文獻(xiàn)[22-23]采用完全氣體模型對(duì)幾類(lèi)再入返回器的氣動(dòng)特性進(jìn)行了對(duì)比分析，需要進(jìn)一步開(kāi)展真實(shí)氣體非平衡效應(yīng)對(duì)超高速再入返回艙跨流域氣動(dòng)特性的影響研究。文獻(xiàn)[15,24-25]采用三維可壓縮層流N-S方程的化學(xué)非平衡流動(dòng)計(jì)算方法，分析比較了再入返回艙(器)在完全氣體模型和化學(xué)非平衡/輻射加熱氣體模型下的流動(dòng)特征與氣動(dòng)特性，顯示出利用當(dāng)前空氣動(dòng)力學(xué)模擬手段對(duì)返回艙(器)以第一、二宇宙速度再入極高超聲速氣動(dòng)熱環(huán)境準(zhǔn)確預(yù)測(cè)還存在較大困難。

返回艙(器)以高超聲速再入大氣層時(shí)，氣流被前體強(qiáng)擾動(dòng)產(chǎn)生嚴(yán)重壓縮激波，波后溫度可能達(dá)10 000 K以上，飛行器周?chē)纬筛哽蕦樱殡S復(fù)雜的熱化學(xué)非平衡過(guò)程。黏性激波層、非平衡態(tài)、化學(xué)反應(yīng)及其對(duì)流動(dòng)的影響等都是再入過(guò)程跨流域空氣動(dòng)力學(xué)必須深入研究和妥善解決的關(guān)鍵基礎(chǔ)問(wèn)題，急待發(fā)展相關(guān)基礎(chǔ)理論、計(jì)算與實(shí)驗(yàn)?zāi)M方法。為此，國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目“航天飛行器跨流域空氣動(dòng)力學(xué)與飛行控制關(guān)鍵基礎(chǔ)問(wèn)題研究(2014CB744100)”頂層設(shè)計(jì)制定了核心前沿基礎(chǔ)研究方向一“返回艙跨流域非平衡氣動(dòng)力熱繞流與姿態(tài)配平模擬”[26-27]，確立了項(xiàng)目核心一子課題群圍繞我國(guó)回收類(lèi)航天器如飛船返回艙、月地高速再入返回器再入地球大氣層過(guò)程，建立氣動(dòng)力、氣動(dòng)熱、彈道計(jì)算分析方法與新的手段途徑，致力返回艙(器)再入跨流域非平衡氣動(dòng)力/熱繞流與姿態(tài)配平模擬，研究揭示返回艙再入跨流域真實(shí)氣體非平衡繞流現(xiàn)象規(guī)律。

項(xiàng)目執(zhí)行五年來(lái)，結(jié)合工程需求開(kāi)展前沿基礎(chǔ)研究與圍繞核心方向子課題群集約式研究實(shí)施、深化發(fā)展。針對(duì)項(xiàng)目核心任務(wù)研究方向一，回收類(lèi)航天器(返回艙)再入大氣層跨流域飛行過(guò)程，在初步建立返回艙(器)再入跨流域氣動(dòng)力/熱繞流模擬新的研究方法、手段途徑基礎(chǔ)上，搭建了求解Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法(GKUA)[28-32]、DSMC[21,33-34]、N-S/DSMC[35-36]、滑移N-S解算器[37]、低密度風(fēng)洞實(shí)驗(yàn)測(cè)試[38]等多種空氣動(dòng)力學(xué)模擬手段相互驗(yàn)證結(jié)合，返回艙再入從外層空間自由分子流到近地面連續(xù)流跨流域空氣動(dòng)力學(xué)一體化模擬平臺(tái)框架，已初步研制形成以彈道計(jì)算為主線，耦合氣動(dòng)特性計(jì)算模塊與動(dòng)態(tài)演示，結(jié)合數(shù)據(jù)庫(kù)技術(shù)與軟件工程，研制返回艙再入跨流域氣動(dòng)模擬軟件系統(tǒng)。本文擬在介紹考慮平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)能激發(fā)熱力學(xué)非平衡效應(yīng)的Boltzmann型速度分布函數(shù)方程可計(jì)算建模氣體動(dòng)理論統(tǒng)一算法取得重要進(jìn)展基礎(chǔ)上，分析驗(yàn)證GKUA與DSMC、連續(xù)/近連續(xù)流區(qū)(滑移)N-S解算器、N-S/DSMC耦合算法、低密度風(fēng)洞實(shí)驗(yàn)?zāi)M手段途徑的正確性，并簡(jiǎn)略介紹由此研發(fā)的返回艙(器)再入跨流域氣動(dòng)力/熱繞流模擬系統(tǒng)特點(diǎn)與工程適應(yīng)性。

1 Boltzmann方程可計(jì)算建模氣體動(dòng)理論統(tǒng)一算法

1.1 方法介紹

將宏觀流體力學(xué)與微觀分子動(dòng)力學(xué)連接起來(lái)的介觀Boltzmann速度分布函數(shù)方程[39-40]，如式(1)，通過(guò)描述氣體流動(dòng)過(guò)程分子速度分布函數(shù)基于位置空間、速度空間任一時(shí)刻由非平衡態(tài)向平衡態(tài)的演化，將各個(gè)流域不同尺度間分子輸運(yùn)現(xiàn)象統(tǒng)一起來(lái)。

如果得到分子速度分布函數(shù)，氣體動(dòng)力學(xué)中感興趣的各種宏觀物理量就可以通過(guò)對(duì)分布函數(shù)乘以分子速度某種形式的函數(shù)再對(duì)整個(gè)速度空間積分而確定。然而，求解該速度分布函數(shù)方程最大的困難在于其碰撞項(xiàng)高度非線性、高維積分微分屬性，屬于復(fù)雜多尺度剛性問(wèn)題，精確求解描述各流域氣體流動(dòng)特征的Boltzmann方程未成現(xiàn)實(shí)。隨著流體力學(xué)發(fā)展，國(guó)際上已有眾多學(xué)者基于質(zhì)量、動(dòng)量、能量守恒定律，由數(shù)學(xué)上較簡(jiǎn)單的統(tǒng)計(jì)碰撞模型代替Boltzmann方程碰撞項(xiàng)，提出許多表征原始Boltzmann方程碰撞松弛和統(tǒng)計(jì)物理特性的氣體動(dòng)理學(xué)理論(氣體分子運(yùn)動(dòng)論)模型方程[41-44]，特別是基于Boltzmann-BGK模型方程逐步發(fā)展起來(lái)的格子Boltzmann方法[45-46]、離散速度有限差分求解[47-50]、BGK型有限體積系列格式如GKS、UGKS、DUGKS[51-57]、線化Boltzmann方程[58-59]、矩分析法[60-63]等氣體動(dòng)理學(xué)理論研究途徑。為了探索跨流域氣體流動(dòng)問(wèn)題一體化模擬方法，過(guò)去近二十年來(lái)，我們從分析氣體分子速度分布函數(shù)與宏觀流動(dòng)量依賴(lài)關(guān)系出發(fā)，使用氣體分子碰撞松弛參數(shù)和當(dāng)?shù)仄胶鈶B(tài)分布函數(shù)，對(duì)Boltzmann方程碰撞積分可計(jì)算建模，先后建立起模擬各流域一維、二維、三維繞流問(wèn)題氣體動(dòng)理論(氣體運(yùn)動(dòng)論)統(tǒng)一算法(GKUA)[64-69]與復(fù)雜飛行器高超聲速氣動(dòng)力/熱問(wèn)題高性能并行計(jì)算應(yīng)用研究[28-32,68-70]。下面擬在此高效并行平臺(tái)基礎(chǔ)上，進(jìn)一步介紹含內(nèi)能激發(fā)熱力學(xué)非平衡效應(yīng)Boltzmann模型方程統(tǒng)一算法與再入各流域繞流問(wèn)題研究進(jìn)展。

針對(duì)航天飛行器再入各流域多尺度繞流過(guò)程存在嚴(yán)重的間斷粒子稀薄非平衡效應(yīng)，為了表征氣體黏性輸運(yùn)與熱傳導(dǎo)屬性，這里基于氣體動(dòng)理論Shakhov模型[43]，也可基于不同的氣體動(dòng)理學(xué)理論模型如Holway的橢球統(tǒng)計(jì)模型[30,42]或雙原子氣體Rykov模型[28-29,44]等，用以Maxwellian分布函數(shù)fM作為權(quán)函數(shù)的Hermite多項(xiàng)式的三階漸近展開(kāi)而得到氣體分子當(dāng)?shù)仄胶鈶B(tài)分布函數(shù)fN。另一方面，為了反映跨流域氣體分子相互作用、稀薄程度、分子模型與內(nèi)部能量變化關(guān)系，分析推導(dǎo)聯(lián)系于氣體溫度、密度、不同流區(qū)流態(tài)控制參數(shù)、分子黏性碰撞截面與擴(kuò)散碰撞截面及分子相互作用規(guī)則冪指數(shù)的氣體分子碰撞松弛參數(shù)ν統(tǒng)一模型。由此，確立可描述各流域全局馬赫數(shù)復(fù)雜流動(dòng)輸運(yùn)現(xiàn)象統(tǒng)一的Boltzmann模型速度分布函數(shù)方程(比原始Boltzmann方程要簡(jiǎn)單得多)，其無(wú)量綱形式為：

fN=fM·[1+(1-Pr)c·q(2c2/T-5)/(5PT/2)](3)

fM=n/(πT)3/2exp(-c2/T)(4)

Kn∞=λ∞/L(5)

這里，Pr是Prandtl數(shù),Pr=μCp/K，Cp是定壓比熱比，μ是描述氣體分子輸運(yùn)能力的黏性系數(shù)，K是熱傳導(dǎo)系數(shù)；n為數(shù)密度，P、T分別為壓力、溫度，q為熱流矢量，c為氣體分子熱運(yùn)動(dòng)速度，且c=c。λ∞是來(lái)流氣體分子平均自由程，Kn∞為基于特征長(zhǎng)度L的來(lái)流當(dāng)?shù)豄nudsen(克努森)數(shù)，是表征飛行器再入各流域流動(dòng)狀態(tài)控制參數(shù)；ω、α是聯(lián)系于氣體分子變徑硬球(VHS)與變徑軟球(VSS)模型指數(shù)[71]；χ是相關(guān)于氣體分子相互作用規(guī)則常數(shù)：χ=(ζ+3)/[2(ζ-1)]，ζ是聯(lián)系于分子間作用力大小F與分子間距離r的負(fù)冪律指數(shù)：F=1/rζ。

方程式(2)通過(guò)描述氣體流動(dòng)過(guò)程中分子速度分布函數(shù)f對(duì)位置空間、速度空間和時(shí)間的變化率關(guān)系，從分子運(yùn)動(dòng)論一級(jí)給出氣體的統(tǒng)計(jì)描述[39-40]，通過(guò)ν和fN將不同流域流態(tài)控制參數(shù)、宏觀流動(dòng)物理量、氣體黏性輸運(yùn)特性、熱力學(xué)效應(yīng)及氣體分子相互作用規(guī)則與分子模型用統(tǒng)一表達(dá)式聯(lián)系起來(lái)。氣體動(dòng)力學(xué)中感興趣的各種宏觀流動(dòng)量，如氣體密度、流動(dòng)速度、溫度、壓力、應(yīng)力張量、熱流矢量等，均可通過(guò)速度分布函數(shù)f乘以分子速度某種形式的函數(shù)對(duì)速度空間積分的辦法得到。

氣體分子速度分布函數(shù)實(shí)質(zhì)上是一種遵從氣體運(yùn)動(dòng)論概率統(tǒng)計(jì)規(guī)律的幾率密度分布函數(shù)，分布函數(shù)對(duì)速度分量Vx、Vy、Vz的函數(shù)依賴(lài)關(guān)系屬于指數(shù)型[73]，通過(guò)發(fā)展和應(yīng)用氣體動(dòng)理論離散速度坐標(biāo)法(DVOM)[30,48,65,68,74]，研制適于跨流域高超聲速繞流問(wèn)題模擬的離散速度坐標(biāo)點(diǎn)自適應(yīng)優(yōu)化選取技術(shù)和基于離散速度分布函數(shù)適應(yīng)大規(guī)模并行分布式求和的宏觀流動(dòng)量與物面氣動(dòng)力/熱動(dòng)態(tài)確定離散速度數(shù)值積分方法。可將單一的速度分布函數(shù)方程式(2)化為在各個(gè)離散速度坐標(biāo)點(diǎn)處基于時(shí)間、位置空間具有非線性源項(xiàng)的雙曲型守恒方程組，在計(jì)算位置空間(ξ,η,ζ)其守恒形式為：

根據(jù)氣體分子碰撞松馳與對(duì)流運(yùn)動(dòng)的非定常特點(diǎn)，將算子分裂有限差分方法拓展應(yīng)用于基于時(shí)間、位置空間與速度空間離散Boltzmann模型方程(11)數(shù)值求解，構(gòu)造直接捕捉分子速度分布函數(shù)演化更新的氣體動(dòng)理論耦合迭代數(shù)值格式：

其中，LS、Lξ、Lη、Lζ分別代表式(13)～式(16)差分算子：

這里，計(jì)算網(wǎng)格步長(zhǎng)只需由物理空間所要求計(jì)算精度決定，因所有宏觀流動(dòng)量通過(guò)離散速度分布函數(shù)積分求和動(dòng)態(tài)確定，使得物理空間網(wǎng)格劃分極為粗糙[28-32,64-69]，也能得到穩(wěn)定收斂解的算法優(yōu)勢(shì)；計(jì)算時(shí)間步長(zhǎng)Δt由所使用格式的差分穩(wěn)定條件確定：

其中，CFL為時(shí)間步長(zhǎng)調(diào)節(jié)系數(shù)，取CFL=0.99。

一旦各個(gè)離散速度坐標(biāo)點(diǎn)處氣體分子速度分布函數(shù)被數(shù)值求解，任一時(shí)刻物理空間各點(diǎn)的宏觀流動(dòng)量如氣體密度、流動(dòng)速度、溫度、應(yīng)力張量和熱流矢量等，需通過(guò)文獻(xiàn)[65,68,74-76]發(fā)展的相應(yīng)離散速度數(shù)值積分方法而被演化更新。由于一般大尺度高超聲速飛行器外形復(fù)雜，各部位幾何尺度差異很大，往返大氣層跨流域高超聲速飛行，會(huì)經(jīng)過(guò)稀薄過(guò)渡流區(qū)和近連續(xù)滑移流區(qū)，并會(huì)在飛行器不同繞流區(qū)域因分子碰撞頻率巨大差別而出現(xiàn)連續(xù)流、過(guò)渡流或高稀薄流共存的混合流動(dòng)現(xiàn)象，導(dǎo)致物面邊界條件的表述方式對(duì)氣動(dòng)力/熱特性影響較大?；贐oltzmann模型方程的氣體動(dòng)理論統(tǒng)一算法(GKUA)需要求解的是物理空間與速度空間各個(gè)離散網(wǎng)格點(diǎn)處的氣體分子速度分布函數(shù)，于是各類(lèi)邊界條件與物面氣動(dòng)熱力學(xué)特性數(shù)學(xué)模型建立與計(jì)算實(shí)現(xiàn)[65,68,74-76,78]，始終通過(guò)求解氣體分子速度分布函數(shù)演化更新來(lái)進(jìn)行。任何固體壁面邊界均具有對(duì)氣體分子的不可浸透性，在任一計(jì)算時(shí)刻物面附近的分子要么撞擊物面，要么未撞擊物面，而所有撞擊物面的分子均會(huì)反射回氣體中，只是氣體與壁面相互作用與固體表面層的狀態(tài)有密切關(guān)系，壁面的溫度、粗糙度及清潔度等直接影響氣體分子與固體壁的相互作用。為了確定氣體分子與固體表面相互作用數(shù)學(xué)模型，便于數(shù)值處理，這里假定撞擊物面的氣體分子緊接著以完全調(diào)節(jié)于壁面溫度和速度的Maxwellian平衡態(tài)速度分布散射[64-69]。為此，物面邊界網(wǎng)格在t時(shí)刻的速度分布函數(shù)可計(jì)算模型為：

若c·n≥0，則表示氣體分子撞擊物面并發(fā)生了反射，經(jīng)物面反射的氣體分子速度分布函數(shù)的數(shù)值離散形式為:

其中，(Uw,Vw,Ww)為物面速度，一般為0。

氣體分子從物面散射的數(shù)密度nw先前是未知的，將隨著入射分子的速度分布及固體壁面情況而變化。本文利用物面上的質(zhì)量平衡條件[65-68,74]，即在壁面處應(yīng)用垂直于物面的零質(zhì)量流量條件決定

若c·n<0，則表示氣體分子未撞擊物面，此時(shí)氣體分子速度分布函數(shù)由流場(chǎng)內(nèi)鄰近邊界的幾排網(wǎng)格單邊二階差分實(shí)時(shí)計(jì)算確定[64-69,74]。

由此可見(jiàn)，統(tǒng)一算法設(shè)計(jì)的壁面邊界條件計(jì)算模型，避免了針對(duì)純粹微觀粒子進(jìn)行隨機(jī)統(tǒng)計(jì)模擬的DSMC方法產(chǎn)生統(tǒng)計(jì)漲落與對(duì)低Knudsen數(shù)近連續(xù)滑移過(guò)渡流難以模擬計(jì)算的問(wèn)題，也避開(kāi)了基于N-S方程等宏觀流體力學(xué)解算器純粹使用宏觀流動(dòng)量進(jìn)行邊界表述而不同位置的稀薄效應(yīng)強(qiáng)弱不同、難于用固定關(guān)系式數(shù)值實(shí)現(xiàn)的不足。

雖然適于航天再入繞流問(wèn)題Boltzmann模型方程統(tǒng)一算法計(jì)算空間，是由離散速度空間和位置空間組成的六維空間，并考慮內(nèi)部能級(jí)分布，形成多相空間，需要大量計(jì)算機(jī)內(nèi)存，使用文獻(xiàn)[30,77-78]發(fā)展的離散速度空間區(qū)域分解并行計(jì)算技術(shù)，結(jié)合MPI+Open MP多級(jí)并行[79]，可獲得求解Boltzmann模型方程高性能可擴(kuò)展大規(guī)模并行算法。根據(jù)氣體動(dòng)理論數(shù)值格式(12)式與基于離散速度空間分布函數(shù)宏觀流動(dòng)取矩離散求和，對(duì)任一時(shí)刻物理空間各點(diǎn)的宏觀流動(dòng)量進(jìn)行動(dòng)態(tài)確定與實(shí)時(shí)演化更新，可獲得位置空間各網(wǎng)格點(diǎn)處的氣動(dòng)力/熱繞流特性。

1.2 算法驗(yàn)證

為了考驗(yàn)求解Boltzmann模型方程統(tǒng)一算法對(duì)可回收類(lèi)航天器再入高稀薄流到近連續(xù)流區(qū)高超聲速流動(dòng)模擬能力，使用來(lái)自文獻(xiàn)[80]相同條件下可重復(fù)使用衛(wèi)星體飛行器(底部半徑R=503.5 mm，飛行器總長(zhǎng)L=1410 mm，錐身段半錐角θ=11.4°)，選取飛行器底部半徑為特征尺度。擬定Ma∞=5兩種繞流狀態(tài)H=70.8 km、Kn∞=0.002、Re∞=4074.37與H=105 km、Kn∞=0.74、Re∞=10.19，在位置空間網(wǎng)格51×25×31和速度空間網(wǎng)格60×40×40設(shè)置下，使用具有160MB/CPU的512CPU并行計(jì)算。圖2、圖3分別繪出上述兩種繞流流場(chǎng)軸對(duì)稱(chēng)面內(nèi)馬赫數(shù)等值線與繞流場(chǎng)、物面流線結(jié)構(gòu)，由圖看出，隨著飛行器飛行高度從H=105 km再入到H=70.8 km，氣體繞流由強(qiáng)烈的稀薄效應(yīng)逐漸過(guò)渡到近連續(xù)流區(qū)脫體激波、膨脹波系，飛行器繞流場(chǎng)受擾動(dòng)區(qū)域由全場(chǎng)過(guò)渡到繞衛(wèi)星體周?chē)^小區(qū)域；對(duì)H=105 km、Kn∞=0.74高稀薄流，飛行器繞流并不存在激波等強(qiáng)間斷現(xiàn)象與背風(fēng)旋渦回流結(jié)構(gòu)，而在飛行器周?chē)纬珊窈竦膹?qiáng)擾動(dòng)區(qū)域，氣流附著物面流動(dòng)；而對(duì)H=70.8 km、Kn∞=0.002近連續(xù)滑移流，飛行器繞流物體前面產(chǎn)生明晰清脆的脫體激波，罩在飛行器前面，包括駐點(diǎn)域、附著激波及飛行器底部出現(xiàn)真空區(qū)、尾渦流場(chǎng)結(jié)構(gòu)和飛行器后部流場(chǎng)再壓縮尾跡流動(dòng)現(xiàn)象均能很好地捕捉；截然不同的兩種流態(tài)反映了稀薄氣體繞流場(chǎng)完全不同于近連續(xù)流區(qū)繞流面貌。

圖2 衛(wèi)星體再入不同高度(Ma∞=5)繞流場(chǎng)馬赫數(shù)等值線Fig.2 Mach number contours over a sphere-cone shaped satellite spacecraft at Ma∞=5

圖3 衛(wèi)星體不同高度繞流(Ma∞=5)場(chǎng)與物面流線結(jié)構(gòu)Fig.3 Streamlines patterns over a sphere-cone shaped satellite spacecraft at Ma∞=5

圖4繪出統(tǒng)一算法計(jì)算該飛行器繞流(H=105 km、Ma∞=5、Tw/T∞=1)物面熱流分布與DSMC模擬值[80]定量化比較，可看出GKUA得到的結(jié)果從球頭前駐點(diǎn)沿物面向后不同物面角熱流計(jì)算值均與DSMC結(jié)果吻合較好，兩種結(jié)果偏差范圍0.27%～8.56%，且GKUA得到的迎風(fēng)物面熱流分布的非線性效應(yīng)更明顯。圖5繪出GKUA計(jì)算Ma∞=10上述衛(wèi)星體飛行器(H=75.9 km、Kn∞=0.004、Re∞=4074.5、Tw/T∞=1.64)繞流駐點(diǎn)線密度、溫度分布與典型文獻(xiàn)[80]結(jié)果定量化比較，可看出GKUA計(jì)算值(帶符號(hào)點(diǎn)劃線表示)與文獻(xiàn)結(jié)果(符號(hào)“○”表示)變化規(guī)律吻合很好，兩種結(jié)果計(jì)算偏差0.39%～3.26%。圖中顯示出由于設(shè)置Tw/T∞=1.64冷壁面，在飛行器前面一定區(qū)域會(huì)出現(xiàn)流動(dòng)高溫區(qū)，整體上看GKUA對(duì)駐點(diǎn)線流動(dòng)參數(shù)計(jì)算分辨率優(yōu)于DSMC模擬值。

圖4 衛(wèi)星體再入流(H=105 km,Kn∞=0.74,Ma∞=5)迎風(fēng)物面熱流分布GKUA與DSMC計(jì)算比較Fig.4 Comparison of heat flux distribution along windward surface around SARA satellite at H=105 km,Kn∞=0.74, Ma∞=5

圖5 衛(wèi)星體再入繞流(H=75.9 km,Ma∞=10,Tw/T∞=1.64)駐點(diǎn)線密度、溫度分布GKUA與參考文獻(xiàn)[80]結(jié)果比較Fig.5 Comparison of stagnation line distribution around SARA satellite shape at H=75.9 km,Ma∞=10,Tw/T∞=1.64

為了將求解Boltzmann模型方程統(tǒng)一算法應(yīng)用于飛船返回艙再入近連續(xù)過(guò)渡流區(qū)繞流問(wèn)題模擬，擬定H=80.9 km,Ma∞=12.7,Kn∞=0.0016,Re∞=12 175.4,Pr=0.72,α=20°繞流狀態(tài)，網(wǎng)格布局為位置空間流向×法向×周向45×16×25和速度空間130×125×100，使用具有75.12MB/CPU分布內(nèi)存的4096CPU并行計(jì)算。圖6繪出H=80.9 km,Ma∞=12.7,Kn∞=0.0016繞流流場(chǎng)壓力、溫度、熱流、流線及迎風(fēng)區(qū)物面熱流與背風(fēng)區(qū)極限流線分布。對(duì)此高超聲速近連續(xù)過(guò)渡流區(qū)繞流，離返回艙駐點(diǎn)前端一定距離出現(xiàn)較清晰的脫體激波，顯示出較明顯的近連續(xù)流動(dòng)特征。因壁溫設(shè)置Tw/T0=0.5435與Tw/T∞=18.0529，從圖6(b)看出在脫體激波與駐點(diǎn)前端出現(xiàn)高溫區(qū)。遠(yuǎn)前方氣流跨越脫體激波改變流動(dòng)方向，繞流物面，在返回艙肩部，以超聲速膨脹往后流去，在返回艙尾部出現(xiàn)流動(dòng)真空區(qū)。從圖6(c)看出在脫體激波與前駐點(diǎn)物面附近溫度劇烈變化區(qū)域出現(xiàn)高熱流區(qū)，而在其他部位熱流影響很小。這也說(shuō)明返回艙再入近連續(xù)過(guò)渡流區(qū)，因出現(xiàn)脫體激波，而將熱流匯聚在脫體激波與駐點(diǎn)區(qū)物面，為此，防熱設(shè)計(jì)需要考慮前體防熱結(jié)構(gòu)。從圖6(e,f)看出對(duì)此80.9 km高度，因空氣較強(qiáng)的稀薄效應(yīng)，回流區(qū)旋渦尚在發(fā)展中，該高超聲速繞流狀態(tài)Kn∞=0.0016、Ma∞=12.7、Re∞=12 175.4，因受到較強(qiáng)的稀薄效應(yīng)影響，氣流自前駐點(diǎn)繞流物面基本上附著物面流動(dòng)，僅在后駐點(diǎn)背風(fēng)真空區(qū)存在一定程度的回流擾動(dòng)，對(duì)背風(fēng)區(qū)的流動(dòng)參數(shù)分布有較強(qiáng)的影響。為進(jìn)一步分析GKUA對(duì)返回艙物面壓力分布計(jì)算的正確性，圖7(a)繪出該返回艙繞流不同子午面φ=0°、90°、180°壁面壓力分布，其中線段表示GKUA并行計(jì)算結(jié)果，符號(hào)“○”為N-S方程解算器計(jì)算得到子午面φ=90°的壁面壓力分布，可看出兩種方法得到的結(jié)果在迎風(fēng)區(qū)吻合較好，跨過(guò)肩部進(jìn)入背風(fēng)區(qū)由于出現(xiàn)真空區(qū)嚴(yán)重稀薄效應(yīng)，N-S解幾乎為常數(shù)，而GKUA結(jié)果逐漸減小到返回艙后端面壓力最小，兩種方法在迎風(fēng)區(qū)計(jì)算結(jié)果的吻合相容、背風(fēng)區(qū)不一致源于在背風(fēng)低壓真空區(qū)出現(xiàn)逆壓梯度、較強(qiáng)稀薄效應(yīng)繼續(xù)用連續(xù)流體力學(xué)N-S方程計(jì)算模型是產(chǎn)生差別的主要原因。相較N-S方程解算器，GKUA直接捕捉流場(chǎng)物面任意網(wǎng)格里的速度分布函數(shù)演化更新，所有宏觀流動(dòng)量均通過(guò)式(6)～式(10)的數(shù)值積分離散求和，實(shí)時(shí)表征更為客觀真實(shí)，也由于在六維空間計(jì)算求解，需要大量計(jì)算內(nèi)存，這是該方法的局限性，如這里使用的位置空間網(wǎng)格數(shù)較少。對(duì)照?qǐng)D7(b)與圖6(d)分別繪出的迎風(fēng)區(qū)物面壓力與熱流分布云圖，表明對(duì)此大鈍頭返回艙繞流迎風(fēng)物面力熱分布，壓力最大值發(fā)生在駐點(diǎn)、熱流最大值集中在曲率最大的肩部，這與已開(kāi)展的風(fēng)洞試驗(yàn)與飛行測(cè)試結(jié)果相吻合一致。

圖6 返回艙再入H=80.9 km, Ma∞=12.7,Kn∞=0.0016, α=20°繞流場(chǎng)壓力、溫度、熱流與流線結(jié)構(gòu)Fig.6 Flow field and surface streamline structures around reentry spacecraft Kn∞=0.0016, Ma∞=12.7, H=80.9 km

圖7 返回艙再入H=80.9 km, Ma∞=12.7, Kn∞=0.0016, α=20°繞流物面壓力分布Fig.7 Surface pressure distribution in different meridian planes of φ=0°，90°，180° past reentry spacecraft (H=80.9 km)

1.3 含轉(zhuǎn)動(dòng)非平衡效應(yīng)的統(tǒng)一算法計(jì)算驗(yàn)證

返回艙再入跨流域極高速、高溫繞流流場(chǎng)中，氣體分子的微觀自由度(平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)和電子態(tài))因受到一定程度的激勵(lì)，會(huì)出現(xiàn)能量彼此傳遞，使分子和原子間發(fā)生化學(xué)和電離反應(yīng)。氣體的宏觀運(yùn)動(dòng)和狀態(tài)變化同相應(yīng)的微觀物理化學(xué)過(guò)程相互影響呈現(xiàn)復(fù)雜的流動(dòng)非平衡效應(yīng)。根據(jù)表征分子微觀自由度之間能量傳遞或組元之間進(jìn)行化學(xué)反應(yīng)的特征弛豫時(shí)間與流動(dòng)特征時(shí)間大小尺度的不同，可將非平衡流動(dòng)分為平動(dòng)和轉(zhuǎn)動(dòng)非平衡流、振動(dòng)和化學(xué)非平衡流以及電離輻射非平衡流。如果特征流動(dòng)時(shí)間極小或流場(chǎng)的物理量變化梯度極大，平動(dòng)與轉(zhuǎn)動(dòng)非平衡效應(yīng)表現(xiàn)為與分子性質(zhì)相關(guān)聯(lián)的氣體介質(zhì)特性，如比熱、黏性系數(shù)、熱導(dǎo)率等不再保持常數(shù)，氣體運(yùn)動(dòng)的基本方程就與通常氣體動(dòng)力學(xué)連續(xù)性質(zhì)量、動(dòng)量、能量方程及狀態(tài)關(guān)系不同，出現(xiàn)黏性、熱傳導(dǎo)和擴(kuò)散的非平衡變化，這是一種最基本而接近高超聲速再入多尺度非平衡流動(dòng)現(xiàn)實(shí)環(huán)境。為此，如何構(gòu)造描述熱力學(xué)非平衡效應(yīng)的Boltzmann模型方程及數(shù)值求解途徑，是氣體動(dòng)理學(xué)理論可計(jì)算建模在極高超聲速再入空氣動(dòng)力學(xué)研究焦點(diǎn)。為了將上述求解Boltzmann模型方程統(tǒng)一算法推廣應(yīng)用于近地空間飛行環(huán)境跨流域非平衡繞流問(wèn)題研究，作為這方面研究探索的第一步，我們基于對(duì)轉(zhuǎn)動(dòng)自由度松弛變化[81]、氣體分子運(yùn)動(dòng)論Rykov模型[44,82]研究，在氣體分子速度分布函數(shù)演化更新求解中考慮轉(zhuǎn)動(dòng)自由度影響，把氣體分子轉(zhuǎn)動(dòng)能作為分布函數(shù)的自變量，提出考慮轉(zhuǎn)動(dòng)非平衡效應(yīng)Boltzmann模型方程數(shù)值算法[28-29,31]。

基于轉(zhuǎn)動(dòng)松弛特性Rykov模型[44.82]，采用轉(zhuǎn)動(dòng)慣量描述氣體分子自旋運(yùn)動(dòng)，利用分子總角動(dòng)量守恒作為一個(gè)新的碰撞不變量。將分子內(nèi)部能量作為氣體分子速度分布函數(shù)自變量，引入能量模式配分函數(shù)將能量在各自由度平均分布?；跉怏w分子速度分布函數(shù)f=f(r,V,t,e)，這里r、V分別是位置空間和速度空間坐標(biāo)、e為內(nèi)能，在求解Boltzmann模型方程統(tǒng)一算法框架[28-32,64-79]下，使用權(quán)因子1和e對(duì)速度分布函數(shù)依賴(lài)的速度空間無(wú)窮積分，引入能級(jí)約化速度分布函數(shù)f0與f1，去掉方程對(duì)內(nèi)部能量連續(xù)依賴(lài)關(guān)系，可確立含轉(zhuǎn)動(dòng)非平衡效應(yīng)各流域統(tǒng)一Boltzmann模型方程，其無(wú)量綱形式為：

分析推導(dǎo)反映黏性與熱傳導(dǎo)屬性的擴(kuò)散時(shí)間與氣體分子平均碰撞時(shí)間，描述含內(nèi)能激發(fā)的多原子分子彈性與非彈性碰撞松弛變化統(tǒng)一表達(dá)式，以此表征氣體分子速度分布函數(shù)趨于當(dāng)?shù)仄胶鈶B(tài)分布的碰撞松弛速率。

φ(t)=0.767+0.233t-1/6exp[-1.17(t-1)],

μt=(Tt)2/3φ(B)/φ(BTt),B=T∞/T*。

這里，T*與T∞一樣都是有量綱值。對(duì)于氮?dú)舛?，T*=91.5 K。

速度分布函數(shù)f0、f1同樣是依賴(lài)于時(shí)間t、位置空間r和速度空間V的函數(shù)，使用上節(jié)介紹的氣體動(dòng)理論離散速度坐標(biāo)法與位置空間數(shù)值求解格式，可得到求解式(20,21)考慮轉(zhuǎn)動(dòng)非平衡效應(yīng)Boltzmann模型方程統(tǒng)一算法。

由于空氣中的N2和O2分子振動(dòng)特征溫度遠(yuǎn)高于室溫，分別為3371 K和2256 K，通常情況下氣體分子振動(dòng)能模態(tài)不易激發(fā)，由此建立考慮轉(zhuǎn)動(dòng)非平衡效應(yīng)的氣體動(dòng)理論統(tǒng)一算法可在增加一個(gè)考慮轉(zhuǎn)動(dòng)能級(jí)約化速度分布函數(shù)f1較完全氣體多一倍計(jì)算量情況下實(shí)現(xiàn)任意雙原子氣體流動(dòng)問(wèn)題的模擬。為了考驗(yàn)該計(jì)算模型實(shí)現(xiàn)的正確性，擬定一類(lèi)返回艙外形體再入不同流域高超聲速繞流問(wèn)題。圖8(a,b)分別繪出了飛行高度H=105 km、Kn∞=0.091 93與H=83 km、Kn∞=0.001 95，馬赫數(shù)Ma∞=8，α=0°，γ=1.4，在流動(dòng)介質(zhì)氮?dú)庵蟹祷嘏撏庑误w繞流流場(chǎng)馬赫數(shù)等值線云圖與流線分布計(jì)算結(jié)果。圖中可見(jiàn)，對(duì)于數(shù)米量級(jí)大尺度真實(shí)飛行器，即使在H=105 km氣體分子平均自由程高達(dá)λ∞=0.3365 m的高度，因Kn∞=0.091 93，飛行器繞流仍表現(xiàn)在前端出現(xiàn)厚的脫體激波強(qiáng)擾動(dòng)，遠(yuǎn)前方來(lái)流跨越脫體激波層，流動(dòng)改變方向，繞流物面由于稀薄效應(yīng)很?chē)?yán)重，氣體繞流完全附著物面流動(dòng)；相較之下，對(duì)于圖8(b)繪出H=83 km的繞流圖譜，在此高度氣體分子平均自由程降低為λ∞=0.0071 m很小，Kn∞=0.001 95表現(xiàn)為近連續(xù)流，氣體繞流出現(xiàn)明顯的脫體激波，遠(yuǎn)前方來(lái)流跨越清脆明晰的脫體激波強(qiáng)間斷，流動(dòng)改變方向向外流去，繞流物面，跨越最高點(diǎn)進(jìn)入倒錐區(qū)，氣流產(chǎn)生膨脹波系，并在物體后部出現(xiàn)流動(dòng)分離、尾渦結(jié)構(gòu)等近連續(xù)繞流面貌。該圖直觀再現(xiàn)了返回艙外形體再入飛行不同流域發(fā)生的繞流現(xiàn)象，揭示了飛行器繞流隨著Kn∞數(shù)減小由稀薄流趨近連續(xù)流的變化過(guò)程。

圖8 再入不同高度H=105 km與H=83 km返回艙體繞流馬赫數(shù)等值線云圖與流線結(jié)構(gòu)分布(Ma∞=8)Fig.8 Mach number contours and streamline distribution around reentry spacecraft with Ma∞=8

圖9(a、b)繪出返回艙H=105 km、Ma∞=8、Kn∞=0.091 93、α=0°、γ=1.4條件下繞流流場(chǎng)壓力p/p∞等值線云圖與遠(yuǎn)場(chǎng)到前端駐點(diǎn)線密度ρ/ρ∞分布，圖9(a)上半部分、圖9(b)實(shí)線“GKUA”為統(tǒng)一算法結(jié)果，而下半部分及符號(hào)則是“DSMC”模擬流場(chǎng)分布?？梢?jiàn)兩種方法計(jì)算結(jié)果從流場(chǎng)結(jié)構(gòu)到壓力等值線、遠(yuǎn)前方來(lái)流到前駐點(diǎn)線密度分布均吻合很好，證實(shí)統(tǒng)一算法與DSMC方法適應(yīng)的再入高超聲速稀薄流模擬在刻畫(huà)Boltzmann方程描述的客觀物理過(guò)程方面具有很好的收斂一致性[83]。

圖9 再入返回艙體繞流流場(chǎng)壓力、駐點(diǎn)線密度分布GKUA與DSMC模擬驗(yàn)證(H=105 km, Ma∞=8)Fig.9 Validation of GKUA and DSMC on pressure and density distribution around reentry spacecraft with H=105 km and Ma∞=8

圖10 再入H=105 km返回艙體繞流Ma∞=8沿軸向不同位置物面無(wú)量綱壓力、熱流分布GKUA與DSMC模擬驗(yàn)證Fig.10 Validation of GKUA and DSMC on surface pressure and heat flux distribution along with the X-axis around reentry spacecraft with H=105 km and Ma∞=8

為了進(jìn)一步驗(yàn)證統(tǒng)一算法用于計(jì)算返回艙再入高超聲速繞流物面氣動(dòng)力/熱特性的準(zhǔn)確性，對(duì)H=105 km、α=0°返回艙體Ma∞=8繞流進(jìn)行DSMC與GKUA結(jié)果比較分析，圖10繪出兩種方法計(jì)算返回艙沿軸線方向的物面無(wú)量綱壓力系數(shù)p/p∞和熱流q分布對(duì)比情況，可看出，兩種方法得到的表面壓力和熱流吻合較好，變化趨勢(shì)也完全一致。分析表明，對(duì)此嚴(yán)重稀薄效應(yīng)繞流，迎風(fēng)區(qū)物面壓力、熱流值均遠(yuǎn)高于背風(fēng)區(qū)，物面壓力最大值出現(xiàn)在返回艙頭部駐點(diǎn)位置；而熱流最大值后移至返回艙肩部最高，過(guò)了拐點(diǎn)進(jìn)入背風(fēng)區(qū)，物面熱流則迅速下降。對(duì)于物面壓力分布，DSMC與GKUA結(jié)果幾乎完全重合；對(duì)于物面熱流分布，總體看出，GKUA在揭示具有巨大曲率過(guò)渡帶的返回艙肩部最高拐點(diǎn)熱流分布的分辨率明顯優(yōu)于DSMC結(jié)果，如GKUA肩部拐點(diǎn)熱流390 W/m2，而DSMC結(jié)果低于此結(jié)果達(dá)13%，這是因?yàn)樵诖饲屎艽蟮募绮抗拯c(diǎn)附近，繞流流場(chǎng)網(wǎng)格劃分需要自適應(yīng)加密，DSMC模擬特點(diǎn)決定較難像求解Boltzmann模型方程GKUA確定論數(shù)值算法易于實(shí)現(xiàn)，尤其在精細(xì)捕捉流動(dòng)信息方面，GKUA基于物理空間與速度空間網(wǎng)格布局劃分，直接對(duì)氣體分子速度分布函數(shù)演化更新數(shù)值求解，更能刻畫(huà)流動(dòng)細(xì)致結(jié)構(gòu)，同樣在跨越肩部進(jìn)入倒錐過(guò)渡背風(fēng)區(qū)，GKUA仍然較DSMC有更強(qiáng)捕捉微弱流動(dòng)傳熱模擬能力，顯示出GKUA結(jié)果較DSMC更好地光滑過(guò)渡到背風(fēng)熱流極小區(qū)，且在整個(gè)倒錐、后柱段，GKUA結(jié)果均高于DSMC模擬值，通過(guò)將GKUA與經(jīng)過(guò)實(shí)踐廣泛檢驗(yàn)較為成熟的DSMC兩種計(jì)算模型結(jié)果的比較驗(yàn)證，證實(shí)兩種方法具有較強(qiáng)的工程實(shí)用互補(bǔ)性。通過(guò)對(duì)圖8～圖10初步分析，證實(shí)求解轉(zhuǎn)動(dòng)非平衡效應(yīng)Boltzmann模型方程統(tǒng)一算法用于計(jì)算大型航天器氣動(dòng)力與氣動(dòng)熱繞流問(wèn)題具有較高可靠性與工程置信度，這為進(jìn)一步發(fā)展復(fù)雜高超聲速非平衡繞流問(wèn)題統(tǒng)一算法解算器，開(kāi)展返回艙往返大氣層跨越飛行各流域繞流問(wèn)題大規(guī)模并行計(jì)算，特別是近連續(xù)、過(guò)渡流區(qū)氣動(dòng)力/熱與姿態(tài)配平繞流問(wèn)題研究奠定了基礎(chǔ)。

1.4 含熱力學(xué)振動(dòng)能激發(fā)的Boltzmann模型方程統(tǒng)一算法計(jì)算驗(yàn)證與分析

返回艙再入高超聲速繞流引起嚴(yán)重的氣動(dòng)加熱高溫非平衡效應(yīng)，尤其是脫體激波與物面邊界層內(nèi)流動(dòng)氣體分子振動(dòng)能將被激發(fā)，出現(xiàn)量子能級(jí)分布，在振動(dòng)能量子效應(yīng)不可忽略情況下，需要考慮如何基于描述各流域氣體分子輸運(yùn)現(xiàn)象的Boltzmann方程出發(fā)，實(shí)施含平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能激發(fā)的熱力學(xué)非平衡效應(yīng)Boltzmann模型方程設(shè)計(jì)與氣動(dòng)力熱影響問(wèn)題模擬。早年王承書(shū)與Uhlenbeck提出了一種處理具有內(nèi)能影響的多原子氣體的半經(jīng)典方法，即平動(dòng)能根據(jù)自由度按經(jīng)典方法處理，而內(nèi)自由度按量子力學(xué)觀點(diǎn)處理，由此得到了Boltzmann方程的推廣形式—WCU方程[84-85]，該方程在分子內(nèi)能為非簡(jiǎn)并態(tài)時(shí)成立。Morse[86]在王承書(shū)與Uhlenbeck研究成果基礎(chǔ)上，將彈性碰撞與非彈性碰撞解耦用于處理能量松弛過(guò)程，依此構(gòu)造了一種類(lèi)似BGK模型的考慮分子內(nèi)能間斷能級(jí)的模型方程。在這些研究中，內(nèi)能作為單一模式處理，沒(méi)有區(qū)分轉(zhuǎn)動(dòng)能和振動(dòng)能。如果按量子力學(xué)觀點(diǎn)處理分子能級(jí)躍遷，每個(gè)能級(jí)均為一個(gè)獨(dú)立的分布函數(shù)，勢(shì)必會(huì)加大數(shù)值處理的難度。

為此，我們?cè)谇蠼釨oltzmann模型方程統(tǒng)一算法框架[64-69,74-79]下，為了模擬高溫高馬赫數(shù)下多原子氣體內(nèi)能激發(fā)對(duì)跨流域非平衡流動(dòng)的影響，將轉(zhuǎn)動(dòng)能、振動(dòng)能分別作為氣體分子速度分布函數(shù)中的自變量，把轉(zhuǎn)動(dòng)能和振動(dòng)能處理為連續(xù)分布的能量模式，將Boltzmann方程碰撞項(xiàng)分解成彈性碰撞和非彈性碰撞，同時(shí)將非彈性碰撞按一定松弛速率分解為平動(dòng)-轉(zhuǎn)動(dòng)能松弛和平動(dòng)-轉(zhuǎn)動(dòng)-振動(dòng)能松弛，構(gòu)造了一類(lèi)考慮振動(dòng)能激發(fā)的Boltzmann模型方程，并證明其守恒性和H定理。在此基礎(chǔ)[27-29,31,87]上，對(duì)分布函數(shù)基于內(nèi)部能量變量無(wú)窮積分，引入三個(gè)約化速度分布函數(shù)，得到一組考慮振動(dòng)能激發(fā)的約化速度分布函數(shù)控制方程組：

其中，i=1，2，3，且有

在確定含振動(dòng)能激發(fā)的氣體動(dòng)理學(xué)模型方程及離散速度坐標(biāo)法可計(jì)算模型后，將其加入到氣體動(dòng)理論統(tǒng)一算法求解體系，可建立同時(shí)考慮平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)能激發(fā)的熱力學(xué)非平衡Boltzmann模型方程統(tǒng)一算法。相較僅考慮轉(zhuǎn)動(dòng)非平衡效應(yīng)Boltzmann模型方程的氣體動(dòng)理論數(shù)值算法，通過(guò)分別相應(yīng)于平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)能松弛變化的三個(gè)約化速度分布函數(shù)f1、f2、f3的可計(jì)算表征實(shí)現(xiàn)，使得計(jì)算量較完全氣體Boltzmann模型方程的統(tǒng)一算法增加了三倍，依靠所發(fā)展的Boltzmann模型方程統(tǒng)一算法離散速度空間區(qū)域分解多級(jí)并行MPI+OpenACC程序架構(gòu)體系，可獲得模擬跨流域高超聲速氣動(dòng)力/熱繞流問(wèn)題的氣體動(dòng)理論統(tǒng)一算法大規(guī)模并行計(jì)算應(yīng)用研究平臺(tái)。

圖12繪出分別使用含平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)非平衡效應(yīng)Boltzmann模型方程GKUA與DSMC計(jì)算物面熱流與壓力分布對(duì)比情況，其中“GKUA-Shakhov”表示將N2作為簡(jiǎn)單氣體采用基于Shakhov模型的GKUA計(jì)算值，“GKUA-ES”表示僅考慮轉(zhuǎn)動(dòng)能激發(fā)而基于多原子氣體ES模型的GKUA結(jié)果，“GKUA-vibration”表示考慮N2振動(dòng)能激發(fā)的Boltzmann模型方程GKUA結(jié)果，“DSMC-translation”表示將N2作為簡(jiǎn)單氣體不考慮內(nèi)能激發(fā)DSMC模擬值，“DSMC-rotation”表示僅考慮N2轉(zhuǎn)動(dòng)能影響的DSMC模擬值，“DSMC-vibration”表示含N2振動(dòng)能激發(fā)DSMC結(jié)果。圖12(a)所示物面熱流分布看出，駐點(diǎn)附近“DSMC-translation”結(jié)果明顯高于其他幾種結(jié)果，且存在嚴(yán)重統(tǒng)計(jì)波動(dòng)，繞流圓柱面中部，“GKUA-vibration”偏離其他幾種結(jié)果，最大偏差不超過(guò)15%；圖12(b)所示物面壓力分布看出，駐點(diǎn)附近“DSMC-translation”結(jié)果最小，“GKUA-Shakhov”次之，其他結(jié)果幾乎重合，說(shuō)明內(nèi)能激發(fā)熱力學(xué)非平衡效應(yīng)對(duì)氣動(dòng)力影響很小。

圖11 含振動(dòng)能激發(fā)的Boltzmann模型方程GKUA計(jì)算圓柱繞流宏觀流動(dòng)量與DSMC方法結(jié)果對(duì)比Fig.11 Macroscopic flow variables circular cylinder simulated by GKUA with comparison of DSMC results

圖12 考慮平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)非平衡效應(yīng)的GKUA、DSMC計(jì)算圓柱繞流物面壓力和熱流分布, Ma∞=5,Kn∞=0.01Fig.12 Surface heat flux and pressure around circular cylinder with different non-equilibrium effects computed by GKUA and DSMC, Ma∞=5,Kn∞=0.01

為了考察含振動(dòng)能激發(fā)的Boltzmann模型方程統(tǒng)一算法對(duì)可回收類(lèi)航天器再入繞流氣動(dòng)力熱特性的影響，擬定與1.2節(jié)相同衛(wèi)星體飛行器以來(lái)流馬赫數(shù)5再入飛行兩種繞流狀態(tài)：來(lái)流Kn∞數(shù)分別為1、0.01。從前面的研究體會(huì)到，只有當(dāng)溫度較高時(shí)振動(dòng)能的激發(fā)才會(huì)對(duì)流動(dòng)產(chǎn)生作用，這里將來(lái)流溫度和物面溫度增加，研究考察振動(dòng)能激發(fā)所致熱力學(xué)非平衡效應(yīng)對(duì)流場(chǎng)結(jié)構(gòu)影響問(wèn)題。為此設(shè)置來(lái)流溫度為500 K、壁溫Tw=2000 K。圖13、圖14分別繪出上述兩種狀態(tài)繞流場(chǎng)振動(dòng)溫度、等效溫度、壓力和馬赫數(shù)等值線分布云圖?？煽闯?，對(duì)較高Kn∞=1稀薄氣體繞流，振動(dòng)溫度與全場(chǎng)等效溫度差別較大，達(dá)700 K，說(shuō)明高Kn∞數(shù)氣體分子數(shù)較少，能量松弛過(guò)程通過(guò)分子間碰撞完成，較少的分子數(shù)致能量松弛過(guò)程變緩，非平衡效應(yīng)越發(fā)顯著；另一方面，若流動(dòng)趨于Kn∞=0.01近連續(xù)流，由于連續(xù)介質(zhì)流中，氣體分子速度分布函數(shù)呈現(xiàn)當(dāng)?shù)豈axwell平衡態(tài)分布，流場(chǎng)中平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)溫度接近或趨于與等效溫度相一致的分布，這與分別考慮平動(dòng)或轉(zhuǎn)動(dòng)非平衡的跨流域繞流現(xiàn)象相一致[28-29]。圖中還清晰反映了從高稀薄流過(guò)渡到近連續(xù)流演變過(guò)程中流場(chǎng)變化規(guī)律，再入飛行器前部一定區(qū)域出現(xiàn)脫體強(qiáng)擾動(dòng)逐漸演變成清晰的激波強(qiáng)間斷、由稀薄低壓過(guò)渡到高壓繞流特征。

圖13 含振動(dòng)能激發(fā)Boltzmann模型方程GKUA計(jì)算衛(wèi)星體再入高稀薄繞流流場(chǎng)分布,Kn∞=1、Ma∞=5Fig.13 Flow distribution in vibrational energy excitation around satellite shape with Kn∞=1,Ma∞=5 computed by GKUA

圖14 含振動(dòng)能激發(fā)Boltzmann模型方程GKUA計(jì)算衛(wèi)星體再入近連續(xù)繞流流場(chǎng)分布,Kn∞=0.01、Ma∞=5Fig.14 Flow distribution in vibrational energy excitation around satellite shape with Kn∞=0.01,Ma∞=5 computed by GKUA

圖15、圖16分別繪出Kn∞=1、0.01兩種再入非平衡繞流狀態(tài)下衛(wèi)星體表面壓力和熱流等值線分布云圖?？煽闯?，隨著來(lái)流Kn∞數(shù)由1減小到0.01，來(lái)流由高稀薄流過(guò)渡到近連續(xù)流，物面壓力、熱流的有量綱值均急劇增大，對(duì)此衛(wèi)星體再入飛行器從Kn∞=1再入到Kn∞=0.01，物面壓力分布就增大了125倍、熱流增大近20倍。

圖15 含振動(dòng)能激發(fā)Boltzmann模型方程GKUA計(jì)算衛(wèi)星體再入物面壓力與熱流分布,Kn∞=1、Ma∞=5Fig.15 Pressure and heat flux distribution in vibrational energy excitation around satellite surface with Kn∞=1,Ma∞=5 computed by GKUA

圖16 含振動(dòng)能激發(fā)Boltzmann模型方程GKUA計(jì)算衛(wèi)星體再入物面壓力與熱流分布,Kn∞=0.01、Ma∞=5Fig.16 Pressure and heat flux distribution in vibrational energy excitation around satellite surface with Kn∞=0.01,Ma∞=5 computed by GKUA

2 跨流域空氣動(dòng)力學(xué)模擬方法

上節(jié)介紹的氣體動(dòng)理論統(tǒng)一算法作為近二十年來(lái)由Boltzmann型速度分布函數(shù)方程對(duì)完全氣體到含內(nèi)能激發(fā)包括轉(zhuǎn)動(dòng)能、振動(dòng)能熱力學(xué)非平衡流動(dòng)一維、二維、三維直至復(fù)雜結(jié)構(gòu)航天器再入過(guò)程可計(jì)算建模，發(fā)展至今，逐漸成為日趨廣泛用于航天軍工[26-32,68-70]、MEMS微尺度[66-67,88-89]行之有效的跨流域空氣動(dòng)力學(xué)模擬方法一種新途徑，該方法將氣體分子速度分布函數(shù)在位置空間和速度空間的信息保存起來(lái)進(jìn)行數(shù)值求解，所需要的內(nèi)存、計(jì)算量較大，仍在進(jìn)一步發(fā)展中。無(wú)論是載人航天器還是跨大氣層飛行器，以及各類(lèi)返回器的再入過(guò)程，其飛行航跡都是采取氣動(dòng)方式來(lái)控制。由于在稀薄氣體流域飛行的時(shí)間較長(zhǎng)，飛行速度非常高，如2014年10月我國(guó)圓滿成功完成試驗(yàn)任務(wù)的探月工程三期再入返回飛行試驗(yàn)器的再入速度達(dá)到10.6 km/s，而美國(guó)彗星探測(cè)器Stardust[5]返回再入速度更是高達(dá)12.8 km/s。如此極端速度條件下，氣體分子間較高的相對(duì)運(yùn)動(dòng)能量足以發(fā)生化學(xué)反應(yīng)和電離過(guò)程，燒蝕和所謂的“通信黑障”問(wèn)題更加嚴(yán)重。流場(chǎng)稀薄氣體效應(yīng)影響很大，飛行器前方繞流流場(chǎng)，特別是前方弓形激波附近流場(chǎng)，氣體分子狀態(tài)遠(yuǎn)離平衡態(tài)，不僅氣體分子的平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度存在顯著差別，且氣體分子間的化學(xué)反應(yīng)、電離輻射效應(yīng)劇烈[34]。能否準(zhǔn)確預(yù)測(cè)超高速再入導(dǎo)致的稀薄過(guò)渡流區(qū)嚴(yán)重的熱化學(xué)電離非平衡效應(yīng)對(duì)飛行器氣動(dòng)性能的影響，一直是對(duì)再入空氣動(dòng)力學(xué)研究的挑戰(zhàn)。為此，本文立足既有空氣動(dòng)力學(xué)研究手段，發(fā)展適于返回艙(器)再入各流域相應(yīng)數(shù)值模擬方法，使之形成驗(yàn)證結(jié)合互為補(bǔ)充綜合分析應(yīng)用研究平臺(tái)。

2.1 稀薄流區(qū)DSMC方法

DSMC方法用一個(gè)仿真分子代表一定數(shù)量FN的真實(shí)分子(包括分子、原子、離子和電子)，計(jì)算機(jī)存儲(chǔ)仿真分子位置坐標(biāo)、速度分量和內(nèi)能信息，隨仿真分子運(yùn)動(dòng)、與邊界相互作用及分子間碰撞而改變，最后統(tǒng)計(jì)網(wǎng)格內(nèi)仿真分子狀態(tài)實(shí)現(xiàn)真實(shí)氣體流動(dòng)的模擬，較容易引入符合物理實(shí)際的模型。方法關(guān)鍵是在小于當(dāng)?shù)胤肿悠骄鲎矔r(shí)間步長(zhǎng)Δt內(nèi)將分子運(yùn)動(dòng)和碰撞解耦，即讓所有分子運(yùn)動(dòng)一定距離并考慮邊界反射，然后計(jì)算此時(shí)間段內(nèi)具有代表性分子間碰撞。由于計(jì)算效率問(wèn)題，DSMC方法主要用于高稀薄過(guò)渡流的模擬，為了擴(kuò)展DSMC方法在近連續(xù)過(guò)渡流區(qū)模擬熱化學(xué)非平衡復(fù)雜流動(dòng)的能力和計(jì)算效率，根據(jù)高超聲速壓縮流動(dòng)中連續(xù)流假設(shè)失效的判斷準(zhǔn)則，將流場(chǎng)區(qū)域分解為氣體分子速度分布函數(shù)處于平衡態(tài)分布的連續(xù)流區(qū)和強(qiáng)熱非平衡的稀薄流區(qū)。在連續(xù)/近連續(xù)流區(qū)域，采用限制碰撞的DSMC方法，使用網(wǎng)格自適應(yīng)和變時(shí)間步長(zhǎng)技術(shù)來(lái)減少數(shù)值擴(kuò)散的誤差。對(duì)一些流動(dòng)問(wèn)題，可以在保證計(jì)算精度的基礎(chǔ)上提高計(jì)算效率。同時(shí)發(fā)展大規(guī)模并行算法[90]，形成了適于模擬高超聲速近連續(xù)流區(qū)氣動(dòng)特性的碰撞限制器DSMC混合方法[72]。文獻(xiàn)[91]通過(guò)計(jì)算不同網(wǎng)格尺度、不同壁面反射模型的返回艙再入條件下表面參數(shù)分布和氣動(dòng)力系數(shù)隨迎角變化、真實(shí)氣體效應(yīng)的影響，驗(yàn)證了算法可靠性。

稀薄氣體電離現(xiàn)象是航天器再入速度持續(xù)增大面臨的新問(wèn)題，直接導(dǎo)致通信黑障等傳統(tǒng)難題向稀薄區(qū)域大幅延伸。傳統(tǒng)的再入物理通過(guò)求解含化學(xué)反應(yīng)的N-S方程組[92]得到再入體繞流電子數(shù)密度分布，但是該方程組對(duì)于稀薄氣體流動(dòng)失效。近年來(lái)發(fā)展的稀薄氣體數(shù)值模擬方法中，基于Boltzmann方程的分析和計(jì)算[28-32,39,40,93,68,87]，目前尚無(wú)法包含電離過(guò)程，只有基于粒子隨機(jī)統(tǒng)計(jì)模擬的DSMC(Direct Simulation Monte Carlo)方法[5,21-22,33-34]是最有可能解決稀薄氣體電離數(shù)值預(yù)測(cè)難題的手段?；谙∮薪M分權(quán)重格式處理稀薄氣體含電離化學(xué)反應(yīng)過(guò)程，對(duì)不同權(quán)重粒子之間的碰撞根據(jù)概率改變內(nèi)能狀態(tài)，對(duì)不同權(quán)重粒子之間的化學(xué)反應(yīng)依概率刪除或保留反應(yīng)物粒子，同時(shí)對(duì)生成物中的稀有組分粒子進(jìn)行復(fù)制?；贛PI并行環(huán)境、直角/非結(jié)構(gòu)混合網(wǎng)格和網(wǎng)格自適應(yīng)技術(shù)，建立了適于三維真實(shí)復(fù)雜外形航天器極高超聲速再入電離熱化學(xué)非平衡流動(dòng)過(guò)程DSMC模擬平臺(tái)[7,21,33-34,72,91]。首次計(jì)算考察了月地高速再入返回器和神舟飛船返回艙再入稀薄氣體電離特性，分析了稀薄弱電離等離子體電子振動(dòng)頻率對(duì)應(yīng)的電磁波頻率，計(jì)算結(jié)果表明月地高速再入返回器和神舟飛船返回艙將分別在約85 km和80 km高度開(kāi)始發(fā)生通信中斷，這與飛行試驗(yàn)觀測(cè)結(jié)果一致，證實(shí)了極高超聲速再入電離熱化學(xué)非平衡流DSMC模擬平臺(tái)對(duì)解決國(guó)家需求的工程實(shí)用性[21,33-34,72,94]。

2.2 連續(xù)/近連續(xù)流區(qū)(滑移)N-S解算器

在返回艙再入近空間連續(xù)/近連續(xù)流區(qū)，為了借助N-S方程數(shù)值求解，揭示高空繞流物面稀薄非平衡效應(yīng)的滑移邊界模型對(duì)流場(chǎng)計(jì)算影響，發(fā)展了含熱解燒蝕與壁面催化、滑移稀薄效應(yīng)的熱化學(xué)N-S方程數(shù)值方法[37,95-97]，并從完全氣體和單一組分逐漸發(fā)展到考慮多組分以及振動(dòng)非平衡影響的情況， 實(shí)現(xiàn)了滑移邊界修正，并用于近連續(xù)流和滑移流區(qū)飛行器氣動(dòng)特性的預(yù)測(cè)和模擬[100]。通過(guò)構(gòu)造高溫多組分混合氣體物面滑移邊界條件數(shù)學(xué)模型，在Gupta等[99]關(guān)于多組分滑移邊界條件研究基礎(chǔ)上，開(kāi)展考慮振動(dòng)能激發(fā)非平衡與物面滑移、稀薄氣體間斷粒子效應(yīng)、黏性干擾效應(yīng)可計(jì)算建模，發(fā)展了一套考慮振動(dòng)非平衡、黏性干擾與稀薄氣體物面滑移效應(yīng)的飛行器再入近連續(xù)滑移流區(qū)高超聲速非平衡繞流物面氣動(dòng)熱特性N-S方程數(shù)值算法[37]。為了比較驗(yàn)證上述N-S方程數(shù)值算法對(duì)返回艙在氮?dú)饬鲃?dòng)介質(zhì)再入60公里不同來(lái)流迎角繞流實(shí)驗(yàn)狀態(tài)：來(lái)流馬赫數(shù)Ma∞=9.5、克努森數(shù)Kn∞,L=1.13×10-5、雷諾數(shù)Re∞,D=1.14×105，表1列出來(lái)流迎角α=10°、15°、18°、20°、22°條件下，常規(guī)N-S方程解算器與考慮壁面滑移效應(yīng)修正的N-S方程數(shù)值算法即表中分別備注“計(jì)算”、“滑移”的計(jì)算結(jié)果軸向力、法向力、俯仰力矩、壓心、升力、阻力系數(shù)及升阻比與N2實(shí)驗(yàn)條件下的低密度風(fēng)洞實(shí)驗(yàn)測(cè)量結(jié)果對(duì)比分析，可看出計(jì)算與實(shí)驗(yàn)總體上吻合較好，表中最后一列給出了各迎角狀態(tài)下兩種途徑計(jì)算得到的升阻比與實(shí)驗(yàn)數(shù)據(jù)相對(duì)偏差范圍11.97%～3.92%，考核驗(yàn)證了計(jì)算模型與數(shù)值方法實(shí)現(xiàn)的正確性。計(jì)算發(fā)現(xiàn)，對(duì)此H=60 km返回艙大鈍頭繞流，在α=20°以?xún)?nèi)，經(jīng)滑移邊界修正的N-S解算器結(jié)果較常規(guī)N-S計(jì)算結(jié)果更好吻合風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)，然而迎角更大，考慮壁面滑移效應(yīng)修正的N-S算法結(jié)果與實(shí)驗(yàn)數(shù)據(jù)偏差增大，精度不如常規(guī)N-S解算器，說(shuō)明滑移N-S解算器對(duì)大鈍頭迎風(fēng)區(qū)較簡(jiǎn)單物形繞流計(jì)算具有較好適應(yīng)性，然而對(duì)大迎角出現(xiàn)流動(dòng)分離復(fù)雜結(jié)構(gòu)體繞流局部物面稀薄滑移效應(yīng)的表征困難而致計(jì)算局限性。從表1反映出，針對(duì)該大鈍頭返回艙外形體繞流，迎角增大，會(huì)致迎風(fēng)物面法向力、升力系數(shù)增大；背風(fēng)區(qū)分離干擾區(qū)增大，致軸向力、阻力系數(shù)有所減小，以致升阻比隨迎角增大而明顯增大，為此大迎角情況下較大的升阻比計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)偏差絕對(duì)值與較小迎角情況下即使相當(dāng)，但會(huì)出現(xiàn)計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的相對(duì)偏差在大迎角情況比小迎角情況減小的特點(diǎn)。圖17(a、b)繪出上述返回艙實(shí)驗(yàn)?zāi)Ｐ驮偃?0 km、α=20°條件下繞流場(chǎng)馬赫數(shù)、壓力等值線云圖結(jié)構(gòu)，可看出在較大迎角情況下迎風(fēng)區(qū)壓力遠(yuǎn)遠(yuǎn)高于背風(fēng)區(qū)，稀薄效應(yīng)與背風(fēng)區(qū)分離多流區(qū)干擾致滑移N-S計(jì)算不適應(yīng)。

表1 返回艙模型繞流計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比Table 1 Comparison of Cal. results and Exp. data forre-entry capsule

圖17 返回艙試驗(yàn)?zāi)Ｐ屠@流場(chǎng)結(jié)構(gòu),Ma∞=9.5,α=20°Fig.17 Hypersonic flow past re-entry capsule (Kn∞,L=1.13×10-5,Ma∞=9.5,H=60 km)

2.3 近連續(xù)過(guò)渡流區(qū)N-S/DSMC耦合算法

為了開(kāi)展不同理論模型與數(shù)值方法間的耦合補(bǔ)充，通過(guò)研究連續(xù)流體力學(xué)N-S方程隨著當(dāng)?shù)乜伺瓟?shù)增大與復(fù)雜非規(guī)則物面稀薄效應(yīng)增強(qiáng)致其失效的判斷準(zhǔn)則為基礎(chǔ)，發(fā)展了適于近連續(xù)過(guò)渡流區(qū)高超聲速繞流流場(chǎng)分區(qū)N-S與DSMC耦合模擬技術(shù)[35]。根據(jù)高超聲速壓縮流動(dòng)特點(diǎn)，選取基于梯度變化的當(dāng)?shù)豄n∞數(shù)作為連續(xù)流方程失效判斷參數(shù)，根據(jù)氣體分子當(dāng)?shù)仄骄杂沙膛c流動(dòng)梯度定義的特征尺度比值，確定流動(dòng)當(dāng)?shù)豄n∞數(shù)，利用Kn∞l≥0.02判斷連續(xù)流方程失效，發(fā)展連續(xù)流方程失效界面選取法。構(gòu)建MPC耦合計(jì)算模型，在統(tǒng)一網(wǎng)格系統(tǒng)下，N-S方程對(duì)整個(gè)流場(chǎng)進(jìn)行計(jì)算，而DSMC只計(jì)算連續(xù)流方程失效部分區(qū)域。在分界面，兩種方法進(jìn)行信息交換，N-S方程結(jié)果為DSMC方法提供邊界條件，DSMC模擬值替換N-S方程在DSMC區(qū)域的計(jì)算結(jié)果。構(gòu)造抑制DSMC統(tǒng)計(jì)波動(dòng)對(duì)N-S計(jì)算影響的亞松弛技術(shù)，對(duì)N-S方程和DSMC方法計(jì)算程序進(jìn)行基于網(wǎng)格與信息交換融合分析設(shè)計(jì)，解決分區(qū)耦合計(jì)算信息傳遞問(wèn)題，DSMC方法計(jì)算區(qū)域由計(jì)算過(guò)程中流場(chǎng)參數(shù)自動(dòng)選取，根據(jù)N-S方程結(jié)果不斷調(diào)整DSMC方法計(jì)算區(qū)域，直到流場(chǎng)達(dá)到穩(wěn)定。建立了一套適于復(fù)雜飛行器繞流氣動(dòng)特性模擬的N-S/DSMC耦合算法。為了驗(yàn)證確認(rèn)該N-S/DSMC耦合模擬技術(shù)與求解Boltzmann模型方程GKUA、DSMC方法在近連續(xù)過(guò)渡流區(qū)復(fù)雜構(gòu)型繞流問(wèn)題求解的可靠性，擬定來(lái)自文獻(xiàn)[100]圖18(a)所示15°尖前緣高超聲速中空柱裙流動(dòng)算例：Ma∞=9.91、T∞=51 K、p∞=6.3 Pa、Tw=293 K。流動(dòng)介質(zhì)N2，Lref=0.1017 m，Kn∞=7.9135×10-4、Re∞,L=18 916，分別使用上述三種方法對(duì)該算例進(jìn)行計(jì)算比較。

圖18(b～d)繪出三種算法對(duì)該中空柱裙繞流流場(chǎng)壓力分布等值線云圖比較情況。在流場(chǎng)結(jié)構(gòu)上，三種結(jié)果基本一致。在流場(chǎng)細(xì)節(jié)上，三種結(jié)果有一定差異，氣流在尖前緣處產(chǎn)生斜激波，經(jīng)過(guò)斜激波在壓縮拐角區(qū)域形成λ波流場(chǎng)結(jié)構(gòu)，三種算法較好地捕捉到了這種復(fù)雜結(jié)構(gòu)。N-S/DSMC耦合算法與DSMC方法的結(jié)果更為接近，在遠(yuǎn)離物面與流場(chǎng)背風(fēng)區(qū)存在一些微弱波動(dòng)，而GKUA結(jié)果具有更高分辨率，更為光滑再現(xiàn)局部繞流結(jié)構(gòu)。圖19(a)繪出該中空柱裙繞流流向位置x/L=0.6處氣流密度沿橫向y方向分布的三種方法計(jì)算結(jié)果比較情況，可看出三種結(jié)果變化規(guī)律總體吻合較好，特別是本項(xiàng)目GKUA與DSMC模擬值從邊界層流動(dòng)到跨越斜激波以至遠(yuǎn)場(chǎng)區(qū)均較為一致，N-S/DSMC耦合算法結(jié)果在壁面邊界層及跨越激波過(guò)渡帶較前者差別或高或低一些，反映了N-S/DSMC耦合算法在模擬稀薄非平衡效應(yīng)嚴(yán)重的流場(chǎng)，存在較大的系統(tǒng)誤差，為此我們將該方法定位于近連續(xù)過(guò)渡流區(qū)一種工程預(yù)測(cè)分析手段，能較好補(bǔ)充與彌補(bǔ)跨流域空氣動(dòng)力學(xué)模擬手段途徑的完備性。圖19(b)繪出該中空柱裙繞流物面壓力分布比較情況，看出三種方法得到的壓力系數(shù)分布較為一致，物面壓力最大值在斜板上部拐角處，僅在極值點(diǎn)上，DSMC方法模擬結(jié)果較高，GKUA與N-S/DSMC耦合算法結(jié)果相當(dāng)一致，有待進(jìn)一步確認(rèn)造成該拐角處DSMC模擬與GKUA偏差更大的原因，是否因DSMC統(tǒng)計(jì)波動(dòng)所致。計(jì)算表明三種方法均能獲得較好的壁面壓力分布，不僅證實(shí)了本文前述Boltzmann方程可計(jì)算建模GKUA用于計(jì)算復(fù)雜構(gòu)型跨流區(qū)繞流非平衡相當(dāng)嚴(yán)重的流場(chǎng)與物面繞流問(wèn)題準(zhǔn)確可靠性，而且反映出本項(xiàng)目DSMC、N-S/DSMC耦合方法用于各自所適應(yīng)的流動(dòng)領(lǐng)域復(fù)雜飛行器繞流問(wèn)題模擬研究的工程可行性。

(a) Configuration of hollow column skirt

(b) Pressure contours from N-S/DSMC

(c) Pressure contours from GKUA

(d) Pressure contours from DSMC

(a) Lateral density x/L=0.6

(b) Surface pressure

2.4 返回艙模型低密度風(fēng)洞測(cè)力技術(shù)

高超聲速低密度風(fēng)洞是模擬航天器高空、高速飛行狀態(tài)的一種地面試驗(yàn)設(shè)備[101]，在中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所新建高超聲速低密度風(fēng)洞噴管出口尺寸為1 m量級(jí)，可模擬帶凸起物飛船返回艙真實(shí)外形，準(zhǔn)確預(yù)測(cè)其氣動(dòng)力特性。為此，針對(duì)返回艙外形模型，開(kāi)展了系列低密度風(fēng)洞測(cè)力試驗(yàn)相關(guān)技術(shù)研究，包括開(kāi)展小升阻比縱向三分量天平設(shè)計(jì)技術(shù)以及天平、模型防隔熱技術(shù)研究，發(fā)展高溫流場(chǎng)下高精度低升阻比測(cè)力天平設(shè)計(jì)技術(shù)，建立大鈍頭模型稀薄過(guò)渡流區(qū)氣動(dòng)力測(cè)試技術(shù)與實(shí)驗(yàn)方法。在此基礎(chǔ)上，開(kāi)展返回艙稀薄過(guò)渡流區(qū)不同高度、不同馬赫數(shù)與迎角下氣動(dòng)力特性研究。設(shè)計(jì)了帶凸起物(穩(wěn)定翼、俯仰發(fā)動(dòng)機(jī)座、分離插座板)外形試驗(yàn)?zāi)Ｐ?，模型總長(zhǎng)150.0 mm、最大截面直徑151.02 mm。試驗(yàn)?zāi)Ｐ筒捎梅侄卧O(shè)計(jì)，主要分為頭部、中段和后部等3個(gè)主要部分，如圖20(a)所示?？紤]到天平元件在縱向尺寸(長(zhǎng)度)上受到限制的現(xiàn)實(shí)，盡量避免天平與模型、支桿的連接錐面占用有限的模型縱向空間。將天平與模型中段的連接部分在不影響模型頭部與中段裝配的情況下盡量伸入模型頭部?jī)?nèi)腔，同時(shí)將天平與支桿的連接錐面置于模型外的支桿內(nèi)腔，使有限的模型內(nèi)腔空間用來(lái)布置天平測(cè)力元件。針對(duì)天平載荷特點(diǎn)，天平測(cè)力元件結(jié)構(gòu)盡量采取簡(jiǎn)單對(duì)稱(chēng)的形式進(jìn)行布局。這樣，一方面可有效減小熱結(jié)構(gòu)變形對(duì)天平測(cè)量產(chǎn)生的干擾影響，另一方面也便于天平機(jī)械加工和減小分量間的相互干擾影響。故在進(jìn)行天平結(jié)構(gòu)初步設(shè)計(jì)時(shí)，將軸向力元件布局在天平設(shè)計(jì)中心處，在天平設(shè)計(jì)中心前后處對(duì)稱(chēng)設(shè)置兩個(gè)矩形截面梁，用來(lái)測(cè)量法向力和俯仰力矩，并盡量縮短矩形梁的長(zhǎng)度。為減小天平的溫度效應(yīng)，在應(yīng)變計(jì)選用上使用溫度自補(bǔ)償?shù)闹袦貞?yīng)變計(jì)。對(duì)于高溫流場(chǎng)環(huán)境下使用的小尺寸應(yīng)變天平，熱對(duì)天平準(zhǔn)確測(cè)量的影響較大。天平測(cè)量梁受熱產(chǎn)生熱應(yīng)力，天平加工材料彈性模量受熱發(fā)生變化及天平測(cè)量梁上電阻應(yīng)變片受熱電阻和靈敏度系數(shù)發(fā)生變化等因素，都會(huì)對(duì)天平輸出信號(hào)產(chǎn)生干擾。同樣，天平結(jié)構(gòu)不同使得天平傳熱途徑發(fā)生變化，從而導(dǎo)致天平各測(cè)量梁溫度分布不同，產(chǎn)生不同的溫度影響。對(duì)高溫流場(chǎng)條件下小尺寸應(yīng)變天平技術(shù)進(jìn)行深入研究時(shí)，需要考慮天平結(jié)構(gòu)和傳熱之間的相互耦合影響。為此，需要開(kāi)展天平防隔熱與天平元件結(jié)構(gòu)布局研究。針對(duì)天平和模型間采用簡(jiǎn)單隔熱套結(jié)構(gòu)的方案，在風(fēng)洞來(lái)流條件下，采用有限元方法，對(duì)模型、天平及隔熱套整體結(jié)構(gòu)進(jìn)行傳熱分析，獲取各部件的溫升分布情況，并將分析結(jié)果同試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證有限元方法對(duì)該類(lèi)問(wèn)題進(jìn)行分析的有效性。針對(duì)不同的天平元件結(jié)構(gòu)布局形式和天平測(cè)量梁幾何尺寸，在天平設(shè)計(jì)載荷下，采用有限元方法，對(duì)天平的強(qiáng)度、剛度、靈敏度進(jìn)行詳細(xì)分析，選取綜合性能較優(yōu)的結(jié)構(gòu)方案，實(shí)現(xiàn)天平結(jié)構(gòu)優(yōu)化設(shè)計(jì)[102]。

在滿足模型外形相似情況下，試驗(yàn)選取克努森數(shù)Kn∞作為主要模擬參數(shù)，在流場(chǎng)名義馬赫數(shù)和前室總溫確定的情況下，按照風(fēng)洞試驗(yàn)流場(chǎng)的克努森數(shù)與高空飛行狀態(tài)的克努森數(shù)相同的原則，得出風(fēng)洞試驗(yàn)時(shí)的總壓，以此總壓值和確定的實(shí)際馬赫數(shù)、總溫作為風(fēng)洞運(yùn)行參數(shù)來(lái)模擬高空相應(yīng)飛行狀態(tài)：高度H=57.96 km、來(lái)流馬赫數(shù)Ma∞=9.95、總壓p0=1.95×106Pa、總溫T0=1070 K、來(lái)流雷諾數(shù)Re∞,D=1.7821×105、克努森數(shù)Kn∞,L=8.117×10-5。試驗(yàn)結(jié)果數(shù)據(jù)如表1相關(guān)欄目所示，圖20(a)繪出試驗(yàn)?zāi)Ｐ图疤炱?、模型防隔熱結(jié)構(gòu)示意圖，圖20(b～f)繪出對(duì)應(yīng)總壓p0=1.95 MPa、總溫T0=1100 K、來(lái)流迎角α=10°、15°、18°、20°、22°風(fēng)洞試驗(yàn)所測(cè)紋影照片，可看出不同來(lái)流迎角所呈現(xiàn)的脫體激波流場(chǎng)結(jié)構(gòu)，其相應(yīng)的氣動(dòng)力系數(shù)與(滑移)N-S解算器比較驗(yàn)證情況見(jiàn)表1所示，證實(shí)計(jì)算與實(shí)驗(yàn)?zāi)M可靠性。

圖20 風(fēng)洞試驗(yàn)?zāi)Ｐ团c測(cè)力天平布局及試驗(yàn)紋影照片F(xiàn)ig.20 Wind tunnel test model and force balance layout and test schlieren photograph for Shenzhou-type re-entry capsule

3 跨流域空氣動(dòng)力學(xué)模擬方法驗(yàn)證分析平臺(tái)

在研究建立求解Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法(GKUA)、DSMC、N-S/DSMC、(滑移)N-S解算器、低密度風(fēng)洞實(shí)驗(yàn)測(cè)試多種空氣動(dòng)力學(xué)模擬手段基礎(chǔ)上，為了驗(yàn)證確認(rèn)跨流域空氣動(dòng)力學(xué)模擬方法計(jì)算精度與不同方法模型特點(diǎn)，擬定圓球再入高度H=110～30 km，Kn∞=4.43～2.48×10-5，Ma∞=0.75～3.2，Re∞=0.25～1.37×105條件下十三組狀態(tài)，分別采用N-S方程解、N-S/DSMC耦合算法、DSMC與GKUA對(duì)全飛行流域不同高度氣動(dòng)阻力系數(shù)計(jì)算比較，表2給出GKUA計(jì)算分別與30～60 km連續(xù)流區(qū)N-S解、65～75 km近連續(xù)過(guò)渡流區(qū)N-S/DSMC耦合算法和80～110 km高稀薄流區(qū)DSMC模擬值定量化比較，可看出兩類(lèi)結(jié)果吻合很好、最大偏差不超過(guò)4.5%，在連續(xù)流區(qū)GKUA與N-S偏差0.01%～4.5%、在過(guò)渡流區(qū)GKUA與N-S/DSMC偏差3.4%～3.5%、在高稀薄流區(qū)GKUA與DSMC偏差0.2%～4.4%，證實(shí)所發(fā)展不同流域計(jì)算方法模型的正確可靠性與統(tǒng)一算法具有全流域計(jì)算一致收斂性?xún)?yōu)勢(shì)[83]。圖21繪出該球體隕落從高稀薄流H=105 km至過(guò)渡流區(qū)H=85 km以至連續(xù)流區(qū)H=30 km不同流域繞流面貌，直觀揭示了氣體繞流從高稀薄近自由分子流的物面附近強(qiáng)擾動(dòng)演變?yōu)橄”〗B續(xù)過(guò)渡流區(qū)的激波過(guò)渡帶，再到連續(xù)流區(qū)明晰清脆的脫體激波變化規(guī)律。

圖21 落球繞流馬赫數(shù)等值線云圖分布, H=105～30 kmFig.21 Mach No. contours around re-entry sphere with H=105～30 km computed by GKUA

將求解Boltzmann模型方程統(tǒng)一算法應(yīng)用于返回艙再入稀薄過(guò)渡流區(qū)配平迎角計(jì)算與實(shí)驗(yàn)驗(yàn)證，擬定相應(yīng)于如下條件的繞流狀態(tài):飛行高度H=88.34 km，Kn∞=0.0063，Ma∞=15.587，Re∞=3729.15，Tw/T0=0.5435，迎角α=15°、18°、20°、22°、26°、30°，在位置空間網(wǎng)格設(shè)置45×16×25，使用4096CPU進(jìn)行并行計(jì)算，圖22(a)繪出α=26°上述返回艙繞流軸對(duì)稱(chēng)面馬赫數(shù)等值線計(jì)算結(jié)果，圖22(b)給出返回艙繞重心俯仰力矩系數(shù)隨α變化關(guān)系，顯示出GKUA計(jì)算返回艙H=88.34 km繞流配平迎角αT,Cal=25.06°，這與高超聲速低密度風(fēng)洞實(shí)驗(yàn)測(cè)得的配平迎角αT,Exp=25.39°相當(dāng)一致，偏差僅1.3%，證實(shí)從介觀層次Boltzmann方程可計(jì)算建模發(fā)展氣體動(dòng)理論統(tǒng)一算法用于計(jì)算返回艙再入氣動(dòng)力/熱問(wèn)題可靠性與置信度，特別是高分辨率捕捉微弱流動(dòng)量模擬能力。這為進(jìn)一步發(fā)展各流域復(fù)雜高超聲速非平衡繞流問(wèn)題統(tǒng)一算法解算器，開(kāi)展新一代飛船返回艙再入各流域真實(shí)氣體繞流問(wèn)題計(jì)算，特別是近連續(xù)、滑移、過(guò)渡流區(qū)氣動(dòng)力/熱與姿態(tài)配平繞流問(wèn)題一體化模擬研究奠定了基礎(chǔ)。

作為跨流域空氣動(dòng)力學(xué)模擬方法集成升華對(duì)接工程需求與應(yīng)用，構(gòu)建了Boltzmann方程可計(jì)算建模氣體動(dòng)理論統(tǒng)一算法、DSMC、N-S/DSMC、滑移N-S解算器、低密度風(fēng)洞實(shí)驗(yàn)測(cè)試多種空氣動(dòng)力學(xué)模擬手段相互驗(yàn)證結(jié)合，返回艙再入從外層空間自由分子流到近地面連續(xù)流跨流域空氣動(dòng)力學(xué)一體化模擬平臺(tái)。以此為基礎(chǔ)，建立跨流域氣動(dòng)力、氣動(dòng)熱、彈道聯(lián)合計(jì)算分析機(jī)制與數(shù)據(jù)庫(kù)、圖書(shū)館式基礎(chǔ)研究結(jié)果共享系統(tǒng)，研制形成返回艙(器)從高稀薄自由分子流到近地面連續(xù)流全飛行流域沿彈道再入姿態(tài)配平繞流一體化仿真與可視化軟件系統(tǒng)，具備三大職能：1) 計(jì)算執(zhí)行功能；2) 演示分析功能；3) 存儲(chǔ)檢索“圖書(shū)館”功能。其中，計(jì)算飛行過(guò)程使用云圖粒子流優(yōu)化3D演示渲染效果；通過(guò)Fortran算法庫(kù)接入，實(shí)時(shí)計(jì)算執(zhí)行功能實(shí)現(xiàn)，優(yōu)化計(jì)算速度等。圖23繪出飛船返回艙再入跨流域氣動(dòng)模擬軟件系統(tǒng)主界面實(shí)時(shí)計(jì)算操作臺(tái)與利用該平臺(tái)計(jì)算得到返回艙再入120～70 km流場(chǎng)壓力分布。

表2 球體再入H=110～30 km繞流阻力系數(shù)GKUA與DSMC、N-S/DSMC、N-S計(jì)算比較Table 2 Comparison of drag coefficients of GKUA, DSMC, N-S/DSMC and N-S solvers for re-entry sphere with H=110～30 km

圖22 GKUA計(jì)算返回艙再入Kn∞=0.0063, Ma∞=15.6流場(chǎng)馬赫數(shù)與俯仰力矩系數(shù)隨迎角變化Fig.22 Mach contours around re-entry capsule with α=26° and pitching moment coefficients vs. α computed by the GKUA

圖23 返回艙再入跨流域氣動(dòng)模擬軟件系統(tǒng)實(shí)時(shí)計(jì)算與姿態(tài)配平繞流壓力分布120～70 kmFig.23 Real-time computing pressure distribution during the spacecraft capsule re-entry with 120～70 km

4 結(jié) 論

本文針對(duì)可回收類(lèi)航天器，如飛船返回艙再入過(guò)程所遇跨流域多尺度非平衡變化過(guò)程繞流問(wèn)題，綜述了基于Boltzman方程碰撞積分物理分析與可計(jì)算建模，發(fā)展分別考慮完全氣體、轉(zhuǎn)動(dòng)非平衡效應(yīng)、含振動(dòng)能激發(fā)熱力學(xué)非平衡效應(yīng)Boltzmann模型方程統(tǒng)一算法大規(guī)模并行計(jì)算研究進(jìn)展與算法檢驗(yàn)，證實(shí)統(tǒng)一算法在高稀薄流區(qū)，與DSMC模擬值吻合很好；在連續(xù)流區(qū)，與(滑移)N-S解算器得到的結(jié)果相一致；在中間過(guò)渡帶，與N-S/DSMC耦合算法預(yù)測(cè)值相容；統(tǒng)一算法具有全飛行流域計(jì)算一致收斂性?xún)?yōu)勢(shì)。

為克服直接對(duì)Boltzmann型速度分布函數(shù)方程在位置空間與速度空間離散數(shù)值求解計(jì)算量較大的問(wèn)題，可發(fā)展適于三維復(fù)雜外形航天器極高超聲速再入稀薄電離熱化學(xué)非平衡流動(dòng)過(guò)程DSMC方法；發(fā)展考慮振動(dòng)非平衡、黏性干擾與稀薄氣體物面滑移效應(yīng)適于航天器再入近連續(xù)滑移流區(qū)高超聲速非平衡繞流物面空氣動(dòng)力特性N-S方程解算器及其與DSMC耦合模擬技術(shù)，使之成為返回艙再入近連續(xù)過(guò)渡流區(qū)空氣動(dòng)力特性工程預(yù)測(cè)分析手段，研制基于小尺寸小載荷內(nèi)式應(yīng)變天平的返回艙再入稀薄過(guò)渡流區(qū)氣動(dòng)力測(cè)試技術(shù)與實(shí)驗(yàn)方法。建立基于GKUA、DSMC、N-S/DSMC、滑移N-S解算器、低密度風(fēng)洞實(shí)驗(yàn)驗(yàn)證結(jié)合互為補(bǔ)充跨流域空氣動(dòng)力學(xué)模擬方法綜合分析應(yīng)用研究平臺(tái)。將此平臺(tái)用于再入H=110～30 km各流域三維球體、高超聲速尖前緣中空柱裙、返回艙稀薄過(guò)渡流以至近連續(xù)流區(qū)氣動(dòng)力/熱與姿態(tài)配平繞流問(wèn)題計(jì)算與實(shí)驗(yàn)?zāi)M驗(yàn)證，一方面證實(shí)所設(shè)計(jì)跨流域空氣動(dòng)力學(xué)不同模型算法實(shí)現(xiàn)的正確可靠性，結(jié)合軟件工程、數(shù)據(jù)庫(kù)，初步研制返回艙再入跨流域氣動(dòng)力/熱沿彈道聯(lián)合計(jì)算分析機(jī)制與模擬軟件系統(tǒng)。另一方面，研究體會(huì)到，從介觀層次Boltzmann方程可計(jì)算建模發(fā)展氣體動(dòng)理論統(tǒng)一算法用于計(jì)算返回艙再入氣動(dòng)力/熱問(wèn)題具有全流域計(jì)算收斂性，直接模擬客觀物理過(guò)程的DSMC方法對(duì)返回艙再入稀薄熱化學(xué)非平衡流動(dòng)具有不可替代優(yōu)越性；N-S/DSMC雜交耦合算法、滑移N-S解算器可用于較低來(lái)流迎角近連續(xù)過(guò)渡流區(qū)氣動(dòng)特征預(yù)測(cè)分析手段；結(jié)合低密度風(fēng)洞實(shí)驗(yàn)，如何將此平臺(tái)應(yīng)用于返回(艙)器以第一、二宇宙速度再入跨流域熱化學(xué)非平衡氣動(dòng)力熱問(wèn)題模擬研究，是未來(lái)進(jìn)一步研究發(fā)展方向。

致謝:本工作得到973、杰青課題組彭傲平等同志的大力幫助，部分計(jì)算在國(guó)家超級(jí)計(jì)算無(wú)錫中心、天津中心、濟(jì)南中心進(jìn)行，特此感謝！