感应加热淬火是机械加工中一项至关重要的工艺。相关研究表明，传动零件失效的主要原因是由于表面损坏造成的，例如齿轮的失效主要源于齿面磨损、齿面点蚀和轮齿折断等。因此，采用合适的工艺对零件进行表面淬火，增加零件表面的硬度和耐磨损程度对提高机械设备的寿命和可靠性具有重要意义。感应加热淬火相比于传统的淬火方式具有操作简便、加热速度快、能源利用率高，更加节能环保等优势，因此广泛用于我国的机床制造、汽车、工程机械和石油机械等领域。随着“碳达峰、碳中和”目标的提出，企业更加重视节能环保，此外国家针对感应加热也提出了相应的电耗定额。为了实现更高的能源利用效率，需要借助于计算机模拟仿真技术对感应加热装置进行设计与优化。通过电磁-热多物理耦合仿真计算得到工件内涡流分布和温度场分布，根据温度场随时间的变化过程可以了解感应加热的效果，以此为依据可以对感应加热频率、激励线圈几何形状、激励线圈与工件的相对位置等参数进行优化设计。有关感应加热的数值模拟方法已有很多研究，张根元等基于ANSYS模拟了感应淬火过程中电流密度、电流频率和加热时间对淬火工艺参数的影响，黄军等利用ANSYS仿真研究了不同感应电流频率的组合对无缝钢管加热速度的影响，张鹏飞等利用Maxwell 3D软件建立了带钢电磁-热耦合模型，研究了加热线圈的电流和频率对带钢温度的影响。以上基于商业软件的仿真分析为感应加热过程的优化提供了指导，然而如何改进多物理耦合计算方法，以便进一步提高计算效率，还需要更多的研究。
本文基于自主开发的电磁及多物理仿真平台 EMPbridge，对齿轮感应淬火中涉及的电磁-热多物理耦合过程进行了仿真，考虑了材料参数随温度的变化，采用频域、时域相结合的方法提高了计算效率，为后续进一步开发感应加热专用仿真和优化软件提供了参考。
一、感应加热过程的电磁-热多物理耦合分析
电磁场数学模型：感应加热中涉及涡流损耗的计算，需要求解磁准静态场，其微分方程形式如下：
式中：E为电场强度(V·m-1) ；H为磁场强度(A·m-1)；D为电位移矢量(C·m-2)；B为磁感应强度(Wb·m-2)；ρe为电荷密度(C・m-3)；J为电流密度(A·m-2)。
由本构关系和电流密度的定义可以得到材料特性与电磁场的关系：
式中：μ为磁导率，σ为电导率，本文中它们是随温度变化的参数，其变化曲线通常通过测量得到。
由于感应加热的激励电流一般为正弦电流，假定材料参数不随磁场和电场变化，上式可以写为复数形式在频域求解：
式中：带上标的
分别表示各场量的复数形式，ω为激励电流的角频率。
结合相应的边界条件，在频域中求解以上磁准静态场问题，即可得到在感应加热频率下激励电流在齿轮工件中产生的涡流分布，进而得到相应的涡流损耗。
热场数学模型：感应加热过程中的热场由瞬态热传导方程描述：
式中：ρ为密度(kg・m-3)；Cp为比热容(J・kg-1・K-1)；T为温度(K)；λ为导热系数(W・m-1・K-1)；Q为热源密度(W・m-3)。热源来自于感应涡流产生的涡流损耗，其表达式为：
结合边界条件和温度的初始条件，在时域中求解以上热扩散方程，可以得到各个时刻对应的温度场分布。
电磁-热多物理耦合过程：感应加热中涉及的电磁-热多物理耦合过程为双向耦合，即交变电磁场在工件中感应出涡流，涡流损耗作为热源引起温度场的变化，温度场变化进一步引起材料参数的变化，而材料参数的变化又导致电磁问题中涡流场分布的变化，其耦合示意图如图 1所示。
图1 感应加热中电磁-热双向耦合示意图
感应加热的激励频率在数十kHz数量级，对应的时间变化周期在数百微秒左右，而加热时间却需要数秒才能完成，两者相差四个数量级，若以适应电磁场变化快慢的时间步进行离散，则总体的时间步数将过于庞大。为解决此多时间尺度问题，在本仿真计算中，电磁场问题在频域中求解，热场问题在时域中进行求解，其耦合求解过程如图2所示，即首先在求解区域给定一个初始温度分布并根据初始温度时的材料参数计算初始时刻。(t=t0)时的频域电磁场，获得涡流分布，然后将涡流损耗作为热源代入，计算时间步长为△t的时域热场问题，得到t1时刻的温度场分布。根据此温度场分布，结合材料参数(μ，σ)随温度变化的曲线，可以得到t1时刻温度变化后的材料参数。由于材料参数的变化，涡流分布也相应变化，因此t1时刻的电磁场需要利用新的材料参数重新计算，并得到新的涡流损耗，以代入热场方程计算下一时刻的温度场分布。之后的过程依此类推，直到到达瞬态热场设置的仿真时长为止。
图2 电磁-热耦合求解过程示意图
值得注意的是，此过程为显式算法，即下一时间步的解可以通过前一时间步的解直接计算得到，因此时间步长的取值不能太大，否则单个时间步内温度变化太大会导致材料参数无法满足线性近似条件。同样，时间步长也不宜选取得过小，否则将增加不必要的计算量和计算时间。可根据某时刻磁导率和电导率的梯度决定时间步长的大小，梯度大(材料参数变化快)时时间步长小，反之，时间步长大。
二、齿轮淬火数值分析算例
齿轮淬火仿真模型的建立：如图3所示，利用EMPbridge平台建立齿轮的几何模型。齿轮的几何尺寸采用模数为2.5的齿轮，齿数为26，可以计算得到齿轮的直径约为65 mm。
图3齿轮几何模型示意图
齿轮的材料是常用的45#钢，其通过测量得到磁导率和电导率参数随温度的变化如图4所示。
图4 45# 钢材料参数随温度的变化
对齿轮进行网格剖分时，需要将涡流透入深度的大小作为网格剖分尺寸的参考。涡流透入深度δ由感应加热频率、磁导率和电导率共同决定，其计算公式为：
仿真中选用的感应加热电流频率为50kHz，在室温20℃时，根据图4查到的相对磁导率和电导率可以计算得到涡流透入深度约为：0.07 mm。当工件升温至900℃时，由于温度已经高于居里温度，工件快速退磁，相对磁导率下降到1左右，且电导率也随温度下降，此时的涡流透入深度增加到2.4mm左右。在透入深度范围内涡流幅值成指数形式下降，因此在齿轮表面透入深度范围内需要对网格进行精细剖分。如图5所示，在齿轮的齿部进行精细剖分，而在远离齿轮边缘的部分采用相对较粗的网格剖分尺寸。
图 5 齿轮网格剖分示意图
涡流与温度场的变化规律：随着感应加热过程的进行，齿轮齿部的温度逐渐上升，而温度的上升将引起材料参数的改变，从而导致涡流场分布产生变化。如图6(a)所示，在初始时刻，温度为室温，涡流透入深度很浅，涡流场集中于齿轮边缘。
当齿轮齿根部位的温度上升到700T时，此时的材料的磁导率到达居里温度附近，根据图4(b)的曲线可知，到达居里点后温度的小幅升高将导致磁导率的大幅下降，此时透入深度增加，涡流场分布向齿轮内部扩散，如图6(b)所示。
随着感应加热的进一步进行，齿轮的整个齿部都已经升温至居里温度以上，此时涡流分布再次沿齿轮外沿分布，但是透入深度增加到2.4mm左右，如图6(c)所示。
图6 涡流分布随温度的变化趋势
由以上结果可见，涡流分布随温度的上升而变化，并反过来影响温度的分布，因此感应加热过程需要考虑电磁-热的双向耦合过程。
齿轮材料参数随温度的变化：由于感应加热过程中，齿轮齿部的温度高，而齿轮内部的温度低，因此齿轮不同部位的材料参数并不一致，而是与温度场的分布相关。如图7所示，在加热过程中，温度高的部位，磁导率和电导率下降明显，而齿轮内部由于温度升高不明显，其材料参数与室温时基本相同。对比图6(b)和图7可以进一步印证1.3节中所阐述的电磁-热双向耦合过程。图7中齿根部位的温度升高导致此区域的磁导率和电导率大幅下降，根据公式(10)可知透入深度将升高，观察图6(b)可以发现齿根部位的电流密度升高，涡流透入齿轮齿根部位。
图 7 齿轮材料参数随温度的分布
三、结论
本文针对齿轮淬火加热的实际工程问题，对感应加热中涉及的多时间尺度电磁-热多物理耦合过程进行了剖析。在仿真中考虑了电导率和磁导率参数随温度的非线性变化，阐述了利用频域-时域混合的方法对电磁-热双向耦合问题进行仿真计算的方法。通过分析仿真结果，给出了涡流场分布受热场分布影响的规律，为齿轮感应加热的优化设计提供了参考。
参考文献略.