壁溫對壓縮拐角流動影響的數(shù)值模擬研究

劉 亮，邱 波，曾 磊，姚 杰，朱目成

(1. 西南科技大學(xué) 制造科學(xué)與工程學(xué)院，綿陽 621010；2. 中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引 言

高超聲速飛行器氣動布局中普遍存在一級或多級壓縮拐角結(jié)構(gòu)，如進氣道、襟翼等。拐角區(qū)域存在分離再附流動、激波與邊界層干擾等復(fù)雜現(xiàn)象，對局部熱環(huán)境的分布會產(chǎn)生較大影響。針對這一結(jié)構(gòu)，國內(nèi)外學(xué)者進行了大量的研究。Davis和Sturvant[1]開展了非平衡真實氣體對壓縮拐角分離長度的研究。童福林[2]采用直接數(shù)值模擬（DNS）方法研究激波對邊界層的影響，結(jié)果表明分離氣泡上方的剪切層對分離流動有顯著影響。Simeonides和Vermeulen[3]提供了二維壓縮角中完全層流和過渡激波/附面層相互作用的實驗數(shù)據(jù)，并對兩個Navier-Stokes解算器進行了驗證。John[4]研究表明前緣鈍度會減小二維流場的分離區(qū)，增大軸對稱流場的分離區(qū)。李素循[5]的研究表明楔面上的壓力峰值隨著楔角的增加迅速增大，且壓力峰值的位置慢慢向楔面根部方向移動。陳蘇宇[6]、趙一龍[7]對相關(guān)研究進行了總結(jié)和進一步研究，表明層流邊界層在側(cè)板激波作用下首先發(fā)生轉(zhuǎn)捩，然后分離，分離后的流動為顯著的半錐形流場結(jié)構(gòu)，雷諾數(shù)對峰值熱流大小和轉(zhuǎn)捩過程影響明顯。

總的來說，上述研究對于壓縮拐角的流動特征已經(jīng)有了較為全面的認識，但大多是在300 K均勻“冷壁”條件下進行的。而在真實的高超聲速飛行中，氣動加熱會使得飛行器表面溫度升高，壁面溫度對于壓縮拐角的影響顯然不可忽略。

21世紀(jì)初，Marini[8]在壓縮拐角的試驗研究中考慮到了壁溫比的變化。2006年，德國的Thomas和Herbert[9]針對雙楔模型開展了風(fēng)洞試驗，選取了三種前緣鈍度以及三種壁溫進行研究，獲取不同前緣鈍度、不同壁溫下的壓力和熱流數(shù)據(jù)。隨后，Reinartz[10]等研究表明壁溫會影響流動分離發(fā)生的位置。在此基礎(chǔ)上，2014年，尚慶[11]等研究表明壁溫會影響鈍雙楔轉(zhuǎn)的轉(zhuǎn)捩與流動分離，按第一級壓縮面層流且第二級壓縮面湍流進行數(shù)值模擬的結(jié)果能夠在一定程度上與Thomas等的試驗結(jié)果相符。Bleilebens和Olivier[12]采用可預(yù)熱雙楔模型進行實驗研究，壓力和熱流測量的結(jié)果表明，主要影響因素不是壁面溫度本身，而是壁溫與自由來流溫度的比值。2016年，代光月[13]對一薄殼狀兩級壓縮楔進行了風(fēng)洞試驗研究，并采用多場耦合計算方法進行了驗證分析，通過熱壁修正公式將壁面溫度作為流固耦合的中間邊界條件。

上述研究讓我們了解了壁面溫度對壓縮拐角影響的基本規(guī)律，但在溫度變化范圍、壁溫對于流場結(jié)構(gòu)的影響等方面仍有很多工作需要進一步深入，特別是在目前廣泛使用的熱壁修正公式的適用范圍方面鮮有文獻進行研究。針對該問題，在前人研究工作的基礎(chǔ)上，本文選用文獻[13]中的鈍雙楔模型并設(shè)計了多個計算工況，采用基于Navier-Stokes方程的自研程序進行數(shù)值模擬，在較大溫度范圍內(nèi)分析了壁面溫度對壓縮拐角的流場結(jié)構(gòu)、熱流分布的影響。此外，本文使用熱壁修正公式對無干擾區(qū)域和干擾區(qū)域的熱流進行修正，通過與變壁溫直接計算熱流進行對比，初步分析了該公式對于壓縮拐角流動的適用性。

1 數(shù)值計算方法簡介

在笛卡爾坐標(biāo)系下，守恒形式的三維非定?？蓧嚎sNavier-Stokes方程無量綱化后可寫成如下形式：

其中，Re是無量綱雷諾數(shù)，Q是守恒狀態(tài)變量，E、F、G是無黏通量向量，Ev、Fv、Gv是黏性通量向量?；谕耆珰怏w假設(shè)，取Sutherland公式所給出的黏性系數(shù)，普朗特數(shù)Pr= 0.72，空氣的比熱比γ= 1.4。

引入 計 算空 間 (ξ,η,ζ)， 使 其與 物 理空 間 (x,y,z)存在唯一的單值坐標(biāo)變換，形式如下：

變換后Navier-Stokes方程可表述為：

本文數(shù)值計算采取有限體積法，將上述方程轉(zhuǎn)化為積分形式：

式中，f是封閉曲面S上的通量矢量，V是S包圍的體積，n是S的單位法向向量。在網(wǎng)格線ξ、η、ζ包圍的網(wǎng)格單元內(nèi)對式（4）積分，可得半離散化方程，如式（5）所示：

式中，在ξ、η、ζ三個方向均采用Van Leer通量矢量分裂方法，以ξ方向為例：

式中，采用MUSCL方法在交界面處插值，限制器取為Van Albada限制器。

對于高雷諾數(shù)的定常黏性流動，為了準(zhǔn)確模擬邊界層，物面附近必須使用很小的網(wǎng)格間距。為避免因顯式格式時間步長過小所導(dǎo)致的收斂時間的大量花費[14-15]，本文中的定常Navier-Stokes方程采用Yoon提出的隱式LU-SGS方法。

湍流選用k-ωSST（Shear-Stress Transport）模型[16-17]，其控制方程為：

式中，k、ω分別代表湍動能和湍動能的比耗散率，μt為湍流渦黏性系數(shù)，β*為模型方程中相關(guān)常數(shù)，Φ為湍動能的產(chǎn)生項。

2 計算方法驗證

為了驗證本文計算方法對壓縮拐角模擬的有效性，選用了Thomas和Herbert[9]的風(fēng)洞試驗?zāi)Ｐ秃蛠砹鳁l件，對比分析中同時加入了Reinartz[10]和尚慶[11]對于該試驗的計算結(jié)果。模型結(jié)構(gòu)和尺寸如圖1（a）所示；計算網(wǎng)格如圖1（b）所示，第一層網(wǎng)格高度y+=0.002 mm，物面周向布點451個，物面法向布點91個；計算來流條件為Ma∞= 8.1，T∞= 106 K，Re∞= 3.8×106/m，p∞= 520 Pa；第一壓縮面采用層流、第二壓縮面采用湍流計算。

圖1 Thomas和Herbert[9]風(fēng)洞試驗?zāi)Ｐ图熬W(wǎng)格Fig. 1 Thomas and Herbert[9] wind tunnel test model and Grids

本文計算流場結(jié)構(gòu)與文獻[9]試驗紋影照片的對比如圖2所示，分離點和再附點的位置、激波高度及激波角度的對比如表1所示?？梢钥吹剑m然本文計算結(jié)果在激波高度和激波角與試驗結(jié)果略有差異，但整體結(jié)構(gòu)與試驗結(jié)果吻合良好。

圖2 壓縮拐角分離區(qū)流場結(jié)構(gòu)對比Fig. 2 Flow structure in separation zone of compression corner

表1 與文獻[9]流場參數(shù)對比Table 1 Validation of Flow Parameters with Ref. [9]

圖3是迎風(fēng)中心線上的壓力系數(shù)對比結(jié)果，由該圖可以看出，本文計算的壓力系數(shù)與文獻[10]和文獻[11]數(shù)值計算結(jié)果對比良好，且與風(fēng)洞試驗結(jié)果的吻合度較高。

按公式將熱流轉(zhuǎn)換熱流斯坦頓數(shù)St，再與文獻[10, 11]進行對比，公式[18]表示為：

式中，Q為當(dāng)?shù)責(zé)崃?，Tw為當(dāng)?shù)乇诿鏈囟?，T0為來流總溫，Cp∞為 來流壓力系數(shù)，ρ∞為 來流密度，u∞為來流速度。

圖3 壓力系數(shù)對比驗證Fig. 3 Validation of pressure coefficient distribution

迎風(fēng)中心線的熱流斯坦頓數(shù)結(jié)果如圖4所示。由該圖可以看出，本文計算的熱流斯坦頓數(shù)與文獻[10, 11]數(shù)值計算結(jié)果對比良好，但計算結(jié)果均比試驗結(jié)果偏大，可能是由于第一壓縮面采用層流、第二壓縮面采用湍流模型，且不同計算方法對壓縮拐角流動的數(shù)值模擬能力不同所導(dǎo)致的。

圖4 直接數(shù)值模擬的熱流斯坦頓數(shù)結(jié)果對比Fig. 4 Validation of St number of DNS results

綜合本文計算結(jié)果與文獻在流場結(jié)構(gòu)、壓力系數(shù)和熱流斯坦頓數(shù)三方面的對比表現(xiàn)，可以看出本文所用計算方法對壓縮拐角的數(shù)值模擬結(jié)果具有較高的可信度。

3 計算模型與工況條件

本文以文獻[13]中的鈍雙楔形模型作為參考，取二維模型，模型包含兩級壓縮面，總長603.4 mm，前緣半徑為3 mm，第一級壓縮面水平長度445.93 mm,與下壁面夾角7.34°，第二級壓縮面水平長度97.47 mm，與下壁面夾角25.52°，肩部長度60 mm，如圖5所示。來流條件為：高度60 km、馬赫數(shù)9、迎角-10°，壁溫條件分別取300、500、700、900、1 100、1 300、1 500 K，來流溫度T∞=247.021 K，則對壁溫進行無量綱化（除以來流溫度）可得壁溫條件（Tw/T∞）分別為1.124 5、2.024、2.834、3.643、4.453、5.263、6.072。

圖5 計算模型(單位:mm)Fig. 5 DNS configuration (unit: mm)

計算網(wǎng)格如圖6所示，第一層網(wǎng)格高度y+= 0.1 mm，物面周向布點581個，物面法向布點111個，網(wǎng)格雷諾數(shù)為8。在盡量保證網(wǎng)格尺度均勻、過渡光滑、物面正交性良好的基礎(chǔ)上，在壓縮拐角區(qū)域進行了適當(dāng)?shù)募用芴幚怼?/p>

圖6 計算網(wǎng)格Fig. 6 Computational grid

4 計算結(jié)果分析

4.1 壁面溫度對分離渦結(jié)構(gòu)的影響

高速來流經(jīng)過壓縮拐角處，由于逆壓梯度的存在會形成分離渦，進而產(chǎn)生分離流動和再附流動[19]。壁面溫度的變化會引起分離渦的變化，進而影響壓縮拐角的熱流分布[20]。

圖7給出了分離點、再附點的具體變化情況，圖8給出了部分壁溫下干擾區(qū)的分離渦結(jié)構(gòu)?？梢钥闯?，不同壁溫下流場結(jié)果并未發(fā)生本質(zhì)改變，但分離點、再附點的位置以及分離渦的大小發(fā)生了變化：隨著壁面溫度升高，壁溫比增大，分離點向前移動、再附點向后移動，分離渦增大。分離渦增大。

圖7 不同壁面溫度下分離點、再附點位置Fig. 7 Separation and reattachment points at different wall temperatures

圖8 不同壁面溫度下分離渦結(jié)構(gòu)與大小Fig. 8 Structure and size of separated vortices at different wall temperatures

圖9是不同壁面溫度下近壁面區(qū)域的氣體密度、黏性系數(shù)、近壁面馬赫數(shù)和壓力變化圖，通過對比可以看出：在高速來流條件下，壁面溫度升高，壁溫比增大，分離渦的整體流場結(jié)構(gòu)基本不變，但近壁面區(qū)域的流體密度降低（圖9（a））、黏性系數(shù)增大（圖9（b））、速度邊界層變厚（圖9（c）），導(dǎo)致了逆壓力梯度減?。▓D9（d）、圖10），分離點向前移動、再附點向后移動，干擾區(qū)域增大。

圖9 不同壁面溫度下近壁面氣體密度、黏性系數(shù)、壓力和馬赫數(shù)變化圖Fig. 9 Variation of near wall gas density, viscosity coefficient, pressure and Mach number at different wall temperatures

圖10 不同壁面溫度下的壓力系數(shù)對比Fig. 10 Comparison of pressure coefficient under different wall temperature

4.2 壁面溫度對熱流的影響

迎風(fēng)中心線的熱流分布如圖11所示,從圖中可以看到：在分離點之前，熱流沿著來流方向降低，同時隨著壁面溫度升高而降低；在分離點之后，熱流出現(xiàn)了急劇的下降，且壁面溫度越高，熱流下降的越緩慢；在干擾區(qū)內(nèi)，熱流逐漸降至最低且在一段區(qū)域內(nèi)基本保持不變；在第二級壓縮面，熱流開始逐漸增大。

圖11 迎風(fēng)中心線熱流密度分布Fig. 11 Heat flux distribution at windward centerline

統(tǒng)計不同壁面溫度下分離點和再附點的熱流，并以該壁面溫度下的駐點熱流為基準(zhǔn)進行無量綱化處理，結(jié)果如圖12所示。由圖可知，隨著壁面溫度升高，壁溫比增大，分離點無量綱熱流增大，再附點無量綱熱流減小。綱熱流減小。

圖12 不同壁面溫度下分離點、再附點的無量綱熱流Fig. 12 Dimensionless heat flux at separation and reattachment points at different wall temperatures

為了進一步量化壁面溫度對壓縮拐角不同位置熱環(huán)境的影響，接下來分別取駐點、a點x= 150 mm、b點x= 315 mm、c點x= 360 mm、d點x= 480 mm、e點x= 541 mm等幾個典型位置來分析熱流的變化，如圖13所示，其中a、e兩點位于無干擾平板上，b、c、d三個點處于分離點和再附點之間的干擾區(qū)范圍內(nèi)。

圖13 各位置取點示意圖Fig. 13 Schematic of different points

統(tǒng)計這些點的熱流，如表2所示。以各壁面溫度下的駐點熱流值為參考量將各個點的熱流進行無量綱化處理，結(jié)果如圖14所示?？梢钥吹剑v點、a、b、d、e各點的熱流均是隨著壁面溫度的升高而減小，其中b點熱流減小的幅度最大，但c點的熱流隨著壁面溫度的升高而逐漸增大。

由此可看出，對于本文計算的壓縮拐角模型，在駐點和無干擾區(qū)域，熱流隨著壁面溫度的升高而降低；在干擾區(qū)內(nèi)，大部分區(qū)域的熱流隨著壁面溫度升高而減小，且減小幅度比無干擾區(qū)更大，但在分離點與拐角之間的部分區(qū)域熱流會隨壁面溫度的升高而增大，這主要是由于分離點前移、干擾區(qū)增大，該部分逐漸遠離分離點所導(dǎo)致的。

表2 不同壁面溫度下各個點的熱流密度值(單位: kW/m2)Table 2 Heat flux density at different wall temperatures (unit: kW/m2)

圖14 幾個典型位置無量綱熱流隨壁面溫度變化Fig. 14 Dimensionless heat flux varies with wall temperature in several typical positions

4.3 熱壁修正公式適用性分析

Chen等[21]基于高超聲速邊界層理論發(fā)展了一種高焓條件下熱壁修正方法，該方法通過恢復(fù)焓和壁面焓對壁面熱流進行修正，較大地簡化了對邊界層各物性參數(shù)的討論，得到了廣泛的應(yīng)用。熱壁修正公式可表示為[22-23]：

式中，Qw為當(dāng)?shù)乇诿鏌崃鳎琎300K為300 K條件下當(dāng)?shù)乇诿鏌崃?，Hre為當(dāng)?shù)貧饬骰謴?fù)焓，Hw為當(dāng)?shù)乇诿骒?，H300K為300 K條件下當(dāng)?shù)乇诿骒省?/p>

恢復(fù)焓Hre表示為[18]：

轉(zhuǎn)化為：

式中，He邊界層外緣氣流焓；r0是恢復(fù)因子，層流狀態(tài)下取0.842 6，湍流狀態(tài)下取0.892 1；ue是邊界層外緣速度；普朗特數(shù)Pr≈0.71。

綜合以上分析，人們對大量工程實踐經(jīng)驗進行總結(jié)，將恢復(fù)焓表示為與總焓的一個關(guān)系式，如下所示：

式中，r為恢復(fù)焓系數(shù)，層流狀態(tài)下一般取0.89，湍流狀態(tài)下一般取0.92，H0為總焓。

這里利用公式（10）將表2中各點的高壁溫?zé)崃餍拚?00 K壁溫下的熱流，前緣駐點熱流及無干擾平板熱流的修正結(jié)果如圖15。由該圖可以看出，修正后的熱流與300 K的計算熱流吻合較好，駐點最大修正誤差在-1%以內(nèi)，無干擾平板上a點和e點最大修正誤差-5%左右，均在合理的范圍之內(nèi)。由此可見，對于本文所計算的壓縮拐角模型，熱壁修正公式對于前緣駐點和無干擾平板區(qū)域具有較好的修正作用。

圖15 駐點及平板上熱流修正結(jié)果Fig. 15 Correction results of heat flux on stationary points and flat plates

干擾區(qū)熱流的修正結(jié)果如圖16所示，修正后b點最大誤差-41%，d點最大誤差-33%，c點最大誤差60%。由此可見，熱壁修正公式對干擾區(qū)的熱流也具有一定程度的修正作用，但效果非常有限，導(dǎo)致修正結(jié)果與計算結(jié)果之間存在較大偏差。

圖16 干擾區(qū)熱流修正結(jié)果Fig. 16 Correction results of heat flux in interference zone

綜合來看，對于本文所計算的壓縮拐角模型，其在圓柱前緣部分和無干擾平面區(qū)域較為適用，但其在拐角處干擾區(qū)的直接使用則存在修正精度降低、適用性不足的問題。

5 結(jié) 論

本文采用數(shù)值模擬的方法，在較大溫度范圍內(nèi)(300 K～1 500 K)分析了壁溫對高超聲速壓縮拐角流場結(jié)構(gòu)和熱流分布的影響規(guī)律，并探究了熱壁修正公式在該結(jié)構(gòu)上的適用性。針對本文所計算的壓縮拐角模型，可以得出如下結(jié)論：

1）壁溫升高，壁溫比增大，近壁面區(qū)域的流體物性參數(shù)變化較大，導(dǎo)致逆壓力梯度減小，拐角處分離渦增大。

2）壓縮拐角大部分區(qū)域的熱流隨著壁溫升高而減小，但熱流并不完全遵循隨壁溫升高而減小的規(guī)律，在分離點與拐角之間的部分區(qū)域，熱流會隨壁溫的升高而增大，這主要是由分離渦增大、分離點位置前移導(dǎo)致的。

3）熱壁修正公式在圓柱前緣部分和無干擾平面區(qū)域較為適用，但其在拐角處干擾區(qū)則存在修正精度降低、適用性不足的問題。