張旭平，趙劍衡，譚福利，王桂吉，羅斌強(qiáng)，莫建軍，種 濤，孫承緯，劉倉理

（1.中國工程物理研究院流體物理研究所，四川 綿陽621999；2.中國工程物理研究院，四川 綿陽621999）

磁驅(qū)動加載技術(shù)是近十年發(fā)展起來的一種新的脈沖載荷加載技術(shù)，在高能量密度物理、沖擊動力學(xué)和航空航天等領(lǐng)域具有重要應(yīng)用［1－3］。對磁驅(qū)動飛片的數(shù)值模擬研究主要是在對飛片加載歷史計(jì)算［1，4］、飛片擊靶前狀態(tài)確定［5］、負(fù)載優(yōu)化設(shè)計(jì)［6］等方面。磁驅(qū)動過程中極板受洛倫茲力作用會產(chǎn)生嚴(yán)重的變形，并且極板加載面受焦耳熱而有很大的溫升，導(dǎo)致極板加載面有燒蝕、相變［7－8］。這些負(fù)載變化都會反饋到回路中影響磁驅(qū)動裝置對負(fù)載的放電電流。如果磁流體計(jì)算軟件不能耦合電路計(jì)算，程序的計(jì)算范圍只能是利用實(shí)驗(yàn)數(shù)據(jù)進(jìn)行計(jì)算和物理分析，則一般的磁流體動力學(xué)模擬計(jì)算只有在給定電流是實(shí)驗(yàn)電流的情況下，得到的結(jié)果才接近于真實(shí)情況，這樣對一些預(yù)測性的計(jì)算分析將非常困難。并且研究磁驅(qū)動裝置負(fù)載參數(shù)對電路的影響和對實(shí)驗(yàn)放電物理過程的認(rèn)識都需要采用耦合電路的磁流體動力學(xué)計(jì)算［9－10］。相對于國外磁流體動力學(xué)軟件的快速研發(fā)，國內(nèi)數(shù)值模擬方面工作還處在起步階段。國內(nèi)的模擬研究都是在磁流體動力學(xué)軟件中直接導(dǎo)入電流曲線進(jìn)行計(jì)算的［11－12］。

鑒于動態(tài)電路分析的重要性，本文中利用LS－DYNA980 MHD計(jì)算軟件和自編的電路分析程序，通過2個(gè)計(jì)算程序計(jì)算結(jié)果的迭代，建立一種等效耦合電路的磁流體動力學(xué)計(jì)算方法，并且利用建立好的計(jì)算方法對實(shí)驗(yàn)放電電流曲線進(jìn)行模擬，并分析磁驅(qū)動相同尺寸飛片時(shí)極板長度的選擇。

1 耦合電路計(jì)算方法和流程

1.1 電路計(jì)算模塊建立和參數(shù)確定

對磁驅(qū)動實(shí)驗(yàn)平臺電路分析時(shí)可以簡化為一個(gè)R－L－C電路［13］。脈沖功率加載裝置CQ－4的等效電路如圖1所示，第1部分為主電容、開關(guān)及第1段傳輸線的總等效電容C1、電感L1和電阻R1，第2部分為峰值電容部分的等效電容C2、電感L2和電阻R2，第3部分為負(fù)載區(qū)的等效電阻R3和電感L3。

根據(jù)電容、開關(guān)等的額定參數(shù)和磁驅(qū)動實(shí)驗(yàn)，可以確定R1、L1、C1、R2、L2、C2的值。R1中包含氣體開關(guān)的火花電阻r1，r1在放電過程中不斷變化，在放電初期r1較高，在放電回路電流很大時(shí)則很低。

根據(jù)等效的電路我們編了電路求解程序，程序主要采用四階龍格－庫塔方法求解微分方程組：

圖1 脈沖加載裝置電路圖Fig.1 Circuit of pulsed power generator

這個(gè)程序的優(yōu)點(diǎn)是電路中的各參數(shù)可以賦固定的值，也可以賦值為時(shí)變的表達(dá)式。

利用程序和擬定的R、L、C參數(shù)對CQ－4短路放電實(shí)驗(yàn) CQ－4 shot－1和 CQ－4 shot－2做了模擬計(jì)算，結(jié)果如圖2所示，分別是在CQ－4充電電壓73、80 k V下的短路電流曲線。程序計(jì)算的電流曲線和實(shí)驗(yàn)值符合較好。

圖2 模擬和實(shí)驗(yàn)電流曲線對比Fig.2 Comparison of simulated and measured load currents

1.2 計(jì)算流程和方法

磁驅(qū)動飛片的計(jì)算是一個(gè)典型的耦合電路的磁流體動力學(xué)計(jì)算問題。計(jì)算中以磁驅(qū)動裝置的充電電壓和負(fù)載區(qū)參數(shù)（極板構(gòu)型、飛片大?。┳鳛槌跏紬l件，得到裝置的放電電流和飛片速度等結(jié)果。一般耦合電路的磁流體動力學(xué)計(jì)算都是按時(shí)間步計(jì)算的（見圖3（a））。以裝置的充電電壓和靜態(tài)電路參數(shù)為初始條件，通過程序的電路模塊計(jì)算出加載到負(fù)載的電流值，利用程序的磁流體動力學(xué)模塊計(jì)算得到該時(shí)間步負(fù)載的電磁－力－熱響應(yīng)結(jié)果和新的負(fù)載電感、電阻值。接著將得到的新的負(fù)載參數(shù)反饋到下一個(gè)時(shí)間步的計(jì)算，作為下一個(gè)時(shí)間步的輸入條件，依次循環(huán)直到計(jì)算結(jié)束。圖中V是充電電壓，T為計(jì)算時(shí)間，R、L是負(fù)載電阻和電感，Lt、Rt、it分別表示t時(shí)刻的電感、電阻和電流值；MHD即電磁－力－熱耦合的磁流體動力學(xué)計(jì)算，加入材料的本構(gòu)關(guān)系、狀態(tài)方程、熱狀態(tài)方程、電導(dǎo)率模型等，通過求解流體動力學(xué)方程組、麥克斯韋方程組和熱傳導(dǎo)方程，得到負(fù)載的結(jié)構(gòu)變形、受力、溫度、電磁場分布等結(jié)果；CIRCUIT指電路分析計(jì)算，已知t時(shí)刻的參數(shù)Lt、Rt等，求解電路方程組（式（1）～（4）），得到電流值it。

由于建立以上耦合計(jì)算方法的磁流體動力學(xué)程序很困難，本文中在磁流體動力學(xué)程序LS－DYNA980 MHD和自建的電路計(jì)算程序的基礎(chǔ)上，利用實(shí)驗(yàn)時(shí)各參數(shù)（電流、電感L（t）、電阻R（t）等）對時(shí)間t的單值性，通過2個(gè)程序計(jì)算結(jié)果的迭代實(shí)現(xiàn)了耦合電路的磁流體動力學(xué)計(jì)算。單值指實(shí)驗(yàn)中任一時(shí)刻的所有參數(shù)（電流、電感、電阻等參數(shù)）都只有唯一值。用于數(shù)值計(jì)算即用t＋Δt時(shí)刻計(jì)算輸入的電流值I（t＋Δt）和用t時(shí)刻的整個(gè)電路的參數(shù)（電感L（t）、電阻R（t）等）求解的電流值是相等的，具體計(jì)算流程如圖3（b）所示。圖中t為計(jì)算時(shí)間，T為總計(jì)算時(shí)間，Δt為時(shí)間步長，n為循環(huán)次數(shù)，N 為最大循環(huán)次數(shù)，MHD是磁流體動力學(xué)程序LS－DYNA980 MHD，CIRCUIT為自編的外部電路計(jì)算程序，Iin是輸入磁流體計(jì)算軟件的電流曲線，Iout是CIRCUIT程序計(jì)算出的電流，In（t＋Δt）是第n次迭代電流曲線的（t＋Δt）時(shí)刻的電流值。

圖3 計(jì)算流程示意圖Fig.3 Schematic diagram of calculation flow

本文的計(jì)算流程（見圖3（b））與以上計(jì)算（見圖3（a））的主要區(qū)別是：第1，初始輸入條件是一條電流曲線；第2，不能在計(jì)算中將磁流體動力學(xué)計(jì)算的負(fù)載電路參數(shù)變化及時(shí)反饋到下一時(shí)間步的計(jì)算，而是磁流體和電路模塊分別計(jì)算后迭代。計(jì)算中以｜In（t＋Δt）out－In（t）in｜≤η作為單值原理的判據(jù)，表示兩曲線上對應(yīng)時(shí)刻的電流值相差小于η為止。否則下一循環(huán)輸入的電流曲線為

這種計(jì)算的收斂速度很快，一般循環(huán)次數(shù)N≤5時(shí)，η就已經(jīng)很小，能得到很好的結(jié)果。

2 負(fù)載動態(tài)參數(shù)計(jì)算

2.1 負(fù)載動態(tài)電感計(jì)算

耦合電路的磁流體動力學(xué)計(jì)算的難點(diǎn)是對負(fù)載時(shí)變電感和電阻的計(jì)算，磁驅(qū)動是一個(gè)電磁－力－熱耦合問題，計(jì)算中涉及到極板變形、磁擴(kuò)散和材料相變等問題。電流、電磁場的分布很難準(zhǔn)確計(jì)算，對電感和電阻計(jì)算也有一定的影響。

磁驅(qū)動飛片實(shí)驗(yàn)負(fù)載的原理示意圖如圖4所示，電流沿極板的內(nèi)表面?zhèn)鞑?，在極板之間形成磁場，

流過極板的電流與極板加載 面磁場相互作用后在兩極 板上產(chǎn)生相互排斥的磁壓 力。圖4中h為極板到短路端的長度，d為極板 之間的距離，a為從極板 電流入口處到飛片處的長 度，b為從飛片到極板短路處距離。在飛片發(fā)射過程 中兩極板之間間隙增大和 極板變形都會影響磁場的 分布和大小，導(dǎo)致負(fù)載電感變化。

圖4 磁驅(qū)動飛片負(fù)載物理模型Fig.4 Physics model of magnetically driven flyer

F.G.Steven等［9］在計(jì)算電感時(shí)采用了2種方法。第1種是加載的時(shí)變電感與加載時(shí)變電流有如下的關(guān)系：

式中：μ為磁導(dǎo)率，I（t）為電流，B（t）為磁場。利用2D－MHD對計(jì)算區(qū)域內(nèi)磁場能量積分的方法計(jì)算時(shí)變電感L（t）。第2種是等效計(jì)算一個(gè)完全電導(dǎo)體（PEC）條片模型的電場，PEC條片模型中條片之間中點(diǎn)的磁場強(qiáng)度是MHD計(jì)算的不同時(shí)間點(diǎn)的磁場強(qiáng)度，通過計(jì)算模型的電場能量計(jì)算電感：

式中：E/V是歸一化電場。以上2種方法電流為零時(shí)計(jì)算都存在奇點(diǎn)。

采用假定電流在極板表面均勻分布，對一端短路的平行條片的電感計(jì)算公式［14］：

式中：W 為極板寬度，h為極板長度，d為極板之間的距離。而其中公式極板之間的時(shí)變距離則通過三維的LS－DYNA980 MHD計(jì)算軟件計(jì)算。電流均勻分布于極板表面的假定，忽略了磁擴(kuò)散的影響和極板寬度方向電流的不均勻性，計(jì)算中需要對極板間隙d（t）做修正。

2.2 負(fù)載動態(tài)電阻計(jì)算

電阻計(jì)算采用下式：

式中：電流流過的橫截面積高度δ的計(jì)算采用磁擴(kuò)散深度代替，假設(shè)在這一磁擴(kuò)散深度內(nèi)阻值均勻。

從歐姆定律和麥克斯韋方程導(dǎo)出導(dǎo)體中的磁通密度的演化方程［14］：

假設(shè)只有寬度方向的磁場且電導(dǎo)率是固定值，寬度方向的磁場表示為

向電極板厚度方向磁擴(kuò)散深度和擴(kuò)散時(shí)間的關(guān)系也可表示為：

式中：σ為電導(dǎo)率，μ為磁導(dǎo)率。

將式（12）代入電阻計(jì)算公式（9）推出

實(shí)驗(yàn)中極板加載面受焦耳熱影響溫度升高，甚至極板受焦耳熱燒蝕發(fā)生相變，導(dǎo)致材料電導(dǎo)率在加載中發(fā)生變化，推導(dǎo)前面的公式是假定電導(dǎo)率不變。模擬中首先取2個(gè)電導(dǎo)率的區(qū)間節(jié)點(diǎn)計(jì)算了時(shí)變電阻（見圖5），一個(gè)值是電導(dǎo)率初值（曲線1），另一個(gè)是MHD計(jì)算時(shí)間終點(diǎn)時(shí)溫度和密度變化后的加載面電導(dǎo)率值（曲線2）。這2個(gè)電導(dǎo)率的值相當(dāng)于實(shí)驗(yàn)中的最大值和最小值，實(shí)驗(yàn)中負(fù)載的電阻值應(yīng)該在圖中計(jì)算的2條時(shí)變電阻的區(qū)間內(nèi)。而對于加載初期電導(dǎo)率接近于初值，對于加載后期電導(dǎo)率接近于下降后的值。計(jì)算的極板中電流流過的等效厚度內(nèi)，能量與溫度Tm成正比，能量與時(shí)間也成正比，溫度與時(shí)間成正比，電導(dǎo)率隨溫度的變化曲線可以擬合為一個(gè)衰減函數(shù)式

即在前面假設(shè)的條件下可以推出電導(dǎo)率隨時(shí)間的變化曲線，將其代入電阻計(jì)算公式（13），在前面的假定下磁驅(qū)動飛片過程電阻按圖中曲線3變化。電阻計(jì)算假定寬度方向電流均勻和在計(jì)算的磁擴(kuò)散厚度內(nèi)各處電導(dǎo)率保持均勻一致。

在前面電感和電阻的計(jì)算假定下，在磁擴(kuò)散較小和極板電流分布均勻的情況下能給出較好的計(jì)算結(jié)果。我們采用建立的耦合電路的磁驅(qū)動飛片的計(jì)算方法，對CQ－4的驅(qū)動飛片實(shí)驗(yàn)shot－13和shot－32的電流結(jié)果進(jìn)行了數(shù)值計(jì)算（見圖6），實(shí)驗(yàn)CQ－4 shot－13中飛片尺寸為10 mm×5 mm×0.7 mm，實(shí)驗(yàn)CQ－4 shot－32中飛片尺寸為?5 mm×0.35 mm。從圖6可以看出，計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果在飛片能測到速度的時(shí)間范圍內(nèi)符合較好。后期飛片脫離極板在鏜孔中飛行和極板加載面受焦耳熱燒蝕后形成等離子體層等，這些原因均可導(dǎo)致實(shí)驗(yàn)中電流不再沿飛片加載面?zhèn)鞑?，而是直接擊穿在等離子體中傳播，這時(shí)計(jì)算的電感和電阻值都出現(xiàn)偏差。

圖5 負(fù)載電阻變化曲線Fig.5 Resistance varied with time

圖6 驅(qū)動飛片實(shí)驗(yàn)電流和計(jì)算電流比較Fig.6 Measured and calculated currents of flyer plate driven experiments

3 負(fù)載極板長度設(shè)計(jì)

磁驅(qū)動飛片極板長度的設(shè)計(jì)主要受2個(gè)因素影響，第1是飛片沿極板長度方向的平面性，第2是回路電感。實(shí)驗(yàn)中電流在極板流入端有一個(gè)匯流過程，電流密度由小變大；在負(fù)載的短路端，電流流過垂直極板的短路板，形成的磁場會疊加到極板之間的磁場；這2種原因都導(dǎo)致極板加載面從電流入口端到短路端磁場逐漸增強(qiáng)，在發(fā)射飛片時(shí)會造成飛片加載面壓力不均勻，從而使飛片平面性變差。通過增加極板長度可以減弱這種不均勻性，但是極板長度增加會造成回路電感增大，導(dǎo)致加載電流減小，影響飛片速度。對只給定加載電流曲線計(jì)算極板長度影響加載面磁壓力均勻性的模擬計(jì)算，可以給出均勻性結(jié)果，但對飛片速度或加載電流的大小不能比較，不能起到優(yōu)化計(jì)算的效果。耦合電路的磁流體動力學(xué)計(jì)算不僅可以給出飛片平面性信息，也可以給出飛片速度歷史。

下面分別對相同飛片長度（l＝10 mm）、極板寬度（w＝6 mm），不同極板長度情況下飛片的平面性和速度進(jìn)行數(shù)值計(jì)算。第1種情況是極板長度h＝20 mm、a＝5 mm、b＝5 mm；第2種情況h＝25 mm、a＝10 mm、b＝5 mm；第3種情況h＝30 mm、a＝15 mm、b＝5 mm。與第1種情況相比后2種情況只是增加了入口端到飛片處的極板長度a（見圖4）。

圖7是3種極板長度在飛片發(fā)射實(shí)驗(yàn)充電電壓73 k V時(shí)發(fā)射飛片的自由面速度比較，圖8是飛片長度方向的平面性比較。由圖7～8可知，在加載時(shí)間t＝0.72μs時(shí)，極板長度h＝25 mm的飛片比h＝20 mm的飛片速度低約1.2 km/s，但是飛片的平面性有明顯的提高。極板長度h增加到30 mm時(shí)，相比于h＝25 mm時(shí)對飛片的平面性提高較小。所以認(rèn)為實(shí)驗(yàn)中極板長度取h＝25 mm比較合適。

圖7 不同極板長度h時(shí)飛片自由面速度比較Fig.7 Comparison of free surface velocities under different strip line lengths

圖8 極板長度h不同時(shí)沿極板長度方向飛片的平面性比較Fig.8 Comparison of the flatness under different strip line lengths

4 結(jié) 論

分析了磁驅(qū)動飛片實(shí)驗(yàn)和數(shù)值計(jì)算的研究現(xiàn)狀。根據(jù)實(shí)驗(yàn)裝置等效電路建立了可帶入?yún)?shù)為表達(dá)式的電路計(jì)算程序。利用建立的電路計(jì)算程序結(jié)合LS－DYNA980 MHD計(jì)算軟件，建立了一種可以耦合電路分析的磁驅(qū)動飛片計(jì)算方法。采用這種方法模擬了CQ－4的發(fā)射飛片實(shí)驗(yàn)，計(jì)算結(jié)果與實(shí)驗(yàn)基本符合。最后計(jì)算了極板長度h分別為20、25、30 mm在充電73 k V時(shí)磁驅(qū)動飛片的速度和平面性，認(rèn)為實(shí)驗(yàn)中極板長度取25 mm比較合適。

磁驅(qū)動飛片涉及到高壓、高溫、等離子體燒蝕等問題，數(shù)值計(jì)算中普遍存在的問題是材料本構(gòu)關(guān)系、狀態(tài)方程、電導(dǎo)率模型的選擇對計(jì)算結(jié)果的影響較大。文中對電感和電阻的計(jì)算由于涉及到磁擴(kuò)散和焦耳熱燒蝕，耦合電路的磁驅(qū)動計(jì)算結(jié)果的誤差也主要來自于電阻和電感的計(jì)算的誤差。計(jì)算流程為耦合電路的磁流體動力學(xué)計(jì)算模擬提供了一種新的耦合計(jì)算思路。

本文工作得到陳學(xué)秒、吳剛、蔡進(jìn)濤、稅榮杰、胥超、馬驍、鄧順益的幫助，在此特別感謝！

［1］Lemke R E，Knudson M D，Davis J P.Magnetically driven hyper－velocity launch capability at the Sandia Z accelerator［J］.International Journal of Impact Engineering，2011，38（6）：480－485.

［2］孫承緯，趙劍衡，王桂吉，等.磁驅(qū)動準(zhǔn)等熵平面壓縮和超高速飛片發(fā)射實(shí)驗(yàn)技術(shù)原理、裝置及應(yīng)用［J］.力學(xué)進(jìn)展，2012，42（2）：206－218.Sun Cheng－wei，Zhao Jian－h(huán)eng，Wang Gui－ji，et al.Progress in magnetic loading techniques for isentropic compression experiments and ultra－h(huán)igh velocity flyer launching［J］.Advances in Mechanics，2012，42（2）：206－218.

［3］趙劍衡，孫承緯，譚福利，等.一維平面磁驅(qū)動等熵加載發(fā)射飛片技術(shù)［J］.爆炸與沖擊，2005，25（4）：303－308.Zhao Jian－h(huán)eng，Sun Cheng－wei，Tan Fu－li，et al.Launch technique for isentropic compression flyer plates magnetically driven by using fast pulsed power［J］.Explosion and Shock Waves，2005，25（4）：303－308.

［4］Lemke R W，Knudson M D，Robinson A C，et al.Self－consistent，two－dimensional，magneto－h(huán)ydrodynamic simulations of magnetically driven flyer plates［J］.Physics of Plasmas，2003，10（5）：1867－1874.

［5］Lemke R W，Knudson M D，Bliss D E，et al.Magnetically accelerated，ultrahigh velocity flyer plates for shock wave experiments［J］.Journal of Applied Physics，2005，98（7）：073730.

［6］Ao T，Asay J R，Chantrenne S，et al.A compact strip－line pulsed power generator for isentropic compression experiments［J］.Review of Scientific Instruments，2008，79（1）：013903.

［7］Gael L B，Petit J，Chanal P Y，et al.Modelling the dynamic magneto－thermomechanical behaviour of materials using a multi－phase EOS［C］∥7th European LS－DYNA Conference.Salzburg，Austrian，2009.

［8］Gael L B，Chanal P Y，Herei P L.Ramp wave compresion in a copper strip line：Comprarison between MHD numerical simulations（LS－DYNA）and experimental results（GEPI device）［C］∥10th International LS－DYNA users conference.Dearborn，MI，USA，2008.

［9］Glover S F，Schneider L X，Reed K W，et al.Genesis：A 5 MA programmable pulsed power driver for isentropic compression experiments［C］∥IEEE Pulsed Power Conference（PPC）.Washington DC，2009：763－767.

［10］Glover S F，Davis J P，Schneider L X，et al.Impact of time－varying loads on the programmable pulsed power driver called genesis［J］.IEEE Transacts on Plasma Science，2012，40（10）：2588－2596.

［11］王剛?cè)A，孫承緯，趙劍衡，等.磁驅(qū)動平面飛片的一維磁流體力學(xué)計(jì)算［J］.爆炸與沖擊，2008，28（3）：261－264.Wang Gang－h(huán)ua，Sun Cheng－wei，Zhao Jian－h(huán)eng，et al.One－demensional，magnetohydrodynamic simulations of magnetically driven flyer plates［J］.Explosion and Shock Waves，2008，28（3）：261－264.

［12］王桂吉，蔣吉昊，孫承緯，等.磁驅(qū)動飛片的一維磁流體動力學(xué)數(shù)值研究［J］.計(jì)算力學(xué)學(xué)報(bào)，2008，25（6）：776－781.Wang Gui－ji，Jiang Ji－h(huán)ao，Sun Cheng－wei，et al.One－demensional，magneto－h(huán)ydrodynamic simulation of magnetically driven flyer plates［J］.Chinese Journal of Computational Mechanics，2008，25（6）：776－781.

［13］王桂吉，譚福利，王剛?cè)A，等.小型脈沖功率發(fā)生器放電電流波形的調(diào)節(jié)模式［J］.高電壓技術(shù)，2007，33（11）：223－226.Wang Gui－ji，Tan Fu－li，Wang Gang－h(huán)ua，et al.Shaping models of discharging current pulse for the compact pulsed power generator［J］.High Voltage Engineering，2007，33（11）：223－226.

［14］孫承緯.磁驅(qū)動等熵壓縮和高速飛片的實(shí)驗(yàn)技術(shù)［J］.爆轟波與沖擊波，2005（3）：125－138.Sun Cheng－wei.The technique of magnetically driven isentropic compression and high－velocity flyer plates［J］.Detonation and Shock Waves，2005（3）：125－138.