基于DDES 模擬的葉頂泄漏流與尾跡非定常干涉機理

李會，黃通，蘇欣榮，袁新

清華大學 熱科學與動力工程教育部重點實驗室，北京 100084

渦輪葉片是航空發動機的重要熱-功轉換部件。葉頂泄漏流是渦輪動葉葉片通道內部的關鍵流動結構之一，其對渦輪葉片的氣動損失、流動非定常性、高溫部件可靠性等有明顯影響。由于葉片兩側存在壓差，導致流體在壓力差驅動下從壓力面側經過葉頂間隙后流向吸力面側，進而形成泄漏流。在黏性作用下，泄漏流在吸力面附近與主流摻混并卷吸形成泄漏渦，其引起的流動損失可以占到整個通道損失的30%甚至更高［1］。Denton［1］認為葉頂泄漏流產生的損失主要由3 部分組成，其中泄漏流在通道內造成的摻混損失為主要的損失來源。因此，深入開展葉頂泄漏損失產生機制方面的研究對于減少泄漏流帶來的負面影響具有十分重要的理論意義。

葉頂泄漏流進入到主流通道后，會受到上游流動結構的影響從而產生非定常效應。Ameri等［2］采用數值方法研究了上游尾跡對下游葉頂間隙的傳熱影響，發現上游尾跡對下游葉片吸力面側的壓力分布產生的影響較小，但是會帶來傳熱的非定常特性。劉火星等［3］采用實驗方法研究了上游尾跡對葉柵通道內泄漏流的影響，發現上游尾跡使得下游通道渦區域損失減少，同時泄漏流損失也相應減少。Zhou K 和Zhou C［4］對比了雷諾時均Navier-Stokes（Reynolds Averaged Navier-Stokes， RANS）和非定常雷諾時均Navier-Stokes（Unsteady Reynolds Averaged Navier-Stokes， URANS）的結果，發現由于上游靜葉產生的周期性尾跡的影響，下游動葉葉頂泄漏流造成的損失反而減小；之后他們通過實驗改變上游尾跡達到了減少下游泄漏流損失的效果，進一步驗證了該結論［5］。楊佃亮和豐鎮平［6］通過施加非定常邊界條件的方法，發現上游尾跡會對葉頂吸力面前緣附近和尾緣附近流場產生強烈的影響。

除了上游流場引起的非定常流動，在葉片通道內部葉頂泄漏流也會失穩并誘發較強的流場脈動。Gao 等［7］發現葉頂泄漏流向下游輸運的過程中，在渦輪葉片吸力面側逆壓梯度的作用下，可能出現泄漏渦破碎現象，破碎后的流動呈現出明顯的非定常特征。Sell 等［8］通過數值研究發現，在葉片近間隙吸力面側后半部分存在回流區域，認為可能是葉頂泄漏渦破碎所致。Huang等［9］利用控制體模型方法研究了渦輪近葉頂負荷分布對泄漏渦動力學特性及損失的影響，研究發現，泄漏渦的穩定性受流向逆壓梯度影響較大，流向逆壓梯度足夠大時會誘發泄漏渦破碎，帶來額外的渦破碎損失。Yang 和Ma［10］使用數值模擬對一典型透平動葉泄漏渦的特征做了研究，結果表明逆壓梯度導致渦核附近的雷諾正應力增強，進而造成泄漏渦破碎。以上研究都說明了泄漏流在向下游流動過程中自身存在著非定常特性，普遍的觀點是由于逆壓梯度的存在導致了泄漏渦破碎，但對于是什么流動結構引起的逆壓梯度卻缺乏分析和說明。

為了準確地模擬泄漏流復雜的流動現象，理論上只有采用直接數值模擬（Direct Numerical Simulation， DNS）方法和大渦數值模擬 （Large Eddy Simulation， LES） 方法。受限于高雷諾數流動中的計算量問題，DNS 和LES 主要用于簡單幾何和較低雷諾數泄漏流的機制研究。Shang等［11］使用DNS 對扁平水翼做了仿真，將泄漏渦分成了3 個階段，即生成段、發展段、分裂破碎段。Gao 和Liu［12］采用了LES 方法對簡化的泄漏流模型進行研究，并且獲得了豐富的流場結構。目前高負荷透平通道內部流動分析主要采用RANS方法，其能夠在耗費較少計算資源的前提下，快速給出流場的時均結果，但是該方法無法實現物理上多尺度湍流的精細模擬［13］。對于葉柵通道內部的三維尾跡渦結構、泄漏流非定常流動等復雜流動結構，RANS 方法存在局限性，進而導致復雜流動的機制研究受限。混合RANS/LES 方法能夠以較小的代價實現對復雜幾何、高雷諾數流動的高精度湍流模擬，可以較好地解決計算精度和計算資源之間的矛盾。Spalart［14］在1997 年發展了脫體渦模擬（Detached Eddy Simulation，DES）方法，DES 方法是一種將LES 和RANS 相結合的方法。DES 方法在近壁面邊界層內采用非定常RANS 模型，在遠離壁面的強剪切區域自動切換到LES 方法。為了改進DES 方法對應力模化的能力，2006 年Spalart 等［15］將DES 方法升級為延遲脫體渦方法 （Delayed Detached Eddy Simulation， DDES）。DDES 方法既能很好地模擬主流的湍流特性，又比LES 和DNS 方法極大地減少了計算時間。目前關于葉輪機械的計算已有很多DDES 算例，Lin 等［16］采用DDES 方法對高壓透平級的非定常流動進行了分析；Bian等［17］對端壁流動的湍流特性分析也采用了DDES 方法；Wang 等［18］采用了DDES 方法對壓氣機的角渦分離和尾跡做了相關研究。然而目前還鮮有學者采用DDES 方法對葉頂泄漏渦和尾跡干涉的流動結構做精細模擬。對于泄漏渦破碎和與尾跡之間干涉的機制尚沒有明確認識。

針對現有研究在葉頂泄漏流與尾跡干涉以及泄漏渦破碎方面研究的不足，采用DDES 方法研究某高負荷透平動葉葉頂間隙泄漏流以及下游尾跡渦的非定常特性。首先從葉頂間隙泄漏渦的演變、泄漏渦與端壁渦和尾跡渦的相互作用現象進行分析，說明泄漏渦破碎的機制；然后具體分析泄漏流與尾跡干涉機制；最后系統地分析泄漏流與尾跡干涉帶來的流場損失特性，以期揭示其在二次流損失中的重要性。

1 數理模型

1.1 計算域與計算網格

所研究的計算模型為某高壓渦輪動葉直列葉柵，葉頂間隙高度為葉片高度的1%。計算域如圖1 所示（Ca為軸向弦長），進口馬赫數為0.38，基于軸向弦長的雷諾數為2.93×105。圖2為葉片表面和端壁上的計算域網格分布，還給出了前緣和尾緣處的網格細節。計算網格采用非結構網格，葉片表面上沿著流向方向有652 個網格點，間隙高度方向上布置了40 個網格點，壁面外第Ⅰ層網格點上的無量綱壁面距離y+≈1。基于RANS 計算結果和DDES 的初步計算結果對網格進行了多輪調整，主要加密了尾跡區和泄漏渦區域的網格點分布，以實現對小尺度流動結構的準確捕捉，最終使用的計算網格包含約1 400 萬個網格單元。

圖1 計算域示意圖Fig.1 Schematic diagram of computation domain

圖2 計算域網格劃分示意圖Fig.2 Schematic diagram of computation domain mesh

1.2 數值方法

DDES 計算采用自主開發的多塊并行有限體積代碼［19-22］，所求解的積分形式Navier-Stokes方程為

式中：τ0為虛擬時間；U為N-S 方程的守恒形式自變量；V為網格單元的體積；t為物理時間；Fc為對流通量；Fv為黏性通量；n為積分面積法向量；S為積分面積分別表示虛擬時間項、物理時間項和空間殘差項。

在計算中，采用雙時間步方法進行非定常計算，且物理時間項采用隱式四點二階方法進行離散。為盡可能減小數值耗散對流動結構模擬結果的影響，數值模擬中采用5 階WENO 格式和自適應低耗散格式［16］，進一步降低數值黏性、提高對小尺度流動結構的解析精度［23］。為了加快每個物理時間步內的收斂速度，程序中對虛擬時間項采用了隱式時間推進方法進行求解。為了進一步加速收斂，還采用了多重網格方法。在之前的研究工作［24］中已對數值方法進行了實驗驗證工作，驗證了該數值方法的有效性。

為了準確捕捉高頻、小尺度非定常流動結構，通過對泄漏渦區域非定常特性的分析，確定非定常計算中的物理時間步長為2.5×10-6s，以確保能準確計算高頻非定常流動結構。圖3 給出了計算過程中連續7 個物理時間步內的收斂情況。從圖3 中可見，每一個物理時間步約需內部迭代30 步后，殘差可降低3 個數量級。從圖4 尾跡區的監測點可以看出，非定常流場計算已經進入統計平均收斂階段，可以對非定常流動計算結果進行分析和統計平均。

圖3 非定常求解收斂歷史Fig.3 Convergence history of unsteady simulation

圖4 監測點壓力分布Fig.4 Pressure distribution at monitoring points

2 泄漏區非定常特性

在DDES 方法中，邊界層內使用RANS 模型，而在遠離壁面區采用LES 方法。除了基于y+等參數來判斷邊界層內部網格的密度外，還需要判斷主流區網格密度是否合適。采用湍動能模化率指標來判斷主流區的網格分辨率是否滿足LES 求解要求。湍動能模化率IQ 的定義為

式中：Eslov為通過直接求解解析得到的湍動能；Esgs為亞格子應力模化得到的湍動能。通常認為大渦模擬需要解析80%以上的湍動能，才能解析得到大部分湍流成分［25］，從圖5 中可以看出，主流區的湍流模化率均在80%以上，說明主流區基本解析到了大尺度湍流結構，滿足LES 對網格分辨率的要求。

圖5 湍動能模化率分布Fig.5 Contour of index of resolution quality

此外，還在中葉展的尾跡渦區域、葉頂間隙內的泄漏流分離泡和端壁渦區域分別設置了監測點，對監測點速度進行傅里葉變換并得到功率譜密度圖（Power Spectral Density， PSD），如圖6所示，可以看出DDES 計算結果在慣性子區很好地滿足了-5/3 次方率，這一結果也可以從側面驗證該數值方法的有效性。圖7 展示了泄漏渦、尾跡渦、尾跡與泄漏渦干涉區域最大幅值對應的頻率，其中U、V和W分別為3 個方向上的速度分量。可見泄漏渦與尾跡渦最大幅值對應的頻率相同，均為1 800 Hz，且尾跡渦的幅值小于泄漏渦的幅值。從圖7（c）中可以看出，在尾跡與泄漏渦相互作用區域，除了1 800 Hz 的頻率外，還出現了更高的倍頻，這是兩種渦系干涉后產生的破碎渦系的頻率，此時干涉區的幅值也有明顯增大，說明泄漏渦與尾跡渦之間具有很強的非定常干涉作用。

圖6 監測點速度功率譜密度分布Fig.6 Velocity power spectral density distributions at monitoring points

圖7 不同位置處最大幅值對應頻率Fig.7 Frequencies corresponding to maximum amplitude at different positions

2.1 葉片表面壓力系數分布

圖8 展示了50%和95%葉高截面上的壓力系數分布。圖8（a）中h為葉片高度；z為軸向位置；Cp為壓力系數。由圖8（a）中可以看出，50%葉高和95%葉高的壓力面側壓力系數分布基本重合，意味著泄漏流對葉片壓力面側的載荷分布幾乎不產生影響。圖8（b）中y為展向位置。從圖8（b）和圖8（c）中可以看出，葉片不同高度下壓力面側的云圖分布也基本一致。但是在不同葉高下，吸力面側的壓力系數分布在約10%～90%軸向弦長范圍內存在差異。在95%葉頂高度處，吸力側受到泄漏流的影響更為嚴重，因此產生了相較于50%葉頂高度處的壓力系數分布差異。順著流向方向，隨著泄漏流不斷增強，吸力面側的壓力系數差異也隨之逐漸增加。在60%軸向弦長以后，由于泄漏流逐漸遠離葉片吸力面側，其對吸力面側壓力的影響也逐漸減小。可見泄漏流對葉片表面載荷的分布具有較大的影響。

圖8 葉片表面和葉柵通道壓力分布Fig.8 Load distributions of blade surface and cascade passage

2.2 流動結構分析

圖9 為葉頂泄漏流時均流場結果，用Q準則表達渦結構，并用無量綱流向渦量系數進行染色。Q準則的定義是速度梯度的第二不變量，當Q大于0時表示流體中旋轉程度大于變形程度，認為該區域中存在渦。無量綱流向渦量系數的定義為

圖9 時均流場渦結構（Q=1×108 s-2）Fig.9 Time averaged flow field vortical structures （Q=1×108 s-2）

式中：Ω為渦量；u為速度；l為特征長度，此處選為軸向弦長；uin為進口速度大小。

為了方便展示流場中的渦結構，視圖選擇了從葉根朝向葉頂方向的視角。從圖9 中可以明顯地分辨出葉頂泄漏渦（TLV）、通道渦（PV）和端壁渦（EV）。在泄漏渦往下游發展過程中，端壁渦開始纏繞在泄漏渦周圍，使得泄漏渦出現不穩定現象，并使得泄漏渦的完整性被破壞。通道渦被泄漏渦擠到葉頂下方，然而其與泄漏渦互不干涉。在時均流場中，由于時間平均的作用，只有大尺度的渦結構保留下來，觀察不到小尺度的渦結構。為了探究葉頂區域的非定常特性，需要對瞬時流場做進一步分析。

選取了某一時刻的瞬時非定常流場進行分析，如圖10 所示。相較于時均流場，瞬時流場中包含更為豐富的小尺度渦結構。在端壁渦與泄漏渦相互作用之前，泄漏渦和通道渦以及端壁渦的結構與時均流場一致，這意味著這些渦結構在相互作用之前都保持穩定和近似定常的狀態。在瞬時流場中可以發現端壁渦和泄漏渦相互作用后泄漏渦被破壞，產生一系列誘導渦（IV）。尾跡渦（WV）從尾緣處脫落并向下游輸運，在輸運的過程中可以發現誘導渦與尾跡渦發生相互干涉。尾跡渦受到誘導渦的干涉后，尾跡渦也破碎成更小尺度的渦結構并往下游輸運。可見泄漏渦、端壁渦和尾跡渦之間相互作用使得葉頂區域的流場出現很強的非定常性。

圖10 瞬時流場渦結構（Q=1×108 s-2）Fig.10 Instantaneous flow field vortical structures （Q=1×108 s-2）

2.3 泄漏渦與尾跡干涉機制

為了更好地說明泄漏渦與端壁渦，以及其和尾跡渦間的作用過程，圖11 給出了尾跡渦脫落周期內（周期為τ）的4 個典型時刻。紅色虛線圓圈表示泄漏渦和端壁渦相互作用產生的誘導渦，從圖11（a）到圖11（d），該渦逐漸向下游輸運。紅色實線表示泄漏渦與尾跡相互干涉區尾跡渦的軌跡，紅色箭頭表示尾跡擺動的方向。在（1/4）τ時刻，泄漏渦和端壁渦相互作用后，誘導渦開始產生。此時從尾緣脫落的尾跡渦也開始向下游移動。到了（2/4）τ和（3/4）τ時刻，誘導渦繼續向下游傳播，并且泄漏渦和端壁渦相互作用產生的新的誘導渦也開始從上游往下游傳播。此時之前上游產生的誘導渦還沒有和尾跡渦發生干涉。到了（4/4）τ時刻，從紋影圖中可以看到誘導渦和尾跡渦相互作用，流場變得更加無序，并且渦結構的完整性遭到破壞。從不同時刻的紋影圖中還可以看出尾跡渦的運動軌跡也受到誘導渦系的影響，發生了明顯的左右擺動。

圖11 95%葉高處密度梯度紋影圖Fig.11 Contour of density gradient magnitude at 95% blade height

通過對非定常計算結果進行后處理，可以計算得到湍動能。圖12（a）為95%葉高截面上泄漏區與尾跡區的湍動能分布。從圖12（a）中可以看出：在此截面上，湍動能最大的區域不是尾跡區和泄漏流區域，而是泄漏流和尾跡渦發生相互干涉的區域。截取從尾緣到下游泄漏流與尾跡干涉區，并作周向平均得到圖12（b）。從圖12（b）中可以看出，從尾緣往下游移動的過程中，湍動能呈現先增加后減小的趨勢。這是因為隨著端壁渦與泄漏渦相互作用產生的擬序結構使得流場中的湍動能逐漸增加，在泄漏渦與尾跡渦發生干涉的區域湍動能達到了最大，隨著干涉后產生的渦破碎，湍動能逐漸耗散，導致湍動能逐漸減小。

圖12 湍動能分布Fig.12 Turbulent kinetic energy distribution

3 泄漏區非定常損失機制

從熱力學的角度來看，流場中的不可逆性是損失的來源，其可以分為兩類：由流體黏性引起的不可逆損失和由傳熱過程引起的不可逆損失。在目前常用的損失分析方法中，總壓損失系數、熵增等指標只能反映定常流場中總的損失情況，無法明確給出流場中損失生成的位置和強度，且其應用于非定常流場時常出現損失系數小于0 或者大于1 等非物理的情況。因此采用熵生成率對葉頂區域的時均和瞬時流場進行分析，非定常的熵輸運方程形式為［26］

式中：ρ為密度；s為熵；uj為速度分量；T為溫度；xj為某方向上的分量；τij為黏性應力張量分量；Sij為應變率張量分量；λ為導熱系數分別表示由流體黏性引起的熵增和由傳熱引起的熵增。在相關研究中［16，27］，發現由傳熱引起的不可逆性損失占比很小，損失主要為黏性引起的不可逆性損失，因此重點分析黏性不可逆性損失。黏性耗散引起的熵生成率為

與圖11 中4 個物理時刻相對應，圖13 給出了相應的黏性損失熵產率分布。從黏性耗散的云圖中可以看出，不同時刻下的黏性耗散除了發生在上游泄漏渦和端壁渦所處的位置，在下游泄漏渦和端壁渦相互作用的位置也出現了很強的黏性耗散。在（4/4）τ時刻，誘導渦和尾跡渦發生相互作用的區域，也可以看到黏性耗散增強。通過比對時均結果和瞬時結果的黏性耗散云圖可以看出，泄漏渦和端壁渦以及尾跡渦相互作用增強了流場的非定常性，并且帶來了黏性耗散損失。

圖13 不同時刻的黏性耗散損失分布Fig.13 Distributions of viscous dissipation loss at different times

圖14 選取了130%Ca截面，得到了3 個周期內的瞬時黏性耗散引起的熵生成率（紅色方塊）、時均流場黏性耗散引起的熵生成率（綠色實線）、時均流場總的熵產率（黑色實線）。從圖14 中可以看出，對于黏性耗散而言，瞬時流場下的結果大于時均流場下的結果，損失偏差最大為71%。這意味著泄漏渦和端壁渦以及尾跡渦相互作用導致的非定常損失不可忽略。從不同時刻的損失分布可以看出，損失出現周期性，這也說明了葉頂泄漏區的流動具有很強的非定常性。時均流場總的熵產率包括由黏性耗散引起的和由小尺度脈動引起的，定義小尺度脈動引起的熵產貢獻為時均流場下總的熵產率和黏性耗散的熵產率的差值，即

圖14 130% Ca軸向截面處的黏性耗散損失分布Fig.14 Distribution of viscous dissipation at 130% Ca axial section

由小尺度脈動引起的熵產如圖14 中粉色部分所示。可以發現，泄漏流與尾跡干涉區產生的小尺度脈動熵產貢獻占到時均流場總熵產的23.3%，占比較大，是不可忽略的一部分。這說明泄漏流與尾跡相互干涉產生的非定常損失是高壓透平葉頂泄漏流損失中的重要來源之一。因此對于葉頂泄漏流損失的研究，如果僅使用RANS 方法，僅能得到時均下的流場損失結果，對實際的損失估計不足，建議采用DDES 方法等可以得到瞬時流場且精度更高的方法。

4 結 論

使用DDES 方法對渦輪動葉葉頂間隙泄漏流進行數值模擬，并分析了泄漏渦與端壁渦、尾跡渦之間的干涉作用，最后分析了泄漏流區域的損失機制，得到如下結論。

1）DDES 方法可以精細地捕捉到泄漏渦和尾跡渦等渦結構，以及各種渦系相互作用產生的更小尺度的渦結構，對于分析葉頂泄漏流等復雜渦系之間的作用機制，可以很好地平衡計算代價和計算精度之間的矛盾。

2）葉頂泄漏渦在發展過程中會首先與端壁渦產生相互作用，泄漏渦出現不穩定現象，繼而發生破碎，并周期性地產生誘導渦；誘導渦在下游會與尾跡渦產生干涉作用，使得尾跡渦的軌跡出現周期性左右擺動，并破碎成更小尺度的渦系結構，并逐漸耗散成無序渦結構。渦系之間的相互作用具有很強的非定常性。

3）基于非定常熵輸運方程，得到了流場黏性耗散帶來的熵增，其中瞬時流場下的損失大于時均流場下的損失，最大偏差可達到71%，說明非定常流動帶來的損失不可忽略。由小尺度脈動所產生的損失，可以占到時均流場總損失的23.3%，是非常重要的損失來源之一。對于高效透平葉柵設計，需要考慮非定常流動帶來的損失，以進一步提高透平性能。