一種耦合電路分析的磁驅動飛片數值計算方法＊

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

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

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

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

1 耦合電路計算方法和流程

1.1 電路計算模塊建立和參數確定

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

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

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

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

這個程序的優點是電路中的各參數可以賦固定的值，也可以賦值為時變的表達式。

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

圖2 模擬和實驗電流曲線對比Fig.2 Comparison of simulated and measured load currents

1.2 計算流程和方法

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

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

圖3 計算流程示意圖Fig.3 Schematic diagram of calculation flow

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

這種計算的收斂速度很快，一般循環次數N≤5時，η就已經很小，能得到很好的結果。

2 負載動態參數計算

2.1 負載動態電感計算

耦合電路的磁流體動力學計算的難點是對負載時變電感和電阻的計算，磁驅動是一個電磁－力－熱耦合問題，計算中涉及到極板變形、磁擴散和材料相變等問題。電流、電磁場的分布很難準確計算，對電感和電阻計算也有一定的影響。

磁驅動飛片實驗負載的原理示意圖如圖4所示，電流沿極板的內表面傳播，在極板之間形成磁場，

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

圖4 磁驅動飛片負載物理模型Fig.4 Physics model of magnetically driven flyer

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

式中：μ為磁導率，I（t）為電流，B（t）為磁場。利用2D－MHD對計算區域內磁場能量積分的方法計算時變電感L（t）。第2種是等效計算一個完全電導體（PEC）條片模型的電場，PEC條片模型中條片之間中點的磁場強度是MHD計算的不同時間點的磁場強度，通過計算模型的電場能量計算電感：

式中：E/V是歸一化電場。以上2種方法電流為零時計算都存在奇點。

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

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

2.2 負載動態電阻計算

電阻計算采用下式：

式中：電流流過的橫截面積高度δ的計算采用磁擴散深度代替，假設在這一磁擴散深度內阻值均勻。

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

假設只有寬度方向的磁場且電導率是固定值，寬度方向的磁場表示為

向電極板厚度方向磁擴散深度和擴散時間的關系也可表示為：

式中：σ為電導率，μ為磁導率。

將式（12）代入電阻計算公式（9）推出

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

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

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

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

圖6 驅動飛片實驗電流和計算電流比較Fig.6 Measured and calculated currents of flyer plate driven experiments

3 負載極板長度設計

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

下面分別對相同飛片長度（l＝10 mm）、極板寬度（w＝6 mm），不同極板長度情況下飛片的平面性和速度進行數值計算。第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種極板長度在飛片發射實驗充電電壓73 k V時發射飛片的自由面速度比較，圖8是飛片長度方向的平面性比較。由圖7～8可知，在加載時間t＝0.72μs時，極板長度h＝25 mm的飛片比h＝20 mm的飛片速度低約1.2 km/s，但是飛片的平面性有明顯的提高。極板長度h增加到30 mm時，相比于h＝25 mm時對飛片的平面性提高較小。所以認為實驗中極板長度取h＝25 mm比較合適。

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

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

4 結 論

分析了磁驅動飛片實驗和數值計算的研究現狀。根據實驗裝置等效電路建立了可帶入參數為表達式的電路計算程序。利用建立的電路計算程序結合LS－DYNA980 MHD計算軟件，建立了一種可以耦合電路分析的磁驅動飛片計算方法。采用這種方法模擬了CQ－4的發射飛片實驗，計算結果與實驗基本符合。最后計算了極板長度h分別為20、25、30 mm在充電73 k V時磁驅動飛片的速度和平面性，認為實驗中極板長度取25 mm比較合適。

磁驅動飛片涉及到高壓、高溫、等離子體燒蝕等問題，數值計算中普遍存在的問題是材料本構關系、狀態方程、電導率模型的選擇對計算結果的影響較大。文中對電感和電阻的計算由于涉及到磁擴散和焦耳熱燒蝕，耦合電路的磁驅動計算結果的誤差也主要來自于電阻和電感的計算的誤差。計算流程為耦合電路的磁流體動力學計算模擬提供了一種新的耦合計算思路。

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

［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］孫承緯，趙劍衡，王桂吉，等.磁驅動準等熵平面壓縮和超高速飛片發射實驗技術原理、裝置及應用［J］.力學進展，2012，42（2）：206－218.Sun Cheng－wei，Zhao Jian－heng，Wang Gui－ji，et al.Progress in magnetic loading techniques for isentropic compression experiments and ultra－high velocity flyer launching［J］.Advances in Mechanics，2012，42（2）：206－218.

［3］趙劍衡，孫承緯，譚福利，等.一維平面磁驅動等熵加載發射飛片技術［J］.爆炸與沖擊，2005，25（4）：303－308.Zhao Jian－heng，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－hydrodynamic 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］王剛華，孫承緯，趙劍衡，等.磁驅動平面飛片的一維磁流體力學計算［J］.爆炸與沖擊，2008，28（3）：261－264.Wang Gang－hua，Sun Cheng－wei，Zhao Jian－heng，et al.One－demensional，magnetohydrodynamic simulations of magnetically driven flyer plates［J］.Explosion and Shock Waves，2008，28（3）：261－264.

［12］王桂吉，蔣吉昊，孫承緯，等.磁驅動飛片的一維磁流體動力學數值研究［J］.計算力學學報，2008，25（6）：776－781.Wang Gui－ji，Jiang Ji－hao，Sun Cheng－wei，et al.One－demensional，magneto－hydrodynamic simulation of magnetically driven flyer plates［J］.Chinese Journal of Computational Mechanics，2008，25（6）：776－781.

［13］王桂吉，譚福利，王剛華，等.小型脈沖功率發生器放電電流波形的調節模式［J］.高電壓技術，2007，33（11）：223－226.Wang Gui－ji，Tan Fu－li，Wang Gang－hua，et al.Shaping models of discharging current pulse for the compact pulsed power generator［J］.High Voltage Engineering，2007，33（11）：223－226.

［14］孫承緯.磁驅動等熵壓縮和高速飛片的實驗技術［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.