電磁場(chǎng)作用下金屬凝固過程的建模方法

殷子豪 張 帆 仲紅剛

（上海大學(xué)先進(jìn)凝固技術(shù)中心，上海 200444）

金屬的凝固是非常復(fù)雜的高溫動(dòng)態(tài)過程，金屬在自然狀態(tài)下的凝固組織往往由于溫度梯度而存在凝固的先后次序，從而實(shí)現(xiàn)了固/液界面的推進(jìn)。這一行為最終導(dǎo)致金屬件心部產(chǎn)生縮松縮孔以及成分不均勻等各種缺陷。

由于高效、非接觸的特點(diǎn)，施加電磁場(chǎng)改善金屬凝固組織已成為凝固領(lǐng)域研究的重點(diǎn)。早在1922年，Mcneil就提出了電磁攪拌的專利；20世紀(jì)50年代初，Junghaus和Schaaber等在德國(guó)Hückingen連鑄機(jī)上首次試驗(yàn)了電磁攪拌[1]，以改善鑄坯質(zhì)量。美國(guó)麻省理工學(xué)院Flemings等最早采用脈沖電流細(xì)化金屬凝固組織[2]，之后大量學(xué)者展開了脈沖電流[3-6]與脈沖磁場(chǎng)[7-13]對(duì)金屬凝固組織影響的一系列研究。上海大學(xué)翟啟杰團(tuán)隊(duì)將脈沖磁致振蕩凝固均質(zhì)化技術(shù)成功運(yùn)用于工業(yè)實(shí)踐[14-17]。

電磁場(chǎng)的引入使凝固過程更加復(fù)雜，需要對(duì)熔體流動(dòng)及組織演變行為進(jìn)行更深入的研究。但由于受實(shí)際環(huán)境的限制，一般很難直接通過試驗(yàn)得到結(jié)果。近幾年計(jì)算機(jī)模擬技術(shù)的飛速發(fā)展，使得數(shù)值模擬成為研究這一過程的重要手段之一。但數(shù)值模擬很難將所有的變量都考慮到模型中，因此合理的假設(shè)是建立數(shù)學(xué)模型的必要步驟。根據(jù)模型所求解的物理量，忽略一部分影響相對(duì)較小的變量，或者將其等效在其他物理量上，以便于方程的求解。通常，凝固過程數(shù)值模擬需要在建立幾何模型后對(duì)不同的幾何域定義不同的物理意義，并根據(jù)網(wǎng)格劃分進(jìn)行離散化處理，最終歸結(jié)為多元代數(shù)方程組，求解后獲得數(shù)值解。通過調(diào)整網(wǎng)格劃分與求解方式的設(shè)置等可以提高數(shù)值解的可靠性，進(jìn)一步通過實(shí)驗(yàn)室或者工業(yè)生產(chǎn)的數(shù)據(jù)進(jìn)行校正擬合與驗(yàn)證，才能得到相對(duì)可靠的數(shù)值模型，進(jìn)而從后續(xù)的參數(shù)化研究中得到接近真實(shí)的定量分析結(jié)論，以深刻理解金屬熔體凝固過程。

本文綜述了電磁場(chǎng)作用下金屬凝固過程的數(shù)值模擬方法及應(yīng)用，闡述了電磁場(chǎng)作用下流場(chǎng)與溫度場(chǎng)的耦合模型，分析比較了不同模型的應(yīng)用環(huán)境，并展望了未來的發(fā)展方向。

1 熔體內(nèi)電磁場(chǎng)數(shù)值模擬

1.1 矢量勢(shì)法

Fujisaki等[23]將電磁場(chǎng)視為準(zhǔn)穩(wěn)態(tài)場(chǎng)，通過A-φ變量有限元法計(jì)算電磁場(chǎng)，對(duì)比了旋轉(zhuǎn)攪拌和平行攪拌的電磁效應(yīng)，發(fā)現(xiàn)旋轉(zhuǎn)磁場(chǎng)下的電磁力存在渦流分布特征，這有利于提高鑄坯的表面質(zhì)量及優(yōu)化夾雜物。Kolesnichenko等[24]通過邊界元法計(jì)算了圓形立式連鑄機(jī)二冷區(qū)的二維電磁場(chǎng)分布，通過邊界元法求解空氣中的磁場(chǎng)分布，獲得的表達(dá)式作為求解金屬熔體內(nèi)部電磁場(chǎng)的邊界條件，再利用A-φ法的有限差分方程求解熔體內(nèi)部電磁場(chǎng)分布。

Zhang等[10]基于ANSYS中的磁場(chǎng)模塊通過磁矢量勢(shì)法得到空間電磁場(chǎng)分布。在模擬引入諧波磁場(chǎng)時(shí)，將體積力分為時(shí)間相關(guān)分量和時(shí)間無關(guān)分量。通過時(shí)均體積力，解決電磁場(chǎng)周期與熔體動(dòng)量響應(yīng)時(shí)間尺度相差較大的問題。對(duì)電磁場(chǎng)和流場(chǎng)采用不同大小的網(wǎng)格，利用線性插值法將電磁力精確插入不同的網(wǎng)格中：

式中：FFLUENT為插值點(diǎn)的時(shí)均電磁力密度，該點(diǎn)為FLUENT單元格的質(zhì)點(diǎn)；FiANSYS為ANSYS中第i個(gè)點(diǎn)的時(shí)均電磁力密度；di為插值點(diǎn)和第i個(gè)點(diǎn)的距離，如圖1所示。

圖1 ANSYS（a）和FLUENT 網(wǎng)格（b）及插值流程示意圖（——ANSYS；···FLUENT）（c）Fig.1 Mesh in ANSYS（a）and FLUENT（b）and schematic diagram of interpolation procedure（——ANSYS；···FLUENT）（c）

1.2 直接求解法

通過直接法[25]也可求解空間磁場(chǎng)分布，即對(duì)磁場(chǎng)強(qiáng)度H直接求解，根據(jù)Biot-Savart定律，有：

式中：V為材料體積的積分；r為離源點(diǎn)的距離；r′為積分用虛擬變量。

Meyer等[25]利用直接求解法對(duì)連鑄方坯和板坯的電磁攪拌過程進(jìn)行數(shù)值模擬，發(fā)現(xiàn)三維模型的模擬結(jié)果與二維模型的不同，即線圈上下邊緣處鑄坯橫截面上電磁力明顯衰減，縱截面上電磁力以軸向分量為主。

R?biger等[26]通過OPERA 3D 程序模擬了脈沖電極從頂部導(dǎo)入GaInSn共晶合金中的電磁場(chǎng)與流場(chǎng)分布規(guī)律。針對(duì)電荷守恒和不同電導(dǎo)率材料界面處的電流密度J連續(xù)性，即Jn＝－σ1n×▽?duì)?＝ －σ2n×▽?duì)?，求解電勢(shì)的Laplace方程▽2φ＝0。其中感應(yīng)磁場(chǎng)利用Biot-Savart定律通過電流密度J＝σ（－▽?duì)眨┯?jì)算，即：

電勢(shì)的邊界條件通過電極末端的電流幅值定義：IDC＝∫Ads×J，其中A為電極橫截面積。模擬發(fā)現(xiàn)感應(yīng)磁場(chǎng)主要集中在電極附近，且最大值出現(xiàn)在電極表面。感應(yīng)焦耳熱主要分布在電極底部，電極底部的電流密度遠(yuǎn)大于其他區(qū)域。電磁力同樣集中于電極底部，方向主要向下。Xu等[27]將電極從側(cè)面導(dǎo)入工業(yè)純鋁中，脈沖電流的表達(dá)式可表示為i（t）＝ I0e－αtsin（2πft），其中I0為電流幅值，α為衰減系數(shù)，f為電流頻率。為了研究電極位置對(duì)組織的影響，分別模擬了電極兩端位于一側(cè)和兩側(cè)的情況。

Zhang等[28]模擬了脈沖電流作用下定向凝固的鋁硅合金熔體流動(dòng)情況，熔體中洛倫茲力分布如圖2所示。硅顆粒非對(duì)稱時(shí)，洛倫茲力的變化主要出現(xiàn)在糊狀區(qū)，且初生硅附近的糊狀區(qū)受到的洛倫茲力減少，其余區(qū)域均增大；V型界面上方，洛倫茲力明顯增強(qiáng)，洛倫茲力的總和將熔體從熔池邊緣推向中心。

圖2 脈沖電流作用下熔體中洛倫茲力分布Fig.2 Distributions of Lorentz force in melt under electric current pulse

式中：Ik和Φk分別為第k匝線圈多相電流的幅值與相移；dlk為線圈的單元長(zhǎng)度；dV為熔體的單元體積；｜r－r′｜為熔體中某一點(diǎn)到各自積分單元的距離。

式（16）等號(hào)右邊第一項(xiàng)表示磁場(chǎng)的旋轉(zhuǎn)部分，第二和第三項(xiàng)則分別表示線圈與感應(yīng)電流對(duì)磁場(chǎng)的影響，因此可以涵蓋任意形式感應(yīng)線圈的電磁場(chǎng)計(jì)算，并且不涉及自由空間的電磁場(chǎng)計(jì)算，大大節(jié)省了計(jì)算機(jī)的運(yùn)算時(shí)間與存儲(chǔ)空間。這是首次嘗試在三維坐標(biāo)系下進(jìn)行方坯和板坯內(nèi)的電磁場(chǎng)與流場(chǎng)耦合計(jì)算，也為今后對(duì)方坯和板坯內(nèi)電磁場(chǎng)、流場(chǎng)和溫度場(chǎng)的復(fù)雜建模提供了理論框架。

2 電磁場(chǎng)作用下流場(chǎng)與溫度場(chǎng)耦合模型

由于電磁場(chǎng)、流場(chǎng)與溫度場(chǎng)的時(shí)間尺度問題，往往先對(duì)電磁場(chǎng)單獨(dú)求解，再對(duì)流場(chǎng)和溫度場(chǎng)進(jìn)行耦合求解。目前常用的流場(chǎng)數(shù)值法為雷諾平均納維-斯托克斯模型（Reynolds-averaged Navier-Stokes equations，RANS）[30]，由于RANS 采用了經(jīng)驗(yàn)壁面定律，可在粗大邊界層上捕捉大速度梯度，因此這種方法具有很高的計(jì)算效率，且可精確估計(jì)穩(wěn)態(tài)流型[31]。

描述熔體流動(dòng)的基本方程如下：

不可壓縮流體的質(zhì)量守恒方程：

不可壓縮流體的動(dòng)量守恒方程：

式中：Sm為電磁力、熱浮力、流動(dòng)阻尼力等體積源項(xiàng)。

熔體流動(dòng)過程中往往伴隨著熱交換，使得熔體溫度不斷變化。因此研究流場(chǎng)的同時(shí)，必須考慮熔體的傳熱。根據(jù)傅里葉傳熱與能量守恒定律，可得到流體傳熱的能量守恒方程。

描述傳熱過程的基本方程：

式中：Q為熱源項(xiàng)。

金屬凝固傳熱通常分為單個(gè)域與多個(gè)域兩種方法。焓-多孔介質(zhì)模型[32-33]是一種能較好地模擬流動(dòng)與凝固傳熱的單域模型，其考慮了糊狀區(qū)內(nèi)潛熱釋放與流動(dòng)，收斂性較快，并且與試驗(yàn)結(jié)果[34]吻合較好。這與Bennon 的連續(xù)模型[35]和Beckermann的體積平均模型[36]內(nèi)容基本一致，為后續(xù)電磁場(chǎng)作用下熔體內(nèi)流場(chǎng)與溫度場(chǎng)的耦合模擬奠定了基礎(chǔ)。多域模型[37]由于復(fù)雜的界面域網(wǎng)格重構(gòu)而使用較少。

金屬凝固過程中，糊狀區(qū)同時(shí)包含柱狀晶、等軸晶和液相，而柱狀晶和等軸晶會(huì)阻礙金屬液在糊狀區(qū)的流動(dòng)，進(jìn)而影響金屬的傳熱與凝固過程。一般采用變黏度法[38-40]和達(dá)西源項(xiàng)法[41]處理糊狀區(qū)的流動(dòng)。

2.1 變黏度法

Roplekar等[42]基于體積平均模型，結(jié)合潛熱回升-混合長(zhǎng)度模型研究了旋轉(zhuǎn)磁場(chǎng)下鋁合金半連鑄過程中的宏觀偏析現(xiàn)象，并進(jìn)行了試驗(yàn)驗(yàn)證。將式（18）中的熔體黏度μ寫成渦流黏度μt與分子黏度μ0之和，其表達(dá)式分別為：

式中：lm表示湍流中液滴混合或碰撞的平均自由程，且lm＝ min（0.435ln，0.09lc），其中l(wèi)n為距最近壁面的距離，lc為特征長(zhǎng)度，本文取圓柱結(jié)晶器的半徑；fs為固相率。μ0曲線如圖3（a）所示。

Kumar等[38]在模擬電磁攪拌條件下鋁合金半連鑄凝固過程的相變時(shí)，采用Flemings[43]測(cè)定的黏度數(shù)據(jù)（圖3（b））代入方程。研究發(fā)現(xiàn)，半固態(tài)熔體在結(jié)晶器內(nèi)受電磁力影響，壁面處熔體沿外壁向上運(yùn)動(dòng)，中心處熔體沿中軸線向下運(yùn)動(dòng)。由于電磁攪拌的作用，結(jié)晶器內(nèi)溫度分布更加均勻。

圖3 黏度隨固相率的變化Fig.3 Variation of viscosity with solid fraction

2.2 達(dá)西源項(xiàng)法

然而，黏度變化曲線往往很難獲取，因此達(dá)西源項(xiàng)法的應(yīng)用更為廣泛。但學(xué)者對(duì)達(dá)西源項(xiàng)的表達(dá)式存在不同的意見[41，44-46]。目前常用的有一次枝晶臂間距法、二次枝晶臂間距法和經(jīng)驗(yàn)值法。其中，大多數(shù)學(xué)者[36，47-49]依據(jù)Carman-Kozeny 方程通過枝晶臂間距確定多孔介質(zhì)中的滲透率K，滲透率K與達(dá)西源項(xiàng)SD的關(guān)系為：

（1）一次枝晶臂間距法

Pfeiler等[50]根據(jù)Blake-Kozney 方程將達(dá)西源項(xiàng)寫成如式（23）所示的形式，其中λ1為一次枝晶臂間距。在此基礎(chǔ)上預(yù)測(cè)了結(jié)晶器內(nèi)鋼液的流動(dòng)狀況、凝固坯殼的生長(zhǎng)和夾雜物在固/液界面的運(yùn)動(dòng)軌跡，并證明了直徑大于一次枝晶臂間距的夾雜物難以被糊狀區(qū)捕獲。

Ramirez等[47]分別對(duì)比了滲透率K兩種不同計(jì)算方法的適用條件，如式（24）、（25）所示。而Gu等[51]則通過式（24）模擬研究了大型鑄錠凝固過程中的流動(dòng)、溫度及物質(zhì)濃度分布規(guī)律：

Aboutalebi等[52]利用式（26）模擬了拉速、鋼種、水口結(jié)構(gòu)對(duì)流型和凝固形態(tài)的影響，并與實(shí)測(cè)結(jié)果吻合較好。Willers等[49]根據(jù)該方法研究了旋轉(zhuǎn)磁場(chǎng)下Al-Si合金在模鑄條件下的凝固過程，發(fā)現(xiàn)提高變向頻率能消除低頻時(shí)熔體中磁場(chǎng)引起的溫度波動(dòng)現(xiàn)象。

（2）二次枝晶臂間距法

Zhang等[10]基于體積平均模型將焦耳熱與凝固潛熱加入能量守恒方程，利用達(dá)西源項(xiàng)處理糊狀區(qū)流動(dòng)問題。設(shè)為枝晶干涉臨界固相率，令初始滲透率K0＝（dSDAS）2/180（dSDAS為二次枝晶臂間距）。表1為凝固不同階段的黏度與滲透率。

表1 不同階段的黏度與滲透率Table 1 Viscosity and permeability in the different stages

表1中，χ是使能量守恒方程收斂的較小正數(shù)；X為需要求解的變量；Xs表示固相中的對(duì)應(yīng)變量。研究發(fā)現(xiàn)[53]，隨著電磁場(chǎng)頻率的增加，熔體流動(dòng)強(qiáng)度先增大后減小，而溫度場(chǎng)和熔池深度變化不大。

Jiang等[54]在多孔介質(zhì)模型的基礎(chǔ)上結(jié)合表觀黏度模型，建立了電磁攪拌作用下連鑄二冷區(qū)凝固行為與溶質(zhì)傳輸模型；并對(duì)柱狀晶區(qū)采用多孔介質(zhì)模型，將等軸晶區(qū)視為漿體區(qū)，采用表觀黏度模型模擬熔體流場(chǎng)。糊狀區(qū)滲透系數(shù)K如式（27）所示。對(duì)等軸晶區(qū)內(nèi)的鋼液流動(dòng)情況，采用μm替代式（18）中的μ，如式（28）所示：

Sun 等[55]為了在式（27）的基礎(chǔ)上增加一個(gè)較小的正數(shù)ξ，給出了二次枝晶臂與局部凝固時(shí)間的關(guān)系，如式（29）、（30）所示，并利用該方法模擬研究了M-EMS和F-EMS對(duì)高碳鋼連鑄方坯凝固過程與元素偏析的影響。通過連鑄坯表面溫度、EMS中心磁通量密度、鑄坯形貌及溶質(zhì)分布，驗(yàn)證了電磁場(chǎng)下連鑄坯的凝固與偏析過程。

（3）經(jīng)驗(yàn)值法

采用經(jīng)驗(yàn)值處理達(dá)西源項(xiàng)時(shí)引入一個(gè)糊狀區(qū)系數(shù)Amush：

國(guó)內(nèi)外學(xué)者[56-57]針對(duì)Amush對(duì)熔體流動(dòng)的影響進(jìn)行了探討。Hietanen等[56]模擬大方坯連鑄結(jié)晶器內(nèi)鋼液流動(dòng)時(shí)發(fā)現(xiàn)，較小的Amush會(huì)導(dǎo)致鋼液溫度下降較快，而較大的Amush會(huì)導(dǎo)致模型計(jì)算不易收斂。而李少翔等[57]模擬大方坯連鑄結(jié)晶器內(nèi)的流動(dòng)凝固行為時(shí)，發(fā)現(xiàn)隨著糊狀區(qū)系數(shù)的增大，糊狀區(qū)寬度變窄，并給出了合理的糊狀區(qū)范圍為（1～5）×108。

Zhang等[58]利用式（31）結(jié)合拉格朗日法求解了小方坯結(jié)晶器內(nèi)夾雜物的運(yùn)動(dòng)規(guī)律，該方法的計(jì)算結(jié)果與實(shí)測(cè)值基本吻合。Wang等[59]研究了板坯連鑄中流動(dòng)控制結(jié)晶器對(duì)坯殼生長(zhǎng)的影響，并通過對(duì)比固液相線與實(shí)測(cè)坯殼厚度的變化曲線，驗(yàn)證了模型的可靠性。Sun等[60]研究了不同水口類型與M-EMS組合對(duì)結(jié)晶器內(nèi)鋼液流動(dòng)及坯殼厚度的影響，其中Amush＝105。He等令A(yù)mush＝108，建立了3種類型的浸入式水口對(duì)結(jié)晶器內(nèi)鋼液流動(dòng)的影響，并通過物理水模型驗(yàn)證了流場(chǎng)的可靠性。近幾年將Amush代入達(dá)西源項(xiàng)模擬研究熔體的凝固過程[44，61-65]已成為一種趨勢(shì)。

除了在連續(xù)鑄鋼中的模擬應(yīng)用外，東北大學(xué)樂啟熾團(tuán)隊(duì)利用該方法建立了脈沖磁場(chǎng)下鎂合金半連續(xù)鑄造的連續(xù)介質(zhì)模型，其中Amush＝105。研究發(fā)現(xiàn)，差相脈沖磁場(chǎng)作用下的z向電磁力與z向流速均大于一般脈沖磁場(chǎng)下的，從而導(dǎo)致了更強(qiáng)的對(duì)流、更加均勻的溫度分布以及更加細(xì)化的組織[66]；同時(shí)也分析了電磁力、電流峰值、電壓占空比及頻率對(duì)熔體的影響，發(fā)現(xiàn)電流增大，熔體流動(dòng)性增強(qiáng)，導(dǎo)致熔體熱量散失，熔體溫度和熔池深度均降低；頻率的增加則不會(huì)引起多大變化；占空比的增加僅對(duì)電磁力有所增益[67]。

3 總結(jié)與展望

本文介紹了電磁場(chǎng)數(shù)學(xué)模型、流場(chǎng)與溫度場(chǎng)耦合模型在電磁場(chǎng)處理金屬凝固過程中的建模與應(yīng)用?？偨Y(jié)了求解電磁場(chǎng)空間分布的部分算法模型，其中矢量勢(shì)求解法主要應(yīng)用于時(shí)變感應(yīng)電磁場(chǎng)，可直接用于解決脈沖、諧波等瞬變電磁參數(shù)問題；直接求解法一般采用Biot-Savart定律對(duì)電流產(chǎn)生的磁場(chǎng)求解，適用于解決直流電等恒定電磁參數(shù)的問題，或者先將瞬變信號(hào)轉(zhuǎn)化為等效的恒定信號(hào)再求解計(jì)算，具有一定的局限性。接著介紹了電磁場(chǎng)作用下流場(chǎng)與溫度場(chǎng)的耦合模型，以基于單域模型的體積平均模型為主，闡述了不同的糊狀區(qū)處理方法，其中添加達(dá)西源項(xiàng)的方法應(yīng)用較為廣泛，但對(duì)達(dá)西源項(xiàng)系數(shù)的取值仍存在不同意見；此外，將糊狀區(qū)黏度處理為不同形式的函數(shù)也是解決方法之一。

目前已能夠相對(duì)準(zhǔn)確地預(yù)測(cè)電磁場(chǎng)作用下的熔體演變過程、溫度分布、坯殼生長(zhǎng)以及金屬液受電磁力驅(qū)動(dòng)的強(qiáng)制對(duì)流流場(chǎng)分布等，但大多仍是基于試驗(yàn)數(shù)據(jù)而展開的。將更多的物理場(chǎng)耦合在一起，建立不同空間尺度與時(shí)間尺度的數(shù)值模擬，并應(yīng)用于實(shí)際工業(yè)環(huán)境，仍然是個(gè)巨大的挑戰(zhàn)。