超聲速槽道流動中的湍動能負產(chǎn)生率研究

周航，陳方

(上海交通大學 航空航天學院，上海200240)

0 引言

湍動能輸運方程中的產(chǎn)生項表示雷諾應(yīng)力通過平均運動的變形率向湍流脈動輸入的平均能量，一般情形下為正值。但是在充分發(fā)展的彎曲槽道流動研究[1]中，由于流場局部出現(xiàn)方向相反的雷諾應(yīng)力和切應(yīng)力，首次發(fā)現(xiàn)湍動能負產(chǎn)生率(negative production of turbulent kinetic energy, NPTKE)現(xiàn)象。后續(xù)研究指出NPTKE也出現(xiàn)在許多其他流動類型中。如射流流動中，射流速度降低后會在軸線處出現(xiàn)NPTKE區(qū)域，雷諾應(yīng)力張量各向異性減弱[2]，脈沖作用會加劇這種現(xiàn)象[3-4]。圓柱組繞流中，相鄰圓柱間的流體加速促進湍動能產(chǎn)生率由正轉(zhuǎn)負，并使該區(qū)域內(nèi)湍流強度變大[5]。環(huán)流截面中的NPTKE現(xiàn)象平衡了速度峰值兩側(cè)的曲率，使流動交匯處的平均速度型面保持連續(xù)[6]。可壓縮均勻剪切流中，條紋結(jié)構(gòu)使NPTKE增大，而傾斜渦管結(jié)構(gòu)使其減小[7]。立方T形管流動[8]中NPTKE區(qū)域集中在回流區(qū)，可能是由于邊界層中壓力梯度分布不均導(dǎo)致。ABE H等[9]觀測到NPTKE總是出現(xiàn)在正流線曲率半徑處。振蕩邊界層流動中，湍流結(jié)構(gòu)受到剪切作用變形會使湍動能產(chǎn)生率變?yōu)樨撝礫10-11]。LIBERZON A等對流實驗[12]說明浮力會導(dǎo)致湍動能負產(chǎn)生率，與湍流標量通量和標量梯度[13-14]之間的正負關(guān)系有關(guān)，流場拉伸特性會促進NPTKE[15]。對于后臺階流動[16]，在較長的流向尺度上發(fā)生基于雷諾應(yīng)力的非線性相互作用才可能出現(xiàn)NPTKE現(xiàn)象。強非對稱流場計算運用傳統(tǒng)的k-ε等梯度假設(shè)湍流模型無法準確計算出NPTKE[2-17]。GRETLER W等[18]和WALLIN S等[19]提出了擴展的雷諾應(yīng)力模型，可以計算出NPTKE的準確數(shù)值。在壓氣機和噴管[20]等重要的航空發(fā)動機部件中，通過控制NPTKE的變化來實現(xiàn)氣流品質(zhì)調(diào)節(jié)，有利于飛行器的高效穩(wěn)定飛行。

以上研究對象為低速不可壓縮流體，而超聲速流動中介質(zhì)的可壓縮性不可忽略，存在的激波也會對NPTKE產(chǎn)生影響。本文分析采用雷諾應(yīng)力模型(reynolds stress model, RSM)對放置翼型的超聲速槽道流數(shù)值模擬，研究超聲速槽道流中的NPTKE現(xiàn)象分布特征及成因、入口馬赫數(shù)的影響規(guī)律。

1 計算方法及數(shù)值驗證

1.1 控制方程

可壓縮流動遵循Reynolds平均連續(xù)方程、動量方程和能量方程，RSM模型通過對湍流渦黏性系數(shù)的各項異性?；?，對復(fù)雜流動具有更好的預(yù)測潛力。雷諾應(yīng)力輸運方程如式(1)所示。

(1)

(2)

式中：Ck是湍動能在平均運動軌跡上的增長率；Π是湍動能產(chǎn)生率項；D是湍動能擴散項，包含壓力速度相關(guān)的擴散、湍流的擴散及分子黏性產(chǎn)生的擴散；ε是湍動能耗散項。

在二維直角坐標系下，湍動能產(chǎn)生率項Π表達式如下：

(3)

1.2 計算模型

根據(jù)劉宇陸等[22]所研究的實驗三維模型簡化為二維計算模型。如圖1所示，槽道流動的物理模型中試驗段總長度530mm，槽道高度25mm，在距入口330mm處的下壁面放置最大厚度9mm的軸對稱翼型。超聲速空氣入口采用均勻來流條件，入口總壓P0=150kPa，總溫T0=300K，馬赫數(shù)Ma=4。

圖1 計算模型幾何參數(shù)

1.3 數(shù)值驗證

圖2結(jié)果表明采用15萬網(wǎng)格節(jié)點數(shù)的算例計算滿足精度要求且計算量較小，最為合理。

圖2 網(wǎng)格無關(guān)性驗證

2 計算結(jié)果分析

2.1 流場NPTKE分布

圖3表示了特征區(qū)域(x=0.32～0.53m)內(nèi)壓力分布云圖和密度拉普拉斯算子數(shù)值陰影圖。翼型位置產(chǎn)生的斜激波與上壁面邊界層相互作用導(dǎo)致流動分離，分離點下游是流線包圍回流泡，在分離點和再附點分別產(chǎn)生分離反射激波和再附反射激波。特征區(qū)域內(nèi)湍動能產(chǎn)生率Π和NPTKE分布如圖4所示，NPTKE出現(xiàn)位置在收縮后的槽道主流中，分布較為離散。邊界層和分離區(qū)中湍動能產(chǎn)生率一般為正值，在兩道反射激波之間NPTKE值較大。

圖3 特征區(qū)域流動結(jié)構(gòu)

圖4 特征區(qū)域參數(shù)分布

利用數(shù)學統(tǒng)計方法，對湍動能產(chǎn)生率作概率密度分布(probability density function, PDF)圖5。從圖中可以看出NPTKE的峰值小于正產(chǎn)生率，流場中NPTKE出現(xiàn)概率也小于正產(chǎn)生率，這也與低速不可壓流動中的結(jié)論相吻合[1]。

圖5 湍動能產(chǎn)生率與概率密度分布圖

2.2 激波處湍動能產(chǎn)生率變化規(guī)律

對微元的運動進行分析(圖7)。在圖7(a)下壁面附近激波位置處，由于渦旋運動激波前的高速流體微元流向激波后的低速流層，該微元從法向速度近乎為0流入激波后存在法向速度的流層中，可以認為法向速度脈動v″<0，而原微元的流向速度大于周圍流體，則流向速度脈動u″>0，

圖7 流體微元運動情況

2.3 NPTKE與應(yīng)變率張量關(guān)系

湍動能產(chǎn)生率Π的計算式還可以轉(zhuǎn)化為兩種形式，第一種形式如式(4)所示，表示雷諾應(yīng)力張量和應(yīng)變率張量的內(nèi)積。對其夾角余弦值cos(τRe,Sij)進行概率密度統(tǒng)計。

(4)

從圖8可以看出在NPTKE區(qū)域該余弦值均保持為負值，由此可知NPTKE與流場的應(yīng)變率張量有密切的關(guān)系。湍動能正產(chǎn)生率區(qū)域中概率峰值點處的夾角余弦值為0.271，對應(yīng)角度約為-74.3°，與Hanjalic的平面通道非對稱流動實驗[23]中根據(jù)應(yīng)力主軸方向角計算式得到的第二主軸方向-73°基本一致。

圖8 不同湍動能產(chǎn)生率區(qū)域cos(τRe,Sij)PDF圖

為探究應(yīng)變率張量對湍動能產(chǎn)生率Π的貢獻，對式(4)進一步分解，得到第二種表達形式如下：

(5)

圖9 不同湍動能產(chǎn)生率區(qū)域ΠΛk與分布

2.4 入口馬赫數(shù)影響

入口馬赫數(shù)的改變會使激波強度和偏轉(zhuǎn)角度發(fā)生改變，湍動能產(chǎn)生分布也會隨之發(fā)生變化。為了研究其影響規(guī)律，保持算例的幾何結(jié)構(gòu)和入口總溫總壓不變，計算入口馬赫數(shù)在[3,5]范圍以0.5為間隔的其他4組算例。圖10表示隨馬赫數(shù)增大、流場中湍動能負產(chǎn)生率的取值范圍減小，分布更集中于0點。氣體壓縮性效應(yīng)增強，流場中湍流脈動逆向傳遞到平均流動的湍動能總量減小，NPTKE效應(yīng)減弱。

圖10 不同馬赫數(shù)工況下Π與PDF圖

3 結(jié)語

本文使用雷諾應(yīng)力湍流模型對超聲速放置翼型的槽道流動進行數(shù)值模擬，探究湍動能負產(chǎn)生率的變化規(guī)律。得到以下主要結(jié)論：

1) 激波位置處湍動能產(chǎn)生率值始終為正，并取得局部最大值，法向正應(yīng)力與梯度乘積項是決定性因素；

2) 當雷諾應(yīng)力與應(yīng)變率張量的夾角余弦值為負時，會導(dǎo)致局部NTKEP出現(xiàn)；

3) 與亞聲速情況類似，超聲速流場應(yīng)變的拉伸特性會促使NTKEP區(qū)域出現(xiàn)，而壓縮特性則保持湍動能正產(chǎn)生率，整體積分下湍動能產(chǎn)生率恒為正值；

4) 隨著馬赫數(shù)增大，湍動能產(chǎn)生負產(chǎn)生率取值范圍減小，氣體壓縮性效應(yīng)增強，NPTKE效應(yīng)減弱。NPTKE的極大值點位于第一道反射激波后靠近上壁面處。