毫米级微型涡轮发动机性能仿真模型

刘传凯, 李艳茹

(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量级，毫米级微型涡轮发动机加速工作线仍比稳态工作线离喘振边界更远。零部件的非稳态传热导致微型涡轮发动机整个启动过程的工作线呈现特有的多拐点现象。