面向電磁發(fā)射的多場耦合仿真分析方法

周東強，嚴 藝，楊 易，高振良，劉奕鑫

（北京空間飛行器總體設計部，北京100094）

0 引言

電磁發(fā)射是一種全新概念的發(fā)射方式。電磁發(fā)射技術將電磁能變換為發(fā)射物體需要的瞬態(tài)動能，可在短距離內將發(fā)射物體加速至高速，能有效突破傳統(tǒng)發(fā)射的速度和能量極限，在軍事和民用領域具有巨大的優(yōu)勢和廣闊的應用前景[1-2]。電磁發(fā)射過程中對發(fā)射裝置加載脈沖大電流，伴隨著磁場擴散、熱傳導、結構變形等物理過程，是具有瞬態(tài)特性的多物理場耦合問題[3]。

電磁發(fā)射的數值仿真是發(fā)射裝置設計、試驗的重要支撐手段。國內外對電磁發(fā)射的數值仿真與計算都非常重視，開展了大量研究工作，使用多種數值計算方法模擬電磁發(fā)射的多物理場相互作用。Powell 和Zielinski[4]提出了計算電磁場、電流密度、導軌和電樞溫度等隨空間、時間分布的數學模型。Ghassemi 和Pasandeh[5]考慮導軌?電樞之間的大電流產生的熱能對導軌?電樞組件的電學、熱學、力學性能影響，在麥克斯韋方程組中引入能量方程，同時考慮導軌和電樞中的電傳導、歐姆熱和摩擦熱，得到非線性控制微分方程組，計算導軌和電樞中的熱和電磁感應分布。隨著有關數值仿真應用需求的增長，研究人員逐漸開發(fā)出用途各異的軟件，如美國的EMAP3D，法國和德國的EMAS，英國的MEGA 等[6-11]，一些大型的商業(yè)軟件也可以進行電磁軌道發(fā)射裝置仿真計算。然而，國外部分自用仿真軟件不對外公開，且其他大型商業(yè)軟件針對軌道電磁發(fā)射在電樞高速運動多物理場耦合計算方面存在局限性，在處理高速滑動電接觸問題上能力不足，無法計算特有的速度趨膚效應等[3,11-15]，因此可以說，對固體電樞電磁發(fā)射過程的建模和仿真尚處在不斷探索與研究之中[16-20]。

本文基于有限元法的三維電磁、熱、力瞬態(tài)耦合算法，根據時間步進順序耦合方法構建電磁、熱、力多物理場耦合集成框架，用以指導相關計算程序的編制，實現對電磁發(fā)射中電樞高速運動狀態(tài)下電磁、熱、力瞬態(tài)耦合計算以及電樞高速滑動等多物理場瞬態(tài)狀態(tài)的定量分析。

1 瞬態(tài)電磁分析與有限元模型

軌道電磁發(fā)射是瞬態(tài)電磁相互作用的過程，伴隨著導軌與電樞間的電流傳導、電磁感應以及趨膚效應、渦流等現象。麥克斯韋方程是表征所有宏觀電磁場現象的基本方程[3,5,18]，電磁分析的實質是求解給定邊界條件下的麥克斯韋方程組問題。

電磁發(fā)射的三維瞬態(tài)電磁計算基于麥克斯韋理論方程和體現發(fā)射裝置材料特性的本構關系，采用伽遼金方法[21]對計算模型進行離散，使用三維高斯積分方法進行數值積分，通過計算機編程實現瞬態(tài)電磁、電位的計算。電磁和電位瞬態(tài)分布的仿真計算是電磁發(fā)射中多物理場耦合分析的重要環(huán)節(jié)。

1.1 麥克斯韋方程微分形式與電磁分析本構關系

三維瞬態(tài)電磁場分析是以麥克斯韋方程組為基礎，根據設定的邊界條件求解場域內電場、磁場、電流密度、渦流等狀態(tài)量隨時間和空間變化的過程[6]。麥克斯韋方程組由4個基本定律組成，分別為安培環(huán)路定律、法拉第電磁感應定律、高斯電通定律和高斯磁通定律。為了方便使用有限元數值算法，選用其微分形式建立三維瞬態(tài)電磁分析的數學模型：

式(1)～式(4)中：H為磁場強度；J為電流密度；D為電位移；E為電場強度；B為磁感應強度；ρ為電荷密度。

不同材質的發(fā)射裝置對發(fā)射中多物理量的作用結果不同，為了反映電磁發(fā)射裝置的材質構成對電場、磁場、感應電流的相互作用的影響，在電磁發(fā)射的數值計算中使用反映材質宏觀性質的本構關系來計算介電常數、磁導率、電導率對電磁場的影響和作用[6]：

式(5)～式(7)中：ε為介電常數、μ為磁導率、σ為電導率，各向同性材料用標量表示，各向異性材料用張量表示；v為電樞發(fā)射速度。

為了求解麥克斯韋方程，通過定義標量電勢Φ與矢量磁勢A，將包含電場與磁場2個場量的一階微分方程轉化為2個獨立的包含1個場量的二階微分方程[5]。將電位標量勢E= -?Φ與本構關系D=εE代入高斯電通定律（式(3)）中推導可得電位標量勢方程

磁位矢量勢滿足高斯磁通定律，由本構關系B=μH知μH= ?×A，代入安培環(huán)路定律（式(1)）中推導可得磁位矢量勢方程

1.2 電位標量與磁位矢量有限元列式

為了使用計算機編程實現電磁發(fā)射中電磁、電位的數值積分計算，本文使用伽遼金方法建立電位標量勢方程的三維有限元格式，構造8節(jié)點六面體等參單元，其形函數N為三維空間參考坐標(ξ,η,ζ)的函數，

2 熱傳導分析與有限元模型

電磁發(fā)射過程中產生熱的途徑主要有電樞與導軌摩擦生熱、脈沖激勵電流生熱及感生渦流生熱等；熱擴散途徑主要為熱傳導和熱對流。本文主要考慮瞬態(tài)熱傳導問題的基本方程與有限元方法在熱傳導中的應用。

2.1 熱傳導與熱對流方程

2.2 瞬態(tài)熱傳導有限元列式

與電磁計算相似，將計算域離散成有限個8節(jié)點六面體單元，單元的溫度場可以用節(jié)點溫度Ti插值得到，即

在瞬態(tài)熱傳導計算中此溫度場Ti是與時間t有關的函數，Ni為8節(jié)點六面體形函數。使用伽遼金積分法建立三維瞬態(tài)熱傳導的等效積分形式[5]

將節(jié)點溫度插值代入熱傳導方程（式(15)），得到有限元方程

由此，熱傳導方程可以轉換為以時間t為獨立變量的線性常微分方程組。式(19)中：C為熱容矩陣；K為熱傳導矩陣；P為溫度載荷列陣；T˙為節(jié)點溫度對時間偏倒列陣。

3 瞬態(tài)多場耦合計算程序

計算程序對時間單元進行離散，以電磁、熱、力收斂為準則進行編制與開發(fā)，構建瞬態(tài)多場耦合計算程序框架如圖1所示。計算數據信息主要包括控制信息與有限元前處理信息：控制信息主要為物理計算時間與計算時間步長以及輸出控制；有限元前處理信息主要包含節(jié)點信息、單元拓撲關系、材料信息和邊界條件等。

圖1 瞬態(tài)多場耦合計算程序框架Fig.1 Framework of the coupled transient multi-field calculation program

計算程序中的各子模塊功能分別為：

1）電場分析子模塊：電流曲線插值，組裝電場剛度矩陣；

2）磁場分析子模塊：組裝磁場剛度矩陣、磁場遠場邊界；

3）熱傳導分析子模塊：組裝熱傳導剛度矩陣、主動冷卻與被動冷卻邊界條件；

4）動態(tài)響應子模塊：結構動態(tài)響應計算，電樞速度與位移計算；

5）電接觸子模塊：電樞與導軌接觸面準匹配網格搜索，瞬變狀態(tài)量插值。

軌道電磁炮發(fā)射激勵源是由外在電流加載到2根導軌端部形成回路，電樞與導軌接觸是電場分布計算的先決條件；當電樞發(fā)射脫離導軌后回路中斷，外在激勵源不再維持回路電磁感應，計算時僅需考慮熱場傳導和結構動態(tài)效應，不再考慮電磁感應，因此電樞出膛后不再調用電場和磁場分析子模塊。

使用FORTRAN90搭建程序主體框架，各子模塊與關鍵部分均采用子程序模式，主體框架與子程序間為調用與被調用關系，主體框架代碼如下：

4 典型軌道電磁發(fā)射計算與分析

基于已完成的多物理場耦合有限元計算程序，以典型軌道電磁炮發(fā)射過程為例，分析計算電流密度、電磁場、熱場以及發(fā)射出口速度等多物理場的狀態(tài)量分布情況。發(fā)射裝置主要由導軌、電樞、聚碳酸酯絕緣體、約束槽等組成，為了同時保證計算精度、速度與效率，對計算模型進行合理簡化。

4.1 軌道電磁發(fā)射模型

軌道電磁發(fā)射仿真模型包含有限電磁傳播域網格、導軌網格、電樞網格和電流激勵片。如圖2所示，電磁傳播域網格單元、導軌網格單元、電樞網格單元為8 節(jié)點六面體體單元，電流激勵片網格單元為四邊形面單元。金屬導軌尺寸為15 mm×20 mm×500 mm，磁場域尺寸為200 mm×200 mm×1000 mm。

圖2 軌道電磁發(fā)射計算模型Fig.2 Model for calculating orbit electromagnetic emission

根據電磁發(fā)射裝置設計通常使用的材質情況，在數值計算中對各計算區(qū)域進行材料參數賦值。如表1所示：導軌采用銅制材料；電樞采用鋁制材料；導軌與電樞以外電磁域媒質選用真空媒質材料參數。

表1 計算所用的相關材料參數Table 1 Material parameters for the electromagnetic calculation

4.2 載荷與邊界條件

根據軌道電磁炮電磁、熱、結構計算的載荷與邊界情況，定義激勵載荷1個（脈沖電流荷載），邊界條件3個（無限電磁域、熱對流邊界、導軌位移約束）。

電流加載方式采用瞬態(tài)脈沖式電流，在模擬電磁發(fā)射過程時選擇歷時3 ms的單脈沖式電流加載。選擇導軌與電樞形成的導體回路的2個端頭分別為電流輸入端和電流輸出端（參見圖3）；在輸入端與輸出端建立加載單元，通過導軌輸入端頭與輸出端頭加載單元將激勵電流分配到導軌加載單元節(jié)點上（參見圖4），使電流均勻加載到激勵單元，按照激勵單元各節(jié)點電流貢獻值分配電流承載值。

圖3 軌道電流加載方向示意Fig.3 Schematic diagram of current loading direction

圖4 電流加載單元與加載節(jié)點示意Fig.4 The current loading unit and theengaging nodes

4.3 多場耦合計算結果與分析

4.3.1 脈沖電流對電樞速度的影響

為了比對脈沖電流峰值對電樞發(fā)射速度、加速度的影響規(guī)律，分別對電流脈沖峰值為3.0 kA、3.2 kA、3.4 kA、3.6 kA 共4種激勵工況進行計算，脈沖加載歷程皆為3 ms，峰值狀態(tài)均持續(xù)0.5 ms。

圖5～圖7分別為電樞的發(fā)射速度曲線、發(fā)射位移曲線和發(fā)射加速度曲線，各圖中曲線C1～C4分別代表電流脈沖峰值為3.0 kA、3.2 kA、3.4 kA、3.6 kA 的4種激勵工況?？梢钥闯觯簭拈_始加載電流到0.6 ms 時刻期間，電樞發(fā)射速度明顯提升，0.6 ms之后速度增長趨緩；曲線C4與曲線C1相比，出口速度提高68.61 m/s，電樞脫離導軌時間縮短0.6 ms；電樞加速度變化與加載電流變化高度相關，電樞脫離導軌時電樞與導軌間接觸有效電阻發(fā)生變化，加速度會出現微小波動。

圖5 計算出的電樞發(fā)射速度曲線Fig.5 Curve of calculated launching speed of the armature

圖6 計算出的電樞位移曲線Fig.6 Curve of calculated launching distance of the armature

圖7 電樞加速度計算曲線Fig.7 Curve of calculated acceleration of the armature

4.3.2 導軌與電樞中電流密度分布

分析電磁發(fā)射全過程導軌和電樞中的電流密度分布情況：電樞在導軌中滑行時，導軌的最高電流密度約為平均電流密度的1.39～1.42倍，電樞的最高電流密度約為平均電流密度的1.53～1.57倍；當電樞即將脫離導軌時，電樞與導軌的接觸面較小，電流集中明顯增強。

圖8顯示了0.4 ms時刻電流密度集中區(qū)域，其中電流密度最大區(qū)域在電樞尾部與導軌接觸位置，體現了電樞高速滑動尾部趨膚效應。為了觀察電流集中區(qū)內電流密度的時間歷程曲線，并且能夠與非電流集中區(qū)內電流密度的時間歷程曲線對比，選取7個檢測點作為電流密度時間歷程取值點（見圖9）：node1為導軌網格單元附屬節(jié)點，處在導軌外載電流加載端頭內側；node2、node3、node5、node6位于電樞與導軌接觸面；node4與node7位于電樞尾部與頭部對稱位置。

圖8 計算出的0.4 ms時刻電流密度分布云圖Fig.8 Cloud map of calculated current density after 0.4 ms of current loading

圖9 電流密度時間歷程曲線取值點Fig.9 The location pointsfor investigating the timeprofile of current density

圖10為導軌各取值點電流密度的時間歷程曲線，可以看到，node2＞node3＞node1＞node5＞node6。表明：導軌上電樞尾角與導軌接觸位置電流密度最大；在電樞與導軌接觸面上沿著電樞發(fā)射方向電流密度逐漸變小，電樞頭部電流密度小于外載電流加載端電流密度。

圖10 導軌測試點電流密度時間歷程曲線Fig.10 Time history of current density at typical positions of the rails

圖11為電樞各取值點電流密度的時間歷程曲線，可以看到，node2＞node3＞node4＞node5＞node7＞node6。表明：電樞尾角處電流密度最大；電樞與導軌接觸面電流密度大于電樞體非接觸面的電流密度；電樞頭部電流密度小于電樞尾部電流密度。

圖11 電樞測試點電流密度時間歷程曲線Fig.11 Time history of current density at typical positions of the armature

4.3.3 空間計算域電磁場分布

圖12為導軌與電樞的電磁力計算結果矢量分布，可以看到，2根導軌主要承受y方向排斥力，并且電磁力主要集中在電流回路，導軌內側電磁力大于外側電磁力。截取磁場域z方向截面磁勢分布云圖（圖13）可以看到，兩導軌之間磁勢大約為導軌兩側磁勢的2倍，電樞x方向尾部磁勢大約為頭部磁勢的2倍。磁勢最大位置在電樞尾部-x方向，最大磁勢區(qū)域與電樞所在區(qū)域不存在交集，鑒于此可以改善電樞結構，優(yōu)化最大磁勢分布區(qū)域。

圖12 導軌與電樞電磁力計算結果矢量分布Fig.12 Distributions of calculated electromagnetic force vectorsof the railsand the armature

圖13 z 方向截面磁勢分布云圖Fig.13 Cloud diagram of magnetic potential on cross section

4.3.4 導軌與電樞中熱場分布

導軌與電樞是發(fā)熱最嚴重的部件，其熱源主要為外在電流生熱，熱源體積密度為q=J2/σ。電磁發(fā)射過程中電流生熱會在導體中累積，且由于電流傳遞不均勻，導軌與電樞中的熱分布梯度非常明顯。圖14 為不同時刻導軌與電樞的熱密度分布云圖。從圖中可以看到：0.3 ms時刻熱密度分布與電流密度分布情況相近，熱密度較高區(qū)域與電流密度集中區(qū)域相似，在導軌與電樞接觸面上熱密度較大，電樞尾部熱密度比電樞頭部熱密度大，最大值為1.5×1012J/m3；隨后逐步出現熱密度累積現象，在0.6 ms時刻電流密度最大的位置為與導軌接觸面的電樞尾部，最大熱密度為2.4×1012J/m3。圖14(c)與圖14(d)分別代表電樞離開初始位置一段距離后的熱密度分布情況，最大值分別為3.0×1012J/m3和2.8×1012J/m3?？梢钥闯?，熱密度集中區(qū)域仍然為電樞在初始位置尾部的導軌處，2.0 ms時刻導軌電流加載端頭整個區(qū)域都成為熱密度集中區(qū)。也就是說，在重復發(fā)射情況下，電樞初始位置將會成為熱燒蝕集中發(fā)生部位。

圖14 不同時刻的熱密度分布云圖Fig.14 Cloud diagram of heat density at different moments

5 結束語

本文基于時間步進耦合算法和瞬態(tài)電磁、熱傳導的有限元算法，研究了適用于電磁發(fā)射的多物理場耦合計算程序設計方法，對電磁軌道炮進行三維瞬態(tài)模擬計算，得到磁感應密度、電流密度、熱密度的時域和空間域分布。計算結果表明：時間步進耦合算法可以較方便地實現電磁發(fā)射中多物理場瞬態(tài)耦合仿真計算；電磁發(fā)射過程中電流密度分布呈動態(tài)擴散特點，在電樞與導軌接觸面存在明顯的速度趨膚效應，最大電流密度分布在電樞尾部與導軌接觸位置；導軌中電樞發(fā)射的初始位置存在熱密度累積現象，該位置也是電磁發(fā)射短時間內連續(xù)使用較容易發(fā)生熱燒蝕和熱損傷之處。

以上研究有助于對電磁發(fā)射中的多物理場耦合演化過程的了解，可為發(fā)射裝置設計提供輔助支撐。但計算程序中暫未考慮接觸摩擦磨損等因素，須對計算模型進一步優(yōu)化和完善。