林慶華，栗保明

(南京理工大學(xué) 瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210094)

0 引言

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

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

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

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

南京理工大學(xué)瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室開展電磁發(fā)射理論與技術(shù)研究已有十余年，逐步建立起涵蓋電磁場(chǎng)、熱場(chǎng)和結(jié)構(gòu)場(chǎng)的瞬態(tài)多物理場(chǎng)模型，編寫了電磁炮的求解器程序。本文簡要介紹了該求解器的組成框架、基本功能、數(shù)學(xué)模型和計(jì)算方法，通過電磁軌道炮和電磁線圈炮的具體算例，展示了該求解器對(duì)電磁發(fā)射過程的仿真計(jì)算能力。

1 求解器框架與功能介紹

求解器的組成框架如圖1所示。它針對(duì)的是包含脈沖電源、發(fā)射器和ILP在內(nèi)的整個(gè)發(fā)射系統(tǒng)，具有電路、電磁場(chǎng)、熱場(chǎng)、結(jié)構(gòu)場(chǎng)4個(gè)功能模塊，其中脈沖電源的放電過程用電路模型描述，發(fā)射器和ILP的工作過程用場(chǎng)模型描述。

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

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

熱場(chǎng)以電磁場(chǎng)產(chǎn)生的焦耳熱和結(jié)構(gòu)場(chǎng)的塑性功為熱源，在導(dǎo)體內(nèi)產(chǎn)生溫升，并使熱量在發(fā)射器和ILP內(nèi)傳導(dǎo)。熱場(chǎng)的計(jì)算結(jié)果可以傳遞到電磁場(chǎng)和結(jié)構(gòu)場(chǎng)，對(duì)電導(dǎo)率、力學(xué)本構(gòu)關(guān)系等材料屬性產(chǎn)生影響。

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

圖1中的4個(gè)功能模塊既可以同時(shí)耦合運(yùn)行，也可以單獨(dú)運(yùn)行，或者選擇其中幾個(gè)組合運(yùn)行，例如在剛體假設(shè)下可以不考慮結(jié)構(gòu)場(chǎng)，只進(jìn)行電磁、電路- 電磁耦合、電磁- 溫度耦合等計(jì)算。

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

2 數(shù)學(xué)模型與計(jì)算方法

2.1 電磁場(chǎng)

用運(yùn)動(dòng)坐標(biāo)下的磁擴(kuò)散方程和電流連續(xù)性方程來描述導(dǎo)體內(nèi)的電磁場(chǎng)，滿足

(1)

(2)

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

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

圖2 發(fā)射系統(tǒng)的電路模型Fig.2 Circuit model of launching system

當(dāng)模塊的電容器電壓uCk>0 V時(shí)，電路的數(shù)學(xué)模型如下：

(3)

當(dāng)uCk<0 V時(shí)，電路的數(shù)學(xué)模型如下：

(4)

式中：

(5)

2.1.1 負(fù)載電阻

負(fù)載電阻和電感通過電磁場(chǎng)模型計(jì)算。通過電磁場(chǎng)可以計(jì)算出電流密度j，將其在導(dǎo)體區(qū)域VC內(nèi)積分后，可以得到負(fù)載內(nèi)的焦耳熱功率為

(6)

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

We=i2RL.

(7)

由于電磁場(chǎng)是由電流驅(qū)動(dòng)的，在負(fù)載電流i已知的情況下，可以將負(fù)載電阻表示為

(8)

2.1.2 負(fù)載電感

磁場(chǎng)能量密度瞬時(shí)值wm與磁通密度B和磁場(chǎng)強(qiáng)度H有關(guān)，定義為

(9)

對(duì)于電磁炮這種不含高磁導(dǎo)率材料的磁路，將磁場(chǎng)能量密度在整個(gè)求解域V內(nèi)積分，得到磁場(chǎng)儲(chǔ)能為

(10)

由于磁場(chǎng)儲(chǔ)能與負(fù)載電感的關(guān)系為

(11)

在已知負(fù)載電流i的情況下，可以根據(jù)磁場(chǎng)儲(chǔ)能計(jì)算出負(fù)載電感為

(12)

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

2.2 熱場(chǎng)

由于電樞是運(yùn)動(dòng)的，熱場(chǎng)模型也被定義在運(yùn)動(dòng)坐標(biāo)下，用傅里葉熱傳導(dǎo)方程描述為

(13)

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

(14)

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

2.3 結(jié)構(gòu)場(chǎng)

由于電磁炮發(fā)射時(shí)產(chǎn)生很高的應(yīng)力，ILP與內(nèi)膛之間存在接觸碰撞，而且發(fā)射周期只有短短的幾毫秒，因此采用顯式有限元方法[17]求解結(jié)構(gòu)場(chǎng)。根據(jù)連續(xù)介質(zhì)力學(xué)理論，結(jié)構(gòu)場(chǎng)的控制方程組中包含質(zhì)量守恒、動(dòng)量守恒和能量守恒方程，滿足

(15)

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

3 電磁軌道炮發(fā)射過程的數(shù)值仿真研究

3.1 計(jì)算模型及參數(shù)

圖3為計(jì)算模型的結(jié)構(gòu)尺寸及有限元網(wǎng)格的示意圖，模型中考慮了發(fā)射器及ILP的基本結(jié)構(gòu)特征。圖3(a)為發(fā)射器的截面，截面中心為ILP，D形軌道在其左右兩側(cè)相向布置，ILP上下兩側(cè)的矩形區(qū)域?yàn)榻^緣體，發(fā)射器的外層為纖維纏繞層。圖3(b)為ILP，由電樞、卡瓣、彈丸(內(nèi)部裝填炸藥)組成。

圖3 模型網(wǎng)格與尺寸Fig.3 Mesh and size of the model

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

脈沖電源的總儲(chǔ)能為60 MJ，由1 200個(gè)50 kJ儲(chǔ)能單元組成，單元參數(shù)如表1所示。

表1 儲(chǔ)能單元參數(shù)表Tab.1 Parameters of energy storage unit

儲(chǔ)能單元分為20組，通過設(shè)置不同的開關(guān)觸發(fā)延時(shí)進(jìn)行時(shí)序放電，其中第1～4組在0 ms時(shí)刻觸發(fā)，其他各組依次延時(shí)300 μs觸發(fā)。

3.2 計(jì)算結(jié)果與討論

3.2.1 脈沖電源放電過程

圖4為脈沖電源的電流和電容器電壓曲線。圖4(b)為電容器兩端的電壓曲線，圖中標(biāo)注的數(shù)字與延時(shí)觸發(fā)的第5～20組單元對(duì)應(yīng)。

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

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

圖5 負(fù)載電壓- 時(shí)間曲線Fig.5 Voltage-time curve of load

3.2.2 電磁場(chǎng)計(jì)算結(jié)果

圖6為軌道、電樞、彈體等導(dǎo)體上的電流密度J序列分布。高電流密度區(qū)主要集中在電樞及其后部的一段軌道上，這與趨膚效應(yīng)及速度趨膚效應(yīng)有關(guān)。另外，從圖6中可以看出，電樞和彈體沿發(fā)射方向運(yùn)動(dòng)，它們?cè)诓煌瑫r(shí)刻所到達(dá)的位置是通過結(jié)構(gòu)場(chǎng)與電磁場(chǎng)的耦合計(jì)算獲得的。

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

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

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

圖8所示彈體表面的磁通密度B分布，呈現(xiàn)出彈底高、彈尖低的分布特征。電樞和軌道組成的導(dǎo)電通路為非軸對(duì)稱結(jié)構(gòu)，決定了彈體上的磁場(chǎng)也是非軸對(duì)稱的。在發(fā)射過程中，盡管磁通密度的幅度會(huì)有變化，但是分布模式?jīng)]有發(fā)生明顯的改變。通過對(duì)瞬態(tài)磁場(chǎng)的計(jì)算，可以為彈上敏感部件安裝位置及方向的設(shè)計(jì)提供參考。

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

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

(16)

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

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

3.2.3 熱場(chǎng)計(jì)算結(jié)果

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

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

由于彈體上的電流密度比電樞和軌道上的電流密度低幾個(gè)數(shù)量級(jí)，相比之下，圖10(b)所示彈體上的溫升要小得多，最大溫升區(qū)域位于彈底，由感應(yīng)電流所導(dǎo)致。

3.2.4 結(jié)構(gòu)場(chǎng)計(jì)算結(jié)果

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

圖11 發(fā)射器和ILP的應(yīng)力分布與演化Fig.11 Distribution and evolution of stress on launcher and ILP

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

圖12 ILP的應(yīng)力分布與演化Fig.12 Distribution and evolution of stress on ILP

ILP及發(fā)射器的各部件組成了一個(gè)多體動(dòng)力學(xué)系統(tǒng)。盡管發(fā)射過程中彈丸會(huì)發(fā)生變形和姿態(tài)改變，通過對(duì)模型網(wǎng)格的積分，仍可以得到其質(zhì)心的運(yùn)動(dòng)規(guī)律。圖13為發(fā)射過程中彈丸的位移x和速度v隨時(shí)間的變化曲線，圖14為彈丸的橫向位移d和繞質(zhì)心的轉(zhuǎn)動(dòng)角度ω的變化曲線。這些曲線表明彈丸在膛內(nèi)存在多個(gè)自由度的運(yùn)動(dòng)，尤其是橫向運(yùn)動(dòng)和擺動(dòng)，盡管幅度很小，但也有可能在發(fā)射階段產(chǎn)生橫向過載，影響發(fā)射穩(wěn)定性和安全性，這些因素在ILP的設(shè)計(jì)階段應(yīng)予以考慮。

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

圖14 彈丸的橫向運(yùn)動(dòng)和擺動(dòng)曲線Fig.14 Lateral motion and balloting curves of projectile

ILP的橫向運(yùn)動(dòng)也反映在表2所示接觸壓力p分布上。表2中列出了4個(gè)時(shí)刻的電樞- 軌道接觸面壓力云圖，可以看出接觸壓力并不是嚴(yán)格對(duì)稱分布的，與發(fā)射過程中各部件的變形以及接觸碰撞有關(guān)。

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

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

4 同步感應(yīng)線圈炮發(fā)射過程數(shù)值仿真研究

4.1 單級(jí)線圈炮

為驗(yàn)證求解器在線圈炮數(shù)值仿真中的可用性，首先針對(duì)單級(jí)線圈炮，與文獻(xiàn)[22]給出的理論和實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。算例的結(jié)構(gòu)尺寸如圖15所示，電樞為鋁筒，其電導(dǎo)率為3.0×107S/m，線圈由60匝導(dǎo)線繞成，驅(qū)動(dòng)電流如圖16所示。

圖15 單級(jí)線圈炮結(jié)構(gòu)示意圖Fig.15 Schematic diagram of single-stage coilgun

圖16 單級(jí)線圈炮的驅(qū)動(dòng)電流波形Fig.16 Current waveform for driving the single-stage coilgun

圖17列出了本文求解器的計(jì)算結(jié)果與文獻(xiàn)[22]結(jié)果的對(duì)比。由圖17可見：3條曲線基本吻合，發(fā)射初期，電樞受到前向推動(dòng)力，速度迅速上升；在電樞將要離開線圈時(shí)則會(huì)受到拖拽力作用，使速度有所下降，之后作勻速運(yùn)動(dòng)。

圖17 計(jì)算和實(shí)驗(yàn)結(jié)果的對(duì)比圖Fig.17 Comparison of calculated and experimental results

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

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

該算例結(jié)果表明，本文的求解器可以用于處理同步感應(yīng)線圈炮這種含有運(yùn)動(dòng)導(dǎo)體的渦流場(chǎng)問題。

4.2 多級(jí)線圈炮

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

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

圖20 各級(jí)線圈的激勵(lì)電流曲線Fig.20 Driving current curves of coils

圖21為第1～4級(jí)線圈逐次放電時(shí)的電流密度分布圖，圖中顯示了線圈和電樞的3/4部分。由于線圈由多匝導(dǎo)線繞成，線圈上的電流密度均勻分布，而電樞上呈現(xiàn)出了密度不均的電流分布，高電流密度區(qū)出現(xiàn)在電樞的外側(cè)、靠近正在放電線圈的位置。

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

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

圖22 應(yīng)力的序列分布Fig.22 Sequential contours of stress

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

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

由于電樞與導(dǎo)向管之間有1 mm的間隙，電樞會(huì)在管內(nèi)擺動(dòng)，并與管壁發(fā)生碰撞。圖24顯示了電樞的橫向位移dy、dz的變化過程。

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

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

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

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

上述算例結(jié)果表明，本文的求解器有能力處理同步感應(yīng)電磁線圈炮問題，它對(duì)電樞在脈沖磁場(chǎng)作用下的運(yùn)動(dòng)、變形過程的刻畫有助于為電磁線圈炮的設(shè)計(jì)提供參考依據(jù)。

5 結(jié)論

本文瞬態(tài)多物理場(chǎng)求解器基于電磁場(chǎng)、熱場(chǎng)、結(jié)構(gòu)場(chǎng)的動(dòng)力學(xué)模型，通過多場(chǎng)耦合、場(chǎng)路耦合，實(shí)現(xiàn)了對(duì)電流擴(kuò)散過程和趨膚效應(yīng)、熱傳導(dǎo)過程和溫升效應(yīng)、結(jié)構(gòu)動(dòng)力學(xué)過程和多部件的力學(xué)接觸與沖擊效應(yīng)的模擬，初步具備了針對(duì)電磁炮全系統(tǒng)、全尺寸、全發(fā)射周期的數(shù)值仿真能力。通過數(shù)值模擬，得出如下主要結(jié)論：

1) 電磁炮的發(fā)射過程是多物理場(chǎng)耦合、多部件相互作用下一個(gè)復(fù)雜的動(dòng)力學(xué)過程。

2) 發(fā)射工況和發(fā)射性能與結(jié)構(gòu)、材料以及激勵(lì)條件之間具有密切的相關(guān)性。

3) 求解器可用于模擬動(dòng)態(tài)發(fā)射工況，確定發(fā)射器、一體化彈的危險(xiǎn)點(diǎn)，對(duì)設(shè)計(jì)方案進(jìn)行強(qiáng)度校核，評(píng)價(jià)設(shè)計(jì)方案的可行性。求解器的計(jì)算功能全部采用Fortran源代碼實(shí)現(xiàn)，因此具有良好的使用靈活性和功能可擴(kuò)展性。

下一步，將結(jié)合電磁發(fā)射理論與技術(shù)的進(jìn)展，繼續(xù)提高多物理場(chǎng)仿真計(jì)算的精細(xì)化程度，另外還將拓展該求解器的應(yīng)用方向，使其能在更多的工程技術(shù)領(lǐng)域發(fā)揮作用。