李佳莉， 王玉璋， 王 星

(上海交通大學燃氣輪機研究院，上海 200240)

隨著航空發(fā)動機及地面燃機的迅速發(fā)展，渦輪高溫葉片入口溫度不斷提高。目前國外新型軍用燃氣渦輪發(fā)動機的燃氣溫度已達1538～1871℃，而現(xiàn)役發(fā)動機的高溫部件材料許用溫度均在1100℃以下［1，2］。在目前先進燃氣渦輪發(fā)動機中，高溫熱防護涂層以及高效冷卻是渦輪葉片防護兩大關鍵技術。先進的熱障涂層能夠在工作環(huán)境下降低熱端部件溫度170℃左右［3］。

熱障涂層目前主要的制備工藝是電子束物理氣相沉淀法(EB-PVD)和等離子噴涂(APS)，這兩種方法制備的熱障涂層都是各向異性的非致密多孔狀介質，其孔隙率在5% ～20%之間［4，5］。熱噴涂涂層呈層狀結構，表面比較粗糙;而物理氣相沉淀涂層呈柱狀結構，表面較為光滑，但其孔隙率較高且垂直于涂層基體表面;不同的制備工藝以及長時間服役導致的涂層微孔和微裂紋的變化，會對葉片表面高溫燃氣與冷卻氣膜流場產(chǎn)生不同的影響。

目前對葉片表面流固耦合流場的研究主要是采用傳統(tǒng)的數(shù)值模擬從宏觀上求解控制方程，探索冷氣射流和高溫衡流相互摻混的漩渦結構和流場分布形成及發(fā)展的影響因素如氣膜孔孔型、孔幾何參數(shù)和氣動參數(shù)等［6～8］。高溫燃氣及冷卻射流會直接沖擊在涂層表面，而涂層微孔結構尺度大約為10 μm［9］，目前從微觀角度揭示熱障涂層微結構對葉片表面流動特性影響的文獻相對較少。

格子Boltzmann方法(LBM)是20世紀80年代基于分子動理論發(fā)展起來的一種新型的數(shù)值模擬，LBM從介觀層面上考慮粒子團的遷移和碰撞，通過粒子分布函數(shù)演化來得到宏觀問題的解，在研究具有復雜表面結構如多孔介質內流體流動與換熱及滲流現(xiàn)象中得到了廣泛的應用［10～13］。本工作基于不同熱障涂層的微結構，研究冷卻氣膜和熱障涂層交界處微觀流動，分析不同涂層微結構對冷卻氣膜穩(wěn)定性的影響。

1 數(shù)學模型

1.1 D2Q9格子Boltzmann方法

本工作采用Boltzmann-BGK碰撞松弛模型中的格子波爾茲曼模型的二維九速度(D2Q9)不可壓格子Boltzmann模型［14］，該模型可用來計算微觀尺度下多孔材料、具有復雜表面結構的流固耦合流場分布;該模型的離散速度模型(圖1)和速度配置(式(1))如下:

圖1 D2Q9模型Fig.1 D2Q9 model

在無相變多孔結構內部及表面流場中，密度分布函數(shù)的演化方程為:

式中:r是坐標矢量，t是時間，τf是流體無量綱松弛時間;fi(r，t)和feqi(r，t)分別為與i方向相對應的粒子密度分布函數(shù)和平衡態(tài)密度分布函數(shù)。平衡態(tài)分布函數(shù)feqi(r，t)和無量綱松弛時間分別為:

其中ωf為密度分布函數(shù)的權系數(shù)，u為流體宏觀流動速度，δt和δx分別為格子單位下對應的時間步長和網(wǎng)格步長，cs為格子聲速，ν為流體黏性系數(shù)。

流體宏觀密度和速率可由下式確定:

式中，u為流體宏觀流動速度，ei為上文中圖1所示二維九速度模型中不同方向的速度。

1.2 邊界條件

真實結構熱障涂層是多孔狀介質，四參數(shù)隨機生長方法在獲取涂層幾何結構模型上取得了較好的結果［15～17］，本工作模擬熱障涂層表面結構與冷卻氣膜耦合作用下的流場，以電鏡掃描獲得的柱狀以及層狀熱障涂層真實結構為基礎，如圖2，3所示，其中黑色部分代表涂層固體材料骨架，灰白色區(qū)域為涂層孔隙及流體流動區(qū)域，經(jīng)過灰度處理及Matlab像素識別功能，獲得了0，1數(shù)值結構網(wǎng)格，局部網(wǎng)格結構如圖2，3所示，經(jīng)像素識別，圖2中像素點為711×357，網(wǎng)格數(shù)為 NX×NY=711×357，其中涂層厚度約為78μm，每一格厚度為0.42μm;圖3中像素點為716×281，網(wǎng)格數(shù)為NX×NY=716×281，其中涂層厚度約為 277μm，每一格厚度為 1.964μm，邊界條件設置如圖中所示:

圖2 柱狀結構涂層邊界設置及網(wǎng)格結構Fig.2 Boundary condition and grid of columnar TBC

圖3 APS噴涂結構涂層邊界設置及網(wǎng)格結構Fig.3 Boundary condition and grid of TBC made by APS

1.2.1 左右邊界條件

本工作模擬的流動區(qū)域為局部小范圍區(qū)域，在穩(wěn)態(tài)流動過程中，流體的流速和密度等宏觀參數(shù)在左右邊界變化很小，故在模擬中將流動區(qū)域的左右邊界條件設置為周期性邊界條件，同時忽略重力對流場的影響，如圖4所示，以進口邊界為例，進口邊界節(jié)點未知分布函數(shù)f1，f5，f8，可由出口處相應節(jié)點粒子分布函數(shù)遷移得到，處理格式如下:

1.2.2 上邊界條件

上邊界為流體流動邊界，模擬中給定上邊界水平方向流動速度為U，豎直方向速度無脈動;密度采用與其相鄰的內推節(jié)點的密度值，即:

圖4 周期性邊界條件示意圖Fig.4 Periodic boundary condition

上邊界的分布函數(shù)采用Guo等［18］提出的非平衡外推格式:

1.2.3 下邊界條件

下邊界為固體邊界條件，該固體邊界為無滑移邊界，由真實結構知，在固體邊界上無流體分布，對于固體邊界有:

1.2.4 不規(guī)則固體骨架邊界

反彈格式［19］常用于靜止固體邊界以及內部與流體接觸的不規(guī)則固體骨架邊界上，常用的處理方法就是對邊界上的粒子作彈回處理，如圖5所示邊界節(jié)點(i，1)，采用標準反彈格式來處理，即:

圖5 標準反彈格式Fig.5 Bounce-back scheme

2 結果分析與討論

基于上述模型，在Microsoft Visual Studio 2010平臺上開發(fā)了微觀尺度下流場分析軟件，實現(xiàn)了對多孔涂層表面及內部氣膜流動的數(shù)值模擬，其收斂條件設置為:

此處取ε為小量即ε=10－6，當兩個相鄰時層速度相對誤差小于設定值ε即error≤ε時，計算收斂。

基于兩種工藝制備的涂層微結構形成的計算網(wǎng)格，模擬流動參數(shù)如下:涂層表面摻冷氣體溫度Tf=1212℃，涂層底面溫度Tc=844℃，上邊界流體水平方向速度設為U=0.02(換算為實際流動速度約為23m/s)，豎直方向速度無脈動;模擬中時間步長δt=1，格子步長δx=δy=1，圖2所示模型流動無量綱松弛時間取為τf=0.602，圖3所示模型流動無量綱松弛時間取為τf=0.523，流動計算中均選取格子單位下的無量綱量。

圖6為不同涂層結構水平方向速度分布。由圖可見，柱狀結構涂層表面水平方向速度變化平穩(wěn)，這是由于柱狀結構表面較為平整;而層狀結構涂層表面水平方向速度出現(xiàn)局部波動，這是由于層狀結構表面比較粗糙。圖7為不同涂層結構豎直方向速度分布，由圖可見，柱狀結構涂層孔隙中存在微弱的氣體流動，降低了涂層的隔熱性能，這是由于柱狀結構涂層內部孔隙垂直于涂層表面，且孔隙往往與外部流場相通，孔隙率較高;層狀結構涂層豎直方向速度波動幅度較大，氣膜穩(wěn)定性較差，對表面流場的黏性底層會產(chǎn)生較強影響。

圖6 不同涂層結構水平方向速度分布(a)柱狀結構;(b)APS噴涂層結構Fig.6 Horizontal velocity distribution of different coatings(a)columnar structure;(b)layered structure by APS

圖7 不同結構豎直方向速度分布圖Fig.7 Vertical velocity distribution of different coatings

圖8中a和b分別為柱狀結構涂層和APS噴涂層狀結構涂層表面的氣膜系統(tǒng)流線圖及局部放大圖。從放大圖中可以清晰地看到流體在涂層表面的微孔間隙中產(chǎn)生了漩渦，緊貼涂層表面處流線明顯波動，對流動產(chǎn)生了擾動。

圖9中a和b分別為柱狀和層狀結構不同橫截面豎直方向速度波動圖，選取的三個截面分別從涂層表面至遠流場邊界，圖中縱坐標取為流體豎直方向速度V與來流速度U的無量綱比。由圖可見兩種結構涂層在靠近表面處速度波動強烈，遠流場區(qū)域速度波動趨于平緩，說明涂層表面微結構改變了流場流動狀況。對于涂層表面處豎直方向上的速度波動，層狀結構的波動幅度明顯強于柱狀結構，最大可達2.5%。涂層表面流體溫度場與流場相互耦合，進而會影響流體溫度場的分布。

圖8 系統(tǒng)流線圖及局部流線圖 (a)柱狀結構;(b)APS噴涂層狀結構Fig.8 System flow diagram and local streamlines (a)columnar structure;(b)layered structure by APS

圖9 不同橫截面豎直方向速度波動圖 (a)柱狀結構;(b)APS噴涂層狀結構Fig.9 Distribution of vertical velocity at different vertical sections (a)columnar structure;(b)layered structure by APS

3 結論

(1)采用LBM方法研究涂層微結構對氣膜流動的影響，能夠較好的反映真實的微觀流場分布;

(2)熱障涂層微結構對冷卻氣膜流動產(chǎn)生較大影響，涂層表面凹陷處產(chǎn)生漩渦，柱狀結構涂層內部豎直孔隙中存在微弱的氣體流動，降低了涂層的隔熱性能;

(3)由于制備工藝不同，層狀結構涂層表面比較粗糙，相比柱狀結構對表面氣膜流動的影響更加劇烈，豎直方向速度波動幅度可達2.5%(與來流速度無量綱比)，氣膜穩(wěn)定性差;但在遠流場區(qū)域兩種結構涂層速度波動趨于平緩。

［1］文生瓊，何愛杰，王皓.熱障涂層在航空發(fā)動機渦輪葉片上的應用［J］.燃氣渦輪試驗與研究，2009，22(1):59－62.(WEN S Q，HE A J，WANG H.Development of TBCs on turbine blade of aero-engine［J］.Gas Turbine Experiment and Research，2009，22(1):59 －62.)

［2］李振軍，張紅松，魏媛.(La0.75Sm0.25)2Zr2O7陶瓷的熱導率［J］.人工晶體學報，2009，38(1):117－121.(LI Z J，ZHANG H S，WEI Y.Thermal conductivity of(La0.75Sm0.25)2Zr2O7ceramic［J］.Journal of Synthetic Crystals，2009，38(1):117 －121.)

［3］GUO H B，VAEN R，STVER D.Atmospheric plasma sprayed thick thermal barrier coatings with high segmentation crack density［J］.Surface and Coatings Technology，2004，186(3):353－363.

［4］魏紹斌，陸峰，何利民，等.熱障涂層制備技術及陶瓷層材料的研究進展［J］.熱噴涂技術，2013，5(1):31－37.(WEI S B，LU F，HE L M，et al.Progress in processing techniques and ceramic materials of thermal barrier coatings［J］.Thermal Spray Technology，2013，5(1):31 －37.)

［5］SCHULZ U，LEYENS C，F(xiàn)RITSCHER K，et al.Some recent trends in research and technology of advanced thermal barrier coatings［J］.Aerospace Science and Technology，2003，7(1):73－80.

［6］郭婷婷，金建國，李少華，等.不同出射角度對氣膜冷卻流場的影響［J］.中國電機工程學報，2006，26(16):117－121.(GUO T T，JIN J G，LI S H，et al.A numerical simulation in film cooling flows with different injection angles［J］.Proceedings of the CSEE，2006，26(16):117－121.)

［7］WALTERS D K，LEYLEK J H.A detailed analysis of filmcooling physics:Part I—streamwise injection with cylindrical holes［C］//ASME 1997 International Gas Turbine and Aeroengine Congress and Exhibition.USA:American Society of Mechanical Engineers，1997:V003T09A052-V003T09A052.

［8］BOHN D，REN J，KUSTERER K.Conjugate heat transfer analysis for film cooling configurations with different hole geometries［C］//ASME Turbo Expo 2003，collocated with the 2003 InternationalJointPowerGeneration Conference.USA:American Society of Mechanical Engineers，2003:247－256.

［9］張紅松，劉振啟，關紹康.等離子噴涂納米與微米 YSZ熱障涂層的孔隙結構比較［J］.表面技術，2010，39(5):4－7.(ZHANG H S，LIU Z Q，GUAN S K.Comparison of pore structure between plasma-sprayed nano-YSZ and micron-YSZ thermal barrier coatings［J］.Surface Technology，2010，39(5):4 －7.)

［10］申林方，王志良，李邵軍.基于格子 Boltzmann方法的飽和土體細觀滲流場［J］.排灌機械工程學報，2014，32(10):883－887.(SHEN L F，WANG Z L，LI S J.Mesoscopic seepage field of saturated soil with Lattice Boltzmann method［J］.Journal of Drainage and Irrigation Machinery Engineering，2014，32(10):883 －887.)

［11］郭照立，鄭楚光.格子Boltzmann方法的原理及應用［M］，北京:科學出版社，2008:8.

［12］MADADI M，SAHIMI M.Lattice Boltzmann simulation of fluid flow in fracture networks with rough，self-affine surfaces［J］.Physical Review(E)，2003，67(2):026309.

［13］WANG J，WANG M，LI Z.A lattice Boltzmann algorithm for fluid-solid conjugate heat transfer［J］.International Journal of Thermal Sciences，2007，46(3):228－234.

［14］何雅玲，王勇，李慶.格子Boltzmann方法的理論及應用［M］，北京:科學出版社，2008:31－56.

［15］凌錫祥，王玉璋，王星，等.層狀熱障涂層孔隙微結構對其隔熱性能影響的數(shù)值研究［J］.中國有色金屬學報，2015，25(2):408－414.(LING X X，WANG Y Z，WANG X，et al.Numerical study of effect of pore microstructure of layered thermal barrier coatings on thermal insulation performance［J］.The Chinese Journal of Nonferrous Metals，2015，25(2):408－414.)

［16］凌錫祥，王玉璋，王星.基于柱狀結構的熱障涂層隔熱性能數(shù)值研究［J］.航空材料學報，2014，34(5):69－74.(LING X X，WANG Y Z，WANG X.Numerical study of thermal insulation properties of thermal barrier coating based on columnar structure［J］.Journal of Aeronautical Materials，2014，34(5):69 －74.)

［17］WANG Y Z，WANG X，WENG Y W.Analysis of effective thermal conductivity of thermal barrier coatings based on microstructure［C］//ASME Turbo Expo 2014:Turbine Technical Conference and Exposition.USA:American Society of Mechanical Engineers，2014: V006T22A015-V006T22A015.

［18］GUO Z L，ZHENG C G，SHI B C.Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method［J］.Chinese Phys-ics，2002，11(4):0366 －0374.

［19］ZIEGLER D P.Boundary conditions for lattice Boltzmann simulations［J］.Journal of Statistical Physics，1993，71(5/6):1171－1177.