大型工件感應(yīng)加熱的有限差分模型分析

王秀蓮，曹 宏，周光遠(yuǎn)，藤金輝

(沈陽(yáng)理工大學(xué) 1.信息科學(xué)與工程學(xué)院;2.經(jīng)濟(jì)管理學(xué)院，遼寧 沈陽(yáng) 110159)

感應(yīng)加熱過(guò)程具有加熱速度快、環(huán)境污染少，易實(shí)現(xiàn)自動(dòng)化生產(chǎn)等特點(diǎn)。感應(yīng)加熱過(guò)程的數(shù)值分析，為實(shí)際加工過(guò)程中工藝參數(shù)的合理選取和工藝線路的確定提供了理論依據(jù)，大大降低了生產(chǎn)成本。

感應(yīng)加熱計(jì)算方法有積分方程法、邊界元法、有限元法和有限差分法等［1－3］。航空工業(yè)部第五設(shè)計(jì)院干肇智［1］用有限差分方法對(duì)感應(yīng)加熱工件溫度分布進(jìn)行模擬，但其有限差分方法求出的電磁場(chǎng)和溫度場(chǎng)的迭代矩陣都是非正定的，程序編制的工作量大;杭國(guó)平［2］等人采用有限元法計(jì)算感應(yīng)加熱溫度分布，但其離散方法復(fù)雜，計(jì)算成本高;張媛媛［3］采用邊界元法進(jìn)行感應(yīng)加熱分析，其基礎(chǔ)是積分方程控制誤差困難。本研究來(lái)源于圓柱體型工件熱處理生產(chǎn)實(shí)際，在工件熱處理過(guò)程中，利用有限差分模型，并基于 VS2008軟件C++計(jì)算程序，計(jì)算感應(yīng)加熱過(guò)程中溫度分布，實(shí)時(shí)控制感應(yīng)加熱過(guò)程，熱處理溫度達(dá)到要求，可及時(shí)停止感應(yīng)加熱處理。程序簡(jiǎn)單有效，計(jì)算速度快，易于對(duì)加熱過(guò)程進(jìn)行控制。

1 模型建立

1.1 感應(yīng)加熱原理

感應(yīng)加熱實(shí)質(zhì)是利用電磁感應(yīng)在導(dǎo)體內(nèi)產(chǎn)生的渦流發(fā)熱加熱工件過(guò)程。當(dāng)感應(yīng)加熱的線圈中通過(guò)交變電流時(shí)，在線圈內(nèi)部和其周?chē)a(chǎn)生一個(gè)交變磁場(chǎng)，感應(yīng)加熱時(shí)，工件被交變的磁力線切割。根據(jù)電磁場(chǎng)理論，變化著的磁場(chǎng)產(chǎn)生感應(yīng)電動(dòng)勢(shì)，感應(yīng)電動(dòng)勢(shì)在工件表面形成封閉的電流回路，通常把這種電流稱(chēng)為渦流。渦流在工件內(nèi)部產(chǎn)生熱量，使工件得到加熱。

1.2 滲透層

由于電磁場(chǎng)的集膚效應(yīng)，感應(yīng)加熱過(guò)程中存在電流透入層，感應(yīng)電流集中于該層，提供感應(yīng)加熱大部分能量，決定了最終的加熱效果，這個(gè)深度便是感應(yīng)加熱的“滲透層深度”。渦流的理論滲透層深度 h(m)為［4］

式中:ρ為工件材料的電阻率(Ω·m);μr為工件材料的相對(duì)導(dǎo)磁率;μ0為真空磁導(dǎo)率，μ0=4π×10－7H/m;f為電流的頻率(Hz)。

1.3 模型建立

大型工件在加熱過(guò)程中其溫度場(chǎng)隨時(shí)間變化而變化。假設(shè)工件沿不同方向的導(dǎo)熱系數(shù)為常數(shù)，采用顯式有限差分方法，導(dǎo)出感應(yīng)加熱的有限差分模型。取平行六面體微元做控制體，應(yīng)用熱力學(xué)第一定律［5］，可得關(guān)于溫度T的通用導(dǎo)熱方程

式中:T為x、y、z和時(shí)間t的函數(shù);λ為導(dǎo)熱系數(shù);ρ為密度;c為比熱容;q為內(nèi)熱源。將其轉(zhuǎn)換至圓柱體坐標(biāo)系下，可得圓柱體坐標(biāo)系下的導(dǎo)熱方程

考慮圓柱體加熱時(shí)的徑向?qū)ΨQ(chēng)性，可按二維(z，r)不穩(wěn)定圓柱體導(dǎo)熱模型處理，得導(dǎo)熱模型如下:

式中α=λ/ρc為熱擴(kuò)散系數(shù)。將圓柱體按軸向和徑向離散化，劃分成不同的離散單元?？紤]T在(β，τ)點(diǎn)處的偏導(dǎo)數(shù)，使用前向有限差分［4］可得有源非穩(wěn)態(tài)二維圓柱體坐標(biāo)下導(dǎo)熱的顯式差分模型:

式中:i=1，2，…，(nz－1)，nz為圓柱體軸向離散化后的分段數(shù);j=1，2，…，(nr－1)，nr為圓柱體徑向離散化后的分層數(shù);n表示當(dāng)前時(shí)刻，(n+1)表示下一時(shí)刻;Δz、Δr代表離散單元的軸向和徑向長(zhǎng)度。

2 模型分析

2.1 模型計(jì)算的邊界處理

在感應(yīng)加熱過(guò)程中其邊界對(duì)流對(duì)溫度計(jì)算的影響遠(yuǎn)小于邊界輻射對(duì)溫度計(jì)算的影響。因此，在模型計(jì)算過(guò)程中忽略對(duì)流對(duì)溫度計(jì)算的影響，通過(guò)調(diào)整輻射系數(shù)的方法進(jìn)行加權(quán)，可減少邊界計(jì)算約束條件，又能保證計(jì)算的準(zhǔn)確性;對(duì)于邊界輻射，則采用斯忒潘－波爾茲曼定律［6］

進(jìn)行計(jì)算，式中的環(huán)境溫度To采用工件表面溫度加權(quán)來(lái)實(shí)時(shí)更新，保證計(jì)算的準(zhǔn)確性。

2.2 滲透層對(duì)渦流的影響

研究表明，在感應(yīng)加熱過(guò)程中，從工件表面到中心，感應(yīng)電流的大小按照指數(shù)規(guī)律下降和變化。距離表面x處的強(qiáng)度為［7］

式中:Io為表面的渦流強(qiáng)度(A/m2);x為距離工件表面的距離(m);h為滲透層深度(m)。分析可知，滲透層深度主要由相對(duì)磁導(dǎo)率和頻率［8］決定。當(dāng)頻率一定的情況下，相對(duì)磁導(dǎo)率成為影響滲透層深度主要因素;當(dāng)離散點(diǎn)溫度高于居里點(diǎn)溫度［9］(750℃)時(shí)，其相對(duì)磁導(dǎo)率發(fā)生突變，滲透層深度變大，渦流影響范圍變大。程序中使用簡(jiǎn)化計(jì)算，當(dāng)離散單元的溫度高于居里點(diǎn)溫度時(shí)，即把該離散單元加入到渦流影響的范圍內(nèi)，統(tǒng)計(jì)全部有渦流的離散單元，將所有能量均分到每個(gè)離散單元中，這樣既能保證計(jì)算的簡(jiǎn)單，又能保證計(jì)算精度。

2.3 模型計(jì)算在加熱過(guò)程中的作用

在感應(yīng)加熱過(guò)程中，模型計(jì)算起主要的控制作用。根據(jù)已知的初始溫度和加熱電源的總功率，計(jì)算出整個(gè)加熱過(guò)程所需數(shù)據(jù)，包括跟隨溫度變化的工件材料電阻率、磁導(dǎo)率、比熱和密度等物性參數(shù)［10］，為準(zhǔn)確計(jì)算下一時(shí)刻的溫度提供了有效保障，直至整個(gè)加熱過(guò)程完畢，得到加熱過(guò)程的溫度分布情況。

3 模型計(jì)算結(jié)果分析與對(duì)比

以尺寸 φ1630mm×2050mm圓柱體工件為例，工頻電源額定功率P=4000kw，頻率為50Hz;磁場(chǎng)感應(yīng)效率η=0.8～0.9;圓柱體初始溫度為300℃;離散化:軸向每段100mm均分，多余部分分別加在兩端;徑向(由表層向圓心)按照外6層每層10mm，中間3層每層20mm，其余為一層的離散方法。滲透層深度隨溫度的變化而變化，在10～50mm之間非線性變化;保證模型計(jì)算穩(wěn)定，取最小時(shí)間步長(zhǎng)Δt=2s。按工藝要求，以徑向距離表面70mm、軸向距離底面 1025mm、溫度達(dá)到800℃時(shí)，程序退出計(jì)算。模型計(jì)算過(guò)程利用VS2008環(huán)境下的C++程序完成，得到的計(jì)算結(jié)果用MATLAB繪圖分析。

在離散化中采用軸向均分、徑向上在表面和目標(biāo)層單元格分多層的方法，因?yàn)橹恍杩刂颇繕?biāo)層單元和表面溫度，對(duì)靠近圓心的溫度計(jì)算精度要求不高，這樣劃分既能滿足工藝要求，也能簡(jiǎn)化計(jì)算，加快計(jì)算速度。

3.1 有限元方法的介紹

以變分原理為基礎(chǔ)建立起來(lái)的有限元法是將整個(gè)區(qū)域分割成很多小的子區(qū)域，并構(gòu)造定義在子區(qū)域上的簡(jiǎn)單嘗試函數(shù)，將整個(gè)區(qū)域中各個(gè)子區(qū)域所對(duì)應(yīng)的嘗試函數(shù)線性組合起來(lái)，便形成近似解的表達(dá)式，完成計(jì)算的方法［2］。

3.2 模型計(jì)算速度分析

在滿足工藝計(jì)算要求的情況下，由前面的分析和離散化可看出，有限差分模型的離散化簡(jiǎn)單，計(jì)算選取的時(shí)間步長(zhǎng)既能保證模型的計(jì)算速度，又能保證模型計(jì)算的穩(wěn)定性;模型計(jì)算過(guò)程中對(duì)邊界條件和感應(yīng)渦流分布的處理方法也簡(jiǎn)化了計(jì)算，提高了模型的計(jì)算速度;模型計(jì)算過(guò)程中將電磁場(chǎng)與溫度場(chǎng)分開(kāi)計(jì)算，減少兩個(gè)場(chǎng)之間的耦合。而有限元方法采用ANSYS軟件進(jìn)行復(fù)雜的電磁場(chǎng)和溫度場(chǎng)計(jì)算，其耦合計(jì)算，計(jì)算過(guò)程復(fù)雜，計(jì)算速度和計(jì)算成本比較高，不能較好地滿足工業(yè)現(xiàn)場(chǎng)快速完成計(jì)算的要求。

3.3 結(jié)果分析與對(duì)比

利用有限差分程序使用MATLAB模擬出整個(gè)圓柱體中每個(gè)離散單元的溫度分布，如圖1所示，縱向表示圓柱體徑向方向上的溫度分布，橫向表示圓柱體軸向方向上的溫度。

圖1 圓柱體整體溫度分布圖

圖2是有限差分圓柱體三個(gè)離散單元在感應(yīng)加熱過(guò)程中的溫度上升情況，圖中說(shuō)明部分表示在圓柱體離散單元中的坐標(biāo)(x，y)，x表示徑向方向，y表示軸向。

圖2 溫度上升過(guò)程圖

圖3 溫度對(duì)比圖

圖3是在相同的離散情況下分別利用有限差分算法、有限元方法和實(shí)際測(cè)溫點(diǎn)測(cè)得圓柱體一個(gè)離散單元在感應(yīng)加熱過(guò)程中溫度變化對(duì)比圖，其中誤差一表示實(shí)際溫度與差分算法的誤差，誤差二表示實(shí)際溫度與有限元算法的誤差。

4 結(jié)論

(1)感應(yīng)加熱過(guò)程中圓柱體表面溫度上升速度較快且幅度大，而內(nèi)層單元格溫度上升速度較慢且幅度小，越靠近圓柱體圓心這種現(xiàn)象越明顯，驗(yàn)證了在感應(yīng)加熱過(guò)程中由于集膚效應(yīng)的存在使感應(yīng)電流集中于圓柱體表面，內(nèi)層則靠導(dǎo)熱來(lái)進(jìn)行溫度傳遞。

(2)在滿足實(shí)際工程要求的前提下，采用有限差分模型計(jì)算感應(yīng)加熱的溫度分布情況比采用有限元方法離散化簡(jiǎn)單，使用程序易于實(shí)現(xiàn)，計(jì)算程序簡(jiǎn)單有效，計(jì)算速度快、成本較低。

(3)在感應(yīng)加熱過(guò)程中可通過(guò)調(diào)整電源的功率使圓柱體表面溫度最終趨于穩(wěn)定，不會(huì)為達(dá)到加熱目的使表面溫度過(guò)高損傷工件屬性。

(4)利用模型的數(shù)值解有效分析了感應(yīng)加熱過(guò)程中溫度場(chǎng)的分布特性，為合理選取工藝參數(shù)提供依據(jù);為類(lèi)似感應(yīng)加熱熱處理過(guò)程計(jì)算提供簡(jiǎn)單有效的計(jì)算方法，并為控制感應(yīng)加熱過(guò)程提供有益參考。

(5)在相同的前提條件下，采用相同的離散化方法，采用有限差分方法比采用有限元方法計(jì)算感應(yīng)加熱溫度分布更接近實(shí)際溫度分布情況，且滿足工程精度要求，可指導(dǎo)現(xiàn)場(chǎng)生產(chǎn)實(shí)踐。

［1］干肇智.感應(yīng)加熱的數(shù)值計(jì)算［J］.電爐，1985(5):1－12.

［2］杭國(guó)平.感應(yīng)加熱問(wèn)題的有限元分析［D］.上海:上海大學(xué)，1991.

［3］張媛媛.基于有限元和邊界元方法的軸類(lèi)感應(yīng)加熱分析及數(shù)值模擬［D］.天津:天津大學(xué)，2006.

［4］趙長(zhǎng)漢，姜士林.感應(yīng)加熱原理與應(yīng)用［M］.天津:天津科技翻譯出版公司，1992.

［5］D.皮茨.L西索姆.全美經(jīng)典－傳熱學(xué)［M］.葛新石等譯.北京:科學(xué)出版社，2002:1 －9，69 －78.

［6］陶文銓.數(shù)值傳熱學(xué)［M］(第二版).西安:西安交通大學(xué)出版社，2001.

［7］金曉昌.感應(yīng)加熱技術(shù)中的趨膚效應(yīng)［J］.武漢化工學(xué)院學(xué)報(bào)，1995，17(4):65 －68.

［8］劉白.雙頻感應(yīng)淬火的計(jì)算與應(yīng)用［J］.貴州工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2001，30(4):20 －24.

［9］Fireteanu，Tudorache T.Electromagntic forces in transverse flux induction heating［J］.IEE Transactions on Magneties，2000，36(4):1792 －1795.

［10］Kawase Y，Miyatake T.Thermal analysis of steel blade quenching by induction heating［J］.IEEE Transactions on Magnetics，2000，36(4):1788 －1791.