地铁车辆牵引逆变器热管散热器的温升试验及热仿真

丁　杰，张　平

(1.湘潭大学 土木工程与力学学院，湖南 湘潭　411105；2.南车株洲电力机车研究所有限公司 南车电气技术与材料工程研究院，湖南 株洲　412001)

牵引逆变器是地铁车辆的关键部件，而绝缘栅双极型晶体管(Insulated Gate Bipolar Transistor，IGBT)模块是牵引逆变器中最为核心的器件。根据IGBT模块的失效机理可知，其在封装时各层材料的热膨胀系数不一致，在长期高温热循环作用下可发生铝键合线断裂或脱落、硅芯片与衬板之间及衬板与基板之间的焊料层老化、栅极氧化层损坏和芯片失效等[1-3]。因此，IGBT模块的散热设计至关重要。

设计地铁车辆牵引逆变器时，大多数采用热管散热器加走行风冷的方式，这样可以省去风道和风机，简化牵引逆变器的柜体结构，减少风机故障隐患和振动噪声的影响[4-5]。但由于地铁线路站点多、站间距短，列车启停、加减速频繁，使牵引逆变器的电气参数不断变化，导致IGBT模块的损耗也实时改变；加之走行风速有快有慢，使热管散热器的散热能力也相应变化，因此仅计算热管散热器和IGBT模块某一工况的稳态温升是无法准确评估热管散热器的性能和IGBT模块的疲劳寿命的，有必要对其瞬态温升进行计算。

目前，研究热管散热器和IGBT模块温升的方法可分为试验分析法[6]和仿真分析法[7]。试验分析法可在实验室或实际运行的车辆中进行，优点是测试的温升数据可信度高，缺点在于无法直接测量封装于绝缘材料中芯片及焊料等部位的温度。仿真分析法大多基于电气、热等学科建立仿真模型并进行单独计算。电气仿真借助Matlab/Simulink，Saber和Simplorer等软件建立主电路仿真模型，得到IGBT模块的电压、电流等信息，进而计算芯片损耗。热仿真主要有热阻抗网络、有限元、有限体积和模型降阶等方法。热阻抗网络法将三维结构的导热物体映射成为热阻和热容串级连接的一维热阻抗网络模型，使用方便、计算量很小，然而结温、壳温和热管散热器温度是虚拟的理论平均值，与实际三维结构不同部位表现出明显的温度梯度有较大差异，且难以考虑同一热管散热器上不同IGBT模块之间以及同一IGBT模块中不同芯片之间的热耦合效应[7]。有限元法需要将IGBT模块和热管散热器的几何结构划分为有限元网格，在热管散热器与冷却介质的接触面设置等效对流换热系数，而且计算时用较小的时间步长才能得到较为准确的瞬态温度场计算结果，因此对计算机资源要求较高且求解速度慢。有限体积法是计算流体动力学(Computational Fluid Dynamic，CFD)中应用最为广泛的方法，需要将IGBT模块、热管散热器和冷却介质划分为有限体积网格，自动计算流固耦合面的对流换热系数，其计算结果比有限元法更符合实际情况，然而边界层对流动与换热的影响使得网格尺寸被限制为较小的值，导致网格数量大幅度增加至数千万，与此同时，流体区域的压力与速度耦合，需要消耗大量的计算机资源进行迭代，才能得到收敛的稳态计算结果；并且，计算瞬态问题时，每一时间子步的迭代计算需要消耗计算机资源且计算结果文件需要数GB的硬盘存储空间，因此，瞬态CFD计算的代价远远超出有限元法，极少应用于实际工程中的复杂对象。模型降阶法[8-11]采用Krylov子空间、平衡截断、正交分解、积分全等变换等算法得到阶次较低的近似模型，可以在较高准确性基础上极大提高计算效率。Simplorer和Caspoc等部分商业软件包含了模型降阶的功能，可以实现多个软件的联合仿真，由于商业软件为了保证计算效率，输出结果是多个探测点的变量变化曲线，很难输出和云图显示整个物理场的数据，因此不能完全满足实际的应用需求。

为快速、准确地分析地铁车辆牵引逆变器热管散热器的性能及其IGBT模块的温升，本文首先搭建单个变流模块的温升试验装置，通过试验得到可用于热管散热器热仿真计算的输入条件及验证仿真结果准确性的试验数据；然后采用FLUENT软件对热管散热器进行稳态CFD仿真，得到其温度场和对流换热系数的分布；再将有限元法与模型降阶法相结合，开发快速计算程序，最后通过Matlab/Simulink软件实现热与电气的联合仿真，以有效解决复杂工况的瞬态计算问题。

1　温升试验

1.1　试验装置

某地铁车辆牵引逆变器采用二电平电压型直—交逆变主电路，其主要参数为：额定电压DC750 V，额定输出容量600 kVA，输出电压0～585 V，输出频率0～160 Hz，开关频率500 Hz，额定输出电流600 A，牵引工况最大输出电流有效值884 A，制动工况最大输出电流有效值1 156 A，额定工作点效率0.98，采用直接转矩控制方式。当车辆处于牵引工况时，直流电经过高压电器箱、滤波电抗器箱分别对2个IBCM60G1型变流模块供电，经牵引逆变器输出三相变压变频的交流电向4台异步牵引电机供电。当车辆处于再生制动工况时，牵引逆变器将异步牵引电机输出的三相交流电变成直流电并反馈回电网，或由制动电阻消耗掉。三相交流电逆变及电阻制动消耗环节的开关管均为FZ1600R17KF6C-B2型IGBT模块。

考虑到IGBT模块的芯片被绝缘材料封装起来，无法对其直接进行温度测量，因此在热管散热器的安装面上布置PT100热电阻进行温度测量。热管散热器上测温点的布置情况如图1所示。图中：V2—V7和V1，V8分别为实现三相逆变和电阻制动功能的IGBT模块；J1—J3为温度继电器；1#—6#为相邻的IGBT模块间的测温点；7#为J2的测温点。

图1　测温点布置情况

在实验室条件下进行温升试验的基本步骤为：首先在变流模块组装之前于热管散热器上布置测温点，然后将其与总装后的牵引逆变器、通风机、风罩、电源和负载(牵引电机)等相连，再将测温点的数据线连接至温度巡检仪和计算机，最后进行试验测试和数据分析。

考虑到实验室条件下不能真实模拟走行风且只有1台牵引电机作为负载，因此重点研究在减小变流模块1的输入电流时热管散热器的温升情况。在变流模块1的翅片部分设计1个风罩，通过轴流通风机进行吹风，连接通风机的风罩部分为圆形，包裹翅片的风罩部分为矩形，风罩中间部分采用渐变的形式并设置整流隔栅，以使通风机的出风尽量均匀地吹向热管散热器。在风罩入口处布置1个热电偶，用于测量牵引逆变器附近的空气温度。搭建的试验装置如图2所示。

图2　试验装置

1.2　试验结果

在风罩入口和出口处分别测得9个点的风速值，如图3所示。由图3可以看出：与入口处其他部位相比，中心点处的风速因通风机的阻挡而偏小，入口上部的风速高于下部，而冷却空气经过通风机、风罩、热管散热器翅片及外罩后，出口上部的风速低于下部，出口处风速与入口处风速不对应的情况说明风罩内的空气流动是非常复杂的。

图3　风罩入口和出口处的风速测量结果(单位： m·s-1)

根据风罩入口处的风速测量结果，可以计算出平均风速约为5.7 m·s-1。

图4为由各测温点得到的温度变化曲线。由图4可以看出：试验开始后0～30 min内温度上升速度很快，30 min后上升速度逐渐变缓；风罩入口处空气的温度为31.1 ℃，持续试验220 min后升为33.8 ℃，这是由于试验场地的空间较为狭窄，变流模块1产生的热量难以有效耗散至实验室外部，且牵引逆变器周围还存在其他工作中的试验设备，导致风罩入口处的空气温度不断上升；随着试验时间的推移，各测温点的温度与风罩入口处空气温度的差值趋于平稳，说明热管散热器的温升达到了平衡状态，因此风罩入口处空气温度的变化在试验误差允许范围内；测温点中4#的温度最高、3#的温度次之，而6#的温度最低，这主要是由于V1和V8未工作且处于冷却空气上游、持续工作的V4和V5处于冷却空气最下游所致。

图4　温升试验得到的温度变化曲线

2　CFD稳态热仿真分析

2.1　CFD仿真模型

热管散热器内部涉及复杂的相变过程，目前还很难对其本身进行数值模拟。较为常用的方法是将热管假设为轴向导热系数很高、径向导热系数为热管材料导热系数的实体棒杆[4-6]。由于IGBT模块中包含AlSiC基板、焊料、铜层、AlN基片、芯片和绝缘材料，因此结合各层材料的厚度和整个模型的复杂程度，将仿真模型中网格的基本尺寸设置为2 mm，远离热管散热器的流体区域网格尺寸设置为5 mm。以风罩为参照，将被风罩包裹的变流模块1的流体外边界设置为壁面边界条件；变流模块2因为没有被外罩包裹，故将其流体区域适度扩大，并将其流体外边界设置为压力出口条件。利用HyperMesh软件建立六面体为主的高质量网格，数量为3 862万个，其中流体区域的网格数量为2 520万个。

为便于热仿真结果与温升试验数据的对比，以温升试验条件作为热仿真计算的输入条件。考虑到风罩内的空气流速分布极为不规则，为简化仿真条件的参数设置，将风罩入口处的风速设为5.7 m·s-1、空气温度设为33.8 ℃。假设冷却空气在风道内的流动状态为完全湍流，采用标准k—ε模型进行模拟。利用FLUENT软件和DELL T7600台式工作站对CFD仿真模型进行稳态计算，可以得到收敛的仿真结果。

2.2　稳态仿真结果

图5为冷却空气在试验装置内流动时不同截面处冷却空气的流速分布云图。由图5可以看出：冷却空气受热管散热器翅片、热管和外罩的阻挡，在不同部位处的流速分布并不相等，在局部截面面积小的区域其流速可达到最大值，为16.61 m·s-1。

图5　不同截面的冷却空气流速分布云图

图6为经过87 h仿真计算得到的变流模块1的热管散热器及IGBT模块温度场分布云图。由图6可以看出：因冷却空气主要沿z坐标轴正方向流动，处于冷却空气最下游的V4和V5这2个IGBT模块的温度最高，可达55.43 ℃；V1和V8模块的温度最低，这是由于它们不发热且位于冷却空气的上游；1#—7#测温点的温度分别为39.65，45.81，49.02，48.95，45.61，39.49和45.33 ℃，其中3#点的温度略高于4#点的，这与温升试验结果有区别的主要原因是试验中通风机入口处冷却空气的风速不均匀，而仿真时假设入口处的冷却空气风速完全相等；对比仿真结果与试验数据可知，各测温点温度的相对误差不超过5%，说明仿真结果的准确性。

图6　热管散热器及IGBT模块的温度场分布云图

图7为变流模块1的热管散热器对流换热系数分布云图。由图7可以看出：冷却空气流经散热片、热管、外罩和底板时，会产生不同程度的对流换热作用，从而表现出对流换热系数值的差异，即对流换热系数的分布并不一致；对流换热系数值大的区域表示对流换热较为充分，反之则不充分。

图7　热管散热器的对流换热系数分布云图

对流换热系数分布处于CFD模型中流体和固体相接触的流固耦合面上，流体与固体之间的热量通过这些面进行传递，因此在有限元模型只需考虑固体和流固耦合面的换热，则对流换热系数分布是有限元模型中非常关键的边界条件。在大多数有限元法和模型降阶法的工程应用中，将边界条件即对流换热系数设为常数，这种处理方式与实际情况有差异，这也是有限元法的计算精度低于有限体积法的根本原因。

3　快速瞬态热仿真分析

3.1　理论基础

图6所示的温度场分布是经过较长计算时间得到的稳态结果，而图4所示的温度变化曲线是通过温升试验得到的瞬态数据，故为了更好地实现仿真结果与温升试验数据的对比，需要进行瞬态热仿真。考虑到CFD仿真的准确性高但计算效率非常低，故在满足计算精度的前提下为进一步实现热管散热器和IGBT模块瞬态导热问题的快速计算，可将图7所示的对流换热系数分布插值至有限元模型中，以保证有限元模型与CFD仿真模型的边界条件一致，再结合有限元法及模型降阶法中最基本和最重要的Krylov子空间法[8-11]进行程序开发。

利用有限元法计算热管散热器和IGBT模块等固体区域的温度场时，瞬态导热问题的微分方程为

(1)

式中：ρ为密度；cp为定压比热容；T为温度；t为时间；λx，λy和λz为沿物体3个主方向的导热系数；q为物体内部的热流密度。

有限元模型的边界条件分为给定温度、给定热流密度和给定对流换热系数3种。由微分方程等效积分形式的伽辽金提法，在空间域离散后得到包含n个有限元模型节点数目的1阶常微分方程组矩阵方程如下。

(2)

式中：Cp为热容矩阵；K为热传导矩阵；P为温度载荷列阵；φ为节点温度列阵。

式(2)涉及到大规模矩阵运算，消耗的计算机资源较大，Krylov子空间法可以将物理模型的大规模状态空间投影到以1组基矢量为特征的低维空间(降阶阶数r通常取5～50，主要取决于时间步长[12]，该数值远小于原有阶数n)，从而有效降低矩阵的维度，并在保证计算结果准确性的基础上大幅度降低计算规模来提高计算效率。

3.2　快速计算程序

文献[13]介绍了模型降阶法在IGBT风冷散热器上进行瞬态问题快速计算的算例，并给出了提取热传导矩阵及热容矩阵、求解降阶模型的部分源程序。针对算例中只能计算单一热源、不同软件间手动操作、计算效率存在提升空间等问题，本文进行了改进并开发了快速计算程序。具体为：针对不同的IGBT模块芯片依次设置不同热源的参数并提取相应的矩阵文件，降阶后的矩阵使用常微分方程求解器进行求解，然后将结果投影到原模型上。通过时间步与循环功能重复热源设置、求解和投影等步骤，自动实现多热源多载荷步的瞬态计算，其中的时间步由固定步长改为变步长，大幅度减少载荷子步的计算量。借助Matlab软件的调用功能实现程序的自动化，大大简化程序使用的复杂度。使用GRUS稀疏矩阵求解器的算法代替Matlab软件自带的矩阵分解算法，从而进一步提高了模型降阶法的计算效率。

3.3　仿真与试验的对比分析

利用HyperMesh软件，网格基本尺寸取8 mm，将变流模块1的热管散热器和IGBT模块划分为如图8所示的有限元粗糙网格，其中单元数目为212 125个，节点数目为165 560个。

将图7所示的对流换热系数分布插值到有限元模型的流固耦合面上，且以试验工况为输入条件，冷却空气的初始温度取33.8 ℃，利用ANSYS软件可以在1 min内计算出如图9所示的稳态温度场分布云图。由图9可以看出：V4模块的芯片温度最高，为55.105 ℃，1#—7#测温点的温度分别为39.46，45.58，48.82，48.72，45.38，39.30和45.11 ℃，与图6所示的CFD仿真结果相比，温度云图分布规律一致且温度值的差异小于1%，说明将由精细网格CFD模型得到的对流换热系数分布插值到粗糙网格有限元模型是非常有效的，这可以说明此步骤使有限元模型的计算规模急剧减小，计算速度提高了数千倍。根据牛顿冷却定律可知，通过对流换热系数的引入，可以将风冷、水冷等所有的对流换热过程以简单的规律进行描述，因此，利用较少数量的有限元粗糙网格和对流换热系数插值的方法进行瞬态热仿真，可以推广应用到其他类型的对流换热应用中。

图8　有限元网格

图9　基于对流换热系数插值的有限元模型温度场分布云图

取降阶阶数为5，利用快速计算程序对图8所示的有限元模型进行时间范围为0～10 000 s的瞬态计算，可用2 min左右计算出如图10所示的热管散热器的温升曲线。由图10可以看出：在0～1 000 s时间范围内温度迅速上升，之后温度上升的速度变缓，10 000 s时刻的温度处于完全稳定状态；将图10所示的仿真结果与图9所示的稳态温度场结果进行比较，发现最高温度和1#—7#测温点的温度计算结果完全一致，说明快速计算程序所用方法的可行性与准确性。

图10　基于快速计算程序的温升曲线

在时间范围相同的条件下，利用ANSYS软件对图8所示的有限元模型进行瞬态计算，则需要耗时约5 h才能得到如图11所示完全一致的温升曲线，说明笔者开发的快速计算程序的计算准确性和计算效率很高。

4　电气—热联合仿真

4.1　电气仿真

由IGBT模块的用户手册可知，其损耗与电流、管压降及结温有关，因此可将导通电压、开通损耗等曲线拟合成二次函数[14-15]。在选用Matlab/Simulink的SimPowerSystems模型库建立地铁车辆主电路的电气仿真模型时，IGBT模块的损耗等复杂工况的计算常使用线性化的参数，而不是实际的随二次函数变化的参数，这不利于反映电流周期及暂态变化对功耗波动的影响，因此下面对IGBT模块的损耗计算方法进行改进。

IGBT模块的损耗来自IGBT芯片和二极管芯片。其中，对于IGBT芯片主要计算其导通状态下的损耗以及开通和关断过程中的开关损耗，而忽略其关断状态下的损耗；对于二极管芯片主要计算其导通状态下的和反向恢复过程中的损耗，而忽略其关断状态下的和开通过程中的损耗。

以单个桥臂上管IGBT芯片为例，当IGBT模块的电流大于0时，通过拟合的二次函数计算出电流通过该模块时的开关损耗，以下管IGBT芯片脉宽调制下降沿判断单次开关是否结束，开关结束时则累加而得IGBT模块的开关损耗。通过拟合的二次函数计算电流流过模块时的导通功率，再经过积分得到IGBT模块导通损耗。将IGBT模块开关损耗和导通损耗进行累加，再计算出1个周期内的损耗平均值，即为IGBT芯片的平均损耗。以单个桥臂上管二极管芯片为例，当IGBT模块的电流小于0时，通过查找导通功耗表得出电流通过该模块时的导通功率值，然后再积分得到导通损耗。当桥臂输入电流大于0时，通过拟合的二次函数计算开关损耗，以下管二极管脉宽调制上升沿判断上管二极管是否反向恢复，二极管反向恢复时则累加而得到二极管反向恢复损耗。将二极管导通功耗和反向恢复损耗进行累加，再计算出1个周期的损耗平均值，即为二极管芯片的平均损耗。其他IGBT模块损耗的计算与上管相同，不做赘述。

图11为在电机的全速范围(电机转速在20 s内从0加速至300 rad·s-1)内，利用主电路电气仿真得到的电机电磁转矩和IGBT模块损耗随时间变化的曲线。

图11　电机全速范围内的主电路电气仿真结果

4.2　瞬态热仿真

以图11(b)所示的IGBT模块损耗数据为快速计算程序的输入，对图8所示的有限元模型进行瞬态计算，可在7 min内计算出如图12所示的温升曲线。由图12可以看出：V2—V7的芯片在不同时刻产生变化的损耗，可以直接影响其温度的响应；热管散热器安装面上的各测温点因稍远离产生损耗的芯片，相应的温度受到的影响较小；V1和V8距离产生损耗的芯片较远，相应的温度受到的影响最小，这除与距离热源的远近有关外，还与热管散热器底板具有较大的热容量有关。

图12　全速范围的温升曲线

以图11(b)所示的IGBT模块损耗数据为ANSYS软件的输入，对图8所示的有限元模型进行瞬态计算，可得到与图12完全一致的计算结果，但计算消耗的时间长达34.7 h，说明笔者开发的快速计算程序具有很高的计算效率和准确性，与电气联合仿真可为解决复杂工况的瞬态热仿真问题提供有力支持。

5　结　语

针对某地铁车辆牵引逆变器的IGBT模块和热管散热器的温升问题，首先在实验室条件下完成了温升测试，然后以试验结果为输入开展了CFD稳态热仿真、快速瞬态热仿真和电气—热联合仿真，试验数据及不同的仿真方法验证了笔者开发的快速计算方法可同时满足计算效率与计算精度的要求，为解决复杂的瞬态热仿真问题提供了1种新的途径。

[1]LUTZ J, SCHLANGENOTTO H, SCHEUERMANN U, et al. Semiconductor Power Devices, Physics, Characteristics, Reliability[M]. Berlin: Springer-Verlag Berlin Heidelberg, 2011, 380-409.

[2]CIAPPA M. Selected Failure Mechanisms of Modern Power Modules[J]. Microelectronics Reliability, 2002 (42): 653-667.

[3]HAMIDI A, KAUFMANN S, HERR E. Increased Lifetime of Wire Bond Connections for IGBT Power Modules[C]//16th IEEE Applied Power Electronic Conference and Exposition(APEC). Anaheim: IEEE, 2001: 1040-1044.

[4]丁杰, 唐玉兔. 地铁牵引变流器的热管散热器的数值模拟[J]. 大功率变流技术, 2012 (5): 31-33.

(DING Jie, TANG Yutu. Numerical Simulation of Heat-Pipe Radiator for Metro Vehicle Converter[J]. High Power Converter Technology, 2012(5): 31-33. in Chinese)

[5]PERPINA X, JORDA X, VELLVEHI M, et al. Long-Term Reliability of Railway Power Inverters Cooled by Heat Pipe-Based Systems[J]. IEEE Tractions on Industrial Electronics, 2011, 58(7): 2662-2672.

[6]唐强. IGBT元件热管冷却传热性能的实验与数值研究[D]. 兰州: 兰州交通大学, 2013.

(TANG Qiang. Experimental and Numerical Analysis on Heat Transfer Performance of Heat Pipe Cooling for IGBT Component[D]. Lanzhou: Lanzhou Jiaotong University, 2013. in Chinese)

[7]YUN Chansu, MALBERTI P, CIAPPA M, et al. Thermal Component Model for Electrothermal Analysis of IGBT Module Systems[J]. IEEE Transactions on Advanced Packaging, 2001, 24(3): 401-406.

[8]蒋耀林. 模型降阶方法[M]. 北京: 科学出版社, 2010.

[9]WILHELMUS H A Schilders, HENK A van der Vorst, JOOST Rommes. Model Order Reduction: Theory, Research Aspects and Applications[M]. Berlin: Springer, 2008.

[10]TAMARA Bechtold, EVGENII B Rudnyi，JAN G Korvink. Dynamic Electro-Thermal Simulation of Microsystems-a Review[J]. Journal of Micromechanics and Microengineering, 2005, 15(11): 17-31.

[11]RUDNYI E B, KORVINK J G. Model Order Reduction for Large Scale Engineering Models Developed in ANSYS[J]. Lecture Notes in Computer Science, 2006, 3732(1): 349-356.

[12]YANG Y J, YU C C, SHEN K Y. MEMS Heat Transfer Arnoldi-Based Macromodels and the Study of Minimum Required Orders[J]. Tamkang Journal of Science and Engineering, 2005, 8(3): 185-190.

[13]丁杰, 唐玉兔. 模型降阶方法在瞬态热仿真中的应用[J]. 机车电传动, 2014(5): 51-55，65.

(DING Jie, TANG Yutu. Application Research of Model Order Reduction Method in Transient Thermal Analysis[J]. Electric Drive for Locomotives, 2014 (5): 51-55，65. in Chinese)

[14]白保东, 陈德志, 王鑫博. 逆变器IGBT损耗计算及冷却装置设计[J]. 电工技术学报, 2013, 28(8): 97-106.

(BAI Baodong, CHEN Dezhi, WANG Xinbo. Loss Calculation of Inverter IGBT and Design of Cooling Device[J]. Transactions of China Electrotechnical Society, 2013, 28(8): 97-106. in Chinese)

[15]杨宁, 宋术全, 李红. 高速动车组辅助变流器箱体的热仿真设计方法[J]. 中国铁道科学, 2013, 34(3): 87-92.

(YANG Ning, SONG Shuquan, LI Hong. Thermal Simulation Design Method for the Auxiliary Converter Case of High-Speed EMU[J]. China Railway Science, 2013, 34(3): 87-92. in Chinese)