基于瞬態多物理場求解器的電磁軌道炮發射過程建模與仿真

林慶華，栗保明

(南京理工大學 瞬態物理國家重點實驗室，江蘇 南京 210094)

0 引言

電磁發射是近年比較活躍的一個前沿技術領域，電磁軌道炮、電磁線圈炮等研究方向受到廣泛關注，并已經在一些關鍵技術上取得突破。發射過程的建模與數值仿真是支撐電磁炮技術發展的重要基礎之一，它利用對發射過程物理圖像及其規律的理論描述，為工程設計和試驗提供參考。

大電流高速滑動電接觸是電磁軌道炮發射工況的典型特征[1-2]。國內外為了研究高速滑動電接觸及其所導致的電磁學、熱學和力學效應，已經或正在開發一些專用計算程序，包括EMAP3D[3]、MEGA[4]、HERB[5]、EFEM3D[6]、RGUN3D[7]等。這些程序的核心功能是求解帶有高速滑動電接觸的電磁場問題，它們普遍采用了有限元方法，但在具體算法上有些差別，例如程序EMAP3D使用了節點元，程序EFEM3D使用的是棱邊元。在電磁場計算的基礎上，Hopkins等[8]將EMAP3D程序與著名的顯式動力學計算程序DYNA3D耦合，用于求解電磁軌道炮在脈沖電磁力作用下的結構動態響應。Leyden等[9]將程序MEGA計算出的溫度載荷傳遞給程序DYNA3D，實現了對電樞熱軟化過程的模擬。

電磁線圈炮中沒有滑動電接觸，其電磁場問題相對容易解決，然而一旦涉及到電樞在局域化感應渦流及電磁力作用下的發熱和變形，以及電樞和導向管間的高速摩擦、碰撞，也會遇到與電磁軌道炮相似的非線性、多物理場耦合等問題。

以往的電磁炮建模與數值仿真以研究發射機理為主要目的，因此模型中只需考慮電樞- 軌道或者電樞- 線圈系統[10-11]。隨著電磁炮技術的發展，數值仿真工作逐漸涉及到先進發射器結構和一體化發射單元(ILP)[12-14]，特別是對于攜帶含能材料和電子部件的ILP，其發射安全性的研究和評估工作需要相關模型和數值仿真程序的支撐，這對電磁炮的建模和數值仿真提出了新的要求。

南京理工大學瞬態物理國家重點實驗室開展電磁發射理論與技術研究已有十余年，逐步建立起涵蓋電磁場、熱場和結構場的瞬態多物理場模型，編寫了電磁炮的求解器程序。本文簡要介紹了該求解器的組成框架、基本功能、數學模型和計算方法，通過電磁軌道炮和電磁線圈炮的具體算例，展示了該求解器對電磁發射過程的仿真計算能力。

1 求解器框架與功能介紹

求解器的組成框架如圖1所示。它針對的是包含脈沖電源、發射器和ILP在內的整個發射系統，具有電路、電磁場、熱場、結構場4個功能模塊，其中脈沖電源的放電過程用電路模型描述，發射器和ILP的工作過程用場模型描述。

圖1 求解器框架Fig.1 Framework of solver

電磁場是整個求解器的核心，它以脈沖電流作為激勵條件，在求解域的導體內傳導或感應出電流，產生焦耳熱和電磁力的分布。激勵電流可以通過數據曲線的形式直接輸入，也可以通過場路耦合由電路模型實時計算。進行場路耦合計算時，發射器和ILP被視作電路的負載，每一時刻的負載電阻和電感參數可通過電磁場模型計算得到。

熱場以電磁場產生的焦耳熱和結構場的塑性功為熱源，在導體內產生溫升，并使熱量在發射器和ILP內傳導。熱場的計算結果可以傳遞到電磁場和結構場，對電導率、力學本構關系等材料屬性產生影響。

結構場以電磁場計算出的電磁力(洛倫茲力)為動力，在發射器和ILP內產生動力學響應，導致部件的變形、運動以及應力波在結構內的傳播。結構場的運動和變形被回饋給電磁場，用于更新導體構型、計算電磁場的演化。為了描述發射過程中材料的力學行為，結構場中定義了如下5種材料模型：1)彈性模型；2)彈塑性模型；3)Johnson-Cook模型；4)正交各向異性彈性模型；5)彈塑性流體動力學模型。這些模型基本能夠涵蓋電磁炮所涉及的金屬、非金屬、纖維纏繞結構、含能裝藥等材料。

圖1中的4個功能模塊既可以同時耦合運行，也可以單獨運行，或者選擇其中幾個組合運行，例如在剛體假設下可以不考慮結構場，只進行電磁、電路- 電磁耦合、電磁- 溫度耦合等計算。

求解器的各物理場共用一套網格。輸入量包括電路參數、有限元網格、接觸、材料屬性等，通過數據文件的形式讀入。由于在求解器中已經完成輸出量的計算和處理，輸出的數據文件可以直接導入作圖軟件來繪制曲線圖和云圖。

2 數學模型與計算方法

2.1 電磁場

用運動坐標下的磁擴散方程和電流連續性方程來描述導體內的電磁場，滿足

(1)

(2)

式中：A為矢量磁位；σ為電導率；μ為磁導率；φ為標量電位。當σ=0 S/m時，(1)式退化為Laplace方程，用于描述導體以外的絕緣體及自由空間。由(1)式和(2)式組成的電磁場方程組，可通過有限元- 邊界元耦合方法求解[15]。

在電磁炮導電路徑的某個截面上，施加具有一定波形的電流作為激勵條件，或者采用場路耦合的方式，獲得激勵電流。以脈沖電容器組驅動下的軌道炮發射系統為例，電路模型如圖2所示。圖2中，電源采用多模塊并聯方式，每個模塊k(k=1,2,…,n)中包含電容Ck、電感Lk、硅堆Dk、開關Kk，并計入雜散電阻Rk和RDk；軌道炮負載可視作可變電阻RL、可變電感LL以及電樞與軌道間接觸電壓Ua的串聯，負載電流i流過軌道和電樞，并在炮尾軌道兩端產生電壓Ub，ik為各模塊電流，R0為連接電源與軌道炮的輸電線電阻，L0為電感。

圖2 發射系統的電路模型Fig.2 Circuit model of launching system

當模塊的電容器電壓uCk>0 V時，電路的數學模型如下：

(3)

當uCk<0 V時，電路的數學模型如下：

(4)

式中：

(5)

2.1.1 負載電阻

負載電阻和電感通過電磁場模型計算。通過電磁場可以計算出電流密度j，將其在導體區域VC內積分后，可以得到負載內的焦耳熱功率為

(6)

式中：Ω為積分的體積微元。從電路的角度來看，焦耳熱功率為

We=i2RL.

(7)

由于電磁場是由電流驅動的，在負載電流i已知的情況下，可以將負載電阻表示為

(8)

2.1.2 負載電感

磁場能量密度瞬時值wm與磁通密度B和磁場強度H有關，定義為

(9)

對于電磁炮這種不含高磁導率材料的磁路，將磁場能量密度在整個求解域V內積分，得到磁場儲能為

(10)

由于磁場儲能與負載電感的關系為

(11)

在已知負載電流i的情況下，可以根據磁場儲能計算出負載電感為

(12)

在每個時間步，通過電磁場計算出負載電阻和負載電感，代入電路的常微分方程組，即(3)式～(5)式，然后用Ronge-Kutta方法求解出下一個時間步的負載電流。

2.2 熱場

由于電樞是運動的，熱場模型也被定義在運動坐標下，用傅里葉熱傳導方程描述為

(13)

式中：ρ為材料密度；c為材料比熱；T為溫度；kc為導熱系數；t為時刻；Q為熱源密度，來自于流過導體的電流所產生的焦耳熱以及塑性變形產生的熱量，

(14)

σs和ε分別為應力和應變，β為塑性功轉化為熱量的比例。熱場模型的求解采用有限元法[16]。

2.3 結構場

由于電磁炮發射時產生很高的應力，ILP與內膛之間存在接觸碰撞，而且發射周期只有短短的幾毫秒，因此采用顯式有限元方法[17]求解結構場。根據連續介質力學理論，結構場的控制方程組中包含質量守恒、動量守恒和能量守恒方程，滿足

(15)

采用罰函數法處理發射器和ILP的各部件之間存在著的接觸。檢查每個從節點是否穿透主面，如果存在穿透，則在從節點和接觸點之間施加界面力，其大小正比于穿透量，相當于在所有從節點和主表面之間布置一系列法向界面彈簧。另外，接觸面的摩擦也利用等效彈塑性彈簧來實現。

3 電磁軌道炮發射過程的數值仿真研究

3.1 計算模型及參數

圖3為計算模型的結構尺寸及有限元網格的示意圖，模型中考慮了發射器及ILP的基本結構特征。圖3(a)為發射器的截面，截面中心為ILP，D形軌道在其左右兩側相向布置，ILP上下兩側的矩形區域為絕緣體，發射器的外層為纖維纏繞層。圖3(b)為ILP，由電樞、卡瓣、彈丸(內部裝填炸藥)組成。

圖3 模型網格與尺寸Fig.3 Mesh and size of the model

發射器長6.0 m，軌道表面之間的最短距離為90 mm.ILP整體質量約5.1 kg，其中彈丸質量3.1 kg，電樞、卡瓣等寄生質量為2.0 kg.在該算例中，軌道、絕緣體、卡瓣為彈性材料，電樞為彈塑性材料，發射器的纏繞結構為正交各向異性彈性材料(材料參數見文獻[18])，裝藥為彈塑性流體動力材料(材料參數見文獻[19])。

脈沖電源的總儲能為60 MJ，由1 200個50 kJ儲能單元組成，單元參數如表1所示。

表1 儲能單元參數表Tab.1 Parameters of energy storage unit

儲能單元分為20組，通過設置不同的開關觸發延時進行時序放電，其中第1～4組在0 ms時刻觸發，其他各組依次延時300 μs觸發。

3.2 計算結果與討論

3.2.1 脈沖電源放電過程

圖4為脈沖電源的電流和電容器電壓曲線。圖4(b)為電容器兩端的電壓曲線，圖中標注的數字與延時觸發的第5～20組單元對應。

圖4 電流和電容器電壓曲線Fig.4 Current and capacitor voltage curves

在圖4(a)電流- 時間曲線的上升沿，即0.75 ms之前，只有第1～5組單元放電，它們的電容器電壓- 時間曲線顯示這部分單元的電容器上存在著剩余電壓，表明它們存儲的電能并未完全釋放。文獻[20]分析了這種現象的產生機理，認為它與晶閘管開關的關斷特性以及發射初期較高的負載電壓有關。圖5所示負載電壓(即炮尾電壓)Ub-時間曲線顯示了發射初期的高幅度電壓以及時序放電所造成的電壓波動。

圖5 負載電壓- 時間曲線Fig.5 Voltage-time curve of load

3.2.2 電磁場計算結果

圖6為軌道、電樞、彈體等導體上的電流密度J序列分布。高電流密度區主要集中在電樞及其后部的一段軌道上，這與趨膚效應及速度趨膚效應有關。另外，從圖6中可以看出，電樞和彈體沿發射方向運動，它們在不同時刻所到達的位置是通過結構場與電磁場的耦合計算獲得的。

圖6 電流密度的序列分布圖Fig.6 Sequential contours of current density

由于ILP會攜帶引信、裝藥等敏感部件，需要考慮電磁場對它們的影響。彈丸的殼體一般為導電材料，在脈沖磁場下會感應出渦流。將圖6的電流密度降低2個數量級，單獨顯示彈體上的電流密度分布，如圖7所示。在剛開始放電時(t=0.005 ms)，彈體上出現了明顯的感應電流，并且在靠近彈底的地方幅度較高。隨著時間推移，彈體上的感應電流逐漸變小。

圖7 彈體上的電流密度序列分布Fig.7 Sequential contours of current density on projectile

圖8所示彈體表面的磁通密度B分布，呈現出彈底高、彈尖低的分布特征。電樞和軌道組成的導電通路為非軸對稱結構，決定了彈體上的磁場也是非軸對稱的。在發射過程中，盡管磁通密度的幅度會有變化，但是分布模式沒有發生明顯的改變。通過對瞬態磁場的計算，可以為彈上敏感部件安裝位置及方向的設計提供參考。

圖8 彈體的磁通密度序列分布圖Fig.8 Sequential contours of flux density on projectile

電磁場計算的核心功能之一是獲得作用于電樞的電磁力。電磁力的三維分布特征在文獻[21]中曾作過討論，不再贅述，這里主要討論電磁力的合力。在電樞上，通過對電磁力密度的積分，可以獲得沿發射方向的合力F.根據

(16)

可以計算出軌道的等效電感梯度L′.另外，根據(9)式～(12)式，從磁場儲能的角度也可以計算出電感梯度L′1.將兩種方法得到的結果列于圖9.由圖9可見，在5.0 ms之前，兩條曲線是基本吻合的，但5.0 ms后，虛線表示的由受力計算出的等效電感梯度L′2有所下降，這段時間恰恰處于電流的下降沿，表明電樞的電磁力分布可能受到電流降低的影響而發生了變化。

圖9 電感梯度隨時間的變化曲線Fig.9 Variation of inductance gradient with time

3.2.3 熱場計算結果

電流在導體內產生的焦耳熱不斷累積，使溫度升高。圖10為ILP到達炮口時的電樞和彈體的溫升分布對比。由于導致溫升的熱源主要來自于電流密度所產生的焦耳熱，而且在短短幾毫秒的發射周期內，熱量來不及充分擴散，因此溫升分布呈現出明顯的三維特征，在圖10(a)中表現得尤其明顯，電樞上的局部溫升超過了600 K，最大溫升區位于電樞的喉部和電樞臂的尾端，與圖6的電流密度分布情況相一致，體現了電磁場對溫度場的作用。圖10(a)的電樞溫升分布表明，電樞喉部和電樞臂是潛在的易損傷部位，在發射試驗中經常會出現電樞喉部的磁鋸損傷和電樞臂磨損現象，與電樞溫升的分布特性有一定關系。

圖10 電樞與彈體的溫度分布Fig.10 Temperature contours of armature and projectile

由于彈體上的電流密度比電樞和軌道上的電流密度低幾個數量級，相比之下，圖10(b)所示彈體上的溫升要小得多，最大溫升區域位于彈底，由感應電流所導致。

3.2.4 結構場計算結果

結構場通過對連續介質力學方程組的求解，一方面追蹤ILP的膛內運動，并將每一時間步的位移和速度及時回饋給電磁場和溫度場，用于更新在ILP運動情況下的有限元網格構型；另一方面計算結構部件的微小變形、運動以及部件之間的相互作用，以獲得對結構設計有重要參考意義的應力參數。圖11為不同時刻的應力分布，圖11中顯示的是發射器和ILP的1/2模型。發射器和ILP內的應力主要源自于電磁力，由于電磁力產生于載流導體，隨著電樞的運動，載流軌道的長度不斷增加，發射器的高應力區也在不斷擴大。另外，當t分別為3.447 ms和5.091 ms時，不僅ILP右側的絕緣體上呈現出不均勻的應力分布，而且ILP左側的一段發射器上也有應力區在發展，這與應力波在發射器內的傳播有關。

圖11 發射器和ILP的應力分布與演化Fig.11 Distribution and evolution of stress on launcher and ILP

圖12為ILP上的應力分布。由圖4的電流曲線可知0.978 ms時電流最大，意味著電樞受到的電磁力最大，ILP承受的過載最大，因此該時刻的應力幅度最高，如圖12(a)所示。在隨后的兩個時刻中，由于電樞的推進力主要作用于彈丸殼體，而且其承力截面比較小，因此最大應力區位于彈丸的殼體上。從圖12(b)的ILP剖面圖上可以看出彈丸內部裝藥上的應力分布，它底部的應力最大，能夠達到百兆帕量級。圖12的應力分布呈現出一定的對稱性，而且三維分布特征在各個時刻表現得都很明顯，其產生原因是多方面的，以電樞為例，既與電磁力的三維分布有關，也是整個系統在電磁力作用下的結構響應結果。

圖12 ILP的應力分布與演化Fig.12 Distribution and evolution of stress on ILP

ILP及發射器的各部件組成了一個多體動力學系統。盡管發射過程中彈丸會發生變形和姿態改變，通過對模型網格的積分，仍可以得到其質心的運動規律。圖13為發射過程中彈丸的位移x和速度v隨時間的變化曲線，圖14為彈丸的橫向位移d和繞質心的轉動角度ω的變化曲線。這些曲線表明彈丸在膛內存在多個自由度的運動，尤其是橫向運動和擺動，盡管幅度很小，但也有可能在發射階段產生橫向過載，影響發射穩定性和安全性，這些因素在ILP的設計階段應予以考慮。

圖13 彈丸的速度和位移曲線Fig.13 Velocity and displacement curves of projectile

圖14 彈丸的橫向運動和擺動曲線Fig.14 Lateral motion and balloting curves of projectile

ILP的橫向運動也反映在表2所示接觸壓力p分布上。表2中列出了4個時刻的電樞- 軌道接觸面壓力云圖，可以看出接觸壓力并不是嚴格對稱分布的，與發射過程中各部件的變形以及接觸碰撞有關。

表2 電樞- 軌道接觸面上的壓力分布Tab.2 Pressure distribution on armature-rail contact surface

上述電磁軌道炮算例表明，求解器通過焦耳熱和電磁力載荷在物理場間的實時傳遞以及網格構型的實時更新，實現了電磁場、溫度場和結構場的耦合計算，不但能夠用電流、電壓、電阻、電感、速度等集總參數從系統的角度描述發射過程，而且能夠用電流密度、磁通密度、溫度、應力等場量從細節處描述各部件的物理場分布及演化，甚至可以捕捉到擺動、接觸碰撞、應力波等因素帶來的發射過程非理想特性，計算所獲得的信息可以為電磁軌道炮的設計提供參考。

4 同步感應線圈炮發射過程數值仿真研究

4.1 單級線圈炮

為驗證求解器在線圈炮數值仿真中的可用性，首先針對單級線圈炮，與文獻[22]給出的理論和實驗結果進行對比。算例的結構尺寸如圖15所示，電樞為鋁筒，其電導率為3.0×107S/m，線圈由60匝導線繞成，驅動電流如圖16所示。

圖15 單級線圈炮結構示意圖Fig.15 Schematic diagram of single-stage coilgun

圖16 單級線圈炮的驅動電流波形Fig.16 Current waveform for driving the single-stage coilgun

圖17列出了本文求解器的計算結果與文獻[22]結果的對比。由圖17可見：3條曲線基本吻合，發射初期，電樞受到前向推動力，速度迅速上升；在電樞將要離開線圈時則會受到拖拽力作用，使速度有所下降，之后作勻速運動。

圖17 計算和實驗結果的對比圖Fig.17 Comparison of calculated and experimental results

圖18從左至右分別給出了發射初期的電流密度、磁通密度，以及電磁力密度的矢量圖。從圖18(a)的電流密度矢量可以看出，電樞局部感應出與線圈電流密度方向相反的電流；從圖18(b)的磁通密度矢量可以看出，磁通密度較大的區域位于電樞和線圈之間；從圖18(c)的電磁力密度矢量圖可以看出，盡管電樞和線圈都受到了電磁力的作用，但電樞上的電磁力主要位于與線圈交疊的地方。

圖18 電流- 磁場- 電磁力的矢量圖Fig.18 Vectors of current density,flux density and EM force density

該算例結果表明，本文的求解器可以用于處理同步感應線圈炮這種含有運動導體的渦流場問題。

4.2 多級線圈炮

仿照文獻[23]的線圈炮尺寸，建立10級線圈炮的計算模型如圖19所示。線圈由20匝導線繞制，電樞為鋁制圓筒，設置為彈塑性材料，電樞與線圈之間為絕緣的導向管，設置為正交各向異性彈性材料，它與線圈緊密接觸，而與電樞之間保留了1 mm的間隙。各級線圈的驅動電流如圖20所示。

圖19 多級線圈炮模型Fig.19 Model of multi-stage coilgun

圖20 各級線圈的激勵電流曲線Fig.20 Driving current curves of coils

圖21為第1～4級線圈逐次放電時的電流密度分布圖，圖中顯示了線圈和電樞的3/4部分。由于線圈由多匝導線繞成，線圈上的電流密度均勻分布，而電樞上呈現出了密度不均的電流分布，高電流密度區出現在電樞的外側、靠近正在放電線圈的位置。

圖21 電流密度的序列分布Fig.21 Sequential contours of current density

圖22為線圈、電樞和導向管上的應力分布。由圖22可見，在初始的0.178 ms，電樞上的高應力區位于其尾部，隨著線圈的逐個放電，電樞上的應力分布也發生了改變。導向管上存在著應力分布，在2.270 ms時幅度比較大，這可能是電樞在較高速度下與導向管發生碰撞而造成的。

圖22 應力的序列分布Fig.22 Sequential contours of stress

圖23為電樞速度- 時間曲線。曲線上有多個鼓包，對應著電樞經過每一級線圈時所經歷的加速和減速過程。另外，該曲線是通過質心位移對時間的差分來計算的，由于電樞在發射過程中存在著變形和姿態的改變，計算出的速度- 時間曲線有些振蕩，特別是速度超過100 m/s后振蕩更明顯。

圖23 電樞的速度- 時間曲線Fig.23 Velocity-time curve of armature

由于電樞與導向管之間有1 mm的間隙，電樞會在管內擺動，并與管壁發生碰撞。圖24顯示了電樞的橫向位移dy、dz的變化過程。

圖24 橫向位移曲線Fig.24 Lateral displacement curves

增加電流是提高電樞速度的途徑之一，但帶來的負面影響是造成電樞的變形。圖25為將第1、2級線圈的放電電流峰值各自增加10 kA時出現的電樞變形。變形發生在電樞的尾部，變形模式與圖26所示文獻[24]的實驗現象非常相似。

圖25 大電流下的電樞塑性變形Fig.25 Plastic deformation of armature subjected to high current

圖26 變形后的電樞照片[24]Fig.26 Photo of deformed armature[24]

上述算例結果表明，本文的求解器有能力處理同步感應電磁線圈炮問題，它對電樞在脈沖磁場作用下的運動、變形過程的刻畫有助于為電磁線圈炮的設計提供參考依據。

5 結論

本文瞬態多物理場求解器基于電磁場、熱場、結構場的動力學模型，通過多場耦合、場路耦合，實現了對電流擴散過程和趨膚效應、熱傳導過程和溫升效應、結構動力學過程和多部件的力學接觸與沖擊效應的模擬，初步具備了針對電磁炮全系統、全尺寸、全發射周期的數值仿真能力。通過數值模擬，得出如下主要結論：

1) 電磁炮的發射過程是多物理場耦合、多部件相互作用下一個復雜的動力學過程。

2) 發射工況和發射性能與結構、材料以及激勵條件之間具有密切的相關性。

3) 求解器可用于模擬動態發射工況，確定發射器、一體化彈的危險點，對設計方案進行強度校核，評價設計方案的可行性。求解器的計算功能全部采用Fortran源代碼實現，因此具有良好的使用靈活性和功能可擴展性。

下一步，將結合電磁發射理論與技術的進展，繼續提高多物理場仿真計算的精細化程度，另外還將拓展該求解器的應用方向，使其能在更多的工程技術領域發揮作用。