王 鵬

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

1 引言

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

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

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

2 數(shù)值方法

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

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

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

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

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

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

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

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

3 物理模型和計算設置

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

實驗中有6個均布方孔，其排列方向和射流方向均與橫流流動方向垂直，可實現(xiàn)多孔流動的周期性。數(shù)值模擬中選取一個周期作為計算域，圖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。

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

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

4 結(jié)果分析

4.1 時均流場分析

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

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

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

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

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

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

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

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

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

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

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

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

圖15為射流孔附近主流壁面的極限流線分布。射流的存在對主流流動帶來阻礙，主流繞射流主體流動，可近似看作方柱擾流流動，在前緣存在馬蹄渦系，在擾流后方存在回流區(qū)。但射流不同于固體，自身可以變形且?guī)碜杂杉羟袑?，使得該區(qū)域渦系結(jié)構(gòu)更加復雜。

4.2 瞬態(tài)流場分析

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

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

5 結(jié)論

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

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

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

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

(4)與實驗結(jié)果相比，在射流影響區(qū)域，SST-DDES方法預測的時均結(jié)果仍存在流向速度虧損、展向擴張程度、湍動能和雷諾切應力偏大等不足，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].航空計算技術(shù)，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，F(xiàn)orsythe 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方法的平板孤立冷卻孔數(shù)值模擬研究[J].工程熱物理學報，2011，32(3)：395—398.

[11]翟麗娜.基于DES方法的氣膜冷卻數(shù)值模擬[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].航空科學技術(shù)，2009，(2)：40—43.