基于克里金插值的连铸坯感应加热三维动态有限元建模与仿真

杭州电子科技大学自动化学院 杨钦塔 方 东 程晓明

基于克里金插值的连铸坯感应加热三维动态有限元建模与仿真

杭州电子科技大学自动化学院 杨钦塔 方 东 程晓明

本文采用有限元法建立了连铸坯感应加热过程的三维动态模型，模拟方形连铸坯以固定速度通过单个中频（1000Hz）感应器加热的过程，预测了连铸坯随时间和空间变化的温度分布。该模型利用kriging插值建立非均匀三维初始温度场模拟钢坯进入感应器前温度分布呈现的“头低尾高”现象。仿真结果表明，在轴向磁通的加热方式下，连铸坯存在很严重的“端部效应”，即加热后连铸坯端部局部温度比周围温度高很多，调整电流能有效改善轴向温度分布不均的问题，但是不能消除端部效应。

克里金；三维动态模型；端部效应

0 引言

感应加热具有很多优点，如精确的加热深度和加热表面积、损耗低、工作环境清洁、加热时间短等，因而能适合生产过程的自动化[1]。为了获得好的加热效果（如一致的温度分布、控制晶体尺寸等），设计者多借助经验和实验数据获得比较理想的加热策略和过程参数，成本高、精度低、耗时长。随着计算机技术的发展，各种电热设备的制造商已经开始利用数值模拟来辅助设计和开发过程[2]。

感应加热过程是电、磁、热三大物理场相互耦合作用的过程，其数学模型的特点是包含一组耦合的偏微分方程，任意瞬时时刻，电磁场分布显式依赖于被加热物体由热传导方程决定的历史温度，因此该方程组不能解耦[2]，直接求解难度很大。因此，很多文章都将感应加热过程的耦合场问题作为研究课题。文献[3]提出了无限长工件感应加热过程的二维数学模型和数值模型。C.Chaboudez, S.Clain等[4][5]在文献[3]的基础上分别研究了长工件和截面不一致的轴对称工件的感应加热过程二维数值模型。文献[5]基于FEM建立了感应加热问题的2D和3D静态模型。K.F Wang[6]考虑线圈与工件的相对运动，通过线圈网格重新划分的方式实现感应加热过程的二维动态模型。Y.Favennec[7]不但实现了二维动态数值模型，而且结合模型对加热过程提出优化策略。这些模型通过对比实验都得到了比较满意的结果，但是也是有缺陷的。首先，它们都以均匀初始温度场（冷坯25℃）作为初始状态，而在实际当中，如本文主要研究的连铸坯感应加热过程，钢坯初始温度场大都是不均匀的，轴向存在着“头低尾高”现象。其次，由于钢坯长度有限，感应加热过程存在明显的“端部效应”，即钢坯两端局部温度比周围温度高出很多，如果模型中忽略这一点，其潜在危险是调整的过程参数导致钢坯局部温度过高而产生淬火裂纹或者柔韧度不够等问题[8]。因此，为了获得理想的控制方法，需要建立3D有限元动态模型，调整过程参数进行仿真分析。

本文利用有限元软件ANSYS建立了轴向磁通动态感应加热模型。模型考虑了工件的实际三维几何结构，工件与感应器的相对运动，材料属性与温度的非线性依赖关系，热辐射边界条件等问题。该模型还引入了kringing插值算法构建工件初始温度场，使模型更加符合实际。

1 数学模型及边界条件

1.1 电磁场数学模型

感应加热过程电磁场可以由麦克斯韦方程微分形式表示：

其中：

为了避免求解复杂的一阶偏微分方程组，往往引入位函数作为间接变量。常用的方法是引入矢量磁位A和标量电位φ，它们的定义为：

得到涡流区的控制方程为：

为了确保矢量磁位A的唯一性，应该对A施加库仑规范即：

1.2 温度场数学模型

热传递问题遵循能量守恒定律：热能或者机械能通过边界进入系统的速率减去能量通过边界离开系统的速率，再加上系统体内产生能量的速率，必须等于系统存储能量的速率，其方程[Book]为：

直角坐标系下微元体应用能量守恒定律，得到热传导方程如下所示：

边界条件：高温下通过辐射损失的热量在低速感应加热形式下是相当可观的。大多数金属在温度为850℃以上时，通常产生的热损失很小。当温度在850℃以上时，热损失变得显著[]。热辐射由下式给出：

其中，ε为表面辐射系数，σ为斯蒂夫—玻尔兹曼常数，T1为表面绝对温度K(℃+273)，T2为环境绝对温度K(℃+273)，Qe为每平方米表面积上的辐射热损失。感应加热的工况空气流动缓慢，热对流损失相对于输入功率很小，在模型中忽略这部分损失。

1.3 Kriging模型

回归多项式是l个可选函数的线性组合。

其中R(θ,xi,xj)是参数为θ的相关模型。

在预测点x处定义：

以及随机过程的相关矩阵R：

x处的响应值y(x)预测估计值的表达式为[]：

其中：

文献[9]提供了一种Matlab的Kriging工具箱。它提供了三种回归模型，分别为0、1和2阶多项式，六种相关模型，常用的有指数、球形和高斯模型。本文中回归模型选择1阶模型regploy1，相关模型选择指数模型，相关参数初值为θ= [1,1,1]。

图1 几何模型

2 有限元模型

问题描述：长度为3m的方形高温钢坯通过1m长的单感应器进行加热，钢坯的传送速度为5m/min。线圈为方孔形，内孔边长为0.46m，外孔边长为0.51m，匝数为30。三维几何模型的1/4模型如图1所示。

实现钢坯的运动：

要实现钢坯的运动，一般有两种方法，第一种是重新划分线圈网格来实现钢坯运动[]，但是这种方法同样十分耗时，更重要的是重新划分的网格可能产生很大的畸变导致结果不准确甚至运算终止。第二种方法是预先将线圈运动路径上的的区域分别定义且划分好，将线圈和空气的电磁材料属性来回更替，相当于线圈在相对于空气运动。本文基于方法二，利用载荷的移动来实现线圈与钢坯的相对运动如图2所示。先将线圈的运动路径全部预先划分网格，在t时刻，第一对箭头包含的线圈区域加载电流载荷，在t +△t时刻第二对箭头包含的线圈区域加载电流载荷。载荷运动的距离恰好为钢坯传送的距离△s。划分的段数越多，△s越小计算的精度就越高，相对的计算时间也越长。分段数可以利用试错法得到，将仿真结果与预分析结果或者实际数据逐渐逼近，在精度与计算时间之间平衡得到一个合适值，本文通过这种方法设定的分段数为25。这里需要对模型作如下假设：（1）忽略线圈每匝之间的间隙；（2）忽略线圈运动路径上线圈材料与空气电阻率和磁导率不同产生的影响。

图2 载荷移动示意图

Ansys提供了多种耦合场分析方法，本文采用的是顺序耦合物理分析，即顺序的求解物理场，并将前一次求解结果作为下一次分析的载荷，FEM求解流程图如图3所示。前面提到初始温度场是不均匀的，指定各个节点的初始温度后，需要一个稳态热分析来建立初始温度场，设置很短的时间长度，如0.001s。各个节点的温度设定，根据实际情况需要频繁的使用和修改，可以利用ANSYS创建宏文件。Kriging的预报函数可以计算钢坯每个节点的温度插值结果，但是在ANSYS中按编号逐个给几万个节点设定初始温度显然不可行。本文采用的方法是，先由MATLAB工具包dace建立Kriging模型，模型中包括式(18)中的矩阵和规范化后的测点矩阵[s]，在宏文件中用数组保存这些矩阵。接着在宏中获得

一个节点ω的坐标，再根据选择的回归模型和相关模型的形式求得f(ω)T和r(ω)T，最后利用公式(18)求得温度值。本文建立的有限元模型如图4所示。

图3 FEM求解流程图

图4 (a) 有限元模型

图4 (b) 网格划分

3 实例与分析

实验模拟的工况为某钢厂连铸坯感应加热过程。选取的钢坯为是45号钢，几何参数如图5所示，其参数与温度的非线性关系如图5所示，加热一次时间为50.4s。

图5 45号钢在不同温度下的物理参数：a电阻率；b相对磁导率；c比热容；d热导率

初始温度场：

实际钢坯的温度分布呈现“头低尾高”的现象，头部温度比尾部温度整体大致高100℃。

初始温度插值结果，如图6所示：

图6 初始温度插值结果

实验一：梯形电流

图7(a)是实际该钢厂钢坯通过两个感应器的电源输出电流曲线（保持端电压不变），可见电源电流整体呈现为梯形，且左右对称。原因是在钢坯进入第一个感应线圈的过程中，系统阻抗从空载到满载，当感应器中完全包含钢坯的时候，虽然钢坯还在前进，但是系统的阻抗保持不变，电流也维持恒定。钢坯进入第二个感应器的时候，由于钢坯的长度大于两个感应器与其间隙的长度和，所以会出现电流曲线两次上升的情况，钢坯出感应器的情况同理。因此，对于单感应器加热的情况，本文设计了如图7(b)所示的的电流载荷曲线。需要指出，图7(b)的电流曲线指的是每匝线圈流过的电流，它与电源电流存在如下关系：

图7 (a) 实际电源输出电流曲线

图7 (b) 实验一感应线圈电流曲

钢坯的温度分布即是空间坐标的函数又是时间的函数，为了更好的观察实验结果，图8中左边定义了1/4模型上的四条路径R1、R2、R3、RZ，R1为钢坯中心到钢坯表面沿y方向的路径，R2为钢坯表面中心到钢坯角部沿x方向的路径，R3为钢坯中心到角部沿对角线方向的路径，RZ表示从钢坯表面中部从头部到尾部的路径。右边定义了1/8截面上的六个坐标点A、B、C、D、E、F。

图8 后处理路径与截面节点

图9将R1、R2、R3路径(z=1.5m截面)和RZ路径在加热前后的温度曲线作了对比，R1主要反映了感应加热主要集中在表面，导致表面温度大幅上升，另外，随着加热时间增加，表面中心和边角的温度差也逐渐增加。从R2可以看出，初始时刻，中心温度比表面温度高，由于感应加热集中在表面，表面温度迅速上升，最终超过中心温度。同时，经过感应加热表芯温差大大缩小，如果均温一段时间，表芯温差能够更小。从R3看出，感应加热让钢坯横截面温度分布更加均匀，钢坯最热的部分在内表面的部分，通过均温可以让温度分布更加均匀。从RZ路径看出，初始温度是呈现前低后高趋势的，在梯形电流作用下，加热头尾温度也同样呈现前低后高的总体趋势，同时还有明显的端部效应。图9为横截面结点温度历程图，图10为不同时刻钢坯感应加热图。

图9 加热前后各路径上的温度曲线图

图10 横截面结点温度历程图

图11 不同时刻钢坯感应加热图

实验二：优化电流

通过实验一，发现实际中感应加热结束后，温度的不均匀性是普遍存在的。实际中一般设置均温段对钢坯进行均热，但是这种方法对于缩小截面温差比较有效，而对于缩小轴向温差效果不明显，原因是钢坯长度很长，轴向高温区到低温区的距离远，由傅里叶定律：

轴向热流密度低，所以均热效率也低。因此在加热过程就应该加入控制手段。本文希望通过调节线圈电流来实现缩小温差的目标，设计了如图12所示的电流曲线。

实验结果如图13所示，与图9比较可以发现R1、R2、R3的变化规律几乎是一致的，只有RZ路径上的温度分布有了明显的变化，钢坯中间段的温度分布更加均匀，这说明钢坯表面温度与线圈电流之间存在某种线性关系。同时，两个实验表明，这种加热方式，轴向感应加热的端部效应不可避免，原因可能是涡流在端部聚集，导致电流密度很高，从而产生很高的热量。所以，实际中应当适当控制电流，不能让钢坯的端部温度过高。

图12 实验二感应线圈电流曲线

图13 加热前后各路径上的温度曲线图

4 结论

本文利用有限元软件ANSYS建立了轴向磁通动态感应加热模型。通过引入了kringing插重建了连铸坯进入感应加热器前的初始温度场。在设定的初始温度条件下，预测了在梯形电流产生的电磁场作用下连铸坯感应加热过程的温度分布情况。结果表明，该模型能够立体的反映了钢坯感应加热的过程的温度变化情况。

通过调整电流曲线，能够消除连铸坯初始的“头低尾高”现象，发现表面温度与电流之间存在某种线性关系，这有利于研究温度的实时控制或者前馈控制。

该模型不但能仿真连铸坯感应加热过程，对于冷坯、热轧中间坯等的加热过程，只要对模型稍做修改，同样可以模拟。

[1]Davies J,Simpson P.Induction heating handbook[M]. London:McGraw-Hill,1979.

[2]Lavers J D.Numerical solution methods for electroheat problems[J].Magnetics,IEEE Transactions on,1983,19(6):2566-2572.

[3]Clain S,Rappaz J,Swierkosz M,et al.Numerical modeling of induction heating for two-dimensional geometries[J].Mathematical models and methods in applied sciences,1993,3(06):805-822.

[4]Chaboudez C,Clain S, Glardon R,et al.Numerical modelling of induction heating of long workpieces[J].IEEE transactions on magnetics,1994,30(6):5028-5037.

[5]Chaboudez C,Clain S,Glardon R,et al.Numerical modeling in induction heating for axisymmetric geometries[J].Magnetics,IEEE Transactions on,1997, 33(1):739-745.]

[6]Wang K F,Chandrasekar S,Yang H T Y.Finite-element simulation of moving induction heat treatment[J].Journal of materials engineering and performance,1995,4(4):460-473.

[7]Favennec Y,Labbé V,Bay F.Induction heating processes optimization a general optimal control approach[J].Journal of computational physics,2003,187(1):68-94.

[8]Sadeghipour K,Dopkin J A,Li K.A computer aided finite element/experimental analysis of induction heating process of steel[J]. Computers in industry,1996,28(3):195-205.

[9]Lophaven S N,Nielsen H B,S ndergaard J.DACE-A Matlab Kriging toolbox,version 2.0[R].2002.

杨钦塔（1991—），浙江温州人，硕士，现就读于杭州电子科技大学。