毛宏霞， 賈居紅， 傅德彬， 姜毅

(1.北京理工大學 宇航學院， 北京 100081； 2.北京環(huán)境特性研究所 光學輻射重點實驗室, 北京 100854)

0　引言

臨近空間高超聲速飛行器橫跨連續(xù)流區(qū)、滑移流區(qū)、過渡流區(qū)和自由分子流區(qū)，受黏性效應、高溫真實氣體效應等多種因素影響，存在層流與湍流轉(zhuǎn)捩、激波與邊界層干擾等復雜現(xiàn)象，需要對高超聲速飛行器氣動力、熱問題進行深入研究[1-2]。

裙體結(jié)構(gòu)是高超聲速飛行器后部結(jié)構(gòu)的重要形式之一。在高超聲速飛行時，氣流經(jīng)過裙體拐角受到壓縮產(chǎn)生的激波和邊界層相互干擾，不僅改變流場狀態(tài)，也對裙體結(jié)構(gòu)氣動加熱產(chǎn)生顯著影響。激波與邊界層干擾由于涉及激波和邊界層兩個原本就非常復雜的問題，因此在相應的機理研究或工程應用中，通常采用試驗或數(shù)值方法，結(jié)合雷諾數(shù)、馬赫數(shù)等流動條件對激波與邊界層流動結(jié)構(gòu)形態(tài)和氣動加熱特性等進行研究[3-5]。高超聲速國際飛行研究試驗(HIFiRE)是近年來澳大利亞國防科技集團(DSTO)、美國空軍研究所(AFRL)聯(lián)合昆士蘭大學和波音公司開展的系列高超聲速飛行試驗項目，主要采用探空火箭發(fā)射進行飛行器的上升和再入試驗，旨在通過低成本系列飛行試驗進一步探索高超聲速飛行器基礎(chǔ)問題和物理現(xiàn)象，為下一代空天飛行器奠定技術(shù)基礎(chǔ)[6]。第1項HIFIRE-1試驗主要研究邊界層轉(zhuǎn)捩和激波與邊界層干擾下的高超聲速氣動熱問題，2010年3月進行了飛行試驗，并獲得上升和再入段壓強、溫度和傳熱數(shù)據(jù)[7]。飛行試驗之外，MacLean等[8]和Wadhams等[9]先期進行了風洞試驗，并研究了氣動熱力學及邊界層激流片尺度的影響，結(jié)果表明裙體附近的流動分離與再附狀態(tài)會引起再附區(qū)氣動熱的顯著增加。Kimmel等[10]公布了HIFiRE-1飛行測試得到的基礎(chǔ)氣動數(shù)據(jù)。近年來，圍繞HIFiRE-1試驗模型對應的鈍頭體邊界層轉(zhuǎn)捩、柱裙模型激波與邊界層干擾等復雜流動問題，一些科研人員利用風洞試驗和數(shù)值計算開展了更為深入的研究。如Stanfield等[11]對邊界層轉(zhuǎn)捩進行風洞試驗，明確了不同雷諾數(shù)下的頭錐邊界層轉(zhuǎn)捩的位置和狀態(tài)，Yentsch等[12]采用不同湍流模型模擬了氣動力和氣動熱，表明k-ω模型模擬分離區(qū)過大，熱流率偏低。這些研究為認識帶裙體結(jié)構(gòu)的激波與邊界層干擾氣動熱問題提供了良好基礎(chǔ)，但對于來流參數(shù)對激波與邊界層干擾作用下?lián)Q熱能力、結(jié)構(gòu)參數(shù)對流動狀態(tài)以及氣動加熱的影響等，還有待進一步地分析和研究。

為更深入地研究高超聲速來流條件下帶裙體結(jié)構(gòu)表面的激波與邊界層干擾現(xiàn)象以及由此引起的氣動熱效應，本文采用數(shù)值計算方法，在相關(guān)文獻結(jié)合來流條件對HIFiRE-1試驗模型表面熱流密度大小進行研究的基礎(chǔ)上，重點從流動形態(tài)、斯坦頓數(shù)等角度出發(fā)，對激波與邊界層干擾條件下的熱流狀態(tài)和換熱能力進行研究；與此同時，同樣結(jié)合流動狀態(tài)，研究飛行器裙體張角、裙體長度等結(jié)構(gòu)參數(shù)對氣動加熱的影響，并分析其機理。

1　理論模型與計算方法

1.1　理論模型

本文采用三維可壓縮Navier-Stokes方程作為氣體流動的基本方程，其基本形式可表示為

(1)

采用雷諾平均(RANS)方法對(1)式進行求解計算，并采用基于渦黏性假設的SST湍流模型和Realizablek-ε湍流模型[12-13]對模型方程進行封閉處理。其中，Realizablek-ε模型主要求解湍動能及其耗散項，并在耗散項中添加源項以改善逆壓梯度流動模擬能力，SST模型在近壁面保留k-ω模型以捕捉到黏性底層的流動，在遠離壁面的區(qū)域使用k-ε模型以彌補k-ω模型對入口湍流參數(shù)過于敏感的問題，具體表達式及參數(shù)見文獻[13]。

封閉方程組采用有限體積法進行離散，黏性通量項用中心差分格式進行離散，無黏通量采用HLLC格式[14]進行離散，時間推進采用隱式方法。

1.2　計算模型

HIFiRE-1試驗模型主要有鼻錐、錐體、柱體、裙體、延長體等部分，本文計算模型為簡化后的錐- 柱- 裙結(jié)合體，具體尺寸見圖1，其中球頭半徑為2.5 mm，錐體半錐角為7°，裙體張角為33°，裙體長度L為103.7 mm，其余尺寸如圖1所示。模型驗證采用文獻[15]中的風洞試驗數(shù)據(jù)，該試驗來流條件為：壓強4.62 kPa，溫度231.7 K，速度2 182.4 m/s，馬赫數(shù)7.16.

考慮到模型的對稱性，本文僅對上半部分模型進行計算，幾何模型以鼻錐前頂點為坐標原點。采用商用軟件ICEM-CFD劃分三維結(jié)構(gòu)化網(wǎng)格，軸向、法向和周向節(jié)點數(shù)為311×101×108，在邊界層內(nèi)布置60層網(wǎng)格，網(wǎng)格總數(shù)330萬，全局和局部網(wǎng)格形式如圖2所示。采用商用計算流體力學軟件CFD++，并使用MPICH2軟件分塊進行并行處理。計算模型的入口采用全參數(shù)超聲速來流入口，下游出口采用超聲速出口，壁面采用等溫無滑移壁面條件，壁面溫度300 K.

1.3　網(wǎng)格無關(guān)性驗證

高超聲速來流通過弓形激波后被強烈壓縮，加之黏性滯止作用，來流速度在壁面處降為0. 根據(jù)能量守恒可知，外界來流巨大的動能轉(zhuǎn)化為熱能，必將在駐點區(qū)產(chǎn)生高溫，這個高溫區(qū)域在冷壁面必然產(chǎn)生很大的溫度梯度。氣動熱數(shù)值模擬對網(wǎng)格異常敏感，壁面網(wǎng)格尺度稍有改變，熱流密度誤差就可能成倍增大[16]，因此必須嚴格限制壁面法向網(wǎng)格尺度。文獻[17-18]研究表明，當壁面網(wǎng)格雷諾數(shù)接近1時，計算較為準確。

依據(jù)以上相關(guān)研究提出的壁面網(wǎng)格條件，初步設定第1層網(wǎng)格高度為5×10-7m. 為了驗證網(wǎng)格準確性，結(jié)合風洞試驗數(shù)據(jù)，采用3套網(wǎng)格方案進行試算，壁面法向網(wǎng)格高度Δn從5.7×10-7m依次遞減至5.7×10-8m，網(wǎng)格延伸率保持在1.05，壁面y+從1.954減小到0.226，如表1所示，表中p為壓強，Q為熱流密度，Re為壁面網(wǎng)格雷諾數(shù)。

表1　網(wǎng)格無關(guān)性驗證Tab.1　Grid independence

模擬得到的激波與邊界層干擾流區(qū)壁面壓強和熱流分布如圖3所示。本文試驗數(shù)據(jù)均來源于Wadhams等[9]的風洞試驗。從圖3中可以看出，對于網(wǎng)格方案2和網(wǎng)格方案3，分離區(qū)內(nèi)法向和軸向網(wǎng)格都足夠密集，壁面最大y+小于0.5，計算分析獲得的駐點及壁面的熱流和壓強吻合良好，并與試驗數(shù)據(jù)誤差較小，可以認為網(wǎng)格已滿足法向無關(guān)性要求。需要補充說明的是，軸向網(wǎng)格對錐面/圓柱過度區(qū)域及拐角區(qū)域的流動和熱流狀態(tài)同樣具有較大影響，研究表明對于本文所采用的模型結(jié)構(gòu)，當拐角或過渡區(qū)域軸向網(wǎng)格尺寸小于0.3 mm、網(wǎng)格縮放比例為1.05、周向網(wǎng)格布置間隔為2.4°時，能夠保證軸向和周向的網(wǎng)格無關(guān)性條件。

1.4　湍流模型校驗

除網(wǎng)格因素之外，氣動熱模擬還必須考慮湍流模型的影響。MacLean等[8]計算表明SA模型預測分離區(qū)過??；Yentsch等[12]對HIFiRE-1試驗模型上升段氣動熱的模擬表明k-ω模型模擬分離區(qū)過大，熱流率偏低。因此，本文對這兩類湍流模型不再進行分析，重點對應用較為廣泛的兩方程k-ε模型和SST模型進行對比驗證。

圖4給出了兩種湍流模型對應的壁面熱流和壓強分布。從圖4中可以看出，兩種湍流模型的模擬分離泡大小相當、位置基本相同，k-ε模型模擬平臺區(qū)壓強稍大，再附區(qū)壁面壓強稍低于SST模型，誤差在3%以內(nèi)。但從圖4(a)可以明顯看出二者模擬得到的壁面熱流差別，SST模型高估了再附點熱流，峰值誤差9.8%，而k-ε模型完全低估了分離區(qū)熱流密度，最大誤差-27.9%. 由此可以看出，湍流模型對激波與邊界層干擾氣動熱模擬準確性也有較大影響。考慮到SST模型對再附點熱流峰值預估誤差小于10%，本文后續(xù)使用SST模型進行計算和分析。

2　結(jié)果與分析

2.1　激波與邊界層干擾氣動熱

結(jié)合前述風洞試驗對應的計算工況，從模擬得到的壁面熱流和壓強分布(見圖3和圖4)可以看出，受柱裙拐角氣流壓縮影響，氣流邊界層在1.55 m處開始分離，此時壁面熱流和壓強上升并保持一段平臺期，而后在1.63 m處開始突增，并于1.66 m處再附于裙體斜面，此時熱流密度和壓強均迅速增加并出現(xiàn)峰值，二者最大值分別為388.4 W/cm2和135.4 kPa. 而后在扇形壓縮波區(qū)域，熱流密度和壓強快速降低并趨于常數(shù)，氣流流過延長體部分時迅速膨脹，熱流密度和壓強都迅速降低。對比文獻[9]試驗數(shù)據(jù)可以看出，計算獲得的分離位置、再附位置均與試驗狀態(tài)吻合良好。圖5(a)為數(shù)值模擬獲得的馬赫數(shù)Ma分布云圖，顯示了拐角處流場特性和激波結(jié)構(gòu)。由圖5(a)可知，來流由左向右，上游流動分離誘導出一個明顯的分離激波，向上偏折后再附到裙體壁面，形成再附激波。拐角處分離激波帶來的大逆壓梯度使壁面邊界層發(fā)生較為嚴重的變形，并向上游和下游傳播，造成拐角附近邊界層膨脹，形成氣泡狀分離區(qū)，并且分離區(qū)內(nèi)為亞聲速流動。再附激波沿裙體向下游發(fā)展并遠離壁面，形成扇形壓縮波，扇形波流經(jīng)裙體延長部分時氣流迅速膨脹，出現(xiàn)明顯的膨脹波。圖5(b)給出了流場壓強分布，圖5(c)給出了流場溫度和壁面熱流分布。由圖5(b)和圖5(c)可知，在分離區(qū)內(nèi)為高溫區(qū)，邊界層再附后出現(xiàn)環(huán)形高熱流帶。此外，在裙體前沿的局部區(qū)域，受邊界層分離引起的流動振蕩現(xiàn)象影響，壓強和熱流呈現(xiàn)周向非均勻分布。本文采用時間推進法進行迭代求解，能夠在一定程度上捕獲這類局部流動振蕩現(xiàn)象。

2.2　來流條件對壁面熱流的影響分析

為了研究來流條件對激波與邊界層干擾氣動熱的影響，選擇不同來流雷諾數(shù)、馬赫數(shù)的工況進行計算分析。計算模型的基本狀態(tài)如表2所示。

表2　來流參數(shù)與裙體尺寸Tab.2　Free stream parameters and structure size

2.2.1來流雷諾數(shù)影響

首先考慮來流馬赫數(shù)、溫度等條件保持不變的前提下，雷諾數(shù)從1.738×107降低到3.280×106和1.640×106，計算獲得的分離區(qū)熱流及壓強分布狀態(tài)。

圖6(a)給出了不同雷諾數(shù)下的壁面熱流分布。從圖6中可以看出，模擬得到的流場結(jié)構(gòu)、分離區(qū)位置基本保持不變，分離點仍保持在1.55 m位置附近，但壁面壓強和熱流發(fā)生了較大變化。當雷諾數(shù)降低到3.28×106時，壁面熱流得到顯著降低，再附點熱流峰值由394.4 W/cm2降低到185.1 W/cm2；當雷諾數(shù)降低到1.64×106時，壁面熱流繼續(xù)變化，再附點熱流峰值降低到108.9 W/cm2. 事實上，隨著雷諾數(shù)的降低，來流密度降低，熱流密度自然也會降低，表明計算結(jié)果與理論預計是相符的。圖6(b)給出了壓強分布情況，從中可以看出，壁面壓強同樣持續(xù)變化，峰值由132.7 kPa分別降低到48.3 kPa和24.8 kPa，平臺區(qū)壓強由10 kPa降低到4 kPa和2 kPa.

2.2.2來流馬赫數(shù)影響

在其余參數(shù)不變情況下，馬赫數(shù)由7.16降低至5.00和3.00，壁面熱流和壓強結(jié)果如圖7(a)和圖7(b)所示。從圖7(a)和圖7(b)中可以看出：馬赫數(shù)降低后，壁面熱流密度大幅度降低；馬赫數(shù)為5.00時熱流峰值為101.2 W/cm2；馬赫數(shù)為3.00時壁面熱流基本保持平緩發(fā)展，峰值熱流密度僅為12.8 W/cm2；壁面壓強也大幅度降低，峰值分別為69.7 kPa和25.8 kPa.

此時，流場結(jié)構(gòu)如圖7(c)和圖7(d)所示。從圖7(c)和圖7(d)可以看到：馬赫數(shù)降低后激波厚度增大，分離泡范圍增加；流場波系結(jié)構(gòu)較為清晰，分離激波向下游發(fā)展，與再附激波相交，形成三叉激波，這也是再附點壓強和熱流急劇上升的主要原因；馬赫數(shù)為3.00時三叉激波距離裙體壁面的位置明顯遠于馬赫數(shù)為5.00時的結(jié)果，這說明來流馬赫數(shù)對再附點熱流大小有重要影響。

2.2.3來流條件影響壁面熱流的綜合分析

為進一步明確來流條件對壁面熱流的影響，下面分析不同來流條件下表征換熱能力的斯坦頓數(shù)St變化情況。St采用來流密度、速度和溫度作為參考值。

圖8給出了不同雷諾數(shù)和不同馬赫數(shù)條件對應的St數(shù)分布情況。從圖8(a)可以看出：在不同來流條件下，受裙體上氣流再附效應影響，裙體再附區(qū)域的St數(shù)均顯著增加，表明換熱能力在這一區(qū)域明顯增強；對于不同來流條件，隨著雷諾數(shù)的增加，氣流再附效應引起的St數(shù)增量幅值減小，表明來流雷諾數(shù)增加雖然使裙體上最大熱流密度絕對值增加，但表征換熱能力的St數(shù)反而降低。從圖8(b)可以看出，隨著來流馬赫數(shù)(對應來流速度)的增加，裙體上St數(shù)顯著增加，表明氣流再附引起的換熱能力增加。結(jié)合具體研究工況，可明確：隨著來流密度的增加，表征換熱能力的St數(shù)減小；隨著來流速度的增加，表征換熱能力的St數(shù)增加。

2.3　裙體結(jié)構(gòu)對壁面熱流的影響分析

2.3.1裙體張角影響

結(jié)合表2所示來流條件，選擇裙體張角從33°變化為27°、30°、37°進行計算，獲得壁面熱流分布與壓強分布如圖9(a)和圖9(b)所示，馬赫數(shù)分布如圖9(c)和圖9(d)所示。由圖9可以看出：當張角縮小至30°時，邊界層分離區(qū)變小，上游分離點和下游再附點都有回縮趨勢，分離點位于1.58 m，再附點位于1.66 m，壁面熱流和壓強相應降低，峰值分別降到318 W/cm2和104.1 kPa；當裙體張角進一步縮小至27°時，分離泡縮小至即將消失，來流在拐角處基本直接沿壁面向上流動，受外側(cè)來流影響，同樣形成扇形壓縮波；當增大張角至37°時，激波與邊界層干擾分離區(qū)增大，分離點向上游發(fā)展至1.52 m，再附點向下游發(fā)展至1.67 m，裙體表面熱流密度和壓強相應上升，再附點熱流和壓強峰值升高到481.7 W/cm2和172 kPa.

2.3.2裙體長度影響

在流動參數(shù)保持不變情況下，通過控制改變?nèi)贵w長度L，計算后得到的壁面熱流和壓強分布如圖10(a)和圖10(b)所示，馬赫數(shù)分布如圖10(c)和圖10(d)所示。由圖10可見，當裙體長度從103.7 mm縮短到82.0 mm時，流動依然完成在裙體壁面的再附，壁面熱流密度和壓強分布基本保持與原模型一致。為了探究流動無法完成再附情況下的壁面熱流變化，進一步將裙體長度縮短至32.0 mm位置，此時流動依然被壓縮產(chǎn)生邊界層分離區(qū)，但剪切層無法在裙體斜面完成再附，表現(xiàn)為分離泡縮小，壁面熱流和壓強無法達到峰值即開始迅速降低。從圖5中流場馬赫數(shù)分布可以看出，分離激波的發(fā)展會受到自由來流的影響，為了驗證以上判斷，把裙體長度增大至202.0 mm，此時飛行器錐體外部來流與裙體上方分離激波相交，分離激波受壓之后沿壁面向下游發(fā)展，直接導致扇形壓縮波內(nèi)壓強增大，進而引發(fā)邊界層分離泡增大，分離點前移，再附點后移。此時，模擬得到的壁面熱流和壁面壓強出現(xiàn)明顯上升趨勢，熱流峰值464.1 W/cm2，壓強峰值159.53 kPa.

3　結(jié)論

本文對HIFiRE-1高超聲速飛行器激波與邊界層干擾氣動熱問題及影響因素進行了研究，得到以下結(jié)論：

1)當來流雷諾數(shù)改變時，分離區(qū)流場結(jié)構(gòu)變化較小，但熱流密度隨著雷諾數(shù)的增加而增加。當來流馬赫數(shù)改變時，激波與邊界層干擾區(qū)的激波結(jié)構(gòu)發(fā)生相應改變，熱流密度同樣隨著馬赫數(shù)的增加而增加。此外，在來流變化過程中，裙體氣流再附區(qū)域換熱能力也會產(chǎn)生顯著變化。隨著來流速度增加，裙體上氣流再附區(qū)表征換熱能力的St數(shù)增加；而隨著來流密度增加，再附區(qū)表征換熱能力的St數(shù)減小。

2)當裙體張角改變時，分離區(qū)壁面熱流分布大幅改變，當裙體張角增大到37°時，分離區(qū)增大，分離點向上游發(fā)展，再附點向下游發(fā)展；當裙體張角減小到27°時，分離區(qū)基本消失，壁面熱流將近降到一半。裙體長度增大會導致分離泡增大、熱流增加，裙體長度減小至分離流無法完成再附時，壁面熱流無法達到峰值并迅速降低。

參考文獻(References)

[1]Persovaa M G, Soloveichika Y G, Belovb V K. Modeling of aerodynamic heat flux and thermoelastic behavior of nose caps of hypersonic vehicles [J].Acta Astronautica, 2017, 136(1): 312-331.

[2]張志豪, 孫得川.飛行器氣動加熱燒蝕工程計算[J]. 兵工學報, 2015, 36(10): 1949-1954.

ZHANG Zhi-hao, SUN De-chuan. Calculation of aerodynamic heating and ablation of multi-layer thermal protection material[J]. Acta Armamentarii, 2015, 36(10): 1949-1954. (in Chinese)

[3]童福林, 李新亮, 唐志共, 等. 轉(zhuǎn)捩對壓縮拐角激波/邊界層干擾分離泡的影響[J]. 航空學報， 2016, 37(10): 2909-2912.

TONG Fu-lin, LI Xin-liang, TANG Zhi-gong, et al. Transition effect on separation bubble of shock wave/boundary layer interaction in a compression ramp[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(10): 2909-2912. (in Chinese)

[4]宋友富, 徐晶磊, 張楊, 等. 壓縮拐角激波/邊界層干擾的可壓縮湍流模型研究[J]. 推進技術(shù), 2017, 38(2): 281-288.

SONG Yong-fu, XU Jing-lei, ZHANG Yang, et al. Research of compressible trubulence model in shock wave/boundary-layer interaction flow at a compression corner[J]. Journal of Propulsion Technology, 2017, 38(2): 281-288. (in Chinese)

[5]董祥瑞, 陳耀慧, 董剛, 等. 基于雙微楔的高超聲速激波與邊界層干擾控制研究[J]. 兵工學報, 2016, 37(9): 1624-1632.

DONG Xiang-rui, CHEN Yao-hui, DONG Gang, et al. Research on control of hypersonic shock wave/boundary layer interactions by double micro-ramps[J]. Acta Armamentarii, 2016, 37(9): 1624-1632. (in Chinese)

[6]Bisek N J. High-fidelity simulations of the HIFiRE-6 flow path, AIAA-2016-1115[R]. Reston, VA, US: AIAA, 2016.

[7]Kimmel R L, Adamczak D W. HIFiRE-1 flight trajectory estimation and initial experimental results, AIAA-2011-2358 [R]. Reston, VA, US: AIAA, 2011.

[8]MacLean M, Wadhams T, Holden M, et al. Ground test studies of the HIFiRE-1 transition experiment—part 2: computational analysis[J]. Journal of Spacecraft and Rockets, 2008, 45(6): 1149-1164.

[9]Wadhams T, Mundy E, MacLean M, et al. Pre-flight ground testing of the full-scale HIFiRE-1vehicle at fully duplicated flight conditions: part II, AIAA-2008-0639[R]. Reston, VA, US: AIAA, 2008.

[10]Kimmel R L, Adamczak D W. HIFiRE-1 preliminary aerothermodynamics measurements, AIAA-2011-3413[R]. Reston, VA, US: AIAA, 2011.

[11]Stanfield S A, Kimmel R L, Adamczak D W, et al. Boundary layer transition experiment during reentry of HIFiRE-1[J]. Journal of Spacecraft and Rockets, 2015, 52(3): 637-649.

[12]Yentsch R J, Gaitonde D V, Kimmel R L. Performance of turbulence modeling in simulation of the HIFiRE-1 flight test[J]. Journal of Spacecraft and Rockets, 2014, 51(1): 117-127.

[13]Rautaheimo P R. Developments in turbulence modeling with Reynolds-averaged Navie Stokes equations[M]. Helsinki, Finland: Finnish Academies Technology Press, 2001: 17-19.

[14]Xu L, Wu Q J, Weng P F. HLLC Riemann solver based on high-order reconstruction for unsteady inviscid compressible flows[C]∥Proceedings of IEEE International Conference on Computer Science and Automation Engineering. Shanghai, China: IEEE, 2011: 618-622.

[15]Holden M, Wadhams T, MacLean M, et al. Experimental studies of shock wave/turbulent boundary layer interaction in high Reynolds number supersonic and hypersonic flows to evaluate the performance of CFD codes, AIAA-2010-4468[R]. Reston, VA, US: AIAA, 2010.

[16]Yang J L, Liu M. A wall grid scale criterion for hypersonic aerodynamic heating calculation[J]. Acta Astronautica, 2017, 136(1): 137-143.

[17]張智超, 高振勛, 蔣崇文, 等. 高超聲速氣動熱數(shù)值計算壁面網(wǎng)格準則[J]. 北京航空航天大學學報, 2015, 41(4): 594-600.

ZHANG Zhi-chao, GAO Zhen-xun, JIANG Chong-wen, et al. Grid generation criterions in hypersonic aeroheating computations[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(4): 594-600. (in Chinese)

[18]Sohag F A, Mohanta L, Cheung F B. CFD analyses of mixed and forced convection in a heated vertical rod bundle[J]. Applied Thermal Engineering, 2017, 117(4): 85-93.