劉景源
(南昌航空大學(xué)飛行器工程學(xué)院,江西 南昌 330063)
高超聲速戰(zhàn)略、戰(zhàn)術(shù)飛行器設(shè)計(jì)是世界各航天航空大國(guó)所瞄準(zhǔn)的一個(gè)重點(diǎn)研究方向[1-2]。為了提高高超聲速導(dǎo)彈、動(dòng)能打擊等武器的機(jī)動(dòng)性,舵面偏角會(huì)在較大范圍內(nèi)發(fā)生變化,因此準(zhǔn)確預(yù)估舵面偏角的氣動(dòng)力、氣動(dòng)熱對(duì)高超聲速武器系統(tǒng)的研制具有重要意義[3]。
應(yīng)用基于雷諾平均的渦粘湍流模型仍然是高超聲速飛行器壓縮拐角類繞流的氣動(dòng)力、氣動(dòng)熱預(yù)估的設(shè)計(jì)工具。湍流熱通量的建模精度是制約氣動(dòng)熱準(zhǔn)確預(yù)估的一個(gè)重要因素[4-5]。
對(duì)湍流熱通量的建模有零方程、一方程、二方程模型,以及直接建立湍流熱湍流的微分方程模型等4種類型。零方程模型相比其余3種模型,因?yàn)闊o(wú)需建立微分方程,具有建模簡(jiǎn)單、計(jì)算效率高、魯棒性強(qiáng)等優(yōu)點(diǎn),在高超聲速武器的氣動(dòng)特性預(yù)估中得到了廣泛的應(yīng)用[6-7]。由于零方程的梯度模型是基于雷諾比擬的假設(shè)建立的,比擬系數(shù)為湍流普朗特?cái)?shù),并且假設(shè)湍流普朗特?cái)?shù)是常數(shù),對(duì)高超聲速激波/湍流邊界層干擾導(dǎo)致的復(fù)雜流動(dòng)的精確模擬存在缺陷[8-9]。文獻(xiàn)[10]提出計(jì)及湍流渦粘性與分子粘性系數(shù)比的變湍流普朗特?cái)?shù)模型,但是該模型僅僅對(duì)不可壓縮流動(dòng)進(jìn)行的。文獻(xiàn)[6,11]發(fā)展了一種湍流普朗特?cái)?shù)代數(shù)模型,高超聲速尖錐、雙楔、HIFiRE-1等外形的氣動(dòng)熱計(jì)算結(jié)果表明,與湍流普朗特?cái)?shù)為常數(shù)的數(shù)值結(jié)果相比,所發(fā)展的模型給出的表面熱流更精確。而文獻(xiàn)[4,5,9]提出的零方程湍流熱通量模型,均基于計(jì)及激波非定常效應(yīng),并在修正了的k-ω模型基礎(chǔ)上進(jìn)行的。由于該k-ω模型眾所周知的缺陷,所提出的湍流熱通量模型聯(lián)合其他渦粘及雷諾應(yīng)力模型的模擬精度需要進(jìn)一步評(píng)估。
對(duì)高超聲速激波/湍流邊界層干擾導(dǎo)致的非平衡流動(dòng),湍流熱通量建模中的湍流普朗特?cái)?shù)并非常數(shù),而應(yīng)為非平衡流動(dòng)特征參數(shù)的函數(shù),從而提出了一種修正的湍流普朗特?cái)?shù)熱通量模型。應(yīng)用數(shù)值模擬及理論分析方法,對(duì)高超聲速平板、壓縮拐角(高超聲速飛行器舵面的簡(jiǎn)化外形)等繞流的數(shù)值分析,評(píng)估了所提出的零方程湍流熱通量模型。
基于RANS方程組及二方程SSTk-ω湍流模型方程組[12]進(jìn)行零方程湍流熱通量模型修正及數(shù)值分析。數(shù)值模擬中,RANS及SSTk-ω方程組的對(duì)流項(xiàng)采用具有低耗散特性的經(jīng)熵修正的二階TVD格式,RANS的黏性項(xiàng),及SSTk-ω方程組的擴(kuò)散項(xiàng),湍動(dòng)能生成項(xiàng)等均采用二階中心差分格式進(jìn)行離散。湍流方程組的負(fù)源項(xiàng)則采用點(diǎn)隱格式以提高計(jì)算效率。具體的數(shù)值方法、邊界條件,以及對(duì)數(shù)值方法的驗(yàn)證等見文獻(xiàn)[13,14]。
文獻(xiàn)[4-9]均表明湍流普朗特?cái)?shù)為常數(shù)的零方程湍流熱通量模型,在模擬高超聲速飛行器激波/湍流邊界層干擾等強(qiáng)非平衡湍流繞流的表面熱流存在缺陷,因此需要發(fā)展計(jì)及高超聲速非平衡流動(dòng)特性的變湍流普朗特?cái)?shù)模型。基于衡量高超聲速非平衡特性的參數(shù)——湍動(dòng)能生成與湍流耗散率比值的無(wú)量綱數(shù),提出一種湍流普朗特?cái)?shù)修正模型如下:
(1)
式中:Pk為湍動(dòng)能生成項(xiàng);ε為湍動(dòng)能耗散率[15]。
對(duì)于SSTk-ω模型,湍動(dòng)能耗散率ε可表示為:
ε=cμρωk+δ
(2)
其中δ是為避免式(1)中的分母為零而出現(xiàn)的小量,取δ=1×10-30。
文獻(xiàn)[15]采用限制器,對(duì)湍動(dòng)能生成與湍流耗散率比值的大小進(jìn)行限制,避免湍動(dòng)能生成項(xiàng)的非物理增長(zhǎng)。式(1)亦采用上述方法對(duì)湍流普朗特?cái)?shù)的大小進(jìn)行限制。由文中分析結(jié)果可見,所采用的限制方法能夠提高湍流模型預(yù)測(cè)精度。
由式(1)可見,當(dāng)非平衡無(wú)量綱參數(shù)(湍動(dòng)能生成與湍流耗散率比值)的比值小于1.2時(shí),即為湍流原湍流普朗特常數(shù)0.9;而當(dāng)非平衡參數(shù)大于1.2時(shí),湍流普朗特?cái)?shù)變小,極小值為0.2。湍流普朗特?cái)?shù)的極小值取值對(duì)高超聲速激波/湍流邊界層干擾導(dǎo)致的局部高熱流具有重要的影響,文獻(xiàn)[16]基于提出的二方程湍流熱通量模型,計(jì)算的高超聲速激波/湍流邊界層干擾流動(dòng)的表面熱流結(jié)果表明,湍流普朗特?cái)?shù)最小可取0.2,因此基于文獻(xiàn)[16],文中給出的最小湍流普朗特?cái)?shù)為0.2。
為了驗(yàn)證所提出的湍流熱通量的普朗特?cái)?shù)修正方法,下面給出了高超聲速平板以及15°,30°,34°壓縮拐角等4個(gè)算例的驗(yàn)證結(jié)果。其中網(wǎng)格、數(shù)值方法、邊界條件及驗(yàn)證等見文獻(xiàn)[13,14]。
選取文獻(xiàn)[17]給出的高超聲速平板繞流的實(shí)驗(yàn)測(cè)量結(jié)果,進(jìn)行修正湍流普朗特?cái)?shù)模型的驗(yàn)證。實(shí)驗(yàn)的來(lái)流馬赫數(shù)為8.18,來(lái)流溫度81 K,基于來(lái)流條件的每米雷諾數(shù)為4.9×106m-1,平板壁面溫度為300 K。
圖1給出了原SSTk-ω湍流模型、在原SSTk-ω湍流模型基礎(chǔ)上分別進(jìn)行湍流普朗特修正的Kays模型[10]、文中提出的式(1)模型,計(jì)算高超聲速平板速度型,并與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比結(jié)果。
圖1 高超聲速平板速度型(x=1.80 m,Ma=8.18)Fig.1 Computed velocity profile for a Mach 8.18 flat plate at x=1.80 m
由圖1可見,文中的修正模型與湍流普朗特?cái)?shù)為常數(shù)的模型計(jì)算結(jié)果相同。原因?yàn)槲闹行拚P蛯?duì)高超聲速平板湍流平衡繞流,式(1)等號(hào)右端第二項(xiàng)設(shè)計(jì)的函數(shù)項(xiàng)不起作用,即對(duì)高超聲速平衡繞流流動(dòng),文中的修正模型,退化為原模型的普朗特?cái)?shù)。而對(duì)于Kays模型,由于該模型是基于不可壓縮流動(dòng)提出的,對(duì)高超聲速可壓縮平板流動(dòng),計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相比有一定差距,該模型還需要進(jìn)一步研究改進(jìn)。
圖2給出了原SSTk-ω湍流模型、在原SSTk-ω湍流模型基礎(chǔ)上分別進(jìn)行湍流普朗特修正的Kays模型[10]、文中提出的式(1)模型,計(jì)算的高超聲速平板壁面熱流并與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比曲線。由圖2可見,文中的修正模型與湍流普朗特?cái)?shù)為常數(shù)的模型計(jì)算結(jié)果相同,其原因與對(duì)圖1的分析相同。Kays模型計(jì)算的平板壁面熱流雖然也與實(shí)驗(yàn)結(jié)果符合較好,但比其他模型稍高。原因?yàn)镵ays模型計(jì)算給出的Kármár常數(shù)κ較大,由文獻(xiàn)[14]的分析,則摩擦系數(shù)cf增大。再由摩擦系數(shù)與壁面熱流的雷諾比擬理論,則Kays模型計(jì)算的壁面熱流比其他兩個(gè)模型稍大。
圖2 高超聲速平板壁面熱流Fig.2 Wall heat flux over the flat plate
下面選擇文獻(xiàn)[18]的實(shí)驗(yàn)條件為計(jì)算條件高超聲速15°壓縮拐角繞流進(jìn)行算例驗(yàn)證。算例的來(lái)流馬赫數(shù)及來(lái)流溫度分別為9.22 K、64.5 K,基于來(lái)流參數(shù)的每米雷諾數(shù)為4.7×107m-1,壁面溫度為295 K。
圖3(a)與圖3(b)分別給出了原SSTk-ω湍流模型、在原SSTk-ω湍流模型基礎(chǔ)上分別進(jìn)行湍流普朗特修正的Kays模型[10]、文中提出的式(1)模型,計(jì)算的高超聲速15°壓縮拐角壁面壓強(qiáng)、壁面熱流并與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比曲線。由圖3可見,對(duì)高超聲速15°壓縮拐角繞流,由于拐角角度較小,對(duì)湍流流動(dòng),由激波/湍流邊界層干擾導(dǎo)致逆壓梯度并未使流動(dòng)分離,流動(dòng)的非平衡性較弱,對(duì)式(1)中湍流普朗特?cái)?shù)的影響較小,因此,修正模型與原模型及Kays模型計(jì)算給出的壁面壓強(qiáng)和壁面熱流均符合較好。各模型給出的壁面壓強(qiáng)均與實(shí)驗(yàn)數(shù)據(jù)有一定差別,可能是實(shí)驗(yàn)有一定誤差,需要進(jìn)一步研究。修正模型計(jì)算的表面熱流與實(shí)驗(yàn)結(jié)果更接近,最大熱流相對(duì)誤差小于5%。
圖3 Ma=9.22,15°壓縮拐角壁面壓強(qiáng)、熱流Fig.3 Skin pressure and the skin heat flux distributions for free stream Mach 9.22 flow into a 15° compression corner
下面選擇30°壓縮拐角進(jìn)行算例驗(yàn)證,其余所有條件均與算例3.2節(jié)相同。
圖4(a)與圖4(b)分別給出了原SSTk-ω湍流模型、在原SSTk-ω湍流模型基礎(chǔ)上分別進(jìn)行湍流普朗特修正的Kays模型[10]、文中提出的式(1)模型,計(jì)算的高超聲速30°壓縮拐角壁面壓強(qiáng)、壁面熱流,及與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比曲線。由圖4(a)可見,由于原SSTk-ω湍流模型與Kays模型并未準(zhǔn)確計(jì)及高超聲速激波/湍流邊界層流動(dòng)的湍流普朗特?cái)?shù),導(dǎo)致數(shù)值模擬結(jié)果出現(xiàn)了流動(dòng)分離及再附,以及由于流動(dòng)分離導(dǎo)致了原SSTk-ω湍流模型與Kays模型,在拐角附近區(qū)域出現(xiàn)分離泡(見圖5(a)),在再附區(qū)出現(xiàn)了壓強(qiáng)極值,及壓強(qiáng)極值后的平臺(tái)效應(yīng),從而使計(jì)算的壁面壓強(qiáng)與實(shí)驗(yàn)數(shù)據(jù)相差較大。由于原SSTk-ω湍流模型與Kays模型模擬給出了流動(dòng)分離及再附,但實(shí)驗(yàn)結(jié)果表面,流動(dòng)基本上并未分離,因此導(dǎo)致圖4(b)中原SSTk-ω湍流模型與Kays模型計(jì)算的拐角附近區(qū)域的熱流較大,并且再附區(qū)的熱流遠(yuǎn)高于實(shí)驗(yàn)數(shù)據(jù)。文中提出的修正模型,即式(1),由于正確預(yù)測(cè)了流動(dòng)基本無(wú)分離(見圖5(b))則給出了與實(shí)驗(yàn)符合良好的壁面熱流結(jié)果,最大熱流相對(duì)誤差小于5%。
圖4 Ma=9.22,30°壓縮拐角壁面壓強(qiáng)、熱流Fig.4 Skin pressure and the skin heat flux distributions for free stream Mach 9.22 flow into a 30° compression corner
圖5 Ma=9.22,30°壓縮拐角流線分布Fig.5 Streamline distributions for 30° compression corner at free stream Mach 9.22
下面選擇34°壓縮拐角進(jìn)行算例驗(yàn)證,其余所有條件均與算例3.2節(jié)、3.3節(jié)相同。
圖6(a)與圖6(b)分別給出了原SSTk-ω湍流模型、在原SSTk-ω湍流模型基礎(chǔ)上分別進(jìn)行湍流普朗特修正的Kays模型[10]、文中提出的式(1)模型,計(jì)算的高超聲速34°壓縮拐角壁面壓強(qiáng)、壁面熱流,及與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比曲線。由圖6(a)可見,雖然原SSTk-ω湍流模型與Kays模型均在壓縮拐角附近計(jì)算出了流動(dòng)分離及再附,但由于原SSTk-ω湍流模型與Kays模型并未準(zhǔn)確計(jì)及高超聲速激波/湍流邊界層時(shí)的湍流普朗特?cái)?shù),與實(shí)驗(yàn)結(jié)果相比,預(yù)測(cè)的分離區(qū)偏大(如圖7(a)所示)、峰值壓強(qiáng)及再附區(qū)后移。文中的修正模型則較好地預(yù)測(cè)了流動(dòng)分離及再附(如圖7(b)所示),因此給出與實(shí)驗(yàn)數(shù)據(jù)符合較好結(jié)果。由于原SSTk-ω湍流模型與Kays模型預(yù)測(cè)的分離及再附區(qū)較大,由圖6(b)可見,所計(jì)算的峰值熱流位置及大小均與實(shí)驗(yàn)結(jié)果有一定差別,而文中提出的修正模型的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合較好,最大熱流相對(duì)誤差小于6%。
圖6 Ma=9.22,34°壓縮拐角壁面壓強(qiáng)、熱流Fig.6 Skin pressure and the skin heat flux distributions for Mach 9.22 flow into a 34° compression corner
圖7 Ma=9.22,34°壓縮拐角流線分布Fig.7 Streamline distributions for 34° compression corner at free stream Mach 9.22
應(yīng)用數(shù)值模擬及理論分析方法,針對(duì)高超聲速壓縮拐角激波/湍流邊界層干擾導(dǎo)致的復(fù)雜流動(dòng),提出了一種湍流熱通量建模中的湍流普朗特?cái)?shù)并非常數(shù),而應(yīng)為非平衡流動(dòng)特征參數(shù)的函數(shù)的湍流普朗特?cái)?shù)模型,并與實(shí)驗(yàn)結(jié)果、原湍流普朗特?cái)?shù)模型、Kays湍流普朗特?cái)?shù)模型的計(jì)算結(jié)果進(jìn)行對(duì)比,主要結(jié)論如下:
1)對(duì)高超聲速平板的平衡湍流流動(dòng)建立的模型不改變湍流普朗特?cái)?shù),修正模型與原模型計(jì)算精度相同。
2)對(duì)高超聲速激波/湍流邊界層干擾導(dǎo)致的復(fù)雜非平衡流動(dòng),湍流普朗特?cái)?shù)并非常數(shù),應(yīng)進(jìn)行修正。
3)應(yīng)用所提出的湍流普朗特?cái)?shù)為非平衡流動(dòng)特征參數(shù)的函數(shù)的湍流熱通量模型,數(shù)值模擬與實(shí)驗(yàn)結(jié)果的對(duì)比表明,該模型計(jì)算氣動(dòng)力、氣動(dòng)熱更精確。
雖然文中所提出的模型,在計(jì)算高超聲速壓縮拐角繞流中取得了較好的結(jié)果,但是由于高超聲速激波/湍流邊界層干擾的復(fù)雜性,進(jìn)一步應(yīng)研究所提出的修正湍流普朗特?cái)?shù)模型對(duì)復(fù)雜流動(dòng)的模擬能力。