楊玉東，付成芳，薛 文，趙環(huán)宇

（1.淮陰工學院 電子與電氣工程學院，江蘇 淮安 223003；2.南京理工大學 電子工程與光電技術學院，江蘇 南京 210094）

在固體電樞電磁發(fā)射器發(fā)射結果檢測中發(fā)現(xiàn)軌道的內(nèi)側及電樞后尾部經(jīng)常會出現(xiàn)燒蝕現(xiàn)象，經(jīng)研究表明，這種現(xiàn)象是由于電樞高速運動時，軌道和電樞中電流密度分布不均勻造成的，電流密度聚集在軌道和電樞的尾部，使得局部區(qū)域的電流密度過大，該局部溫度急劇上升，造成軌道材料的燒蝕［1－4］。描述軌道和電樞相對運動的電磁方程在數(shù)學上屬于對流擴散方程，目前有多種方法求解對流擴散方程，如有限差分法、有限元法和邊界元法等，但每種算法均容易出現(xiàn)數(shù)值擴散或解的振蕩解等問題，尤其是對流項占優(yōu)情況下。如文獻［5］采用有限差分法對對流項占優(yōu)的對流擴散方程求解，離散過程中對離散時間步長采用細化，用于提高速度項，實質上是采用弱化對流項方式進行求解，一旦時間過長或速度項過高，如速度超過800m／s時就會出現(xiàn)擴散及振蕩現(xiàn)象，因此如何恰當處理對流項成為解對流擴散方程的關鍵。有限元或迎風差分格式雖然能在一定程度上消除解的振蕩，但數(shù)值擴散現(xiàn)象很嚴重，即使采取更高的精度格式也難于解決此類問題［6］。目前一種求解對流擴散反應型方程的算子分裂格式DF－AD－REA （擴散－對流－反應）被廣泛應用于對流擴散方程的求解。在規(guī)則網(wǎng)格上，處理對流擴散算子和反應項時，可以分別采用不同的各自最適應其算子的數(shù)學物理特征的格式［7］。筆者采用算子分裂格式處理軌道和電樞間運動電磁場方程，計算出在激勵電流作用下運動導體內(nèi)部電磁場參數(shù)的分布，揭示電磁發(fā)射過程中出現(xiàn)的速度趨膚效應現(xiàn)象，為選擇合適的軌道、電樞材料和形狀奠定一定的理論基礎。

1 控制方程

固體電樞與軌道之間聯(lián)系是電接觸式，電樞固定，電樞與軌道間滑動接觸，激勵源從軌道的兩端加入，其結構及其二維模型如圖1所示。

由Maxwell方程組可推導出軌道和滑動電樞內(nèi)部電磁場分布方程：

式中：B為導體中磁感應強度；μ為導體的磁導率；σ為導體的電導率；t為磁場擴散時間；vx為沿著x軸方向的速度。在二維坐標中，式（1）可寫成如下偏微分方程形式：

2 計算方法

目前算子分裂法用于求解對流擴散方程已成為有效的方法之一，其優(yōu)點是算子分裂后的方程簡單、容易求解且格式靈活，穩(wěn)定性好。算子分裂方法主要有：標準的Lie分裂、Strang分裂和Sourec分裂［8］，筆者將采用Strang分裂算法，把對流項和擴散項分別分解為兩個算子。

2.1 算子分裂格式

對于式（2）采用Strang算子分裂.在每一個時間步［tn，tn＋1］采用如下A、B分裂格式：

式中，A 格式為擴散項表達式，B格式為對流項表達式。式（3）中第二個表達式是A 格式的邊界值和初始值，其含義是A 格式中的初始值和邊界值，采用式（2）的初始值和邊界值；式（4）中第二個表達式是B格式的邊界值和初始值，其采用的值為A 格式計算后得到的值。A、B 兩格式順序計算，B 格式計算的結果就是一個時間步結束后，最終得到的結果。即B（x，y，tn＋1）≈B（2）（x，y，tn＋1）。

2.2 激勵源設置

激勵源采用8 個單模塊脈沖成形網(wǎng)絡同時觸發(fā)，電路參數(shù)為：軌道長度2m，電樞質量為45g，平均電感梯度為0.46μH／m，每個模塊中儲能電容容量為15mF，調波電感為60μH，電容充電電壓為4 kV，總能量為12kJ，驅動電流波形如圖2所示。

激勵電流密度可用如下函數(shù)近似描述：

式中：ω＝，S為軌道x軸方向橫截面面積。

經(jīng)仿真計算可得到電樞運行速度與時間的關系如圖3所示。

2.3 初始條件和邊界條件

設軌道和電樞內(nèi)部的初始磁感應強度和電流密度為0；以0時刻為計算時間，此時在軌道兩端施加激勵電流，當時間步長k＝1時，電樞的左邊界磁感應強度為兩軌道產(chǎn)生的磁感應強度的疊加，B＝2μj；軌道的左下邊界、電樞左下側邊界為B＝μj；右側下邊界B＝0，在后續(xù)的計算過程中，電樞左側邊界條件一直保持為B＝2μj，右側邊界為B＝0；上下邊界條件為?B／?y＝0；軌道邊界條件左側?B／?x＝0；右側和上邊界B＝0［9－10］。

2.4 中間條件

考慮到A、B計算格式的順序性，先計算電磁場的擴散運動，再計算電磁場的對流運動。A 格式反映的是軌道和電樞中外部激勵電流在導體內(nèi)部引起的磁感應強度隨著時間由外至內(nèi)逐步擴散的過程；B格式反映的是電樞運動過程中切割磁力線引起的感應電流在導體內(nèi)部產(chǎn)生的磁感應強度的變化，其中間量是有物理意義的。電樞運動產(chǎn)生的在導體內(nèi)部的感生電流屬于渦流，渦流的方向與外加電流方向相反，實際測量的軌道激勵電流中已經(jīng)包含了感生電流，因此，B格式中采用的邊界條件應與A 格式中邊界條件相同，B 格式計算的初始條件采用A 格式計算結果。

2.5 格式計算方法

為了與文獻［5］相比較，筆者采取與文獻［5］一致的格式剖分，即考慮到軌道和電樞外形是對稱的，可選取軌道上半部分和電樞的上半部分作為計算區(qū)域，一半導軌被剖分為5 000×25個網(wǎng)格，一半電樞被剖分為25×25個網(wǎng)格［5］。為編程方便，設電樞相對于軌道靜止，軌道相對于電樞作反向運動。因此用式（3）、（4）分裂格式計算時，軌道計算取速度為－vx；電樞計算取速度vx＝0。即計算軌道磁感應強度分布采用A、B格式，而計算電樞中磁感應強度僅使用A 格式即可。

式（3）方程為二維拋物型方程，選擇P－R格式進行離散。P－R 差分格式結構較為簡單，無條件穩(wěn)定，且具有二階精度［11］。設離散格式x、y方向的空間步長為h，時間步長為τ，進一步推導可得式（3）的具體差分格式為

式（4）是一維一階雙曲型方程，選擇如下Euler差分格式進行離散：

此差分格式具有二階精度［11］。

當求得某時刻電樞和軌道內(nèi)磁感應強度B后，根據(jù)▽×B＝μJ，可相應地求出電流密度J在x、y軸方向的分量。

則J＝

3 計算結果及分析

式（3）～式（4）的初始條件確定后，采用Turbo C 2.0編程計算得到0－3.5ms內(nèi)不同時刻的軌道及電樞內(nèi)部磁感強度和電流密度分布數(shù)據(jù)，部分數(shù)據(jù)經(jīng)過Surfer處理后顯示結果如圖4和圖5所示。1.02ms時刻，導軌、電樞內(nèi)磁感強度最大值為0.7 T。此時電樞的運動速度為600m／s。最大磁感強度分布在軌道的下邊緣和電樞左邊緣，并由軌道下邊緣向上、電樞左邊緣向右遞減擴散，在電樞的尾部集中；在2.78ms時刻，導軌、電樞內(nèi)磁感強度最大值為1.2T。此時電樞的運動速度為972 m／s。隨著電樞運行速度增加，磁感應強度在軌道內(nèi)部向外面擴散的深度減小，而在電樞內(nèi)部，向右深度擴散增強，并有擴散至電樞整個區(qū)域的趨勢。

此計算結果、電磁場分布云圖與參考文獻［5］呈一致性，區(qū)別在于文獻［5］在計算過程中，當速度項超過600m／s時，計算結果容易發(fā)散或出現(xiàn)振蕩現(xiàn)象，即使對時間步長進一步減小到10－6，也只能使得速度限制在800m／s左右，計算過程耗時長且結果不穩(wěn)定。而采用算子分裂格式時，當計算時長達到1.5ms以上，即時速度達到700m／s，時間步長僅為10－5時，計算數(shù)據(jù)依然穩(wěn)定收斂。在不減小時間步長情況下，當計算至2.78 ms、即時速度為972m／s時，仍能得到穩(wěn)定的計算結果。事實上式（3）為拋物型方程，其中不含速度項，采用P－R 格式是采用交替隱形格式進行求解，無條件穩(wěn)定，具有二階精度；不穩(wěn)定的因素存在于式（4）中，式（4）為一維一階雙曲型方程，其穩(wěn)定和收斂的條件為

式中λ＝τ／h。則根據(jù)式（10）的條件，設時間步長為10－5，空間步長為10－2，則最大收斂速度可達到1km／s。欲提高速度項，可通過減小時間步長或增大空間步長來獲得，如空間步長不變，時間步長減小至10－6，則速度可提高至10km／s。這是理論上分析的結果，實際上電樞運動的速度并不是線性增加的，當電樞運行至激勵電流下降期，其加速度增加緩慢，速度變化需要更多的時間，反映在程序運行上，迭代的次數(shù)將成倍增加，嚴重耗費機時和內(nèi)存。

圖5為瞬時電流密度在軌道和電樞內(nèi)部的分布云圖。1.02ms時刻軌道和電樞內(nèi)的電流密度最大值約為8.8×108A／m2；2.78 ms時約為2.3×108A／m2，與圖4對比可知，軌道和電樞內(nèi)電流密度與磁感應強度分布規(guī)律是一致的，電流密度由邊緣至導體內(nèi)部擴散呈遞減趨勢，在軌道和電樞拐角處最大；隨著速度逐漸增加，電流密度越來越趨近于軌道和電樞的內(nèi)側。這是因為：在靜止情況時，導體中電流密度在短時間內(nèi)會迅速擴散；而當電樞處于高速運動狀態(tài)時，在軌道和電樞接觸面的電流密度會被強行牽引至其尾部，使得電樞和軌道的邊緣處，尤其電樞與軌道相接觸局部區(qū)域電流密度最大［12］。尤其需要注意的是，當軌道的內(nèi)表面及電樞的尾部的電流密度過于集中時，會引起接觸區(qū)局部區(qū)域導體產(chǎn)生過高的歐姆熱，從而引起局部區(qū)域溫度短時間內(nèi)迅速升高，當溫度達到電樞及軌道材料的熔點時，會引起該區(qū)域的導體燒蝕，上述計算結果與文 獻［5］、［9］、［10］、［12］結論和數(shù)據(jù)基本吻合。

4 結論

在高速運動的固體電磁發(fā)射過程中出現(xiàn)的速度趨膚效應是運動導體在強磁場中呈現(xiàn)的一種現(xiàn)象，是電磁場和熱傳遞共同作用的結果，尤其在大電流、高速度情況下顯得尤其突出。其結果會造成電樞和導軌局部區(qū)域出現(xiàn)燒蝕，使得難以進行多次發(fā)射，會大大降低電磁發(fā)射效率。建立了軌道和電樞內(nèi)部電磁場分布的二維對流擴散方程，并采用算子分裂格式將方程分解為擴散方程和對流方程。因為擴散方程是拋物型方程，對其求解有穩(wěn)定、精度較高的成熟算法，對流方程為一階雙曲方程，對其求解也比較簡單。求解過程避免了對流項的存在造成方程數(shù)值振蕩和擴散問題。數(shù)值計算過程表明，當速度達到972 m／s時，此時對流項絕對占優(yōu)，計算數(shù)據(jù)仍然穩(wěn)定收斂，避免了文獻［5］中當電樞運動速度高時導致數(shù)據(jù)發(fā)散的現(xiàn)象。計算過程存在以下不足：模型中假設電樞與軌道間接觸良好，不存在接觸電阻，而在實際試驗中發(fā)現(xiàn)兩者的接觸電阻是存在的，且是造成熱損耗的一個重要因素；實際的電磁發(fā)射過程中，需要考慮運動、電磁場、熱場之間的耦合關系，而這種耦合關系是極為復雜的，筆者對這種耦合關系處理較為簡單，不可避免帶來誤差。還有其他若干不足地方不能一一列舉，這將在今后的研究工作中加以逐步改進。

（References）

［1］MARSHALL R A.The distributed store railgun，its efficiency，and its energy store implications［J］.IEEE Transaction on Magnetics，1997，23（1）：582－583.

［2］KAMRAN D，RAHIMZADEH M.Dynamic response and critical velocity studies in an electromagnetic railgun［J］.IEEE Transaction on Magnetics，2007，43（1）：126－127.

［3］THOMAS G，ENGEL Jese M.Efficiency and scaling of constant inductance gradient DC electromagnetic launchers［J］.IEEE Transaction on Magnetics，2006，42（8）：2044－2047.

［4］DAVID A H.Analysis of startup behavior in a Cshaped armature using linked EMAP3D LDYNA3D finite element codes［J］.IEEE Transaction on Magnetics，1999，35（1）：60－62.

［5］楊玉東，王建新，薛文.軌道炮速度趨膚效應的分析與仿真［J］.強激光與粒子束，2011，23（7）：1965－1968.YANG Yudong，WANG Jianxin，XUE Wen.Simulation and analysis for velocity skin effect of railgun［J］.High Power Laser and Particle Beams.2011，23（7）：1965－1968.（in Chinese）

［6］POWELL J D，ZIELINSKI A E.Ohmic heating in a double－taper sabot－armature［J］.IEEE Transaction on Magnetics，2003，39（1）：153－157.

［7］賈宏恩，李開泰，鐘賀.算子分裂法求解對流－擴散－反應方程［J］.工程數(shù)學學報，2012，29（1）：89－91.JIA Hongen，LI Kaitai，ZHONG He.An operator splitting scheme for advection－diffusion－reaction equation［J］.Chinese Journal of Engineering Mathematics，2012，29（1）：89－91.（in Chinese）

［8］曹志先，魏良琰.對流擴散反應型方程的一種穩(wěn)定的算子分裂格式［J］.水動力學研究與進展A 輯，1991，6（1）：60－62.CAO Zhixian，WEI Liangyan.An operator splitting scheme with stability for advection－diffusion－reaction equations［J］.Journal of Hydrodynamics，1991，6（1）：60－62.（in Chinese）

［9］李昕，翁春生.固體電樞電磁導軌炮非穩(wěn)態(tài)電磁效應［J］.南京理工大學學報，2009，2（1）：108－109.LI Xin，WENG Chunshen.Unsteady electromagnetic effection solid armature railguns［J］.Journal of Nanjing University of Science and Technology，2009，2（1）：108－109.（in Chinese）

［10］李昕，翁春生.U 型電樞非穩(wěn)態(tài)電磁場二維數(shù)值模擬［J］.火炮發(fā)射與控制學報，2009，（1）：1－3.LI Xin，WENG Chunshen.Two dimension numericalsimulation of unsteady electromagnetic in U shaped solid armature［J］.Journal of Gun Launch＆Control，2009，（1）：1－3.（in Chinese）

［11］孫志忠.偏微分方程數(shù)值解法［M］.北京：科學出版社，2005：120－137.SUN Zhizhong.Numer methods partial differential eq［M］.Beijing：Science Press，2005：120－137.（in Chinese）

［12］POWELL J D，ZIELINSKI A E.Observation and simulation of solid－armature railgun performance［J］.IEEE Transaction on Magnetics，1999，35（1）：84－89.