粉末注射成形材料微波燒結過程的有限元模擬

石建軍，成志強，GELIN J C，柳葆生，BARRIERE T

（1.西南交通大學力學與工程學院，成都610031；2.Femto-ST研究所應用力學實驗室，法國貝桑松25000）

0 引 言

微波加熱是將材料自身吸收的微波能轉化為材料內部分子動能和熱能的過程，具有整體加熱均勻、升溫快速、燒結效率高等優(yōu)點［1］。在微波加熱過程中，電磁波以波的形式滲透到介質內部，引起介電損耗發(fā)熱，熱量從材料內部產生，加之微波的高頻振動，使燒結過程中材料內部的溫度梯度可以減小到最小程度。同時在微波電磁能作用下，材料內部分子或離子的動能增大，擴散程度增加。這些特性對于提高材料的微觀、亞微觀均勻性是非常有益的［2］。微波加熱還具有選擇性加熱［3］和無污染性［4］等特點，因此得到了廣泛應用。

微波燒結工藝結合了粉末冶金與微波加熱技術，近幾年被廣泛研究。如，晉勇［5］等采用微波燒結新技術研究了納米金屬陶瓷材料的燒結工藝與性能；王美娜等［6］利用微波燒結技術制備了鈦酸鎂基介電陶瓷；Roy等［7］利用微波技術成功燒結制備了金屬粉體，伺候微波燒結技術被廣泛應用于制造多種標準金屬粉末成形件。作者團隊和法國Femto-ST實驗室合作，在粉末注射成形和燒結致密化的模擬和試驗方面進行了長期研究，并致力于粉末喂料的配方試驗［8］，標定PIM（粉末注射成形）過程中的物理參數［9］，建立雙向流填充的數學模型［10］并進行傳統(tǒng)熱燒結試驗和數值模擬分析［11］；最近研究團隊配置了微波燒結試驗設備，對PIM的微波燒結進行進一步理論分析和試驗探討。

微波燒結是一種多物理場的耦合，包括電磁場、熱傳遞和粉末燒結過程的耦合，只有綜合考慮這些復雜的耦合現象才能真實地模擬微波燒結過程?，F有的研究成果大多集中在對微波加熱過程的研究，如，文獻［12］運用FDTD方法模擬多模腔的微波燒結試驗，文獻［13］通過擬和手段分段逼近材料的物理性能，提出了簡單的半解析模型，文獻［14］建立了一種半解析模型研究微波燒結過程中的能量轉換。這些研究耦合了電磁場和熱傳遞的宏觀模型，但是都沒有涉及材料的致密化過程。鑒于此，作者應用微波加熱的相關理論，借助Zhang等［15］提出的熱粘塑性本構模型，引入研究團隊成員利用熱膨脹計和傳統(tǒng)熱燒結爐標定的適用于多孔體注射成形坯件的宏觀燒結模型［16-17］，進一步完善了描述非鐵磁性陶瓷粉末成形件微波燒結全過程的數學模型和模擬方法；此外，還借助有限元求解器COMSOL Multiphysics實現了微波電場、熱傳遞和材料致密化過程的耦合計算，展示了單模腔燒結爐微波燒結過程的三維數值模擬結果。

1 微波燒結的數學模型

微波是通過電磁場能量在傳播過程中產生損耗的方式對物質進行加熱，由于其內部加熱的特性，微波損耗作為內熱源加熱被燒結物體。建立微波燒結過程的數學模型主要需要解決以下五個問題。

1.1 求解Maxwell方程得到加熱腔內的電場分布

電磁場問題的宏觀控制方程是Maxwell方程。

相對于電場效應而言，微波場中的磁效應對非鐵磁性陶瓷粉末成形件的加熱作用甚微［18］，故暫不考慮磁效應的影響。加熱腔內的電場分布通過式（1）［18］求取。

式中：μr為相對磁導率；E為電場強度；ω為微波源的角頻率；εr為相對復介電常數（加熱腔內的積分區(qū)域分為兩部分：燒結件占據部分和其它區(qū)域，這兩個區(qū)域具有不同的相對復介電常數）；為介電常數；為介電損耗因子，表征介質的損耗特性，是決定電場分布的重要參數。

在燒結過程中，燒結件內的電磁參數隨其密度和溫度的變化而變，因此，εr，μr等參數是燒結件密度和溫度的函數。這些參數的變化規(guī)律，目前尚無標定研究工作，根據文獻［19］取參照值進行模擬。

1.2 求解微波燒結過程中熱量的分布

材料在電磁場中生成的電阻熱即微波加熱過程中的損耗功率P，可按式（3）計算［19］。

式中：H 為磁場強度；σEC為電導率為被加熱物質的磁損耗因子（相對磁導率的虛部）。

由于可以忽略磁效應的作用，故式（3）右端第二項在求解中可被消除。

1.3 求解傳熱方程得到燒結件內部的溫度分布

微波加熱的原理是將被加熱物體作為微波傳輸的有耗介質，將有耗介質對微波的損耗變?yōu)闊崮?，使加熱物體的溫度升高。在加熱過程中，由于電磁場分布的不均勻性導致?lián)p耗功率分布不均勻，繼而造成了溫度分布的不均勻，所以在微波加熱過程中存在著熱量的傳遞過程，可通過求解傳熱方程（Flourier方程）得到溫度場的分布：

式中：ρ為密度；CP為常壓熱容；k為熱傳導系數；T為溫度；t為時間。

熱力學參數CP和k隨燒結過程中燒結件密度和溫度的變化而變，需要通過試驗測定，模擬中取文獻［19］中的參照值。

1.4 求解燒結過程的控制方程得到燒結件的結構響應

PIM材料注射成形及脫脂后的坯件由粉體和孔隙組成，其致密化過程是由高溫燒結過程中材料晶粒和晶界的擴散造成的。宏觀力學模型基于連續(xù)介質力學原理，以類粘塑性本構關系描述粉末的燒結致密化行為，同時考慮材料的彈性效應和熱膨脹效應。

燒結過程是一個熱力耦合現象，溫度場的演變控制燒結件的致密化進程。力學模型的控制方程組包括傳統(tǒng)的平衡方程、協(xié)調方程，這里不再一一列出，可參照參考文獻［20］。燒結過程模擬的特有方程還包括質量守恒方程和燒結本構方程。其中質量守恒方程為：式中為總應變率張量；tr（）為總應變率的跡；為密度變化率；為彈性應變率張量為熱應變率張量；為粘塑性應變率張量。

式中：σ為燒結件中分布的應力張量；σ′為應力偏張量；tr（σ）為應力張量的跡；σs為燒結應力；Ⅰ為二階單位張量；GP和KP分別為燒結過程中多孔體的剪切黏度模量和體積黏度模量。

GP、KP通過試驗標定燒結材料的單軸黏度和黏度泊松比確定；燒結應力σs通過不同的燒結溫度進程標定；燒結應力張量σ為密度和燒結粉末表面能和界面能的函數，通過熱膨脹儀對小試樣進行燒結標定［17］。在微波燒結模擬建模的框架下，這些標定參數還是溫度進程和密度的函數。

進行模擬還需確定燒結件的熱膨脹系數α，以及與彈性變形部分有關的材料參數。這些關系可以示意性地表達為：

α＝α（ρ，T）（8）

Ee＝Ee（ρ，T）（9）

ν＝ν（ρ，T）（10）

式中：Ee為多孔材料的彈性模量；ν為泊松比。

1.5 Maxwell方程、傳熱方程和材料本構方程的耦合

微波燒結的過程包含三個性質不同的物理過程，其一是微波能在加熱腔內的傳播過程，在傳播過程中由于被燒結材料的損耗，一部分微波能量轉化為熱能；其二是傳熱過程，由于加熱腔內電場分布的不均勻性，導致被燒結材料內部生熱不均勻，存在熱傳導的過程；其三是材料的燒結致密化過程，PIM材料顆粒在高溫狀態(tài)下會相互連接在一起，顆粒間的空隙逐漸減少或消失，而宏觀上則表現為燒結后產品外形收縮，體積減小，密度增大。同時，燒結密度和溫度的變化進程也將導致各物理方程的參數變化，形成復雜的耦合關系。

2 微波燒結過程的有限元模擬

2.1 模型定義

微波燒結是微波加熱實現燒結致密化的過程。頻域電磁分析和分布熱源涉及穩(wěn)態(tài)方程的求解；瞬態(tài)熱傳導模擬可獲得溫度在燒結件內的分布；燒結過程的力學模擬可獲得粉末成形件在燒結過程中密度和收縮變形的演化。

PIM材料微波燒結過程的數學模型在時間積分上采用增量算法。在同一時間增量內，分別對代表以上三種現象的控制方程進行分步求解。三種現象的耦合通過時間步的遞增實現，在每一時間步內，采用直接法求解該時刻的電場分布，溫度場和材料的致密演化需要進行迭代求解。在求解過程中，損耗功率的分布作為內熱源引入傳熱式（4），求解加熱腔內的溫度分布。將該溫度場導入PIM材料的致密化方程和本構方程中，求解得到PIM材料燒結過程的致密演化與變形值。

由于介質的電磁參數和熱力學參數隨材料密度和溫度的變化而變化，在增量求解中可利用更新后的數值求解當前時間步的各變量場。如此循環(huán)遞增，直到完成預定的燒結過程。所以整個微波燒結過程是Maxwell方程、傳熱方程和結構平衡方程/材料致密化方程/本構方程的交替數值求解過程，如圖1所示。

將微波爐內腔看成是一個金屬盒，其與一個2.45GHz、輸出功率連續(xù)可調的微波源相連。由于忽略微波中的磁效應，故可將微波源設置為一橫電（TE）波。爐腔和波導的壁面看成是完美的電導體，用邊界條件n×E＝0（n為垂直于邊界平面的單位矢量；E為電場強度矢量）表示。利用對稱性截取一半模型，由于電場關于對稱面鏡像對稱，所以對稱面可看成是完美的磁導體，用邊界條件n×H＝0（H為磁場強度矢量）表示。燒結模擬過程中要考慮重力效應。

模擬使用的陶瓷粉末成形件的物理參數取自參考文獻［19］中氧化鋯粉體的近似參數。微波源輸入功率設為恒定值1kW，初始溫度T0＝293K，初始相對密度ρr0＝0.65。燒結材料的復介電常數εr＝10-i0.1e0.0017（T-293），磁導率k＝1，密度 ρ＝600 0ρrkg·m-3（ρr為相對密度），熱傳導率k＝30ρrW·（m·K-1），熱容CP＝900ρrJ·kg-1·K-1，熱膨脹系數α＝1.2×10-2K-1，重力系數g＝9.8N·kg-1，體積黏度模量GP＝1×1011Pa，剪切黏度模量KP＝1×1011Pa，泊松比ν＝0.33。需要強調的是，由于試驗測試技術的局限性，類似于材料的相對復介電常數、電導率等電介質參數以及控制材料致密化過程的熱粘彈塑性本構關系，目前仍是材料科學領域需要研究的課題，從文獻查閱的相關數據也存在很大差異。因此，上述部分參數設置為初始參數，在目前的數值算法研究中暫時作為常量處理。為了準確真實地模擬微波燒結過程，提供真實的材料參數是至關重要的。而作者側重于微波燒結致密化過程數值模擬方法的研究，在下一步的研究工作中將通過試驗手段獲得符合微波燒結過程的盡可能真實的材料參數。

圖1 微波燒結過程模擬多物理場耦合的示意Fig.1 Abridged general view of coupling of multi-physics in simulation of microwave sintering

2.2 數值計算

對PIM材料微波燒結過程的數值模擬，需要聯(lián)合求解上述建立的數學模型中涉及的五個問題。這里借助有限元軟件COMSOL Multiphysics作為求解器，使用該軟件自帶的4個計算模塊：利用射頻模塊中的電磁波（“rfw”）求解電磁場的分布；利用廣義熱傳遞模塊的瞬態(tài)分析（“htgh”）求解溫度場分布；利用結構力學模塊中的實體、應力-應變分析（“smsld”）求解燒結件的結構響應；利用通用形式偏微分方程（PDE模式）的瞬態(tài)分析（“g”）求解質量守恒方程（5），更新瞬時密度衍變值。

為確定放置粉末成形件的適當位置，首先對微波燒結爐的空爐腔做電磁場分布分析，如圖2所示。由圖2可知，為了讓粉末成形件盡可能多地吸收電磁波的能量，應該將其置于電場峰值附近的區(qū)域，如圖3所示。這里采用形狀規(guī)則的長方體，可直接使用COMSOL軟件自帶的四面體自由網格劃分模型，通過調節(jié)網格尺寸參數，可控制不同區(qū)域的網格大小。燒結區(qū)域的有限元模型如圖4所示。為了簡化模型，假設粉末成形件外表面與燒結爐內腔環(huán)境熱絕緣。在進行熱傳遞、結構應力-應變和相對密度演化的求解時，抑制微波爐內爐腔部分的計算，而只計算粉末成形件內部材料的反應。

圖2 微波燒結爐空爐腔內的電場分布Fig.2 Electric field in empty cavity of microwave sintering furnace

圖3 粉末成形件放置的位置Fig.3 Placement of powder molding sample

圖4 COMSOL中的有限元模型Fig.4 FEM in COMSOL

從圖5中可以看到，經過1 500s之后，PIM燒結件內部溫度達到1 650K左右，峰值溫度位于燒結件內部，且其整體溫度梯度不大。從圖6，7中可以看到，由于熱效應和重力效應的共同作用，此時，PIM燒結件向中心和底部方向產生了比較明顯的收縮變形，這就是燒結的致密化現象。

圖5 加熱1 500s后PIM燒結件內部溫度場的縱向截面云圖Fig.5 Longitudinal section of temperature field inside PIM sintered sample after heating for 1 500s

圖6 加熱1 500s后PIM燒結件x方向的線位移云圖Fig.6 Displacement along x-direction of PIM sintered sample after heating for 1 500s

圖7 加熱1 500s后PIM燒結件的總位移云圖Fig.7 Total displacement of PIM sintered sample after heating for 1 500s

采用COMSOL非線性求解器進行求解，整個微波燒結過程耗時1 500s（25min），從0～1 500s均等截取6個時間值作為計數點保存對應時刻的結果。從圖8，9可以看到，由于復介電常數是溫度的函數，其損耗因子隨溫度的升高而增大，故在加熱的初始階段，PIM材料內部產生的電阻熱很小，加熱緩慢；在600s之前，材料內部的溫度和材料的相對密度無明顯變化，600s之后，PIM材料內部溫度持續(xù)急劇升高并達到峰值，其相對密度也隨著溫度的升高而不斷增加，最后接近完全致密材料。由于模擬采用固定的輸入功率（1kW），所以材料內部溫度會不斷升高，如此時調整微波源輸入功率，即可控制微波燒結進入保溫階段，隨后進行冷卻處理。將模擬結果與Charmond等［21］對納米級氧化鋯粉末成形件微波燒結的試驗結果進行對比分析。Charmond等［21］發(fā)現氧化鋯粉體在加熱溫度低于400℃時與微波的耦合作用較微弱，當溫度超過400℃時，微波的加熱效應逐步增強，這與模擬結果（圖8，9）十分吻合。由圖10可見，有限元模擬結果比微波燒結的試驗結果略小，這是因為采用的燒結模型并非經過標定的確切本構關系，所以得到的相對密度及其變化趨勢只具有參照意義。

圖8 粉末成形件中心產生的電阻熱-時間曲線Fig.8 Evolution of resistance heat generation in the center of PIM sample

圖9 PIM燒結件中心溫度和相對密度的變化曲線Fig.9 Variation curves of center temperature and relative density of PIM sintered sample

圖10 在1 360℃下采用不同方法制得燒結件相對密度與有限元模擬結果的比較Fig.10 Comparison of relative density of sintered samples prepared at 1 360 ℃ by different sintering processes with FEM result

3 結 論

（1）利用COMSOL Multi-physics有限元軟件模擬氧化鋯粉末成形件的微波燒結過程，模擬結果表明整個微波燒結過程耗時1 500s，初始階段由于粉末成形件內部溫度較低，粉體材料吸波能力較弱，加熱緩慢；當溫度升至400℃后，粉末成形件內部溫度持續(xù)急劇升高；當加熱至1 360℃時，燒結件的相對密度達到92%，滿足粉末燒結工藝要求。

（2）建立的數學模型和有限元手段能有效模擬微波燒結過程中PIM材料內部的電場、溫度場分布以及密度的衍化過程，為深入分析和認知微波燒結過程提供了有效工具。

（3）為了更準確模擬微波燒結過程，需要通過試驗手段標定粉末材料真實的物理參數與本構關系，同時應具體考慮邊界條件、重力以及摩擦力等外部影響的因素。

［1］陳鼎，李林，陳振華.金屬材料微波燒結的研究現狀［J］.機械工程材料，2012，36（4）：7-10.

［2］GUO Y，YI J，LUO S，et al.Fabrication of W-Cu composition by microwave infiltration［J］.Journal of Alloys and Compounds，2010，492：L75-L78.

［3］SHON I J，KIM B R.Initial stage sintering of binderless tungsten carbide powder under microwave radiation［J］.Ceramics International，2011，37（2）：505-512

［4］TOMPSETT G A，CONNER W C，YNGVESSON K S.Microwave synthesis of nanoporous materials［J］.Chemphyschem，2006，7（2）：296-319.

［5］晉勇，薛屺，湯小文，等.納米金屬陶瓷材料的微波燒結工藝研究［J］.機械工程材料，2004，28（12）：49-51.

［6］王美娜，丘泰，沈春英.鈦酸鎂基介電陶瓷的微波燒結［J］.機械工程材料，2009，33（1）：23-25.

［7］ROY R，AGRAWAL D，CHENG J，et al.Full sintering of powdered-metal bodies in a microwave field［J］.Nature，1999，399：668-670.

［8］KONG X，QUINARD C，BARRIERE T，et al.Mixing and characterisation of stainless steel 316Lfeedstock［J］.International Journal of Material Forming，2009，2（1）：709-712.

［9］QUINARD C，BARRIERE T，GELIN J C.Development and property identification of 316Lstainless steel feedstock for PIM andμPIM［J］.Powder Technology，2009，190（1/2）：123-128.

［10］LARSEN G，CHENG Z Q，BARRIERE T，et al.Modeling and numerical simulation of biphasic fluid flow［J］.Steel Research International，2010，81（9）：1458-1461.

［11］SONG J，BARRIERE T，LIU B，et al.Experimental and numerical analysis on the sintering behaviours of injection moulded components in 316Lstainless steel powder［J］.Powder Metallurgy，2010，53（4）：295-304.

［12］ISKANDER M D，ANDRADE A O N M.FDTD simulation of microwave sintering of ceramics in multimode cavities［J］.IEEE Trans Microwave Theory and Techniques，1994，42（5）：793-799.

［13］KOZLOV P V，RAFATOV I R，KULUMBAEV E B，et al.On modeling of microwave heating of a ceramic material［J］.Journal of Physics D：Applied Physics，2007，40（9）：2927-2935.

［14］LASRI J，RAMESH P D，SCHACHTER L.Energy conversion during microwave sintering of a multiphase ceramic surrounded by a susceptor［J］.Journal of American Ceramic Society，2000，83（6）：1465-1468.

［15］ZHANG B，GASIK M.Stress evolution in graded materials during densification by sintering processes［J］.Computational Materials Science，2002，25（1/2）：264-271.

［16］QUINARD C，SONG J，BARRIERE T，et al.Elaboration of PIM feedstocks with 316Lfine stainless steel powders for the processing of micro-components［J］.Powder Technology，2011，208（2）：383-389.

［17］SONG J，BARRIERE T，GELIN J C，et al.Powder injection molding of metallic and ceramic hip implants［J］.International Journal of Powder Metallurgy，2009，45（3）：25-34.

［18］SANTOS T，VALENTE M A，MONTERIO J，et al.Electromagnetic and thermal history during microwave heating［J］.Applied Thermal Engineering，2011，31（16）：3255-3261.

［19］BOUVARD D，CHARMOND S，CARRY C P.Mutiphysics simulation of microwave sintering in a monomode cavity［C］//12th Seminar Computer Modeling in Microwave Engineering＆ Applications.Grenoble，France：［s.n］，2010.

［20］GERMAN R M.Sintering theory and practice［M］.［S.l.］：Wiley-Interscience Publication，1996.

［21］CHARMOND S，CARRY C P，BOUVARD D.Densification and microstructure evolution of Y-tetragonal zirconia polycrystal powder during direct and hybrid microwave sintering in a single-mode cavity［J］.Journal of the European Ceramic Society，2010，30（6）：1211-1221.