基于SST-DDES方法的孤立氣膜孔流動研究

王 鵬

(中國航發沈陽發動機研究所，沈陽110015)

1 引言

航空發動機氣冷渦輪中的冷氣/主流摻混過程直接影響渦輪的流場結構、流動損失和冷卻效果。其中的復雜機理及規律已成為先進氣冷渦輪設計中的一個重要問題和技術瓶頸，同時也是進一步深入挖掘氣冷渦輪效率潛能的重要研究方向。

目前，研究手段和方法是制約冷氣摻混研究的關鍵因素。對于真實復雜氣動熱力環境下的氣膜冷卻研究和設計優化，參數繁多且相互制約，相應流動機理的揭示還需要有詳細可靠的流場數據支撐。用于設計優化結果驗證的實驗手段常常受到諸多限制，而數值手段的適用性則主要取決于其精度，目前制約冷氣摻混數值計算精度的主要因素之一是湍流的預測。氣膜冷卻的各向異性湍流流場中存在著復雜的時空多尺度渦系結構，數值模擬必須考慮到所有相干結構的綜合作用才能更精確地預測摻混流動，而這恰恰是RANS(雷諾平均)方法的力不能及之處。LES(大渦模擬)雖然可以較好地解決這一問題，但當雷諾數較大或幾何較復雜時，對計算資源要求很高，且對壁面附近流動的模擬存在固有困難，要在氣冷渦輪的系統性研究中應用存在很大困難。

在這種情況下，DES[1](分離渦模擬)類RANS/LES混合計算方法受到了廣泛重視。其通過對流場特征的判斷，在以耗散為主的區域采用RANS，以大渦輸運為主的區域采用LES，并實現兩者的自動過渡和切換，兼顧了計算的精度和效率[2]。然而DES類方法自建立以來主要用于外流中的大尺度分離及噪聲等問題[3-6]，對內流中小尺度分離尤其是渦輪葉片氣膜孔附近復雜流動的DES類研究相對較少[7-11]，且研究中采用的大都為DES97等原始的DES類方法。這些方法易產生雷諾應力模化不足即MSD(模型雷諾應力損耗)問題，進而帶來網格誘導分離等非物理問題，而近年來改進的DES類方法如延遲的DES方法(DDES)等還少有應用。為此，本文采用基于SST湍流模型的DDES方法，從渦輪葉片氣膜冷卻孔附近流動的本質即橫流射流問題入手，開展射流與主流摻混過程中的流動機理和損失機理研究。

2 數值方法

湍流模擬中采用基于SST湍流模型的DDES方法[12-13]。該DDES方法的構造基礎是Strelets[3]將DES的基本思想與SST湍流模型相結合所得到的SST-DES方法，即將SST模型的湍動能輸運方程耗散項改寫為與網格尺度Δ和湍流長度尺度Lt相關的DES形式：

該SST-DES方法與早期的SA-DES等方法一樣，存在著MSD問題和相應的網格誘導分離問題。為減小此類風險，Menter等基于SST模型的第一混合函數F1和第二混合函數F2對其進行改進：

式中：系數CDES取0.61。本文取FSST=F2，以最為嚴格地避免發生網格誘導分離。圖1為改進前后FDES函數中第一項隨混合函數變化特性對比。可見，改進后原SST-DES方法具有了延遲特征，成為一種DDES方法。

在基于有限體積法的數值計算中，空間離散格式的精度是難點之一。對于DES類方法，由于在不同流場區域同時存在著RANS和LES模式，若全場采用中心型格式則人工粘性系數將引入較多不確定性，若全場采用迎風型格式則耗散較大進而抑制小尺度結構的生成和發展。因此，引入了具有加權混合思想的混合型空間離散方法，在不同流動區域采用不同的空間離散格式[14-15]。對于編號為ip的空間某點處的通量Φip可表示為：

式中：σ為與流場有關的混合函數，Φip,U為中心型格式計算的通量，Φip,C為迎風型格式計算的通量。

該混合型空間離散方法具有以下特點：在LES模式流動區域趨近于0，表現為低耗散的中心型格式，利于解析小尺度流動結構；在壁面及遠場無旋區域趨近于1，表現為耗散較大的迎風型格式，利于抑制數值振蕩[15]。為進一步避免DES計算中發生網格誘導分離風險，對σ作如下附加限制：

式中：BF1和BF2為混合因子，本文選擇第二混合函數，即BF1=0，BF2=1。為避免在中高CFL數情況下中心型格式帶來數值振蕩，對σ又施加了基于CFL數的限制器：

式中：CFLmax、CFLEXP 分別取5.0和1.0。

3 物理模型和計算設置

Ajersch等[16]對射流雷諾數Rejet為4 700的方孔橫流射流開展了大量實驗研究，獲取了0.5、1.0和1.5三種吹風比[17]下的平均流參數和湍流脈動參數在流場不同位置處的空間分布，可用于數值模擬方法在此類流動中適用性的校驗。本文以其物理模型為計算模型，采用SST-DDES方法研究了真實跨聲氣冷渦輪氣膜冷卻中較典型的0.5吹風比狀態，同時采用相同的物理模型開展了基于SST湍流模型的RANS和URANS計算，以進行對比分析。

實驗中有6個均布方孔，其排列方向和射流方向均與橫流流動方向垂直，可實現多孔流動的周期性。數值模擬中選取一個周期作為計算域，圖2為具體的計算幾何模型。定義流向為x方向，展向為y方向，法向為z方向，射流孔中心為坐標原點。令射流孔邊長D為特征長度，與實驗模型一致取為12.7 mm，計算域流向長度共為51.0D，法向為25.0D，展向為3.0D。射流孔位置為 x=(-0.5～0.5)D，y=(-0.5～0.5)D，射流孔入口位置為 z=-5.0D，主流入口位置為x=-10.5D，出口位置為x=40.5D。

主流進口邊界條件為給定總溫、速度分布和湍流度分布，具體分布形式按照實驗中的實測結果給出[18]，見圖3；出口邊界給定大氣壓作為背壓。主流端壁為無滑移壁面，兩側面為對稱邊界，上邊界為滑移壁面。射流進口給定總溫、速度和中等湍流度，側壁為無滑移壁面。主流和射流進口總溫均為293.15 K。

計算域網格總數為300萬，分27塊，全場均為結構化網格。圖4為計算域網格分布示意圖，各處壁面的 y+均保證在1以下。DDES計算和URANS計算中的物理時間步長取0.000 15 s，虛擬時間步數取10步。非定常時間離散采用二階歐拉后差，空間離散以準二階迎風格式為基礎，能量方程求解中采用熱能方程形式。

4 結果分析

4.1 時均流場分析

圖5為 y/D=0截面上x/D=0、1.0、8.0位置處流向速度U沿法向的分布，其中DDES和URANS的結果為時間平均值。速度采用射流出口速度Vjet進行無量綱處理，法向距離采用射流孔邊長進行無量綱處理。由于復雜流動集中在壁面附近，因此主要針對法向約4.0D范圍內的結果進行分析。

在x/D=0即射流孔中心位置，不同層次的數值模擬方法所預測的壁面附近流向速度分布幾乎無差異，且都與實驗結果吻合較好。這主要是由于該位置處于射流與主流相互作用的起始階段，平均流的壓力梯度和速度梯度占主導，流動的非定常性和摻混程度都較弱，RANS就可以實現較好的預測。從分布趨勢看，射流對主流的影響主要集中在法向2.0D以下。值得注意的是，在法向2.0D以上區域，計算值大于實驗值。這主要是由于考慮到計算資源的限制，對非研究重點但實際中存在的主流上壁面邊界層不予模擬，進而使得主流的流向速度在該位置偏高，但這并不影響對下壁面附近流動特征的研究。在法向(0.2～0.4)D范圍內，流向速度沿法向的梯度突然減小，實驗由于測量點有限未能捕捉到該現象，該現象對應主流與射流之間的強剪切過程。

在x/D=1.0位置，流動的非定常性和摻混程度加強，不同數值模擬方法預測結果間的差異開始顯現。相比而言，DDES方法預測結果與實驗結果吻合得更好。主要差異集中在法向0.6D以下范圍內，URANS和RANS方法均過高地預測了射流尾跡中回流流動的程度，特別是RANS方法。在法向約0.8D位置，也存在一個流向速度法向梯度突變的現象，這與射流內側邊緣的剪切流動相對應。

在x/D=8.0位置，經過一定距離的流動摻混后，DDES的預測精度明顯更高，主要體現在法向1.2D以下射流主體區域，URANS和RANS預測的流向速度均明顯偏小。

圖 6為 y/D=-0.5截面上 x/D=0、3.0、5.0位置處展向速度W的分布。在x/D=0位置，各方法的預測結果仍基本無差異，且均與實驗結果吻合較好。該位置處于射流孔側面邊緣，射流從孔內噴出后，孔內的邊界層突然變為自由剪切層，在K-H不穩定性的作用下，剪切層失穩進而卷起一對反向旋轉的旋渦。圖中展向速度在近壁面附近的速度峰即對應該旋渦的起始發展過程。

在x/D=3.0和x/D=5.0位置，反旋渦對進一步發展，強度減弱且相互融合，DDES預測結果與實驗結果吻合得更好，URANS和RANS均過高預測了反旋渦對上半區的強度。在反旋渦對下半區，各方法均預測出了近壁面附近的速度峰，表征了反旋渦對的強烈下洗過程，但實驗結果受測量手段影響未能捕捉到該速度峰。值得注意的是，在x/D=3.0位置，實驗結果在法向0.8D附近存在兩個較弱的速度峰，URANS和RANS均未捕捉到該特征，而DDES實現了對該特征的捕捉。整體而言，DDES在展向速度方面的預測精度相對較高，但在局部也與實驗結果存在一定差異。這除了表明DDES方法自身還需進一步改進外，也與計算域受計算資源限制僅截取了單個流動周期有關。

圖7為y/D=-1.0截面x/D=0位置處法向速度V的分布。在射流影響下，展向1.0D位置處的主流也具有顯著的法向速度，在法向1.0D附近達到峰值；DDES的預測結果同樣與實驗結果吻合得較好，其他兩種方法的預測結果均偏小。

圖8為y/D=0截面x/D=8.0位置處無量綱湍動能的分布。由前文可知，在法向1.2D附近區域，?U/?z較大，URANS和RANS高估了該區域的湍動能峰值，DDES則較好地預測了湍動能的峰值水平。但應注意，實驗測量結果中在壁面附近還存在另一個局部湍動能峰值，而各方法計算結果均未捕捉到該峰值，這和計算預測的速度場等與實驗結果存在一定差異有關。

圖9為y/D=0截面x/D=3.0、8.0位置處雷諾切應力的分布。在法向1.0D附近即射流迎風面與主流發生強剪切的區域，URANS和RANS捕捉的局部應力峰值顯著過高，這與圖8中計算渦粘系數的湍動能局部峰值預測過高相對應，而DDES則實現了相對更為合理的捕捉。對于該現象，Andreopoulos等[19]的研究認為，RANS方法的渦粘模型在與?U/?z正負相反的流動區域理應具有一定的預測能力，但當存在流線匯聚或發散以及?W/?y的影響較顯著等情況時其預測能力較弱，該區域附近的流動恰好由于主流繞射流的流動而存在這些情況。

圖10進一步給出了DDES方法預測的射流孔出口法向速度(以Vjet進行無量綱)等值線與實驗結果對比。在所研究的0.5吹風比狀態下，射流孔出口流動受主流強烈影響，從射流迎風側到背風側法向速度逐漸增加，最高值約1.7Vjet，上游半區的梯度相對較大。圖11為DDES方法預測的射流孔出口平面流線分布與實驗結果中速度矢量分布的對比，絕對速度沿流向逐漸增加反映了射流在流向的偏轉程度，同時流動也向兩側偏轉，且其偏轉程度隨距射流孔中心距離的增加而增加。從圖10和圖11的對比結果看，DDES方法預測的主要特征均與實驗結果吻合較好。

圖12為y/D=0截面射流孔附近DDES方法預測的時均流場流線分布。在主流作用下，射流還未到射流孔出口便開始發生彎曲和流管收縮，流速相應增加，射流主體逐漸偏轉，最終與主流方向趨于平行，在射流迎風側和背風側均存在回流區。由于吹風比較小，隨著向下游的發展，射流速度的法向分量迅速減小，即射流的穿透能力較弱。

圖13為射流孔背風側邊緣即x/D=0.5截面DDES方法預測的時均流線分布。在y/D=±0.5的射流孔兩側區域卷起了兩個旋向相反且具有良好對稱性的旋渦，其核心區域的流速高于主流流速，且該反旋渦對在該位置上游便已形成并逐漸變大。圖14為x/D=2.0截面時均流場流線分布。反旋渦對離開射流孔區域后向下游繼續發展，尺寸不斷增加，相互之間的影響逐漸加強，各自的強度則逐漸減弱，旋渦核心區的流速已顯著低于主流區。

圖15為射流孔附近主流壁面的極限流線分布。射流的存在對主流流動帶來阻礙，主流繞射流主體流動，可近似看作方柱擾流流動，在前緣存在馬蹄渦系，在擾流后方存在回流區。但射流不同于固體，自身可以變形且帶來自由剪切層，使得該區域渦系結構更加復雜。

4.2 瞬態流場分析

針對某時刻瞬態流場中射流孔附近的流場結構進行了基于Q法則的旋渦識別，Q值取320 000，同時采用流向渦量對渦系結構進行渲染，見圖16。由圖可知，瞬態流場中的大渦擬序結構，主要包括射流前緣穩定的馬蹄渦、射流側壁產生的反旋渦對以及射流后逐漸形成并發展的發卡渦結構等。在本文較低吹風比狀態下未見高吹風比狀態下存在的從射流前緣開始隨射流向下游移動的反向渦頭[20]。從大渦結構上的渲染狀態，可清楚看到反旋渦對的反向旋轉特征及前緣馬蹄渦兩個分支的反向旋轉特征，下游的發卡渦無論是平行渦腿還是垂直渦腿也都具有兩側分支旋轉方向相反的特征，同一側的馬蹄渦分支和反旋渦對分支的旋轉方向也彼此相反。流向渦較強的區域主要集中在反旋渦對和發卡渦的平行渦腿等位置。

圖17進一步給出了采用展向渦量和法向渦量渲染的渦系結構。馬蹄渦和發卡渦的渦頭等位置具有相對較高的展向渦量，且方向相同。法向渦量較強的區域則主要集中在發卡渦的垂直渦腿等位置，且也具有兩側方向相反的特征。圖18為某瞬時x/D=1.0截面的熵分布，可看到流動損失主要來源于包含反旋渦對在內的射流背風側下游區域。

5 結論

采用基于SST湍流模型的DDES法，對典型吹風比下的孤立氣膜孔流動進行了數值模擬，并與實驗測量結果及RANS/URANS結果進行了詳細對比。分析了包括各速度分量、湍動能、雷諾應力和大渦擬序結構等在內的流場信息。研究結果表明：

(1)在吹風比0.5的條件下，橫流射流的入射能力較弱，進入主流后迅速偏轉，對下游壁面形成覆蓋，其影響區域主要集中在法向2.0D范圍內；同時，射流還存在一定的展向擴展，但在方孔邊長為相鄰孔距離1/3的條件下，相鄰射流孔流動不存在明顯干涉。

(2)摻混流場中存在著射流前緣馬蹄渦、產生于側壁且具有強烈下洗特征的反旋渦對以及包絡其上的發卡渦等大渦擬序結構，且這些大渦結構直接影響平均流特征和氣膜覆蓋效果，同時也是摻混損失的主要來源。

(3)對于橫流射流類復雜摻混流動，與傳統的RANS/URANS方法相比，基于RANS/LES混合的SST-DDES方法不僅對時均流場的預測更合理、更準確，還可有效捕捉瞬態大渦擬序結構及其演化過程，更利于流動機理研究和精細化設計。

(4)與實驗結果相比，在射流影響區域，SST-DDES方法預測的時均結果仍存在流向速度虧損、展向擴張程度、湍動能和雷諾切應力偏大等不足，SST-DDES方法仍需適當改進以進一步提高對橫流射流流動的預測精度。

參考文獻：

[1]Spalart P R，Jou W，Strelets M，et al.Comments on the fea?sibility of LES for wings and on a hybrid RANS/LES ap?proach[J].Advances in DNS/LES，1997，1：4—8.

[2]孫明波，汪洪波，梁寒劍，等.復雜湍流流動的混合RANS/LES 方法研究[J].航空計算技術，2011，41(1)：24—33.

[3]Strelets M.Detached eddy simulation of massively separat?ed flows[R].AIAA 2001-0879，2001.

[4]Squires K D，Forsythe J R，Morton S A，et al.Progress on detached-eddy simulation of massively separated flows[R].AIAA 2002-1021，2002.

[5]Forsythe J R，Squires K D，Wurtzler K E，et al.De?tached-eddy simulation of fighter aircraft at high alpha[R].AIAA 2002-0591，2002.

[6]Rulik S，Dykas S，Wroblewski W.Modeling of aerodynam?ic noise using hybrid SAS and DES methods[R].ASME GT2010-22696，2010.

[7]Roy S，Kapadia S，Heidmann J D.Film cooling analysis us?ing DES turbulence model[R].ASME GT2003-38140，2003.

[8]Kim S I，Hassan I G，Zhang X Z.Unsteady simulation of film cooling flow from an inclined cylindrical jet[R].ASME HT2007-32401，2007.

[9]Funazaki K，Kawabata H，Takahashi D，et al.Experi?mental and numerical studies on leading edge film cooling performance:effects of hole exit shape and free-stream turbulence[R].ASME GT2012-68217，2012.

[10]梁俊宇，康 順.基于DES方法的平板孤立冷卻孔數值模擬研究[J].工程熱物理學報，2011，32(3)：395—398.

[11]翟麗娜.基于DES方法的氣膜冷卻數值模擬[D].北京：華北電力大學，2012.

[12]Menter F R，Kuntz M.Development and application of a zonal DES turbulence model for CFX-5[R].CFX Valida?tion Report CFX-VAL17/0503.

[13]Menter F，Kuntz M.A zonal SST-DES formulation[R].Russia：St.Petersburg Technical University，2003.

[14]Travin A，Shur M，Strelets M，et al.Physical and numeri?cal upgrades in the detached-eddy simulation of complex turbulent flows[J].Advances in LES of Complex Flows，2004，65：239—254.

[15]劉 健.采用DES類方法研究基本起落架非定常大范圍分離流動[D].北京：清華大學，2011.

[16]Ajersch P，Zhou J M，Ketler S，et al.Multiple-jets in a cross-flow:detailed measurements and numerical simula?tions[J].Journal of Turbomachinery，1997，119：330—342.

[17]Holdeman J D，Walker R E.Mixing of a row of jets with a confined cross-flow[J].AIAA Journal，1977，15(2)：243—249.

[18]Zhou J M.A computational and experimental investigation of gas turbine blade film cooling effectiveness[D].Univer?sity of British Columbia，1994.

[19]Andreopoulos J，Rodi W.Experimental investigation of jets in a cross-flow[J].Journal of Fluid Mechanics，1984，138：93—127.

[20]劉 斌，任麗蕓，葉 建，等.方孔橫向射流的大渦模擬研究[J].航空科學技術，2009，(2)：40—43.