導(dǎo)軌式電磁發(fā)射裝置電樞熔化波有限元計(jì)算

譚賽,魯軍勇,張曉,關(guān)曉存,龍?chǎng)瘟?/p>

(海軍工程大學(xué)艦船綜合電力技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室,430033,武漢)

?

導(dǎo)軌式電磁發(fā)射裝置電樞熔化波有限元計(jì)算

譚賽,魯軍勇,張曉,關(guān)曉存,龍?chǎng)瘟?/p>

(海軍工程大學(xué)艦船綜合電力技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室,430033,武漢)

針對(duì)單元燒蝕算法精度較低的問(wèn)題,基于導(dǎo)軌式電磁發(fā)射裝置的二維控制方程,推導(dǎo)出采用求解量移動(dòng)法處理運(yùn)動(dòng)問(wèn)題的伽遼金有限元離散方程,建立了電磁-溫度-運(yùn)動(dòng)耦合場(chǎng)有限元計(jì)算模型。以單元節(jié)點(diǎn)為燒蝕判斷對(duì)象,采用節(jié)點(diǎn)燒蝕算法處理電樞材料燒蝕問(wèn)題,建立了電樞熔化波模型,對(duì)運(yùn)動(dòng)電樞的燒蝕速度進(jìn)行了數(shù)值計(jì)算。模型的計(jì)算結(jié)果表明:在激勵(lì)電磁感應(yīng)強(qiáng)度為40 T、電樞運(yùn)動(dòng)速度為150 m/s的情況下,采用節(jié)點(diǎn)燒蝕法得到的熔化波燒蝕速度為Barber理論模型計(jì)算值的94%,而相關(guān)文獻(xiàn)采用單元燒蝕法的計(jì)算值為Barber理論模型計(jì)算值的73%。因此,與采用單元燒蝕法相比,采用節(jié)點(diǎn)燒蝕法的計(jì)算值與Parks和Barber兩個(gè)經(jīng)典理論模型的計(jì)算值更加相近,驗(yàn)證了該方法的正確性。

導(dǎo)軌式電磁發(fā)射裝置;運(yùn)動(dòng)邊界;電樞;熔化波;燒蝕速度;耦合場(chǎng);相變

在脈沖電流的驅(qū)動(dòng)下,導(dǎo)軌式電磁發(fā)射裝置的電樞在兩軌道間做加速運(yùn)動(dòng),電樞與軌道存在著滑動(dòng)電接觸,電樞的高速運(yùn)動(dòng)對(duì)軌道及電樞內(nèi)部的磁場(chǎng)、電流分布具有顯著的影響[1-2],這種現(xiàn)象稱(chēng)為速度趨膚效應(yīng)。文獻(xiàn)[3-4]采用有限差分法建立了二維電磁場(chǎng)、溫度場(chǎng)數(shù)學(xué)離散模型,文獻(xiàn)[5-6]則采用有限元法建立計(jì)算模型,得到了電樞運(yùn)動(dòng)下的電流分布特點(diǎn)。

由于電流集中于電樞尾翼,所以在焦耳熱的作用下,電樞材料最先在尾翼處發(fā)生熔化,熔化部分在電磁力的作用下迅速?gòu)碾姌斜倔w剝離,隨后熔化區(qū)域在新的電樞邊界處重新形成,這種熔化形式稱(chēng)為熔化波。對(duì)于電樞熔化燒蝕問(wèn)題,國(guó)外學(xué)者Parks和Barber分別從理論上推導(dǎo)出了二維情況下的電樞熔化燒蝕速度[7-8]。此后,國(guó)內(nèi)外學(xué)者分別從一維、二維模型出發(fā),采用有限差分法、有限元法建立了基于單元燒蝕算法的電樞燒蝕模型,得出了不同電樞速度下的熔化波燒蝕速度[9-12]。從公開(kāi)文獻(xiàn)中的分析結(jié)果可知:采用單元燒蝕算法時(shí),若電樞運(yùn)動(dòng)速度較大,則其計(jì)算結(jié)果與理論值偏差過(guò)大。本文認(rèn)為:選取單元作為燒蝕判斷對(duì)象時(shí),實(shí)際上是取單元內(nèi)節(jié)點(diǎn)的溫度算術(shù)平均值作為單元的溫度,在單元內(nèi)溫度梯度較大的情況下,勢(shì)必會(huì)產(chǎn)生較大的誤差,甚至可導(dǎo)致相變過(guò)程的數(shù)值計(jì)算方法發(fā)生失效。

為此,本文首先建立了考慮電樞運(yùn)動(dòng)及材料相變的電磁-溫度耦合場(chǎng)的二維有限元模型,并針對(duì)單元燒蝕算法應(yīng)用在電樞熔化燒蝕計(jì)算中時(shí)精度較差的問(wèn)題,采用節(jié)點(diǎn)燒蝕法建立電樞熔化波模型,對(duì)運(yùn)動(dòng)電樞的燒蝕速度進(jìn)行數(shù)值計(jì)算,得到了與經(jīng)典理論計(jì)算值更加接近的電樞熔化燒蝕速度。

1 導(dǎo)軌式電磁發(fā)射裝置控制方程

假設(shè)電樞與軌道間的接觸為理想接觸,不考慮電樞與軌道間的接觸電阻,當(dāng)電樞材料發(fā)生熔化時(shí),電樞材料立即在電磁力的作用下脫離電樞,電樞的實(shí)際物理邊界發(fā)生相應(yīng)變化,且忽略電樞非熔化區(qū)域的形變。

1.1 電磁場(chǎng)控制方程

求解域如圖1所示,采用靜止坐標(biāo)系,二維下電樞與軌道內(nèi)的磁擴(kuò)散方程[2-3]為

(1)

求解域邊界條件為

(2)

式中:σ、μ、vx分別為材料的電導(dǎo)率、磁導(dǎo)率及電樞的運(yùn)動(dòng)速度,均是與區(qū)域相關(guān)的參數(shù),假設(shè)軌道沿x軸的負(fù)方向運(yùn)動(dòng),則vx為軌道運(yùn)動(dòng)速度,而如果電樞靜止不動(dòng),則vx為0;Bz為電磁感應(yīng)強(qiáng)度,對(duì)于二維問(wèn)題,僅存在z方向的分量;J為施加的電流密度,其方向平行于求解域平面,對(duì)于邊界S2、S3,J僅存在x方向的分量,而對(duì)于邊界S1、S7,J僅存在y方向分量。模型中邊界S1、S2處施加的電磁感應(yīng)強(qiáng)度用于等效引入電流密度激勵(lì)。

圖1 導(dǎo)軌式電磁發(fā)射裝置二維計(jì)算模型

1.2 溫度場(chǎng)控制方程

如圖1所示,在二維情況下,電樞與軌道內(nèi)的熱傳導(dǎo)方程滿足

(3)

式中:κ、c、ρ、EL分別為材料的熱導(dǎo)率、比熱容、質(zhì)量密度和熔化潛熱。熔化潛熱即材料發(fā)生熔化相變過(guò)程中單位質(zhì)量所吸收的熱量。

在溫度場(chǎng)中,求解域的邊界條件為

(4)

1.3 網(wǎng)格移動(dòng)法處理運(yùn)動(dòng)問(wèn)題

對(duì)于運(yùn)動(dòng)電磁場(chǎng)問(wèn)題,為消除運(yùn)動(dòng)項(xiàng)系數(shù)矩陣對(duì)稱(chēng)性的影響,通常采用獨(dú)立坐標(biāo)系。求解域的運(yùn)動(dòng)部分與靜止部分分別采用不同的坐標(biāo)系,將控制方程中的運(yùn)動(dòng)項(xiàng)消去,通過(guò)網(wǎng)格的實(shí)際移動(dòng)來(lái)考慮運(yùn)動(dòng)項(xiàng)的影響[13]。

1.4 求解量移動(dòng)法處理運(yùn)動(dòng)問(wèn)題

假設(shè)電樞靜止,軌道沿x負(fù)方向運(yùn)動(dòng),通過(guò)移動(dòng)求解量考慮運(yùn)動(dòng)項(xiàng)的影響,其等效性說(shuō)明如下。

矢量B是時(shí)間和空間上的函數(shù),在二維情況下

(5)

由式(1)、式(5)可得

(6)

對(duì)于運(yùn)動(dòng)部分,運(yùn)動(dòng)坐標(biāo)系下的時(shí)間離散表達(dá)式為

(7)

因此,在t時(shí)刻的節(jié)點(diǎn)求解量Bz計(jì)算完成后,運(yùn)動(dòng)部分運(yùn)動(dòng)的距離為ds。在進(jìn)行t+Δt時(shí)刻的計(jì)算時(shí),需用到節(jié)點(diǎn)在t時(shí)刻時(shí)所處位置的求解量。采取的處理方法為:網(wǎng)格不發(fā)生移動(dòng),t時(shí)刻計(jì)算所得的節(jié)點(diǎn)求解量沿運(yùn)動(dòng)方向平移ds后作為t+Δt時(shí)刻計(jì)算時(shí)所用的求解量數(shù)值,即通過(guò)求解量的移動(dòng)來(lái)模擬電樞的運(yùn)動(dòng)。

1.5 有限元離散方程

根據(jù)求解區(qū)域的電磁場(chǎng)和溫度控制方程,采用移動(dòng)求解量法處理運(yùn)動(dòng)問(wèn)題,最終得到導(dǎo)軌式電磁發(fā)射裝置二維耦合場(chǎng)計(jì)算模型的伽遼金有限元離散方程為

(8)

(9)

采用四節(jié)點(diǎn)四邊形單元,則系數(shù)矩陣和右端向量在局部坐標(biāo)系下的形式為

(10)

(11)

(12)

(13)

(14)

式中

(15)

下角標(biāo)i、j為單元節(jié)點(diǎn)編號(hào);ξ、η為局部坐標(biāo)系下的坐標(biāo)軸;J為雅克比矩陣,表示整體坐標(biāo)系與局部坐標(biāo)系之間的轉(zhuǎn)換關(guān)系。

對(duì)于形如式(10)～式(14)中的數(shù)值積分,采用高斯積分法進(jìn)行計(jì)算。對(duì)于瞬態(tài)場(chǎng)問(wèn)題,采用Crank-Nicholson法對(duì)時(shí)間進(jìn)行離散。

2 電樞材料燒蝕的處理

電樞在加速運(yùn)動(dòng)過(guò)程中,由于電流在電樞尾翼處的集中,可導(dǎo)致此處的材料發(fā)生熔化,造成電樞的燒蝕,因此在溫度場(chǎng)計(jì)算過(guò)程中,需考慮材料的相變過(guò)程。本文采用熱源法,為避免單元溫度梯度較大導(dǎo)致的計(jì)算精度降低、甚至相變過(guò)程的數(shù)值計(jì)算失敗,以節(jié)點(diǎn)為判斷對(duì)象:當(dāng)單元內(nèi)的節(jié)點(diǎn)達(dá)到材料熔點(diǎn)時(shí),節(jié)點(diǎn)溫度不再發(fā)生變化;當(dāng)其吸收的熱量達(dá)到熔化潛熱時(shí),節(jié)點(diǎn)處發(fā)生相變,認(rèn)為節(jié)點(diǎn)及其附近發(fā)生了燒蝕。當(dāng)節(jié)點(diǎn)發(fā)生熔化時(shí),其對(duì)應(yīng)位置的電樞材料在電磁力的作用下迅速脫離電樞本體,磁場(chǎng)迅速擴(kuò)散至電樞的非熔化區(qū)域邊界,即電磁感應(yīng)強(qiáng)度的邊界條件發(fā)生了改變。由式(16)可知:磁場(chǎng)在材料中的滲透深度δ與材料電阻率ρ的1/2次方成正比。

(16)

在程序中增大熔化位置的電阻率,可增加磁場(chǎng)的滲透深度,使磁場(chǎng)較快地滲透至電樞的非熔化邊界,相當(dāng)于更新了Bz的邊界條件。在程序編寫(xiě)過(guò)程中發(fā)現(xiàn):將熔化位置的電阻率增大4個(gè)數(shù)量級(jí)時(shí)即可得到較為準(zhǔn)確的結(jié)果。因此,對(duì)于存在燒蝕節(jié)點(diǎn)的單元,采用非均勻材質(zhì)單元的處理方法,在高斯積分的過(guò)程中,將熔化節(jié)點(diǎn)的電阻率增大4個(gè)數(shù)量級(jí),這種方法稱(chēng)為節(jié)點(diǎn)燒蝕法。

具體的處理過(guò)程描述如下:圖2為正處于相變過(guò)程的單元,其中節(jié)點(diǎn)1、4為已發(fā)生燒蝕的節(jié)點(diǎn),節(jié)點(diǎn)2、3為非燒蝕節(jié)點(diǎn),在電磁場(chǎng)計(jì)算過(guò)程中,其單元系數(shù)矩陣ke的元素為

(17)

在溫度場(chǎng)計(jì)算過(guò)程中,其單元的右端向量Fe中的元素為

(18)

式中:n為每一維的高斯積分點(diǎn)數(shù),本文取n=2;σe(ξp,ηq)為單元內(nèi)高斯積分點(diǎn)(p,q)處的電導(dǎo)率,其具體數(shù)值與單元內(nèi)節(jié)點(diǎn)的相態(tài)有關(guān),對(duì)于燒蝕節(jié)點(diǎn)1、4,其對(duì)應(yīng)的高斯積分點(diǎn)g1、g4的電導(dǎo)率減小4個(gè)數(shù)量級(jí)。至此,完成了電樞燒蝕的數(shù)值模擬計(jì)算。

圖2 局部坐標(biāo)系下的單元高斯積分點(diǎn)示意圖

考慮材料相變的二維電樞熔化波有限元計(jì)算流程如圖3所示。

圖3 電磁-溫度-運(yùn)動(dòng)耦合場(chǎng)有限元計(jì)算流程圖

3 有限元計(jì)算結(jié)果分析與驗(yàn)證

為方便對(duì)比分析,本文模型所用的材料參數(shù)及求解域均與文獻(xiàn)[12]相同。求解域如圖1所示,所用材料的參數(shù)見(jiàn)表1。本文模型與文獻(xiàn)[12]中Stefani模型的求解區(qū)域單元剖分密度相同(每個(gè)單元邊長(zhǎng)為0.2 mm)。

在激勵(lì)磁感應(yīng)強(qiáng)度Bz為40 T(對(duì)應(yīng)的區(qū)域?yàn)镾1、S2邊界)、電樞運(yùn)動(dòng)速度vx為100 m/s的情況下,電樞的熔化波燒蝕演變?nèi)鐖D4所示,對(duì)應(yīng)的電流密度云圖如圖5所示。當(dāng)電樞材料發(fā)生熔化燒蝕時(shí),電流集中點(diǎn)從電樞尾翼處移至熔化波前端,形成了燒蝕間隙在電樞軌道接觸面由電樞尾翼至電樞前端的擴(kuò)展過(guò)程,燒蝕間隙的擴(kuò)展速度即為熔化波燒蝕速度。

表1 軌道電樞材料的參數(shù)

圖4 電樞熔化波燒蝕演變圖

圖5 熔化波燒蝕過(guò)程的電流密度

模型vm/m·s-1Bz=32T,v=100m/sBz=32T,v=150m/sBz=40T,v=100m/sBz=40T,v=150m/sParks10.115.215.823.7Barber8.612.813.420.0Stefani7.09.011.714.6本文7.59.412.218.8

由圖6a可知:對(duì)于特定的某一工況,本文模型計(jì)算得出的熔化波電樞燒蝕深度與時(shí)間近似成線性關(guān)系,即不同時(shí)刻的電樞熔化波燒蝕速度基本不變,這與Parks和Barber的理論公式相吻合。由圖6b可知:用Stefani模型計(jì)算出的燒蝕深度與時(shí)間成明顯的非線性關(guān)系,且在計(jì)算前期,電樞熔化波燒蝕速度計(jì)算值明顯偏小。

由表2的對(duì)比可知:當(dāng)Bz=40T、vx=100m/s時(shí),本文模型計(jì)算出的熔化波燒蝕速度略大于Stefani模型的熔化波燒蝕速度計(jì)算值,但小于Parks、Barber理論公式的計(jì)算值;當(dāng)Bz=40T、vx=150m/s時(shí),本文模型計(jì)算出的熔化波燒蝕速度為Barber理論公式計(jì)算值的94%,而Stefani模型的熔化波燒蝕速度計(jì)算值為Barber理論公式計(jì)算值的73%。

(a)本文模型計(jì)算結(jié)果

(b)Stefani模型計(jì)算結(jié)果[12]圖6 熔化波燒蝕深度與時(shí)間的關(guān)系

顯然,在燒蝕計(jì)算的整個(gè)過(guò)程中,Stefani模型的計(jì)算結(jié)果均偏小,主要原因在于:Stefani模型選取單元作為燒蝕判斷對(duì)象,實(shí)際上是取單元內(nèi)節(jié)點(diǎn)的溫度算術(shù)平均值作為單元的溫度,當(dāng)單元內(nèi)節(jié)點(diǎn)溫度差異較大時(shí),會(huì)導(dǎo)致其確定的單元燒蝕時(shí)刻滯后,最終造成計(jì)算前期材料燒蝕時(shí)刻延遲,且電樞熔化波燒蝕速度計(jì)算值偏小。

圖7 不同電樞速度下的熔化波燒蝕速度

取Bz=40T,并增加求解域的單元剖分密度(每個(gè)單元邊長(zhǎng)為0.05mm),計(jì)算不同電樞速度下的熔化波燒蝕速度,并與Parks和Barber理論公式的計(jì)算值進(jìn)行對(duì)比,結(jié)果如圖7所示。由圖7可知:在不同電樞速度下,本文模型計(jì)算出的熔化波燒蝕速度值介于Parks與Barber理論公式計(jì)算值之間。

由以上對(duì)比分析可知:與單元燒蝕法相比,在激勵(lì)磁感應(yīng)強(qiáng)度場(chǎng)較大、電樞運(yùn)動(dòng)速度較高的情況下,采用本文提出的節(jié)點(diǎn)燒蝕法可得到穩(wěn)定的、且與經(jīng)典理論公式計(jì)算值較吻合的燒蝕速度。

4 結(jié) 語(yǔ)

本文通過(guò)求解量移動(dòng)的方法處理運(yùn)動(dòng)邊界,采用節(jié)點(diǎn)燒蝕法處理電樞的燒蝕問(wèn)題,建立了考慮材料相變的導(dǎo)軌式電磁發(fā)射裝置二維電磁-溫度-運(yùn)動(dòng)耦合場(chǎng)計(jì)算模型,并進(jìn)行了電樞熔化波燒蝕速度的數(shù)值計(jì)算。計(jì)算結(jié)果表明:與單元燒蝕法相比,采用本文模型得到的熔化波燒蝕速度與Parks和Barber兩個(gè)經(jīng)典模型得到的理論計(jì)算結(jié)果更接近,間接驗(yàn)證了本文方法的正確性。

[1] ENGEL T G, NERI J M, VERACKA M J. Characterization of the velocity skin effect in the surface layer of a railgun sliding contact [J]. IEEE Transactions on Magnetics, 2007, 44(17): 1837-1844.

[2] GHASSEMI M, MOLAYI Y, HAMEDI M H. Analysis of force distribution acting upon the rails and the armature and prediction of velocity with time in an electromagnetic launcher with new method [J]. IEEE Transactions on Magnetics, 2007, 43(1): 132-136.

[3] 李聽(tīng), 翁春生. 電磁導(dǎo)軌炮固體電樞非穩(wěn)態(tài)電磁效應(yīng)研究 [J]. 南京理工大學(xué)學(xué)報(bào): 自然科學(xué)版, 2009, 33(1): 108-111. LI Xin, WENG Chunsheng. Unsteady electromagnetic effect on solid armature railguns [J]. Journal of Nanjing University of Science and Technology: Natural Science, 2009, 33(1): 108-111.

[4] HSIEH K T. Lagrangian formulation for mechanically, thermally coupled electromagnetic diffusive processes with moving conductors [J]. IEEE Transactions on Magnetics, 1995, 31(1): 604-609.

[5] LIU Hsingpang, LEWIS M C. 3-D electromagnetic analysis of armatures and rails for high launch energy applications [J]. IEEE Transactions on Magnetics, 2009, 45(1): 322-326.

[6] 王剛?cè)A, 謝龍, 王強(qiáng), 等. 電磁軌道炮電磁力學(xué)分析 [J]. 火炮發(fā)射與控制學(xué)報(bào), 2011(1): 69-71 WANG Ganghua, XIE Long, WANG Qiang, et al. Analysis on electromagnetic mechanics in electromagnetic railgun [J]. Journal of Gun Launch & Control, 2011(1): 69-71.

[7] PARKS P. Current melt wave model for transitioning solid armatures [J]. Journal of Applied Physics, 1990, 67(7): 3511-3516.

[8] BARBER J, DREIZIN Y. Model of contact transitioning with realistic armature rail interface [J]. IEEE Transactions on Magnetics, 1995, 31(1): 96-100.

[9] BENTON T, STEFANI F, SATAPATHY S, et al. Numerical modeling of melt wave erosion in conductors [J]. IEEE Transactions on Magnetics, 2003, 39(1): 129-133.

[10]MERRILL R, STEFANI F. Electrodynamics of the current melt wave erosion boundary in a conducting half-space [J]. IEEE Transactions on Magnetics, 2003, 39(1): 66-71.

[11]鞏飛, 翁春生. 固體電樞熔化波燒蝕的二維數(shù)值模擬 [J]. 南京理工大學(xué)學(xué)報(bào), 2012, 36(3): 487-491. GONG Fei, WENG Chunsheng. Two dimensional numerical simulation of melt wave erosion in solid armatures [J]. Journal of Nanjing University of Science and Technology, 2012, 36(3): 487-491.

[12]STEFANI F, MERRILL R, WATT T. Numerical modeling of melt wave erosion in two-dimensional block armatures [J]. IEEE Transactions on Magnetics, 2005, 41(1): 437-441.

[13]黃濤, 阮江軍, 張宇嬌, 等. 瞬態(tài)運(yùn)動(dòng)電磁問(wèn)題的時(shí)步有限元方法研究 [J]. 中國(guó)電機(jī)工程學(xué)報(bào), 2013, 33(6): 168-175. HUANG Tao, RUAN Jiangjun, ZHANG Yujiao, et al. Research on the time stepping element method for solving the transient motion electromagnetic problems [J]. Proceedings of the CSEE, 2013, 33(6): 168-175.

[14]WATT T, STEFANI F. The effect of current and speed on perimeter erosion in recovered armatures [J]. IEEE Transactions on Magnetics, 2005, 41(1): 448-452.

(編輯 葛趙青)

Finite Element Analysis of Melt Wave Ablation in Electromagnetic Rail Launcher Armatures

TAN Sai,LU Junyong,ZHANG Xiao,GUAN Xiaocun,LONG Xinlin

(National Key Laboratory for Vessel Integrated Power System Technology, Naval University of Engineering, Wuhan 430033, China)

Aiming at the low accuracy of element melt method in predicting the wave melt ablation (WMA) velocity, a new finite element method of melt wave ablation is proposed. Based on a two-dimensional governing equation of electromagnetic rail launcher, a Galerkin-form finite element discrete formulation for dealing with moving boundary by potential moving method is derived, then the finite element model of transient electromagnetic-thermal-moving coupling problem in two dimensions is established. Using node element as the criterion, the armature WMA model is built. Computational results show that when the armature moving speed is 150 m/s and the applied magnetic induction intensity is 40 T, the WMA velocity predicted by the node melt method (NMM) is 94% of that predicted by Barber’s models, while the WMA velocity predicted by Stefani using element melt method(EMM) is 73% of that predicted by Barber’s models. Thus, compared with EMM, the WMA velocity predicted by NMM is closer to the results of Parks’ and Barber’s classic models, which verifies the correctness of the proposed NMM method.

electromagnetic rail launcher; moving boundary; armature; wave melt ablation velocity; coupling field; phase change

10.7652/xjtuxb201603017

2015-08-19。 作者簡(jiǎn)介:譚賽(1988—),男,博士生;魯軍勇(通信作者),男,研究員,博士生導(dǎo)師。 基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(51207162,51522706,51407191);國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃資助項(xiàng)目(613262)。

TM315

:A

:0253-987X(2016)03-0106-06