杭州電子科技大學(xué)自動化學(xué)院 楊欽塔 方 東 程曉明

基于克里金插值的連鑄坯感應(yīng)加熱三維動態(tài)有限元建模與仿真

杭州電子科技大學(xué)自動化學(xué)院 楊欽塔 方 東 程曉明

本文采用有限元法建立了連鑄坯感應(yīng)加熱過程的三維動態(tài)模型，模擬方形連鑄坯以固定速度通過單個中頻（1000Hz）感應(yīng)器加熱的過程，預(yù)測了連鑄坯隨時間和空間變化的溫度分布。該模型利用kriging插值建立非均勻三維初始溫度場模擬鋼坯進(jìn)入感應(yīng)器前溫度分布呈現(xiàn)的“頭低尾高”現(xiàn)象。仿真結(jié)果表明，在軸向磁通的加熱方式下，連鑄坯存在很嚴(yán)重的“端部效應(yīng)”，即加熱后連鑄坯端部局部溫度比周圍溫度高很多，調(diào)整電流能有效改善軸向溫度分布不均的問題，但是不能消除端部效應(yīng)。

克里金；三維動態(tài)模型；端部效應(yīng)

0 引言

感應(yīng)加熱具有很多優(yōu)點，如精確的加熱深度和加熱表面積、損耗低、工作環(huán)境清潔、加熱時間短等，因而能適合生產(chǎn)過程的自動化[1]。為了獲得好的加熱效果（如一致的溫度分布、控制晶體尺寸等），設(shè)計者多借助經(jīng)驗和實驗數(shù)據(jù)獲得比較理想的加熱策略和過程參數(shù)，成本高、精度低、耗時長。隨著計算機(jī)技術(shù)的發(fā)展，各種電熱設(shè)備的制造商已經(jīng)開始利用數(shù)值模擬來輔助設(shè)計和開發(fā)過程[2]。

感應(yīng)加熱過程是電、磁、熱三大物理場相互耦合作用的過程，其數(shù)學(xué)模型的特點是包含一組耦合的偏微分方程，任意瞬時時刻，電磁場分布顯式依賴于被加熱物體由熱傳導(dǎo)方程決定的歷史溫度，因此該方程組不能解耦[2]，直接求解難度很大。因此，很多文章都將感應(yīng)加熱過程的耦合場問題作為研究課題。文獻(xiàn)[3]提出了無限長工件感應(yīng)加熱過程的二維數(shù)學(xué)模型和數(shù)值模型。C.Chaboudez, S.Clain等[4][5]在文獻(xiàn)[3]的基礎(chǔ)上分別研究了長工件和截面不一致的軸對稱工件的感應(yīng)加熱過程二維數(shù)值模型。文獻(xiàn)[5]基于FEM建立了感應(yīng)加熱問題的2D和3D靜態(tài)模型。K.F Wang[6]考慮線圈與工件的相對運動，通過線圈網(wǎng)格重新劃分的方式實現(xiàn)感應(yīng)加熱過程的二維動態(tài)模型。Y.Favennec[7]不但實現(xiàn)了二維動態(tài)數(shù)值模型，而且結(jié)合模型對加熱過程提出優(yōu)化策略。這些模型通過對比實驗都得到了比較滿意的結(jié)果，但是也是有缺陷的。首先，它們都以均勻初始溫度場（冷坯25℃）作為初始狀態(tài)，而在實際當(dāng)中，如本文主要研究的連鑄坯感應(yīng)加熱過程，鋼坯初始溫度場大都是不均勻的，軸向存在著“頭低尾高”現(xiàn)象。其次，由于鋼坯長度有限，感應(yīng)加熱過程存在明顯的“端部效應(yīng)”，即鋼坯兩端局部溫度比周圍溫度高出很多，如果模型中忽略這一點，其潛在危險是調(diào)整的過程參數(shù)導(dǎo)致鋼坯局部溫度過高而產(chǎn)生淬火裂紋或者柔韌度不夠等問題[8]。因此，為了獲得理想的控制方法，需要建立3D有限元動態(tài)模型，調(diào)整過程參數(shù)進(jìn)行仿真分析。

本文利用有限元軟件ANSYS建立了軸向磁通動態(tài)感應(yīng)加熱模型。模型考慮了工件的實際三維幾何結(jié)構(gòu)，工件與感應(yīng)器的相對運動，材料屬性與溫度的非線性依賴關(guān)系，熱輻射邊界條件等問題。該模型還引入了kringing插值算法構(gòu)建工件初始溫度場，使模型更加符合實際。

1 數(shù)學(xué)模型及邊界條件

1.1 電磁場數(shù)學(xué)模型

感應(yīng)加熱過程電磁場可以由麥克斯韋方程微分形式表示：

其中：

為了避免求解復(fù)雜的一階偏微分方程組，往往引入位函數(shù)作為間接變量。常用的方法是引入矢量磁位A和標(biāo)量電位φ，它們的定義為：

得到渦流區(qū)的控制方程為：

為了確保矢量磁位A的唯一性，應(yīng)該對A施加庫侖規(guī)范即：

1.2 溫度場數(shù)學(xué)模型

熱傳遞問題遵循能量守恒定律：熱能或者機(jī)械能通過邊界進(jìn)入系統(tǒng)的速率減去能量通過邊界離開系統(tǒng)的速率，再加上系統(tǒng)體內(nèi)產(chǎn)生能量的速率，必須等于系統(tǒng)存儲能量的速率，其方程[Book]為：

直角坐標(biāo)系下微元體應(yīng)用能量守恒定律，得到熱傳導(dǎo)方程如下所示：

邊界條件：高溫下通過輻射損失的熱量在低速感應(yīng)加熱形式下是相當(dāng)可觀的。大多數(shù)金屬在溫度為850℃以上時，通常產(chǎn)生的熱損失很小。當(dāng)溫度在850℃以上時，熱損失變得顯著[]。熱輻射由下式給出：

其中，ε為表面輻射系數(shù)，σ為斯蒂夫—玻爾茲曼常數(shù)，T1為表面絕對溫度K(℃+273)，T2為環(huán)境絕對溫度K(℃+273)，Qe為每平方米表面積上的輻射熱損失。感應(yīng)加熱的工況空氣流動緩慢，熱對流損失相對于輸入功率很小，在模型中忽略這部分損失。

1.3 Kriging模型

回歸多項式是l個可選函數(shù)的線性組合。

其中R(θ,xi,xj)是參數(shù)為θ的相關(guān)模型。

在預(yù)測點x處定義：

以及隨機(jī)過程的相關(guān)矩陣R：

x處的響應(yīng)值y(x)預(yù)測估計值的表達(dá)式為[]：

其中：

文獻(xiàn)[9]提供了一種Matlab的Kriging工具箱。它提供了三種回歸模型，分別為0、1和2階多項式，六種相關(guān)模型，常用的有指數(shù)、球形和高斯模型。本文中回歸模型選擇1階模型regploy1，相關(guān)模型選擇指數(shù)模型，相關(guān)參數(shù)初值為θ= [1,1,1]。

圖1 幾何模型

2 有限元模型

問題描述：長度為3m的方形高溫鋼坯通過1m長的單感應(yīng)器進(jìn)行加熱，鋼坯的傳送速度為5m/min。線圈為方孔形，內(nèi)孔邊長為0.46m，外孔邊長為0.51m，匝數(shù)為30。三維幾何模型的1/4模型如圖1所示。

實現(xiàn)鋼坯的運動：

要實現(xiàn)鋼坯的運動，一般有兩種方法，第一種是重新劃分線圈網(wǎng)格來實現(xiàn)鋼坯運動[]，但是這種方法同樣十分耗時，更重要的是重新劃分的網(wǎng)格可能產(chǎn)生很大的畸變導(dǎo)致結(jié)果不準(zhǔn)確甚至運算終止。第二種方法是預(yù)先將線圈運動路徑上的的區(qū)域分別定義且劃分好，將線圈和空氣的電磁材料屬性來回更替，相當(dāng)于線圈在相對于空氣運動。本文基于方法二，利用載荷的移動來實現(xiàn)線圈與鋼坯的相對運動如圖2所示。先將線圈的運動路徑全部預(yù)先劃分網(wǎng)格，在t時刻，第一對箭頭包含的線圈區(qū)域加載電流載荷，在t +△t時刻第二對箭頭包含的線圈區(qū)域加載電流載荷。載荷運動的距離恰好為鋼坯傳送的距離△s。劃分的段數(shù)越多，△s越小計算的精度就越高，相對的計算時間也越長。分段數(shù)可以利用試錯法得到，將仿真結(jié)果與預(yù)分析結(jié)果或者實際數(shù)據(jù)逐漸逼近，在精度與計算時間之間平衡得到一個合適值，本文通過這種方法設(shè)定的分段數(shù)為25。這里需要對模型作如下假設(shè)：（1）忽略線圈每匝之間的間隙；（2）忽略線圈運動路徑上線圈材料與空氣電阻率和磁導(dǎo)率不同產(chǎn)生的影響。

圖2 載荷移動示意圖

Ansys提供了多種耦合場分析方法，本文采用的是順序耦合物理分析，即順序的求解物理場，并將前一次求解結(jié)果作為下一次分析的載荷，F(xiàn)EM求解流程圖如圖3所示。前面提到初始溫度場是不均勻的，指定各個節(jié)點的初始溫度后，需要一個穩(wěn)態(tài)熱分析來建立初始溫度場，設(shè)置很短的時間長度，如0.001s。各個節(jié)點的溫度設(shè)定，根據(jù)實際情況需要頻繁的使用和修改，可以利用ANSYS創(chuàng)建宏文件。Kriging的預(yù)報函數(shù)可以計算鋼坯每個節(jié)點的溫度插值結(jié)果，但是在ANSYS中按編號逐個給幾萬個節(jié)點設(shè)定初始溫度顯然不可行。本文采用的方法是，先由MATLAB工具包dace建立Kriging模型，模型中包括式(18)中的矩陣和規(guī)范化后的測點矩陣[s]，在宏文件中用數(shù)組保存這些矩陣。接著在宏中獲得

一個節(jié)點ω的坐標(biāo)，再根據(jù)選擇的回歸模型和相關(guān)模型的形式求得f(ω)T和r(ω)T，最后利用公式(18)求得溫度值。本文建立的有限元模型如圖4所示。

圖3 FEM求解流程圖

圖4 (a) 有限元模型

圖4 (b) 網(wǎng)格劃分

3 實例與分析

實驗?zāi)M的工況為某鋼廠連鑄坯感應(yīng)加熱過程。選取的鋼坯為是45號鋼，幾何參數(shù)如圖5所示，其參數(shù)與溫度的非線性關(guān)系如圖5所示，加熱一次時間為50.4s。

圖5 45號鋼在不同溫度下的物理參數(shù)：a電阻率；b相對磁導(dǎo)率；c比熱容；d熱導(dǎo)率

初始溫度場：

實際鋼坯的溫度分布呈現(xiàn)“頭低尾高”的現(xiàn)象，頭部溫度比尾部溫度整體大致高100℃。

初始溫度插值結(jié)果，如圖6所示：

圖6 初始溫度插值結(jié)果

實驗一：梯形電流

圖7(a)是實際該鋼廠鋼坯通過兩個感應(yīng)器的電源輸出電流曲線（保持端電壓不變），可見電源電流整體呈現(xiàn)為梯形，且左右對稱。原因是在鋼坯進(jìn)入第一個感應(yīng)線圈的過程中，系統(tǒng)阻抗從空載到滿載，當(dāng)感應(yīng)器中完全包含鋼坯的時候，雖然鋼坯還在前進(jìn)，但是系統(tǒng)的阻抗保持不變，電流也維持恒定。鋼坯進(jìn)入第二個感應(yīng)器的時候，由于鋼坯的長度大于兩個感應(yīng)器與其間隙的長度和，所以會出現(xiàn)電流曲線兩次上升的情況，鋼坯出感應(yīng)器的情況同理。因此，對于單感應(yīng)器加熱的情況，本文設(shè)計了如圖7(b)所示的的電流載荷曲線。需要指出，圖7(b)的電流曲線指的是每匝線圈流過的電流，它與電源電流存在如下關(guān)系：

圖7 (a) 實際電源輸出電流曲線

圖7 (b) 實驗一感應(yīng)線圈電流曲

鋼坯的溫度分布即是空間坐標(biāo)的函數(shù)又是時間的函數(shù)，為了更好的觀察實驗結(jié)果，圖8中左邊定義了1/4模型上的四條路徑R1、R2、R3、RZ，R1為鋼坯中心到鋼坯表面沿y方向的路徑，R2為鋼坯表面中心到鋼坯角部沿x方向的路徑，R3為鋼坯中心到角部沿對角線方向的路徑，RZ表示從鋼坯表面中部從頭部到尾部的路徑。右邊定義了1/8截面上的六個坐標(biāo)點A、B、C、D、E、F。

圖8 后處理路徑與截面節(jié)點

圖9將R1、R2、R3路徑(z=1.5m截面)和RZ路徑在加熱前后的溫度曲線作了對比，R1主要反映了感應(yīng)加熱主要集中在表面，導(dǎo)致表面溫度大幅上升，另外，隨著加熱時間增加，表面中心和邊角的溫度差也逐漸增加。從R2可以看出，初始時刻，中心溫度比表面溫度高，由于感應(yīng)加熱集中在表面，表面溫度迅速上升，最終超過中心溫度。同時，經(jīng)過感應(yīng)加熱表芯溫差大大縮小，如果均溫一段時間，表芯溫差能夠更小。從R3看出，感應(yīng)加熱讓鋼坯橫截面溫度分布更加均勻，鋼坯最熱的部分在內(nèi)表面的部分，通過均溫可以讓溫度分布更加均勻。從RZ路徑看出，初始溫度是呈現(xiàn)前低后高趨勢的，在梯形電流作用下，加熱頭尾溫度也同樣呈現(xiàn)前低后高的總體趨勢，同時還有明顯的端部效應(yīng)。圖9為橫截面結(jié)點溫度歷程圖，圖10為不同時刻鋼坯感應(yīng)加熱圖。

圖9 加熱前后各路徑上的溫度曲線圖

圖10 橫截面結(jié)點溫度歷程圖

圖11 不同時刻鋼坯感應(yīng)加熱圖

實驗二：優(yōu)化電流

通過實驗一，發(fā)現(xiàn)實際中感應(yīng)加熱結(jié)束后，溫度的不均勻性是普遍存在的。實際中一般設(shè)置均溫段對鋼坯進(jìn)行均熱，但是這種方法對于縮小截面溫差比較有效，而對于縮小軸向溫差效果不明顯，原因是鋼坯長度很長，軸向高溫區(qū)到低溫區(qū)的距離遠(yuǎn)，由傅里葉定律：

軸向熱流密度低，所以均熱效率也低。因此在加熱過程就應(yīng)該加入控制手段。本文希望通過調(diào)節(jié)線圈電流來實現(xiàn)縮小溫差的目標(biāo)，設(shè)計了如圖12所示的電流曲線。

實驗結(jié)果如圖13所示，與圖9比較可以發(fā)現(xiàn)R1、R2、R3的變化規(guī)律幾乎是一致的，只有RZ路徑上的溫度分布有了明顯的變化，鋼坯中間段的溫度分布更加均勻，這說明鋼坯表面溫度與線圈電流之間存在某種線性關(guān)系。同時，兩個實驗表明，這種加熱方式，軸向感應(yīng)加熱的端部效應(yīng)不可避免，原因可能是渦流在端部聚集，導(dǎo)致電流密度很高，從而產(chǎn)生很高的熱量。所以，實際中應(yīng)當(dāng)適當(dāng)控制電流，不能讓鋼坯的端部溫度過高。

圖12 實驗二感應(yīng)線圈電流曲線

圖13 加熱前后各路徑上的溫度曲線圖

4 結(jié)論

本文利用有限元軟件ANSYS建立了軸向磁通動態(tài)感應(yīng)加熱模型。通過引入了kringing插重建了連鑄坯進(jìn)入感應(yīng)加熱器前的初始溫度場。在設(shè)定的初始溫度條件下，預(yù)測了在梯形電流產(chǎn)生的電磁場作用下連鑄坯感應(yīng)加熱過程的溫度分布情況。結(jié)果表明，該模型能夠立體的反映了鋼坯感應(yīng)加熱的過程的溫度變化情況。

通過調(diào)整電流曲線，能夠消除連鑄坯初始的“頭低尾高”現(xiàn)象，發(fā)現(xiàn)表面溫度與電流之間存在某種線性關(guān)系，這有利于研究溫度的實時控制或者前饋控制。

該模型不但能仿真連鑄坯感應(yīng)加熱過程，對于冷坯、熱軋中間坯等的加熱過程，只要對模型稍做修改，同樣可以模擬。

[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—），浙江溫州人，碩士，現(xiàn)就讀于杭州電子科技大學(xué)。