劉朋欣 孫 東 李 辰 郭啟龍 袁先旭
(中國空氣動力研究與發(fā)展中心空氣動力學國家重點實驗室,四川綿陽,621000)
(中國空氣動力研究與發(fā)展中心計算空氣動力研究所,四川綿陽,621000)
高超聲速飛行器在以極高馬赫數(shù)飛行時,飛行器周圍會產(chǎn)生較強的激波并壓縮來流空氣,使其溫度急劇升高.溫度的升高會激發(fā)復雜的熱力學過程,比如振動能激發(fā)和氣體分子離解,使飛行器表面的氣體處于非平衡狀態(tài)[1].當飛行器在低空稠密大氣高速飛行時,飛行器表面流動會發(fā)生轉(zhuǎn)捩并形成湍流.湍流會使得壁面摩阻升高數(shù)倍[2],損害飛行器飛行性能.湍流與非平衡效應耦合在一起,使得表面摩阻的預測更加困難,高溫非平衡效應對壁面摩阻產(chǎn)生機制的影響也尚不清楚.
目前同時考慮高溫化學非平衡效應和湍流流動的研究尚處于起步階段,相關(guān)報道還十分有限,且均是采用直接數(shù)值模擬方法(direct numerical simulation,DNS)進行分析研究.Duan 和 Martin[3-4]使用5組分空氣反應模型進行了DNS,結(jié)果表明高焓化學非平衡條件下基于Morkovin 假設的大部分相關(guān)尺度關(guān)系仍然成立,并評估了湍流對化學反應速率的影響.Kim[5]使用單一組分的O2或N2,并考慮熱化學非平衡過程,發(fā)現(xiàn)半當?shù)爻叨瓤梢暂^好地消除熱化學非平衡效應的影響.劉朋欣等[6-7]發(fā)現(xiàn)在邊界層的對數(shù)區(qū),各流動參數(shù)的拓展自相似理論的相對標度指數(shù)基本符合標度規(guī)律.吳正園等[8]的研究表明湍流邊界層中的高溫氣體效應對壁面平均壓強和脈動壓強有顯著的增加效果.Renzo 和Urzay[9-10]研究了含化學非平衡的邊界層轉(zhuǎn)捩過程,并基于模擬結(jié)果分析了壁面摩阻、熱流的流向演化.Passiatore 等[11]對比了兩種不同氣體模型(凍結(jié)流動和化學非平衡流動)下的湍流邊界層,發(fā)現(xiàn)吸熱反應會使得溫度脈動和密度脈動低于凍結(jié)流模型,但化學非平衡效應對速度的一階和二階統(tǒng)計量影響較小,并簡單地采用RD (Renard &Deck)方法對摩阻進行了分解研究.Volpiani[12]開發(fā)了一套模擬高超聲速激波/化學非平衡湍流邊界層相互的直接數(shù)值模擬程序,發(fā)現(xiàn)化學非平衡效應使得激波干擾區(qū)的流動特性發(fā)生了顯著變化.
上述研究結(jié)果表明,高溫非平衡效應會顯著影響流場熱力學特性分布.溫度分布的改變又會直接影響?zhàn)ば韵禂?shù)的大小,繼而對摩阻分布產(chǎn)生影響.因此有必要研究高溫非平衡效應下的壁面摩阻產(chǎn)生機制.目前高溫非平衡湍流邊界層的研究大都關(guān)注于湍流邊界層宏觀統(tǒng)計特性規(guī)律.而對工程中較為關(guān)心的摩阻特性缺乏深入分析.
摩阻分解技術(shù)將壁面摩阻的產(chǎn)生與流動機制結(jié)合起來,是分析摩阻產(chǎn)生機理的一種有效手段.目前常見的有FIK 分解[13]和RD 分解[14].Fukagata 等[13]通過對平均動量方程的連續(xù)三次積分,建立了摩阻系數(shù)與雷諾應力之間的關(guān)系,將摩阻分解為層流部分的貢獻和湍流部分的貢獻.后續(xù)的改進工作提升了公式的易用性[15-16],并將該公式推廣至任意復雜外形[17].Renard 和Deck 指出[14]FIK 公式?jīng)]有明確的物理意義,并提出了另外一種摩阻分解方法,即RD 分解該方法在假設遠離壁面附近流體保持靜止的絕對參考坐標系下,通過對流向動量方程做變換和積分,將摩阻的產(chǎn)生分為分子黏性耗散作用和湍動能生成作用,具有更明確的物理意義.此方法在旋渦結(jié)構(gòu)輸運過程[18]、大尺度湍流運動[19]和破碎[20]對摩阻的影響研究得到了應用.Li 等[21-22]將其推廣至可壓縮湍流中,并研究了雷諾數(shù)效應對摩阻產(chǎn)生的影響.
本文利用RD 摩阻分解技術(shù),分析了高溫化學非平衡湍流邊界層摩阻生成機制,給出了對摩阻生成起主要作用的流動過程;通過設置低焓完全氣體模型對比算例,從摩阻空間分布特點、熱力學和流動特性方面分析了高焓效應對主要摩阻產(chǎn)生項的影響.
低焓算例工質(zhì)為空氣,采用完全氣體模型,控制方程為可壓縮Naiver-Stokes (N-S);高焓算例由于需要考慮空氣的化學非平衡效應,采用含化學反應的N-S 方程作為控制方程.后者的形式如下
其中 ρ和p是密度和壓力;uj為分速度;fk是第k種組分的質(zhì)量分數(shù);E為單位體積的總內(nèi)能;hk是比焓,通過溫度的多項式擬合得到;Sk是第k種組分的生成速率; τij為黏性應力;K和Dk分別為熱傳導系數(shù)和組分質(zhì)量擴散系數(shù).在完全氣體計算中,通常取Pr=0.72,比熱比 γ=1.4,黏性系數(shù)采用Sutherland 公式計算.而對于多組分氣體,則首先通過Lennard-Jones 分子黏性系數(shù)計算公式求得各組分的黏性系數(shù),并通過Wilke 公式計算混合物黏性系數(shù); γ 也不再為常數(shù);通過Pr數(shù)和S c數(shù)與黏性系數(shù)的關(guān)系式來計算K和Dk.多組分輸運系數(shù)(黏性系數(shù)、熱傳導系數(shù)和組分擴散系數(shù))具體計算公式可參考文獻[23].
計算采用課題組發(fā)展的Inhouse 程序.該計算程序已經(jīng)成功開展大量的可壓縮湍流[24-25]、激波/邊界層干擾[26-27]、旋轉(zhuǎn)爆震精細流場結(jié)構(gòu)[28]、高溫非平衡湍流邊界層[6-7]等復雜流動的數(shù)值模擬,計算精度和魯棒性均得到了充分驗證.在本文計算中,無黏通量采用了7 階有限差分格式WENO-Z[29],黏性項采用4 階中心差分格式,時間推進采用具有TVD 性質(zhì)的3 階Runge-Kutta 方法.高焓算例計算不考慮溫度非平衡效應,且只考慮空氣組分的離解過程,化學反應源項求解采用Gupta 等[30]提出的5 組分(N2,O2,N,O,NO) 6 基元反應模型.
高焓算例的計算狀態(tài)示意圖見圖1.為了激發(fā)空氣的化學非平衡態(tài),選擇30 km 高空以馬赫數(shù)20 飛行的楔形體,其頭部產(chǎn)生一道斜激波.斜激波后的溫度較高且已經(jīng)足夠激發(fā)空氣發(fā)生化學反應,選擇此狀態(tài)作為湍流邊界層外緣流動狀態(tài).高焓算例來流組分質(zhì)量分數(shù)設置為f(O2)=0.267 5,f(N2)=0.732 5,并使用空氣化學反應模型模擬化學非平衡過程,記為TH 算例.作為對比,設置同等邊界層參數(shù)下的低焓算例,以評估高溫非平衡效應對湍流邊界層統(tǒng)計特性和摩阻的影響,并記為TL 算例.兩算例的來流狀態(tài)及壁溫設置見表1.
表1 來流狀態(tài)及壁溫設置Table 1 Inflow condition and wall temperature
圖1 計算模型狀態(tài)示意圖Fig.1 Schematic of computational model
計算初場采用相同來流狀態(tài)下RANS[31]計算結(jié)果的剖面生成.入口的脈動通過數(shù)字濾波合成湍流方法[32]產(chǎn)生;出口采用特征邊界條件;上邊界固定為初始值;壁面采用黏性等溫條件;展向采用周期性邊界條件.
計算域設置為20 δ × 4 δ × 2 δ (流向 × 法向 × 展向,δ 為邊界層名義厚度),相應的網(wǎng)格為901 × 161 ×301.流向和展向網(wǎng)格均勻,法向網(wǎng)格在壁面附近加密.文獻[6]進行了同等工況下的網(wǎng)格無關(guān)性驗證,本文所采用的網(wǎng)格可以滿足直接數(shù)值模擬的要求.湍流充分發(fā)展段的網(wǎng)格分辨率及邊界層厚度和雷諾數(shù)見表2.此時的邊界層厚度 δ 約為5 mm.保持兩個算例的Reδ2一致,且可以發(fā)現(xiàn)Reθ和Reτ也相差不大.分析所采用流向位置如圖2(b)中黑色實線所示.且如無特殊說明,下文中分析默認基于此流向截面處的流場.在下文的表述中約定,代表變量q的Reynolds 平均,表示Favre 平均.兩種平均方式所對應的脈動量分別為:.
圖2 瞬時流場結(jié)構(gòu)圖Fig.2 Instantaneous flow structures
表2 邊界層厚度、雷諾數(shù)和網(wǎng)格尺度Table 2 Thickness,Reynold number and grid resolution
圖2 給出了瞬時流場旋渦結(jié)構(gòu)和密度梯度分布圖.從圖中可以看出,在入口添加的擾動,經(jīng)過一段距離的發(fā)展,形成了充分發(fā)展的湍流,流動結(jié)構(gòu)變得紊亂無序.從整體上看,高焓算例和低焓算例的流動結(jié)構(gòu)沒有明顯的變化,說明化學非平衡過程對流動運動過程影響較小,并不會顯著改變流場結(jié)構(gòu)特性.
圖3 和圖4 分別給出了平均流向速度和雷諾應力分布.Van Driest 變換的定義為:.從圖3 中可以看出,此時速度邊界層仍然存在明顯的線性區(qū)、緩沖區(qū)和對數(shù)區(qū),只不過截距有所增大,與文獻中高焓算例[3]結(jié)果符合較好.除了在緩沖區(qū),高焓算例的平均流向速度略高于低焓算例之外,焓值并不會明顯改變湍流邊界層速度剖面的整體分布特性.從圖4 中可以看出,經(jīng)過密度加權(quán)的雷諾應力分布具有相似的分布趨勢,與文獻完全氣體結(jié)果[33]及高焓算例結(jié)果[3]符合較好.高焓效應對雷諾應力分布的影響不大,并沒有對湍流速度脈動產(chǎn)生明顯影響,這與圖2 中的瞬時流場結(jié)構(gòu)對比結(jié)果是一致的.
圖3 Van Driest 變換后的速度分布Fig.3 Profile of Van Driest transformed velocity
圖4 雷諾應力分布Fig.4 Reynolds stress distribution
本節(jié)首先對比高低焓值下摩阻的空間分布,分析高溫非平衡效應對摩阻分布的整體影響.圖5 中給出了不同焓值工況下的局部瞬時摩阻系數(shù)在壁面上的分布,從圖中可以看出瞬時摩阻的分布是相似的,都呈現(xiàn)出局部高摩阻斑點,并且存在一些摩阻較弱的條帶結(jié)構(gòu).另外,高焓情況下的條帶結(jié)構(gòu)尺度有減小趨勢.
圖5 瞬時壁面摩阻系數(shù)分布Fig.5 Instantaneous distribution of skin friction coefficients
圖6 給出了摩阻脈動的概率分布函數(shù)(probability distribution function,PDF),使用其均方根進行歸一化,并與高斯分布函數(shù)進行對比.結(jié)果表明,高焓非平衡流動和低焓流動的PDF 曲線均呈現(xiàn)出傾斜趨勢;兩者分布基本一致且均偏離高斯分布的中心位置,峰值位置約在-0.5 處;高焓算例的峰值略低于低焓算例,但兩者的峰值概率密度略高于高斯分布;計算結(jié)果與文獻[34]完全氣體模型(Reθ=2300)較為一致.這說明摩阻脈動的PDF 受化學非平衡效應的影響有限.PDF 曲線向負方向一側(cè)傾斜,說明當?shù)氐退贄l帶發(fā)生的概率較大;>0 對應于高速條帶流動事件,且摩阻脈動的幅值達到4 左右,表明高速條帶產(chǎn)生的脈動更強.
圖6 摩阻脈動PDF 分布Fig.6 PDF profile of skin friction fluctuations
圖5 的結(jié)果表明,摩阻脈動的高低值條帶在尺寸上存在差異.進一步可以通過摩阻脈動的空間兩點相關(guān)系數(shù)來定量分析高溫非平衡效應對摩阻脈動條帶尺度的影響.壁面上相關(guān)系數(shù)定義為
其中代表摩阻脈動,Δx,Δz分別代表流場和展向間距,(x0,z0) 在所選流向站位的展向中間位置.式(3)是指摩阻脈動的空間兩點相關(guān)系數(shù)值等于0.2 的流向尺度和展向尺度.
圖7 給出了相關(guān)系數(shù)的分布云圖.不同焓值下的相關(guān)系數(shù)均呈現(xiàn)出橢圓形分布,且流向尺度遠大于展向尺度.高焓流動的分布范圍在流向和展向兩個方向上均小于低焓流動.這與圖5 中瞬時摩阻脈動的條帶尺寸一致.沿兩個方向進行長度尺度積分[34],可以得到定量的尺寸,具體見式(3).低焓流動的流向尺度為1.41δ,展向尺度為0.12δ;高焓流動的流向尺度為1.20δ,展向尺度為0.11δ.可見高焓算例的流向尺度減小了14.9%,展向尺度減小了8.3%.
圖7 壁面摩阻脈動兩點相關(guān)系數(shù)空間分布
進一步,采用壁面黏性尺度δv進行對比.低焓流動的流向尺度為1130.6δv,展向尺度為96.0δv;高焓流動的流向尺度為993.5δv,展向尺度為95.5δv.可見高焓算例的流向尺度減小了12.1%,展向尺寸基本一致.
高焓條件下摩阻條帶結(jié)構(gòu)尺度減小的原因可能來自于高溫非平衡效應.熱力學量的脈動引起局部流向速度梯度增大,如圖3 中所示在緩沖區(qū)較高,條帶結(jié)構(gòu)的能量更高.相較低焓流動,輸運相同能量,高焓邊界層條帶結(jié)構(gòu)更小.具體原因還需要進一步驗證.
研究摩阻的產(chǎn)生機理以及與之相關(guān)的近壁區(qū)動力學過程,對于基礎研究和工程實際都有重要意義,可以為減阻方法的設計提供指導.對于可壓縮壁湍流,假設展向均勻,可將壁面摩阻寫成如下分解積分關(guān)系式[22]
其中 ρ∞和U∞分別為來流的密度和速度,本算例中分別等于 ρe和Ue.Cf,V代表分子黏性耗散項;Cf,T代表湍動能生成項;Cf,G1表示平均流向動能的變化項;Cf,G2表示流向黏性正應力不均勻性項,Cf,G3表示流向雷諾正應力不均勻性項,兩者合起來稱為流向不均勻項.其中.
在對湍流邊界層摩阻空間分布分析的基礎上,本節(jié)繼續(xù)采用摩阻分解式(4)來分析邊界層中摩阻生成機制.表3 比較了采用摩阻分解技術(shù)得到的摩阻平均值(Cf,decomp)、直接采用時間平均得到的摩阻值(Cf,avg)以及兩者之間的相對誤差.相對誤差的定義為error=(Cf,decomp-Cf,avg)/Cf,avg.結(jié)果顯示,本文實現(xiàn)的摩阻分解具有較高的精度,相對誤差均在1.5%以下.
表3 摩阻分解公式相對誤差Table 3 Relative error of decomposition formula
表4 給出了式(4)中不同流動和輸運過程對摩阻產(chǎn)生的貢獻比例.可以看出,分子黏性耗散項(Cf,V)和湍動能生成項(Cf,T) 所占比例最大,基本上在40%以上;其次是平均流向動能的變化項(Cf,G1),比例在10%;流向黏性正應力不均勻性項(Cf,G2)和流向雷諾正應力不均勻性項(Cf,G3)占比很小,基本上可以忽略.對比不同焓值算例可以發(fā)現(xiàn),高焓值情況下,Cf,V和Cf,T的占比略有下降,減少的這部分貢獻基本上由Cf,G1補上.摩阻分解分析表明,湍流邊界層中對摩阻生成起主要作用的流動過程是Cf,V和Cf,T.下面將進一步分析兩者產(chǎn)生的原因以及高低焓值邊界層摩阻產(chǎn)生差異的原因.
表4 摩阻分解公式中各項貢獻Table 4 Contributions of different terms to skin friction
分子黏性耗散過程是壁面摩阻產(chǎn)生的主要流動過程之一.圖8 給出了分子黏性耗散項預乘積分函數(shù)沿法向的分布.兩算例的分布趨勢較為一致,約在y≈0.015δ(y+≈10)處達到峰值.且分子黏性耗散項的主要作用區(qū)域在y≤0.1δ(y+≤30)的近壁區(qū),而在外區(qū)基本沒有貢獻.這與文獻中基于槽道、平板湍流[21-22]的分析結(jié)果一致.高焓會使得y<0.015δ 以內(nèi)的近壁區(qū)分子黏性耗散略大于低焓流動;而在y>0.015δ的區(qū)域中略小于低焓流動.
圖8 黏性耗散項預乘積分函數(shù)Fig.8 Pre-multiplied integrand of viscous dissipation terms
進一步可以將分子黏性耗散項分解為平均部分Cf,V,m和脈動部分Cf,V,f[21]
表4 中給出了兩部分的占比.可以發(fā)現(xiàn),Cf,V,m占據(jù)了分子黏性耗散項的絕大部分,占比高達95%以上.這就說明焓值對分子黏性耗散項的影響主要通過Cf,V,m項體現(xiàn).
將Cf,V,m的積分函數(shù)分為兩部分:和,分別體現(xiàn)了流動特性和熱力學特性的作用.圖9給出了兩者沿壁面法向的分布,并給出了溫度的分布,均采用無窮遠處參數(shù)作為無量綱參考值.從流動特性可以看出,高焓使得近壁區(qū)(y<0.015δ)的?u/?y增大(?v/?x較小,影響可以忽略);而對外區(qū)的流動特性影響很小.從熱力學特性可以看出,高焓對整個邊界層內(nèi)的黏性系數(shù)都有較大影響,基本上均低于低焓算例的黏性系數(shù).黏性系數(shù)主要受到溫度的影響,兩者的分布具有明顯的一致性.當高焓時,邊界層內(nèi)發(fā)生空氣化學非平衡過程,離解反應會吸收一部分熱量,使得高焓邊界層內(nèi)的黏性系數(shù)低于低焓邊界層.綜合流動特性和熱力學特性的分布,可以發(fā)現(xiàn),在y<0.015δ 時,流動特性起到主導作用,導致高焓流動黏性耗散項貢獻略大于低焓流動;而在y>0.015δ 時,熱力學特性起到主要作用,導致低焓流動黏性耗散項貢獻略大于高焓流動.
圖9 Cf,V,m積分函數(shù)分解及溫度分布Fig.9 Integrand parts of Cf,V,mand the temperature distribution
湍動能生成項是壁面摩阻產(chǎn)生的另一個重要部分.圖10 給出了湍動能生成項預乘積分函數(shù)沿法向的分布.兩算例的分布趨勢較為一致;但低焓流動時,峰值位置更加遠離壁面.結(jié)合圖4 中雷諾應力分布基本一致的情況來看,峰值位置的外移應該來自于的分布不同;結(jié)合圖9 來看,高焓流動在預乘積分函數(shù)附近的要大于低焓流動.兩算例在邊界層外區(qū)均出現(xiàn)了二次峰值,這與較高雷諾數(shù)下的槽道流動[21]一致.二次峰值的產(chǎn)生是由于較高雷諾數(shù)下,邊界層外區(qū)存在大尺度湍流運動導致.大尺度湍流運動會引起外區(qū)產(chǎn)生大量湍動能,從而通過摩阻的湍動能生成項影響摩阻;且外區(qū)大尺度湍流流動對內(nèi)區(qū)小尺度湍流運動具有調(diào)制作用,促發(fā)湍動能在內(nèi)外區(qū)的重新分配[35].
圖10 湍動能生成項預乘積分函數(shù)Fig.10 Pre-multiplied integrand of TKE production terms
為了進一步分析湍動能生成項中起主導作用的流動事件,本節(jié)采用象限分析技術(shù)[36-37]將邊界層中的主導流動事件與摩阻分解相結(jié)合以分析不同流動事件對摩阻中湍動能生成項的影響.
在邊界層中,可以按照脈動的正負將流向速度脈動和法向速度脈動分解到4 個象限中.第1 象限:u′′>0,v′′>0表示外向運動(outward motion);第2 象限:u′′<0,v′′>0表示上拋運動;第3 象限:u′′<0,v′′<0表示內(nèi)向運動(inward motion);第4 象限:u′′>0,v′′<0表示下掃運動.圖11 給出了Q1~Q4 事件在邊界層內(nèi)發(fā)生的概率.不同焓值下不同事件的概率分布曲線是相似的.Q2 和Q4 發(fā)生的概率顯著大于Q1 和Q3.
圖11 Q1~Q4 事件發(fā)生概率Fig.11 Probability distribution of Q1~Q4 events
圖12 中給出了不同事件對應的雷諾應力與總雷諾切應力的比值.可以看出,不同焓值的比值曲線很好的重合在一起,說明該比值對焓值變化不敏感;Q2 和Q4 事件是雷諾切應力的主要貢獻項,并且在y≈0.03δ以內(nèi),Q4 事件的比值要大于Q2 事件;在y≈0.03δ以外,Q2 事件的比值大于Q4 事件,且這一趨勢會在整個邊界層范圍內(nèi)得到保持.
表5 定量給出了不同事件的湍動能生成項貢獻的摩阻與整體摩阻的比值,表中以Cf,T,Qi表示Qi事件對湍動能生成項的貢獻.可以看出,Q2 和Q4 事件起到絕對主導作用,且Q2 貢獻要略大于Q4 事件,這與圖12 中所指出的Q2 對雷諾應力的貢獻大于Q4 一致.且Q1 和Q3 對壁面摩阻的貢獻為負,兩者的值基本相同,且都遠小于Q2 和Q4.
表5 Q1~Q4 湍流生成項對壁面摩阻的貢獻Table 5 Contributions of TKE production terms produced by Q1~Q4 to the time-averaged value
圖12 Q1~Q4 事件對湍動能生成項的貢獻Fig.12 Contributions of Q1~Q4 events to TKE production terms
圖13 給出了不同事件的被積函數(shù)沿法向的分布,以考察主要起作用的邊界層區(qū)域范圍.可以看出,不同焓值下,被積函數(shù)的曲線分布趨勢都是相似的,高焓效應僅會稍微改變幅值;Q4 事件的峰值位置要比Q2 事件更為靠內(nèi);Q1 和Q3 在整個邊界層范圍內(nèi)的值都是負值.
圖13 Q1~Q4 湍流生成項預乘積分函數(shù)Fig.13 Pre-multiplied integrand of TKE production terms produced by Q1~Q4
本文通過采用直接數(shù)值模擬方法,研究了高焓湍流邊界層壁面摩阻產(chǎn)生機制,并通過設置同等邊界層參數(shù)下的低焓算例,對比分析了高溫化學非平衡效應對摩阻產(chǎn)生的影響,可以得到以下結(jié)論.
(1)不同焓值下的摩阻脈動條帶尺寸存在差異.高焓流動的條帶尺寸在流向和展向兩個方向上均小于低焓流動.
(2)湍流邊界層中對摩阻生成起主要作用的流動過程是分子黏性耗散項Cf,V和湍動能生成項Cf,T.
(3) 分子黏性耗散項Cf,V的主要作用區(qū)域在y≤0.1δ的近壁區(qū).焓值的影響主要通過Cf,V的平均部分Cf,V,m項體現(xiàn).且在y<0.015δ 時,流動特性(速度梯度)起到主導作用,導致高焓流動Cf,V,m略大于低焓流動;而在y>0.015δ 時,熱力學特性(黏性系數(shù))起到主要作用,導致低焓流動Cf,V,m略大于高焓流動.
(4)邊界層中的上拋和下掃運動是影響摩阻分解湍動能生成項Cf,T的主導事件,外向運動和內(nèi)向運動對Cf,T的影響很小,且起到負貢獻.