可壓縮湍流邊界層近壁區(qū)馬蹄渦的演化

史萬里，葛 寧，陳 林，唐登斌

(1.南京航空航天大學能源與動力學院，江蘇南京 210016;2.南京航空航天大學航空宇航學院，江蘇南京 210016)

0 引言

在流體力學的研究領(lǐng)域中，湍流邊界層問題是一個至今尚未解決的重大難題，而擬序結(jié)構(gòu)問題則對于湍流邊界層的進一步深入研究和探索至關(guān)重要［1］，尤其在湍流邊界層近壁區(qū)域，流動的許多宏觀特性，如湍流的生成和維持、湍流能量的輸運和耗散、阻力減小、熱量傳遞等都與之密切相關(guān)［2］。大量實驗研究與數(shù)值模擬結(jié)果表明，流向渦結(jié)構(gòu)是擬序結(jié)構(gòu)的一種重要典型特征［3］，馬蹄渦(或發(fā)卡渦)更受到了廣泛的重視［4－5］，其生成與演化規(guī)律的研究，能夠為湍流邊界層內(nèi)的流動現(xiàn)象提供很多合理的解釋，并有助于深入理解和解決工程中遇到的流動問題。

在湍流邊界層中，馬蹄渦與發(fā)卡渦在結(jié)構(gòu)上稍有不同，通常在動量厚度雷諾數(shù)大于800時稱為發(fā)卡渦，低于此值時渦頭部分通常較寬，因此稱為馬蹄渦［6－7］。以前人們認為這種馬蹄渦(或發(fā)卡渦)結(jié)構(gòu)可用對稱的一對反向旋轉(zhuǎn)的流向渦對來表示［8］。最近的研究發(fā)現(xiàn)［2，6－7，9］，湍流近壁區(qū)的擬序結(jié)構(gòu)不僅存在對稱的流向渦對，而且絕大多數(shù)流向渦對更是以非對稱形式存在，這種結(jié)構(gòu)有時甚至是以一個流向渦的形式出現(xiàn)，這種現(xiàn)象引起人們的很大興趣。

在不同的數(shù)值方法中，大渦模擬(LES)方法相比于雷諾平均方法(RANS)，能夠顯示出流場的細節(jié)和瞬態(tài)過程［10］，而對計算資源的需求大大小于直接數(shù)值模擬(DNS)方法，所以本文是采用LES方法研究三維馬蹄渦的結(jié)構(gòu)。

對于采用LES方法的湍流模擬，如果直接從一般層流開始轉(zhuǎn)捩到湍流，那么計算湍流場則需要附加很大的計算量。因此本文采用進口給定三維擾動波方法［11－12］，作為大渦模擬的湍流進口邊界的給定方法，模擬了可壓縮三維槽道湍流流動。根據(jù)計算結(jié)果，與DNS解的數(shù)據(jù)進行了比較，并對近壁湍流的馬蹄渦結(jié)構(gòu)的演化進行了深入的研究和分析。

1 控制方程

對預處理下［13－14］的可壓縮流 Navier-Stokes方程，采用Faver提出的密度加權(quán)過濾方法，得到曲線坐標系下，無量綱通量形式的LES控制方程為:

式中，上標“－”表示系統(tǒng)平均，“～”表示細網(wǎng)格過濾，轉(zhuǎn)換矩陣 M= ?Q/?q，原始變量為為守恒變量。Re為特征雷諾數(shù)表示亞格子粘性應力為亞格子熱流通量，雅克比矩陣轉(zhuǎn)換矩陣

需要注意的是，Γ=diag{1 1 1 1 β(Mr)}為預處理矩陣［13－14］，預處理參數(shù) β 的選取:若 Mr＜1.0，β(Mr)=Mr2;若Mr≥1.0，β(Mr)=1.0。其中為特征馬赫數(shù)。這里采用預處理技術(shù)，主要是為了改善邊界層近壁區(qū)低馬赫數(shù)流動情況下的計算格式的穩(wěn)定性以及加快計算結(jié)果的收斂性。

2 數(shù)值方法

將對控制方程使用有限體積法進行離散，粘性項采用二階中心格式，對流項采用六階精度對稱WENO格式，時間推進采用三步Runge-Kutta法進行精確求解。其亞格子模型、計算域和進口邊界條件分述如下:

2.1 亞格子模型

建立合理的亞格子模型是大渦模擬的關(guān)鍵，Smagorinsky模型是廣泛使用的一種亞格子模型，表達式如下［15］:

相應的亞格子熱流通量為:

但是該模型存在著某些缺陷，如系數(shù)C實際最佳值取決于運動流體的性質(zhì)與狀態(tài)而不是常數(shù)，且模型未體現(xiàn)流體壓縮性的影響等，因此本文采用Germano提出的動態(tài)亞格子模型，具體如下［15］:

動態(tài)模型是動態(tài)確定亞格子渦粘系數(shù)C，對湍流場做粗細兩次過濾，并假設(shè)粗網(wǎng)格上最小脈動產(chǎn)生的應力等于粗細網(wǎng)格分別過濾產(chǎn)生的亞格子應力之差，在此基礎(chǔ)上用最小二乘法推導出亞格子模型中的系數(shù)C如下:

其中:

式中，上標“^”表示粗網(wǎng)格過濾，“＜ ＞”表示在統(tǒng)計方向上平均量。

2.2 計算域

計算域見圖1(a)，計算中選取進口湍流邊界層的位移厚度δ*作為特征長度，計算域無量綱尺度取為50×10×8，分別對應流向x、法向y和展向z的長度。三個方向上計算網(wǎng)格點數(shù)分別取為192×96×64，在流向和展向上計算網(wǎng)格平均分布，而在法向上的近壁區(qū)域內(nèi)進行加密。按照進口邊界條件確定的網(wǎng)格尺度為△x+=13.0，△z+=6.25，第一層網(wǎng)格的y+=0.23，粘性底層區(qū)域(即 y+=10)內(nèi)有22個網(wǎng)格點，x-y面局部網(wǎng)格如圖1(b)。為了方便并行計算將網(wǎng)格平均分為4份，每份在一個CPU上計算，編號依次從0到3，如圖1(a)，利用MPI來傳遞各區(qū)之間信息交換，并行計算。

流場初始化在邊界層內(nèi)采用Spalding湍流平均速度剖面。在計算域的出口給定反壓，兩邊的展向位置給定周期性邊界條件，在上邊界采用對稱邊界條件，平板壁面處采用絕熱無滑移邊界條件。因此圖1中所給出的僅是法向y上的半個槽道區(qū)域。

圖1 計算域及網(wǎng)格示意圖Fig.1 Schematic of the parallel computational domain decomposition

2.3 進口邊界條件

在計算域進口采用了進口給定三維擾動波方法，這是基于目前湍流研究者普遍認同的湍流是擬序結(jié)構(gòu)的組合這一假設(shè)，在進口產(chǎn)生包括時間和空間變化的擬序結(jié)構(gòu)，即擾動波。這種方法相比于其它常用的湍流進口給定方法，是簡單可行的，其進口發(fā)展段更短，更適用于工程上使用。

進口邊界條件，主要是確定進口處的瞬時速度。進口處瞬時速度通常由三部分組成:瞬時值=平均值+擾動值+隨機噪音。平均值是采用Spalding湍流平均速度剖面，隨機噪音的最大幅值不超過平均自由流速度的1%。

對于擾動值的確定，本文直接給定三維波。通常湍流邊界層按動態(tài)特征可分為內(nèi)層和外層，這種方法是認為內(nèi)層中存在低速條帶，用一道三維波表示，而把外層中的擾動波用三維基波和亞諧波組合的渦結(jié)構(gòu)［11－12］來表示。其公式［11］及參數(shù)的選擇如下:

式中，展向速度w按照自由散度條件由u和v解出，基本參數(shù)參照表1。表中的ωj是脈動頻率，βj是展向波數(shù)，φj是相位。

表1 湍流邊界層進口湍流脈動參數(shù)Table 1 Parameters of inflow turbulent fluctuations for the turbulent boundary-layer

3 結(jié)果與討論

湍流邊界層的結(jié)構(gòu)是三維的，邊界層外緣具有弱的自由剪切層的性質(zhì)，內(nèi)外層之間有相互作用，內(nèi)層又受到壁面的抑制，因此湍流邊界層內(nèi)的流動結(jié)構(gòu)非常復雜。為了檢驗計算結(jié)果的可靠性，這里首先把計算結(jié)果與Spalart的DNS解［16］等進行了比較，然后對LES計算得到的馬蹄渦結(jié)構(gòu)進行分析。文中算例的馬赫數(shù)Mr=0.6，基于動量厚度和自由流速度的雷諾數(shù)為670(對應于基于位移厚度和自由流速度的雷諾數(shù)為1000)。

3.1 時間平均速度和表面摩擦系數(shù)

圖2(a)是LES計算得到的平均速度剖面圖，同時給出了Spalart的DNS的結(jié)果［16］。由圖可以看出，本文的LES計算結(jié)果與相同雷諾數(shù)下的DNS精確解符合的很好，說明了所采用的進口剖面在經(jīng)過很短的發(fā)展段之后，已發(fā)展成穩(wěn)定的湍流。

圖2(b)是湍流脈動速度，即湍流強度分布圖，其縱坐標是脈動速度的絕對值(已用自由流速度無量綱化)。該圖也將LES計算值與Spalart的湍流DNS解進行了比較，發(fā)現(xiàn)各個方向上的湍流強度的幅值與DNS解基本一致，僅在流向脈動速度的幅值稍有偏離，其原因可能是DNS解為不可壓縮流計算的結(jié)果，而本文的LES解是預處理下可壓縮流計算的結(jié)果。

圖2(c)是時間平均和展向平均的表面摩擦系數(shù)Cf。可以看到從計算域入口僅經(jīng)過很短的入口過渡距離，大約在流向Rex為139000時，計算的Cf就符合了湍流的經(jīng)驗值。在圖2中，湍流摩擦系數(shù)Cf并不是完全符合湍流的經(jīng)驗值，而是存在波動，這是由于湍流的非定常性引起的。

3.2 湍流邊界層的馬蹄渦結(jié)構(gòu)

本節(jié)對計算得到的湍流邊界層內(nèi)馬蹄渦結(jié)構(gòu)的演化進行分析，其中給出的時間是以進口標準周期T為單位。因為渦核處存在強的渦量及最小壓力值，所以采用拉普拉斯壓力［17］P，kk=(ωiωi/2 － SijSij)的正的等值面是很容易識別擬序結(jié)構(gòu)的，以下圖中渦結(jié)構(gòu)顯示均為P，kk=0.6 的等值面。

3.2.1 馬蹄渦的形成和演化

圖3是t=6.62T時的速度云圖和流場渦結(jié)構(gòu)，圖中深藍色區(qū)域為低速區(qū)域。圖中可以清楚地看到流場存在很多的流向渦結(jié)構(gòu)，這些流向渦結(jié)構(gòu)與低速條帶是相對應的，流向渦在x-z平面內(nèi)是左右擺動、交替出現(xiàn)的，因而高低速條紋是左右擺動而不單純是流向的。這和文獻［9］和文獻［18］中的結(jié)論一致。

圖4中顯示了對稱流向渦對形成馬蹄渦結(jié)構(gòu)的過程，這個結(jié)構(gòu)在后面的圖8中的標號為E1。正如圖3中提到的，流向渦結(jié)構(gòu)左右擺動，在圖4中可以看到，當兩條流向渦相互靠近時，兩渦腿的頭部之間會形成橋狀的展向渦結(jié)構(gòu)，進而合并融合而成的，這種對稱產(chǎn)生的馬蹄渦結(jié)構(gòu)形成過程與文獻［5］的DNS解一致。

圖2 LES計算結(jié)果與DNS解和經(jīng)驗公式的對比Fig.2 Comparison of result by LES，DNS and experiential formula

圖3 流場渦結(jié)構(gòu)與低速條帶俯視圖Fig.3 Top view of vortex structures and low-speed streaks in flow field

圖4 對稱流向渦形成的馬蹄渦結(jié)構(gòu)Fig.4 Generation of the horseshoe-shaped structures by symmetric streamwise vortex

兩條渦腿相互靠近會形成帶兩條長渦腿的馬蹄渦結(jié)構(gòu)，文中的計算結(jié)果也清楚地展示了這一點，但是如果兩條渦腿距離較遠，則馬蹄渦結(jié)構(gòu)的渦頭不會封閉，出現(xiàn)中間開口的現(xiàn)象，這樣最終形成由單側(cè)流向渦腿形成的馬蹄渦結(jié)構(gòu)。圖5顯示了F1結(jié)構(gòu)由沿流向傾斜的流向渦發(fā)展為馬蹄渦的過程，F(xiàn)1馬蹄渦的渦腿是在右側(cè)。圖5(a)中，F(xiàn)1為流向渦，其圖中顯示的流向渦的渦頭部分在自誘導的作用下向上抬頭，同時向一側(cè)的展向傾斜。圖5(b)、圖5(c)中可以看到這個向展向傾向的渦頭沿展向拉伸，因此渦頭部分的渦管變細，沿展向的長度變長。在圖5(c)、圖5(d)中可以看到，渦頭的展向部分在拉伸的同時，在擾動的作用下，形成小的馬蹄渦。在圖5(e)～圖5(h)可以看到這個單側(cè)渦腿的馬蹄渦較短渦腿的形成過程，它正如Robinson所述的帶有單側(cè)渦腿的馬蹄渦形成過程那樣［6］(如圖6(a))，馬蹄渦的頸部向壁面移動，逐漸下降并拉伸成一個加長的渦腿。在圖5(e)、圖5(f)，可以看到在馬蹄渦的短渦腿側(cè)出現(xiàn)了細碎的渦管a1，這是在頸部的移動、下降和拉伸時，這個被加長的渦腿的渦量較弱，在圖中不能完全顯示出來，因此看上去這個渦腿是斷開的。圖6(b)是Delo等人［7］在實驗中得到了類似馬蹄渦和流向渦的結(jié)果。圖中標號為作者對流場中可能的渦結(jié)構(gòu)的注釋，其中V6是指流向渦，V4是指可能的單腿馬蹄渦，這個馬蹄渦的渦腿在左側(cè)。

圖5 單側(cè)流向渦形成的較小馬蹄渦對結(jié)構(gòu)Fig.5 Generation of the horseshoe-shaped vortex pair by the one-side streamwise vortex

圖6 湍流邊界層中馬蹄渦Fig.6 The horseshoe-shaped vortex in the turbulent boundary layer

由于文中進口邊界條件所產(chǎn)生的是反向旋轉(zhuǎn)的流向渦對，因此計算結(jié)果中這種單側(cè)渦腿形成的馬蹄渦結(jié)構(gòu)總是成對出現(xiàn)的，如在后面的圖8中標識的F1和F2及F3和F4結(jié)構(gòu)。圖7(a)是 t=5.94T時刻，F(xiàn)1和F2馬蹄渦的俯視圖和側(cè)視圖，俯視圖中的背景為速度云圖，圖中深藍色區(qū)為低速區(qū)。在俯視圖中可以看到F1和F2結(jié)構(gòu)都騎跨在低速條帶上，在側(cè)視圖中可以看到馬蹄渦頭部的抬起。

圖7(b)是t=5.98T時刻的典型馬蹄渦結(jié)構(gòu)渦腿處速度矢量場?？梢钥吹?，F(xiàn)1結(jié)構(gòu)中環(huán)繞馬蹄渦結(jié)構(gòu)兩條流向渦腿的流體是反向旋轉(zhuǎn)的，因此在兩條流向渦腿之間的流體存在強烈的上升，而在兩條渦腿外側(cè)，流體下掃。同時還注意到，在兩條渦腿之間還存在一些碎的流向渦結(jié)構(gòu)，而流體在流向渦兩側(cè)同樣也形成上升和下掃運動。

3.2.2 亞諧波影響下的馬蹄渦結(jié)構(gòu)

圖7 單腿馬蹄渦細節(jié)圖Fig.7 Particular view of one-side horseshoe vortices

圖8顯示的是流場中存在的擬序結(jié)構(gòu)，為更清楚起見，把計算結(jié)果沿展向擴展了一個周期。我們注意到，圖8中計算域的前面部分(即圖中點劃線之前)為進口發(fā)展段，這里的擬序結(jié)構(gòu)只是研究在這之后(即點劃線之后)的結(jié)果。圖中顯示了在湍流場中存在很多馬蹄渦，且這些馬蹄渦結(jié)構(gòu)的形狀并不是規(guī)則一致的，也不存在嚴格的對稱性。文獻［2］、文獻［9］指出，這是由于湍流的非定常性影響，即湍流脈動會扭曲渦管的發(fā)展，從而很難在湍流邊界層中得到完整一致的馬蹄渦。

圖8(a)、圖8(b)分別顯示了兩個不同時間的流場擬序結(jié)構(gòu)圖，圖中標示的大寫英文字母序列表示在亞諧波影響下出現(xiàn)的馬蹄渦結(jié)構(gòu)。如圖8中所示，可以發(fā)現(xiàn)近壁區(qū)存在幾種對稱流向渦對形成的馬蹄渦的復雜結(jié)構(gòu)(如E1到E3所標識的擬序結(jié)構(gòu)，這種結(jié)構(gòu)的渦頭部分的形成如圖4中所示)，但是基本的流動特征還是表示為亞諧波失穩(wěn)引起的對稱產(chǎn)生的馬蹄渦交錯出現(xiàn)。即沿流向x方向看，前一個對稱產(chǎn)生的馬蹄渦結(jié)構(gòu)中心位置與后一個對稱產(chǎn)生的馬蹄渦結(jié)構(gòu)中心位置沿流向并不是在一條直線上的，前排結(jié)構(gòu)與后排結(jié)構(gòu)的展向位置呈叉排交錯分布。另外，在圖8中還可以發(fā)現(xiàn)，在兩組對稱產(chǎn)生的馬蹄渦結(jié)構(gòu)之間通常存在著由單側(cè)渦腿產(chǎn)生的較小馬蹄渦結(jié)構(gòu)(如F1到F4所標識的擬序結(jié)構(gòu)，這種結(jié)構(gòu)的渦頭部分的形成如圖5中所示)。此外，還可以看到在兩組結(jié)構(gòu)之間沿展向不均勻分布著流向渦結(jié)構(gòu)，這些流向渦結(jié)構(gòu)總是左右擺動，而不是單純的沿流向向前發(fā)展。

圖8 湍流場中存在的渦結(jié)構(gòu)Fig.8 Vortex structures at different time instants

由上面的論述可知，流向渦結(jié)構(gòu)在向下游發(fā)展過程中，既能單獨演化，又有相互作用。單獨演化時往往會出現(xiàn)單側(cè)渦腿產(chǎn)生的較小馬蹄渦，相互作用則可發(fā)現(xiàn)兩個流向渦結(jié)構(gòu)相互靠近合并成較大馬蹄渦結(jié)構(gòu)。同時這些馬蹄渦形狀結(jié)構(gòu)并非完全對稱的。

4 結(jié)論

本文采用動態(tài)亞格子模型和預處理的大渦模擬方法對槽道可壓縮流湍流進行了精確的數(shù)值模擬，得到了滿意的結(jié)果。通過對結(jié)果的分析比較，可以得出以下結(jié)論:

(1)文中LES計算得到的湍流平均速度剖面、湍流強度及平均表面摩擦系數(shù)，與精確的DNS解和經(jīng)驗公式符合得很好，驗證了計算結(jié)果的可靠性;

(2)計算結(jié)果表明，湍流場內(nèi)的流向渦結(jié)構(gòu)與低速條帶是相對應的，流向渦在x-z平面內(nèi)會左右擺動、交替出現(xiàn)，因而快慢條紋會左右擺動、而不單純是流向的;

(3)結(jié)果清楚地展示了湍流場中的流向渦演化成馬蹄渦的過程，這些流向渦向下游的發(fā)展過程中，既可以相互作用，由對稱渦腿產(chǎn)生馬蹄渦結(jié)構(gòu)，又能單獨演化成單側(cè)渦腿的馬蹄渦結(jié)構(gòu)。特別是，在亞諧波的作用下，湍流場中可能會出現(xiàn)的復雜馬蹄渦結(jié)構(gòu)，這些對稱渦腿產(chǎn)生的馬蹄渦相互之間是交錯出現(xiàn)的，這些結(jié)構(gòu)由于湍流擾動，會出現(xiàn)形狀也不是嚴格對稱的現(xiàn)象。

［1］RINGUETTE M J，WU M，MARTIN M P.Coherent structures in direct numerical simulation of turbulent boundary layers at Mach 3［J］.Journal of Fluid Mechanics，2008，594:59-69.

［2］連祺祥.湍流邊界層擬序結(jié)構(gòu)的實驗研究［J］.力學進展，2006，36(3):373-388.

［3］張楠，陸利蓬，段真真，等.湍流邊界層內(nèi)準流向發(fā)卡渦生成的數(shù)值模擬［J］.應用數(shù)學和力學，2008，29(1):13-20.

［4］ZHOU J，ADRIAN R J，BALACHANDAR S，et al.Mechanisms for generating coherent packets of hairpin vortices in channel flow［J］.Journal of Fluid Mechanics，1999，387:353-396.

［5］陳林，唐登斌，劉小兵，等.邊界層轉(zhuǎn)捩過程中環(huán)狀渦和尖峰結(jié)構(gòu)的演化［J］.中國科學G輯，2009，39(10):1520-1526.

［6］ROBINSON S K.Coherent motions in the turbulent boundary layer［J］.Ann Rev Fluid Mech，1991，23:601-639.

［7］DELO C J，KELSO R M，SMITS A J.Three-dimensional structure of a low-Reynolds-number turbulent boundary layer［J］.Journal of Fluid Mechanics，2004，512:47-83.

［8］CANTWELL B J.Organized motions in turbulent flows［J］.Annual Review of Fluid Mechanics，1981，13:457-515.

［9］周恒，熊忠民.湍流邊界層近壁區(qū)擬序結(jié)構(gòu)起因的研究［J］.中國科學 A 輯，1994，24(9):941-948

［10］MOIN P.Advances in large eddy simulation methodology for complex flows［J］.International Journal of Heat and Fluid Flow，2002，23:710-720

［11］SANDHAM N D，YAO Y F，LAWAL A A.Large-eddy simulation of transonic turbulent flow over a bump［J］.Inernational Journal of Heat and Fluid Flow，2003，24:584-595.

［12］周恒，平板湍流邊界層底層的不穩(wěn)定波［J］.力學學報，1988，20(6):481-488.

［13］BRILEY W R，TALTLOR L K，WHITFIELD D L.Highresolution viscous flow simulations at arbitrary Mach umber［J］.Journal of Computational Physics，2003，184(1):79-105.

［14］史萬里，葛寧，盛春華.計算全馬赫數(shù)下粘性流動的隱式算法［J］.空氣動力學學報，2009，27(5):566-571.

［15］GERMANO M，PIOMELLI U，MOIN P，et al.A dynamic subgrid-scale eddy viscosity model［J］.Physics of Fluids，1991，3:1760-1765.

［16］SPALART P R.Direct simulation of a turbulent boundary layer up to Reθ=1410［J］.Journal of Fluid Mechanics，1988，187:61-98.

［17］TYAGI M，ACHARYA S.Large eddy simulation of film cooling flow from an inclined cylindrical jet［J］.ASME，2003，125:734-742.

［18］SCHOPPA W，HUSSAIN F.Numerical study of near-wall coherent structures and their control in turbulent boundary layers［A］.6th International Conference on Numerical Methods in Fluid Dynamics，Volume 515［C］.ISBN 978-3-540-65153-6.Springer-Verlag，1998.