影響冷氣摻混數值模擬精度的若干問題分析

何宏疆，綦蕾，王鵬，劉火星

影響冷氣摻混數值模擬精度的若干問題分析

何宏疆，綦蕾，王鵬，劉火星

(北京航空航天大學能源與動力工程學院，航空發動機氣動熱力國家級重點實驗室，北京100191)

對冷氣摻混數值模擬中的計算網格、湍流模型、射流邊界條件等影響精度的若干問題進行分析。在吹風比0.5下，對不同網格分布、湍流模型、射流邊界條件進行數值模擬，所得結果與Ajersch的實驗結果進行對比。結果表明，出口下游的網格分布可適當稀疏以減少計算量；在所研究的幾種湍流模型中，k-ε模型所得結果與實驗結果吻合得最好；考慮射流通道內流動能提高精度，在射流出口給定流場分布也可保證計算結果的精度。

航空發動機；氣膜冷卻；方孔射流；冷氣摻混；數值模擬精度

numerical simulation precision

1 引言

隨著航空發動機效率的提高，渦輪部件進口溫度也不斷提高。而渦輪部件進口溫度的提高，直接威脅到發動機的安全性、可靠性和使用壽命。因此，探索成熟、高效的冷卻方法日益重要。自上世紀70年代以來，氣膜冷卻因其良好的冷卻效果被廣泛應用于航空發動機渦輪部件。冷氣摻混流動十分復雜，且有很強的三維特性，準確預測其特性有助于對冷氣再附位置、入射深度做出判斷。為此，很多學者對冷氣摻混問題進行了大量實驗和數值研究。

Zhou等[1，2]采用簡化的渦輪葉片氣膜冷卻模型，對吹風比0.5、1.0、1.5三種情況進行了實驗和數值模擬。結果表明，射流下游中線上的速度分布和雷諾應力分布相差較大，并認為采用的湍流模型均不能準確預測氣膜冷卻的流動特征。Hoda等[3]采用七種湍流模型來預測冷氣摻混問題，表明所有模型均高估了射流下游尾跡區域內的速度。Hassan等[4]對比了k-ε、RNGk-ε、realizablek-ε、k-ω四種兩方程湍流模型在解決冷氣摻混問題上的效果，并分析了其局限性，認為標準k-ε模型得到的模擬效果最好。Keimasi等[5]采用雷諾平均方法研究了三維湍流垂直入射主流，分別采用了標準k-ε模型加壁面函數模型和SST模型，最終得到的平均速度分布與實驗結果吻合得很好，但湍動能分布相差較大。Garg等[6]在文獻中展示了射流入口邊界條件對準確預測射流出口下游流動的重要性，表明不同射流入口邊界條件可導致射流出口下游高達60%的換熱系數差異。隨著計算能力的提高，大渦模擬(LES)方法也被逐漸應用到對這一問題的研究中來。Guo等[7]采用LES方法，得到了比雷諾平均方法吻合更好的結果。劉斌[8]采用LES方法對方孔射流問題進行了模擬，其計算模型與Ajersch的實驗模型相同，通過模擬吹風比1.5情況，加深了對各種渦結構發展過程的認識。

本文通過數值結果與實驗結果的對比，來研究計算網格分布、湍流模型、邊界條件等對冷氣摻混數值模擬精度的影響。

2 計算域與邊界條件

計算模型與Ajersch[2]的實驗模型相同。實驗中采用了6排垂直入射的方孔，本文選擇一個周期作為計算域。實驗研究了吹風比0.5、1.0和1.5三種情況，本文只針對吹風比0.5的情況。以實驗測得參數作為邊界條件，主流進口邊界條件為給定速度分布、湍流度分布和總溫。主流邊界層內的速度和湍流度分布見圖1，主流中各個方向的速度U=U∞=11 m/s，V=0，W=0，總溫293.15 K。射流進口邊界條件：射流進口湍流度為0.05，速度U=0，V=5.5 m/s，W=0，射流總溫293.15 K。壁面兩側采用對稱邊界條件，壁面采用絕熱無滑移壁面，出口給定大氣壓力，具體邊界條件參見文獻[1]、[2]。與圖2的坐標系相對應，定義X方向為流向，Y方向為法向，Z方向為展向。取射流孔邊長D為特征長度，且D=12.7 mm。主流流向長度為51D，法向長度為25D，展向寬度為3D。取射流孔中心為坐標原點，射流孔位置為X= -0.5D～0.5D，Z=-0.5D～0.5D，其入口位于Y=-5.0D處。主流入口位置為X=-10.5D，出口位于X=40.5D。

圖1 主流進口邊界層內的速度和湍流度分布Fig.1 Velocity and turbulence intensity distribution at mainstream inlet

圖2 計算域示意圖Fig.2 Computational domain

3 結果分析

將數值結果與Ajersch等所做實驗結果進行對比，包括流向速度、展向速度、法向速度、湍動能，且均以射流進口速度Vj進行無量綱化，距離以射流孔邊長進行無量綱化。

3.1網格分布對數值模擬精度的影響

對表1中前四種網格分別進行計算。各種網格在射流孔內均勻分布，射流孔上、下游流向網格的膨脹比分別為1.12和1.06，法向膨脹比為1.06，射流孔中法向膨脹比為1.20，應用SST湍流模型。四種網格下射流出口(Z/D=0)位置法向速度對比如圖3所示。可見，隨著網格數目的增加，不同網格的射流孔出口法向速度的差異逐漸減小。其中第1種網格過高給出了射流孔上游的入射速度，但隨著網格數目的增大，該位置的速度逐漸減小；網格3與網格4的結果相差很小。文獻[1]中，在距射流出口較遠的下游(X/D=3，5)位置，速度分量和湍動能與實驗結果均相差很大。網格3在X/D=10位置的流向尺寸為0.48D，在主流出口為1.66D，網格過于稀疏可能是導致下游流場精度過低的原因。對網格3進行加密，其中主流進口至射流出口下游X/D=1.5位置的網格分布不變，X/D=1.5至主流出口的流向進行加密，主流通道中展向和法向的網格分布不變，網格數目由361 854增加為573 032，加密后的網格稱為網格5。由圖4可見，網格加密后，主流中X/D=10位置的網格尺寸由0.48D縮小為0.17D，主流出口附近最大網格尺寸由1.66D縮小為0.45D。射流通道內網格分布保持不變。但網格加密前后，出口下游各個位置的流向速度分量完全重合。下文中如無特殊說明，均以網格3作為計算網格。

圖3 射流出口中線法向速度對比Fig.3 Comparison of vertical velocity at the center line of jet exit

圖4 網格3與網格5在Z/D=0位置上的流向速度對比Fig.4 Comparison of streamwise velocity forZ/D=0

3.2湍流模型對數值模擬精度的影響

對幾種常用湍流模型(層流模型、零方程模型、四種兩方程模型(k-ε、SST、k-ω、RNGk-ε)、Base?line湍流模型)進行驗證。圖5給出了射流孔中線不同流向位置(Z/D=0)的流向速度對比。在射流孔中心(X/D=0)位置，相對于其它湍流模型，層流模型在Y/D=0～0.5高度位置的流向速度偏小，而零方程模型在Y/D=0.5～2.0位置預測的流向速度偏小。這是由于零方程模型和層流模型只能給定速度分布與總溫，不能給定進口湍流度，導致射流出口附近流場與實驗偏差很大。兩方程模型預測的射流孔出口處流向速度很接近，與實驗結果也吻合得很好。在距射流出口下游1D位置，不同湍流模型預測的流向速度出現了明顯差異。其中層流和零方程模型在壁面附近的趨勢與實驗結果相反，均在壁面附近向正流向運動；其它湍流模型在壁面附近的趨勢與實驗相同，均在壁面附近出現倒流。k-ε模型預測的倒流峰值低于其它湍流模型，SST、k-ω、RNGk-ε模型均給出了偏高的回流速度。隨著X/D的增大，不同湍流模型之間的差異擴大。X/D=3位置，RNGk-ε模型得到的峰值為負數，而實驗結果表明均為正向流動。RNGk-ε模型、零方程模型和層流模型加速偏離實驗值，其中零方程預測速度嚴重偏離實驗值，層流方程在X/D=5位置速度峰值仍是負值。k-ε、k-ω、SST湍流模型結果偏離實驗值較小，但均偏低地預測了壁面附近的速度，k-ε模型結果與實驗值最為接近。

表1 網格分布Table 1 Grid distribution

圖5 不同湍流模型在Z/D=0位置上的流向速度對比Fig.5 Comparison between different turbulence models of streamwise velocity forZ/D=0

圖6不同湍流模型在Z/D=-0.5位置上的展向速度對比Fig.6 Comparison between different turbulence models of cross-tunnel velocity forZ/D=-0.5

圖6 給出了射流孔側面邊緣(Z/D=-0.5)不同湍流模型得到的展向速度與實驗結果對比。在射流孔側邊(X/D=0)位置，各模型結果與實驗結果的展向速度均為負值，表明流向均為遠離射流孔方向，且展向速度在距壁面很短一段距離內就上升到峰值。零方程在射流出口就偏差實驗值較遠，其預測峰值為0.50Vj，其它模型則在0.65Vj～0.70Vj之間。在X/D=1位置，除零方程模型外，其余模型所得結果較為吻合，均預測到近壁面的正峰值和距壁面0.50D的負峰值。在壁面附近，對轉渦受主流作用，展向速度為正值，朝射流孔流動；在對轉渦上半區，展向速度為負，遠離射流孔流動。兩個峰值的位置表明了對轉渦的范圍，而實驗結果并沒有表現出近壁面峰值的存在，計算的對轉渦區域與實驗結果吻合較好。隨著距射流出口距離的增加，對轉渦強度逐漸減小，各數值結果與實驗結果的偏差增大。

圖7給出了Z/D=-1位置上的法向速度分布。在X/D=0位置，零方程模型和層流模型均在壁面附近得到了負值，即指向壁面流動，該位置為對轉渦側邊緣；實驗結果表明法向速度為正值，其流向為遠離壁面方向，其它兩方程模型均與實驗結果吻合得很好。在X/D=1位置，壁面附近法向速度出現低峰值，對應對轉渦外側向下流動；在距壁面1D位置則出現正峰值，為主流在對轉渦位置以上加速，兩個峰值所在高度與圖6的一致。隨著流體向下游流動，展向速度變小，對轉渦強度逐漸減弱，兩方程湍流模型預測趨勢相同，RNGk-ε模型預測的低峰值過低，k-ε模型和k-ω模型與實驗結果吻合得很好。

圖7不同湍流模型在Z/D=-1位置上的法向速度對比Fig.7 Comparison between different turbulence models of vertical velocity for Z/D=-1

圖8 給出了射流出口中線(Z/D=0)位置上的湍動能(k/Vj)分布。在射流出口位置，層流模型高估了主流的湍動能，其它模型均低估了湍動能的峰值，這與流向速度沿法向的速度梯度(?U/?Y)相關。在距射流出口1D位置，各模型均得到了三個峰值，分別對應壁面射流、射流尾跡和剪切層流動，這些位置的?U/?Y較大，k-ε模型高估了峰值，SST、RNGk-ε模型低估了峰值，k-ω模型與實驗結果吻合較好。隨著流體向下游流動，湍動能峰值減小，但峰值所在高度提高。文獻[5]認為，數值結果的不準確與射流孔出口流場有關，在Ajersch的實驗中，射流氣體由儲氣室進入射流通道時，會在銳角處加速，從而導致射流進口氣體湍動能不均勻。因此在下文研究中，將采用實驗得到的射流通道出口截面參數作為數值模擬的射流邊界條件，從而對射流邊界條件的影響進行討論。綜上，相對于其它湍流模型，k-ε模型能較為準確地預測各速度分量和湍動能分布。

圖8 不同湍流模型在Z/D=0位置上的湍動能對比Fig.8 Comparison between different turbulence models of turbulence kinetic energy for Z/D=0

3.3射流邊界條件的影響

對兩種射流邊界條件進行計算：僅把射流進口位置參數作為射流出口邊界條件(CASE2)，和以實驗測得的射流出口速度分布、平均湍動能作為射流邊界條件(CASE3)。所采用的計算域不包括射流通道，并將結果與考慮射流通道的計算結果、實驗結果進行對比。計算網格均采用網格3，湍流模型為k-ε模型。

圖9給出了射流出口位置法向速度(V/Vj)分布。在射流孔迎風側至背風側，法向速度均在0.2Vj～1.7 Vj之間，表明實驗結果作為射流出口位置邊界條件可行。

圖9 射流出口位置法向速度分布Fig.9 Vertical velocity distribution at jet exit

圖10 不同射流邊界條件在Z=0位置上的流向速度對比Fig.10 Comparison between different jet conditions of streamwise velocity for Z/D=0

圖10 給出了Z/D=0位置上的流向速度對比，CASE1結果為考慮射流通道時的計算結果。在射流孔出口中心，CASE2的進口流向速度分量為零，CASE1中射流在冷氣通道內受主流擠壓而產生流向分量，在射流出口流向速度約為0.5Vj，CASE1、CASE3的分布與實驗結果吻合得很好。在X/D=1位置，實驗、CASE1、CASE3結果幾乎重合，由于CASE2的射流出口流向速度分布與實驗偏差較大，在X/D=1位置尾跡區的速度要低于實驗值。在壁面到0.25D高度內，CASE3的低峰值比CASE1的小，說明采用實驗結果作為邊界條件下，在距射流出口很短一段距離內的倒流加強。在下游發展過程中，CASE3在壁面附近的流向速度加速明顯，在X/D=3位置，CASE3尾跡區的流向速度已高于考慮射流通道的流向速度。在距射流出口較長位置上，CASE1和CASE2的結果與實驗結果的偏差逐漸擴大。然而在壁面附近，采用實驗結果作為進口邊界條件所得流速與實驗結果最為接近，CASE1與CASE2所得的流向速度均偏小。在X/D=3～5范圍內，CASE2結果在距壁面1D高度內的速度分布較實驗結果偏低。

圖11給出了方孔兩側(Z/D=-0.5)位置上的展向速度對比。在X/D=0位置，展向速度與實驗結果吻合較好，采用射流出口測量值為邊界條件后，在0.2D高度處的展向速度峰值為0.8Vj，高于考慮射流通道算例中的0.7Vj。對比X/D=0與X/D=1流向位置，所測低峰值絕對值由0.5Vj減小為0.3Vj，而峰值所在高度由0.2D升至0.5D。在X/D=1位置，壁面附近的正峰值高達0.6Vj，CASE1與CASE3在Y/D=0.5高度存在負峰值，CASE2所預測的負峰值位于偏高的Y/D=0.8位置，高估了對轉渦對的范圍。

圖12給出了Z/D=-1位置上的法向速度分布。對比CASE2、CASE3，雖然均不考慮射流通道，但CASE2在X/D=0位置得到的低峰值為負值，在靠近壁面有向下的流動，說明在射流出口對轉渦已形成，這與實驗結果不符；受射流擠壓，主流沿法向加速至0.15Vj。相對于CASE2，CASE3與實驗結果吻合得更好，但在主流區，CASE2、CASE3均低于實驗結果，表明射流出口速度分布對對轉渦的形成至關重要。在X/D=1位置，由前文分析已知，該位置位于對轉渦邊緣，其法向速度方向指向壁面，在對轉渦位置以上，主流向上繞過對轉渦加速流動，法向速度迅速增加，在距壁面4D范圍內，加速后大于0.05Vj。對比考慮射流通道結果，CASE2的正峰值所在高度最高。

圖11 不同射流邊界條件在Z/D=-0.5位置上的展向速度對比Fig.11 Comparison between different jet conditions of cross-tunnel velocity for Z/D=-0.5

圖12 不同射流邊界條件在Z=-1D位置上的法向速度對比Fig.12 Comparison between different jet conditions of vertical velocity for Z/D=-1

圖13 為不同射流邊界條件得到的湍動能對比。在射流孔出口位置，CASE2與CASE3預測的湍動能峰值大小接近，但CASE2預測峰值位置偏高，CASE3所得高度與實驗值吻合較好。在下游位置，各數值結果均與實驗結果吻合不好，在X/D=1位置，CASE2的湍動能比CASE1的低，這是由于主流與射流在射流通道內已開始相互作用，湍動能會在射流通道內部提高，在射流孔下游表現明顯，不考慮射流通道就不能對通道內的流動進行預測，湍動能偏低；湍動能與尾跡區內的流向速度梯度相關，所以CASE3的湍動能偏大。

圖13 不同射流邊界條件在Z=0位置上的湍動能對比Fig.13 Comparison between different jet conditions of turbulence kinetic energy for Z/D=0

4 結論

(1)在射流出口遠下游，各數值結果與實驗結果偏差均較大，加密下游網格，所得流向速度分布保持不變。射流孔下游網格分布可適當稀疏，有利于減少總體網格數目，從而減少計算量。

(2)解決冷氣摻混問題的常用湍流模型中，兩方程模型普遍比層流方程模型和零方程模型表現好。相對于其它湍流模型，k-ε模型結果與實驗值吻合最好，SST、k-ω、RNGk-ε模型均給出了偏高的回流速度；在距射流孔較遠的流向位置，所采用的湍流模型均高估了對轉渦強度；k-ε模型能較為準確地預測射流與主流相互作用產生的射流尾跡、對轉渦等流動特征。

(3)在氣膜冷卻數值研究中，應考慮射流通道，對射流通道內流動進行模擬，有助于得到準確的射流出口流場，從而得到冷氣覆蓋范圍、冷卻效果等參數。在冷卻通道未知的情況下，根據經驗在射流出口給定速度分布和湍流度作為射流的邊界條件，也有助于提高氣膜冷卻問題數值模擬精度。

[1]Zhou J M.A Computational and Experimental Investiga?tion of Gas Turbine Blade Film Cooling[D].Vancouver：University of British Columbia，1994.

[2]Ajersch P，Ketler S，Zhou J M，et al.Multiple Jets in a Cross-Flow:Detailed Measurements and Numerical Simu?lations[J].ASME Journal of Turbomachinery，1997，119：330—342.

[3]Hoda A，Acharya S.Predictions of a Film Coolant Jet in Cross-Flow with Different Turbulence Models[J].ASME Journal of Turbomachinery，2000，122：558—569.

[4]Hassan J S，Yavuzkurt S.Comparison of Four Different Two-Equation Models of Turbulence in Predicting Film Cooling Performance[R].ASME GT2006-90860，2006.

[5]Keimasi M R，Taeibi-Rahni M.Numerical Simulation of Jets in a Cross-Flow Using Different Turbulence Models [J].AIAA J，2001，39：2268—2277.

[6]Garg V K，Gaugler R E.Effect of Velocity and Tempera?ture Distribution at the Hole Exit on Film Cooling of Tur?bine Blades[R].ASME 95-GT-2，1995.

[7]Guo X，Schr?der W.Large Eddy Simulation of Film Cool?ing Flows[J].Comput.Fluids，2006，35：587—606.

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

Parameters Affecting Numerical Simulation Precision for Cooling Air Mixing

HE Hong-jiang，QI Lei，WANG Peng，LIU Huo-xing，ZOU Zheng-ping
(National Key Laboratory of Science and Technology on Aero-Engine Aero-Thermodynamics，School of Energy and Power Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

Several parameters which affect the precision of numerical simulation for cooling air mixing were investigated,including mesh,turbulence model and boundary conditions.The comparison between the results with these parameters and measurement at a blowing ratio of 0.5 was made.The meshs at the far downstream from the exit of injection could be coarser to reduce the amount of computation.Thek-εturbu?lence model produced the result which matched the measurement best.The distribution of velocity and tur?bulence kinetic energy at the exit of injection is non-uniform because of the flow within the jet channel,so the geometry of jet channel should be taken into consideration to elevate the precision of numerical simula?tion.Also,distribution of flow-field at the jet exit could be fixed to ensure the precision of computation.

aero-engine；film cooling；square jet in cross-flow；cooling air in cross-flow；

V231.3

A

1672-2620(2013)03-00 21-08

2012-09-10；

2013-03-26

武器裝備預研基金(9140C410101110C4101)，國家自然科學基金(51106004)

何宏疆(1987-)，男，新疆昌吉人，碩士研究生，主要從事葉輪機內部流動機理研究。