毫米級微型渦輪發動機性能仿真模型

劉傳凱, 李艷茹

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

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

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

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

1 微型葉輪的相似參數

對于常規尺度的發動機，壓氣機和渦輪(本節統稱為葉輪)與氣體工質的換熱熱流遠小于其軸功率，傳熱效應通常可以忽略，并且在大多數情況下，葉輪均可被認為處于雷諾自模區。因此，常規尺度葉輪的工作狀態可由相似轉速和相似流量2個參數唯一確定。

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

(1)

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

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

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

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

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

(2)

(3)

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

2 微型轉子的換熱模型

2.1 轉子結構傳熱分析

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

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

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

圖1 微型渦輪發動機的典型結構示意圖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。依據壓氣機和渦輪的熱平衡關系有

(4)

(5)

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

(6)

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

2.2 微型轉子共同工作方程

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

1) 轉子功率平衡殘量：

(7)

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

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

εT_W=WCC_out-WT_in

(8)

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

3) 尾噴管流量連續殘量：

εN_W=WT_out-WN_cr

(9)

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

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

(10)

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

(11)

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

3 微型渦輪發動機靜子換熱模型

3.1 靜子結構換熱路徑

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

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

3.2 靜子結構換熱模型

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

3.2.1 熱網絡基本元件

1) 固體熱節點

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

(12)

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

2) 氣體熱節點

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

(13)

εgas_Q=hgas_outWgas_out-hgas_inWgas_in-

(14)

εgas_P=Pgas_out-Pgas_in

(15)

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

3) 導熱熱阻元件

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

Qcond=k(Tsolid_1-Tsolid_2)

(16)

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

4) 對流換熱熱阻元件

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

Qconv=hconv(Tsolid-Tgas)

(17)

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

3.2.2 熱網絡的求解

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

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

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

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

4 仿真分析

4.1 仿真對象

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

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

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

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

4.2 仿真計算條件

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

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

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

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

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

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

4.3 結果分析

4.3.1 主要工作參數變化過程分析

1) 轉速變化過程

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

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

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

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

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

4.3.2 啟動過程工作線分析

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

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

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

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

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

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

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

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

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

4.3.3 與國外文獻的對比

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

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

5 結 論

毫米級微型渦輪發動機的低雷諾數效應和傳熱效應對整機工作特性產生巨大的影響。本文提出了微型渦輪發動機轉子共同工作方程，建立了微型渦輪發動機靜子結構的熱網絡計算方法，開發了微型渦輪發動機性能仿真計算軟件，正確模擬了上述影響因素對毫米級微型渦輪發動機工作特性的影響規律，為開展毫米級微型渦輪發動機的總體方案設計和技術研究奠定了技術基礎。基于本文的仿真結果，可以得到如下結論：

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

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

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