毫米級微型渦輪發(fā)動機性能仿真模型

劉傳凱, 李艷茹

(1. 北京航空航天大學 能源與動力工程學院 航空發(fā)動機氣動熱力國家級重點實驗室, 北京 100083；2. 先進航空發(fā)動機協(xié)同創(chuàng)新中心, 北京 100083； 3. 北京航空航天大學 交通科學與工程學院, 北京 100083)

基于微機電系統(tǒng)(MEMS)工藝的毫米級微型渦輪發(fā)動機是一種新型燃氣渦輪動力裝置。具有重量輕和功率密度高等優(yōu)點，在各種微型無人機、微小型單兵飛行器和高能量密度分布式能源等新興領域具有潛在的應用價值，近年來得到了國內外研究機構的普遍關注[1-2]。目前，國外的微型渦輪發(fā)動機的基礎研究和設計技術仍在不斷的探索中[3-4]，尚未發(fā)展成熟。而國內的相關研究才剛剛起步，缺乏相關技術基礎和設計經(jīng)驗。

由于微型渦輪發(fā)動機尺寸十分微小，其低雷諾數(shù)效應和傳熱效應都對整機工作特性產(chǎn)生極其顯著的影響[5-6]，基于發(fā)動機部件雷諾自模假設和絕熱假設的常規(guī)渦輪發(fā)動機性能仿真模型[7]無法準確模擬微型渦輪發(fā)動機的工作特性，甚至可能在微型渦輪發(fā)動機方案設計和性能分析中得到顛覆性的計算結果。因此，建立考慮低雷諾數(shù)效應和傳熱效應的微型渦輪發(fā)動機性能仿真模型，是開展微型渦輪發(fā)動機設計技術研究的必要前提。

本文基于對微型渦輪發(fā)動機工作原理的分析，提出了考慮低雷諾數(shù)效應和轉子傳熱效應的微型渦輪發(fā)動機共同工作方程，并采用局部集總參數(shù)假設建立了微型渦輪發(fā)動機結構件的熱網(wǎng)絡模型，基于面向對象技術開發(fā)了仿真計算程序，實現(xiàn)了微型渦輪發(fā)動機工作特性和部件傳熱的動態(tài)模擬，分析了微型渦輪發(fā)動機典型過渡態(tài)特性及其變化機理。研究結果可為深入理解微型渦輪發(fā)動機工作原理，開展相關基礎研究和方案設計提供參考。

1 微型葉輪的相似參數(shù)

對于常規(guī)尺度的發(fā)動機，壓氣機和渦輪(本節(jié)統(tǒng)稱為葉輪)與氣體工質的換熱熱流遠小于其軸功率，傳熱效應通?？梢院雎裕⑶以诖蠖鄶?shù)情況下，葉輪均可被認為處于雷諾自模區(qū)。因此，常規(guī)尺度葉輪的工作狀態(tài)可由相似轉速和相似流量2個參數(shù)唯一確定。

微型葉輪的葉片一般為硅材料刻蝕的二維葉形，葉片高度不超過0.2～0.4 mm，低雷諾數(shù)效應和傳熱效應的影響非常顯著。文獻[8]采用如式(1)的相似參數(shù)表示微型葉輪的特性：

(1)

為了便于微型渦輪發(fā)動機性能建模，在保證微型葉輪主要物理過程合理性的基礎上，建立以下近似條件或簡化假設：

1) 微型葉輪葉尖間隙的相對影響雖然比常規(guī)尺度發(fā)動機大很多，但葉尖間隙的絕對變化幅度(小于10 μm量級)仍顯著小于葉片高度(約200～400 μm)。葉尖間隙變化對氣動性能的影響遠不及傳熱效應顯著，近似認為葉輪在工作中仍然能保持幾何相似條件。

2) 微型葉輪盤直徑只有毫米量級，畢渥數(shù)非常小。因此輪盤內部溫差遠小于輪盤材料與氣體工質之間的溫差，近似認為輪盤表面為等溫面。

3) 由于氣體工質與葉輪機匣之間的換熱對氣動特性影響相對較小，在葉輪特性分析時不考慮機匣換熱的影響，只考慮轉子換熱的影響。而機匣換熱將作為靜子部件換熱的一部分，在靜子結構的熱網(wǎng)絡模型計算中予以考慮。

基于上述假設，并依據(jù)雷諾數(shù)和努塞爾數(shù)的定義，葉輪相似特性可表示為

(2)

(3)

基于上述通用相似特性的表述形式，可采用部件試驗或CFD方法得到微型轉子的特性，為微型渦輪發(fā)動機的性能建模提供特性數(shù)據(jù)。

2 微型轉子的換熱模型

2.1 轉子結構傳熱分析

微型渦輪發(fā)動機內部輻射換熱的熱流相對較小[8]，在本文的轉子結構和靜子結構的換熱分析中僅考慮固體導熱和對流換熱2種換熱形式。

微型渦輪發(fā)動機的典型結構如圖 1所示，轉子構件由離心壓氣機、向心渦輪及轉子連接軸組成，采用空氣軸承支撐。靜子構件采用機匣、燃燒室及支撐結構的一體化設計，內部設有空氣軸承氣體流路。轉子和靜子結構均采用單晶硅作為主要材料，通過MEMS加工工藝制造而成；燃料通常采用氫氣或甲烷等氣體燃料。

轉子換熱路徑如圖2所示，由于壓氣機和渦輪之間為封閉腔，對流換熱溫差和熱流均相對較小，因此忽略轉子軸向封閉腔的散熱。

圖1 微型渦輪發(fā)動機的典型結構示意圖Fig.1 Schematic of typical structure of micro turbine engine

圖2 轉子換熱路徑示意圖Fig.2 Schematic of heat transfer routes of rotor

設渦輪盤面與主流燃氣的換熱熱流為QFT；渦輪盤緣與空氣軸承氣流的換熱熱流為QBT；壓氣機盤緣與空氣軸承氣流的換熱熱流為QBC；壓氣機盤面與主流氣流的換熱熱流為QFC。設渦輪盤向轉子軸的導熱熱流為QST；轉子軸向壓氣機軸的導熱熱流為QSC。依據(jù)壓氣機和渦輪的熱平衡關系有

(4)

(5)

考慮到轉子軸的質量遠小于壓氣機和渦輪盤的質量，轉子軸的熱容可忽略。壓氣機、渦輪與轉子軸之間的導熱熱流也可由式(6)近似計算：

(6)

式中：λS為轉子軸的導熱系數(shù)；AS為轉子軸的橫截面積；lS為轉子軸的長度。

2.2 微型轉子共同工作方程

微型渦輪發(fā)動機性能仿真模型的示意圖如圖3所示。該性能仿真模型不僅需要考慮傳統(tǒng)的大型渦輪發(fā)動機部件法性能仿真模型[7]的流量連續(xù)和功平衡等共同工作條件，還需要考慮轉子的熱平衡。

1) 轉子功率平衡殘量：

(7)

式中：LC和LT分別為壓氣機功率和渦輪功率，根據(jù)式(3)插值計算得到；ηm為轉子的機械效率；J為轉子的轉動慣量。

2) 燃燒室與渦輪流量連續(xù)殘量：

εT_W=WCC_out-WT_in

(8)

式中：WCC_out為燃燒室出口燃氣流量；WT_in為渦輪入口燃氣流量，根據(jù)式(3)由渦輪特性插值計算得到。

3) 尾噴管流量連續(xù)殘量：

εN_W=WT_out-WN_cr

(9)

式中：WT_out為渦輪出口燃氣流量；WN_cr為依據(jù)尾噴口喉道面積計算的燃氣流量。

4) 壓氣機盤與轉子軸的導熱平衡：

(10)

5) 渦輪盤與轉子軸的導熱平衡：

(11)

圖3 考慮轉子熱平衡的微型渦輪發(fā)動機性能仿真模型示意圖Fig.3 Schematic of micro turbine engine performance simulation model considering thermal balance of rotor

3 微型渦輪發(fā)動機靜子換熱模型

3.1 靜子結構換熱路徑

靜子換熱的主要路徑如圖4所示，其對流換熱包括流道內強制對流和自然對流2種，后者僅發(fā)生于發(fā)動機殼體外壁面，在換熱熱流總量中占比較小。強制對流換熱主要分為2部分。一部分是主流道的氣體與流道壁面的換熱，氣體進入燃燒室之前溫度較低，氣流被流道內壁面加熱；而燃燒室出口氣流溫度遠高于壁面溫度，熱流由氣體流向流道壁面。另一部分發(fā)生于空氣軸承氣流與其流道壁面?？諝廨S承氣體為外部氣源提供的高壓冷空氣，氣體通道全部位于支撐結構內部，熱流方向顯然是由流道壁面指向空氣軸承氣流，故在圖中沒有標注。

圖4 靜子換熱路徑示意圖Fig.4 Schematic of heat transfer routes of stator

3.2 靜子結構換熱模型

如圖5和圖6所示，基于局部集總參數(shù)法的思想，將靜子結構抽象為由兩類熱節(jié)點和兩類熱阻元件組成的熱網(wǎng)絡模型，并與主流道性能模型進行耦合求解。

3.2.1 熱網(wǎng)絡基本元件

1) 固體熱節(jié)點

可集總的固體域被抽象為固體熱節(jié)點。在計算中，試給各熱節(jié)點平均溫度Tsolid，根據(jù)能量守恒建立殘量方程：

(12)

式中：csolid和msolid分別為固體節(jié)點的質量比熱容和質量；Qconv和Qcond分別為對流換熱熱流和固體導熱熱流，吸熱為正，放熱為負。

2) 氣體熱節(jié)點

氣體熱節(jié)點為主流道或空氣軸承中可集總的氣體域。試給各節(jié)點溫度Tgas、質量流量Wgas和壓力Pgas，則根據(jù)質量連續(xù)、能量守恒和壓力平衡，可建立3個殘量方程：

(13)

εgas_Q=hgas_outWgas_out-hgas_inWgas_in-

(14)

εgas_P=Pgas_out-Pgas_in

(15)

式中：W、ρgas、Vgas、P、hgas和u分別為氣體熱節(jié)點的質量流量、密度、體積、壓強、焓和熱力學能；下標gas_in表示進口氣流參數(shù)，gas_out表示出口氣流參數(shù)。

3) 導熱熱阻元件

固體導熱熱阻元件連接發(fā)生換熱的2個固體熱節(jié)點元件，其導熱熱流由式(16)計算：

Qcond=k(Tsolid_1-Tsolid_2)

(16)

式中：Tsolid_1和Tsolid_2分別為相鄰2個固體熱節(jié)點的節(jié)點溫度；k為傳熱系數(shù)，與換熱單元的幾何結構、材料導熱系數(shù)等因素有關[9]。

4) 對流換熱熱阻元件

對流換熱熱阻元件連接一個固體熱節(jié)點元件和一個氣體熱節(jié)點元件，其對流換熱熱流由式(17)計算：

Qconv=hconv(Tsolid-Tgas)

(17)

式中：hconv為對流換熱系數(shù)，由換熱準則式計算確定[8-10]。

3.2.2 熱網(wǎng)絡的求解

以燃燒室結構熱網(wǎng)絡模型為例, 根據(jù)其幾何結構和換熱形式，將其劃分為圖5所示的換熱單元，進而建立圖6所示的熱網(wǎng)絡模型。

圖5 燃燒室結構換熱單元劃分Fig.5 Heat transfer element division of combustor structure

圖6 燃燒室結構的熱網(wǎng)絡模型Fig.6 Heat network model of combustor structure

由于熱阻元件中無試給參數(shù)，而熱節(jié)點元件中試給參數(shù)的數(shù)量與殘量方程的數(shù)量相同。因此由任意數(shù)量的熱阻元件和熱節(jié)點元件組合的熱網(wǎng)絡方程組與微型渦輪發(fā)動機轉子共同工作方程組相聯(lián)立，均可構成封閉可解的耦合方程組。

4 仿真分析

4.1 仿真對象

基于面向對象的程序設計技術和變物性計算方法[11]，在Visual Studio 2010 平臺下自主開發(fā)毫米級微型渦輪發(fā)動機性能仿真軟件。仿真對象為以甲烷為燃料的毫米級微型渦輪發(fā)動機。

由于結構尺寸的限制，微型渦輪發(fā)動機無法采用傳統(tǒng)啟動機進行啟動，風車啟動是最可行的啟動方法[12]。

風車啟動過程主要包括以下4個步驟：①接通空氣軸承氣源，使發(fā)動機轉子處于懸浮狀態(tài)，且能夠自由旋轉。②接通發(fā)動機進口高壓啟動氣源，吹動轉子旋轉。③轉子轉速達到給定的風車轉速時，保持發(fā)動機進口壓力不變，燃燒室點火，甲烷流量迅速增大至設計值并保持不變，轉子轉速迅速增加。④發(fā)動機逐步達到穩(wěn)定的工作狀態(tài)后，撤除高壓啟動氣源，發(fā)動機將工作于自維持狀態(tài)，啟動完畢。

本文對微型渦輪發(fā)動機從穩(wěn)定的風車狀態(tài)開始，到最終達到穩(wěn)定自維持狀態(tài)為止的整個過渡過程開展仿真計算，最小仿真時間步長為0.001 s。

4.2 仿真計算條件

假定發(fā)動機工作于標準大氣條件下，發(fā)動機燃氣流量為被控量，進口總壓(啟動氣源總壓)為被控環(huán)境參數(shù)，啟動過程具體計算參數(shù)如下：

1)t=0 s時，進口總壓為115 kPa，甲烷流量為0。發(fā)動機處于穩(wěn)定的風車狀態(tài)，轉子轉速約為4×105r/min。

2)t=0～0.2 s為點火過程。進口總壓保持115 kPa不變，燃燒室點火并將甲烷流量迅速增加到設計值。

3)t=0.2～100 s為熱穩(wěn)定過程。甲烷流量保持不變，發(fā)動機內部經(jīng)歷復雜的非穩(wěn)態(tài)換熱過程后，逐步達到平衡狀態(tài)。

4)t=100 s時，撤去進口高壓，進口總壓由115 kPa突降至101 kPa。

5)t=100～200 s為自維持穩(wěn)定過程。甲烷流量保持不變，發(fā)動機達到最終的穩(wěn)定工作狀態(tài)。

4.3 結果分析

4.3.1 主要工作參數(shù)變化過程分析

1) 轉速變化過程

圖7給出了啟動過程中微型渦輪發(fā)動機轉速變化過程的模擬結果。由于轉子的轉動慣量非常小，在最初的0.2 s內，轉速即可快速增加至1.23×106r/min，其相對增長過程幾乎能夠與甲烷流量的相對增長過程保持同步。此后，雖然甲烷流量始終穩(wěn)定不變，但隨著零部件平均溫度的持續(xù)升高和總體吸熱熱流的下降，轉速仍會繼續(xù)緩慢增加。50 s以后，轉速能基本穩(wěn)定在1.65×106r/min的設計轉速附近。100 s時撤去高壓氣源，轉速出現(xiàn)了小幅度的突降，之后緩慢回升至設計轉速。

2) 換熱熱流及渦輪軸功率的變化過程

如圖8(a)所示，劇烈的非穩(wěn)態(tài)傳熱效應對微型渦輪發(fā)動機過渡態(tài)性能產(chǎn)生重要的影響。當發(fā)動機剛點火時(t<0.2 s)，整個發(fā)動機轉子仍處于相對較低的溫度，渦輪表面與工質相對溫差很大，渦輪表面的吸熱熱流很快達到約70 W的峰值，顯著高于渦輪輸出的軸功率。此后隨渦輪轉子溫度的快速升高，吸熱熱流在短時間內出現(xiàn)快速下降，而渦輪軸功率則隨靜子溫度升高而緩慢上升。在0.7 s時，渦輪軸功率才超過渦輪表面吸熱熱流。轉子冷端(壓氣機和空氣軸承)表面的換熱熱流則隨著其表面溫度的升高而快速上升。在1 s以后，轉子冷端向工質的總散熱熱流與渦輪從工質的吸熱熱流基本持平。由于靜子結構的總換熱面積和熱容都遠大于轉子，靜子從工質的總吸熱熱流在點火后快速達到220 W的峰值，此后隨靜子溫度的升高，吸熱熱流逐漸減少，50 s后基本達到穩(wěn)定狀態(tài)。第100 s時撤去高壓啟動氣源，氣體流量減小會導致燃燒室余氣系數(shù)下降和燃燒室出口溫度的小幅度突增，轉子與靜子的換熱量也隨之發(fā)生小幅度增加。經(jīng)過一定的時間后，發(fā)動機靜子與氣體工質的總換熱熱流最終穩(wěn)定在15 W左右，這與機匣結構在穩(wěn)定狀態(tài)通過自然對流向環(huán)境的散熱熱流持平。

圖7 轉速及甲烷流量的變化過程Fig.7 Evolution of rotational speed and methane flow

圖8 換熱熱流及軸功率的變化過程Fig.8 Evolution of heat flux and shaft power

4.3.2 啟動過程工作線分析

圖9給出了微型渦輪發(fā)動機啟動過程的工作線?？梢钥闯?，在0～0.2 s的加速過程中，微型渦輪發(fā)動機加速工作線基本上都在穩(wěn)態(tài)工作線(進口壓力為115 kPa時)之下，這與常規(guī)尺寸渦輪發(fā)動機加速工作線的軌跡完全不同。其根本的原因在于微型渦輪發(fā)動機轉子轉動慣量極小加速時間極短。加速過程中，各部件仍處于相對低溫狀態(tài)，整個流道的相對流通能力和壓氣機效率比對應轉速下的穩(wěn)態(tài)值更高。

圖9 微型渦輪發(fā)動機啟動過程共同工作線Fig.9 Startup operating line of micro-engine

在0.2 s時，發(fā)動機“供油(甲烷)量”已經(jīng)達到最高值并保持穩(wěn)定，但發(fā)動機轉子仍需一定的時間(約1 s)才能基本達到熱平衡。此時轉子溫度的快速增加不可避免地降低了壓氣機工作效率，發(fā)動機匹配的工作點因此快速向喘振線方向移動，導致啟動工作線出現(xiàn)第1個拐點。

1 s以后，轉子基本達到熱平衡，發(fā)動機工作點向喘振邊界移動速度變慢，靜子部件換熱成為影響發(fā)動機工作點的主導因素。隨著靜子結構溫度的逐漸提高，靜子結構從氣體工質吸熱的功率也在逐漸減少，發(fā)動機的轉速因此緩慢增加，工作點向右上方緩慢移動，導致啟動工作線出現(xiàn)第2個拐點。

100 s時，發(fā)動機高壓啟動氣源撤除，使得渦輪的膨脹比以及壓氣機的相對流通能力下降，發(fā)動機匹配的工作點再次向喘振邊界移動，導致啟動工作線出現(xiàn)第3個拐點。

綜上，微型渦輪發(fā)動機的微小尺度導致影響其加速特性的物理因素與常規(guī)尺度渦輪發(fā)動機存在顯著不同。

對于常規(guī)大尺度渦輪發(fā)動機，零部件的換熱熱流相對于渦輪功率的比例很小，一般可以忽略。發(fā)動機轉動慣量是影響其加速性能的主導因素，其加速過程的工作線會因此比穩(wěn)態(tài)工作線更靠近喘振邊界[13-14]。

而對于毫米級微型渦輪發(fā)動機，其零部件換熱熱流與渦輪軸功率處于同等量級。微型渦輪發(fā)動機的空氣流量大約為常規(guī)尺度發(fā)動機10-5量級，而轉子轉動慣量只有常規(guī)尺度發(fā)動機10-11量級。因此，在微型渦輪發(fā)動機加速過程中，轉動慣量的相對影響非常微小，而零部件的換熱成為影響發(fā)動機加速工作線的最主要因素。

不僅如此，微型渦輪發(fā)動機轉子總質量只有靜子的1/50左右。這使其在啟動的過程中，“油門”響應時間、轉子達到熱穩(wěn)定的時間、靜子達到熱穩(wěn)定的時間，分別處于10-1、100、102s量級。這是微型渦輪發(fā)動機啟動過程工作線出現(xiàn)若干拐點的根本原因。

4.3.3 與國外文獻的對比

由于微型渦輪發(fā)動機內部參數(shù)極難測量，也極少公開，本文采用文獻[15]的仿真結果作為對比數(shù)據(jù)。圖10給出了油門響應時間在0.02 s的條件下，轉子冷端散熱和熱端吸熱熱流的相對變化過程(以渦輪吸熱熱流在50 s時的穩(wěn)定值作為基準)，并與文獻[15]仿真結果進行對比?？梢钥闯?，本文仿真結果與文獻[15]的變化趨勢一致，能夠反映毫米級微型渦輪發(fā)動機啟動過程中的基本氣動和換熱過程。導致少量相對差異的可能原因是本文模擬的微型渦輪發(fā)動機的幾何結構與文獻[15]并不完全一致、葉輪的氣動結構和傳熱特性存在差異、燃料(文獻[15]為氫氣)存在差異。另外，在仿真算法上，文獻[15]將工質在葉輪內的傳熱和流動抽象為相對獨立的串聯(lián)過程，同時把整機靜子部件集總為單個節(jié)點，不分辨結構件內部溫度分布，這也是導致仿真結果出現(xiàn)相對差異的原因之一。理論上，本文的仿真算法比文獻[15]的算法的具有更高的空間分辨率，能夠更真實地模擬微型渦輪發(fā)動機內部換熱及其對整機性能的影響過程。

圖10 轉子換熱熱流相對值的變化過程Fig.10 Evolution of rotor relative heat flux

5 結 論

毫米級微型渦輪發(fā)動機的低雷諾數(shù)效應和傳熱效應對整機工作特性產(chǎn)生巨大的影響。本文提出了微型渦輪發(fā)動機轉子共同工作方程，建立了微型渦輪發(fā)動機靜子結構的熱網(wǎng)絡計算方法，開發(fā)了微型渦輪發(fā)動機性能仿真計算軟件，正確模擬了上述影響因素對毫米級微型渦輪發(fā)動機工作特性的影響規(guī)律，為開展毫米級微型渦輪發(fā)動機的總體方案設計和技術研究奠定了技術基礎?；诒疚牡姆抡娼Y果，可以得到如下結論：

1) 微型渦輪發(fā)動機轉動慣量對其加速性能的影響微小， 非穩(wěn)態(tài)傳熱效應是影響其過渡態(tài)特性的主要因素。在毫米級微型渦輪發(fā)動機的典型啟動過程中，轉子熱平衡時間要比“油門”響應時間多一個數(shù)量級以上，靜子熱平衡時間比轉子熱平衡時間多2個數(shù)量級。

2) 微型渦輪發(fā)動機渦輪軸功率、渦輪吸熱熱流、轉子冷端散熱熱流、整機靜子吸熱熱流均處于同等數(shù)量級。在點火啟動后的0.7 s內，渦輪軸功率均低于渦輪吸熱熱流。靜子部件的瞬時吸熱峰值功率可達最大渦輪軸功率4倍以上。

3) 即使“油門”響應時間只有10-1s量級，毫米級微型渦輪發(fā)動機加速工作線仍比穩(wěn)態(tài)工作線離喘振邊界更遠。零部件的非穩(wěn)態(tài)傳熱導致微型渦輪發(fā)動機整個啟動過程的工作線呈現(xiàn)特有的多拐點現(xiàn)象。