鋼板移動式感應加熱的多場耦合數(shù)值分析

張雪彪，陳誠，劉玉君，張立衛(wèi)

(大連理工大學船舶工程學院工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧大連116024;2.大連理工大學數(shù)學科學學院，遼寧大連116024)

在船廠，水火彎板工藝廣泛應用于雙曲率板的成形過程中。眾所周知，火焰熱源的水火彎板工藝強烈依賴熟練工人的熱源控制技術(shù)，以達到控制鋼板變形形狀的目的。然而，火焰熱源是乙炔氣體和氧氣發(fā)生化學反應的熱源，精確控制火焰作用在鋼板表面的熱通量是困難的;同時，整個燃燒過程釋放大量的CO2氣體和噪音，對環(huán)境造成污染。為了克服火焰熱源的這些缺陷，人們開始研究采用感應熱源的水火彎板工藝:感應加熱直接利用渦流對鋼板進行加熱，熱效率高;通過調(diào)節(jié)頻率、功率等電學參數(shù)可以實現(xiàn)對熱輸入量的精確控制;最后，感應加熱工作環(huán)境好、噪音小，有利于環(huán)保。

盡管感應加熱具有許多優(yōu)點，由于缺乏足夠的加熱機理的理論分析以及變形特點的深入探索，此項工藝技術(shù)在船體板成形領(lǐng)域仍未被廣泛采用。目前，國內(nèi)外已經(jīng)有學者做了初步研究。例如，LEE J H［1］研究了感應加熱在板材成型中的可行性，指出火焰加熱和感應加熱具有很大的相似性。帥克剛、羅宇［2-3］依據(jù)數(shù)值與實驗結(jié)果的比較建立了感應加熱的感應熱模型，在數(shù)值計算時直接加載該熱源模型進行熱-結(jié)構(gòu)的計算。這一方法缺少電磁-熱耦合的分析過程，忽略了感應加熱過程中的集膚效應、端部效應等感應加熱的基本特點。LEE K S［4］利用商業(yè)軟件 JMAG-Studio開發(fā)了鋼板移動加熱的電磁-熱-結(jié)構(gòu)耦合模型，分析了加熱過程中的溫度分布和彎曲行為。楊玉龍和張雪彪［5-6］也利用COMSOL軟件開發(fā)了鋼板移動式感應加熱的磁-熱耦合的數(shù)值模型，對鋼板加熱過程中的電磁場和溫度的耦合特點進行了研究。雖然人們對鋼板感應加熱已經(jīng)做了一些工作，但是對于鋼板在感應加熱過程中的溫度機理和變形機理仍需仔細研究。因為鋼板感應加熱依靠金屬內(nèi)部的感應渦流進行加熱，然后進行熱傳導，與主要依靠熱傳導的火焰熱源有本質(zhì)的不同?；诖?，本文在前面工作的基礎(chǔ)上建立了移動式感應加熱的電磁-熱-結(jié)構(gòu)耦合的數(shù)值模型，研究多場耦合條件下鋼板的溫度和變形機理，并利用鋼板感應加熱實驗驗證了數(shù)值模型的有效性。

1 數(shù)學模型

如圖1所示，鋼板移動式感應加熱系統(tǒng)主要由3部分組成:鋼板、感應器、Π型導磁體。感應加熱過程中，感應器在鋼板上方以一定速度由板內(nèi)向板邊進行移動加熱。加熱過程中保持鋼板與感應器之間的空氣間隙不變。

圖1 鋼板感應加熱示意圖Fig.1 Schematic of steel plate induction heating

1.1 感應加熱彎板成型的理論

描述電磁感應現(xiàn)象的控制方程是麥克斯韋方程組。因為感應加熱的電磁場是時諧電磁場［7-8］，所以其麥克斯韋方程組為

式中:H是磁場強度，J是電流密度，E是電場強度，B是磁感應強度，ω是正弦變化角頻率(ω=2πf，f是感應加熱的電流頻率)，D是電位移，ρ是電荷密度。

引入磁場矢勢A，滿足:

此外，B、H、J、E兩兩之間滿足下列關(guān)系式:

式中:μ是磁導率，σ是電導率。電場E可以表示為

其中，φ是標量電勢。將式(5)～(8)代入式(1)整理得電磁場的控制方程:

其中，Js是源電流密度。本文中采用的數(shù)值模型是磁矢量勢法的三維有限元模型?；谌S磁矢量勢法的邊界條件如圖2所示。通量垂直邊界條件:

通量平行邊界條件:

遠場邊界條件:

電勢邊界條件:感應器尾部截面電勢φ為零。

圖2 電磁邊界條件Fig.2 Electromagnetic boundary conditions

感應加熱過程中，渦流集中分布于工件表層，這種現(xiàn)象稱之為集膚效應，集膚層深度δ的計算公式:

受到渦流的作用后，依據(jù)焦耳-楞次定律，鋼板加熱區(qū)域內(nèi)感應渦流產(chǎn)生的熱量q。

式中:Je是感應渦流密度。鋼板內(nèi)部的熱傳導現(xiàn)象由微分控制方程(15)來描述:

式中:ρ是鋼板密度，Cp是鋼板的比熱，T是鋼板的溫度，k是鋼板的導熱系數(shù)。

鋼板感應加熱中需要考慮空氣冷卻、水冷卻及輻射換熱等換熱過程，所以鋼板的熱邊界條件是:

與空氣接觸的鋼板:

與水接觸的鋼板:

式中:qa、qw分別為空氣冷卻和水冷卻的熱損量，ε是放射率，σb是玻爾茲曼常數(shù)，T∞是環(huán)境溫度，T是鋼板壁面的溫度，A是在溫度T范圍內(nèi)的面積，h1是空氣的對流換熱系數(shù)，利用平板外部自然對流模型求解，h2是水的對流換熱系數(shù)，利用流動沸騰理論求解，h1、h2的計算過程參考文獻［10］。加熱時鋼板各表面施加空氣冷卻邊界條件(16)，加熱結(jié)束后進行正面水冷，在水冷區(qū)域施加水冷邊界條件(17)。

鋼板的塑性應變過程通過下述方程來描述［11］:

式中:σij，i=?σij/?xi，ui，j=?ui/?xj，uj，i=?uj/?xi，σij是應力張量，εij是應變張量，F(xiàn)i是外力矢量(本分析值中為零)，E是彈性模量，υ是泊松比，δij是克羅內(nèi)克函數(shù)，e是應變張量的對角線部分。

應變張量εij由3部分組成疊加而成:

1.2 數(shù)值模擬方法

鋼板感應加熱過程中，電磁場與溫度場的耦合是一個交互式的耦合過程。從式(16)可以看出，在焦耳-楞次定理的作用下，電磁分析得到的渦流值轉(zhuǎn)化為熱載荷q進行熱分析。隨著鋼板溫度的升高，鋼板的磁導率μ和電導率σ會發(fā)生改變，則在下一載荷步的電磁計算時必須更新這些電磁學參數(shù)，同理在進行下一載荷步的熱分析時，材料的熱物性參數(shù)也要更新;此外，感應器移動時，鋼板、感應器、空氣組成的求解域時刻在發(fā)生變化。上述兩點是建立電磁-熱-結(jié)構(gòu)耦合數(shù)值模型所需解決的問題。

本文采用ANSYS軟件的物理環(huán)境法進行鋼板電磁感應加熱的多場耦合數(shù)值計算。圖3是電磁-熱-結(jié)構(gòu)耦合計算的流程圖，首先設置電磁場計算的初始條件以及初始模型的約束方程，進行電磁諧波分析，計算得到的渦流值和初始鋼板溫度一起作為熱輸入量進行熱分析。熱分析后，判斷是否進行下一載荷步的電磁求解?如果需要進行下一載荷步的電磁場求解，則移動感應器模型，重新建立感應器和鋼板之間的約束方程，更新鋼板的材料屬性，進行新的電磁分析，從而實現(xiàn)電磁-熱物理場之間的交互式耦合;如果仍然需要進行溫度場的求解，則根據(jù)鋼板表面的溫度，更新熱物性參數(shù)，再次進行溫度場計算，直到熱分析結(jié)束。當電磁-熱耦合分析結(jié)束后，開始進行熱-結(jié)構(gòu)的計算，將電磁-熱耦合的溫度計算結(jié)果作為熱載荷作用于結(jié)構(gòu)模型，從而完成對鋼板變形場的求解。電磁和熱分析的邊界條件如前所述，結(jié)構(gòu)變形的邊界條件為四端簡支。

圖3 電磁－熱－結(jié)構(gòu)耦合計算流程圖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單元，結(jié)構(gòu)變形求解中鋼板模型使用SOLID 185單元。圖4(a)是鋼板感應加熱的整體模型，感應器處于空氣包中，圖4(b)是感應器和鋼板結(jié)合部分的局部網(wǎng)格圖。

為了建立感應器和鋼板之間的移動關(guān)系，感應器和鋼板之間的空氣間隙被分成兩層，空氣間隙的上層是移動的，空氣間隙的下層是靜止的，在兩層網(wǎng)格的結(jié)合處設置約束方程以進行感應器移動時的耦合計算。由于磁場強度在感應器周圍較高，所以在感應器附近區(qū)域的網(wǎng)格劃分較密集?？紤]到集膚效應，鋼板表層也劃分了密集的網(wǎng)格。

鋼板的電磁學參數(shù)和熱物性參數(shù)均與溫度相關(guān)［6］，力學性能數(shù)據(jù)來自參考文獻［12］。

2 鋼板感應加熱實驗

為了驗證數(shù)值模型的正確性，課題組在大連理工大學造船工藝實驗室采用鋼板感應加熱設備SPZ-70進行了實驗。

圖5是溫度和變形測量點的布置圖，P1點是溫度測量點，P1(620，150，18)點位于鋼板表面，加熱線的縱向中間位置，距離加熱線20 mm，該點的溫度利用Raytek公司的自動紅外測溫儀進行測量。橫向收縮量的測量則通過兩點距離的變化得出，即在加熱線兩側(cè)確定11對測量點，通過數(shù)顯游標卡尺(精0.02 mm)測量它們在加熱前后的距離值，從而得出橫向收縮量值。加熱時的工藝參數(shù)如表1所示，感應加熱現(xiàn)場如圖6所示，鋼板焰道在加熱后進行正面水冷。

表1 感應加熱實驗的工藝參數(shù)Table 1 Technical parameters of induction heating

圖5 實驗布置Fig.5 Experiment setup

圖6 感應加熱實驗圖Fig.6 Induction heating experiments

3 結(jié)果與討論

3.1 數(shù)值實驗數(shù)據(jù)對比

圖7是P1點溫度的實驗值和計算值的對比圖，在加熱階段兩者的平均相對誤差是6.86%。在34～70 s的時間段內(nèi)感應器經(jīng)過P1點區(qū)域，鋼板金屬受到感應渦流作用，溫度急速上升至550℃左右。當感應器移過P1點后，P1點溫度開始下降，下降速率較慢。111 s后進入冷卻階段，P1點的溫度在水冷作用下急速下降。P1點溫度計算值和實驗值的誤差在可接受的范圍內(nèi)，鋼板溫度的計算結(jié)果是正確的。

圖8是橫向收縮量計算值和實驗值的對比圖。兩者的變化趨勢相同，平均相對誤差為10.77%。圖11中在板邊處的收縮變形偏低，這主要是由感應器加熱到板邊后出現(xiàn)的端部效應導致的。

圖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 鋼板移動式感應加熱的電磁-熱-結(jié)構(gòu)分析

鋼板在移動式感應加熱過程中，鋼板的電磁場場分布、溫度分布是不斷變化的，這與靜止式感應加熱有較大的區(qū)別。根據(jù)感應器的移動位置，本文從初始加熱位置開始確定了6個典型截面，對應截面的位置和時刻如圖9所示。

圖10是1、41、111 s 3個時刻鋼板上表面的磁通量密度云圖。如圖10(a)，t=1 s時鋼板表面上的磁通密度前后對稱分布;隨著加熱的進行，感應器不斷前移，如圖10(b)，磁通密度分布較初始狀態(tài)發(fā)生較大變化。由圖11可以看出，t=41 s時刻，感應器前端鋼板的渦流值(2.8×108A/m2)和磁通量密度(1.5 T)遠大于感應器中后方鋼板的渦流值(0.25×108A/m2)和磁通量密度(0.04 T)，由式(6)可知，B與μ是正比關(guān)系，感應器前端鋼板溫度較感應器中后方鋼板溫度低，相對導磁率大，所以此區(qū)域的磁通密度大，感應渦流大。綜上所述:鋼板移動式感應加熱過程中，感應渦流分布集中于感應器前端的鋼板區(qū)域內(nèi)，有利于感應器前側(cè)鋼板金屬的快速升溫。

圖9 電磁－熱－結(jié)構(gòu)分析中的鋼板截面示意圖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)是加熱到板邊處時的磁通密度云圖，由圖中可以看出在板邊處磁通密度擴張，呈現(xiàn)出較強的端部效應。圖12是4個時刻下鋼板上表面及厚度方向(如圖9所示)的溫度等值線圖，其中厚度方向的等值線圖做了放大處理，這些圖反映了鋼板溫度的動態(tài)變化過程。圖12(a)是1 s時刻的溫度等值線圖，鋼板表面溫度已經(jīng)達到331℃，厚度方向的加熱深度較淺;鋼板上表面的溫度上升很快，從圖12(b)的11 s時刻，鋼板最高溫度已升至900℃以上;自11 s以后，如圖12(c)的41 s時刻，鋼板最高溫度趨于穩(wěn)定，鋼板進入了準穩(wěn)態(tài)加熱階段;在圖12(d)的111 s時刻，鋼板上表面溫度等值線有明顯的外擴趨勢，與圖10(c)中的感應磁場相對應，形成了感應加熱特有的端部效應現(xiàn)象。

圖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

感應加熱過程中，鋼板能快速升溫并進入準穩(wěn)態(tài)加熱狀態(tài)的原因是:感應加熱起始的磁通密度較大，渦流值很大，產(chǎn)生大量焦耳熱，鋼板溫度快速上升，當鋼板表面溫度超過居里點溫度(720℃)后，集膚層內(nèi)的金屬失去磁性(相對導磁率趨于1)，磁通密度減少，渦流值減少，同時感應渦流向鋼板內(nèi)部滲透，鋼板表面的熱輸入量減小，鋼板的熱輸入量與鋼板內(nèi)部的熱傳導外加空氣換熱損失的熱量總和達到平衡，所以鋼板表面溫度上升速率減緩，趨于平穩(wěn)。

圖13 是Y分別為0.18、0.06、0 m(截面位置如圖9所示)3個橫截面上表面節(jié)點分別在41、81、111 s時刻的X方向應力分布圖。

圖13 Sx應力分布圖Fig.13 Distribution of stress of Sx

各截面節(jié)點的Sx值均是負值，即加熱線兩側(cè)的節(jié)點均承受橫向壓應力，其中Y分別為0.18、0.06 m兩截面的Sx分布趨勢相似;然而，Y=0 m橫截面上的應力在遠離加熱線處呈現(xiàn)出較大的壓應力，與前面的截面溫度分布圖相對應，Y分別為0.18、0.06 m的兩截面均已處于準穩(wěn)態(tài)加熱，如圖12(c)所示，高溫區(qū)集中在加熱線附近;而Y=0 m截面，由于感應加熱端部效應的存在，如圖12(d)所示，鋼板端部高溫區(qū)外擴，應力分布也發(fā)生改變，應力最大值出現(xiàn)在遠離加熱線處。

4 結(jié)束語

本文基于Ansys有限元軟件采用耦合約束方程的方法開發(fā)了鋼板移動式感應加熱的多場耦合數(shù)值模型，實驗數(shù)據(jù)與計算結(jié)果吻合，證明了數(shù)值模型的正確性。利用電磁-熱-結(jié)構(gòu)耦合的數(shù)值模型分析了鋼板在線加熱過程中的渦流分布、鋼板的溫度分布和結(jié)構(gòu)變形;鋼板感應加熱速度快，迅速進入準穩(wěn)態(tài)加熱狀態(tài)，加熱時感應加熱的渦流集中于感應器前端，有利于冷卻金屬的迅速升溫;在感應器加熱至板邊時會產(chǎn)生端部效應，影響板邊的收縮變形。

［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］.鍋爐技術(shù)，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］楊玉龍，張雪彪，劉玉君.船體板高頻感應加熱的多場耦合數(shù)值分析［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］張雪彪，楊玉龍，劉玉君，等.鋼板彎曲工藝中的高頻感應加熱過程數(shù)值分析［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版.北京:機械工業(yè)出版社，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］汪驥.水火彎板自動化加工工藝的關(guān)鍵技術(shù)研究［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，F(xiàn)ELIACHI 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.