周航，陳方

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

0 引言

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

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

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

1.1 控制方程

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

(1)

(2)

式中：Ck是湍動能在平均運動軌跡上的增長率；Π是湍動能產(chǎn)生率項；D是湍動能擴散項，包含壓力速度相關的擴散、湍流的擴散及分子黏性產(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)格無關性驗證

2 計算結(jié)果分析

2.1 流場NPTKE分布

圖3表示了特征區(qū)域(x=0.32～0.53m)內(nèi)壓力分布云圖和密度拉普拉斯算子數(shù)值陰影圖。翼型位置產(chǎn)生的斜激波與上壁面邊界層相互作用導致流動分離，分離點下游是流線包圍回流泡，在分離點和再附點分別產(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與應變率張量關系

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

(4)

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

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

為探究應變率張量對湍動能產(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點。氣體壓縮性效應增強，流場中湍流脈動逆向傳遞到平均流動的湍動能總量減小，NPTKE效應減弱。

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

3 結(jié)語

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

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

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

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

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