冷卻氣膜與熱障涂層交界處流場(chǎng)分析

李佳莉， 王玉璋， 王 星

(上海交通大學(xué)燃?xì)廨啓C(jī)研究院，上海 200240)

隨著航空發(fā)動(dòng)機(jī)及地面燃機(jī)的迅速發(fā)展，渦輪高溫葉片入口溫度不斷提高。目前國(guó)外新型軍用燃?xì)鉁u輪發(fā)動(dòng)機(jī)的燃?xì)鉁囟纫堰_(dá)1538～1871℃，而現(xiàn)役發(fā)動(dòng)機(jī)的高溫部件材料許用溫度均在1100℃以下［1，2］。在目前先進(jìn)燃?xì)鉁u輪發(fā)動(dòng)機(jī)中，高溫?zé)岱雷o(hù)涂層以及高效冷卻是渦輪葉片防護(hù)兩大關(guān)鍵技術(shù)。先進(jìn)的熱障涂層能夠在工作環(huán)境下降低熱端部件溫度170℃左右［3］。

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

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

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

1 數(shù)學(xué)模型

1.1 D2Q9格子Boltzmann方法

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

圖1 D2Q9模型Fig.1 D2Q9 model

在無(wú)相變多孔結(jié)構(gòu)內(nèi)部及表面流場(chǎng)中，密度分布函數(shù)的演化方程為:

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

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

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

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

1.2 邊界條件

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

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

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

1.2.1 左右邊界條件

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

1.2.2 上邊界條件

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

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

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

1.2.3 下邊界條件

下邊界為固體邊界條件，該固體邊界為無(wú)滑移邊界，由真實(shí)結(jié)構(gòu)知，在固體邊界上無(wú)流體分布，對(duì)于固體邊界有:

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

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

圖5 標(biāo)準(zhǔn)反彈格式Fig.5 Bounce-back scheme

2 結(jié)果分析與討論

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

此處取ε為小量即ε=10－6，當(dāng)兩個(gè)相鄰時(shí)層速度相對(duì)誤差小于設(shè)定值ε即error≤ε時(shí)，計(jì)算收斂。

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

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

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

圖7 不同結(jié)構(gòu)豎直方向速度分布圖Fig.7 Vertical velocity distribution of different coatings

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

圖9中a和b分別為柱狀和層狀結(jié)構(gòu)不同橫截面豎直方向速度波動(dòng)圖，選取的三個(gè)截面分別從涂層表面至遠(yuǎn)流場(chǎng)邊界，圖中縱坐標(biāo)取為流體豎直方向速度V與來(lái)流速度U的無(wú)量綱比。由圖可見(jiàn)兩種結(jié)構(gòu)涂層在靠近表面處速度波動(dòng)強(qiáng)烈，遠(yuǎn)流場(chǎng)區(qū)域速度波動(dòng)趨于平緩，說(shuō)明涂層表面微結(jié)構(gòu)改變了流場(chǎng)流動(dòng)狀況。對(duì)于涂層表面處豎直方向上的速度波動(dòng)，層狀結(jié)構(gòu)的波動(dòng)幅度明顯強(qiáng)于柱狀結(jié)構(gòu)，最大可達(dá)2.5%。涂層表面流體溫度場(chǎng)與流場(chǎng)相互耦合，進(jìn)而會(huì)影響流體溫度場(chǎng)的分布。

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

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

3 結(jié)論

(1)采用LBM方法研究涂層微結(jié)構(gòu)對(duì)氣膜流動(dòng)的影響，能夠較好的反映真實(shí)的微觀流場(chǎng)分布;

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

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

［1］文生瓊，何愛(ài)杰，王皓.熱障涂層在航空發(fā)動(dòng)機(jī)渦輪葉片上的應(yīng)用［J］.燃?xì)鉁u輪試驗(yàn)與研究，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陶瓷的熱導(dǎo)率［J］.人工晶體學(xué)報(bào)，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］魏紹斌，陸峰，何利民，等.熱障涂層制備技術(shù)及陶瓷層材料的研究進(jìn)展［J］.熱噴涂技術(shù)，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］郭婷婷，金建國(guó)，李少華，等.不同出射角度對(duì)氣膜冷卻流場(chǎng)的影響［J］.中國(guó)電機(jī)工程學(xué)報(bào)，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］張紅松，劉振啟，關(guān)紹康.等離子噴涂納米與微米 YSZ熱障涂層的孔隙結(jié)構(gòu)比較［J］.表面技術(shù)，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方法的飽和土體細(xì)觀滲流場(chǎng)［J］.排灌機(jī)械工程學(xué)報(bào)，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方法的原理及應(yīng)用［M］，北京:科學(xué)出版社，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方法的理論及應(yīng)用［M］，北京:科學(xué)出版社，2008:31－56.

［15］凌錫祥，王玉璋，王星，等.層狀熱障涂層孔隙微結(jié)構(gòu)對(duì)其隔熱性能影響的數(shù)值研究［J］.中國(guó)有色金屬學(xué)報(bào)，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］凌錫祥，王玉璋，王星.基于柱狀結(jié)構(gòu)的熱障涂層隔熱性能數(shù)值研究［J］.航空材料學(xué)報(bào)，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.