高馬赫數(shù)層流摩阻數(shù)值計(jì)算精度

范月華，段毅，周乃楨，楊攀

中國運(yùn)載火箭技術(shù)研究院 空間物理重點(diǎn)實(shí)驗(yàn)室，北京 100076

飛行器在高空高速飛行的顯著特征是高馬赫數(shù)和低雷諾數(shù)，此時(shí)黏性與黏性干擾效應(yīng)十分顯著，極端工況下摩擦阻力在總阻中的占比甚至超過90%[1]，對(duì)飛行器升阻特性等關(guān)鍵氣動(dòng)性能影響顯著。因此，準(zhǔn)確預(yù)估飛行器的摩擦阻力具有重要意義。

目前，地面試驗(yàn)還無法覆蓋實(shí)際飛行工況[2-4]，而飛行試驗(yàn)中還缺乏可行的科學(xué)測量手段獲得有效的摩阻數(shù)據(jù)[5]，因此計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)技術(shù)是高空高速狀態(tài)下摩阻預(yù)示最為有效的手段。近幾十年來，隨著CFD技術(shù)的不斷突破和發(fā)展[6]，航空航天氣動(dòng)設(shè)計(jì)的基本問題越來越多地依賴這一技術(shù)開展研究[7]，但與飛行器阻力特性密切相關(guān)的摩擦阻力的精確預(yù)示一直是CFD的難題和熱點(diǎn)之一。飛行試驗(yàn)結(jié)果分析發(fā)現(xiàn)，通過數(shù)值計(jì)算確定的軸向力系數(shù)與試驗(yàn)辨識(shí)結(jié)果存在較大差異，主要是由摩阻預(yù)示的不準(zhǔn)確造成的?；诂F(xiàn)有的認(rèn)識(shí)，即使考慮高溫氣體效應(yīng)和稀薄氣體效應(yīng)的影響，也無法充分解釋上述摩阻預(yù)示的差異。結(jié)合地面風(fēng)洞試驗(yàn)的結(jié)果，該差異主要源于數(shù)值計(jì)算的摩擦阻力結(jié)果偏高。因此，需緊密結(jié)合高空高速飛行特點(diǎn)，系統(tǒng)研究數(shù)值方法對(duì)摩阻計(jì)算精度的影響，為黏性效應(yīng)的精細(xì)化模擬提供支撐。

在前期工作中，周丹等[8]研究了網(wǎng)格分布、通量格式、限制器對(duì)低速平板層流摩阻計(jì)算精度的影響，鄭世超[9]研究了插值格式精度、通量格式和網(wǎng)格雷諾數(shù)對(duì)于高超平板層流摩阻計(jì)算的影響，張培紅等[10]研究了提高混合網(wǎng)格摩阻預(yù)測精度的熵修正方法。但上述研究均以平板或者低速流動(dòng)為研究對(duì)象，對(duì)于實(shí)際高速飛行器外形摩阻計(jì)算精度的相關(guān)研究仍較少，無法為高馬赫數(shù)流動(dòng)中黏性效應(yīng)的精細(xì)化模擬提供支撐。

本文以具有高速飛行器典型部件特征的球錐、大后掠角三角翼為對(duì)象，結(jié)合風(fēng)洞試驗(yàn)?zāi)ψ铚y量結(jié)果，研究了數(shù)值計(jì)算中影響摩阻計(jì)算的數(shù)值耗散及邊界條件等重要因素，并基于分析結(jié)果提出了對(duì)CFD技術(shù)的發(fā)展需求。第1節(jié)首先介紹表面摩阻計(jì)算與試驗(yàn)測量結(jié)果存在差異的現(xiàn)象，第2節(jié)分析數(shù)值計(jì)算中影響表面摩阻的幾個(gè)重要因素，第3節(jié)基于分析結(jié)果提出CFD技術(shù)的發(fā)展需求，第4節(jié)給出了研究結(jié)論。

1 摩阻計(jì)算與試驗(yàn)測量間的差異現(xiàn)象

1.1 圓錐表面摩阻計(jì)算與試驗(yàn)對(duì)比

Meritt等[11]基于AEDC(Arnold Engineering Development Complex)9號(hào)風(fēng)洞開展了圓錐模型的表面摩阻測量試驗(yàn)，試驗(yàn)狀態(tài)包括層流、轉(zhuǎn)捩、湍流3種不同流態(tài)。模型長度為1 551.5 mm(從理論頂點(diǎn)算起)，半錐角為7°，底端面直徑為381 mm，端頭半徑為0.15 mm；摩阻測量傳感器安裝位置距離理論頂點(diǎn)1 350 mm，位于迎風(fēng)面中心子午線一側(cè)，與其夾角為10°。模型及傳感器的位置如圖1所示，攻角(Angle of Attack，AOA)在xy平面內(nèi)。

圖1 Meritt試驗(yàn)?zāi)Ｐ图皞鞲衅魑恢檬疽鈭D[11]Fig.1 Sketch of cone model and sensor location of Meritt[11]

分別使用國家數(shù)值風(fēng)洞(NNW)工程研發(fā)的計(jì)算流體力學(xué)軟件Flowstar[12](V1.0版)、商業(yè)軟件CFD++和自研程序HyperCFD，針對(duì)其中的層流狀態(tài)試驗(yàn)條件開展數(shù)值模擬，來流參數(shù)如表1所示。其中，F(xiàn)lowstar基于非結(jié)構(gòu)網(wǎng)格，使用T-Rex技術(shù)保證邊界層網(wǎng)格密度，邊界層第1層網(wǎng)格高度為0.01 mm，無黏通量空間離散使用二階HLLE++格式；CFD++和HyperCFD基于結(jié)構(gòu)網(wǎng)格，計(jì)算網(wǎng)格邊界層第1層高度為0.01 mm，無黏通量空間離散CFD++使用HLLC格式[13]，HyperCFD使用Roe格式[14]。壁面溫度均設(shè)為300 K。

表1 圓錐模型層流狀態(tài)風(fēng)洞試驗(yàn)參數(shù)

層流狀態(tài)下測點(diǎn)處摩擦阻力系數(shù)數(shù)值計(jì)算結(jié)果與試驗(yàn)的對(duì)比如圖2所示，可以看到，F(xiàn)lowstar的模擬精度與其他CFD軟件接近，攻角α=10°時(shí)Flowstar相比于CFD++和HyperCFD分別偏大約4%和7%，不同軟件的計(jì)算結(jié)果與試驗(yàn)測量結(jié)果偏差較大，相比試驗(yàn)偏大24%～30%。

圖2 圓錐模型表面摩阻系數(shù)數(shù)值計(jì)算與風(fēng)洞試驗(yàn)[11]差異Fig.2 Difference of skin friction coefficients of cone model between numerical calculation and wind tunnel test[11]

1.2 三角翼表面摩阻計(jì)算與試驗(yàn)對(duì)比

試驗(yàn)用三角翼標(biāo)準(zhǔn)模型全長363.34 mm，前緣后掠角為76.38°，分別采用摩阻天平和液晶涂層測量技術(shù)測量單點(diǎn)的摩擦應(yīng)力和物面的摩擦應(yīng)力分布，風(fēng)洞來流參數(shù)如表2所示。摩阻天平測點(diǎn)位置在下表面中心線上，距頭部頂點(diǎn)140 mm。

表2 三角翼模型風(fēng)洞試驗(yàn)參數(shù)Table 2 With tunnel test parameters of delta-wing model

分別使用Flowstar、CFD++和HyperCFD對(duì)上述試驗(yàn)工況進(jìn)行數(shù)值模擬。其中，F(xiàn)lowstar使用基于T-Rex技術(shù)的非結(jié)構(gòu)網(wǎng)格，邊界層第1層網(wǎng)格高度為0.01 mm；CFD++和HyperCFD同時(shí)采用2套結(jié)構(gòu)化網(wǎng)格進(jìn)行計(jì)算以評(píng)估網(wǎng)格密度對(duì)摩阻計(jì)算精度的影響：粗網(wǎng)格流向、法向和周向網(wǎng)格點(diǎn)數(shù)分別為147、81、327，邊界層第1層高度為0.01 mm；密網(wǎng)格流向、法向和周向網(wǎng)格點(diǎn)數(shù)分別為301、301、327，邊界層第1層高度為0.005 mm。壁面溫度均設(shè)為300 K。

來流馬赫數(shù)8、0°攻角下表面中心線摩阻系數(shù)層流計(jì)算結(jié)果與試驗(yàn)對(duì)比如圖3所示。Flowstar獲得的層流狀態(tài)摩阻系數(shù)結(jié)果介于CFD++和HyperCFD之間，流向位置140 mm處Flowstar計(jì)算結(jié)果與HyperCFD結(jié)果相比約偏大15%；不同軟件計(jì)算的摩阻系數(shù)都明顯高于試驗(yàn)結(jié)果，流向位置140 mm處約為試驗(yàn)值的1.5～2.0倍。

圖3 Ma∞=8、0°攻角下表面中心線摩阻系數(shù)層流計(jì)算結(jié)果與試驗(yàn)對(duì)比Fig.3 Comparison of friction coefficients along center windward ray between laminar calculation and wind tunnel test under conditions of Ma∞=8 and α=0°

通過以上2個(gè)計(jì)算與試驗(yàn)結(jié)果的對(duì)比，可以發(fā)現(xiàn)，即便不考慮流動(dòng)轉(zhuǎn)捩的影響，高馬赫數(shù)層流的摩擦阻力計(jì)算仍然整體上高于試驗(yàn)測量結(jié)果。

2 摩阻數(shù)值計(jì)算精度的影響因素

計(jì)算網(wǎng)格是影響CFD計(jì)算結(jié)果的一個(gè)重要因素，在計(jì)算網(wǎng)格相同的情況下，不同的計(jì)算方法與設(shè)置也會(huì)對(duì)摩阻產(chǎn)生非常明顯的影響。網(wǎng)格的影響不在本文的討論范疇，本節(jié)僅關(guān)注CFD軟件本身的計(jì)算方法與設(shè)置對(duì)摩阻計(jì)算的影響。如無特殊說明，本節(jié)計(jì)算均基于1.1節(jié)的圓錐模型開展，來流條件保持不變。

2.1 數(shù)值耗散的影響

使用自研的HyperCFD結(jié)構(gòu)網(wǎng)格求解器對(duì)表1所列的層流風(fēng)洞試驗(yàn)狀態(tài)進(jìn)行數(shù)值模擬，對(duì)比不同空間格式和熵修正方法對(duì)摩阻計(jì)算的影響，壁面條件為等溫壁，壁溫為300 K。

2.1.1 無黏通量格式

對(duì)比Roe、AUSM+[15-16]、Van Leer[17]3種不同空間離散格式對(duì)摩阻計(jì)算的影響，其中Roe和AUSM+格式使用MUSCL(Monotone Upstream-centered Schemes for Conservation Laws)插值得到空間二階精度[18]，同時(shí)使用minmod限制器防止激波振蕩。此外，Roe格式需要采用熵修正來抑制非物理解的產(chǎn)生并保持計(jì)算穩(wěn)定，常用的熵修正形式如下[19]:

(1)

式中：λ為Navier-Stokes方程無黏通量Jacobi矩陣特征值；δ為熵修正閾值，當(dāng)特征值的絕對(duì)值低于該值時(shí)進(jìn)行熵修正。δ的值可以直接取常數(shù)，也可以表示為特征值的函數(shù)，根據(jù)δ的不同形成了各類熵修正公式[20]。

Van Leer格式屬于矢通量分裂類格式，本身的數(shù)值耗散較大，黏性模擬精度較低，從而導(dǎo)致摩阻計(jì)算偏大，不適合進(jìn)行邊界層等強(qiáng)剪切流的模擬。Roe格式屬于通量差分裂類格式，可以得到線化Riemann問題的精確解，對(duì)線性波的分辨率較高，因而對(duì)黏性剪切層的模擬也相對(duì)精確。圖4 給出了圓錐模型不同空間格式測點(diǎn)摩阻系數(shù)的計(jì)算對(duì)比，其中Grid 1對(duì)應(yīng)的網(wǎng)格邊界層第1層高度為0.01 mm，Grid 2對(duì)應(yīng)的網(wǎng)格邊界層第1層高度為0.002 mm，Roe格式熵修正中δ的計(jì)算使用Müller提出的方法[21], eps0.1表示熵修正比例系數(shù)為0.1?？梢钥吹剑赗oe和AUSM+格式的2套網(wǎng)格獲得的摩阻系數(shù)基本一致，最大相差2%；Van Leer格式對(duì)空間網(wǎng)格的敏感度更高，20°攻角時(shí)2套網(wǎng)格的摩阻系數(shù)相差約8%。不同空間格式間的摩阻計(jì)算結(jié)果差異明顯，而且隨著攻角增加差異趨于明顯；Roe格式計(jì)算的摩阻最小，AUSM+格式比Roe格式略微偏大2%～4%，Van Leer格式的摩阻最大，相比于Roe格式偏大20%～30%。結(jié)合相關(guān)理論及上述計(jì)算結(jié)果，可知黏性分辨率越高、數(shù)值耗散越小，黏性摩擦阻力的計(jì)算結(jié)果越精確。

圖4 不同空間格式測點(diǎn)摩阻系數(shù)對(duì)比Fig.4 Comparison of measurement location friction coefficients among different space schemes

2.1.2 熵修正方法

基于二階Roe格式，對(duì)比Müller[21]和Harten-Yee[22]2種熵修正方法以及不同熵修正比例系數(shù)對(duì)摩阻計(jì)算的影響。2種熵修正方法均基于式(1)，但是δ的取值不同，Harten-Yee熵修正方法δ的表達(dá)式為

(2)

式中：δ*為熵修正比例系數(shù)；V為速度矢量；a為當(dāng)?shù)芈曀?ξ、η、ζ為計(jì)算坐標(biāo)系的3個(gè)方向。Müller熵修正方法δ的表達(dá)式為(以ξ方向?yàn)槔?

(3)

圖5給出了圓錐模型不同熵修正方法和熵修正比例系數(shù)測點(diǎn)摩阻系數(shù)的計(jì)算對(duì)比?？梢钥吹剑?種熵修正方法的測點(diǎn)摩阻系數(shù)計(jì)算結(jié)果相差4%～8%，Müller型熵修正對(duì)黏性邊界層的模擬精度更高。相比于Harten-Yee熵修正，Müller熵修正方法考慮了速度和網(wǎng)格的各向異性，對(duì)網(wǎng)格長細(xì)比較大的邊界層區(qū)域具有更高的分辨率。此外，2種熵修正方法均表現(xiàn)為熵修正比例系數(shù)越大，Roe格式的數(shù)值耗散越大，計(jì)算的摩阻系數(shù)也相應(yīng)越大。

從圖5中可以得到的結(jié)論是，表面摩阻計(jì)算結(jié)果與格式的數(shù)值耗散特性密切相關(guān)，對(duì)于Roe格式而言，熵修正比例系數(shù)越小，計(jì)算摩阻系數(shù)越小?；诖?，本文嘗試在邊界層內(nèi)近壁區(qū)關(guān)閉熵修正，圖6給出了這種處理方式下10°攻角時(shí)圓錐模型迎風(fēng)面子午線(迎風(fēng)面中心線)的摩阻系數(shù)，并與常規(guī)熵修正(Harten-Yee eps0.1)的摩阻結(jié)果進(jìn)行了對(duì)比。

圖5 不同熵修正方法測點(diǎn)摩阻系數(shù)對(duì)比Fig.5 Comparison of measurement location friction coefficients among different entropy fix methods

圖6 近壁區(qū)無熵修正時(shí)10°攻角迎風(fēng)面子午線摩阻系數(shù)Fig.6 Friction coefficient along center windward ray while closing entropy fix near wall under conditions of α=10°

可以看到，邊界層內(nèi)近壁區(qū)關(guān)閉熵修正得到的摩阻系數(shù)在120 mm后小于常規(guī)熵修正，在1 350 mm處偏小7%。然而由于處理方式比較簡單，導(dǎo)致圓錐頭部區(qū)域(120 mm之前)的摩阻系數(shù)有一定程度的增加。

2.2 壁面溫度邊界條件的影響

2.2.1 壁面溫度

從低速到低超聲速流動(dòng)，來流總溫較低，飛行器表面溫度變化不明顯，壁溫對(duì)摩阻計(jì)算的影響問題不突出；而高馬赫數(shù)飛行時(shí)，氣動(dòng)加熱現(xiàn)象顯著，計(jì)算時(shí)壁面溫度的選取對(duì)于摩阻計(jì)算會(huì)產(chǎn)生較大影響。

第一，根據(jù)皖河流域山區(qū)環(huán)境的特點(diǎn)，構(gòu)建低耗、優(yōu)質(zhì)、高產(chǎn)、高效的農(nóng)田生態(tài)系統(tǒng)。主要內(nèi)容有：山區(qū)土地整治和土壤改良，保持耕地的綠色覆蓋，建設(shè)生態(tài)水系和現(xiàn)代灌溉系統(tǒng)，山坡耕地保護(hù)或退耕等。

對(duì)表1所列的層流風(fēng)洞試驗(yàn)狀態(tài)進(jìn)行數(shù)值模擬，對(duì)比不同壁面溫度條件對(duì)表面摩阻的影響，無黏通量空間離散使用二階Roe格式，黏性通量使用二階中心差分格式。圖7給出了不同壁溫條件下，圓錐模型在0°攻角(圖7(a))和20°攻角(圖7(b)) 迎風(fēng)面子午線摩阻系數(shù)的對(duì)比情況，其中Tw為壁面溫度。

圖7 不同壁溫下圓錐外形迎風(fēng)面子午線摩阻系數(shù)對(duì)比Fig.7 Comparison of friction coefficients along center windward ray of cone model at different wall temperatures

0°攻角時(shí)，邊界層流動(dòng)相對(duì)簡單，壁面摩阻對(duì)溫度邊界條件不敏感，流向位置60 mm之前，絕熱壁的當(dāng)?shù)啬ψ韪?；隨著邊界層發(fā)展，流向位置200 mm之后，絕熱壁的當(dāng)?shù)啬ψ璧陀诘葴乇凇?/p>

20°攻角時(shí)，壁面摩阻隨壁面溫度的變化相對(duì)復(fù)雜。流向位置200～900 mm，壁溫250 K(低壁溫)對(duì)應(yīng)的迎風(fēng)面子午線的壁面摩阻小于壁溫300 K的，二者在460 mm處最多相差13%；900 mm 之后，壁溫250 K對(duì)應(yīng)的迎風(fēng)面子午線的壁面摩阻要更大，二者在1 220 mm處最多相差9%。這說明，當(dāng)圓錐有攻角時(shí)，出現(xiàn)繞圓錐的橫向流動(dòng)，邊界層內(nèi)流動(dòng)更加復(fù)雜，壁溫對(duì)流場的影響較大；此時(shí)，不僅壁面附近網(wǎng)格上的耗散特性影響表面摩阻，整個(gè)邊界層流動(dòng)的模擬精度都會(huì)對(duì)摩阻產(chǎn)生較大影響。高馬赫數(shù)流動(dòng)時(shí)，不同物面位置處，壁溫對(duì)表面摩阻的影響也不相同，這種現(xiàn)象對(duì)于復(fù)雜的飛行器外形更加突出，因此在進(jìn)行數(shù)值模擬時(shí)真實(shí)壁面溫度的選取對(duì)于計(jì)算結(jié)果的精度至關(guān)重要。

圖8給出了三角翼模型在Ma∞=8來流中進(jìn)行摩阻測量時(shí)紅外測得的表面溫度場，局部區(qū)域(尤其是頭部和側(cè)緣)氣動(dòng)加熱明顯?？紤]到壁面溫度設(shè)置對(duì)摩阻計(jì)算結(jié)果的影響，如果CFD軟件能夠?qū)崿F(xiàn)基于邊界層當(dāng)?shù)亓鲌鲎兞康谋诿鏈囟茸赃m應(yīng)調(diào)整，將有助于提高現(xiàn)階段數(shù)值計(jì)算的摩阻預(yù)示精度。

圖8 三角翼模型摩阻試驗(yàn)表面溫度分布(Ma∞=8，T0=749 K)Fig.8 Surface temperature distribution of delta-wing model in friction test (Ma∞=8，T0=749 K)

2.2.2 壁面溫度邊界條件的表征方式

目前，高速流動(dòng)問題的數(shù)值模擬通常采用等溫壁面，數(shù)值計(jì)算中2種常用的等溫壁面邊界條件表征方式如圖9所示，1代表壁面內(nèi)第1層網(wǎng)格單元中心，w代表固體壁面，-1代表固體壁面外第1層虛網(wǎng)格中心。方式1第1層虛網(wǎng)格中心的溫度T-1滿足：

圖9 2種常用的等溫壁面邊界條件處理方式Fig.9 Two common ways of isothermal wall boundary in simulation

T-1=2Tw-T1

(4)

方式2第1層虛網(wǎng)格中心的溫度T-1滿足：

T-1=Tw

(5)

圖10 高空高馬赫數(shù)2種等溫壁面處理方式飛行器摩阻系數(shù)對(duì)比Fig.10 Aircraft friction coefficients comparison between two ways of isothermal wall boundary at high altitude and high Mach number

2種處理方式對(duì)摩阻的影響可以通過雷諾比擬來解釋。對(duì)于不可壓平板層流流動(dòng)，可以通過雷諾比擬建立摩阻系數(shù)Cf和斯坦頓數(shù)St間的聯(lián)系[23]：

(6)

(7)

式中：Pr為普朗特?cái)?shù)，層流一般取0.7～0.72；ρe和ue分別為邊界層外緣的密度和速度；hw為壁面焓；haw為絕熱壁焓；qw為壁面熱流密度，其計(jì)算式為

(8)

其中：k為熱傳導(dǎo)率；T為流場溫度；n為壁面法向方向。在數(shù)值求解Navier-Stokes方程時(shí)，基于網(wǎng)格中心的有限體積方法的壁面溫度梯度的計(jì)算表達(dá)式為

(9)

其中：Δd為第1層網(wǎng)格中心與壁面的距離。由式(9)可知，圖9中的表征方式1可以保證壁面處溫度梯度不變，而表征方式2則會(huì)導(dǎo)致壁面處溫度梯度偏小。結(jié)合式(6)～式(8)可知，使用方式2 進(jìn)行等溫壁面的溫度條件處理會(huì)導(dǎo)致斯坦頓數(shù)St和摩阻系數(shù)Cf都相應(yīng)偏小。摩阻和熱流相關(guān)的計(jì)算主要涉及到速度和溫度梯度的精確計(jì)算，目前的CFD方法中多以一階離散為主，需要在現(xiàn)有基礎(chǔ)上發(fā)展流場變量梯度的高階表征方法。

3 摩阻精確預(yù)示的研究需求

通過本文對(duì)摩阻計(jì)算精度影響因素的分析發(fā)現(xiàn)，無黏通量離散格式的耗散特性和壁面溫度條件對(duì)高馬赫數(shù)流動(dòng)的摩阻計(jì)算有重要影響。面向未來的工程應(yīng)用需求，亟需依托于國家數(shù)值風(fēng)洞工程發(fā)展高精度數(shù)值方法以精確模擬邊界層黏性流動(dòng)及其與激波、分離相互干擾作用中的流動(dòng)問題，從而提高摩擦阻力的數(shù)值預(yù)示精度，具體需求包括：

1) 低耗散數(shù)值方法

計(jì)算的表面摩擦黏性應(yīng)力與邊界層近壁區(qū)數(shù)值方法的耗散特性密切相關(guān)，鑒于此，對(duì)整個(gè)流場可以采用不同區(qū)域不同耗散水平的處理方式——即分區(qū)低耗散數(shù)值方法，同時(shí)保證激波區(qū)域的計(jì)算穩(wěn)定性和近壁區(qū)域的低耗散性。以Roe格式為例，可以采用基于當(dāng)?shù)亓鲌鰠?shù)的自適應(yīng)熵修正方法，在近壁薄層內(nèi)減小甚至關(guān)閉熵修正。但當(dāng)邊界層流動(dòng)比較復(fù)雜，尤其流場變量梯度較大時(shí)，僅僅在近壁薄層內(nèi)保證數(shù)值低耗散并不足以保證摩阻計(jì)算的準(zhǔn)確性，需要發(fā)展高精度數(shù)值方法，提高邊界層流動(dòng)的整體模擬精度。

2) 基于邊界層當(dāng)?shù)亓鲌鲎兞康谋诿鏈囟茸赃m應(yīng)調(diào)整技術(shù)

高速飛行時(shí)，來流的總溫較高，氣動(dòng)加熱現(xiàn)象往往比較嚴(yán)重，但是飛行器表面熱流一般會(huì)隨不同部位表現(xiàn)出巨大差異性，導(dǎo)致不同部位的表面溫度差異明顯。而高速流動(dòng)時(shí)壁面溫度與流場的相互耦合影響大，導(dǎo)致計(jì)算的摩擦阻力表現(xiàn)出很強(qiáng)的壁面溫度相關(guān)性。現(xiàn)階段數(shù)值計(jì)算中一般使用等溫壁并將全局設(shè)為同一溫度，其無法反映壁面溫度的空間不均勻性。因此，需要根據(jù)壁面不同部位氣動(dòng)加熱的程度，在給定基準(zhǔn)壁面溫度的基礎(chǔ)上，發(fā)展可實(shí)現(xiàn)所有位置壁溫的高效自適應(yīng)調(diào)整技術(shù)。

3) 真實(shí)表面邊界的建模和數(shù)值模擬方法

首先，材料與工藝水平會(huì)導(dǎo)致真實(shí)飛行器表面并不光滑；此外，高馬赫數(shù)、長時(shí)間飛行時(shí)氣動(dòng)加熱現(xiàn)象顯著，飛行器表面材料出現(xiàn)熔解、燒蝕，產(chǎn)生質(zhì)量引射并形成粗糙表面，這些效應(yīng)與高速流動(dòng)固有的高溫、轉(zhuǎn)捩等復(fù)雜流動(dòng)效應(yīng)相互耦合，對(duì)飛行器流場及摩阻等產(chǎn)生顯著影響。因此，在重點(diǎn)發(fā)展的高精度數(shù)值方法方面，首先需要考慮層流狀態(tài)下的粗糙表面、質(zhì)量引射等邊界的建模及數(shù)值模擬方法問題。

4) 高精度摩阻數(shù)值計(jì)算方法的驗(yàn)證

目前，可用于摩阻計(jì)算精度驗(yàn)證的試驗(yàn)數(shù)據(jù)仍十分缺乏，相關(guān)地面測量技術(shù)尚不成熟。相關(guān)研究已經(jīng)獲取的少量地面試驗(yàn)數(shù)據(jù)其精度和可靠性不足以支撐高精度摩阻數(shù)值計(jì)算方法的驗(yàn)證與改進(jìn)。從計(jì)算方法的精度驗(yàn)證方面，需要開展能夠反映典型流動(dòng)特征的標(biāo)準(zhǔn)試驗(yàn)?zāi)Ｐ驮O(shè)計(jì)，系統(tǒng)開展精細(xì)化摩阻測量方法研究，并針對(duì)典型外形開展多種摩阻測量技術(shù)的對(duì)比驗(yàn)證風(fēng)洞試驗(yàn)，研究并探索飛行試驗(yàn)?zāi)ψ杩茖W(xué)測量技術(shù)，積累可靠的摩阻試驗(yàn)數(shù)據(jù)，為發(fā)展高精度數(shù)值模擬方法的驗(yàn)證提供數(shù)據(jù)支撐。

4 結(jié) 論

目前高馬赫數(shù)層流的摩阻數(shù)值計(jì)算，相比于風(fēng)洞試驗(yàn)測量結(jié)果仍然偏大，本文通過對(duì)無黏通量空間格式數(shù)值耗散和壁面溫度邊界條件對(duì)表面摩阻影響的計(jì)算和分析，得到了以下結(jié)論：

1) NNW-Flowstar在高馬赫數(shù)(Ma=8～10)范圍內(nèi)，摩阻的計(jì)算精度與常用CFD軟件相仿。

2) 表面黏性摩擦應(yīng)力的計(jì)算精度與近壁區(qū)空間格式的耗散密切相關(guān)，數(shù)值耗散越小，表面摩阻的計(jì)算精度越高；在速度較低的邊界層近壁區(qū)內(nèi)關(guān)閉熵修正，將有助于提高表面摩阻的預(yù)示精度。

3) 高馬赫數(shù)流動(dòng)不同部位的壁溫變化明顯，進(jìn)行數(shù)值模擬時(shí)壁面溫度邊界條件對(duì)表面摩阻的計(jì)算有重要影響。

4) 結(jié)合工程需要提出了高精度摩阻數(shù)值預(yù)示的研究需求，主要包括低耗散數(shù)值方法、基于邊界層當(dāng)?shù)亓鲌鲎兞康谋诿鏈囟茸赃m應(yīng)調(diào)整技術(shù)、真實(shí)表面邊界的建模和數(shù)值模擬方法、高精度摩阻數(shù)值計(jì)算方法的驗(yàn)證等。

致 謝

感謝國家數(shù)值風(fēng)洞工程提供的網(wǎng)格劃分軟件NNW-Gridstar和數(shù)值計(jì)算軟件NNW-Flowstar。