三維電磁場(chǎng)有限元數(shù)據(jù)到有限差格式的轉(zhuǎn)換算法

龔海軍，徐達(dá)鳴，傅恒志

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

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

為有效地將高效率凝固耦合傳輸統(tǒng)一模型［13－15］應(yīng)用于電磁凝固傳輸計(jì)算中，需采用FEM與FVM相結(jié)合的耦合計(jì)算方法，即應(yīng)用通用有限元軟件ANSYS進(jìn)行電磁場(chǎng)量的計(jì)算，采用基于有限差法的凝固傳輸統(tǒng)一數(shù)學(xué)模型計(jì)算電磁凝固耦合傳輸過(guò)程.為此，本文提出一種將電磁場(chǎng)FEM計(jì)算結(jié)果轉(zhuǎn)換成凝固傳輸耦合計(jì)算所需的FDM/FVM數(shù)據(jù)格式的有效方法.

1 三維電磁場(chǎng)計(jì)算模型

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

圖1 冷坩堝電磁定向凝固計(jì)算模型結(jié)構(gòu)及FEM網(wǎng)格剖分

2 三維有限元計(jì)算結(jié)果的有限差轉(zhuǎn)換

2.1 ANSYS有限元計(jì)算結(jié)果的存儲(chǔ)形式

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

圖2 ANSYS三維電磁場(chǎng)有限元計(jì)算所用的4種單元形狀

2.2 有限元計(jì)算結(jié)果的有限差轉(zhuǎn)換

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

式中:a1，a2，…，a8分別為待定系數(shù).六面體單元內(nèi)任意一點(diǎn)A在x方向上的電磁場(chǎng)分量u可寫(xiě)成

其中:

式中:［N］為形函數(shù)，是坐標(biāo)x，y，z的函數(shù)，{φ}為六面體單元8個(gè)節(jié)點(diǎn)在x方向的電磁場(chǎng)分量構(gòu)成的列向量，由 ANSYS計(jì)算并導(dǎo)出.|Δi|，…，|Δq|分別為將|Δ|內(nèi)的坐標(biāo)值xi，yi，zi，…，xq，yq，zq替換為xA，yA，zA.同理，y，z方向的電磁場(chǎng)分量v，w的插值計(jì)算可類似地寫(xiě)出.

對(duì)于在如圖3(c)所示的四面體一次單元，同樣設(shè)A為四面體內(nèi)部任意一點(diǎn)，4個(gè)節(jié)點(diǎn)按右手定則順序依次對(duì)其編號(hào)為i，j，k，l.則A點(diǎn)處某一電磁場(chǎng)量值為

其中:

式中:V為由i，j，k，l 4點(diǎn)構(gòu)成的四面體體積，Vi，Vj，Vk，Vl分別為四面體內(nèi)的A點(diǎn)取代i，j，k，l構(gòu)成的四面體體積.

其它兩種形狀單元的形函數(shù)可類似導(dǎo)出.顯而易見(jiàn)，采用四面體及其形函數(shù)進(jìn)行線性插值計(jì)算相對(duì)其它單元形狀要簡(jiǎn)單和方便，故本文將六面體形、楔形和金字塔形等非四面體單元都分解為便于插值計(jì)算的四面體單元來(lái)處理.

根據(jù)空間幾何關(guān)系，任意六面體形單元可劃分為5個(gè)或6個(gè)四面體，且分別有兩種劃分方式，以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.為減少計(jì)算工作量，本文將六面體單元?jiǎng)澐譃?個(gè)四面體，楔形和金字塔形也做類似的分解處理，所有分解出的四面體子單元及其節(jié)點(diǎn)將攜帶原母單元及其節(jié)點(diǎn)上的所有信息.

圖3 六面體FEM單元分解為四面的兩種方式及FDM中心點(diǎn)與分解的四面體FEM的對(duì)應(yīng)關(guān)系

確定了非四面體單元分解次序，便可通過(guò)體積判斷來(lái)對(duì)FDM中心點(diǎn)與FEM單元進(jìn)行查找和對(duì)應(yīng).具體方法是:如果FDM中心點(diǎn)A落入某非四面體FEM單元分解出的一個(gè)小四面體內(nèi)(圖3(c))，則A點(diǎn)與此小四面體4個(gè)頂點(diǎn)組成的4個(gè)新小四面體的體積和與原小四面體體積將相等(＜ε，ε為一誤差限)，那么便認(rèn)為該FDM中心點(diǎn)落入此四面體FEM單元中，隨后以其對(duì)應(yīng)的四面體的4個(gè)節(jié)點(diǎn)來(lái)線性插值;如果FDM中心點(diǎn)A落入小四面體外(圖3(d))，那么根據(jù)四面體體積式(6)，在頂點(diǎn)次序?yàn)榉怯沂侄▌t順序的情況下，A點(diǎn)與此小四面體4個(gè)頂點(diǎn)組成的4個(gè)新小四面體的體積和將為負(fù)值或一小于原小四面體體積的值.ANSYS網(wǎng)格原本為四面體形的單元無(wú)需再分解，可直接進(jìn)行體積判斷來(lái)與FDM單元查找對(duì)應(yīng)，確認(rèn)對(duì)應(yīng)關(guān)系后便可對(duì)場(chǎng)量進(jìn)行插值計(jì)算.根據(jù)線性插值理論，上述插值計(jì)算方法對(duì)FDM中心點(diǎn)落在FEM單元表面、邊界或單元節(jié)點(diǎn)上的情況同樣適用.

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

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

圖4 FEM→FDM數(shù)據(jù)轉(zhuǎn)換計(jì)算程序流程

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

表1 各計(jì)算模型FEM→FDM數(shù)據(jù)轉(zhuǎn)換計(jì)算結(jié)果對(duì)比

3 三維矢量場(chǎng)后處理顯示對(duì)比

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

圖5為冷坩堝電磁定向凝固系統(tǒng)內(nèi)TiAl合金錠的各電磁場(chǎng)量對(duì)比，其中，圖5(a)、圖5(c)、圖5(e)是ANSYS計(jì)算顯示的結(jié)果，圖5(b)、圖5(d)、圖5(f)是本文程序?qū)NSYS計(jì)算結(jié)果轉(zhuǎn)為有限差后的顯示結(jié)果.由圖5(a)可見(jiàn)，磁感應(yīng)強(qiáng)度在線圈所纏繞的中部區(qū)域最強(qiáng)，矢量方向沿合金鑄錠向下，圖5(b)是數(shù)據(jù)換后的等軸測(cè)顯示結(jié)果，可見(jiàn)，轉(zhuǎn)換后磁感應(yīng)強(qiáng)度的大小和方向與轉(zhuǎn)換前是一致的.圖5(c)、圖5(e)所示分別為ANSYS計(jì)算輸出的上下錠中焦耳熱和電流密度分布，相應(yīng)的轉(zhuǎn)換結(jié)果如圖5(d)、圖5(f)所示，可見(jiàn)轉(zhuǎn)換后的電磁場(chǎng)量分布與轉(zhuǎn)換前吻合.

圖5 上、下鈦合金鑄錠中電磁場(chǎng)的ANSYS模擬結(jié)果與本文FEM→FDM數(shù)據(jù)轉(zhuǎn)換后結(jié)果的對(duì)比

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

4 結(jié)論

1)將六面體有限元剖分單元及其退化單元統(tǒng)一地分解為簡(jiǎn)單四面體單元，然后進(jìn)行FEM→FDM三維空間對(duì)應(yīng)和數(shù)據(jù)插值計(jì)算這一處理方法是簡(jiǎn)便、可行的.

2)對(duì)三維電磁場(chǎng)有限元計(jì)算實(shí)例的數(shù)據(jù)處理及結(jié)果顯示表明，本文提出的FEM→FDM數(shù)據(jù)轉(zhuǎn)換計(jì)算方法及開(kāi)發(fā)的顯示手段是成功和有效的.

圖6 不同電磁計(jì)算模型中磁感應(yīng)強(qiáng)度的ANSYS－FEM模擬結(jié)果與本文FEM→FDM數(shù)據(jù)轉(zhuǎn)換后結(jié)果的對(duì)比

3)本文FEM→FDM數(shù)據(jù)轉(zhuǎn)換算法和程序適用于由有限元法計(jì)算的任意三維矢/標(biāo)量場(chǎng)數(shù)據(jù)向有限差格式轉(zhuǎn)換計(jì)算，這一算法及程序與自行開(kāi)發(fā)的3D圖形顯示軟件為實(shí)現(xiàn)任意三維電磁凝固傳輸FEM－FDM耦合數(shù)值計(jì)算及工藝優(yōu)化提供了必要的技術(shù)手段.

［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］包燕平，張濤，蔣偉，等.板坯連鑄結(jié)晶器內(nèi)鋼液流場(chǎng)的三維數(shù)學(xué)模型［J］.北京科技大學(xué)學(xué)報(bào)，2001，23(2):106－110.

［5］李嬌，文光華，祝明妹，等.寬板坯結(jié)晶器流場(chǎng)和溫度場(chǎng)數(shù)值模擬優(yōu)化分析［J］.重慶大學(xué)學(xué)報(bào)，2009，32(2):173－176.

［6］陳瑞潤(rùn)，丁宏升，郭景杰，等.冷坩堝連續(xù)熔鑄與定向凝固Ti6Al4V合金的溫度場(chǎng)計(jì)算［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ìn)武，黃天佑.鑄造過(guò)程溫度場(chǎng)/應(yīng)力場(chǎng)雙向耦合的數(shù)值模擬［J］.清華大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，48(5):769－772.

［10］李輝，李志強(qiáng).鑄造過(guò)程中的應(yīng)力場(chǎng)數(shù)值模擬［J］.中國(guó)鑄造裝備與技術(shù)，2007(6):13－15.

［11］李培鋒，孫立斌，康進(jìn)武，等.基于微機(jī)的鑄件凝固過(guò)程應(yīng)力數(shù)值模擬及工程應(yīng)用［J］.鑄造技術(shù)，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.