網(wǎng)格對(duì)高超聲速鈍頭體表面熱流數(shù)值模擬結(jié)果的影響

馬崇立，劉景源

南昌航空大學(xué) 飛行器工程學(xué)院，南昌 330063

發(fā)展具有一定高超聲速巡航能力、可重復(fù)使用的新型天地往返運(yùn)輸系統(tǒng)及在大氣層內(nèi)飛行的高超聲速飛行器已成為航天航空強(qiáng)國(guó)技術(shù)儲(chǔ)備中需要研究的關(guān)鍵課題，其中高超聲速?gòu)?fù)雜外形飛行器表面熱流精確預(yù)估是其中的一個(gè)核心問(wèn)題。

對(duì)高超聲速飛行器表面熱流精確預(yù)估有試驗(yàn)、理論分析及數(shù)值模擬等3 種方法。但是無(wú)論高超聲速飛行器的地面試驗(yàn)還是飛行試驗(yàn)，均存在成本高、周期長(zhǎng)等問(wèn)題。并且地面試驗(yàn)條件的限制無(wú)法模擬真實(shí)飛行環(huán)境；而飛行試驗(yàn)雖然可以模擬飛行環(huán)境，但有投入巨大、可重復(fù)性弱等缺點(diǎn)。由于高超聲速流動(dòng)方程的非線性性質(zhì)以及高超聲速飛行器流場(chǎng)的復(fù)雜性，理論分析無(wú)法給出可用于指導(dǎo)工程設(shè)計(jì)的結(jié)果。因此，數(shù)值模擬技術(shù)成為高超聲速飛行器表面熱流精確估計(jì)的重要手段。

高超聲速飛行器頭部激波強(qiáng)、流動(dòng)減速劇烈，氣動(dòng)加熱嚴(yán)重。為了有效防熱，高超聲速飛行器采用鈍頭型。因此，數(shù)值模擬精確預(yù)估高超聲速鈍頭體飛行器表面熱流即成為高超聲速?gòu)?fù)雜外形表面熱流必須先行解決的問(wèn)題。

數(shù)值模擬預(yù)估高超聲速鈍頭體的精度受模擬所用網(wǎng)格、選取的離散方法、邊界條件及數(shù)值格式等的影響［1-2］。而鈍體表面法向第1 層網(wǎng)格尺度被認(rèn)為是重要的影響因素［3］。文獻(xiàn)［4-5］給出了基于壁面參數(shù)的壁面法向網(wǎng)格雷諾數(shù)，認(rèn)為壁面法向網(wǎng)格雷諾數(shù)對(duì)氣動(dòng)熱的預(yù)測(cè)具有重要的影響，并經(jīng)數(shù)值分析指出，為了精確預(yù)估壁面熱流，壁面法向網(wǎng)格雷諾數(shù)必須小于一定的數(shù)值。Klopfer 等［6］提出了基于來(lái)流參數(shù)的法向網(wǎng)格雷諾數(shù)，數(shù)值模擬結(jié)果表明，為了精確預(yù)估高超聲速鈍體表面熱流，基于來(lái)流參數(shù)的法向網(wǎng)格雷諾數(shù)亦必須小于一定的數(shù)值。潘沙等［7］研究了基于來(lái)流參數(shù)的法向網(wǎng)格雷諾數(shù)對(duì)高超聲速鈍頭體熱流的影響，并給出了法向網(wǎng)格劃分的建議。張翔等［8］應(yīng)用數(shù)值模擬方法研究了高超聲速二維半圓柱繞流，并分析了法向網(wǎng)格雷諾數(shù)與激波位置網(wǎng)格尺度對(duì)壁面熱流的影響。閆文輝等［9］應(yīng)用數(shù)值模擬方法分析了二維半圓柱超聲速繞流的摩阻及熱流，結(jié)果表明，為準(zhǔn)確計(jì)算摩阻及熱流法向第1 層網(wǎng)格間距、網(wǎng)格長(zhǎng)寬比與擴(kuò)張度等需要取合適的值。Qu［10］系統(tǒng)地研究了幾種二階迎風(fēng)數(shù)值格式在計(jì)算高超聲速壁面熱流時(shí)，對(duì)限制器及網(wǎng)格雷諾數(shù)的依賴(lài)，結(jié)果表明，SLAU數(shù)值格式計(jì)算高超聲速熱流具有高精度及強(qiáng)魯棒性。Zhao 等［11］基于數(shù)值模擬方法研究了不同網(wǎng)格雷諾數(shù)對(duì)高超聲速化學(xué)非平衡繞流流動(dòng)的飛行器表面及駐點(diǎn)熱流的影響，結(jié)果表明，隨網(wǎng)格雷諾數(shù)減小，所計(jì)算的飛行器表面熱流及駐點(diǎn)熱流均增大并最終計(jì)算的熱流具有網(wǎng)格無(wú)關(guān)性。Qu［12］研 究 了Roe、AUSM+及AUSMPW+這3種數(shù)值格式計(jì)算高超聲速真實(shí)飛行器表面熱流的模擬精度，結(jié)果表明，為精確計(jì)算飛行器表面熱流，基于來(lái)流參數(shù)的網(wǎng)格雷諾數(shù)不應(yīng)大于20。

由于缺少理論上的支撐，上述方法確定的壁面法向網(wǎng)格尺度具有一定的經(jīng)驗(yàn)性及不確定性。程曉麗等［1］則提出將分子平均自由程作為壁面法向網(wǎng)格尺度，數(shù)值模擬結(jié)果表明，可將該尺度作為法向網(wǎng)格劃分的最優(yōu)下限。張智超等［2］則進(jìn)一步發(fā)展了程曉麗等提出的方法，改善了基于分子平均自由程準(zhǔn)則需要通過(guò)粗網(wǎng)格試算獲得駐點(diǎn)位 置 壁 面 密 度 的 缺 點(diǎn)。Yang 與Liu［13］亦 基 于 分子運(yùn)動(dòng)論，引入平均自由程給出了一種駐點(diǎn)附近法向網(wǎng)格劃分方法，并應(yīng)用所給出的方法數(shù)值模擬了高超聲速半圓柱繞流熱流流場(chǎng)，結(jié)果表明，為準(zhǔn)確計(jì)算高超聲速表面熱流，壁面法向第1 層網(wǎng)格的尺度應(yīng)為基于壁面參數(shù)的分子平均自由程的大小，并且所提出的方法不但適用于介質(zhì)為空氣的高超聲速繞流，對(duì)氮?dú)庖噙m用。Ren 等［14］根據(jù)一維流動(dòng)理論，給出了基于來(lái)流及壁面參數(shù)的法向網(wǎng)格雷諾數(shù)及基于壁面參數(shù)的分子平均自由程三者在計(jì)算高超聲速連續(xù)流及過(guò)渡區(qū)時(shí)的關(guān)系，并指出為精確計(jì)算壁面熱流，壁面法向第1 層網(wǎng)格尺度必須嚴(yán)格滿(mǎn)足一定的網(wǎng)格準(zhǔn)則。Huang 與Hu［15］基于Kemp-Riddell 經(jīng)驗(yàn)公式，提出了計(jì)算高超聲速飛行器表面熱流的壁面法向第1 層網(wǎng)格劃分的新準(zhǔn)則，并基于新準(zhǔn)則計(jì)算了球頭高超聲速繞流表面熱流，結(jié)果表明，所提出的新方法，比基于網(wǎng)格雷諾數(shù)的方法給出的駐點(diǎn)熱流更精確。張亮等［16］則給出了一種適用于高超聲速高冷壁下表面熱流法向網(wǎng)格尺度上限確定方法。

對(duì)高超聲速鈍頭體表面熱流進(jìn)行精確預(yù)估的絕大多數(shù)的研究，均關(guān)注表面法向網(wǎng)格的劃分方法，而流向網(wǎng)格尺度對(duì)氣動(dòng)熱的影響研究有所欠缺。文獻(xiàn)［17-18］雖然基于差分格式誤差匹配原則，提出了網(wǎng)格劃分依據(jù)，但該方法僅基于Navier-Stokes 方程組的流向動(dòng)量方程，忽略了對(duì)高超聲速流動(dòng)的能量方程的分析，而能量方程對(duì)高超聲速熱流預(yù)估具有重要影響。另外，該方法需要事先給定基于自由來(lái)流參數(shù)的法向網(wǎng)格雷諾數(shù)，而法向網(wǎng)格雷諾數(shù)的選擇具有一定的不確定性的缺陷。

本文基于法向網(wǎng)格劃分方法、高超聲速能量方程及動(dòng)量方程修正方程各誤差項(xiàng)的誤差匹配，提出了基于壁面參數(shù)的法向網(wǎng)格雷諾數(shù)及基于鈍頭體特征長(zhǎng)度的來(lái)流雷諾數(shù)作為網(wǎng)格劃分的兩個(gè)依據(jù)，從而給出了一種對(duì)高超聲速鈍頭體表面繞流的較為完整的網(wǎng)格劃分方法，然后應(yīng)用所提出的方法對(duì)二維半圓柱、球頭等高超聲速鈍頭體繞流進(jìn)行了數(shù)值分析及討論，并指出本文方法不足之處及進(jìn)一步的研究方向。

1 數(shù)值格式的誤差匹配

應(yīng)用數(shù)值模擬方法對(duì)高超聲速鈍頭體表面熱流進(jìn)行預(yù)估，必然與流場(chǎng)控制方程的精確解存在一定的誤差，而合理控制數(shù)值模擬的誤差對(duì)精確預(yù)估飛行器表面熱流具有重要意義。由于高超聲速飛行器表面存在邊界層，所以在流場(chǎng)數(shù)值計(jì)算時(shí)邊界層內(nèi)計(jì)算網(wǎng)格的長(zhǎng)寬比較大，而計(jì)算誤差與網(wǎng)格的長(zhǎng)寬比密切相關(guān)。保證控制方程的差分格式的修正方程各誤差項(xiàng)為同一量級(jí)［18］，是預(yù)估飛行器表面熱流的前提。

1.1 能量方程誤差匹配

二維定常無(wú)量綱高超聲速能量方程［19］的差分格式的修正方程為

式中：為密度；為比定壓熱容；、分別為流向及法向速度分量；為溫度；γ為比熱比；Ma∞為來(lái)流馬赫數(shù)；為壓強(qiáng)；Pr=0.71 為普朗特?cái)?shù)；為熱傳導(dǎo)系數(shù)；、、、分別為流向及法向與溫度相關(guān)的誤差項(xiàng)；為能量方程中壓力離散誤差項(xiàng)；為能量方程耗散函數(shù)離散誤差項(xiàng)；Re為基于飛行器特征長(zhǎng)度的自由來(lái)流雷諾數(shù)。為耗散函數(shù)，表達(dá)式為

為方便起見(jiàn)，以下分析及論述忽略上述方程符號(hào)上標(biāo)“-”，而用無(wú)上標(biāo)量表示無(wú)量綱量。

不失一般性，將無(wú)量綱的熱傳導(dǎo)系數(shù)k當(dāng)作近似常數(shù)處理，式（1）中流向與法向溫度相關(guān)項(xiàng)的誤差項(xiàng)分別為

式中：a1、a2、b1、b2、c1、c2、d1、d2為差分格式的系數(shù)。

式（1） 中RΦ=Rux+Rvy+Rvx+Ruy+Ruxvy+Rvxuy，等號(hào)右端各項(xiàng)的具體形式如下：

式中：g1、g2、h1、h2、r1、r2、s1、s2、e1、f1為差分格式的系數(shù)。

用δ表示邊界層內(nèi)流動(dòng)的特征尺度，根據(jù)邊界層理論，可得無(wú)量綱參數(shù)量級(jí)關(guān)系為［18，20］

式中：上標(biāo)p、q、m、n為各物理量求導(dǎo)階次。

1.1.1 能量方程溫度項(xiàng)誤差匹配

本節(jié)首先對(duì)能量方程式（1）中溫度相關(guān)項(xiàng)進(jìn)行分析。在數(shù)值計(jì)算中，為了能夠準(zhǔn)確捕捉到壁面附近的流動(dòng)參數(shù)，需要將邊界層內(nèi)沿法向方向進(jìn)行加密，即邊界層內(nèi)法向網(wǎng)格尺度Δy要小于δ；所以Δy δ應(yīng)為小量，從而與溫度相關(guān)的誤差項(xiàng)可簡(jiǎn)化為

為了保證數(shù)值計(jì)算的精度，并使法向與流向2 個(gè)方向所選擇的網(wǎng)格尺度最優(yōu)，各誤差項(xiàng)應(yīng)保持為同一量級(jí)［18］，即式（7）～式（10）左側(cè)誤差項(xiàng)的比值量級(jí)為1。對(duì)以上各誤差項(xiàng)進(jìn)行匹配，可得如下結(jié)論：

1）RTx1與RTy1匹配，則

2）RTx2與RTy1匹配，則

3）RTx1與RTy2匹配，則

4）RTx1與RTy2匹配，則

在數(shù)值計(jì)算中，若流向與法向差分格式采取相同的階數(shù)，即a=b=2c=2d；故網(wǎng)格長(zhǎng)寬比量級(jí)應(yīng)取式（7）～式（10）匹配結(jié)果中較小量以保證邊界層內(nèi)網(wǎng)格尺度滿(mǎn)足數(shù)值計(jì)算需求，即

1.1.2 能量方程耗散項(xiàng)誤差匹配

耗散項(xiàng)與能量方程的溫度項(xiàng)相比，對(duì)飛行表面熱流的影響亦重要，因此為了保證數(shù)值計(jì)算的精度，能量方程耗散項(xiàng)中流向與法向的誤差項(xiàng)也應(yīng)為同一量級(jí)。

根據(jù)邊界層理論以及Δy δ應(yīng)為小量，結(jié)合式（5），耗散項(xiàng)各誤差項(xiàng)可記為

同1.1.1節(jié)，式（12）各誤差項(xiàng)亦應(yīng)保持為同一量級(jí)。對(duì)以上各誤差項(xiàng)進(jìn)行匹配，可得如下結(jié)論：

1）Rux與Ruy匹配，則

2）Rvx與Rvy匹配，則

3）Rux與Rvy匹配，則

4）Rvx與Ruy匹配，則

若耗散項(xiàng)流向與法向差分格式采用相同的階數(shù)，即g=h=r=s，為保證數(shù)值計(jì)算精度，網(wǎng)格長(zhǎng)寬比量級(jí)應(yīng)取上述匹配結(jié)果中較小量以保證邊界層內(nèi)網(wǎng)格尺度滿(mǎn)足數(shù)值計(jì)算需求即匹配式（15），其長(zhǎng)寬比量級(jí)關(guān)系式同式（11）。

不失一般性，對(duì)耗散函數(shù)各導(dǎo)數(shù)項(xiàng)的離散均采用相同的差分精度，因此式（12）最后兩項(xiàng)量級(jí)相同，并且與式（12）其他4 項(xiàng)進(jìn)行誤差比較時(shí)，所得結(jié)論包含在式（13）～式（16）內(nèi)，限于篇幅，此處不再贅述。

1.2 動(dòng)量方程誤差匹配

對(duì)動(dòng)量方程的分析與文獻(xiàn)［18］類(lèi)似，二維定常無(wú)量綱高超聲速邊界層流動(dòng)的動(dòng)量方程差分格式的修正方程為

式 中：、分 別 為 流 向 與 法 向 對(duì) 流 項(xiàng) 誤 差項(xiàng)；、分別為流 向及法向 黏性項(xiàng)誤 差項(xiàng)。

同1.1 節(jié)，以下分析及論述忽略上述方程符號(hào)上標(biāo)“-”，用無(wú)上標(biāo)量表示無(wú)量綱量，則式（17）中誤差項(xiàng)分別為

式中：a*1、b*1、c*1、c*2、d*1、d*2為動(dòng)量方程差分格式的系數(shù)。

根據(jù)邊界層理論與式（5），各誤差項(xiàng)可簡(jiǎn)化為

為保證數(shù)值計(jì)算精度，對(duì)流項(xiàng)與黏性項(xiàng)誤差項(xiàng)比值量級(jí)為1。對(duì)以上各誤差項(xiàng)進(jìn)行匹配可得：

1）Rux1與Ruy1匹配，則

2）Rux2與Ruy1匹配，則

3）Rux1與Ruy2匹配，則

4）Rux2與Ruy2匹配，則

若流向與法向差分格式采用相同階數(shù)，即a*=b*=2c*=2d*，網(wǎng)格長(zhǎng)寬比量級(jí)應(yīng)取式（20）～式（23）匹配結(jié)果中較小量，即

1.3 能量方程與動(dòng)量方程誤差匹配結(jié)果

通過(guò)對(duì)能量方程與動(dòng)量方程進(jìn)行誤差分析可得，數(shù)值計(jì)算中為保證數(shù)值計(jì)算精度，網(wǎng)格長(zhǎng)寬比量級(jí)上應(yīng)取較小量以保證邊界層內(nèi)網(wǎng)格尺度滿(mǎn)足數(shù)值耗散較小的要求。分別對(duì)溫度項(xiàng)與耗散項(xiàng)誤差匹配中長(zhǎng)寬比量級(jí)較小項(xiàng)，即式（7）、式（9）、式（15）、式（20）及式（22）等的分析表明：①當(dāng)流向與法向?qū)?shù)各項(xiàng)采用相同差分精度時(shí)，網(wǎng)格長(zhǎng)寬比應(yīng)與飛行器特征長(zhǎng)度的自由來(lái)流雷諾數(shù)平方根同量級(jí)；②當(dāng)流向各導(dǎo)數(shù)采用高階差分精度時(shí)，如法向第1 層網(wǎng)格間距確定，為了滿(mǎn)足誤差匹配，則流向網(wǎng)格尺度應(yīng)增大即放寬對(duì)流向網(wǎng)格尺度的要求；③當(dāng)法向采用高階精度時(shí)，流向網(wǎng)格尺度應(yīng)降低。

當(dāng)對(duì)流項(xiàng)、溫度二階空間導(dǎo)數(shù)及耗散函數(shù)項(xiàng)等均取相同差分精度時(shí)，網(wǎng)格長(zhǎng)寬比應(yīng)與基于鈍頭體特征長(zhǎng)度的來(lái)流雷諾數(shù)的平方根同量級(jí)，即式（11）或式（24）。當(dāng)能量及動(dòng)量方程的對(duì)流項(xiàng)、溫度二階空間導(dǎo)數(shù)及耗散函數(shù)項(xiàng)等取不同差分精度時(shí)，由于網(wǎng)格長(zhǎng)寬比量級(jí)應(yīng)取上述匹配結(jié)果中較小量，網(wǎng)格長(zhǎng)寬比的結(jié)論亦為式（11）或式（24）。

對(duì)于三維問(wèn)題，由于2 個(gè)切向方向的截?cái)嗾`差表達(dá)式類(lèi)似，所以只需要將上述匹配結(jié)果中的流向網(wǎng)格尺度Δx用切向網(wǎng)格尺度hτ=max(Δx， Δz)代替即可把本節(jié)的結(jié)論推廣到三維數(shù)值模擬中。

另 外，應(yīng) 用 式（5）的 量 級(jí) 估 計(jì)，式（1）及式（17）與壓強(qiáng)相關(guān)項(xiàng)的誤差項(xiàng)進(jìn)行誤差匹配的結(jié)果表明，流向與法向網(wǎng)格尺度比為來(lái)流雷諾數(shù)。即所得結(jié)果不對(duì)本節(jié)結(jié)論產(chǎn)生影響，具體分析此處不再贅述。

2 高超聲速鈍頭體繞流流場(chǎng)網(wǎng)格劃分

第1 節(jié)的分析雖然給出了流場(chǎng)流向與法向網(wǎng)格劃分時(shí)，應(yīng)考慮流向與法向網(wǎng)格尺度的比，但數(shù)值模擬中還不能確定流向及法向等方向的網(wǎng)格尺度的具體數(shù)值。若能確定流向或者壁面法向的網(wǎng)格尺度，則高超聲速繞流流場(chǎng)的網(wǎng)格劃分尺度即可確定。

由壁面網(wǎng)格雷諾數(shù)的定義及壁面分子平均自由程的公式［20］可知

式中：aw為基于壁面參數(shù)的聲速。

由式（25）及式（26）可知

由程曉麗等［1］給出的壁面法向第1 層網(wǎng)格為壁面上分子平均自由程，因此基于壁面參數(shù)的網(wǎng)格雷諾數(shù)為1.5。

根據(jù)文獻(xiàn)［21］對(duì)流域劃分，則為了使Navier-Stokes（N-S）方程組成立，努森數(shù)Kn＜0.1。根據(jù)文獻(xiàn)［22］，上述流域劃分方法可推廣到數(shù)值模擬的壁面附近的網(wǎng)格劃分中，則對(duì)于數(shù)值計(jì)算：

由式（27）及式（28）可知

但考慮到數(shù)值計(jì)算中并非直接求解N-S 方程組，而是求解其修正方程組，因此實(shí)際計(jì)算中，應(yīng)結(jié)合數(shù)值格式的精度確定基于壁面參數(shù)的網(wǎng)格雷諾數(shù)選取的具體數(shù)值。

因此，本文提出基于壁面法向網(wǎng)格雷諾數(shù)Recell，w，與基于繞流流場(chǎng)特征長(zhǎng)度的自由來(lái)流雷諾數(shù)相結(jié)合確定繞流流場(chǎng)網(wǎng)格法向及流向網(wǎng)格尺度，用于數(shù)值模擬前網(wǎng)格劃分的參考。具體為，首先由來(lái)流及壁面參數(shù)，根據(jù)指定的Recell，w具體數(shù)值，確定法向網(wǎng)格尺度，然后再根據(jù)流向與法向網(wǎng)格尺度比確定流向網(wǎng)格尺度。周向網(wǎng)格尺度應(yīng)用1.3 節(jié)的三維推廣即可。

3 算例驗(yàn)證

為驗(yàn)證第2 節(jié)所提出的網(wǎng)格劃分方法，本節(jié)應(yīng)用完全氣體模型，對(duì)高超聲速半圓柱、球頭等進(jìn)行數(shù)值模擬，對(duì)鈍錐及升力體等外形引用文獻(xiàn)給出的經(jīng)驗(yàn)網(wǎng)格劃分，對(duì)本文提出的網(wǎng)格劃分方法進(jìn)行驗(yàn)證。數(shù)值模擬均采用三維完全氣體NS 控制方程組；對(duì)流項(xiàng)離散采用效率較高且相對(duì)穩(wěn)定的van Leer 提出的迎風(fēng)數(shù)值格式，應(yīng)用MUSCL 方法使格式精度為3 階，并采用文獻(xiàn)［23］給出的熵修正方法減小該格式的數(shù)值黏性耗散，黏性項(xiàng)采用二階中心差分格式，時(shí)間導(dǎo)數(shù)項(xiàng)采用隱式推進(jìn)，以提高計(jì)算效率。

3.1 二維半圓柱繞流

3.1.1Ma∞=8 半圓柱繞流

首先選擇文獻(xiàn)［24］的圓柱繞流實(shí)驗(yàn)作為驗(yàn)證算例。圓柱直徑D=0.076 m，自由來(lái)流馬赫數(shù)Ma∞=8，靜壓P∞=855 Pa，靜溫T∞=125 K，壁面溫度Twall=294 K，基于圓柱直徑的來(lái)流雷諾數(shù)為Re=3.76×105，限于篇幅，圓柱繞流幾何模型可見(jiàn)文獻(xiàn)［7， 13］。根據(jù)第2 節(jié)給出的流場(chǎng)網(wǎng)格劃分方法，可得周向與法向第1 層網(wǎng)格尺度比為ΔxΔy～O(613)；再基于壁面參數(shù)的網(wǎng)格雷諾數(shù)Recell，w=15，代入自由來(lái)流參數(shù)，可給出Δnw=7.16×10-7m，由第2 節(jié)給出的方法可得周向網(wǎng)格尺度為Δx=4.38×10-4m，對(duì)應(yīng)的周向網(wǎng)格節(jié)點(diǎn)數(shù)約為270。

根據(jù)文獻(xiàn)已有的數(shù)值模擬對(duì)網(wǎng)格劃分方法的建議，可知用第2 節(jié)流場(chǎng)網(wǎng)格劃分方法，給出的半圓柱周向網(wǎng)格點(diǎn)數(shù)過(guò)多。究其原因，第2 節(jié)給出的方法，首先需要確定法向網(wǎng)格尺度，才能根據(jù)周向與法向網(wǎng)格尺度比，再確定周向網(wǎng)格尺度。若確定的法向網(wǎng)格尺度過(guò)于保守（網(wǎng)格尺度過(guò)小），必然導(dǎo)致給出的周向網(wǎng)格尺度也偏于保守，從而致使周向網(wǎng)格點(diǎn)數(shù)過(guò)多。

為驗(yàn)證上述的論述及分析，進(jìn)行以下數(shù)值模擬。第1 步，首先給定周向足夠多的網(wǎng)格點(diǎn)數(shù)，然后改變不同的法向網(wǎng)格尺度，應(yīng)用數(shù)值模擬與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，用數(shù)值方法找到此算例的最大法向網(wǎng)格尺度，確保此法向網(wǎng)格尺度的網(wǎng)格無(wú)關(guān)性及與實(shí)驗(yàn)結(jié)果的一致性。第2 步，再改變周向網(wǎng)格尺度（網(wǎng)格點(diǎn)數(shù)），驗(yàn)證本文提出的周向網(wǎng)格尺度與法向網(wǎng)格尺度的比是否適用。

圖1 為不同法向網(wǎng)格尺度下，半圓柱表面斯坦頓數(shù)St周向分布圖。由圖可見(jiàn)，在確保周向網(wǎng)格尺度足夠小時(shí)，當(dāng)法向網(wǎng)格尺度小于8×10-6m時(shí)，圓柱的斯坦頓數(shù)幾乎不隨法向網(wǎng)格變化而變化，根據(jù)文獻(xiàn)［1］的方法得到的Δnw=1.05×10-7m 過(guò)于保守，即應(yīng)用分子平均自由程確定法向網(wǎng)格尺度偏保守。而本文基于壁面參數(shù)的網(wǎng)格雷諾數(shù)Recell，w=15 給出的法向網(wǎng)格尺度雖然有所改進(jìn)，但還是過(guò)小。

圖1 Ma∞=8 時(shí)法向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響（半圓柱）Fig.1 Influence of normal grid scale over wall on Stanton number at Ma∞=8 （semi-cylindrical）

根據(jù)數(shù)值試驗(yàn)找到的法向網(wǎng)格尺度為Δnw=8×10-6m，然后再根據(jù)第2 節(jié)給出的周向與法向網(wǎng)格尺度比可得，Δx=4.9×10-3m，即半圓柱周向節(jié)點(diǎn)數(shù)為25。為驗(yàn)證本文給出的周向網(wǎng)格是否合理，分別選用壁面法向第1 層網(wǎng)格尺度為Δnw=8×10-6， 4×10-6m，通過(guò)改變周向節(jié)點(diǎn)數(shù)的斯坦頓數(shù)分布結(jié)果如圖2 所示。由圖可見(jiàn)，當(dāng)法向網(wǎng)格取值合理時(shí)，用本文的方法給出的周向網(wǎng)格的取值亦合適。并且從圖2（b）可見(jiàn)，當(dāng)法向網(wǎng)格取值過(guò)于保守時(shí)，與實(shí)驗(yàn)結(jié)果的對(duì)比可知，周向劃分25 個(gè)網(wǎng)格點(diǎn)，對(duì)計(jì)算結(jié)果幾乎無(wú)影響。但根據(jù)法向網(wǎng)格與周向網(wǎng)格的比可知，此時(shí)要求周向50 個(gè)網(wǎng)格點(diǎn)。因此，若確定的法向網(wǎng)格尺度過(guò)于保守，必然導(dǎo)致給出的周向網(wǎng)格尺度也偏于保守，從而致使周向網(wǎng)格點(diǎn)數(shù)過(guò)多。

圖2 Ma∞=8 時(shí)周向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.2 Influence of circumferential grid scale on Stanton number at Ma∞=8

3.1.2Ma∞=16.34 半圓柱繞流

為進(jìn)一步驗(yàn)證第2 節(jié)的網(wǎng)格劃分方法，以及3.1.1 節(jié)的初步結(jié)論。本節(jié)選擇文獻(xiàn)［24］給出的實(shí)驗(yàn)結(jié)果進(jìn)行驗(yàn)證。圓柱直徑D=0.076 m，自由來(lái)流的馬赫數(shù)、靜溫、靜壓分別為Ma∞=16.34、T∞=52 K、P∞=82.95 Pa；壁 面 溫 度Twall=294.4 K，基于圓柱直徑的來(lái)流雷諾數(shù)為Re=2.97×105。 根 據(jù) 第2 節(jié) 確 定 的Δnw=1.77×10-6m、周向網(wǎng)格尺度Δx=9.61×10-4m，對(duì)應(yīng)的周向節(jié)點(diǎn)數(shù)約為125。與3.1.1 節(jié)類(lèi)似，當(dāng)給出的法向網(wǎng)格尺度過(guò)小時(shí)，周向網(wǎng)格尺度亦過(guò)小。

基于3.1.1 節(jié)給出的方法，應(yīng)用數(shù)值模擬方法，根據(jù)圖3 可知數(shù)值計(jì)算確定合理的法向網(wǎng)格尺度Δnw=2.0×10-5m，可得周向網(wǎng)格尺度上限為Δx=1×10-2m，對(duì)應(yīng)的周向網(wǎng)格節(jié)點(diǎn)數(shù)為11。由圖4 亦可得，當(dāng)法向網(wǎng)格取值合理時(shí)，用本文的方法給出的流向網(wǎng)格的取值亦合適。

圖3 Ma∞=16.34 時(shí)法向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.3 Influence of normal grid scale over wall on Stanton number at Ma∞=16.34

圖4 Ma∞=16.34 時(shí)周向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.4 Influence of circumferential grid scale on Stanton number at Ma∞=16.34

3.2 三維球頭繞流

3.2.1Ma∞=8 球頭繞流

為驗(yàn)證第2 節(jié)的網(wǎng)格劃分方法以及3.1.1 節(jié)的初步結(jié)論，本節(jié)選用文獻(xiàn)［17］給出的三維球頭算例進(jìn)行數(shù)值驗(yàn)證。球頭直徑D=15 mm，自由來(lái)流參數(shù)為Ma∞=8，T0=892 K，P0=7.8 MPa，Twall=294 K，基于球頭直徑的特征雷諾數(shù)為Re=1.7×105。根據(jù)第2 節(jié)確定的Δnw=7.5×10-7m，對(duì) 應(yīng) 的 切 向 網(wǎng) 格 尺 度 為hτ=3.09×10-4m，對(duì)應(yīng)的節(jié)點(diǎn)數(shù)約為38。與3.1.1 節(jié)類(lèi)似，即當(dāng)給出的法向網(wǎng)格尺度過(guò)小時(shí)，周向網(wǎng)格尺度亦過(guò)小。

采用與3.1 節(jié)相同方法，應(yīng)用數(shù)值模擬方法確定法向網(wǎng)格尺度上限。根據(jù)圖5 可知，數(shù)值計(jì)算確定合理的法向網(wǎng)格尺度Δnw=4.3×10-6m，對(duì)應(yīng)的切向網(wǎng)格節(jié)點(diǎn)數(shù)為9。由圖6 可得，當(dāng)法向網(wǎng)格尺度取值合理時(shí)，用本文的方法給出的切向網(wǎng)格尺度的取值亦合適。

圖5 Ma∞=8 時(shí)法向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響（球頭）Fig.5 Influence of normal grid scale over wall on Stanton number at Ma∞=8 （sphere）

圖6 Ma∞=8 時(shí)切向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.6 Influence of circumferential grid scale on Stanton number at Ma∞=8

圖7 給出了球頭表面斯坦頓數(shù)分布云圖，由圖可見(jiàn)，數(shù)值方法較好地模擬了球頭表面的斯坦頓數(shù)分布。

圖7 Ma∞=8 球頭壁面斯坦頓數(shù)分布Fig.7 Stanton number distribution over sphere at Ma∞=8

3.2.2Ma∞=13 球頭繞流

為驗(yàn)證第2 節(jié)的網(wǎng)格劃分方法在來(lái)流馬赫數(shù)較高的球頭算例適用性，本節(jié)選用文獻(xiàn)［25］中歐洲風(fēng)洞F4 試驗(yàn)的標(biāo)準(zhǔn)球頭進(jìn)行了數(shù)值計(jì)算，其直徑D=0.35 m，試驗(yàn)來(lái)流參數(shù)為U∞=4 320 m/s，T∞=265 K，P∞=52.81 Pa，Twall=294.4 K，可得基于球頭直徑的特征雷諾數(shù)為Re=6.26×104。由第2 節(jié)的方法則法向網(wǎng)格尺度為Δnw=4.23×10-6m，切向網(wǎng)格與法向網(wǎng)格尺度比約為250，可得切向網(wǎng)格尺度為hτ=1.06×10-3m。

由前述可知，該尺度偏于保守，采用同前兩節(jié)的方法，采用數(shù)值模擬方法確定網(wǎng)格尺度上限。由圖8 可知，數(shù)值計(jì)算確定的Δnw=3.9×10-5m，可得切網(wǎng)格尺度hτ=9.75×10-3m，所對(duì)應(yīng)的切向網(wǎng)格節(jié)點(diǎn)數(shù)為28。由圖9 可得，當(dāng)法向網(wǎng)格尺度取值合理時(shí)，用本文的方法給出的切向網(wǎng)格尺度的取值亦合適。

圖8 Ma∞=13 時(shí)法向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.8 Influence of normal grid scale over wall on Stanton number at Ma∞=13

圖9 Ma∞=13 時(shí)切向網(wǎng)格尺度對(duì)斯坦頓數(shù)的影響Fig.9 Influence of circumferential grid scale on Stanton number at Ma∞=13

圖10 給出了球頭表面斯坦頓數(shù)分布云圖。由圖可見(jiàn)，數(shù)值方法較好地模擬了球頭表面的斯坦頓數(shù)分布。

圖10 Ma∞=13 球頭壁面斯坦頓數(shù)分布Fig.10 Stanton number distribution over sphere at Ma∞=13

3.3 文獻(xiàn)算例驗(yàn)證

為進(jìn)一步驗(yàn)證高超聲速流動(dòng)中誤差匹配原則的適用性，本節(jié)引用已發(fā)表文獻(xiàn)中相關(guān)數(shù)值計(jì)算算例，驗(yàn)證復(fù)雜外形高超聲速數(shù)值計(jì)算中誤差匹配原則的適用性。

3.3.1Ma∞=10 鈍錐繞流

本節(jié)通過(guò)文獻(xiàn)［18］中鈍頭圓錐外形驗(yàn)證第2 節(jié)中給出的網(wǎng)格劃分方法的適用性，其流場(chǎng)的數(shù)值模擬通過(guò)隱式TVD 格式和LU 分解算法進(jìn)行求解三維N-S 方程，鈍錐由直徑D=0.055 88 m 的球頭和半錐角為15°的錐體組成，總長(zhǎng)為0.447 04 m，幾何外形可見(jiàn)文獻(xiàn)［18］。其來(lái)流參數(shù)為Ma∞=10.6，T∞=47.3 K，P∞=109.9 Pa，Twall=294.4 K，攻角為0°。可得基于球頭直徑的特征雷諾數(shù)為Re=2.2×105，由第2 節(jié)的方法，可得法向網(wǎng)格尺度為3.18×10-6m，切向網(wǎng)格與法向網(wǎng)格尺度比約為460，切向網(wǎng)格尺度為hτ=1.46×10-3m。文獻(xiàn)［18］給出的經(jīng)網(wǎng)格無(wú)關(guān)性驗(yàn)證的壁面法向第1層網(wǎng)格尺度為Δnw=4.2×10-6m，切向網(wǎng)格尺度hτ＜1.397×10-3m，與應(yīng)用第2 節(jié)方法給出的切向網(wǎng)格尺度與法向網(wǎng)格尺度比相近。

3.3.2Ma∞=15 升力體繞流

本節(jié)通過(guò)文獻(xiàn)［26］中升力體外形驗(yàn)證第2 節(jié)給出的網(wǎng)格劃分方法的適用性，采用有限體積法對(duì)控制方程進(jìn)行離散，數(shù)值格式為van Leer，應(yīng)用3 階精度的MUSCL 方法使數(shù)值精度為3 階，時(shí)間推進(jìn)采用隱式LU-SGS 方法。其升力體頭部直徑D=140 mm，模型總長(zhǎng)為3.3 m，幾何外形見(jiàn)文獻(xiàn)［26］。來(lái)流參數(shù)為Ma∞=15，T∞=250.35 K，P∞=287.4 Pa，Twall=300 K；基于升力體頭部直徑的自由來(lái)流Re=1.68×105。由第2節(jié)的方法，可得法向網(wǎng)格尺度為6.2×10-7m，切向網(wǎng)格與法向網(wǎng)格尺度比約為410、切向網(wǎng)格尺度為hτ=2.54×10-4m。文獻(xiàn)［26］給出的經(jīng)網(wǎng)格無(wú)關(guān)性驗(yàn)證的壁面法向第1 層網(wǎng)格尺度為Δnw=8.4×10-6m，切向網(wǎng)格尺度hτ＜4.1×10-3，與應(yīng)用第2節(jié)方法給出的切向網(wǎng)格尺度與法向網(wǎng)格尺度比一致。

4 討論及分析

4.1 表面熱流收斂性

高超聲速數(shù)值模擬中，表面熱流密度的收斂性較慢。本文采用的表面熱流收斂要求是，首先計(jì)算至殘差不再下降，并穩(wěn)定在一定的數(shù)值范圍后，再取2 個(gè)相鄰5 000 步的熱流數(shù)值不再變化作為收斂判斷準(zhǔn)則。3.1 節(jié)中Ma∞=8，16.34 圓柱繞流與3.2 節(jié)中Ma∞=8，13 球頭繞流算例殘差歷程如圖11 所示。

圖11 數(shù)值計(jì)算殘差歷程Fig.11 Residual history for different numerical simulations

以上方法確保了本文高超聲速熱流計(jì)算的收斂性。

4.2 數(shù)值格式對(duì)表面熱流的影響

表1 給出了本文及文獻(xiàn)［17-18， 26］計(jì)算鈍頭體表面熱流的網(wǎng)格劃分尺度。由表1 可見(jiàn)，切向網(wǎng)格尺度對(duì)高超聲速鈍頭體表面熱流的影響不如法向網(wǎng)格影響顯著，并且只要鈍頭體壁面法向網(wǎng)格尺度合適，根據(jù)第2 節(jié)給出的切向網(wǎng)格劃分尺度亦合理。

表1 不同算例網(wǎng)格劃分尺度匯總Table 1 Collections of grid scales for different numerical cases

表1 中文獻(xiàn)［17］及文獻(xiàn)［26］雖然采用了改進(jìn)了熵修正函數(shù)的二階TVD 格式，但是數(shù)值模擬給出的壁面網(wǎng)格雷諾數(shù)相比其他算例還是過(guò)小。而文獻(xiàn)［6］由于采用未改進(jìn)二階TVD 格式，為了能準(zhǔn)確計(jì)算表面熱流，所要求的法向網(wǎng)格尺度為壁面分子自由程的大小［1］（此時(shí)基于壁面參數(shù)的網(wǎng)格雷諾數(shù)為1.5）。其原因在于，上述2 個(gè)文獻(xiàn)采用的二階精度TVD 格式在熱流計(jì)算上由于格式的高耗散型而存在困難［27］，并且文獻(xiàn)［28］發(fā)現(xiàn)即使物面網(wǎng)格很小，原TVD 格式計(jì)算駐點(diǎn)熱流仍較實(shí)驗(yàn)和文獻(xiàn)給出的結(jié)果偏低。而文獻(xiàn)［17， 26］由于采用了焓守恒修正的van Leer 格式數(shù)值離散，并用MUSCL 插值到3 階精度，從而有效抑制了數(shù)值格式的黏性效應(yīng)，所以對(duì)壁面網(wǎng)格雷諾數(shù)放寬了要求。文獻(xiàn)［11］研究結(jié)果表明，當(dāng)采用4 種差分格式，并用MUSCL 插值到二階精度時(shí)，Roe 格式的壁面雷諾數(shù)可放寬到約50（此處根據(jù)該文獻(xiàn)基于來(lái)流參數(shù)的網(wǎng)格雷諾數(shù)，換算為基于壁面參數(shù)的網(wǎng)格雷諾數(shù)），而SLAM格式由于有效抑制了數(shù)值格式的黏性，壁面雷諾數(shù)甚至可放寬到約90。文獻(xiàn)［10］的研究結(jié)果則表明，不同迎風(fēng)型數(shù)值格式對(duì)網(wǎng)格雷諾數(shù)的要求不同，應(yīng)選擇數(shù)值黏性小，則可放寬對(duì)網(wǎng)格雷諾數(shù)的要求。

4.3 壁溫及駐點(diǎn)熱流密度的影響

由于高超聲速飛行器表面熱流影響參量較多，以下研究圓柱及球頭不同表面溫度，以進(jìn)一步驗(yàn)證本文結(jié)論的適應(yīng)性。圖12 給出了不同外形、不同壁面溫度時(shí)，法向網(wǎng)格雷諾數(shù)對(duì)斯坦頓數(shù)的影響。由圖可見(jiàn)，壁面溫度不同時(shí)，網(wǎng)格雷諾數(shù)上限范圍幾乎不變，即根據(jù)本文提出的方法確定的網(wǎng)格尺度對(duì)壁面溫度不敏感。

圖12 不同壁溫下法向網(wǎng)格尺度對(duì)表面斯坦頓數(shù)的影響Fig.12 Influence of normal grid scale on Stanton number under different skin temperature conditions

表2 給出了本文及參考文獻(xiàn)［18， 26］中算例流場(chǎng)參數(shù)及駐點(diǎn)熱流密度等數(shù)值，其中Tt為總溫。由表2 可見(jiàn)，本文得出的結(jié)論對(duì)不同來(lái)流總溫與表面溫度比TtTw、變化范圍較大（約50 倍）的駐點(diǎn)熱流密度qs均具有一定的適應(yīng)性。

表2 不同算例流場(chǎng)參數(shù)匯總Table 2 Collections of flow field parameters for different numerical cases

綜上，對(duì)采用數(shù)值模擬進(jìn)行高超聲速飛行器表面熱流預(yù)估時(shí)，若所選用的數(shù)值格式未進(jìn)行數(shù)值黏性影響研究，那么在數(shù)值模擬的具體實(shí)施中需進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。對(duì)于數(shù)值黏性小的格式，可取的壁面網(wǎng)格雷諾數(shù)分別為160、80、40 進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證；若采用數(shù)值黏性較大的格式，則需要采用40、20、10 的壁面網(wǎng)格雷諾數(shù)進(jìn)行驗(yàn)證。而流向網(wǎng)格的劃分，則可采用3 套網(wǎng)格中間的壁面網(wǎng)格雷諾數(shù)給定的壁面法向第1 層網(wǎng)格尺度，應(yīng)用流向和法向網(wǎng)格尺度比的方法，確定流向網(wǎng)格尺度，并采用等間距的流向網(wǎng)格劃分。若所選擇的數(shù)值格式黏性較小，壁面法向網(wǎng)格雷諾數(shù)可增大到40 以上。

5 結(jié) 論

應(yīng)用理論分析及數(shù)值模擬，提出了一種基于壁面網(wǎng)格雷諾數(shù)及基于鈍頭體特征長(zhǎng)度的來(lái)流雷諾數(shù)的網(wǎng)格劃分方法，并給出了壁面網(wǎng)格雷諾數(shù)的取值范圍。應(yīng)用所提出的網(wǎng)格劃分方法對(duì)不同來(lái)流條件的高超聲速半圓柱及球頭表面熱流進(jìn)行了數(shù)值模擬，結(jié)論如下：

1）法向方向網(wǎng)格尺度對(duì)高超聲速鈍頭體表面熱流的影響大于切向網(wǎng)格尺度的影響。

2）基于能量方程的修正方程各誤差項(xiàng)匹配分析表明，為使各誤差項(xiàng)為同一量級(jí)從而減小數(shù)值誤差，流向與法向第1 層網(wǎng)格尺度比應(yīng)為基于鈍體特征長(zhǎng)度的來(lái)流雷諾數(shù)的平方根。

3）基于壁面網(wǎng)格雷諾數(shù)確定的法向網(wǎng)格尺度的取值受數(shù)值格式影響，數(shù)值黏性小的格式可取的壁面網(wǎng)格雷諾數(shù)大。對(duì)于工程上采用的能有效抑制數(shù)值黏性的格式，壁面網(wǎng)格雷諾數(shù)可供參考的取值范圍為20～80。

應(yīng)該注意，本文給出的方法受所選擇的數(shù)值格式的耗散性影響較大，在應(yīng)用本文給出的網(wǎng)格劃分方法前，應(yīng)選擇經(jīng)評(píng)估的、數(shù)值耗散性小的數(shù)值格式。對(duì)數(shù)值黏性大的格式則不能進(jìn)一步放寬網(wǎng)格尺度的要求。

另外，雖然本文給出了一種計(jì)算高超聲速鈍體表面熱流的網(wǎng)格劃分方法，但高超聲速真實(shí)飛行器表面熱流的準(zhǔn)確計(jì)算，受數(shù)值方法、網(wǎng)格質(zhì)量、邊界條件等多種因素的影響，因此真實(shí)高超聲速飛行器表面熱流的精確預(yù)估還需要深入研究，以滿(mǎn)足工程應(yīng)用的需求。