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

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

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

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

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

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

0 引言

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

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

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

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

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

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

其中：

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

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

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

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

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

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

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

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

1.3 Kriging模型

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

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

在預(yù)測(cè)點(diǎn)x處定義：

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

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

其中：

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

圖1 幾何模型

2 有限元模型

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

實(shí)現(xiàn)鋼坯的運(yùn)動(dòng)：

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

圖2 載荷移動(dòng)示意圖

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

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

圖3 FEM求解流程圖

圖4 (a) 有限元模型

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

3 實(shí)例與分析

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

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

初始溫度場(chǎng)：

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

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

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

實(shí)驗(yàn)一：梯形電流

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

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

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

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

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

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

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

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

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

實(shí)驗(yàn)二：優(yōu)化電流

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

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

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

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

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

4 結(jié)論

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

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

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

[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é)。