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

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

式中，在ξ、η、ζ三個(gè)方向均采用Van Leer通量矢量分裂方法，以ξ方向?yàn)槔?/p>

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

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

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

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

2 計(jì)算方法驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

3 計(jì)算模型與工況條件

本文以文獻(xiàn)[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，則對壁溫進(jìn)行無量綱化（除以來流溫度）可得壁溫條件（Tw/T∞）分別為1.124 5、2.024、2.834、3.643、4.453、5.263、6.072。

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

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

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

4 計(jì)算結(jié)果分析

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

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

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

圖7 不同壁面溫度下分離點(diǎn)、再附點(diǎn)位置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），分離點(diǎn)向前移動(dòng)、再附點(diǎn)向后移動(dòng)，干擾區(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所示,從圖中可以看到：在分離點(diǎn)之前，熱流沿著來流方向降低，同時(shí)隨著壁面溫度升高而降低；在分離點(diǎn)之后，熱流出現(xiàn)了急劇的下降，且壁面溫度越高，熱流下降的越緩慢；在干擾區(qū)內(nèi)，熱流逐漸降至最低且在一段區(qū)域內(nèi)基本保持不變；在第二級壓縮面，熱流開始逐漸增大。

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

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

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

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

圖13 各位置取點(diǎn)示意圖Fig. 13 Schematic of different points

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

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

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

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

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

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

式中，Qw為當(dāng)?shù)乇诿鏌崃?，Q300K為300 K條件下當(dāng)?shù)乇诿鏌崃?，Hre為當(dāng)?shù)貧饬骰謴?fù)焓，Hw為當(dāng)?shù)乇诿骒剩琀300K為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是邊界層外緣速度；普朗特?cái)?shù)Pr≈0.71。

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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