基于超大渦模擬的翼端間隙流湍流特性與損失機理分析

陳為升 黎耀軍,2 劉竹青,2 楊 魏,2

(1.中國農業大學水利與土木工程學院， 北京 100083； 2.北京市供水管網系統安全與節能工程技術研究中心， 北京 100083)

0 引言

在葉片正背面壓差作用下，軸流式葉輪機械轉子葉頂間隙產生間隙流動，間隙流與葉輪通道內主流相互摻混，在輪緣區形成間隙泄漏渦(Tip leakage vortex, TLV)、間隙分離渦(Tip separation vortex, TSV)和誘導渦(Induced vortex, IV)等復雜的渦旋結構，導致間隙流動損失，影響機組性能[1-2]。在采用小展弦比葉片的葉輪中，間隙流動損失在總流動損失中的占比可達60%～70%[3]，探明間隙區湍流特性和流動損失機理對高性能軸流式葉輪機械研發具有重要意義。

間隙區流動損失包含粘性損失和湍流損失兩部分，前者為流體粘性導致的能量耗散，從湍流能量級串的角度看，后者則為大尺度湍流結構從平均流動中獲取的能量，該部分能量最終經由小尺度湍流結構耗散為熱能。大量實驗和數值模擬研究結果表明，間隙區湍流損失與間隙渦系結構的形成和消散相關，間隙區湍流特性是影響間隙湍流損失的最直接因素。文獻[4]的研究認為湍流損失與湍動能成正比，并建立了間隙湍流損失的計算模型。文獻[5]的研究結果表明，渦輪的葉頂間隙湍流損失集中在TLV區域，而線性葉柵間隙流動試驗結果[6]表明，葉片端部流動分離形成的TSV和葉片吸力面附近的TLV是間隙湍流損失的主要來源，間隙湍流損失中的40%源于TSV的影響，且總的間隙流動損失在葉片端部間隙區沿流向逐漸增大，在葉片尾緣達到最大值后逐漸降低。在間隙湍流損失的產生機理和影響因素方面，文獻[7]發現間隙湍流損失由展向速度脈動主導；而線性葉柵繞流的研究表明，與端壁平行的速度脈動、流向及展向的速度梯度是產生間隙湍流損失的主要原因[2,8]，間隙區強剪切流動中的速度梯度是湍動能產生的重要原因[9]；文獻[10]的研究考慮了間隙寬度的影響，發現間隙寬度影響間隙渦系的生成位置和移動軌跡，進而改變間隙湍流損失；在線性葉柵中，翼型端部適當的開槽處理可以有效降低間隙流動損失[11]。

為優化轉子性能，軸流式旋轉機械葉片端部邊緣常采用倒圓處理，盡管現有研究已證實葉端邊緣倒圓對間隙渦系和泄漏流量產生一定影響[12-13]，但不同葉端間隙形狀下的間隙湍流損失特性及機理并不明確。由于旋轉機械間隙渦結構與翼型端部間隙渦結構具有相似性[14]，為探明不同間隙形狀下的間隙區湍流特性，本文采用可直接求解大尺度湍流結構的超大渦模擬方法(Very large eddy simulation, VLES)對NACA0009型鈍尾緣水翼進行數值模擬，對比分析不同翼端間隙寬度τ和不同翼端邊緣倒圓半徑r對應的翼端間隙渦結構、雷諾應力和湍動能等湍流特征，探明不同翼端間隙形狀的間隙湍流損失機理。

1 計算模型

1.1 NACA0009型水翼及計算域

研究對象為NACA0009型鈍尾緣對稱水翼[15]，計算域如圖1所示。水翼弦長c，計算域長10c，橫斷面為1.5c×1.5c的正方形截面，水翼位于計算域中部，計算域入口距水翼中心5.1c，來流攻角α=10°，基于水翼弦長c和來流速度u∞的雷諾數Re為1.0×106。

圖1 計算域示意圖Fig.1 Schematic of computational domain and coordinate system

本文選用兩種翼端間隙寬度(τ=0.1c，τ=0.02c)和3種翼端邊緣倒圓半徑(r=0，r=0.005c，r=0.01c)，3種不同的倒圓半徑分別計為Case-R0、Case-R0.5和Case-R1，如圖2所示，圖中SS表示水翼吸力面，PS表示水翼壓力面。

圖2 不同翼端倒圓半徑示意圖Fig.2 Schematics of tip rounding with different rounding radius

為便于湍流統計分析，采用“〈 〉”和“ ′ ”分別代表時間平均統計量和瞬時脈動量。瞬時速度u、v、w分別為x(流向)、y(法向)和z(展向)3個方向的速度分量，u=〈u〉+u′,v=〈v〉+v′,w=〈w〉+w′。

如圖3所示，定義垂直于翼弦方向的速度為間隙泄漏流速Utip，公式為

圖3 流場的統計平面和速度坐標定義Fig.3 Statistical planes for investigating flow field and definitions of velocity coordinate

Utip=〈u〉sinα-〈v〉cosα

(1)

1.2 湍流模型

為準確解析翼端間隙區流場結構，采用可直接求解部分湍流脈動特征的VLES方法[16]進行流場計算。該方法是介于RANS(雷諾平均方程)方法和DNS(直接數值模擬)方法的混合方法，相比RANS方法可以獲取流場中更為豐富的湍流信息，近年來在復雜工程湍流的模擬中得到成功應用[17]。

VLES方法的基本思想是通過引入湍流尺度求解控制函數Fr，實現流場中部分湍流的直接求解，降低對湍流進行模化的比例，即減小模化的雷諾應力，公式為

(2)

控制函數Fr將實現VLES從RANS到DNS間的過渡，當網格分辨率足以求解全部湍流時，即Fr趨近于0，VLES類似于DNS；在粗網格或無限雷諾數條件下，Fr趨近于1，VLES等效于RANS方法；Fr介于0與1之間時，VLES則類似于LES的亞格子尺度模型(SGS)，可以直接求解大尺度湍流結構。本文采用文獻[18]提出的控制函數Fr，即

(3)

其中

Lc=Cx(ΔxΔyΔz)1/3

Li=k3/2/εLk=ν0.75/ε0.25

式中Lc——網格長度尺度

Li——湍流積分長度尺度

Lk——Kolmogorov長度尺度

k——湍動能β——系數

ε——湍流耗散率

ν——流體運動粘度

β=2.0×10-3，n=2，Cx=0.61。

由于RANS方法的SSTk-ω模型考慮了逆壓邊界層中湍流剪切應力傳輸效應，可以有效預測逆壓梯度條件下的流動分離，在復雜工程湍流求解中應用廣泛[19-20]，基于該模型的VLES方法也被成功用于軸流壓氣機葉端間隙泄漏渦的求解[17]，獲得了豐富的間隙區湍流結構。因此，本文以SSTk-ω模型為基礎構造VLES湍流模型。根據Boussinesq假定，對雷諾應力張量的調整可以近似等效為對湍流粘度的調整。因此VLES模型的控制方程與SSTk-ω模型相同，僅需對SSTk-ω模型的渦粘系數進行修正，即

(4)

ρ——流體密度

ω——湍流比耗散率

S——應變率不變量

F2——SSTk-ω模型混合函數[21]

系數a1=0.31。

1.3 計算域網格劃分

采用六面體結構化網格離散計算域，如圖4a所示，共劃分網格總數分別為2.1×106、4.6×106和1.01×107的3套網格進行網格無關性驗證(詳見2.1節)。網格總數為1.01×107的方案，在τ=0.1c和τ=0.02c兩種間隙下，間隙區域內網格層數分別為65層和55層，近壁面第1層網格對應y+<1.5，τ=0.1c間隙下x/c=0位置的斷面內網格分布如圖4b所示。

圖4 網格分布(τ=0.1c)Fig.4 Computational mesh

1.4 邊界條件和求解方案

采用CEL語言將VLES湍流模型集成到ANSYS CFX軟件中進行流場計算。計算域入口給定平均流速u∞，湍流強度設為2%，出口設置靜壓；計算域側壁和水翼表面均設為無滑移壁面。對流項離散采用高階精度離散格式，非定常求解的時間推進采用二階后向歐拉方法，時間步長t=1.5×10-3c/u∞，不同網格對應的庫朗數最大值均小于2。

2 結果與討論

2.1 計算模型驗證

選取間隙寬度τ=0.1c無翼端倒圓條件下的計算結果與文獻[15]的實驗結果進行對比，驗證計算模型的合理性和計算結果的可靠性。圖5為3種網格所得水翼下游x/c=1位置處TLV對應的時均流向渦量〈ωx〉等值線。從圖中可以看到，不同網格計算所得翼端間隙渦范圍差異并不明顯，〈ωx〉c/u∞=4的等值線均近似呈圓形分布，但與實驗結果相比，間隙渦范圍略小且略向下偏移。

圖5 時均流向渦量〈ωx〉等值線(x/c=1, 〈ωx〉c/u∞=4)Fig.5 Iso-line of non-dimensional time-averaged vorticity 〈ωx〉c/u∞=4 at x/c=1

在VLES等可直接求解大尺度湍流結構的湍流模型求解的流場中，總湍動能ktotal由直接求解的湍動能kres(kres=0.5(〈u′u′〉+〈v′v′〉+〈w′w′〉))和模化的湍動能kunres兩部分組成。kunres表示湍流模型模化的小尺度湍流結構包含的湍動能，kunres在總湍動能中的占比越小，表明流場中越多的大尺度湍流結構被直接求解，計算結果越準確。圖6和圖7分別展示了不同網格數量下流向x/c=-0.1位置處模化的湍動能占比云圖和間隙區湍動能沿流向的分布，隨網格精細度的增加，模化的湍動能占比kunres/ktotal由0.5減至約0.2，可見網格總數為1.01×107的網格方案直接求解了間隙區大部分湍流結構，后續以該網格方案的計算結果進行間隙湍流特性及流動損失分析。

圖6 不同網格數量下模化的湍動能占比Fig.6 Ratio of modeled to total kinetic energy for different mesh resolutions

圖7 不同網格數量下翼端間隙區kunres/ktotal沿流向的分布Fig.7 Distribution of kunres/ktotal along flow direction in tip clearance region for different mesh resolutions

2.2 翼端間隙區流場結構

圖8 時均流場的間隙渦結構及直接求解的湍動能分布Fig.8 Structures of time-averaged tip clearance vortices and contour plots of resolved turbulent kinetic energy

翼端邊緣倒圓半徑對TSV的形態和范圍有明顯影響。間隙寬度τ=0.1c時，由于翼端倒圓減小了間隙流與水翼端面的夾角，導致水翼翼端壓力面邊緣處的分離流動速度梯度增加、TSV徑向尺度減小，翼端倒圓半徑r=0.5%c時翼端間隙內形成了多個分離渦結構(圖8b)，且沿流向向下游延伸，倒圓半徑r=1%c時壓力面邊緣的分離流動匯合為單一的TSV。間隙寬度τ=0.02c時，隨倒圓半徑增加，TLV尺度和范圍增大，TSV的影響范圍則呈下降趨勢。大間隙下，kres主要分布在TSV流動區域，且隨翼端倒圓半徑增加而減小；小間隙下，kres主要分布在TLV和IV流動區域內，隨翼端倒圓半徑增加而增大。

圖9展示了翼端間隙區y-z平面內的渦系簡化結構。在水翼端部壓力面與吸力面壓差作用下，水翼壓力面邊緣產生的流動分離，在水翼頂端發展形成一個或多個TSV，翼端間隙流流經水翼吸力面邊緣時形成強剪切流動，在水翼吸力側發展為TLV。小間隙下，TLV對端壁邊界層的強卷吸作用，誘導近端壁區流動反向旋轉，形成誘導渦(IV)。

圖9 翼端間隙區渦系結構示意圖(y-z平面)Fig.9 Schematic of flow structures in tip clearance region (y-z plane)

2.3 間隙渦湍動能和湍流損失

圖10為不同翼端間隙條件下間隙區湍動能沿流向的分布，不同位置處模化的湍動能占總湍動能的比例kunres/ktotal均小于0.20，表明本文采用的VLES方法和求解模型可以直接求解翼端間隙區附近的大部分湍流結構。

圖10 不同間隙條件下的翼端間隙區kunres/ktotal分布Fig.10 Distribution of kunres/ktotal in tip clearance region for different tip clearance schemes

已有研究表明，間隙湍流損失與間隙渦旋結構的形成和消散導致的湍流混合直接相關，文獻[4]認為湍流損失與湍動能成正比，湍流損失系數ξ計算公式為

(5)

式中A——y-z平面內渦結構所占面積

由圖9可知，IV旋轉方向與TSV和TLV的旋轉方向相反，TSV集中在間隙內，TLV主要分布于水翼吸力面側且展向(z向)低于水翼端面，因此在式(5)的計算中，翼端間隙內流向渦量〈ωx〉小于0的區域定義為IV區，間隙內〈ωx〉大于0的區域定義為TSV區，展向坐標小于水翼頂端所在位置且〈ωx〉大于0的區域定義為TLV區。

圖11展示了不同間隙條件下，翼端間隙區湍流損失系數ξ沿流向的分布及不同渦結構對應的湍流損失占比。由圖11可知，兩種間隙尺寸下間隙湍流總損失沿流向均呈增加趨勢。大間隙下(τ=0.1c)，湍流損失隨翼端倒圓半徑增加而減小，對比間隙湍流損失的組成可以發現，其主要由TSV產生，但隨倒圓半徑增加，TSV產生的湍流損失占比顯著降低，無翼端倒圓(Case-R0)時TSV產生的湍流損失最小，占比約為70%，而倒圓半徑r=1%c(Case-R1)時TSV產生的湍流損失占比僅約為30%，TLV對應的湍流損失占比則隨倒圓半徑增加明顯增大。小間隙下(τ=0.02c)，間隙區湍流損失主要由TLV和IV產生，間隙區不同渦結構產生的總湍流損失隨翼端倒圓半徑增大均略有增加，間隙內TSV產生的湍流損失占比隨翼端倒圓半徑增加而減小，無倒圓時其占比約為10%，倒圓半徑r=1%c時其占比約為5%。

圖11 間隙湍流損失分布及不同漩渦對應的損失占比Fig.11 Tip clearance turbulent loss and loss proportions of different tip clearance vortices

2.4 間隙區湍流雷諾應力

為進一步分析翼端間隙流動湍流特性，圖12顯示了不同間隙形態下流向x/c=0.1位置處的雷諾應力分布云圖和流向渦量〈ωx〉等值線，圖中標記了TLV、TSV和IV等間隙渦旋結構。從圖中可以看到，不同翼端形狀下翼端間隙流動的雷諾正應力整體均大于雷諾切應力，且以法向正應力〈v′v′〉和展向正應力〈w′w′〉為主，雷諾應力分布區域與間隙渦系分布區域趨于一致[2]。在大間隙下(圖12a)，雷諾正應力主要產生于TSV區域，由于翼端倒圓減小了水翼壓力面邊緣間隙流動的沖角[22]，抑制了TSV的形成和發展，因此雷諾應力隨著倒圓半徑增加而減小。小間隙下(圖12b)，水翼吸力面邊界層內流向正應力〈u′u′〉占主導，法向正應力〈v′v′〉則主要集中于IV區域，而在TLV區域，展向正應力〈w′w′〉顯著大于其他應力分量；小間隙下翼端倒圓導致間隙泄漏速度Utip(見式(1))和泄漏流量顯著增加，致使間隙流與端壁邊界層的剪切作用增強、圍繞TLV中心的速度梯度增加，因此雷諾應力隨翼端倒圓半徑增加而增大。

圖12 翼端間隙區湍流雷諾應力分布(x/c=0.1)Fig.12 Contour plots of resolved Reynolds stresses in tip clearance region in a y-z plane (x/c=0.1)

2.5 間隙渦系湍流損失機理

翼端間隙渦系增加了水翼間隙區附近的湍動能(圖8)，湍動能的生成與速度梯度和雷諾應力相關[23]，其控制方程為

(6)

其中

式中Pk——湍動能生成項，表示大尺度湍流渦通過速度脈動與平均流動相互作用而獲取能量

Dk——擴散項，表征速度脈動、壓力脈動和粘性等引起的能量擴散

p′——脈動壓力

角標i和j用于表示方向(x、y或z)。

了解翼端間隙渦的湍動能生產機制，有助于揭示翼端流動的湍流損失機理。由于間隙湍流損失與湍動能成正比(式(5))，為分析不同翼端間隙形態下湍動能生成及湍流損失的差異，將湍動能生成項Pk分解為

Pk=P1+P2+P3+P4+P5+P6

(7)

其中

圖13為流向x/c=0.1位置處y-z平面內P1～P6的分布。大間隙下(圖13a)，水翼端部間隙內TSV的湍動能生成以P3和P6兩項為主，且P3和P6均隨翼端倒圓半徑增加而減小；對比圖11可知，翼端間隙內的雷諾正應力顯著大于雷諾切應力，因此法向速度的展向梯度?〈v〉/?z和展向速度的法向梯度?〈w〉/?y是湍動能生成項P6的主要來源。

小間隙下(圖13b)，IV區域的湍動能生成項以P2占主導，在TLV區域湍動能生成項則以P6為主，P2和P6均隨翼端倒圓半徑增加而增大；P1在靠近水翼吸力面存在明顯的負值區域，對比圖12的流向雷諾應力〈u′u′〉，表明靠近翼端間隙的水翼吸力面存在明顯的流向速度梯度?〈u〉/?x；而P4在水翼吸力面附近為正值，由于切應力〈u′v′〉相對較小，因此靠近水翼吸力面的區域內具有較大的流向速度的法向梯度?〈u〉/?y和法向速度的流向梯度?〈v〉/?x。

圖13 湍動能生成項Pk不同分量的分布圖(x/c=0.1)Fig.13 Contour plots of turbulent-kinetic-energy production components in Pk at x/c=0.1

為了定量對比雷諾應力和速度梯度對湍動能生成項的影響，在圖13中，選取line1、line2和line3所示3個特征位置，分析τ=0.1c的TSV流動區域以及τ=0.02c的IV和TLV流動區域的雷諾應力和速度梯度變化規律。line1～line3的位置坐標如表1所示。

表1 特征線位置坐標Tab.1 Characteristic line position coordinates

圖14為大間隙(τ=0.1c)下不同間隙形態的line1位置處湍動能生成項P3和P6及對應的雷諾應力和速度梯度分布。由圖可知，P3和P6均隨翼端倒圓半徑增加而減小，雷諾應力〈w′w′〉和〈v′w′〉亦隨倒圓半徑增加而減小，但不同方向的速度梯度隨翼端倒圓半徑增加整體呈增大的趨勢。因此，大間隙下TSV流動區域內湍動能隨倒圓半徑增加而減小(圖11)主要是由雷諾應力〈w′w′〉和〈v′w′〉的減小導致的。對比不同方向的速度梯度還可以發現(圖14c、14g、14h)，法向速度的展向梯度?〈v〉/?z顯著大于其他速度梯度分量，表明大間隙下TSV的湍動能生成由?〈v〉/?z主導。

圖14 間隙分離渦區域特征位置的雷諾應力與速度梯度(τ=0.1c)Fig.14 Reynolds stresses and velocity gradients in a characteristic position of tip separation vortices (τ=0.1c)

圖15為小間隙(τ=0.02c)時不同翼端形狀下line2和line3位置處湍動能生成項、雷諾應力和速度梯度的分布。在IV區域(line2)，雷諾應力〈v′v′〉和速度梯度?〈v〉/?y的量值均隨翼端倒圓半徑增加呈增大趨勢，對應的湍動能生成項分量P2隨倒圓半徑增加而增大。在TLV區域(line3)，P6隨翼端倒圓半徑增加而增大，其原因是與泄漏渦對應的雷諾應力〈v′w′〉和速度梯度(?〈v〉/?z+?〈w〉/?y)均隨倒圓半徑增加而增大，倒圓半徑r=0.01c條件下增幅尤為明顯。

圖15 誘導渦和間隙泄漏渦特征位置的雷諾應力與速度梯度(τ=0.02c)Fig.15 Reynolds stresses and velocity gradients in characteristic positions of induced vortex and tip separation vortex (τ=0.02c)

3 結論

(1)大翼端間隙下，翼端水翼壓力面邊緣產生的流動分離在水翼頂端發展形成包含一個或多個渦旋結構的間隙分離渦(TSV)，間隙流在水翼吸力面邊緣產生的強剪切流動在水翼吸力側發展為間隙泄漏渦(TLV)；隨翼端間隙寬度減小，間隙泄漏渦對端壁邊界層的卷吸作用增強，誘導端壁區流動反向旋轉形成誘導渦(IV)。

(2)翼端間隙附近的雷諾應力分布與間隙渦系的分布趨于一致，以法向正應力〈v′v′〉和展向正應力〈w′w′〉為主。大間隙下，雷諾應力集中分布于TSV區域，翼端倒圓半徑的增加將減小水翼壓力面邊緣間隙流動與水翼端面的夾角，進而抑制TSV的形成和發展，減小間隙內的雷諾應力。小間隙下，〈v′v′〉和〈w′w′〉分別集中于IV和TLV區域，翼端倒圓將顯著增加泄漏流量，致使間隙流與端壁之間的速度梯度和圍繞TLV中心的速度梯度均大幅增加，雷諾應力也隨翼端倒圓半徑增加而增大。

(3)間隙區湍流損失與湍動能成正比，不同翼端形狀下的間隙湍流損失機理存在差異。大間隙下，TSV是產生湍流損失的主要區域，其湍動能控制方程的湍動能生成項以P3和P6分量為主，速度梯度?〈v〉/?z和雷諾正應力〈w′w′〉分別為影響P3和P6的主要因素；隨翼端倒圓半徑增加，間隙湍流總損失因間隙區〈w′w′〉和〈v′w′〉的顯著減小而減小。小間隙下，間隙區湍流損失主要由TLV和IV產生，IV區域內P2是影響湍動能生成的主要因素，隨翼端倒圓半徑增加，其雷諾正應力〈v′v′〉和速度梯度?〈v〉/?y

增大，對應的湍流損失增加；TLV區域內的湍動能生產主要受P6影響，其主導因素雷諾應力〈v′w′〉和速度梯度(?〈v〉/?z+?〈w〉/?y)均隨倒圓半徑增加而增大，導致湍流損失增加。