電磁軌道發(fā)射器動態(tài)發(fā)射過程的數(shù)值模擬*

譚 賽，魯軍勇，張 曉，龍鑫林

(海軍工程大學(xué) 艦船綜合電力技術(shù)國防科技重點實驗室， 湖北 武漢 430033)

電磁軌道發(fā)射器動態(tài)發(fā)射過程的數(shù)值模擬*

譚 賽，魯軍勇，張 曉，龍鑫林

(海軍工程大學(xué) 艦船綜合電力技術(shù)國防科技重點實驗室， 湖北 武漢 430033)

針對電磁軌道發(fā)射器動態(tài)發(fā)射過程的數(shù)值模擬問題，基于矢量磁位A和時間積分標(biāo)量電位V，采用節(jié)點單元法，并選擇運動坐標(biāo)系描述運動問題，推導(dǎo)出動態(tài)發(fā)射下的電磁軌道發(fā)射器三維渦流場有限元離散方程。結(jié)合溫度場控制方程，建立電磁軌道發(fā)射器的三維電磁-溫度耦合有限元模型。針對一個電磁軌道發(fā)射器動態(tài)發(fā)射問題，對模型進行數(shù)值模擬，得到了動態(tài)發(fā)射下軌道及電樞的溫度場分布特點及發(fā)射器電感梯度隨時間的變化規(guī)律。計算結(jié)果表明，模型計算出的脈沖電流峰值、出口電流大小、出口速度等參數(shù)均與試驗結(jié)果吻合較好，驗證了所開發(fā)的有限元程序代碼的正確性。

動態(tài)發(fā)射；電樞；電磁軌道發(fā)射裝置；電感梯度；有限元

在脈沖電流的驅(qū)動下，電磁軌道發(fā)射器中的電樞在兩軌道間作加速運動，電樞與軌道存在著滑動電接觸，電樞的高速運動對軌道及電樞內(nèi)部的磁場、電流分布具有顯著的影響，這種現(xiàn)象稱之為速度趨膚效應(yīng)[1-4]。這種效應(yīng)將導(dǎo)致轉(zhuǎn)捩等一系列現(xiàn)象的發(fā)生，嚴重影響軌道的壽命。因此，電磁軌道發(fā)射器動態(tài)發(fā)射過程的數(shù)值計算一直是國內(nèi)外學(xué)者研究的熱點。然而，滑動電接觸問題具有其自身的特殊性，主要表現(xiàn)為：運動分界面處需同時滿足磁場強度H、磁感應(yīng)強度B、電場強度E、電流密度J等四個物理量的分界面條件。

為解決這一類滑動電接觸問題，一部分學(xué)者采用基于有限差分法建立了二維狀態(tài)下的電磁場、溫度場數(shù)學(xué)離散模型，分析了電樞運動下的電流分布特點[5-6]。然而，由于傳統(tǒng)有限差分法的固有缺點，難以對復(fù)雜結(jié)構(gòu)的電樞進行精確的分析。

一部分學(xué)者基于節(jié)點元法，在兩個不同導(dǎo)體的分界面處采用雙節(jié)點法，選擇非連續(xù)的矢量磁位A處理電樞的運動問題[7-8]。但這種方法并未在三維運動渦流問題中得到驗證；另一部分學(xué)者利用棱邊元法，選取H作為求解量，建立了電磁軌道發(fā)射器的三維有限元模型[9-10]，并引入迎風(fēng)因子處理電樞運動問題，在考慮電樞運動的同時，保持了網(wǎng)格拓撲結(jié)構(gòu)不變。然而，由于迎風(fēng)因子的引入，電樞速度較大時導(dǎo)致的Peclet數(shù)增大會嚴重影響模型的收斂性和準(zhǔn)確性。美國Hsieh等學(xué)者開發(fā)了具備進行電磁軌道發(fā)射器動態(tài)發(fā)射下的多物理場數(shù)值計算能力的三維機電分析程序代碼(Electro-Mechanical Analysis Program in three Dimensions， EMAP3D)有限元程序代碼[11]。此后，EMAP3D分別通過與三維顯式動力學(xué)程序代碼(DYNAmic in three Dimensions, DYNA3D)耦合、引入彈塑性本構(gòu)模塊，彌補了電樞的剛體假設(shè)不足[12-15]，針對計算規(guī)模需求不斷增加的需求，EMAP3D又推出了并行計算版本，使其仿真效率大幅提升[16-17]。目前，EMAP3D已經(jīng)比較完備，是迄今為止較為成熟的一款電磁軌道發(fā)射器動態(tài)發(fā)射數(shù)值模擬軟件。

但是，從文獻[18]可以看出，EMAP3D在進行電磁軌道發(fā)射器動態(tài)發(fā)射的數(shù)值模擬過程中，仍存在電樞側(cè)總電流計算值與實際電流激勵值誤差偏大的問題。因此，EMAP3D在處理樞軌接觸面處的分界面條件上仍存在一定的不足。本文基于矢量磁位與標(biāo)量磁位結(jié)合方法[19]，推導(dǎo)出了電磁軌道發(fā)射器的三維運動渦流場有限元方程，最終建立了動態(tài)發(fā)射下的電磁軌道發(fā)射器的電磁-溫度耦合模型。并針對一個電磁軌道發(fā)射器的動態(tài)發(fā)射問題進行了數(shù)值計算，將計算結(jié)果與EMAP3D程序計算結(jié)果、試驗結(jié)果進行了對比分析。

1 電磁軌道發(fā)射器的控制方程

在有限元離散方程的推導(dǎo)及數(shù)值計算過程中，做如下基本假設(shè)：

1)電樞與軌道為理想接觸，不考慮兩者之間的接觸電阻、摩擦力及材料的磨損；

2)電樞與軌道均視為剛體，忽略兩者在電磁力作用下的形變。

1.1 三維運動渦流場控制方程

對于如圖1所示的電磁軌道發(fā)射器三維求解域模型，考慮電樞運動的情況下，將求解域劃分為運動部分和靜止部分兩大子區(qū)域。其中運動部分包括電樞及同側(cè)的空氣區(qū)域；靜止部分包括軌道及同側(cè)的空氣區(qū)域。兩大子區(qū)域分別采用獨立的坐標(biāo)參考系：對于運動部分，采用運動坐標(biāo)系描述；對于靜止部分，采用靜止坐標(biāo)系描述。采用A和標(biāo)量電位φ作為變量，則整個求解域?qū)?yīng)的控制方程可統(tǒng)一寫為：

×(v×A)+σ

(1)

(2)

式中：v為磁阻率；σ為電導(dǎo)率。以上參數(shù)均與區(qū)域有關(guān)，比如σ在空氣區(qū)域為零，在導(dǎo)體區(qū)域不為零。

圖1 電磁軌道發(fā)射器的三維求解域示意圖Fig.1 Computational domains of electromagnetic rail launcher in three dimensions

對于運動區(qū)域，A對時間t的全導(dǎo)數(shù)與A對t的偏導(dǎo)數(shù)這兩者之間的關(guān)系為：

(3)

式中，v為導(dǎo)體的運動速度。顯然，對于靜止部分，A對t的全導(dǎo)數(shù)與A對t的偏導(dǎo)數(shù)相等。

引入時間積分標(biāo)量電位V，則φ與V的關(guān)系為：

(4)

在如圖1所示的求解區(qū)域中，分界面條件及邊界條件采用A與V的描述形式為：

在Γi,j上

(5)

在Si上

(6)

式中：nij表示區(qū)域Ωi與Ωj之間的分界面Γi，j的單位法向矢量；ni為邊界面Si上的單位法向矢量；下標(biāo)i，j表示求解域的編號，i，j均為整數(shù)，且滿足i

(7)

對于式(5)中的分界面條件，在分界面處采用雙節(jié)點并通過罰函數(shù)法保證上述所有分界面條件的成立。

1.2 三維溫度場控制方程

對于如圖1所示的求解域，其溫度場控制方程為：

(8)

其中，κ，c，ρ，T分別為材料的熱傳導(dǎo)系數(shù)、比熱容、質(zhì)量密度及溫升。

對于溫度場問題，在Si上的邊界條件為：

(9)

在Γi,j上

(10)

2 有限元離散方程

對于三維渦流場控制方程，采用加權(quán)余量法可對控制方程(1)和方程(2)進行離散如下：

(11)

∫Ωeσ

(12)

式中，Wi,Wi分別為節(jié)點i的矢量形狀函數(shù)和標(biāo)量形狀函數(shù)，下標(biāo)i代表單元內(nèi)節(jié)點的編號。本文研究中采用8節(jié)點六面體單元，則i為1～8的整數(shù)。

進行單元分析后，可得求解域的離散控制方程為：

(13)

式中，K，M，F(xiàn)分別為總體的剛度矩陣、質(zhì)量矩陣和右端向量，U為節(jié)點待求的物理量。

采用歐拉向后差分格式進行時間的離散，式(13)可變?yōu)椋?/p>

(14)

式中Δt為時間步長。

類似地，可對三維溫度場控制方程進行空間和時間上的離散，得到與式(14)類似的離散控制方程。

在發(fā)射過程中，電樞因受到電磁驅(qū)動力不斷向前加速，假設(shè)電樞在每一個時間步長內(nèi)做勻加速運動，則在Δt的時間內(nèi)

(15)

式中，Δx為每一時間步下的電樞位移。為保證運動部分與靜止部分的網(wǎng)格不發(fā)生交錯，運動邊界兩側(cè)的單元在運動方向的尺寸相同，且每一時間步長由式(15)及式(16)確定。

(16)

通過以上離散方法，建立了動態(tài)發(fā)射下的電磁軌道發(fā)射器的三維渦流場有限元模型；結(jié)合三維溫度場的有限元模型，最終建立了相應(yīng)的電磁-溫度耦合模型。同時，采用運動坐標(biāo)系處理電樞運動問題可保證最終形成的系數(shù)矩陣的正定對稱性，在有限元方程數(shù)值求解過程中，模型具有無條件穩(wěn)定性。

3 模型驗證與分析

為驗證所建模型的正確性，以奧卡盧薩電樞測試發(fā)射裝置(Okaloosa Armature Tester, OAT)發(fā)射裝置為研究對象[20-21]，將模型仿真計算結(jié)果與試驗編號為SLK 0018的試驗結(jié)果及EMAP3D仿真結(jié)果進行了對比分析。其中OAT發(fā)射裝置口徑為25 mm × 25 mm，試驗用的電樞型號為MCA 103，其質(zhì)量約為249.72 g。

圖2為OTA發(fā)射裝置的四分之一有限元模型。其中空氣區(qū)域的單元未畫出。發(fā)射裝置的矩形軌道尺寸為31.75 mm × 6.35 mm，軌道與電樞的材料分別為電解韌(Electrolytic Tough Pitch, ETP)銅合金和T7075 鋁合金，兩種材料電導(dǎo)率的溫度特性如圖3所示。

圖2 軌道及電樞的有限元模型Fig.2 Finite element models of rail and armature

圖3 材料電導(dǎo)率的溫度特性Fig.3 Thermal properties of material conductivities

在每一時間步長進行計算過程中，模型選擇迭代因子w=0.15時，實際的最大迭代次數(shù)則為6，具有良好的收斂性。

在仿真過程中，模型與EMAP3D軟件采用相同的輸入條件。圖4為模型計算得出的電流波形與試驗電流波形的對比圖，圖5為模型計算得出的電樞速度曲線與EMAP3D計算值及試驗值的對比圖。其中，電樞速度的試驗擬合曲線是根據(jù)試驗編號為SLK 0018中的電樞出口速度測量值及試驗電流波形，并通過選取合適的軌道電感梯度值推算出的。

由圖5可知：當(dāng)模型采用如圖4所示的脈沖電流波形時，電樞側(cè)與軌道側(cè)的實際電流計算值吻合較好。與EMAP3D計算結(jié)果相比，有限元模型的電樞速度計算結(jié)果與試驗測試結(jié)果吻合更好。

圖4 激勵電流試驗及仿真波形Fig.4 Current profiles of FEM code and measured value

圖5 速度計算值與測試值的對比Fig.5 Comparison of calculated velocity with measured data

如表1、表2所示：在試驗過程中，脈沖電流峰值為677 kA，電樞在膛內(nèi)加速1.0 ms后出口，電樞出口速度為247.0 m/s，出口電流為461.5 kA；模型計算得出的軌道側(cè)峰值電流為676.2 kA，電樞側(cè)峰值電流為669.8 kA，分別比試驗值小0.12%和0.11%，計算至1.0 ms時，模型計算得出的軌道側(cè)出口電流為461.8 kA，電樞側(cè)出口電流為465.0 kA，分別比試驗值大0.06%和0.76%，電樞末速度計算值約為257.1 m/s，比試驗值大4.09%；而EMAP3D程序計算得出的電樞末速度為268.2 m/s，比實驗值大8.58%。由此可見：所建模型的數(shù)值計算值與試驗值吻合較好，間接驗證了模型的正確性。

表1 模型的電流計算值與試驗值的對比Tab.1 Comparison of calculated currents data with measured data

表2 模型的速度計算值與試驗值對比Tab.2 Comparison of calculated velocity data with measured data and EMAP3D’s data

圖6與圖7分別為動態(tài)發(fā)射情況下，電樞運動1.0 ms后的軌道及電樞溫度云圖。由圖6可以看出：由于模型并未考慮摩擦熱的作用，僅在焦耳熱的作用下，由于軌道起始段的通流時間最長，因此，軌道的溫度最大值出現(xiàn)在起始段，具體數(shù)值約為150 ℃。由圖7可知：電樞喉部最高溫度約為953 ℃。

圖6 t=1.0 ms時，軌道溫度云圖Fig.6 Temperature contour of rail at t=1.0 ms

圖7 t=1.0 ms時，電樞溫度云圖Fig.7 Temperature contour of armature at t=1.0 ms

圖8為電樞加速運動過程中，電磁軌道發(fā)射器的電感梯度隨時間的變化圖，從圖中曲線可知：發(fā)射器的電感梯度在電流的峰值時刻附近達到最大值，此時的電感梯度大小約為0.46 mH/m，隨著時間的推移，電感梯度逐漸減小，在仿真結(jié)束時刻，發(fā)射器的電感梯度減小至0.42 mH/m。

圖8 電感梯度隨時間的變化Fig.8 Inductance gradient varies with time

4 結(jié)論

電磁-溫度耦合模型可進行電磁軌道發(fā)射器動態(tài)發(fā)射過程的數(shù)值模擬。通過模型計算值與EMAP3D程序計算值、試驗測試值的對比分析可知：模型的軌道側(cè)及電樞側(cè)電流計算值、電樞出口速度等關(guān)鍵參數(shù)均與試驗值吻合較好，間接驗證了模型的正確性。所開發(fā)的有限元程序可為電磁軌道發(fā)射器動態(tài)發(fā)射過程中的磁、熱分析及發(fā)射器的設(shè)計提供一種有效的仿真工具，為軌道燒蝕機理研究提供技術(shù)支撐，在一定程度上促進電磁軌道發(fā)射技術(shù)的發(fā)展。同時，其研究成果為電磁發(fā)射仿真軟件的國產(chǎn)化和未來軟件的并行化奠定了理論基礎(chǔ)。

References)

[1] Hsieh K T, Kim B K. 3D modeling of sliding electrical contact[J]. IEEE Transactions on Magnetics, 1997, 33(1): 237-239.

[2] 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, 2008, 27(7)：1837-1844.

[3] 王剛?cè)A, 謝龍, 王強, 等. 電磁軌道炮電磁力學(xué)分析[J].火炮發(fā)射與控制學(xué)報, 2011, 1(3): 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(3): 69-71. (in Chinese)

[4] Liu H P, 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.

[5] 李昕, 翁春生. 電磁導(dǎo)軌炮固體電樞非穩(wěn)態(tài)電磁效應(yīng)研究[J]. 南京理工大學(xué)學(xué)報:自然科學(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. (in Chinese)

[6] 李昕， 翁春生. 塊狀固體電樞非穩(wěn)態(tài)電磁效應(yīng)的三維數(shù)值模擬[J]. 彈道學(xué)報, 2009, 21(1): 103-106. LI Xin, WENG Chunsheng. Three dimensional numerical simulation of unsteady electromagnetic effect in block solid armature[J]. Journal of Ballistics, 2009, 21(1): 103-106. (in Chinese)

[7] Rodger D, Leonard P J, Eastham J F. Modeling electromagnetic rail launchers at speed using 3D finite elements[J].IEEE Transactions on Magnetics, 1991, 27(1): 314-317.

[8] Rodger D, Leonard P J. Modeling electromagnetic performance of moving rail gun launchers using finite elements[J]. IEEE Transactions on Magnetics, 1993, 29(1): 496-498.

[9] 余海濤, 邵可然, 周克定．迎風(fēng)-線性棱邊單元法[J]. 中國電機工程學(xué)報, 1996, 16(2): 92-94. YU Haitao, SHAO Keran, ZHOU Keding. Upwind-linear edge element method[J]. Proceedings of the CSEE, 1996, 16(2): 92-94. (in Chinese)

[10] Yu H T, Shao K R, Zhou K D, et al. Upwind-linear edge elements for 3D moving conductor eddy current problems[J].IEEE Transactions on Magnetics, 1996, 32(3): 760-763.

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

[12] Hopkins D A, Stefani F, Hsieh K T. Analysis of startup behavior in a “C-Shaped” armature using linked EMAP3D/DYNA3D finite element codes[J]. IEEE Transactions on Magnetics, 1999, 35(1): 59-64.

[13] Newill J F, Powell J D, Zielinski A E. Coupled finite element codes for armature design[J]. IEEE Transactions on Magnetics, 2003, 39(1): 148-152.

[14] Watt T J, Bryant M D. Cracking and dominant stresses in the throat region of C-shaped solid armatures[J]. IEEE Transactions on Magnetics, 2007, 43(1): 418-421.

[15] Watt T J, Bryant M D. Modeling assumptions for railguns[J]. IEEE Transactions on Magnetics, 2007, 43(1): 380-383.

[16] Thiagarajan V, Milfeld K, Hsieh K T. Adapting EMAP3D to parallel processing[J]. IEEE Transactions on Magnetics, 2001, 37(1): 143-146.

[17] Hsieh K T. Parallelization of EMAP3D based on element-by-element Jacobi preconditioned conjugate gradient method[J]. IEEE Transactions on Magnetics, 2003, 39(1): 139-147.

[18] Vanicek H, Satapathy S. Thermal characteristics of a laboratory electromagnetic launcher[J]. IEEE Transactions on Magnetics, 2005, 41(1): 251-255.

[19] 謝德馨, 楊仕友. 工程電磁場數(shù)值分析與綜合[M]. 北京: 機械工業(yè)出版社, 2009. XIE Dexin，YANG Shiyou. Electromagnetic field numerical analysis and synthesis of engineering[M]. Beijing : China Machine Press, 2009. (in Chinese)

[20] Price J H, Yun H D. Design and testing of integrated metal armature sabots for launch of armor penetrating projectiles from electric guns[J]. IEEE Transactions on Magnetics, 1995, 33(1): 219-224.

[21] Hsieh K T, Kim B K. International railgun modeling effort[J]. Transactions on Magnetics, 1997, 33(1): 245-248.

Numerical simulation of dynamic launching in electromagnetic rail launcher

TAN Sai, LU Junyong, ZHANG Xiao, LONG Xinlin

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

For the numerical simulation method of the dynamic launching in electromagnetic rail launcher, a discrete formulation of three dimensional transient eddy problem for electromagnetic rail launcher considering dynamic launching was established, based on the magnetic vectorAand the time integrated electrical scalar potentialV, combined by using nodal element method and the moving coordinate frame to deal with the motion. The corresponding three dimensional magnetic-thermal coupled model was formed finally by using the governing equation for thermal field. The dynamic launching process of an electromagnetic rail launcher was simulated by the finite element model. Numerical analysis results of FEM (finite element method) code and the experiment results were compared and analyzed. The characteristic of thermal field distribution for launcher and the variation of inductance gradient with time were obtained by the model. Results show that the magnitudes of the pulsed current, muzzle current, and muzzle velocity have good correspondence with measure dates respectively, which has verified the correctness of the FEM code.

dynamic launching; armature; electromagnetic rail launcher; inductance gradient; finite element method

10.11887/j.cn.201606007

2016-03-29

國家自然科學(xué)基金資助項目(51207162，51522706)；國家部委基金資助項目(613262)

譚賽(1988—)，男，河南永城人，博士研究生，E-mail：tansai208@126.com； 魯軍勇(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：ljy041@sohu.com

TM315

A

1001-2486(2016)06-043-06

http://journal.nudt.edu.cn