高焓湍流邊界層壁面摩阻產生機制分析1)

劉朋欣 孫 東 李 辰 郭啟龍 袁先旭

(中國空氣動力研究與發展中心空氣動力學國家重點實驗室,四川綿陽,621000)

(中國空氣動力研究與發展中心計算空氣動力研究所,四川綿陽,621000)

引言

高超聲速飛行器在以極高馬赫數飛行時,飛行器周圍會產生較強的激波并壓縮來流空氣,使其溫度急劇升高.溫度的升高會激發復雜的熱力學過程,比如振動能激發和氣體分子離解,使飛行器表面的氣體處于非平衡狀態[1].當飛行器在低空稠密大氣高速飛行時,飛行器表面流動會發生轉捩并形成湍流.湍流會使得壁面摩阻升高數倍[2],損害飛行器飛行性能.湍流與非平衡效應耦合在一起,使得表面摩阻的預測更加困難,高溫非平衡效應對壁面摩阻產生機制的影響也尚不清楚.

目前同時考慮高溫化學非平衡效應和湍流流動的研究尚處于起步階段,相關報道還十分有限,且均是采用直接數值模擬方法(direct numerical simulation,DNS)進行分析研究.Duan 和 Martin[3-4]使用5組分空氣反應模型進行了DNS,結果表明高焓化學非平衡條件下基于Morkovin 假設的大部分相關尺度關系仍然成立,并評估了湍流對化學反應速率的影響.Kim[5]使用單一組分的O2或N2,并考慮熱化學非平衡過程,發現半當地尺度可以較好地消除熱化學非平衡效應的影響.劉朋欣等[6-7]發現在邊界層的對數區,各流動參數的拓展自相似理論的相對標度指數基本符合標度規律.吳正園等[8]的研究表明湍流邊界層中的高溫氣體效應對壁面平均壓強和脈動壓強有顯著的增加效果.Renzo 和Urzay[9-10]研究了含化學非平衡的邊界層轉捩過程,并基于模擬結果分析了壁面摩阻、熱流的流向演化.Passiatore 等[11]對比了兩種不同氣體模型(凍結流動和化學非平衡流動)下的湍流邊界層,發現吸熱反應會使得溫度脈動和密度脈動低于凍結流模型,但化學非平衡效應對速度的一階和二階統計量影響較小,并簡單地采用RD (Renard &Deck)方法對摩阻進行了分解研究.Volpiani[12]開發了一套模擬高超聲速激波/化學非平衡湍流邊界層相互的直接數值模擬程序,發現化學非平衡效應使得激波干擾區的流動特性發生了顯著變化.

上述研究結果表明,高溫非平衡效應會顯著影響流場熱力學特性分布.溫度分布的改變又會直接影響黏性系數的大小,繼而對摩阻分布產生影響.因此有必要研究高溫非平衡效應下的壁面摩阻產生機制.目前高溫非平衡湍流邊界層的研究大都關注于湍流邊界層宏觀統計特性規律.而對工程中較為關心的摩阻特性缺乏深入分析.

摩阻分解技術將壁面摩阻的產生與流動機制結合起來,是分析摩阻產生機理的一種有效手段.目前常見的有FIK 分解[13]和RD 分解[14].Fukagata 等[13]通過對平均動量方程的連續三次積分,建立了摩阻系數與雷諾應力之間的關系,將摩阻分解為層流部分的貢獻和湍流部分的貢獻.后續的改進工作提升了公式的易用性[15-16],并將該公式推廣至任意復雜外形[17].Renard 和Deck 指出[14]FIK 公式沒有明確的物理意義,并提出了另外一種摩阻分解方法,即RD 分解該方法在假設遠離壁面附近流體保持靜止的絕對參考坐標系下,通過對流向動量方程做變換和積分,將摩阻的產生分為分子黏性耗散作用和湍動能生成作用,具有更明確的物理意義.此方法在旋渦結構輸運過程[18]、大尺度湍流運動[19]和破碎[20]對摩阻的影響研究得到了應用.Li 等[21-22]將其推廣至可壓縮湍流中,并研究了雷諾數效應對摩阻產生的影響.

本文利用RD 摩阻分解技術,分析了高溫化學非平衡湍流邊界層摩阻生成機制,給出了對摩阻生成起主要作用的流動過程;通過設置低焓完全氣體模型對比算例,從摩阻空間分布特點、熱力學和流動特性方面分析了高焓效應對主要摩阻產生項的影響.

1 數值方法與算例設置

1.1 控制方程

低焓算例工質為空氣,采用完全氣體模型,控制方程為可壓縮Naiver-Stokes (N-S);高焓算例由于需要考慮空氣的化學非平衡效應,采用含化學反應的N-S 方程作為控制方程.后者的形式如下

其中 ρ和p是密度和壓力;uj為分速度;fk是第k種組分的質量分數;E為單位體積的總內能;hk是比焓,通過溫度的多項式擬合得到;Sk是第k種組分的生成速率; τij為黏性應力;K和Dk分別為熱傳導系數和組分質量擴散系數.在完全氣體計算中,通常取Pr=0.72,比熱比 γ=1.4,黏性系數采用Sutherland 公式計算.而對于多組分氣體,則首先通過Lennard-Jones 分子黏性系數計算公式求得各組分的黏性系數,并通過Wilke 公式計算混合物黏性系數; γ 也不再為常數;通過Pr數和S c數與黏性系數的關系式來計算K和Dk.多組分輸運系數(黏性系數、熱傳導系數和組分擴散系數)具體計算公式可參考文獻[23].

1.2 數值方法

計算采用課題組發展的Inhouse 程序.該計算程序已經成功開展大量的可壓縮湍流[24-25]、激波/邊界層干擾[26-27]、旋轉爆震精細流場結構[28]、高溫非平衡湍流邊界層[6-7]等復雜流動的數值模擬,計算精度和魯棒性均得到了充分驗證.在本文計算中,無黏通量采用了7 階有限差分格式WENO-Z[29],黏性項采用4 階中心差分格式,時間推進采用具有TVD 性質的3 階Runge-Kutta 方法.高焓算例計算不考慮溫度非平衡效應,且只考慮空氣組分的離解過程,化學反應源項求解采用Gupta 等[30]提出的5 組分(N2,O2,N,O,NO) 6 基元反應模型.

1.3 算例設置

高焓算例的計算狀態示意圖見圖1.為了激發空氣的化學非平衡態,選擇30 km 高空以馬赫數20 飛行的楔形體,其頭部產生一道斜激波.斜激波后的溫度較高且已經足夠激發空氣發生化學反應,選擇此狀態作為湍流邊界層外緣流動狀態.高焓算例來流組分質量分數設置為f(O2)=0.267 5,f(N2)=0.732 5,并使用空氣化學反應模型模擬化學非平衡過程,記為TH 算例.作為對比,設置同等邊界層參數下的低焓算例,以評估高溫非平衡效應對湍流邊界層統計特性和摩阻的影響,并記為TL 算例.兩算例的來流狀態及壁溫設置見表1.

表1 來流狀態及壁溫設置Table 1 Inflow condition and wall temperature

圖1 計算模型狀態示意圖Fig.1 Schematic of computational model

計算初場采用相同來流狀態下RANS[31]計算結果的剖面生成.入口的脈動通過數字濾波合成湍流方法[32]產生;出口采用特征邊界條件;上邊界固定為初始值;壁面采用黏性等溫條件;展向采用周期性邊界條件.

計算域設置為20 δ × 4 δ × 2 δ (流向 × 法向 × 展向,δ 為邊界層名義厚度),相應的網格為901 × 161 ×301.流向和展向網格均勻,法向網格在壁面附近加密.文獻[6]進行了同等工況下的網格無關性驗證,本文所采用的網格可以滿足直接數值模擬的要求.湍流充分發展段的網格分辨率及邊界層厚度和雷諾數見表2.此時的邊界層厚度 δ 約為5 mm.保持兩個算例的Reδ2一致,且可以發現Reθ和Reτ也相差不大.分析所采用流向位置如圖2(b)中黑色實線所示.且如無特殊說明,下文中分析默認基于此流向截面處的流場.在下文的表述中約定,代表變量q的Reynolds 平均,表示Favre 平均.兩種平均方式所對應的脈動量分別為:.

圖2 瞬時流場結構圖Fig.2 Instantaneous flow structures

表2 邊界層厚度、雷諾數和網格尺度Table 2 Thickness,Reynold number and grid resolution

2 結果分析與討論

2.1 湍流邊界層瞬時與統計結果

圖2 給出了瞬時流場旋渦結構和密度梯度分布圖.從圖中可以看出,在入口添加的擾動,經過一段距離的發展,形成了充分發展的湍流,流動結構變得紊亂無序.從整體上看,高焓算例和低焓算例的流動結構沒有明顯的變化,說明化學非平衡過程對流動運動過程影響較小,并不會顯著改變流場結構特性.

圖3 和圖4 分別給出了平均流向速度和雷諾應力分布.Van Driest 變換的定義為:.從圖3 中可以看出,此時速度邊界層仍然存在明顯的線性區、緩沖區和對數區,只不過截距有所增大,與文獻中高焓算例[3]結果符合較好.除了在緩沖區,高焓算例的平均流向速度略高于低焓算例之外,焓值并不會明顯改變湍流邊界層速度剖面的整體分布特性.從圖4 中可以看出,經過密度加權的雷諾應力分布具有相似的分布趨勢,與文獻完全氣體結果[33]及高焓算例結果[3]符合較好.高焓效應對雷諾應力分布的影響不大,并沒有對湍流速度脈動產生明顯影響,這與圖2 中的瞬時流場結構對比結果是一致的.

圖3 Van Driest 變換后的速度分布Fig.3 Profile of Van Driest transformed velocity

圖4 雷諾應力分布Fig.4 Reynolds stress distribution

2.2 湍流邊界層摩阻空間分布

本節首先對比高低焓值下摩阻的空間分布,分析高溫非平衡效應對摩阻分布的整體影響.圖5 中給出了不同焓值工況下的局部瞬時摩阻系數在壁面上的分布,從圖中可以看出瞬時摩阻的分布是相似的,都呈現出局部高摩阻斑點,并且存在一些摩阻較弱的條帶結構.另外,高焓情況下的條帶結構尺度有減小趨勢.

圖5 瞬時壁面摩阻系數分布Fig.5 Instantaneous distribution of skin friction coefficients

圖6 給出了摩阻脈動的概率分布函數(probability distribution function,PDF),使用其均方根進行歸一化,并與高斯分布函數進行對比.結果表明,高焓非平衡流動和低焓流動的PDF 曲線均呈現出傾斜趨勢;兩者分布基本一致且均偏離高斯分布的中心位置,峰值位置約在-0.5 處;高焓算例的峰值略低于低焓算例,但兩者的峰值概率密度略高于高斯分布;計算結果與文獻[34]完全氣體模型(Reθ=2300)較為一致.這說明摩阻脈動的PDF 受化學非平衡效應的影響有限.PDF 曲線向負方向一側傾斜,說明當地低速條帶發生的概率較大;＞0 對應于高速條帶流動事件,且摩阻脈動的幅值達到4 左右,表明高速條帶產生的脈動更強.

圖6 摩阻脈動PDF 分布Fig.6 PDF profile of skin friction fluctuations

圖5 的結果表明,摩阻脈動的高低值條帶在尺寸上存在差異.進一步可以通過摩阻脈動的空間兩點相關系數來定量分析高溫非平衡效應對摩阻脈動條帶尺度的影響.壁面上相關系數定義為

其中代表摩阻脈動,Δx,Δz分別代表流場和展向間距,(x0,z0) 在所選流向站位的展向中間位置.式(3)是指摩阻脈動的空間兩點相關系數值等于0.2 的流向尺度和展向尺度.

圖7 給出了相關系數的分布云圖.不同焓值下的相關系數均呈現出橢圓形分布,且流向尺度遠大于展向尺度.高焓流動的分布范圍在流向和展向兩個方向上均小于低焓流動.這與圖5 中瞬時摩阻脈動的條帶尺寸一致.沿兩個方向進行長度尺度積分[34],可以得到定量的尺寸,具體見式(3).低焓流動的流向尺度為1.41δ,展向尺度為0.12δ;高焓流動的流向尺度為1.20δ,展向尺度為0.11δ.可見高焓算例的流向尺度減小了14.9%,展向尺度減小了8.3%.

圖7 壁面摩阻脈動兩點相關系數空間分布

進一步,采用壁面黏性尺度δv進行對比.低焓流動的流向尺度為1130.6δv,展向尺度為96.0δv;高焓流動的流向尺度為993.5δv,展向尺度為95.5δv.可見高焓算例的流向尺度減小了12.1%,展向尺寸基本一致.

高焓條件下摩阻條帶結構尺度減小的原因可能來自于高溫非平衡效應.熱力學量的脈動引起局部流向速度梯度增大,如圖3 中所示在緩沖區較高,條帶結構的能量更高.相較低焓流動,輸運相同能量,高焓邊界層條帶結構更小.具體原因還需要進一步驗證.

2.3 湍流邊界層摩阻分解

研究摩阻的產生機理以及與之相關的近壁區動力學過程,對于基礎研究和工程實際都有重要意義,可以為減阻方法的設計提供指導.對于可壓縮壁湍流,假設展向均勻,可將壁面摩阻寫成如下分解積分關系式[22]

其中 ρ∞和U∞分別為來流的密度和速度,本算例中分別等于 ρe和Ue.Cf,V代表分子黏性耗散項;Cf,T代表湍動能生成項;Cf,G1表示平均流向動能的變化項;Cf,G2表示流向黏性正應力不均勻性項,Cf,G3表示流向雷諾正應力不均勻性項,兩者合起來稱為流向不均勻項.其中.

在對湍流邊界層摩阻空間分布分析的基礎上,本節繼續采用摩阻分解式(4)來分析邊界層中摩阻生成機制.表3 比較了采用摩阻分解技術得到的摩阻平均值(Cf,decomp)、直接采用時間平均得到的摩阻值(Cf,avg)以及兩者之間的相對誤差.相對誤差的定義為error=(Cf,decomp-Cf,avg)/Cf,avg.結果顯示,本文實現的摩阻分解具有較高的精度,相對誤差均在1.5%以下.

表3 摩阻分解公式相對誤差Table 3 Relative error of decomposition formula

表4 給出了式(4)中不同流動和輸運過程對摩阻產生的貢獻比例.可以看出,分子黏性耗散項(Cf,V)和湍動能生成項(Cf,T) 所占比例最大,基本上在40%以上;其次是平均流向動能的變化項(Cf,G1),比例在10%;流向黏性正應力不均勻性項(Cf,G2)和流向雷諾正應力不均勻性項(Cf,G3)占比很小,基本上可以忽略.對比不同焓值算例可以發現,高焓值情況下,Cf,V和Cf,T的占比略有下降,減少的這部分貢獻基本上由Cf,G1補上.摩阻分解分析表明,湍流邊界層中對摩阻生成起主要作用的流動過程是Cf,V和Cf,T.下面將進一步分析兩者產生的原因以及高低焓值邊界層摩阻產生差異的原因.

表4 摩阻分解公式中各項貢獻Table 4 Contributions of different terms to skin friction

2.4 分子黏性耗散項對壁面摩阻的貢獻分析

分子黏性耗散過程是壁面摩阻產生的主要流動過程之一.圖8 給出了分子黏性耗散項預乘積分函數沿法向的分布.兩算例的分布趨勢較為一致,約在y≈0.015δ(y+≈10)處達到峰值.且分子黏性耗散項的主要作用區域在y≤0.1δ(y+≤30)的近壁區,而在外區基本沒有貢獻.這與文獻中基于槽道、平板湍流[21-22]的分析結果一致.高焓會使得y＜0.015δ 以內的近壁區分子黏性耗散略大于低焓流動;而在y＞0.015δ的區域中略小于低焓流動.

圖8 黏性耗散項預乘積分函數Fig.8 Pre-multiplied integrand of viscous dissipation terms

進一步可以將分子黏性耗散項分解為平均部分Cf,V,m和脈動部分Cf,V,f[21]

表4 中給出了兩部分的占比.可以發現,Cf,V,m占據了分子黏性耗散項的絕大部分,占比高達95%以上.這就說明焓值對分子黏性耗散項的影響主要通過Cf,V,m項體現.

將Cf,V,m的積分函數分為兩部分:和,分別體現了流動特性和熱力學特性的作用.圖9給出了兩者沿壁面法向的分布,并給出了溫度的分布,均采用無窮遠處參數作為無量綱參考值.從流動特性可以看出,高焓使得近壁區(y＜0.015δ)的?u/?y增大(?v/?x較小,影響可以忽略);而對外區的流動特性影響很小.從熱力學特性可以看出,高焓對整個邊界層內的黏性系數都有較大影響,基本上均低于低焓算例的黏性系數.黏性系數主要受到溫度的影響,兩者的分布具有明顯的一致性.當高焓時,邊界層內發生空氣化學非平衡過程,離解反應會吸收一部分熱量,使得高焓邊界層內的黏性系數低于低焓邊界層.綜合流動特性和熱力學特性的分布,可以發現,在y＜0.015δ 時,流動特性起到主導作用,導致高焓流動黏性耗散項貢獻略大于低焓流動;而在y＞0.015δ 時,熱力學特性起到主要作用,導致低焓流動黏性耗散項貢獻略大于高焓流動.

圖9 Cf,V,m積分函數分解及溫度分布Fig.9 Integrand parts of Cf,V,mand the temperature distribution

2.5 湍動能生成項對壁面摩阻的貢獻分析

湍動能生成項是壁面摩阻產生的另一個重要部分.圖10 給出了湍動能生成項預乘積分函數沿法向的分布.兩算例的分布趨勢較為一致;但低焓流動時,峰值位置更加遠離壁面.結合圖4 中雷諾應力分布基本一致的情況來看,峰值位置的外移應該來自于的分布不同;結合圖9 來看,高焓流動在預乘積分函數附近的要大于低焓流動.兩算例在邊界層外區均出現了二次峰值,這與較高雷諾數下的槽道流動[21]一致.二次峰值的產生是由于較高雷諾數下,邊界層外區存在大尺度湍流運動導致.大尺度湍流運動會引起外區產生大量湍動能,從而通過摩阻的湍動能生成項影響摩阻;且外區大尺度湍流流動對內區小尺度湍流運動具有調制作用,促發湍動能在內外區的重新分配[35].

圖10 湍動能生成項預乘積分函數Fig.10 Pre-multiplied integrand of TKE production terms

為了進一步分析湍動能生成項中起主導作用的流動事件,本節采用象限分析技術[36-37]將邊界層中的主導流動事件與摩阻分解相結合以分析不同流動事件對摩阻中湍動能生成項的影響.

在邊界層中,可以按照脈動的正負將流向速度脈動和法向速度脈動分解到4 個象限中.第1 象限:u′′＞0,v′′＞0表示外向運動(outward motion);第2 象限:u′′＜0,v′′＞0表示上拋運動;第3 象限:u′′＜0,v′′＜0表示內向運動(inward motion);第4 象限:u′′＞0,v′′＜0表示下掃運動.圖11 給出了Q1～Q4 事件在邊界層內發生的概率.不同焓值下不同事件的概率分布曲線是相似的.Q2 和Q4 發生的概率顯著大于Q1 和Q3.

圖11 Q1～Q4 事件發生概率Fig.11 Probability distribution of Q1～Q4 events

圖12 中給出了不同事件對應的雷諾應力與總雷諾切應力的比值.可以看出,不同焓值的比值曲線很好的重合在一起,說明該比值對焓值變化不敏感;Q2 和Q4 事件是雷諾切應力的主要貢獻項,并且在y≈0.03δ以內,Q4 事件的比值要大于Q2 事件;在y≈0.03δ以外,Q2 事件的比值大于Q4 事件,且這一趨勢會在整個邊界層范圍內得到保持.

表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 給出了不同事件的被積函數沿法向的分布,以考察主要起作用的邊界層區域范圍.可以看出,不同焓值下,被積函數的曲線分布趨勢都是相似的,高焓效應僅會稍微改變幅值;Q4 事件的峰值位置要比Q2 事件更為靠內;Q1 和Q3 在整個邊界層范圍內的值都是負值.

圖13 Q1～Q4 湍流生成項預乘積分函數Fig.13 Pre-multiplied integrand of TKE production terms produced by Q1～Q4

3 結論

本文通過采用直接數值模擬方法,研究了高焓湍流邊界層壁面摩阻產生機制,并通過設置同等邊界層參數下的低焓算例,對比分析了高溫化學非平衡效應對摩阻產生的影響,可以得到以下結論.

(1)不同焓值下的摩阻脈動條帶尺寸存在差異.高焓流動的條帶尺寸在流向和展向兩個方向上均小于低焓流動.

(2)湍流邊界層中對摩阻生成起主要作用的流動過程是分子黏性耗散項Cf,V和湍動能生成項Cf,T.

(3) 分子黏性耗散項Cf,V的主要作用區域在y≤0.1δ的近壁區.焓值的影響主要通過Cf,V的平均部分Cf,V,m項體現.且在y＜0.015δ 時,流動特性(速度梯度)起到主導作用,導致高焓流動Cf,V,m略大于低焓流動;而在y＞0.015δ 時,熱力學特性(黏性系數)起到主要作用,導致低焓流動Cf,V,m略大于高焓流動.

(4)邊界層中的上拋和下掃運動是影響摩阻分解湍動能生成項Cf,T的主導事件,外向運動和內向運動對Cf,T的影響很小,且起到負貢獻.