渦輪葉柵通道內(nèi)顆粒物沉積過(guò)程的數(shù)值模擬

楊曉軍， 祝佳雄

1.中國(guó)民航大學(xué) 航空工程學(xué)院， 天津 300300 2.中國(guó)民航大學(xué) 中歐航空工程師學(xué)院， 天津 300300

渦輪葉柵通道內(nèi)顆粒物沉積過(guò)程的數(shù)值模擬

楊曉軍1,2，*， 祝佳雄2

1.中國(guó)民航大學(xué) 航空工程學(xué)院， 天津 300300 2.中國(guó)民航大學(xué) 中歐航空工程師學(xué)院， 天津 300300

為了更加準(zhǔn)確獲得顆粒物在渦輪中的沉積分布，以某渦輪葉片為模型，選用最接近航空發(fā)動(dòng)機(jī)內(nèi)部顆粒組成的Jim Bridger Power Station (JBPS)顆粒為污染物，同時(shí)，利用C++編寫(xiě)合適的User Defined Function (UDF)經(jīng)過(guò)調(diào)試來(lái)分析顆粒沉積后葉片邊界的復(fù)雜變形和邊界網(wǎng)格依賴(lài)于時(shí)間變化的重構(gòu)生成，在考慮每個(gè)時(shí)間步長(zhǎng)顆粒沉積在葉片上從而改變?nèi)~片幾何特性和換熱特性的情況下，來(lái)深入研究顆粒物沉積在葉片的整個(gè)過(guò)程，最終得出了沉積的分布情況，并且通過(guò)數(shù)值研究結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比，驗(yàn)證了網(wǎng)格重構(gòu)與融合程序的合理性、準(zhǔn)確性。根據(jù)葉片變形情況預(yù)測(cè)腐蝕的發(fā)生情況。隨后，改變主流溫度、顆粒直徑來(lái)研究顆粒沉積特性。結(jié)果表明：顆粒主要沉積在葉片壓力面中部，但會(huì)使得葉片前緣和壓力面中部均產(chǎn)生明顯變形；葉片前緣由于顆粒沉積使得粗糙度增加形成鋸齒形，最先遭受腐蝕；顆粒直徑影響顆粒沉積的分布與沉積生長(zhǎng)速度；只影響沉積速度，并不改變沉積分布。

沉積； 變形； UDF； 腐蝕； 渦輪葉片； 粒徑； 數(shù)值模擬

隨著燃?xì)廨啓C(jī)在能源供應(yīng)上的提高以及天然氣資源陸續(xù)被耗盡，尋找替代燃料開(kāi)始變得非常重要。然而，替代燃料含有微量的煤灰和其他雜質(zhì)。這些雜質(zhì)會(huì)沉積在葉片和渦輪表面從而影響它們的熱交換特性，并且會(huì)降低氣膜冷卻的效率。同時(shí)伴隨著世界工業(yè)化的進(jìn)步，相當(dāng)多的地區(qū)包含了更多雜質(zhì)顆粒，從而引發(fā)嚴(yán)重的安全以及性能問(wèn)題。因此，顆粒物在航空發(fā)動(dòng)機(jī)內(nèi)部的沉積以及由此帶來(lái)的發(fā)動(dòng)機(jī)性能的損失相關(guān)問(wèn)題和所引發(fā)的環(huán)境問(wèn)題[1]已經(jīng)成為國(guó)際上高度重視的一個(gè)研究課題[2-5]。

為了提高渦輪的工作效率，預(yù)測(cè)渦輪葉片的壽命，也為了優(yōu)化葉片的設(shè)計(jì)，必須深入研究葉片表面的換熱特性，在眾多影響葉片換熱特性的因素當(dāng)中，顆粒的沉積是唯一隨著葉片工作時(shí)間的增長(zhǎng)而對(duì)葉片換熱影響愈加重要的關(guān)鍵因素，這是因?yàn)轭w粒的沉積主要影響葉片表面的粗糙度，甚至?xí)斐蓪?duì)葉片的熱腐蝕，因此研究顆粒沉積過(guò)程對(duì)于葉片設(shè)計(jì)與維護(hù)都有重要的意義[6]。

顆粒沉積與其影響在國(guó)外研究得非常早。在沉積機(jī)理的研究方面，David和Karen[7]進(jìn)行實(shí)驗(yàn)研究發(fā)現(xiàn)，熔融程度不同的粒子沉積對(duì)于冷卻效果的影響可以忽略。Wright等[8]發(fā)現(xiàn)顆粒很有可能通過(guò)相互碰撞或者撞擊葉片而重新分配并且沉積在葉片上。同時(shí)，Wright通過(guò)實(shí)驗(yàn)證明了處于端壁處的前緣氣膜冷卻區(qū)域更易吸引顆粒沉積。對(duì)于顆粒黏附的主要影響因素的實(shí)驗(yàn)研究方面，David和Kaven[7]主要研究了冷卻孔的結(jié)構(gòu)對(duì)于粒子沉積帶來(lái)影響的抵御作用。根據(jù)實(shí)驗(yàn)結(jié)果，槽式孔可以有效地降低由于顆粒沉積而損失的冷卻效率。同時(shí)，他們通過(guò)對(duì)前緣氣膜冷卻與端壁冷卻進(jìn)行實(shí)驗(yàn)，得出橫向槽可以更有效地改善端壁冷卻的效率。Ahluwalia等[9]研究了煤灰在渦輪葉片上面的沉積和黏附現(xiàn)象。他們發(fā)現(xiàn)壓力面的黏著系數(shù)在0.000 3～0.11范圍內(nèi)變化并且與主燃?xì)夂捅诿鏈囟扔兄軓?qiáng)的聯(lián)系，但是沖擊角的變化對(duì)于黏著系數(shù)的影響是可以忽略的。對(duì)于沉積與腐蝕的關(guān)系以及影響沉積因素的研究方面，Crosby等[10]利用渦輪加速沉積設(shè)備(Turbine Accelerated Deposition Facility，TADF)研究顆粒的大小、氣體的溫度和沉積面的溫度如何影響沉積過(guò)程。他們發(fā)現(xiàn)當(dāng)氣體溫度跨過(guò)1 230 K 這個(gè)門(mén)檻，沉積增長(zhǎng)會(huì)越來(lái)越快。同時(shí)發(fā)現(xiàn)，保持主流氣體溫度不變，隨著沉積面溫度的降低，沉積的粗糙度會(huì)下降(即沉積更加均勻)。

國(guó)內(nèi)方面航空發(fā)動(dòng)機(jī)制造起步較晚，渦輪進(jìn)口溫度不夠高，因此研究顆粒物沉積起步較晚，關(guān)于沉積物的數(shù)值與實(shí)驗(yàn)研究很少。火箭固沖發(fā)動(dòng)機(jī)研究者研究了顆粒組分對(duì)于沉積情況的影響，主要針對(duì)液體燃料。王德全等[11]研究了固沖發(fā)動(dòng)機(jī)中顆粒沉積的主要分布位置，并無(wú)其他細(xì)致研究。關(guān)于航空發(fā)動(dòng)機(jī)方面，對(duì)顆粒沉積生長(zhǎng)研究甚少，國(guó)內(nèi)大多數(shù)學(xué)者研究與顆粒沉積相似的機(jī)翼結(jié)冰生長(zhǎng)過(guò)程，以及結(jié)冰對(duì)于機(jī)翼的影響等，其對(duì)于顆粒沉積的研究有著借鑒作用。2004年，張大林和陳維建[12]應(yīng)用Navier-Stokes方程和紊流模型計(jì)算結(jié)冰表面外的空氣流場(chǎng)。在研究過(guò)程中發(fā)現(xiàn)數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果總是有相對(duì)較大的誤差，經(jīng)過(guò)仔細(xì)研究之后，發(fā)現(xiàn)是由于網(wǎng)格的重構(gòu)與融合方面沒(méi)有定義。隨后2005年，陳維建和張大林[13]提出了一種可考慮粗糙度影響的瘤狀冰結(jié)冰過(guò)程數(shù)值模擬方法。并且與國(guó)外文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù)對(duì)比，結(jié)果吻合較好。2011年，孫志國(guó)等[14]專(zhuān)門(mén)對(duì)于機(jī)翼結(jié)冰過(guò)程中二維以及三維的網(wǎng)格重構(gòu)進(jìn)行了研究，提出了網(wǎng)格重構(gòu)與融合的思想，并且將網(wǎng)格模塊NI-GRID集成于結(jié)冰計(jì)算軟件NUAA-ICE3D中，并對(duì)二維翼型和三維機(jī)翼的網(wǎng)格生成與重構(gòu)進(jìn)行了驗(yàn)證。

這些研究為顆粒物沉積生長(zhǎng)方面的探索提供了重要的理論依據(jù)和實(shí)驗(yàn)數(shù)據(jù)，但是這些研究對(duì)于全面探索顆粒物沉積的過(guò)程機(jī)理還有一定的不足，有必要進(jìn)一步展開(kāi)數(shù)值模擬研究來(lái)彌補(bǔ)實(shí)驗(yàn)研究過(guò)程中所不能覆蓋的部分。本文嘗試?yán)糜袔缀螖?shù)據(jù)以及沉積實(shí)驗(yàn)數(shù)據(jù)的葉片來(lái)進(jìn)行顆粒物在渦輪葉片方面的沉積數(shù)值研究，與實(shí)驗(yàn)結(jié)果相結(jié)合檢驗(yàn)數(shù)值研究的可靠性，并分析影響顆粒物沉積的因素與顆粒物沉積對(duì)于渦輪葉片的影響。

1 顆粒沉積數(shù)值模擬方法

1.1 控制方程

關(guān)于流體方面的三大方程，根據(jù)流體力學(xué)中的無(wú)量綱化結(jié)果，可以得知，在航空發(fā)動(dòng)機(jī)工作的情況下，可以得到

(1)

由于渦輪內(nèi)部流動(dòng)的復(fù)雜性和高速度性，在整個(gè)渦輪葉片表面溫度的梯度的分布，在微米以及亞微米級(jí)別的顆粒，其主要受4個(gè)力的作用: 阻力，Saffman升力、熱泳力(熱遷移力)以及布朗力。

這其中，在一定時(shí)間內(nèi)起主導(dǎo)作用的是阻力，隨著時(shí)間延長(zhǎng)，Saffman升力占的比重會(huì)越來(lái)越大，而熱泳力只在溫度梯度大的地方稍微對(duì)固體的運(yùn)動(dòng)有些影響，就葉柵通道全局來(lái)看，影響不大。布朗力是分子隨機(jī)運(yùn)動(dòng)而產(chǎn)生的力，只會(huì)對(duì)整體的運(yùn)動(dòng)有輕微的波動(dòng)性影響，并不會(huì)影響固體顆粒宏觀的運(yùn)動(dòng)趨勢(shì)。進(jìn)一步簡(jiǎn)化之后，對(duì)固體顆粒列三大方程最終無(wú)量綱化的結(jié)果為

(2)

(3)

(4)

式中：St為Stokes數(shù)，其表示顆粒松弛時(shí)間和流體特征時(shí)間的比；下標(biāo)p表示顆粒，下標(biāo)conv表示對(duì)流，下標(biāo)k表示動(dòng)量。

由式(2)～式(4)可以看出，固體顆粒的運(yùn)動(dòng)受到Stokes數(shù)、主流溫度和主流速度的影響。其中，Stokes數(shù)主要是由顆粒的直徑來(lái)決定的，主流溫度與主流速度都是流場(chǎng)特性，通過(guò)改變進(jìn)口溫度與進(jìn)出口壓強(qiáng)比實(shí)現(xiàn)。

1.2 沉積模型

Tafti和Sreedharan[15]基于顆粒的黏性發(fā)展了一種沉積模型。顆粒的黏性隨著溫度改變并且可以根據(jù)煤灰顆粒的性質(zhì)進(jìn)行預(yù)測(cè)。煤灰顆粒的軟化溫度被認(rèn)為是臨界的沉積溫度，記為T(mén)s。高于臨界溫度的顆粒100%會(huì)沉積，低于臨界溫度的顆粒有根據(jù)以下函數(shù)的沉積概率

(5)

式中：μcrit為顆粒在臨界溫度時(shí)候的黏度；μTp為熔融顆粒在目前溫度的黏度。

Senior和Srinivasachar[16]提出了一種根據(jù)煤灰顆粒的組分來(lái)計(jì)算煤灰顆粒黏性的方法。隨后，N’Dala等[17]顯示了溫度和黏性之間的關(guān)系式

(6)

式中：A和B是依賴(lài)于化學(xué)組分的常數(shù)，Tp為顆粒溫度。Senior和Srinivasachar通過(guò)實(shí)驗(yàn)做出了A和B的變化曲線。更詳細(xì)的情況以及機(jī)理在文獻(xiàn)[15]中。

由于研究對(duì)象是第一級(jí)渦輪，渦輪進(jìn)口溫度較高，并且渦輪溫度場(chǎng)分布較速度場(chǎng)梯度變化更加劇烈，因此，數(shù)值模擬基于此沉積模型進(jìn)行模擬研究，并且與真實(shí)沉積結(jié)果進(jìn)行對(duì)比，來(lái)驗(yàn)證在渦輪工作條件下，溫度場(chǎng)對(duì)于沉積的影響是主要的。

2 網(wǎng)格的重構(gòu)與融合

2.1 網(wǎng)格生成思想

首先，利用橢圓型偏微分方程的方法生成網(wǎng)格。因?yàn)闄E圓形偏微分方程生成的網(wǎng)格有優(yōu)秀的光滑性，并且在橢圓偏微分方程的右側(cè)可以加入源項(xiàng)來(lái)控制調(diào)節(jié)生成網(wǎng)格的正交性。同時(shí)，橢圓型偏微分方程生成的網(wǎng)格質(zhì)量更加容易調(diào)節(jié)，利用Hilgenstock法源項(xiàng)可以有效地提高網(wǎng)格的質(zhì)量，從而保證迭代的穩(wěn)定性。

ξi=Pi

(7)

另一種方法是利用求解擴(kuò)散方程來(lái)進(jìn)行網(wǎng)格邊界的調(diào)整，擴(kuò)散光順調(diào)整即通過(guò)求解擴(kuò)散房產(chǎn)成對(duì)于網(wǎng)格的重構(gòu)進(jìn)行光順處理，從而優(yōu)化網(wǎng)格的質(zhì)量，減少產(chǎn)生負(fù)體積的危險(xiǎn)，其擴(kuò)散方程為

(8)

式中：

Υ=1/Vα

(9)

其中：α為擴(kuò)散參數(shù)；V為正則體積，其分別代表了取決于壁面距離的擴(kuò)散光順與取決于單元體積的擴(kuò)散光順。

2.2 網(wǎng)格的重構(gòu)與融合思想

網(wǎng)格重構(gòu)和沉積模型是構(gòu)成User Defined Function (UDF)的主要部分，其編寫(xiě)和選取方法對(duì)結(jié)果的影響非常大。目前發(fā)現(xiàn)主流的重構(gòu)方法主要包括以下2種：

1) 法線重構(gòu)法

將顆粒的體積分別投影到3個(gè)正方向上面，分別計(jì)算3個(gè)正方向上面的位移量，然后再將其合成為總的變形量，以二維來(lái)舉例，則可利用以下關(guān)系式進(jìn)行網(wǎng)格重構(gòu)與融合：

(10)

(11)

(12)

(13)

三維的情況只是多了z方向，此種方法依然需要對(duì)網(wǎng)格進(jìn)行光順處理，否則依然會(huì)很快出現(xiàn)負(fù)體積。

2) 順序重構(gòu)法[13]

首先，假設(shè)不存在顆粒沉積之后由于形狀而產(chǎn)生的滾轉(zhuǎn)以及其余位移形態(tài)選取貼近葉片表面的兩層網(wǎng)格為例，粗線表示葉片邊緣。虛線表示該時(shí)間步長(zhǎng)之內(nèi)顆粒沉積的總量。對(duì)于計(jì)算網(wǎng)格邊界而言，只需要確定陰影邊界即圖1中虛線與法線部分的交點(diǎn)坐標(biāo)即可，這是由于顆粒沉積沿著貼面網(wǎng)格的法線方向，而結(jié)構(gòu)網(wǎng)格即六面體網(wǎng)格的邊恰好是貼面的法線方向，因此極大簡(jiǎn)化了網(wǎng)格的重構(gòu)過(guò)程。具體的計(jì)算過(guò)程如下：

① 根據(jù)顆粒的沉積模型來(lái)設(shè)定顆粒沉積的條件。

② 計(jì)算顆粒單位時(shí)間步長(zhǎng)內(nèi)的沉積數(shù)量和沉積質(zhì)量，隨后根據(jù)顆粒的密度求解變形量。

③ 尋找最近的沒(méi)有顆粒沉積的控制體，如圖1(b)中的控制體8。

④ 從控制體7開(kāi)始計(jì)算各控制體內(nèi)顆粒沉積的量與形狀。具體方法為：根據(jù)控制體7內(nèi)的沉積量計(jì)算出顆粒在控制體7內(nèi)沉積的體積，然后從a點(diǎn)引射線到控制體6與7的交界處，即d點(diǎn)，使得Δacd與顆粒沉積的體積相等，即可求解出d點(diǎn)坐標(biāo)；之后的方法與計(jì)算d點(diǎn)相同，在控制體6中以d作射線與線段ef的交點(diǎn)記為g，將其類(lèi)比為梯形面計(jì)算，則可計(jì)算出梯形面積，從來(lái)可以得出控制體6內(nèi)的沉積體積；

⑤ 控制體1～5內(nèi)的顆粒沉積體積與控制體6中的計(jì)算方法相同。

圖1 等體積順序重構(gòu)法 Fig.1 Equal-volume sequence method

2.3 重構(gòu)方法的比較與選擇

順序重構(gòu)方法的優(yōu)點(diǎn)在于網(wǎng)格的重構(gòu)有著順序性，因此在編輯UDF的工作量上面要小很多，邏輯光面的思考量比較小，同時(shí)，這種順序方法重構(gòu)的網(wǎng)格對(duì)于負(fù)體積的接納能力要優(yōu)于法線重構(gòu)方法。但是它有著一個(gè)致命的缺陷，經(jīng)過(guò)計(jì)算發(fā)現(xiàn)，此種方法的計(jì)算速度慢，并且對(duì)處理器和內(nèi)存的要求高。

同順序重構(gòu)方法相比，法線重構(gòu)方法屬于非結(jié)構(gòu)方法，F(xiàn)luent本身就是非結(jié)構(gòu)計(jì)算軟件，不適用規(guī)定順序的結(jié)構(gòu)化計(jì)算，這種方法首先不需要尋找出現(xiàn)沉積量較少的網(wǎng)格，更不需要再去尋找臨近網(wǎng)格進(jìn)行計(jì)算，從而節(jié)省了給沉積網(wǎng)格排序的時(shí)間，在網(wǎng)格數(shù)量達(dá)到一定程度的情況下，這2種方法的計(jì)算結(jié)果是一致的，并且法線重構(gòu)方法比較節(jié)約時(shí)間和計(jì)算內(nèi)存。它也有自身的缺陷，即法線重構(gòu)法對(duì)于負(fù)體積的接納能力要差于順序重構(gòu)方法，不適合長(zhǎng)時(shí)間的或者濃度很大的顆粒沉積的數(shù)值模擬。

綜合起來(lái)，根據(jù)航空發(fā)動(dòng)機(jī)的工作狀況，顆粒直徑和顆粒濃度都非常小，而且每經(jīng)過(guò)1 000飛行小時(shí)基本上就會(huì)進(jìn)行洗發(fā)工作，因此適合使用法線重構(gòu)方法來(lái)模擬顆粒的沉積過(guò)程。

3 沉積過(guò)程的模擬與驗(yàn)證

3.1 模型與網(wǎng)格

首先，建立三維葉片模型，如圖2所示。進(jìn)行邊界條件設(shè)置，隨后設(shè)置并且編譯UDF進(jìn)入Fluent，進(jìn)行模擬并且與實(shí)驗(yàn)結(jié)果對(duì)比從而驗(yàn)證數(shù)值模擬和UDF的正確性。葉片弦長(zhǎng)為120 mm，葉高為30 mm，葉片材料使用鋁合金，其密度為 2 719 kg/m3，比熱容為871 J/(kg·K)，其物理屬性見(jiàn)文獻(xiàn)[19]。計(jì)算域包括葉柵通道，前后延長(zhǎng)大于一倍弦長(zhǎng)。

圖2 葉片三維模型圖 Fig.2 3D modeling schema of blade

由于葉片表面是個(gè)曲面，因此用近似六面體的網(wǎng)格，并不完全是六面體，在壁面邊界(左右壁面)第1層使用非結(jié)構(gòu)化網(wǎng)格，其余網(wǎng)格均使用六面體網(wǎng)格，主要保證葉片表面的網(wǎng)格正交性比較好就可以了，計(jì)算網(wǎng)格如圖3所示。

圖3 葉片計(jì)算網(wǎng)格 Fig.3 Computing grid for blade

網(wǎng)格的單元總數(shù)為170 880，為了延長(zhǎng)負(fù)體積的時(shí)間，可以進(jìn)一步進(jìn)行優(yōu)化，通過(guò)mesh/repair-improve/improve-quality指令進(jìn)行進(jìn)一步的優(yōu)化，此指令在優(yōu)化正交性的同時(shí)會(huì)降低扭曲度。

采用歐拉-拉格朗日方法進(jìn)行描述，選取離散相模型(Discrete Phase Model，DPM)，因?yàn)轭w粒的Stokes數(shù)變化范圍比較大，歐拉模型和混合物模型都只適用于Stokes數(shù)一定范圍內(nèi)的多相模擬。關(guān)于網(wǎng)格的重構(gòu)與融合部分，利用法線重構(gòu)法，在UDF里面對(duì)于變形進(jìn)行第一次光順，完成UDF輪詢(xún)之后，在時(shí)間步長(zhǎng)結(jié)尾對(duì)網(wǎng)格進(jìn)行第2次光順，此次光順通過(guò)求解擴(kuò)散方程來(lái)進(jìn)行光順。

顆粒選擇為Jim Bridger Power Station(JBPS)顆粒，其組成成分如表1所示，其顆粒的特性如表2所示。

表1 JBPS顆粒組成Table 1 JBPS composition

表2 JBPS顆粒特性Table 2 JBPS properties

3.2 軌跡模擬

在Fluent中，可以畫(huà)出顆粒的軌跡，在同樣進(jìn)出口溫度和進(jìn)出口壓比的條件下，分別取射入直徑d為1，10和100 μm的顆粒，觀察它們的曲線，如圖4所示。

圖4中，紅色軌跡代表100 μm顆粒，綠色代表10 μm顆粒，藍(lán)色代表1 μm顆粒，由圖中可以看出，隨著顆粒直徑的增大，在遇到葉片壓力面的情況下，主流體沿著壓力面表面流走，而顆粒卻由于慣性使得運(yùn)動(dòng)狀態(tài)不能及時(shí)改變，因此撞擊在壓力面上面。同時(shí)，由于主流體與顆粒之間有了相對(duì)速度，這種相對(duì)速度相對(duì)于顆粒來(lái)說(shuō)是向下的，根據(jù)動(dòng)量定理，此時(shí)顆粒會(huì)受到向下的作用力，因此，這也在一定程度上促使顆粒撞擊在葉片表面上。比如紅色顆粒由于直徑過(guò)大，因此撞擊區(qū)域在葉片前緣，進(jìn)行一次反彈后再次撞擊到相鄰葉片的前緣和中部，再次進(jìn)行反射，由于葉片幾何構(gòu)型的原因，才使得其順利流走，同時(shí)，10 μm的綠色顆粒同樣撞擊在葉片尾緣區(qū)域，而1 μm的藍(lán)色顆粒基本上就跟隨著主流運(yùn)動(dòng)，不與葉片進(jìn)行接觸。隨著顆粒直徑的減小，慣性力的作用就會(huì)越來(lái)越小，顆粒傾向于與主流體一起順著壓力面流走，因此，小直徑的顆粒更不容易撞擊在葉片上，連葉片都接觸不到，沉積下來(lái)的小直徑顆粒自然更低。

圖4 顆粒軌跡示意圖 Fig.4 Schematic of particle tracking

3.3 網(wǎng)格變形模擬

進(jìn)口壓強(qiáng)約選為1.8 MPa，出口為1 MPa，進(jìn)口溫度選為1 650 K，顆粒的入射速度為81 m/s，沉積模型選取臨界黏度模型，A與B如表2所示，模擬沉積生長(zhǎng)情況從而得到葉片的變形情況，如圖5所示。

從圖5中可以看出，編譯過(guò)UDF之后的Fluent可以模擬在沉積情況下的葉片變形情況，前緣的沉積大大增加了葉片的粗糙度，使得葉片形成鋸齒狀，而壓力面中部均勻沉積，因此粗糙度改變并不大，只是增加了厚度。這種現(xiàn)象與葉片的幾何構(gòu)型和溫場(chǎng)分布有關(guān)，在前緣溫場(chǎng)梯度較密，因此沉積量不均勻，而到了葉片壓力面中部，溫度場(chǎng)分布較為均勻，葉片開(kāi)始變得接近平板，因此更接近均勻沉積。

圖5 葉片變形情況 Fig.5 Deformation of blade

3.4 葉片沉積生長(zhǎng)模擬與驗(yàn)證

進(jìn)口壓強(qiáng)約選為1.8 MPa，出口為1 MPa，進(jìn)口溫度選為1 650 K，顆粒的入射速度為81 m/s，沉積模型選取臨界黏度模型，A與B如表2所示，圖6和圖7表示壓力面顆粒的沉積變形量與沉積質(zhì)量。

由圖6與圖7可以看出，葉片發(fā)生了明顯變形，沉積變形量分布基本上與沉積質(zhì)量保持一致，對(duì)于葉片壓力面來(lái)說(shuō)，沉積使得葉片呈現(xiàn)橫向條紋式的變形，在壓力面中部，也就是直徑10～100 μm 顆粒的主要沉積區(qū)域，沉積使得縱向呈條紋式的變形，并且其隨著時(shí)間的推進(jìn)程度呈現(xiàn)拋物線形式從葉根與葉冠向著葉片中部發(fā)展。

圖7中，對(duì)于直徑在10～100 μm之間的顆粒，其沉積主要從葉片壓力面的中部進(jìn)行生長(zhǎng)，逐漸向葉片前緣和尾緣生長(zhǎng)，并且由于顆粒直徑的原因，前緣和尾緣很難造成大面積的沉積，這是因?yàn)樵谖簿墔^(qū)域葉片壓力面法線與顆粒速度方向角度變小，使得大直徑的顆粒更加容易撞擊進(jìn)而反彈，前緣法線幾乎與顆粒入射速度角度重合，因此更加容易造成顆粒的反彈，這是由于顆粒的直徑變大，慣性力和動(dòng)量也會(huì)增大，其黏性不足以平衡反作用力，從而導(dǎo)致顆粒無(wú)法沉積，只能反射。

圖6 葉片變形量 Fig.6 Blade deformation

圖7 沉積質(zhì)量 Fig.7 Deposition mass

圖8 實(shí)驗(yàn)結(jié)果與數(shù)值模擬對(duì)比 Fig.8 Comparison of experiment and CFD simulation

圖8中的沉積實(shí)驗(yàn)是在TADF裝置中進(jìn)行的快速沉積實(shí)驗(yàn)，實(shí)驗(yàn)葉片為CFM56系列高壓渦輪葉片，其入射顆粒的流量是200 μg/s，基本上接近于北京地面的PM2.5的濃度，入射顆粒為JBPS顆粒，入射粒徑估算為1～30 μm，進(jìn)口壓強(qiáng)為1.8 MPa，出口壓強(qiáng)為1 MPa，入射溫度為1 500 K，作用3 h。

實(shí)驗(yàn)與數(shù)值研究所選用的沉積物顆粒組分完全相同，但是其質(zhì)量流量增大了萬(wàn)倍；另一方面，由于真實(shí)的入射顆粒直徑無(wú)法控制，因此只能通過(guò)統(tǒng)計(jì)學(xué)理論來(lái)估算平均顆粒直徑為10 μm，其真實(shí)入射顆粒直徑為正態(tài)分布，期望值為10 μm。

至于沉積量的測(cè)量方面，在實(shí)驗(yàn)中測(cè)量沉積量和沉積速度是非常難辦的，往往采用估算的方法，誤差非常大，因此主要對(duì)比實(shí)驗(yàn)與數(shù)值研究的沉積分布，只要分布一致，即以相同條件下數(shù)值研究量化分析來(lái)彌補(bǔ)實(shí)驗(yàn)不足。

從圖8可以看出，數(shù)值模擬和實(shí)驗(yàn)結(jié)果[20]的沉積分布量非常一致，均分布在葉片壓力面中部(數(shù)值研究沉積程度為：藍(lán)<綠<紅)，其次是前緣，尾緣沉積非常少。并且實(shí)驗(yàn)與數(shù)值研究壓力面中部的沉積變形均為縱向擴(kuò)展，而前緣部分變形為橫向鋸齒。在前緣變形的體現(xiàn)上面，實(shí)驗(yàn)與數(shù)值計(jì)算均不能明顯體現(xiàn)出變形，這說(shuō)明前緣沉積量至少與壓力面中部的沉積量積累差一個(gè)等級(jí)。

數(shù)值計(jì)算相對(duì)于實(shí)驗(yàn)的優(yōu)勢(shì)即在于可以完美控制顆粒直徑，并不需要利用統(tǒng)計(jì)學(xué)估計(jì)出平均直徑；另外，數(shù)值研究可以無(wú)休止地計(jì)算下去，而不僅僅是實(shí)驗(yàn)所能進(jìn)行的幾小時(shí)的沉積實(shí)驗(yàn)，這避免了實(shí)驗(yàn)單方面增大太多濃度而造成的誤差。

圖8驗(yàn)證了在氣體剛剛從燃燒室進(jìn)入渦輪中，其氣動(dòng)速度比較低，并沒(méi)有經(jīng)過(guò)多級(jí)渦輪加速，并且同時(shí)發(fā)動(dòng)機(jī)內(nèi)部的雜質(zhì)顆粒直徑比較小，根據(jù)Tafti[15]的研究，在較小顆粒直徑的沉積模擬方面，臨界速度模型與臨界黏度模型的沉積效率相差不多，因此為研究燃燒室沉積而利用的臨界黏度模型完全可以勝任前幾級(jí)渦輪的沉積生長(zhǎng)與分布數(shù)值研究，并且與真實(shí)的沉積情況有著驚人的一致性。

3.5 腐蝕預(yù)測(cè)

首先，由圖9可以看出，前緣的沉積呈現(xiàn)鋸齒狀，這樣會(huì)大大增大葉片前緣的粗糙度，從而增大顆粒的碰撞概率，最終加劇顆粒的磨蝕作用，使得葉片前緣最先遭受腐蝕。

在圖8實(shí)驗(yàn)的基礎(chǔ)上，利用同樣的葉片進(jìn)行長(zhǎng)時(shí)間的沉積實(shí)驗(yàn)，將顆粒入射流量提高到200 mg/s，

溫度提高至1 700 K，其條件仍然與圖8相同，提高顆粒流量與主流溫度可以使得顆粒盡快沉積，造成腐蝕的時(shí)間縮短，實(shí)驗(yàn)6 h之后的結(jié)果如圖10 所示。

圖9 葉片前緣變形 Fig.9 Deformation of blade leading edge

圖10 沉積分布與腐蝕實(shí)驗(yàn)結(jié)果 Fig.10 Experimenal results of deposition and errosion

下面將葉片沉積的模擬與葉片腐蝕的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比：

從圖9中可以看出，X和Y分別為基準(zhǔn)坐標(biāo)軸，沉積會(huì)使得葉片表面尤其是前緣粗糙度增加，而腐蝕最嚴(yán)重的區(qū)域恰恰是沉積變形量分布最不均勻的區(qū)域，即粗糙度最大的區(qū)域葉片前緣，這是與腐蝕的形成機(jī)理有著重要關(guān)系的。正是因?yàn)榇植诙鹊脑黾樱沟酶g速度大大加快，因此前緣和尾緣角落處最先受到腐蝕。腐蝕形成的原因是多方面的，其中比較重要的兩個(gè)方面是熔融顆粒物沉積并且滲透到葉片涂層中和顆粒的磨蝕撞擊，這2個(gè)腐蝕的主要形成原因造成了這種腐蝕與多步長(zhǎng)法計(jì)算的沉積分布有著驚人相似這種結(jié)果。

從沉積變形圖與腐蝕實(shí)驗(yàn)圖的對(duì)比來(lái)看，沉積不均勻的區(qū)域預(yù)示著腐蝕嚴(yán)重的區(qū)域，因此，沉積變形可以預(yù)測(cè)腐蝕的區(qū)域。

4 顆粒沉積特性研究

4.1 粒徑對(duì)于沉積的影響

顆粒直徑影響粒子的Stokes數(shù)，粒子的直徑越大，其對(duì)應(yīng)的Stokes數(shù)也越大，而Stokes數(shù)是影響沉積分布與沉積速度的重要因素。其余邊界條件同3.4節(jié)相同，改變粒徑，從而研究顆粒直徑對(duì)于沉積分布與沉積生長(zhǎng)速率的影響。對(duì)于各直徑區(qū)間內(nèi)部顆粒在葉柵通道內(nèi)的流動(dòng)特性3.2節(jié)已經(jīng)進(jìn)行了分析，選用直徑1～10 μm與直徑10～100 μm的顆粒進(jìn)行數(shù)值模擬。

圖11表示進(jìn)行360 s沉積之后，直徑為10 μm 的顆粒與直徑為100 μm的顆粒沉積情況的對(duì)比，橫坐標(biāo)表示從壓力面前緣至尾緣(0～1)，從圖11中可以看出，小直徑的沉積分布更傾向于葉片尾緣，而大直徑顆粒的沉積分布有向葉片前緣移動(dòng)的趨勢(shì)，在直徑增大10倍的情況下已經(jīng)移動(dòng)到了葉片中部。如果使用更大直徑的顆粒進(jìn)行沉積生長(zhǎng)的數(shù)值與實(shí)驗(yàn)研究會(huì)發(fā)現(xiàn)，隨著顆粒直徑的增大，沉積區(qū)域會(huì)從葉片尾緣向葉片前緣移動(dòng)，這正是顆粒的運(yùn)動(dòng)機(jī)理造成的結(jié)果。

圖12還表明了顆粒直徑對(duì)于沉積速度也是有很大影響的，同樣進(jìn)行360 s的沉積，10 μm的顆粒沉積速度為5×10-10kg/s，而100 μm顆粒沉積速度為2×10-6kg/s，并且100 μm顆粒沉積密度遠(yuǎn)大于10 μm顆粒。從1～100 μm選取多組不同粒徑進(jìn)行模擬。

以Stokes數(shù)為100的顆粒沉積速度為標(biāo)準(zhǔn)進(jìn)行無(wú)量綱化，即可得到一個(gè)無(wú)量綱數(shù)，表征為Stokes數(shù)小于100的顆粒沉積效率，發(fā)動(dòng)機(jī)的內(nèi)的顆粒Stokes一般不超過(guò)10。從圖12可以看出，直徑的增大會(huì)使得沉積速度增大，從而導(dǎo)致沉積效率增加，當(dāng)顆粒直徑增大到一定程度，由于慣性力的作用，顆粒會(huì)100%沉積下來(lái)，達(dá)到最大沉積速度，之后不再改變。

圖11 粒徑對(duì)沉積的影響 Fig.11 Effects of diameter on deposition

圖12 Stocks數(shù)對(duì)沉積效率的影響 Fig.12 Effects of Stocks number on deposition efficiency

4.2 主流溫度對(duì)沉積的影響

保持其余邊界條件不變，與3.4節(jié)相同，改變主流溫度，從而研究主流溫度對(duì)于沉積的影響，將主流溫度1 200 K提高1 600 K，每隔100 K進(jìn)行一次研究，將各個(gè)溫度點(diǎn)對(duì)應(yīng)的沉積速度均除以1 800 K的沉積速度得到無(wú)量綱沉積速度，得到如下結(jié)果：

主流的改變不影響沉積的分布，顆粒沉積的分布如圖13所示，其變化的只是沉積速度的大小，換而言之，主流溫度只影響沉積的生長(zhǎng)速度，主流溫度在一定范圍內(nèi)升高，沉積速度加快；根據(jù)顆粒的性質(zhì)以及入射速度的不同，存在一個(gè)臨界主流溫度，超過(guò)臨界溫度之后主流溫度的提高對(duì)沉積無(wú)影響。

圖13 主流溫度對(duì)沉積的影響 Fig.13 Effects of mainflow temperature on deposition

5 結(jié) 論

1) 基于臨界黏度模型與法線網(wǎng)格重構(gòu)法編寫(xiě)了顆粒沉積計(jì)算與葉片變形方面的程序模塊，將其用UDF嵌入Fluent數(shù)值模擬，并且與國(guó)外已有實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證本文計(jì)算方法的合理性。

2) 由沉積使得葉片變形的結(jié)果來(lái)預(yù)測(cè)了葉片腐蝕的情況，同時(shí)進(jìn)行關(guān)聯(lián)原因分析，并且與國(guó)外已有腐蝕方面的實(shí)驗(yàn)進(jìn)行對(duì)比，從而驗(yàn)證了預(yù)測(cè)方案的可行性。

3) 顆粒直徑既影響顆粒分布區(qū)域也影響顆粒沉積速度，這是由于顆粒直徑影響其受力狀態(tài)，慣性力在顆粒沉積方面有著很大影響。

4) 主流溫度對(duì)顆粒沉積分布無(wú)影響，與沉積速度正相關(guān)，其通過(guò)改變顆粒的物性從而對(duì)沉積速度與效率進(jìn)行影響。

[1] DONKELAAR A, MARTIN R, BRAUER M, et al. Global estimates of ambient fine particulate matter concentrations from satellite-based aerosol optical depth: Development and application[J]. Environmental Health Perspectives， 2010, 118(6): 847-855.

[2] BRUN K, NORED M， KURZ R. Particle transport analysis of sand ingestion in gas turbine engines[J]. Journal of Engineering for Gas Turbines and Power, 2012, 134(1): 012402.

[3] DUNN M G, PADOVA J, MOLLER J E， et al. Performance deterioration of a turbofan and a turbojet engine upon exposure to a dust environment[J]. Journal of Engineering for Gas Turbines and Power, 1987, 109(3): 336-343.

[4] KIM J, DUNN M G， BARAN A J， et al. Deposition of volcanic materials in the hot sections of two gas turbine engines[J]. Journal of Engineering for Gas Turbines and Power, 1993, 115(3): 641-651.

[5] SCHNEIDER O, DOHMEN H J, BENRA F K， et al. Investigations of dust separation in the internal cooling air system of gas turbines： GT-2003-38293[R]. New York: ASME, 2003.

[6] EKKAD S, HAN J C. Detailed heat transfer distributions on a cylindrical model with simulated TBC spallation: AIAA-1997-1595[R]. Reston: AIAA, 1997.

[7] DAVID G, KAREN A. Simulating particle deposition and mitigating deposition degradation effects in film cooled turbine sections: DE-FC21-92MC29061[R]. Washington, D.C.: UTSR, 2010.

[8] WRIGHT I G, LEVENS C, PINT B A. An analysis of the potential for deposition, erosion, or corrosion in gas turbines fueled by the products of biomass gasification or combustion: GT-2000-0019[R]. New York: ASME, 2000.

[9] AHLUWALIA R K, IM K H, WENGLARZ R A. Flyash adhesion in simulated coal-fired gas turbine environment[J]. Journal of Engineering for Gas Turbines and Power, 1989, 111(4): 672-678.

[10] CROSBY J M, LEWIS S, BONS J P, et al. Effects of temperature and particle size on deposition in land based turbines[J]. Journal of Engineering for Gas Turbines and Power, 2008, 130(5): 819-825.

[11] 王德全, 夏智勛, 胡建新, 等. 固沖發(fā)動(dòng)機(jī)沉積數(shù)值模擬與試驗(yàn)研究[J]. 固體火箭技術(shù)， 2009， 32(1)： 38-42.

WANG D Q, XIA Z X, HU J X, et al. Numerical simulation and experimental study of deposition for ducted rockets[J]. Journal of Solid Rocket Technology, 2009, 32(1): 38-42 (in Chinese).

[12] 張大林, 陳維建. 飛機(jī)機(jī)翼表面霜狀冰結(jié)冰過(guò)程的數(shù)值模擬[J]. 航空動(dòng)力學(xué)報(bào), 2004, 19(1): 138-141.

ZHANG D L, CHEN W J. Numerical simulation of the rime ice accretion process on airfoil[J]. Journal of Aerospace Power, 2004, 19(1): 138-141 (in Chinese).

[13] 陳維建, 張大林. 飛機(jī)機(jī)翼結(jié)冰過(guò)程的數(shù)值模擬[J]. 航空動(dòng)力學(xué)報(bào), 2005, 20(6): 1010-1017.

CHEN W J, ZHANG D L. Numerical simulation of ice accretion on airfoils[J]. Journal of Aerospace Power, 2005, 20(6): 1010-1017 (in Chinese).

[14] 孫志國(guó), 朱程香, 付斌, 等. 結(jié)冰計(jì)算中翼面網(wǎng)格重構(gòu)研究[J]. 航空學(xué)報(bào), 2011, 32(2): 231-241.

SUN Z G, ZHU C X, FU B,et al. Research of grid reconstruction on wing surface for icing calculation[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(2): 231-241 (in Chinese).

[15] TAFTI D K, SREEDHARAN S S. Composition dependent model for the prediction of syngas ash deposition with the application to a leading edge turbine vane, GT-20110-23655[C]//ASME Turbo Expo 2010: Power for Land, Sea, and Air. New York: ASME, 2010: 615-626.

[16] SENIOR C J, SRINIVASACHAR S. Viscosity of ash particles in combustion systems for prediction of particle sticking[J]. Energy and Fuels, 1995, 9(2): 277-283.

[17] N’DALA I, CAMBIER F, ANSEAU M R, et al. Viscosity of liquid feldspars. Part I: Viscosity measurements[J]. Transactions & Journal of the British Ceramic Society, 1984, 83(4)： 105-107.

[18] SHIH T I P, BAILEY R T. GRID2D/3—A computer program for generating grid systems in complex-shaped two- and three-dimensional spatial domains: NASA-TM-102453[R]. Washington, D.C.: NASA,1990.

[19] Allegheny Ludlum Corporation Company. Stainless steel types 309 and types 310[S]. Pittsburgh, PA: Allegheny Lulum Corporation Company, 2002.

[20] WEBB J B, CASSADAY B. Coal ash deposition on nozzle guide vanes: Part II—Computational modeling[J]. Journal of Turbomachinery, 2011, 135(1): 1757-1767.

(責(zé)任編輯： 張晗)

URL：www.cnki.net/kcms/detail/11.1929.V.20161221.1520.004.html

Numericalsimulationofparticledepositionprocessinsideturbinecascade

YANGXiaojun1,2,*,ZHUJiaxiong2

1.CollegeofAeronauticalEngineering,CivilAviationUniversityofChina,Tianjin300300,China2.Sino-EuropeanInstituteofAviationEngineering,CivilAviationUniversityofChina,Tianjin300300,China

AnumericalsimulationisconductedtoobtainmoreaccuratedistributionofparticledepositiononturbinebladewithJimBridgerPowerStation(JBPS).C++isusedtocodesuitableUserDefinedFunction(UDF),whichreconstructsgridstoaccomdatecomplexboundarydeformationineverytimestep.Insuchawaythatdepositioncanchangegeometryfeaturesandheattransfercharacteristics,andthewholedepositionprocessonbladecanbeanalyzedandtheexactdepositiondistrbutioncanbeobtained.Afaircomparisonofnumericalresultswithexistingexperimentalresultsshowstherationalityandaccuracyofgridreconstructioncode.Onthebasisofdepositiondistributionandgeometrydeformationofturbineblade,theerosionareaonthebladeispredicted.Thereafter,theeffectofmainstreamtemperatureandparticlediameterisanalyzed.Theresearchresultsshowthatthemaindepositionoccursonthecentralareaofbladepressureside,whereasthereisobviousdeformationinleadingedgeandcentralarea.Theobviousdeformationrevealsthatroughnessisincreasedanderosionismoresevereonleadingedge.Thediameterofparticlehasimportantinfluenceonparticledepositionanddepositionrate,whereastemperatureofmainstreamonlyhasinfluenceondepositionrate.

deposition;deformation;UDF;erosion;turbineblade;particlediameter;numericalsimulation

2016-06-14;Revised2016-07-18;Accepted2016-08-15;Publishedonline2016-12-211520

s：JointFundsoftheNationalNaturalScienceFoundationofChinaandCivilAviationAdministrationofChina(U1633113);TianjinResearchProgramofApplicationandAdvancedTechnology(14JCQNJC06800);theFundamentalResearchFundsfortheCentralUniversities(ZXH2012H004)

.E-mailxiaojunyoung@hotmail.com

2016-06-14；退修日期2016-07-18；錄用日期2016-08-15； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-12-211520

www.cnki.net/kcms/detail/11.1929.V.20161221.1520.004.html

國(guó)家自然科學(xué)基金委員會(huì)與中國(guó)民用航空局聯(lián)合資助項(xiàng)目 (U1633113)； 天津市應(yīng)用基礎(chǔ)與前沿技術(shù)研究計(jì)劃 (14JCQNJC06800)； 中央高校基本科研業(yè)務(wù)費(fèi) (ZXH2012H004)

.E-mailxiaojunyoung@hotmail.com

楊曉軍， 祝佳雄. 渦輪葉柵通道內(nèi)顆粒物沉積過(guò)程的數(shù)值模擬J. 航空學(xué)報(bào),2017,38(5):120530.YANGXJ,ZHUJX.NumericalsimulationofparticledepositionprocessinsideturbinecascadeJ.ActaAeronauticaetAstronauticaSinica,2017,38(5):120530.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0237

V231.1

A

1000-6893(2017)05-120530-12