高馬赫數層流摩阻數值計算精度

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

中國運載火箭技術研究院 空間物理重點實驗室，北京 100076

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

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

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

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

1 摩阻計算與試驗測量間的差異現象

1.1 圓錐表面摩阻計算與試驗對比

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

圖1 Meritt試驗模型及傳感器位置示意圖[11]Fig.1 Sketch of cone model and sensor location of Meritt[11]

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

表1 圓錐模型層流狀態風洞試驗參數

層流狀態下測點處摩擦阻力系數數值計算結果與試驗的對比如圖2所示，可以看到，Flowstar的模擬精度與其他CFD軟件接近，攻角α=10°時Flowstar相比于CFD++和HyperCFD分別偏大約4%和7%，不同軟件的計算結果與試驗測量結果偏差較大，相比試驗偏大24%～30%。

圖2 圓錐模型表面摩阻系數數值計算與風洞試驗[11]差異Fig.2 Difference of skin friction coefficients of cone model between numerical calculation and wind tunnel test[11]

1.2 三角翼表面摩阻計算與試驗對比

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

表2 三角翼模型風洞試驗參數Table 2 With tunnel test parameters of delta-wing model

分別使用Flowstar、CFD++和HyperCFD對上述試驗工況進行數值模擬。其中，Flowstar使用基于T-Rex技術的非結構網格，邊界層第1層網格高度為0.01 mm；CFD++和HyperCFD同時采用2套結構化網格進行計算以評估網格密度對摩阻計算精度的影響：粗網格流向、法向和周向網格點數分別為147、81、327，邊界層第1層高度為0.01 mm；密網格流向、法向和周向網格點數分別為301、301、327，邊界層第1層高度為0.005 mm。壁面溫度均設為300 K。

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

圖3 Ma∞=8、0°攻角下表面中心線摩阻系數層流計算結果與試驗對比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個計算與試驗結果的對比，可以發現，即便不考慮流動轉捩的影響，高馬赫數層流的摩擦阻力計算仍然整體上高于試驗測量結果。

2 摩阻數值計算精度的影響因素

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

2.1 數值耗散的影響

使用自研的HyperCFD結構網格求解器對表1所列的層流風洞試驗狀態進行數值模擬，對比不同空間格式和熵修正方法對摩阻計算的影響，壁面條件為等溫壁，壁溫為300 K。

2.1.1 無黏通量格式

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

(1)

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

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

圖4 不同空間格式測點摩阻系數對比Fig.4 Comparison of measurement location friction coefficients among different space schemes

2.1.2 熵修正方法

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

(2)

式中：δ*為熵修正比例系數；V為速度矢量；a為當地聲速;ξ、η、ζ為計算坐標系的3個方向。Müller熵修正方法δ的表達式為(以ξ方向為例)

(3)

圖5給出了圓錐模型不同熵修正方法和熵修正比例系數測點摩阻系數的計算對比。可以看到，2種熵修正方法的測點摩阻系數計算結果相差4%～8%，Müller型熵修正對黏性邊界層的模擬精度更高。相比于Harten-Yee熵修正，Müller熵修正方法考慮了速度和網格的各向異性，對網格長細比較大的邊界層區域具有更高的分辨率。此外，2種熵修正方法均表現為熵修正比例系數越大，Roe格式的數值耗散越大，計算的摩阻系數也相應越大。

從圖5中可以得到的結論是，表面摩阻計算結果與格式的數值耗散特性密切相關，對于Roe格式而言，熵修正比例系數越小，計算摩阻系數越小。基于此，本文嘗試在邊界層內近壁區關閉熵修正，圖6給出了這種處理方式下10°攻角時圓錐模型迎風面子午線(迎風面中心線)的摩阻系數，并與常規熵修正(Harten-Yee eps0.1)的摩阻結果進行了對比。

圖5 不同熵修正方法測點摩阻系數對比Fig.5 Comparison of measurement location friction coefficients among different entropy fix methods

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

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

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

2.2.1 壁面溫度

從低速到低超聲速流動，來流總溫較低，飛行器表面溫度變化不明顯，壁溫對摩阻計算的影響問題不突出；而高馬赫數飛行時，氣動加熱現象顯著，計算時壁面溫度的選取對于摩阻計算會產生較大影響。

第一，根據皖河流域山區環境的特點，構建低耗、優質、高產、高效的農田生態系統。主要內容有：山區土地整治和土壤改良，保持耕地的綠色覆蓋，建設生態水系和現代灌溉系統，山坡耕地保護或退耕等。

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

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

0°攻角時，邊界層流動相對簡單，壁面摩阻對溫度邊界條件不敏感，流向位置60 mm之前，絕熱壁的當地摩阻更高；隨著邊界層發展，流向位置200 mm之后，絕熱壁的當地摩阻低于等溫壁。

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

圖8給出了三角翼模型在Ma∞=8來流中進行摩阻測量時紅外測得的表面溫度場，局部區域(尤其是頭部和側緣)氣動加熱明顯。考慮到壁面溫度設置對摩阻計算結果的影響，如果CFD軟件能夠實現基于邊界層當地流場變量的壁面溫度自適應調整，將有助于提高現階段數值計算的摩阻預示精度。

圖8 三角翼模型摩阻試驗表面溫度分布(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 壁面溫度邊界條件的表征方式

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

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

T-1=2Tw-T1

(4)

方式2第1層虛網格中心的溫度T-1滿足：

T-1=Tw

(5)

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

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

(6)

(7)

式中：Pr為普朗特數，層流一般取0.7～0.72；ρe和ue分別為邊界層外緣的密度和速度；hw為壁面焓；haw為絕熱壁焓；qw為壁面熱流密度，其計算式為

(8)

其中：k為熱傳導率；T為流場溫度；n為壁面法向方向。在數值求解Navier-Stokes方程時，基于網格中心的有限體積方法的壁面溫度梯度的計算表達式為

(9)

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

3 摩阻精確預示的研究需求

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

1) 低耗散數值方法

計算的表面摩擦黏性應力與邊界層近壁區數值方法的耗散特性密切相關，鑒于此，對整個流場可以采用不同區域不同耗散水平的處理方式——即分區低耗散數值方法，同時保證激波區域的計算穩定性和近壁區域的低耗散性。以Roe格式為例，可以采用基于當地流場參數的自適應熵修正方法，在近壁薄層內減小甚至關閉熵修正。但當邊界層流動比較復雜，尤其流場變量梯度較大時，僅僅在近壁薄層內保證數值低耗散并不足以保證摩阻計算的準確性，需要發展高精度數值方法，提高邊界層流動的整體模擬精度。

2) 基于邊界層當地流場變量的壁面溫度自適應調整技術

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

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

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

4) 高精度摩阻數值計算方法的驗證

目前，可用于摩阻計算精度驗證的試驗數據仍十分缺乏，相關地面測量技術尚不成熟。相關研究已經獲取的少量地面試驗數據其精度和可靠性不足以支撐高精度摩阻數值計算方法的驗證與改進。從計算方法的精度驗證方面，需要開展能夠反映典型流動特征的標準試驗模型設計，系統開展精細化摩阻測量方法研究，并針對典型外形開展多種摩阻測量技術的對比驗證風洞試驗，研究并探索飛行試驗摩阻科學測量技術，積累可靠的摩阻試驗數據，為發展高精度數值模擬方法的驗證提供數據支撐。

4 結 論

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

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

2) 表面黏性摩擦應力的計算精度與近壁區空間格式的耗散密切相關，數值耗散越小，表面摩阻的計算精度越高；在速度較低的邊界層近壁區內關閉熵修正，將有助于提高表面摩阻的預示精度。

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

4) 結合工程需要提出了高精度摩阻數值預示的研究需求，主要包括低耗散數值方法、基于邊界層當地流場變量的壁面溫度自適應調整技術、真實表面邊界的建模和數值模擬方法、高精度摩阻數值計算方法的驗證等。

致 謝

感謝國家數值風洞工程提供的網格劃分軟件NNW-Gridstar和數值計算軟件NNW-Flowstar。