三維電磁場有限元數據到有限差格式的轉換算法

龔海軍，徐達鳴，傅恒志

(哈爾濱工業大學材料科學與工程學院，哈爾濱150001，331ghj@163.com)

電磁約束凝固成形和電磁連鑄等電磁材料加工及其數值模擬技術是目前研究較為活躍的領域［1－3］.電磁場作用下的合金凝固是一個多場作用下的熱量、質量及流體動量強耦合傳輸過程［4］，如冷坩堝電磁定向凝固過程中，高溫金屬熔體與冷坩堝壁之間存在劇烈的傳熱，同時電磁場對冷坩堝中的金屬熔體有高密度的加熱、攪拌及約束作用［5－8］.對這些發生在四維時空中的復雜多物理場耦合問題，通常需要采用有限元(FEM)、有限差(FDM)或有限體積(FVM)法等進行數值解析.許多多物理場相互作用下的復雜場量耦合計算常聯合使用FEM和FDM/FVM來進行，并已成為一種趨勢［9－12］.由于FEM和FDM/ FVM剖分網格空間結構差異較大，高效率的多場FEM/FDM耦合計算須有效解決兩種數據格式的匹配，特別是三維FEM→FDM的數據格式轉換問題.

為有效地將高效率凝固耦合傳輸統一模型［13－15］應用于電磁凝固傳輸計算中，需采用FEM與FVM相結合的耦合計算方法，即應用通用有限元軟件ANSYS進行電磁場量的計算，采用基于有限差法的凝固傳輸統一數學模型計算電磁凝固耦合傳輸過程.為此，本文提出一種將電磁場FEM計算結果轉換成凝固傳輸耦合計算所需的FDM/FVM數據格式的有效方法.

1 三維電磁場計算模型

為檢驗本文數據轉換方法和程序的普適性，首先對冷坩堝電磁定向凝固系統、電磁連鑄系統、電磁熔煉和攪拌系統中的電磁場進行計算，然后將有限元格式的數據轉換為限差格式并顯示.其中，冷坩堝電磁定向凝固系統造型和剖分如圖1所示，模型比例為1∶1.本文計算模型采用ANSYS基于節點的方法，諧波模型除遠場空氣采用INFIN111單元外，其余均采用SOLID97單元，空氣單元定義磁矢位自由度，導體單元定義磁矢位和電壓自由度.

圖1 冷坩堝電磁定向凝固計算模型結構及FEM網格剖分

2 三維有限元計算結果的有限差轉換

2.1 ANSYS有限元計算結果的存儲形式

ANSYS節點法求解電磁場得到的是節點解，對于靜磁場計算結果的插值轉換，需要輸出結點坐標信息(*c.lis)、單元 －材質 －節點信息(*mn.lis)、節點形式的磁感應強度結果(*b.lis)、有單元信息的磁感應強度結果(*e.lis)4個文件;對諧波和行波磁場，則還需另外輸出焦耳熱(*j.lis)和電流密度(*d.lis)兩個結果文件.實際電磁計算中由于模型形狀復雜，無法全部用六面體單元剖分網格，在某些區域會以退化的單元形式出現.通常情況下，8節點空間等參數單元合并其中1個或幾個節點，便可以退化為4～7節點的單元.在ANSYS模型中，六面體形單元只以楔形、金字塔形和四面體形3種退化單元形狀來協調網格劃分，如圖2所示，不會出現7節點六面體單元(即圖2(a)中六面體合并任意相鄰兩點).故一般情況下，電磁場的計算值存儲于4種不同形狀和不同有效節點數目的單元或單元節點上.

圖2 ANSYS三維電磁場有限元計算所用的4種單元形狀

2.2 有限元計算結果的有限差轉換

為將3D－EM的FEM計算結果應用于電磁凝固耦合傳輸計算，需將FEM數據向FDM形式轉換.對于如圖2(a)所示的三維一次8節點六面體單元，設:A為其內部任意一點，8個節點按右手定則順序依次對其編號為i，j，k，l，m，n，p，q.根據有限元線性插值理論，電磁場量節點值在單元內沿3個坐標軸(x，y，z)方向上將呈線性變化，x方向分量u=u(x，y，z)可寫成

式中:a1，a2，…，a8分別為待定系數.六面體單元內任意一點A在x方向上的電磁場分量u可寫成

其中:

式中:［N］為形函數，是坐標x，y，z的函數，{φ}為六面體單元8個節點在x方向的電磁場分量構成的列向量，由 ANSYS計算并導出.|Δi|，…，|Δq|分別為將|Δ|內的坐標值xi，yi，zi，…，xq，yq，zq替換為xA，yA，zA.同理，y，z方向的電磁場分量v，w的插值計算可類似地寫出.

對于在如圖3(c)所示的四面體一次單元，同樣設A為四面體內部任意一點，4個節點按右手定則順序依次對其編號為i，j，k，l.則A點處某一電磁場量值為

其中:

式中:V為由i，j，k，l 4點構成的四面體體積，Vi，Vj，Vk，Vl分別為四面體內的A點取代i，j，k，l構成的四面體體積.

其它兩種形狀單元的形函數可類似導出.顯而易見，采用四面體及其形函數進行線性插值計算相對其它單元形狀要簡單和方便，故本文將六面體形、楔形和金字塔形等非四面體單元都分解為便于插值計算的四面體單元來處理.

根據空間幾何關系，任意六面體形單元可劃分為5個或6個四面體，且分別有兩種劃分方式，以A5、B5和A6、B6表示，如圖3(a)、圖3(b)所示.A5型劃分為6 873，6 581，6 831，6 321，8 314，B5型劃分為7 652，8 754，2 574，7 324，1 245;A6型劃分7 652，2 573，5 321，8 753，8 513，1 483，B6型劃分7 652，8 752，5 218，7 832，1 238，1 348.為減少計算工作量，本文將六面體單元劃分為5個四面體，楔形和金字塔形也做類似的分解處理，所有分解出的四面體子單元及其節點將攜帶原母單元及其節點上的所有信息.

圖3 六面體FEM單元分解為四面的兩種方式及FDM中心點與分解的四面體FEM的對應關系

確定了非四面體單元分解次序，便可通過體積判斷來對FDM中心點與FEM單元進行查找和對應.具體方法是:如果FDM中心點A落入某非四面體FEM單元分解出的一個小四面體內(圖3(c))，則A點與此小四面體4個頂點組成的4個新小四面體的體積和與原小四面體體積將相等(＜ε，ε為一誤差限)，那么便認為該FDM中心點落入此四面體FEM單元中，隨后以其對應的四面體的4個節點來線性插值;如果FDM中心點A落入小四面體外(圖3(d))，那么根據四面體體積式(6)，在頂點次序為非右手定則順序的情況下，A點與此小四面體4個頂點組成的4個新小四面體的體積和將為負值或一小于原小四面體體積的值.ANSYS網格原本為四面體形的單元無需再分解，可直接進行體積判斷來與FDM單元查找對應，確認對應關系后便可對場量進行插值計算.根據線性插值理論，上述插值計算方法對FDM中心點落在FEM單元表面、邊界或單元節點上的情況同樣適用.

需指出的是，ANSYS計算的焦耳熱和電流密度結果是按單元輸出的，不能按上述方式直接插值計算，此時FDM中心點與FEM單元對應獲取焦耳熱和電流密度有兩種方案:1)當某FDM中心代表點落入某單元時，就以該單元的焦耳熱和電流密度值作為FDM中心點的相應值，其缺點是當多個FDM中心代表點落入同一個FEM單元時取值都相同，不能很好顯示焦耳熱和電流密度值沿空間位置連續變化的趨勢，且在靠近表面處焦耳熱和電流密度值的插值結果將偏低;2)將FEM單元的焦耳熱和電流密度值向單元的節點進行“節點化”處理.即對于所有節點，先統計同一節點屬于多少相鄰單元共用，然后將此節點上的焦耳熱和電流密度值取共用該節點的單元焦耳熱和電流密度值的算術平均，最后再按上述處理磁感應強度的方法來查找對應和插值計算.即本文選用第2套方案，這樣處理的結果，也不可避免地會導致緊鄰樣件表面處焦耳熱和電流密度值的插值結果偏低，但可以通過在場量梯度較大區域加細網格，保證足夠的精度.

基于上述有限元的線性插值原理和所提出的實施方案，本文采用Fortran95語言編寫了三維FEM→FDM數據轉換程序，具體流程如圖4所示.

圖4 FEM→FDM數據轉換計算程序流程

對于一些幾何形狀復雜的模型，經過一次完整循環，有時會出現某些有限差中心點找不到與之對應有限元單元的情況，也就是說沒有得到插值結果，這時就需適當調整誤差限ε，進行二次查找對應和插值，一般經過二次查尋都能找到對應單元;幾何形狀規則的模型，不會或很少出現需要二次查尋和對應的情況.通常情況下，為了提高耦合計算精度和適應模型中變尺寸區域，FDM網格比有限元單元小而多，FEM/FDM轉換計算的單元對應過程中可能有多個FDM中心點落入同一FEM單元，如按六面體和五面體形函數進行插值計算，則這些落入同一FEM單元的多個FDM中心點的插值結果相近，這樣便不能很好反映出因坐標變化而帶來的場量變化;如將非四面體FEM單元分解，當 FDM中心點落入某一非四面體FEM單元分解出的小四面體后，以此中心點所在的小四面體節點值進行插值計算，而當下一個FDM中心點再次落入同一非四面體FEM單元時，則以其所對應單元的另一分解出的小四面體節點值進行插值計算.也就是說，不同時以非四面體單元的所有節點來插值，這樣便增加了插值結果對坐標的敏感性，這是分解非四面體單元進行插值計算的另一優點.本文計算和轉換工作在1臺CPU主頻為3.0 GHz、內存為3.0 GB的PC上運行，各計算模型數據及轉換計算用時對比如表1所示.

表1 各計算模型FEM→FDM數據轉換計算結果對比

3 三維矢量場后處理顯示對比

目前有許多商用軟件自備后處理模塊，可顯示3D標量和矢量場，但是這些軟件的顯示模塊通常不通用.在凝固傳輸過程的研究中，尤其是對于涉及凝固的流動過程，純液相區中的流速與固液兩相區中的速度相比一般要大1～2個數量級以上，而這兩個區域的流動行為對于凝固缺陷如宏觀偏析等的定量預測極其重要，需要區分純液相區和固液兩相區采用雙速度標尺才能同時有效顯示整個凝固鑄錠/鑄件中的流場分布.為此，本文通過Visual Fortran 6.6A的QuickWin應用程序平臺，用Fortran 95語言自行編寫了集FEM/FDM數據轉換計算和鑄件凝固傳輸數值模擬數據3D圖形顯示功能于一體的后處理程序.

圖5為冷坩堝電磁定向凝固系統內TiAl合金錠的各電磁場量對比，其中，圖5(a)、圖5(c)、圖5(e)是ANSYS計算顯示的結果，圖5(b)、圖5(d)、圖5(f)是本文程序將ANSYS計算結果轉為有限差后的顯示結果.由圖5(a)可見，磁感應強度在線圈所纏繞的中部區域最強，矢量方向沿合金鑄錠向下，圖5(b)是數據換后的等軸測顯示結果，可見，轉換后磁感應強度的大小和方向與轉換前是一致的.圖5(c)、圖5(e)所示分別為ANSYS計算輸出的上下錠中焦耳熱和電流密度分布，相應的轉換結果如圖5(d)、圖5(f)所示，可見轉換后的電磁場量分布與轉換前吻合.

圖5 上、下鈦合金鑄錠中電磁場的ANSYS模擬結果與本文FEM→FDM數據轉換后結果的對比

本文編寫的后處理程序可根據實際需要，選擇等軸測或斜二測顯示二維和三維場量圖像，如果必要，3個坐標軸方向還可以任意互換，使圖像在屏幕上以不同視角顯現.圖6(b)即為連鑄熔池內靜磁場的斜二測分層顯示，其分布規律和數值大小與轉換前(圖6(a))一致，且水嘴出口處的液流和磁場值都顯示正常.圖6(c)、圖6(d)所示為某一碳化硅坩堝內鋁合金熔體的磁感應強度［16］轉換前后對比，結果表明:對于雙變曲率的曲面結構，本文提出與開發的FEM→FDM轉換算法和計算、顯示程序同樣適用.

4 結論

1)將六面體有限元剖分單元及其退化單元統一地分解為簡單四面體單元，然后進行FEM→FDM三維空間對應和數據插值計算這一處理方法是簡便、可行的.

2)對三維電磁場有限元計算實例的數據處理及結果顯示表明，本文提出的FEM→FDM數據轉換計算方法及開發的顯示手段是成功和有效的.

圖6 不同電磁計算模型中磁感應強度的ANSYS－FEM模擬結果與本文FEM→FDM數據轉換后結果的對比

3)本文FEM→FDM數據轉換算法和程序適用于由有限元法計算的任意三維矢/標量場數據向有限差格式轉換計算，這一算法及程序與自行開發的3D圖形顯示軟件為實現任意三維電磁凝固傳輸FEM－FDM耦合數值計算及工藝優化提供了必要的技術手段.

［1］ASAI S.Recent development and prospect of electromagnetic processing［J］.Science and Technology of Advanced Materials，2000，1(4):191－200.

［2］FORT J，KLYMYSHYN N，GARNICH M.Electromagnetic and thermal－flow modeling of a cold－wall crucible induction melter［J］.Metallurgical and Materials Transactions B，2005，36(1):141－152.

［3］DING Hongsheng，CHEN Ruirun，GUO Jingjie，et al.Directional solidification of titanium alloys by electromagnetic confinement in cold crucible［J］.Materials Letters，2005，59(7):741－745.

［4］包燕平，張濤，蔣偉，等.板坯連鑄結晶器內鋼液流場的三維數學模型［J］.北京科技大學學報，2001，23(2):106－110.

［5］李嬌，文光華，祝明妹，等.寬板坯結晶器流場和溫度場數值模擬優化分析［J］.重慶大學學報，2009，32(2):173－176.

［6］陳瑞潤，丁宏升，郭景杰，等.冷坩堝連續熔鑄與定向凝固Ti6Al4V合金的溫度場計算［J］.稀有金屬材料與工程，2007，36(10):1722－1727.

［7］TANAKA T，YASHIRO N，SHIRAI Y，et al.LPE growth of AlN single crystal using cold crucible under atmospheric nitrogen gas pressure［J］.Physical Status Solid C－Current Topics in Solid State Physics，2007，4(7):2227－2230.

［8］WU Shiping，LIU Dongrong，SU Yanqing，et al.Modeling of microstructure formation of Ti－6Al－4V alloy in a cold crucible under electromagnetic field［J］.Journal of Alloys and Compounds，2008，456(1/2):85－95.

［9］徐艷，康進武，黃天佑.鑄造過程溫度場/應力場雙向耦合的數值模擬［J］.清華大學學報(自然科學版)，2008，48(5):769－772.

［10］李輝，李志強.鑄造過程中的應力場數值模擬［J］.中國鑄造裝備與技術，2007(6):13－15.

［11］李培鋒，孫立斌，康進武，等.基于微機的鑄件凝固過程應力數值模擬及工程應用［J］.鑄造技術，2001 (3):5－8.

［12］BAI Yunfeng，XU Daming，MAO Lihe，et al.FEM/FDM－joint simulation for transport phenomena in directionally solidifying TiAl casting under electromagnetic field［J］.ISIJ International，2004，44(7):1173－1179.

［13］XU Daming，BAI Yunfeng，GUO Jingjie，et al.Numerical simulation of heat，mass and momentum transport behaviors in directionally solidifying alloy castings under electromagnetic fields using an extended direct－simple scheme［J］.Int J Numerical Methods in Fluids，2004，46(7):767－791.

［14］XU Daming，NI Jun.A numerical approach of direct－SIMPLE deduced pressure equations to simulations of transport phenomena during shaped casting［C］//Proceedings of the First International Multi－Symposium on Computer and Computational Sciences.Washington: IEEE Computer Society，2006:815－821.

［15］XU Daming，LI Qingchun.Numerical method for solution of strongly coupled binary alloy solidification problems［J］.Numerical Heat Transfer，Part A，1991，20(2):181－201.

［16］VIVèS Ch，RICOU R.Fluid flow phenomena in a single phase coreless induction furnace［J］.Metallurgical Transactions B，1985，16(2):227－235.