钢板移动式感应加热的多场耦合数值分析

张雪彪，陈诚，刘玉君，张立卫

(大连理工大学船舶工程学院工业装备结构分析国家重点实验室，辽宁大连116024;2.大连理工大学数学科学学院，辽宁大连116024)

在船厂，水火弯板工艺广泛应用于双曲率板的成形过程中。众所周知，火焰热源的水火弯板工艺强烈依赖熟练工人的热源控制技术，以达到控制钢板变形形状的目的。然而，火焰热源是乙炔气体和氧气发生化学反应的热源，精确控制火焰作用在钢板表面的热通量是困难的;同时，整个燃烧过程释放大量的CO2气体和噪音，对环境造成污染。为了克服火焰热源的这些缺陷，人们开始研究采用感应热源的水火弯板工艺:感应加热直接利用涡流对钢板进行加热，热效率高;通过调节频率、功率等电学参数可以实现对热输入量的精确控制;最后，感应加热工作环境好、噪音小，有利于环保。

尽管感应加热具有许多优点，由于缺乏足够的加热机理的理论分析以及变形特点的深入探索，此项工艺技术在船体板成形领域仍未被广泛采用。目前，国内外已经有学者做了初步研究。例如，LEE J H［1］研究了感应加热在板材成型中的可行性，指出火焰加热和感应加热具有很大的相似性。帅克刚、罗宇［2-3］依据数值与实验结果的比较建立了感应加热的感应热模型，在数值计算时直接加载该热源模型进行热-结构的计算。这一方法缺少电磁-热耦合的分析过程，忽略了感应加热过程中的集肤效应、端部效应等感应加热的基本特点。LEE K S［4］利用商业软件 JMAG-Studio开发了钢板移动加热的电磁-热-结构耦合模型，分析了加热过程中的温度分布和弯曲行为。杨玉龙和张雪彪［5-6］也利用COMSOL软件开发了钢板移动式感应加热的磁-热耦合的数值模型，对钢板加热过程中的电磁场和温度的耦合特点进行了研究。虽然人们对钢板感应加热已经做了一些工作，但是对于钢板在感应加热过程中的温度机理和变形机理仍需仔细研究。因为钢板感应加热依靠金属内部的感应涡流进行加热，然后进行热传导，与主要依靠热传导的火焰热源有本质的不同。基于此，本文在前面工作的基础上建立了移动式感应加热的电磁-热-结构耦合的数值模型，研究多场耦合条件下钢板的温度和变形机理，并利用钢板感应加热实验验证了数值模型的有效性。

1 数学模型

如图1所示，钢板移动式感应加热系统主要由3部分组成:钢板、感应器、Π型导磁体。感应加热过程中，感应器在钢板上方以一定速度由板内向板边进行移动加热。加热过程中保持钢板与感应器之间的空气间隙不变。

图1 钢板感应加热示意图Fig.1 Schematic of steel plate induction heating

1.1 感应加热弯板成型的理论

描述电磁感应现象的控制方程是麦克斯韦方程组。因为感应加热的电磁场是时谐电磁场［7-8］，所以其麦克斯韦方程组为

式中:H是磁场强度，J是电流密度，E是电场强度，B是磁感应强度，ω是正弦变化角频率(ω=2πf，f是感应加热的电流频率)，D是电位移，ρ是电荷密度。

引入磁场矢势A，满足:

此外，B、H、J、E两两之间满足下列关系式:

式中:μ是磁导率，σ是电导率。电场E可以表示为

其中，φ是标量电势。将式(5)～(8)代入式(1)整理得电磁场的控制方程:

其中，Js是源电流密度。本文中采用的数值模型是磁矢量势法的三维有限元模型。基于三维磁矢量势法的边界条件如图2所示。通量垂直边界条件:

通量平行边界条件:

远场边界条件:

电势边界条件:感应器尾部截面电势φ为零。

图2 电磁边界条件Fig.2 Electromagnetic boundary conditions

感应加热过程中，涡流集中分布于工件表层，这种现象称之为集肤效应，集肤层深度δ的计算公式:

受到涡流的作用后，依据焦耳-楞次定律，钢板加热区域内感应涡流产生的热量q。

式中:Je是感应涡流密度。钢板内部的热传导现象由微分控制方程(15)来描述:

式中:ρ是钢板密度，Cp是钢板的比热，T是钢板的温度，k是钢板的导热系数。

钢板感应加热中需要考虑空气冷却、水冷却及辐射换热等换热过程，所以钢板的热边界条件是:

与空气接触的钢板:

与水接触的钢板:

式中:qa、qw分别为空气冷却和水冷却的热损量，ε是放射率，σb是玻尔兹曼常数，T∞是环境温度，T是钢板壁面的温度，A是在温度T范围内的面积，h1是空气的对流换热系数，利用平板外部自然对流模型求解，h2是水的对流换热系数，利用流动沸腾理论求解，h1、h2的计算过程参考文献［10］。加热时钢板各表面施加空气冷却边界条件(16)，加热结束后进行正面水冷，在水冷区域施加水冷边界条件(17)。

钢板的塑性应变过程通过下述方程来描述［11］:

式中:σij，i=∂σij/∂xi，ui，j=∂ui/∂xj，uj，i=∂uj/∂xi，σij是应力张量，εij是应变张量，Fi是外力矢量(本分析值中为零)，E是弹性模量，υ是泊松比，δij是克罗内克函数，e是应变张量的对角线部分。

应变张量εij由3部分组成叠加而成:

1.2 数值模拟方法

钢板感应加热过程中，电磁场与温度场的耦合是一个交互式的耦合过程。从式(16)可以看出，在焦耳-楞次定理的作用下，电磁分析得到的涡流值转化为热载荷q进行热分析。随着钢板温度的升高，钢板的磁导率μ和电导率σ会发生改变，则在下一载荷步的电磁计算时必须更新这些电磁学参数，同理在进行下一载荷步的热分析时，材料的热物性参数也要更新;此外，感应器移动时，钢板、感应器、空气组成的求解域时刻在发生变化。上述两点是建立电磁-热-结构耦合数值模型所需解决的问题。

本文采用ANSYS软件的物理环境法进行钢板电磁感应加热的多场耦合数值计算。图3是电磁-热-结构耦合计算的流程图，首先设置电磁场计算的初始条件以及初始模型的约束方程，进行电磁谐波分析，计算得到的涡流值和初始钢板温度一起作为热输入量进行热分析。热分析后，判断是否进行下一载荷步的电磁求解?如果需要进行下一载荷步的电磁场求解，则移动感应器模型，重新建立感应器和钢板之间的约束方程，更新钢板的材料属性，进行新的电磁分析，从而实现电磁-热物理场之间的交互式耦合;如果仍然需要进行温度场的求解，则根据钢板表面的温度，更新热物性参数，再次进行温度场计算，直到热分析结束。当电磁-热耦合分析结束后，开始进行热-结构的计算，将电磁-热耦合的温度计算结果作为热载荷作用于结构模型，从而完成对钢板变形场的求解。电磁和热分析的边界条件如前所述，结构变形的边界条件为四端简支。

图3 电磁－热－结构耦合计算流程图Fig.3 Flowchart of coupled electromagnetic-thermal-structural calculation

图4 钢板感应加热的有限元模型Fig.4 Finite element model of steel plate during induction heating

在电磁场计算中钢板、感应器与空气模型使用六面体SOLID 97单元，温度场分析中钢板模型使用SOLID 70单元，结构变形求解中钢板模型使用SOLID 185单元。图4(a)是钢板感应加热的整体模型，感应器处于空气包中，图4(b)是感应器和钢板结合部分的局部网格图。

为了建立感应器和钢板之间的移动关系，感应器和钢板之间的空气间隙被分成两层，空气间隙的上层是移动的，空气间隙的下层是静止的，在两层网格的结合处设置约束方程以进行感应器移动时的耦合计算。由于磁场强度在感应器周围较高，所以在感应器附近区域的网格划分较密集。考虑到集肤效应，钢板表层也划分了密集的网格。

钢板的电磁学参数和热物性参数均与温度相关［6］，力学性能数据来自参考文献［12］。

2 钢板感应加热实验

为了验证数值模型的正确性，课题组在大连理工大学造船工艺实验室采用钢板感应加热设备SPZ-70进行了实验。

图5是温度和变形测量点的布置图，P1点是温度测量点，P1(620，150，18)点位于钢板表面，加热线的纵向中间位置，距离加热线20 mm，该点的温度利用Raytek公司的自动红外测温仪进行测量。横向收缩量的测量则通过两点距离的变化得出，即在加热线两侧确定11对测量点，通过数显游标卡尺(精0.02 mm)测量它们在加热前后的距离值，从而得出横向收缩量值。加热时的工艺参数如表1所示，感应加热现场如图6所示，钢板焰道在加热后进行正面水冷。

表1 感应加热实验的工艺参数Table 1 Technical parameters of induction heating

图5 实验布置Fig.5 Experiment setup

图6 感应加热实验图Fig.6 Induction heating experiments

3 结果与讨论

3.1 数值实验数据对比

图7是P1点温度的实验值和计算值的对比图，在加热阶段两者的平均相对误差是6.86%。在34～70 s的时间段内感应器经过P1点区域，钢板金属受到感应涡流作用，温度急速上升至550℃左右。当感应器移过P1点后，P1点温度开始下降，下降速率较慢。111 s后进入冷却阶段，P1点的温度在水冷作用下急速下降。P1点温度计算值和实验值的误差在可接受的范围内，钢板温度的计算结果是正确的。

图8是横向收缩量计算值和实验值的对比图。两者的变化趋势相同，平均相对误差为10.77%。图11中在板边处的收缩变形偏低，这主要是由感应器加热到板边后出现的端部效应导致的。

图7 P1点温度计算值实验值对比图Fig.7 Comparison between FEM results and measurement data of temperature at P1

图8 计算与实验横向收缩量对比图Fig.8 Comparison between FEM results and measurement data of transverse shrinkage

3.2 钢板移动式感应加热的电磁-热-结构分析

钢板在移动式感应加热过程中，钢板的电磁场场分布、温度分布是不断变化的，这与静止式感应加热有较大的区别。根据感应器的移动位置，本文从初始加热位置开始确定了6个典型截面，对应截面的位置和时刻如图9所示。

图10是1、41、111 s 3个时刻钢板上表面的磁通量密度云图。如图10(a)，t=1 s时钢板表面上的磁通密度前后对称分布;随着加热的进行，感应器不断前移，如图10(b)，磁通密度分布较初始状态发生较大变化。由图11可以看出，t=41 s时刻，感应器前端钢板的涡流值(2.8×108A/m2)和磁通量密度(1.5 T)远大于感应器中后方钢板的涡流值(0.25×108A/m2)和磁通量密度(0.04 T)，由式(6)可知，B与μ是正比关系，感应器前端钢板温度较感应器中后方钢板温度低，相对导磁率大，所以此区域的磁通密度大，感应涡流大。综上所述:钢板移动式感应加热过程中，感应涡流分布集中于感应器前端的钢板区域内，有利于感应器前侧钢板金属的快速升温。

图9 电磁－热－结构分析中的钢板截面示意图Fig.9 Sectional schematic diagram of steel plate during the electromagnetic-thermal-structural analysis

图10 钢板表面磁通密度云图Fig.10 Contours of the magnetic flux density on the surface of the steel plate

图10(c)是加热到板边处时的磁通密度云图，由图中可以看出在板边处磁通密度扩张，呈现出较强的端部效应。图12是4个时刻下钢板上表面及厚度方向(如图9所示)的温度等值线图，其中厚度方向的等值线图做了放大处理，这些图反映了钢板温度的动态变化过程。图12(a)是1 s时刻的温度等值线图，钢板表面温度已经达到331℃，厚度方向的加热深度较浅;钢板上表面的温度上升很快，从图12(b)的11 s时刻，钢板最高温度已升至900℃以上;自11 s以后，如图12(c)的41 s时刻，钢板最高温度趋于稳定，钢板进入了准稳态加热阶段;在图12(d)的111 s时刻，钢板上表面温度等值线有明显的外扩趋势，与图10(c)中的感应磁场相对应，形成了感应加热特有的端部效应现象。

图11 41 s时刻加热线的涡电流与磁通密度分布图Fig.11 Distribution of eddy current and magnetic flux density in the heating line at time 41 s

图12 钢板表面温度和横截面温度等值线图Fig.12 Contour map of the surface temperature and sectional temperature of steel plate

感应加热过程中，钢板能快速升温并进入准稳态加热状态的原因是:感应加热起始的磁通密度较大，涡流值很大，产生大量焦耳热，钢板温度快速上升，当钢板表面温度超过居里点温度(720℃)后，集肤层内的金属失去磁性(相对导磁率趋于1)，磁通密度减少，涡流值减少，同时感应涡流向钢板内部渗透，钢板表面的热输入量减小，钢板的热输入量与钢板内部的热传导外加空气换热损失的热量总和达到平衡，所以钢板表面温度上升速率减缓，趋于平稳。

图13 是Y分别为0.18、0.06、0 m(截面位置如图9所示)3个横截面上表面节点分别在41、81、111 s时刻的X方向应力分布图。

图13 Sx应力分布图Fig.13 Distribution of stress of Sx

各截面节点的Sx值均是负值，即加热线两侧的节点均承受横向压应力，其中Y分别为0.18、0.06 m两截面的Sx分布趋势相似;然而，Y=0 m横截面上的应力在远离加热线处呈现出较大的压应力，与前面的截面温度分布图相对应，Y分别为0.18、0.06 m的两截面均已处于准稳态加热，如图12(c)所示，高温区集中在加热线附近;而Y=0 m截面，由于感应加热端部效应的存在，如图12(d)所示，钢板端部高温区外扩，应力分布也发生改变，应力最大值出现在远离加热线处。

4 结束语

本文基于Ansys有限元软件采用耦合约束方程的方法开发了钢板移动式感应加热的多场耦合数值模型，实验数据与计算结果吻合，证明了数值模型的正确性。利用电磁-热-结构耦合的数值模型分析了钢板在线加热过程中的涡流分布、钢板的温度分布和结构变形;钢板感应加热速度快，迅速进入准稳态加热状态，加热时感应加热的涡流集中于感应器前端，有利于冷却金属的迅速升温;在感应器加热至板边时会产生端部效应，影响板边的收缩变形。

［1］LEE J H，LEE K H，YUN J S.An electromagnetic and thermo-mechanical analysis of high frequency induction heating for steel plate bending［J］.Key Engineering Materials，2006，326:1283-1286.

［2］帅克刚，罗宇，沙万华，等.高频感应加热弯板成型温度场的有限元分析［J］.锅炉技术，2004，35(3):52-54.SHUAI Kegang，LUO Yu，SHA Wanhua，et al.FME analysis of the temperature field of high frequency induction heating in plate bend molding［J］.Boiler Technology，2004，35(3):52-54.

［3］罗宇，江晓玲，邓德安.材料力学性能对高频感应弯板成形的影响［J］.材料工程，2005(8):35-38.LUO Yu，JIANG Xiaoling，DENG Dean.Influence of mechanical properties of materials on high frequency inducting plate bending［J］.Journal of Materials Engineering，2005(8):35-38.

［4］LEE K S，KIM S W，EOM D H.Temperature distribution and bending behaviour of thick metal plate by high frequency induction heating［J］.Material Research Innovation，2011，15(1):283-287.

［5］杨玉龙，张雪彪，刘玉君.船体板高频感应加热的多场耦合数值分析［J］.船海工程，2011，40(6):24-27.YANG Yulong，ZHANG Xuebiao，LIU Yujun.Numerical analysis of multi-physics coupling for high frequency induction heating of hull plate［J］.Ship ＆ Ocean Engineering，2011，40(6):24-27.

［6］张雪彪，杨玉龙，刘玉君，等.钢板弯曲工艺中的高频感应加热过程数值分析［J］.中国造船，2011，52(2):108-116.ZHANG Xuebiao，YANG Yulong，LIU Yujun，et al.Numerical analysis of the high-frequency induction heating process applied to steel plate bending technology［J］.Shipbuilding of China，2011，52(2):108-116.

［7］RENHART W，STOGNER H，PREIS K.Calculation of 3D eddy current problems by finite element method using either an electric or a magnetic vector potential［J］.IEEE Transactions on Magnetics，1988，24(1):122-125.

［8］谢德馨.三维涡流场的有限元分析［M］.2版.北京:机械工业出版社，2007:8-146.

［9］CHO K.Coupled electro-magneto-thermal model for induction heating process of a moving billet［J］.International Journal of Thermal Sciences，2012，60:195-204.

［10］汪骥.水火弯板自动化加工工艺的关键技术研究［D］.大连:大连理工大学，2007:32-42.WANG Ji.Research on pivotal techniques for automatic line-heating process［D］.Dalian:Dalian University of Technology，2007:32-42.

［11］PANTELYAT M G，FELIACHI M.Magneto-thermo-elasticplastic simulation of inductive heating of metals［J］.EPJ Applied Physics，2002，17(1):29-34.

［12］BIRK-SØRENSEN M，JENSEN J J.Simulation of welding distortion in ship section［D］.Lyngby:Technical University of Denmark，1999:95-96.