電渣重熔過程流場和溫度場的數值模擬

韓麗輝,　于春梅,　曲明磊

(北京科技大學 冶金與生態工程學院，北京 100083)

0　引　言

電渣重熔(Electro Slag Remelting，ESR)是利用電加熱渣來重熔自耗電極進行精煉鑄錠的一種方法，它能夠有效消除偏析、縮孔、疏松等缺陷。由于存在熔滴和熔渣的冶金反應，經電渣重熔的鋼具有純度高、含硫低、非金屬雜質少、鋼錠表面光潔、組織均勻致密、金相組織和化學成分均勻等優點，因此電渣重熔技術越來越多地應用到高溫合金、精密合金、模具鋼、航天軍用鋼等特殊鋼以及鋁、銅、銀等有色金屬合金生產領域[1,2]。在電渣重熔生產過程中，熔池形狀和深度影響著結晶的方向[3]；渣池溫度場及流場與自耗電極末端形狀、尺寸及熔化速度之間相互影響；供電方式、電流大小、填充比及電極浸入深度同樣影響渣池的溫度分布、金屬熔池的形成和渣皮的厚度等重要的工藝指標[4]，因此分析電渣重熔過程中的溫度場和流場對于控制相應的工藝參數提高重熔鑄錠質量有重要的意義。由于ESR過程處于高溫大電流環境，通過實驗測量手段對其流場溫度場進行研究存在一定困難，因此近年來國內外許多學者通過數值模擬方法在ESR研究方面做了很多工作[5-10]。

從數值模擬方法上看，景馨等使用有限元軟件ANSYS計算電磁場和溫度場[11-14]；文獻[15-17]中利用ANSYS軟件計算出焦耳熱和洛倫茲力，將二者導入流體計算軟件FLUENT中，分別作為能量方程和動量方程的源項來計算溫度場和流場,沒有考慮流體運動對磁場的影響；文獻[18-19]中利用Meltflow、CAE軟件計算電磁場和溫度場；文獻[5，7，20-22]中采用有限差分方法模擬具體條件下電渣重熔過程的計算程序；文獻[7，23]中求解外加磁場二維傳輸方程，也沒有考慮流動對磁場的影響；電渣重熔體系內流體的流動屬于磁流體流動范疇，本文采用流體計算軟件FLUENT中的磁流體模塊(Magneto hydro dynamics，MHD)，求解由于流動引起的附加磁場傳輸方程，因交變電流產生的外部磁場由有限元計算軟件ANSYS求解，然后通過C語言編寫的轉換程序將外部磁場數據導入到FLUENT中，由外部磁場和附加磁場共同產生的焦耳熱和電磁力作為源項加入到能量方程和動量方程中，最終由FLUENT耦合求解流場及溫度場。

1　電渣重熔過程

電渣重熔的基本原理如圖1所示。在銅制水冷結晶器內盛有熔融的爐渣，自耗電極一端插入熔渣內。自耗電極、渣池、金屬熔池、鋼錠、底水箱通過短網導線和變壓器形成回路。在通電過程中，渣池放出焦耳熱，將自耗電極端頭逐漸熔化，熔融金屬匯聚成液滴，穿過渣池，落入結晶器，形成金屬熔池，受水冷作用,迅速凝固形成鋼錠。在電極端頭液滴形成階段,以及液滴穿過渣池滴落階段,鋼-渣充分接觸，鋼中非金屬夾雜物及有害元素為爐渣所吸收，從而達到凈化金屬熔液的目的，同時渣池可以隔絕空氣起到精煉作用。上升的熔池在鋼錠外層形成一層薄渣殼，不僅使鋼錠表面光潔，還起到絕緣和隔熱作用，使更多的熱量向下部傳導，有利于鋼錠自下而上的定向結晶[23-24]，等到結晶器內的液相完全凝固后，可以進行脫錠，完成整個生產過程。

圖1　電渣重熔工作原理圖

2　數學模型

數值模擬就是數值求解描述物理問題的微分或積分方程組，而數學模型就是描述物理問題的控制方程及邊界條件。描述電渣重熔過程所用基本控制方程包括麥克斯韋方程組、質量守恒方程(連續方程)、動量守恒方程和能量守恒方程。電渣重熔與普通流動傳熱問題最大的不同在于導電流體內部產生電磁場，流體受到電磁力并產生焦耳熱，流動的驅動力就是電磁力以及焦耳熱差異導致的浮升力，所以電磁力和浮升力以源項形式加到動量方程中，焦耳熱也以源項形式加到能量方程中。為了簡化計算過程，假設：① 忽略電極運動，電極端頭采用平頭形狀；② 側面的渣殼起到電絕緣作用；③ 渣和鋼液的密度、渣的電導率隨溫度變化，其他屬性均為常數[17]；④ 不考慮磁擴散的影響；⑤ 不考慮鋼錠周圍空氣隙引起的接觸熱阻；⑥ 重熔過程溫度超過居里點，鋼錠及鋼液均從鐵磁體變為順磁體，相對磁導率為1。

2.1　電磁場控制方程

麥克斯韋方程組由如下4個方程組成[25]：

(1)

(2)

(3)

(4)

電流連續方程(電荷守恒)：

(5)

電磁場本構關系：

(6)

D=εE

(7)

歐姆定律：

J=σE

(8)

式中：H為磁場強度，A/m；B為磁感應強度，T；E為電場強度，V/m；D為電位移，C/m2；J為電流密度，A/m2；μ為磁導率，H/m；ε為電容率，F/m；σ為電導率，1/(Ω·m)；q為體電荷密度，C/m3；t為時間，s。

電渣重熔使用交流電，交變電流產生的磁場和電場相互影響，可用麥克斯韋方程組描述，由于渣、金屬熔池和鑄錠充分導電且為低頻狀態，可忽略體電荷密度和位移電流[25],式(3)忽略，式(1)變為：

(9)

式(5)變為：

(10)

考慮到帶電流體在磁場中的速度效應，式(8)變為：

J=σ(E+U×B)

(11)

式中：U為流體速度矢量，m/s。

由麥克斯韋方程組、電流連續方程及歐姆定律得到導電流體電磁場的磁感應強度傳輸方程：

(12)

求解出磁感應強度后，根據式(9)(6)求得電流密度。

如何求解磁感應強度傳輸方程是解決電磁場問題的關鍵，考慮到磁場由外加磁場和流體運動產生的感應磁場共同作用而成，得到：

B=B0+b

(13)

式中：B0為外加感應磁場，與流體運動無關；b為由于帶電流體運動而產生的感應磁場,也就是磁流體的速度效應。

B0的傳輸方程為：

(14)

b的傳輸方程為：

(15)

外加感應磁場B0由有限元計算軟件ANSYS通過矢量磁位法求解，磁場和電場與矢量磁位的關系[26]：

(16)

(17)

式中：A為矢量磁位，(Vs)/m；φ為標量電位，V。

利用SOLID97單元求得矢量磁位和標量電位，再由式(16)得到磁感應強度B0。將B0的實部和虛部按照FLUENT MHD模型要求的格式導入到流體計算軟件FLUENT中，由FLUENT求解b的傳輸方程。在FLUENT中通過編寫UDF，求解電流密度，電磁力(洛倫茲力)及焦耳熱：

(18)

Fe=J×B

(19)

(20)

式中：Fe為電磁力，N/m3；Q為焦耳熱，W/m3。

2.2　流場控制方程

描述流體運動的連續方程和動量守恒方程在電渣重熔問題中同樣適用，只是流動主要由電磁力和浮升力引起的，屬于湍流流動，但湍流強度不大；假設粘度為常數，密度使用Boussinesq假設，該假設只是用于計算浮升力Fb，方程中其他項中的密度為常數。另外由于鋼液發生凝固，在固液兩相區流動受到阻力，作為多孔介質區域處理(也稱焓-多孔介質法)，所以在動量方程及湍流方程中要加上相應的源項。

連續方程：

(21)

動量方程：

(22)

Fb=-ρ0β(T-T0)g

(23)

(24)

式中：ρ0為T0時的密度，kg/m3；T0、T分別為參考溫度和流體溫度，K；β為熱膨脹系數；μeff為有效黏性系數，kg/(m·s)；p為壓強，N/m2；q為體電荷密度，c/m3，對于高電導率的低頻條件，體電荷密度可忽略不計；Fe為電磁力，N/m3；g為重力加速度，m/s2；Sp為由于凝固引起的流動阻力，N/m3；fl為固液兩相區液相比例,Amush為兩相區常數，為了使液固界面清晰設為107[15]；Up為鑄錠拉速,m/s，實為渣層上升速度，該選項利于連續方程收斂，但由于數量級很小可以忽略[27]。

湍流模型主要用來計算有效粘性系數，選擇RNGk-ε湍流模型及提高型近壁面函數，該模型適用于低雷諾數的湍流流動。計算域中包括渣和鋼液屬于多相流，選擇VOF模型，該模型適合于兩相間有明顯的分界面，針對第二相求解體積比傳輸方程，兩相的體積比之和為1，設鋼液為第二相。表面張力影響著液滴的形成和鋼渣界面形狀，采用連續性表面張力模型。

2.3　能量方程

通過求解能量方程得到溫度場。由于發生相變，焓中要考慮凝固釋放的潛熱，固液模糊區的焓值通過液相比例設定。焦耳熱作為源項加到能量方程中，

(25)

H=h+flL

(26)

(26)

fl=(T-TS)/(TL-TS)

(27)

式中：keff為有效導熱系數，W/(m·K)，href為參考溫度下的參考焓值，J/kg；L為鋼液凝固釋放出來的潛熱，J/kg；Cp為鋼液比熱，J/(kg·K)；TS鋼液固相線溫度，K；TL鋼液液相線溫度，K。

2.4　邊界條件

ANSYS求解外加電磁場，利用矢量磁位法將xyz3個方向上的磁位和電位作為未知量，電極、渣、鋼以及外圍空間為計算域。電極一端為電流邊界條件，耦合該面上電壓，鋼錠底部電壓約束為0，周圍空間外表面為磁力線平行邊界條件。

用FLUENT進行多物理場耦合計算，計算域不包含電極部分，只包含渣池和金屬域。以電極與渣池接觸面為液態金屬速度入口邊界，溫度高出液相線溫度20 K。雖然金屬熔池的焦耳熱是周期性變化的，但因其周期很短(為0.02 s)，可以認為金屬熔滴是勻速滴落的；取鋼錠下表面為流量出口邊界；渣池與空氣接觸面采用剪應力為零邊界條件，換熱方式為對流和輻射混合模式；其他面為墻壁，采用無滑移壁面邊界條件，換熱采用對流換熱。在FLUENT中通過自定義標量(UDS)求解由于流體運動產生的感應磁場，壁面為絕緣壁。邊界條件參數見表1。

表1　邊界條件具體參數

3　模型調試

3.1　幾何模型及網格劃分

ANSYS計算域包括電極、渣池、金屬鑄錠及其周圍空間，如圖2、3所示。取空間半徑約為渣池半徑的10倍左右，電極高度對計算結果影響不大。電極、渣池及金屬鑄錠為核心計算域，為了減少計算量，核心域網格尺寸比外圍空間小得多。核心域采用SWEEP方式劃分六面體(HEX)網格，網格尺寸10 mm；外圍空間采用四面體網格，網格尺寸0.1 m。FLUENT計算域只包括渣池和鑄錠，全部為六面體網格。幾何模型尺寸：電極直徑為55 mm；鑄錠直徑為0.12 m；外圍空間直徑為1.2 m；渣池高度為60 mm；鑄錠高度為0.215 m；電極浸入深度為2 mm；電極高度為0.1 m。

圖2　ANSYS計算域幾何模型圖3　FLUEN計算域幾何模型

3.2　MHD模型及外加磁場導入轉換

FLUENT中的MHD模型包括兩種求解電磁場方法，一是磁感應法，求解由于流體運動引起的感應磁場；二是電位法，忽略流體運動引起的感應磁場，且外加磁場是由直流電場產生的靜態磁場。電渣重熔使用交流電，且考慮磁流體的速度效應，屬于流動導體介質內的磁感應問題。通過定義自定義標量uds-0(bx)、uds-1(by)、uds-2(bz)，求解其傳輸方程，再通過自定義函數UDF求解電場強度、電磁力以及焦耳熱等并存儲在自定義內存(UDM)中。FLUENT耦合求解多個控制方程，其求解機理見文獻[28]。

外加磁場由ANSYS計算得到，MHD模型需要的實部虛部數據，所需數據域為一個立方體，該立方體的大小要把FLUENT計算域包含在內。立方體需要劃分多個網格，其網格尺寸與ANSYS網格尺寸差不多。立方體大小的定義、網格劃分以及得到每個網格節點的磁場數據都由C語言編寫的數據轉換程序完成。圖4中C為立方體中任意網格，其磁場數據可由加權法或插值法得到，兩種方法得到的磁場數據基本相同，見圖5、6。

圖4　數據轉換示意圖

加權法認為一個網格內包含有多個ANSYS節點數據，網格中心磁場由這些點疊加得到，且離中心越遠的貢獻的磁場越小，權重與距離倒數成正比：

(28)

(29)

式中：BC為網格中心磁場數據；Bi為Pi點磁場數據，由ANSYS計算得到；vi為Pi點權重，Li為Pi點到網格中心的距離，m。

圖5　加權法得到的磁場數據圖6　插值法得到的磁場數據

插值法即定義好網格C的中心點坐標后，找到其周圍8個方向上離中心最近的8個點進行插值。如點pmp的順序為xyz，p表示增加方向(plus)，m表示減小(minus)，pmp表示x方向增加、y方向減小、z方向增加。找到8個點后進行插值，插值順序從邊長小的方向開始，如△x<△y<△z,插值則先從x方向插值，得到1、2、3、4個點的值；再進行y方向插值，由1、2得到5點值，由3、4得到6點值；最后z方向插值，由5、6得到網格中心的值。1點的插值式為：

(30)

其他點雷同。

式中Bmpm及Bppm為mpm點和ppm點的磁場數據，由ANSYS計算得到。

3.3　關鍵參數

渣池內熱浮升力的計算采用Boussinesq模型，需要設置參考溫度以及該溫度下熔渣的參考密度，參考溫度設置的是否合理將影響流場分布[29]。時間步長影響液滴形成及低落過程，一般采用較小的時間步長，如0.001。金屬的導熱系數影響溫度場的分布，最終會影響金屬熔池的形貌，文中采用piecewise-polynomial方法將固、液兩相區域分開設置。渣池焦耳熱與電渣的電導率有關，電導率隨溫度變化而變化，最終會影響溫度場的分布。具體物性參數見表2。

3.4　初始化

合理的初場設置有利于控制方程的收斂并縮短計算時間，采用hybrid方式初始化有利于得到合理的速度初場。因外加磁場產生的焦耳熱是周期性變化的，沒有必要進行溫度初場的計算，直接采用較小時間步長(步長可為周期的1/20)同時計算所有方程。

4　結果分析

4.1　MHD模型

MHD模型的計算結果中包括外加磁場、由于運動產生的磁場、電流密度、電場及洛倫茲力及焦耳熱，見圖7～12。從結果中可以看出電流密度和外加磁場方向符合右手螺旋定則，洛倫茲力與電流方向及磁場方向之間符合右手定則；電場強度方向與電流方向一致，由于金屬的電導率要比電渣的電導率大得多，所以金屬域的電場強度要遠遠小于渣池的電場強度；由于某點磁場是電荷在該點產生的磁場疊加合成，從圖7中可以看出沿半徑方向磁感應強度增加且呈線性變化，由于集膚效應進入到熔池的電流密度在徑向分量可以忽略，另外由于渦流效應渣池部分的磁感應強度大于鋼錠部分的磁感應強度[12]；從圖8可以看出，由于流體在磁場中的運動而產生的附加磁場同外加磁場相比很小很小，完全可以忽略不計，其磁雷諾數為4.16×10-5，符合文獻[7]中當磁雷諾數遠遠小于1時不考慮運動附加磁場的結論；從圖9可以看出，電流進入到熔池后，流經渣池、金屬熔池、鋼錠后返回到變壓器，在渣域由于金屬熔滴電阻遠小于渣池電阻，所以熔滴電流密度較大；從圖11可以看出，洛倫茲力方向指向對稱軸，且在電極底部周圍附近最大，越靠近對稱軸越小，有利于渣池形成中間向下四周向上的漩渦流動；從圖12可以看出，在渣池內電極下方位置的焦耳熱最大，從而有利于自耗電極的熔化滴落，另外金屬域的焦耳熱要比電渣域的小的多，這是因為金屬的電阻遠遠小于電渣的電阻；文中MHD模型計算結果與文獻[23]基本一致。

表2　物性參數表

圖7　外加磁場矢量圖圖8　流動附加磁場矢量圖

圖9　電流密度矢量圖圖10　電場強度矢量圖

圖11　洛倫茲力矢量圖圖12　焦耳熱云圖

4.2　流場及熔滴

渣池內的流動主要由洛倫茲力和浮升力引起，圖13為渣池流動的速度場。從圖中可以看出，由洛倫茲力引起的渣的流動呈漩渦狀，中間部分向下四周向上，

圖13　渣池速度場

且在電極底部周圍附近速度較大，洛倫茲力的攪拌效果增強，有利于熔池加深及渣池溫度均勻化[18]；金屬熔池內也有漩渦狀流動，流動與其上部渣池漩渦流動方向相反，是交界面的渣和滴落的鋼液從內向外的流動引起的熔池漩渦，洛倫茲力對金屬熔池影響很小；在靠近壁面處存在流速向下的漩渦流動，這是由于近壁面處溫度較低內外部密度差形成的浮升力所引起的漩渦流動，該渦流方向與洛倫茲力引起的渦流流動方向相反；在液滴的帶動下，在靠近中心處的渣速變大。

圖14為金屬熔滴在不同時刻形成長大滴落過程。在焦耳熱的作用下，自耗電極端頭開始熔化，液態金屬從端頭周圍向中心聚集，在重力、表面張力以及洛倫茲力的的共同作用下形成熔滴。熔滴經過長大、脫離母體、分散成小熔滴、緩慢穿透渣層等過程進入到金屬熔池，并在渣金界面引起擾動。

t=1st=1．2st=1．3st=1．4st=1．5s

圖14熔滴滴落過程

4.3　溫度場及凝固

從計算中可知傳熱過程在1 840 s之后基本趨于穩定，圖15為1 848 s時的溫度場。由于焦耳熱的不同，造成金屬域和渣域溫度場差距很大，在電極底端周圍電渣溫度較高，另外由于洛倫茲力引起的漩渦流動使得渣池上部電極下方保持一個較高溫度區域，這樣有利于電極熔化。在金屬熔池下方鋼錠區域，金屬內部傳熱為導熱，壁面與外界傳熱為對流，等溫線呈階梯狀變化，鋼錠底部周邊溫度最低。圖16為1 848 s時金屬熔池剖面形貌圖，熔池深度約40 mm，從中心到壁面凝固邊界由兩端弧線組成，在靠近壁面處弧線緩慢降低，拐點之后弧度變大，中心部位較平幾乎沒有弧度。

圖15　熔池溫度場圖16　金屬熔池剖面形貌圖

5　結　語

電渣重熔的流場和溫度場數值計算，其重點是如何解決由電磁場引起的洛倫茲力和焦耳熱對流動和傳熱的影響。文中應用FLUENT軟件的MHD模型中在流動導體內外加磁場法，求解由于磁流體運動而產生的附加磁場傳輸方程，并將洛倫茲力和焦耳熱分別作為源項加入的動量方程和能量方程中。其外部磁場由ANSYS軟件通過矢量磁位法計算得到，使用插值法或加權法編寫磁場數據轉換程序導入到MHD模型中。使用VOF模型進行渣金界面跟蹤，渣池熱浮升力采用Boussinesq模型，用焓-多孔介質法處理金屬凝固模糊區，最后進行多方程同時耦合求解。

計算結果表明，由于磁流體運動引起的附加磁場很小很小，可以忽略不計，所以計算電渣重熔的電磁場時可以不考慮速度效應；渣池內流動主要受洛倫茲力的影響，熱浮升力只對近壁面處流動產生影響；電磁場產生的焦耳熱在渣域遠遠大于金屬域，渣池大量的焦耳熱足以熔化自耗電極，形成穩定的金屬熔滴不斷滴落；非穩態計算在1 840 s之后在溫度場接近穩定，金屬熔池形成穩定的凝固邊界。

參考文獻(References)：

[1]姜周華，陳旭，董艷伍，等.電渣冶金技術的創新與發展[J].世界金屬導報，2015，11(3)：1-5.

[2]李小孟，劉立，趙俊學，等. 電渣重熔冶金過程爐渣粘度性能的研究[J].金屬世界，2015(5)：56-59.

[3]Eremin E N, Zherebtsov S N, Gritsevich V N. Usingcast electroslag blanks instead of forging in flange production [J]. Chemicaland Petroleum Engineering，2003，39：7-8.

[4]王春光，葛鋒，張玉碧，等. 工藝參數對電渣重熔凝固過程的影響綜述[J]. 鑄造技術，2013，34(10)：1321-1323.

[5]張華，關洋，陳瑞，等. 電渣熔鑄導葉數值模擬技術研究[J].鑄造，2013,62(7)：637-640.

[6]梁強，陳希春，付銳，等. 電流對GH4169合金電渣重熔過程影響的數值模擬研究[J].鑄造技術，2011,32(9)：1252-1256.

[7]DILAWARI A H, SZEKELY J. A mathematical model of slag and metal flow in the ESR process [J]. Metallurgical Transactions B，1977，8B (6)：227-236.

[8]Kelkarl K M, Patankarl S V, Srivatsa S K.Computational modeling of electro slag remelting (ESR) process used for the production of high-preformance alloys[C]//Proceedings of the 2013 International Symposium on Liquid Metal Processing and Casting. USA:Wiley, 2013:1-10.

[9]Fezi K, Yanke J, Krane M J M. Modeling macrosegregation during electrslag remelting of alloy 625 [C]//Proceedings of the 2013 International Symposium on Liquid Metal Processing and Casting. USA:Wiley, 2013:151-158.

[10]張赫，雷洪，耿佃橋，等. 電渣重熔過程傳熱及凝固組織的數值模擬[J].工業加熱，2013,42(6)：42-46.

[11]景馨，姜周華，劉福斌，等. 導電結晶器電渣重熔空心鋼錠電場和溫度場的數值模擬[J].材料與冶金學報，2013,12(3)：203-208.

[12]王芳，李寶寬. 電渣重熔過程中的電磁場和Joule熱分析[J].金屬學報，2010,46(7)：794-799.

[13]段吉超，董艷伍，姜周華，等. 工藝參數對電渣重熔過程影響的數值模擬研究[J].江西冶金，2014,34(1)：5-8.

[14]劉福斌，李永旺，姜周華，等. 實驗室規模和工業規模電渣重熔過程中電磁行為[J].中南大學學報，2015,46(10)：3580-3585.

[15]賀鑄，蔡輝，劉雙，等. 不同熔化速率下的電渣重熔過程數值模擬分析[J].重慶大學學報，2015,46(10)：104-110.

[16]王曉花，厲英. 電渣重熔過程渣池流場數值模擬[J].材料與冶金學報，2014,13(2)：133-136.

[17]劉艷賀，賀鑄，劉雙，等. 電渣重熔過程中電磁與流動及溫度場的數值模擬[J].過程工程學報，2014,14(1)：16-21.

[18]梁強，陳希春，郭漢杰. 熔速對GH4169合金電渣重熔凝固過程參數影響的數值模擬研究[J].鑄造技術，2013,34(8)：1016-1020.

[19]唐建軍，耿茂鵬，饒磊. 電渣重熔中電極熔化過程的模擬研究[J].鑄造技術，2010,31(8)：1062-1065.

[20]劉喜海，王君卿，賈維國，等. 核電用大型鋼錠電渣重熔工藝模擬及生產應用[J].鑄造，2010,59(12)：1315-1319.

[21]白海軍，厲英，李寶寬. 電磁攪拌對電渣重熔鋼錠溫度場的影響[J].特殊鋼，2008,29(1)：7-9.

[22]王長周，宋錦春. 電渣重熔爐電極熔化速度對重熔過程的影響[J].冶金設備，2012,4(2)：5-8.

[23]Antje RüCKERT, Herbert PFEIFER. Numerical modelling of the electroslag remelting process [J]. Metal，2007，5：22-24.

[24]耿茂鵬，孫達昕.電渣熔鑄過程控制與模擬仿真[M]．北京:冶金工業出版社，2008．

[25]Moreau R. Magnetohydrodynamics[M]. Kluwer Academic Publishers. 1990.

[26]閻照文.ANSYS工程電磁分析技術與實例詳解[M].北京:中國水利水電出版社，2008．

[27]Jardy A, Ablitzer D, Wadier J F. Magnetohydrodynamic and thermal behavior of electroslag remelting slags[J]. Metall Mater Trans, B 22B(1991)：111-120.

[28]韓麗輝，于春梅，馮根生. 數值模擬質子交換膜燃料電池的方法[J].實驗技術與管理，2014,31(8)：42-43.

[29]胡進，孫本良，李成威，等.鋼液密度對連鑄中間包流場影響的數值模擬[J].冶金能源，2008,27(6)：12-14.