孫承緯，陸 禹，趙繼波，羅斌強(qiáng)，谷卓偉，王桂吉，張旭平，陳學(xué)秒，周中玉，李 牧，袁 紅，張紅平，王剛?cè)A，孫奇志，文尚剛，譚福利，趙劍衡，莫建軍，蔡進(jìn)濤，金云聲，賀 佳，種 濤，趙小明，劉倉理

（1. 中國(guó)工程物理研究院流體物理研究所，四川 綿陽 621999；2. 上海激光等離子體研究所，上海 201800；3. 深圳技術(shù)大學(xué)先進(jìn)材料測(cè)試技術(shù)研究中心，廣東 深圳 518118；4. 深圳技術(shù)大學(xué)大數(shù)據(jù)與互聯(lián)網(wǎng)學(xué)院，廣東 深圳 518118；5. 中國(guó)工程物理研究院化工材料研究所，四川 綿陽 621999；6. 中國(guó)工程物理研究院應(yīng)用電子學(xué)研究所，四川 綿陽 621999；7. 中國(guó)工程物理研究院，四川 綿陽 621999）

極端物理學(xué)顧名思義是極端條件下的物理學(xué)，但按當(dāng)今科學(xué)前沿的實(shí)際理解是高能量密度狀態(tài)下的物理學(xué)，即高能量密度物理[1]。高能量密度狀態(tài)的公認(rèn)條件是受載過程中物質(zhì)內(nèi)能的增量大于0.1 MJ/cm3，即其熱力學(xué)（靜水）壓力超過100 GPa[2]，或者波長(zhǎng)約1 μm 的輻射能束對(duì)物質(zhì)的輻照度大于3 PW/cm2。數(shù)十年前高能量密度狀態(tài)顯然只能在天體、行星、地球數(shù)千千米深部等環(huán)境以及核武器物理過程中存在，但產(chǎn)生高功率短脈沖激光的啁啾技術(shù)和美國(guó)Sandia 實(shí)驗(yàn)室100 TW 級(jí)功率低阻抗強(qiáng)電流驅(qū)動(dòng)器Z 機(jī)器的問世，使這種狀態(tài)已可在不少實(shí)驗(yàn)室內(nèi)實(shí)現(xiàn)。研究高能量密度狀態(tài)下物質(zhì)、輻射及其相互作用的多學(xué)科交叉的極端物理學(xué)前沿領(lǐng)域，對(duì)慣性聚變、天體物理、行星物理和核武器物理研究具有重要現(xiàn)實(shí)意義。

電磁驅(qū)動(dòng)下宏觀物體的高能量密度動(dòng)力學(xué)實(shí)驗(yàn)，包括極端材料動(dòng)力學(xué)、流體動(dòng)力學(xué)行為和高溫高密度等離子體產(chǎn)生等，均具有重要的技術(shù)背景和應(yīng)用價(jià)值[3-4]。從驅(qū)動(dòng)原理上看，這種實(shí)驗(yàn)可分為2 類：磁驅(qū)動(dòng)和磁壓縮。磁驅(qū)動(dòng)系利用導(dǎo)體構(gòu)形表面上強(qiáng)驅(qū)動(dòng)電流脈沖與自生磁場(chǎng)作用的電磁力加載，實(shí)現(xiàn)材料準(zhǔn)等熵（斜波）壓縮或發(fā)射超高速固體飛片（統(tǒng)稱isentropic compression experiment, ICE），可用于沖擊壓縮、驅(qū)動(dòng)金屬套筒內(nèi)爆、形成高溫高壓等離子體等。磁壓縮系利用炸藥爆轟或電磁驅(qū)動(dòng)金屬套筒的低速內(nèi)爆運(yùn)動(dòng)，把較大體積中的磁通量壓縮到小區(qū)域中，形成很高磁場(chǎng)和磁壓，把該小區(qū)域內(nèi)的材料（尤其是體積較大的低密度材料）準(zhǔn)等熵壓縮到高密度。利用炸藥爆轟為動(dòng)力的磁通量壓縮裝置，又稱為磁通量壓縮（MC-1）發(fā)生器。這2 類實(shí)驗(yàn)均涉及極端條件下實(shí)驗(yàn)構(gòu)形的力學(xué)與物理學(xué)耦合問題，難以從理論方程組上完善考慮，但能通過數(shù)值模擬途徑用計(jì)算解決。

1980 年代，美國(guó)Los Alamos 國(guó)家實(shí)驗(yàn)室研制了一維輻射磁流體力學(xué)編碼RAVEN[5]，具有一維任意拉格朗日-歐拉（arbitrary Lagrangian-Eulerian, ALE）功能，可以處理帶真空磁腔的負(fù)載構(gòu)形。RAVEN 編碼利用算子分裂方式，解決了流體彈塑性材料、雙流雙溫等離子體和高溫輻射物質(zhì)等強(qiáng)電磁加載下多物理場(chǎng)耦合模擬問題。目前，公認(rèn)的最先進(jìn)磁流體力學(xué)（magnetohydrodynamics, MHD）構(gòu)形計(jì)算編碼是美國(guó)Sandia 國(guó)家實(shí)驗(yàn)室的多介質(zhì)多物理場(chǎng)ALE 編碼ALEGRA 系列。1990 年代初，Sandia 實(shí)驗(yàn)室將高精度沖擊動(dòng)力學(xué)有限差分（歐拉）編碼CTH 與改進(jìn)的拉格朗日有限元編碼結(jié)合，得到在無結(jié)構(gòu)單元有限元網(wǎng)格上表述的ALE 編碼，即ALEGRA 最初的版本。經(jīng)過30 多年的開發(fā)，ALEGRA 已成為功能強(qiáng)大的計(jì)算輻射磁流體-沖擊動(dòng)力學(xué)-多物理場(chǎng)耦合的ALE 有限元編碼，得到了廣泛應(yīng)用[6]。

本文中敘述的一維沖擊爆轟和磁流體力學(xué)多物理場(chǎng)計(jì)算編碼SSS-MHD，可對(duì)電磁驅(qū)動(dòng)高能量密度動(dòng)力學(xué)實(shí)驗(yàn)開展磁流體力學(xué)多物理場(chǎng)數(shù)值模擬，分析強(qiáng)電流或者炸藥驅(qū)動(dòng)的瞬態(tài)電磁場(chǎng)加載、受載實(shí)驗(yàn)構(gòu)形和材料樣品中發(fā)生的彈塑性磁流體力學(xué)運(yùn)動(dòng)和熱力學(xué)變化等過程。SSS-MHD 編碼的基礎(chǔ)是1980 年代流體物理研究所研發(fā)的沖擊爆轟動(dòng)力學(xué)編碼SSS[7]，后又?jǐn)U展了激光輻照效應(yīng)[8-9]、等離子體電磁驅(qū)動(dòng)的計(jì)算功能[10-11]。隨著磁驅(qū)動(dòng)和磁壓縮實(shí)驗(yàn)的開展[12-13]，SSS 編碼也向MHD 擴(kuò)展，并于2012 年提出SSS-MHD 編碼的初稿[14-16]，10 年來經(jīng)過許多實(shí)例計(jì)算核對(duì)和修改完善后基本定稿[17-22]。圖1 所示的各類實(shí)驗(yàn)裝置基本概括了具有重要技術(shù)應(yīng)用背景的電磁驅(qū)動(dòng)實(shí)驗(yàn)類型，均可用SSS-MHD 編碼進(jìn)行不同程度的模擬。

圖1 適合SSS-MHD 編碼模擬的各種類型電磁驅(qū)動(dòng)高能量密度物理實(shí)驗(yàn)Fig. 1 Kinds of the magnetically-driven high-energy-density physics experiments suitably simulated with the SSS-MHD code

多物理場(chǎng)耦合的數(shù)值研究，主要體現(xiàn)在多個(gè)物理系統(tǒng)方程組在離散化計(jì)算中進(jìn)行恰當(dāng)?shù)模〝?shù)學(xué)上完備的）相互耦合計(jì)算。SSS-MHD 編碼主要理論框架是彈塑性磁流體力學(xué)方程組，自身就是力學(xué)與電磁學(xué)兩大系統(tǒng)的耦合，進(jìn)一步還能擴(kuò)展到與等離子體物理和輻射流體力學(xué)等學(xué)科的耦合。由于使用場(chǎng)景和目的不同，磁流體力學(xué)方程組形式多樣。圖1 的各類實(shí)驗(yàn)涉及低頻脈沖強(qiáng)電流和強(qiáng)磁場(chǎng)與導(dǎo)電體的短時(shí)間非相對(duì)論性相互作用，在電磁力和能量表達(dá)式中磁場(chǎng)的份額遠(yuǎn)高于電場(chǎng)。這種意義下的磁流體力學(xué)近似，只需要2 個(gè)電磁學(xué)變量，如電流密度和磁感應(yīng)強(qiáng)度，相關(guān)的方程式為磁擴(kuò)散方程和歐姆定律。上述耦合體現(xiàn)于運(yùn)動(dòng)方程中的電磁力、能量方程中的焦耳熱，同時(shí)磁擴(kuò)散方程與介質(zhì)運(yùn)動(dòng)速度有關(guān)，這種架構(gòu)為SSS-MHD 編碼結(jié)構(gòu)的簡(jiǎn)化創(chuàng)造了有利條件。

本文組成如下。第1 節(jié)，敘述彈塑性磁流體力學(xué)多物理場(chǎng)耦合基本方程組；第2 節(jié)，建立作為SSSMHD 編碼基礎(chǔ)的拉格朗日一維方程組的具體形式；第3 節(jié)，敘述SSS-MHD 編碼的創(chuàng)新結(jié)構(gòu)和功能；第4 節(jié)，展示平面構(gòu)形ICE 實(shí)驗(yàn)的SSS-MHD 模擬，包括鋁、鉭、PBX 炸藥的準(zhǔn)等熵斜波壓縮和鋁飛片高速發(fā)射實(shí)驗(yàn)；第5 節(jié)，展示圓柱面構(gòu)形的固體套筒電磁內(nèi)爆實(shí)驗(yàn)及炸藥內(nèi)爆磁通量壓縮（MC-1）發(fā)生器實(shí)驗(yàn)的模擬計(jì)算；第4～5 節(jié)將SSS-MHD 編碼計(jì)算與Sandia 實(shí)驗(yàn)室Z 裝置的3 個(gè)高水平實(shí)驗(yàn)和計(jì)算實(shí)例進(jìn)行比較，而且模擬流體物理研究所CQ、CJ 系列裝置的4 個(gè)典型實(shí)驗(yàn)；第6 節(jié)的結(jié)束語中指出，上述算例表明SSS-MHD 模擬與實(shí)驗(yàn)數(shù)據(jù)以及ALEGRA 計(jì)算數(shù)據(jù)的相對(duì)偏差基本不超過5%。本文中的數(shù)值模擬還以豐富的物理量剖面圖顯示實(shí)驗(yàn)過程全面的物理圖像。編碼SSS-MHD 的創(chuàng)新特色以及原有的沖擊動(dòng)力學(xué)編碼的牢固基礎(chǔ)，使其模擬通用性強(qiáng)、功能靈活全面、計(jì)算速度快、精度高。編碼SSS-MHD 的開發(fā)不但能為極端物理學(xué)、極端材料動(dòng)力學(xué)實(shí)驗(yàn)提供有力的數(shù)值模擬平臺(tái)，而且有助于多維MHD 多物理場(chǎng)編碼的研發(fā)。

1 彈塑性磁流體力學(xué)基本方程組

連續(xù)介質(zhì)力學(xué)運(yùn)動(dòng)的基本規(guī)律為質(zhì)量、動(dòng)量和能量守恒定律，從一般的連續(xù)介質(zhì)動(dòng)力學(xué)方程組可以方便地導(dǎo)出磁流體力學(xué)的基本方程組[26]。首先，在動(dòng)量方程和能量方程中分別增加與電磁力和電阻產(chǎn)生熱量（焦耳熱）相關(guān)的項(xiàng)，把向量形式的連續(xù)介質(zhì)動(dòng)力學(xué)方程組改寫如下，依次為連續(xù)性方程（質(zhì)量守恒）、運(yùn)動(dòng)方程（動(dòng)量守恒）和能量方程（能量守恒）：

將式(2)全部點(diǎn)乘以u(píng)，再與式(3)相加，得到以總比能E=e+u·u/2 表達(dá)的能量方程：

式(4)是SSS-MHD 編碼實(shí)際使用的能量方程，其優(yōu)點(diǎn)是簡(jiǎn)化了張量運(yùn)算。后文中，能量方程中無外體力F項(xiàng)，外源項(xiàng)Q寫成WL即單位質(zhì)量介質(zhì)吸收的激光輻射功率。

在沖擊動(dòng)力學(xué)或爆炸力學(xué)中，對(duì)連續(xù)介質(zhì)通常采用流體彈塑性的力學(xué)模型加以描述。將應(yīng)力張量Σ分解為 Σ=-pI+T，寫成分量式：

式中：比容v=1/ρ 。SSS-MHD 編碼計(jì)算中材料EOS 基本均采用完全物態(tài)方程數(shù)據(jù)庫SESAME[27]。SESAME 庫中包括熱力學(xué)和低溫等離子體性質(zhì)，為SSS-MHD 編碼多物理計(jì)算提供了重要基礎(chǔ)。由式(5)～(6)和未寫出的本構(gòu)關(guān)系，溫度以及應(yīng)力張量各個(gè)分量都可以通過微元的比內(nèi)能、比容、粒子速度分量及其導(dǎo)數(shù)或積分表達(dá)，此時(shí)未知的力學(xué)函數(shù)尚有e、 ρ 和u，式(1)～(3)恰好是2 個(gè)標(biāo)量和1 個(gè)矢量方程，只要其中的電磁量J和B能借助于電動(dòng)力學(xué)方程組耦合求解，即可實(shí)現(xiàn)方程組數(shù)學(xué)上的封閉。

考慮電動(dòng)力學(xué)之始，先假定介質(zhì)中不存在自由電荷，無需考慮位移電流和電介質(zhì)的極化或磁化，并得到電場(chǎng)強(qiáng)度E和磁感應(yīng)強(qiáng)度B都是無源的： ?·E=0 ， ?·B=0 。麥克斯韋方程組還有如下定律。

法拉第定律：

安培定律：

再增加廣義歐姆定律：

前面方程組已含有粒子速度u為未知函數(shù)，式(7)～(9)提供了求解E、J和B等3 個(gè)電磁量的方程，其中 μ 為介質(zhì)的磁導(dǎo)率。由于安培定律，式(1)～(3)中僅出現(xiàn)磁場(chǎng)B，在該意義下系統(tǒng)被稱為磁流體力學(xué)方程組。進(jìn)一步化簡(jiǎn)，可從式(7)～(9)中消去E和J，得到?jīng)Q定B的磁擴(kuò)散方程：

至此，式(1)～(6)、(8)、(10)以及未寫明的介質(zhì)材料本構(gòu)關(guān)系，組成了一般情況下多維彈塑性磁流體力學(xué)的封閉方程組，而且是用矢量、張量和微分算子表述的一般形式，提供了寫出任意正交曲線坐標(biāo)系中分量形式微分方程組的方便途徑。

2 一維拉格朗日形式彈塑性磁流體力學(xué)方程組

根據(jù)第1 節(jié)中的基本方程組，可寫出一維情形下平面、柱面、球面坐標(biāo)系中歐拉形式的彈塑性磁流體力學(xué)方程組（分量式），如下。連續(xù)性方程：

運(yùn)動(dòng)方程：

能量方程：

在一維磁流體力學(xué)情形中，幾何指數(shù)N只能取0 或1，SSS-MHD 編碼的電磁學(xué)方程組（式(9)～(10)）的分量形式如下。安培定律：

磁擴(kuò)散方程組：

式中：Cθ=rNBθ稱為磁感應(yīng)強(qiáng)度的廣義 θ 分量，在柱對(duì)稱情形中有重要意義。

在歐拉-拉格朗日形式的轉(zhuǎn)換中，時(shí)間t保持不變，拉格朗日空間坐標(biāo)采用質(zhì)量坐標(biāo)M，即計(jì)算構(gòu)形自最左邊界點(diǎn)到格點(diǎn)處的累計(jì)質(zhì)量。歐拉空間坐標(biāo)r成為格點(diǎn)拉格朗日坐標(biāo)M的位置函數(shù)，位于r處格點(diǎn)的質(zhì)量 ?M則是厚度為dr、底面積分別是單位面積（平面）、r× 單位寬度（柱面）和r2（球面）的（單位弧度和立體角）體積中包含的介質(zhì)質(zhì)量。因此，2 個(gè)坐標(biāo)系之間的轉(zhuǎn)換關(guān)系是：

為了加以區(qū)分，將拉格朗日形式方程組中的r改寫為R，粒子速度u改寫為U，徑（縱）向應(yīng)力 σrr簡(jiǎn)寫為 σ 。每個(gè)拉格朗日格點(diǎn)的質(zhì)量保持不變體現(xiàn)了質(zhì)量守恒，也省去了連續(xù)性方程。為此，需要添加一個(gè)速度方程，即U等于R的時(shí)間偏導(dǎo)數(shù)。

由于高速度、高壓力問題參數(shù)范圍的特點(diǎn)，沖擊爆轟編碼SSS 采用了“沖擊波單位制”，其基本單位是質(zhì)量g、長(zhǎng)度cm 和時(shí)間μs，因而速度單位是cm/μs 即10 km/s，壓力單位是g/(cm · μs2)=100 GPa。前面推導(dǎo)的電磁學(xué)方程組采用國(guó)際單位制（SI），因此在耦合方程組中電磁量相關(guān)項(xiàng)需添加修正系數(shù)，達(dá)到與力學(xué)量單位的自洽。磁導(dǎo)率 μ=μ′μ0，其中真空磁導(dǎo)率 μ0=4π×10-7H/m ， μ′為介質(zhì)的相對(duì)磁導(dǎo)率。

經(jīng)過上述推導(dǎo)和轉(zhuǎn)換，得到SSS-MHD 編碼使用的拉格朗日形式一維彈塑性磁流體力學(xué)方程組如下。

連續(xù)性方程和速度方程：

運(yùn)動(dòng)方程：

能量方程：

安培定律（Cθ≡rNBθ）：

磁擴(kuò)散方程組：

對(duì)式(18)～(23)做離散時(shí)，時(shí)間離散格式是簡(jiǎn)單的前向差分，空間離散中格點(diǎn)序號(hào)自左向右排列，除了每一格的位置、速度、電流密度定義在格子中點(diǎn)以外，其他力學(xué)量和電磁量均定義在格子右點(diǎn)。因此，速度方程和安培定律是中心差分格式，其余方程都是前向差分格式。整個(gè)差分方程組的具體形式可參看文獻(xiàn)[5]。

3 SSS-MHD 編碼的結(jié)構(gòu)和功能

3.1 計(jì)算編碼的結(jié)構(gòu)

SSS-MHD 編碼的結(jié)構(gòu)立足于SSS 程序原有的框架[7]，保持彈塑性流體動(dòng)力學(xué)計(jì)算的各個(gè)子程序模塊構(gòu)架，進(jìn)行適應(yīng)MHD 計(jì)算的改造。首先，在導(dǎo)程序中增加MHD 建模的內(nèi)容，在總卡（GENERL）和區(qū)卡（COMPNT）中引入與MHD 計(jì)算相關(guān)的數(shù)組、變量和指數(shù)等，然后，在專門的MHD 建模子程序（MHDZON）里再對(duì)有關(guān)區(qū)卡賦以MHD 負(fù)載構(gòu)形的類型指數(shù)、構(gòu)形參數(shù)以及初始條件等，完成MHD 計(jì)算建模。接著，在相關(guān)子程內(nèi)部和計(jì)算循環(huán)節(jié)點(diǎn)處加入控制MHD 計(jì)算的語句，如在計(jì)算循環(huán)開始處設(shè)定邊界條件狀態(tài)及參數(shù)的同時(shí)，檢查和調(diào)整空腔和磁腔的狀態(tài)及參數(shù)。

隨著力學(xué)計(jì)算流程的進(jìn)行，在運(yùn)動(dòng)方程子程序（VELOC）、能量方程子程序（ENERGY）和熱輸運(yùn)計(jì)算子程序（HETCON）中，使用上時(shí)刻電磁變量值分別做洛倫茲力和焦耳熱計(jì)算。力學(xué)計(jì)算完成之后，流程進(jìn)入專設(shè)的電磁計(jì)算程序包（MHDBLK，后詳）。最后是數(shù)據(jù)更新（UPDATE）、后處理等，進(jìn)入下一個(gè)時(shí)間循環(huán)，直至計(jì)算時(shí)間t達(dá)到設(shè)定時(shí)刻ts。編碼的計(jì)算流程見圖2。

圖2 SSS-MHD 編碼的計(jì)算流程Fig. 2 Flowchart of the SSS-MHD code

原編碼SSS 是一個(gè)框架式的通用程序，具有多種物理場(chǎng)耦合計(jì)算的功能，可以包容或擴(kuò)展各種物態(tài)方程和物性數(shù)據(jù)庫，加入有關(guān)局部性物理過程的計(jì)算（例如激光與物質(zhì)的熱和沖量耦合、平衡電離的Saha 方程、動(dòng)態(tài)斷裂的NAG 模型和含能材料的燃燒爆炸反應(yīng)等），適合平、柱、球?qū)ΨQ幾何的動(dòng)力學(xué)問題，設(shè)有多種邊界條件和初始條件可供選擇，能方便地實(shí)現(xiàn)沖擊或爆轟的加載模擬[7]，這些特點(diǎn)已全部為SSS-MHD 編碼沿用。下面各小節(jié)將具體闡述SSS-MHD 計(jì)算編碼的關(guān)鍵內(nèi)容和特色。

3.2 MHD 負(fù)載構(gòu)形的統(tǒng)一建模

在磁流體力學(xué)計(jì)算模型的建模中，由于力學(xué)與電磁量具有不同方向的坐標(biāo)分量，首先需要設(shè)定力學(xué)構(gòu)形和相應(yīng)的加載磁場(chǎng)類型，然后賦予各個(gè)MHD 計(jì)算區(qū)域及格點(diǎn)處的電磁學(xué)物性參數(shù)和相應(yīng)初始值，必要時(shí)在整體模型邊界處補(bǔ)充電磁學(xué)邊界條件。其中的實(shí)質(zhì)性工作是負(fù)載構(gòu)形及磁場(chǎng)的分類和表述。

圖3 是SSS-MHD 編碼規(guī)定的實(shí)驗(yàn)負(fù)載構(gòu)形種類，N是一維連續(xù)性方程式(11)中的幾何指數(shù)。這樣的一維幾何中介質(zhì)速度只有徑（縱）向分量，但電流和磁場(chǎng)可以有2 個(gè)方向，即軸向（z）和角向（θ）分量。磁擴(kuò)散相應(yīng)地有2 個(gè)方向的方程式，還可以用2 個(gè)驅(qū)動(dòng)電路分別提供z、θ 方向的加載電流。一維磁流體力學(xué)原則上可以計(jì)算z、θ 及其組合的螺旋磁力線相關(guān)問題。負(fù)載構(gòu)形是各種材料的多層平行平板或者同軸圓柱殼層組成的幾何結(jié)構(gòu)，其間可以包括空腔（有磁場(chǎng)的空腔稱為磁腔）。建模依據(jù)主要的電磁驅(qū)動(dòng)機(jī)制，對(duì)構(gòu)形加載或被其壓縮的磁場(chǎng)（磁力線）有軸向（z）或角向（θ）以及兩者組合（z、θ 螺旋）等3 種形式。螺旋磁場(chǎng)中軸向與角向分量各自成體系、并行不悖，雖然圖中未畫出這個(gè)類別。從計(jì)算體系來說，只有圖3 的4 種基本構(gòu)形，其中平面z和θ 等2 種構(gòu)形實(shí)際相同，構(gòu)形的名稱表示其加載電流方向。構(gòu)形的板、殼可以是多層組合，也沒有具體畫出。

圖3 SSS-MHD 編碼中實(shí)驗(yàn)負(fù)載構(gòu)形的種類（紅線表示回流導(dǎo)體）[29]Fig. 3 Types of experimental configurations in the SSS-MHD code(red curves standing for return conductors)[29]

圖1 中各類實(shí)驗(yàn)的構(gòu)形和磁場(chǎng)均可抽象為圖3 的類型。負(fù)載構(gòu)形中的磁場(chǎng)可以是加載電流感生的（持續(xù)勵(lì)磁），也可以是初始時(shí)由外界給定的（瞬態(tài)種子磁場(chǎng)或定常勵(lì)磁）；推動(dòng)構(gòu)形運(yùn)動(dòng)的可以是加載電流自身形成的洛倫茲力（磁驅(qū)動(dòng)），也可以利用炸藥爆轟或沖擊的力量（磁壓縮）。

3.3 空腔和磁腔

SSS-MHD 編碼沿用了SSS 編碼的空腔功能，并擴(kuò)展到帶磁場(chǎng)的磁腔情形。計(jì)算模型通常為順序排列的介質(zhì)分區(qū)所組成的單連通區(qū)，如果隊(duì)列中存在真空間隔并不妨礙歐拉形式計(jì)算，但對(duì)于拉格朗日計(jì)算則發(fā)生“零質(zhì)量格點(diǎn)”的困難，通常需采取分區(qū)分階段計(jì)算合成或?qū)嵭幸痪SALE 策略，這樣可能使得編碼復(fù)雜、機(jī)時(shí)增多，計(jì)算精度也受影響。

SSS 編碼采取物質(zhì)區(qū)與空腔區(qū)一體化拉格朗日計(jì)算的方法。圖4 表示把空腔（區(qū)I）處理為只有3 個(gè)零質(zhì)量格點(diǎn)（J至J+2）的特殊區(qū)，不參加物質(zhì)區(qū)運(yùn)算，但其左右兩岸點(diǎn)J和J+3 分別具有左右物質(zhì)區(qū)岸點(diǎn)位置坐標(biāo)R的數(shù)據(jù)，空腔間隙等于R(J+3)-R(J) 。若空腔間隙大于零，左右兩岸點(diǎn)J和J+3 處均取自由邊界條件，空腔中間兩點(diǎn)J+1 和J+2 分別是J和J+3 點(diǎn)的鬼點(diǎn)（影子點(diǎn)）。如果空腔間隙為零或負(fù)值，意味著空腔閉合甚至兩岸碰撞過沖，需要退到前一時(shí)刻，適當(dāng)減小時(shí)間步長(zhǎng)重算，使得重算的空腔間隙在規(guī)定誤差下為零，并使點(diǎn)J和點(diǎn)J+3 的R值保持相等，成為同一個(gè)右岸介質(zhì)內(nèi)點(diǎn)。如果空腔原來是閉合的，但表觀重合的左、右岸點(diǎn)相互作用的拉力已變?yōu)槌^一定閾值的斥力，則空腔將被自動(dòng)拉開。對(duì)于空腔開、閉變動(dòng)的暫態(tài)振蕩過程，程序給予致穩(wěn)處理。大量算例表明，SSS 編碼的空腔一體化計(jì)算功能比較成功，步驟簡(jiǎn)便，精度較高。

圖4 SSS 編碼中空腔區(qū)及鄰區(qū)中格點(diǎn)的設(shè)置Fig. 4 Meshes in the cavity and nearby in the code SSS

磁流體力學(xué)計(jì)算中磁空腔較常見，其力學(xué)運(yùn)動(dòng)可沿用SSS 編碼的處理方法，但需增加電磁學(xué)計(jì)算，并保持兩岸點(diǎn)磁場(chǎng)及磁力對(duì)稱相等的處理。經(jīng)過改進(jìn)，SSS-MHD 編碼實(shí)現(xiàn)了一維磁流體力學(xué)多連通區(qū)的拉格朗日一體化計(jì)算，具備了對(duì)于較復(fù)雜實(shí)驗(yàn)及裝置進(jìn)行數(shù)值模擬的能力。

磁腔處理方法基于電磁學(xué)一維問題的特點(diǎn)并利用了SSS 編碼的基礎(chǔ)。如圖3 所示，運(yùn)用積分形式的安培定律和法拉第定律可以證明，實(shí)驗(yàn)構(gòu)形中磁腔界面處的磁場(chǎng)總是與相應(yīng)電流分量正交的，磁腔左右岸導(dǎo)體表面的電流只能是面電流?？涨唬ɑ螂娊橘|(zhì)）中磁場(chǎng)的生成，無外乎饋入加載電流或初始種子磁場(chǎng)以及構(gòu)形導(dǎo)體中產(chǎn)生感應(yīng)電流等幾類原因，都可由磁擴(kuò)散方程組計(jì)算確定。由于不考慮靜電電荷和位移電流，電介質(zhì)的電磁學(xué)性質(zhì)等同于磁腔。如果最右端的回流導(dǎo)體質(zhì)量較大或可不加考慮，則可簡(jiǎn)化為最右端的單純磁腔（力學(xué)上是自由邊界，不是空腔）。此時(shí)的計(jì)算模型稱為單邊驅(qū)動(dòng)，加載電流線只要一條，還需設(shè)定右端磁腔條件來計(jì)算該腔的磁通量。圓柱面θ 構(gòu)形可能存在軸線處的中心磁腔，但只能是環(huán)形電流驅(qū)動(dòng)的均勻軸向磁場(chǎng)，如MC-1 裝置的情形。

3.4 磁腔磁場(chǎng)的計(jì)算

圖3 表明在平面一維情形（N= 0）、以一對(duì)無限大平行電極板為界的平磁腔（即陽陰極間隙AK-gap）中，加載電流可以有2 個(gè)方向，分別產(chǎn)生均勻磁場(chǎng)Bθ或Bz，磁力線均為直線；圓柱面一維（N= 1）幾何具有θ 或z等2 種構(gòu)形，同軸的無限長(zhǎng)雙圓柱筒的間隙或單柱筒內(nèi)部均為圓（環(huán)）磁腔，可近似推廣到雙多邊形柱筒之間的環(huán)形磁腔（環(huán)腔）情形。因加載電流有θ 或z方向，環(huán)腔中磁力線則為軸向直線或同軸環(huán)線。按照?qǐng)D3 畫出的磁力線圍線，由安培定律給出平直腔和圓（環(huán)）腔的磁場(chǎng)與構(gòu)形加載電流的關(guān)系，得到平直磁腔和圓（環(huán)）磁腔中的磁場(chǎng)分別為：

式中：I為構(gòu)形的總加載電流，i為該電流經(jīng)過單位寬度表面的平均電流（線密度），w為腔寬度或環(huán)腔的周線長(zhǎng)度，下標(biāo)表示電流或磁場(chǎng)的方向。按照一維幾何的原意，平直腔是無限寬的，總電流I無限大，式(24)是用有限寬度平均電流近似無限寬的理想狀況，帶來的誤差修正下面說明。式(25)適用于圓腔，由于非圓形環(huán)腔也是封閉的，可以近似采用圓環(huán)腔來處理，從而環(huán)腔的直線段部分相當(dāng)于理想化的平直腔，雖然其電流密度因需按周線長(zhǎng)度平均計(jì)算而大為下降。

柱筒構(gòu)形按式(25)的近似很接近實(shí)際的柱面二維情況，誤差不大。然而，理想平面一維的式(24)與實(shí)際構(gòu)形的三維電磁場(chǎng)以及實(shí)驗(yàn)結(jié)果的差異就較顯著。式(24)第一個(gè)等號(hào)右邊的因子 (I/w) 即以磁腔實(shí)際寬度w計(jì)算的平均加載電流密度，相當(dāng)于用電流密度 (I/w) 、寬度無限的理想磁腔做近似。而如圖3所示，實(shí)際的平面一維“單條片”構(gòu)形（磁驅(qū)動(dòng)實(shí)驗(yàn)的一對(duì)加載電極板）是“折疊”后供加載電流來回的狹長(zhǎng)金屬片，其間的AK-gap 即是磁腔。從能量角度看，式(24)的近似顯然是對(duì)磁腔中電流密度、磁場(chǎng)及磁壓的一種高估。陸禹[30]利用商業(yè)電磁學(xué)軟件[31]進(jìn)行了二維平面電磁場(chǎng)計(jì)算，以頻率0.35 MHz 的正弦電流加載于一對(duì)銅電極板（截面寬2 mm、厚1 mm），兩板間的磁腔截面寬度（w）為2 mm、間隙（g）為0.1 mm。此無限長(zhǎng)“單條片構(gòu)形”電磁響應(yīng)的二維平面高頻計(jì)算結(jié)果表示，加載電流（幅值）密度分布明顯不均勻，集中于有限尺度磁腔截面的4 個(gè)角點(diǎn)，證實(shí)磁腔中部電流密度低于邊部。這種計(jì)算提供了對(duì)有限尺度平直磁腔中電流密度分布不均勻性的定量估計(jì)。定義K為等效電流因數(shù)，即利用商業(yè)電磁學(xué)軟件模擬計(jì)算得到的平面磁腔中部幅值電流密度與全寬度平均幅值電流密度的比值。將式(24)（不計(jì)下標(biāo)）修正為：

式中：修正后的電流密度iK=KI/w。平面一維磁腔的等效電流因數(shù)K是一個(gè)小于1 的正數(shù)，與磁腔幾何參數(shù)w、g以及加載電流參數(shù)有關(guān)，如果w→∞,g→0 ，則K→1 。

嚴(yán)格說來，上述等效電流因數(shù)只反映實(shí)驗(yàn)開始時(shí)的情況，由于運(yùn)動(dòng)中構(gòu)形動(dòng)態(tài)電感的變化需要進(jìn)行難以做到的至少是二維的磁擴(kuò)散計(jì)算，才能對(duì)電流分布作更適當(dāng)?shù)目紤]。動(dòng)態(tài)電感的增大降低了加載電流的幅度，也可看作為運(yùn)動(dòng)過程中因數(shù)K的繼續(xù)下降，然而實(shí)際上MHD 計(jì)算包括了這部分的相互作用（詳見3.6 節(jié)），使得K動(dòng)態(tài)變化的影響大為減弱，絕大部分計(jì)算中只需沿用初始的K值，即可得到全程與實(shí)驗(yàn)數(shù)據(jù)相符的良好結(jié)果。這就是SSS-MHD 編碼處理平面磁腔的方法。

本節(jié)中論述的平面磁腔磁場(chǎng)近似計(jì)算方法對(duì)于SSS-MHD 編碼的磁擴(kuò)散計(jì)算具有重要意義，因?yàn)榇徘恢械拇艌?chǎng)Bz或Cθ都是常值，根據(jù)式(24)～(26)可由構(gòu)形加載總電流的即時(shí)值直接確定。該總電流（時(shí)間函數(shù)）可以作為輸入數(shù)據(jù)給定，也可從電路方程組計(jì)算獲得。因而，這些磁腔的磁場(chǎng)值可以用作為磁擴(kuò)散方程計(jì)算中可靠的定標(biāo)值，避免了隱式計(jì)算的麻煩和誤差。該問題在理論上有更嚴(yán)謹(jǐn)?shù)姆椒?，如Lemke 等[32]利用Poisson 方程計(jì)算了含有四邊形環(huán)形磁腔的多連通區(qū)平面二維磁流體力學(xué)問題，但并不能得到有普遍意義的結(jié)果。

3.5 顯式計(jì)算格式

SSS-MHD 編碼主要用來模擬MHD 實(shí)驗(yàn)構(gòu)形的力學(xué)和電磁運(yùn)動(dòng)，進(jìn)而研究樣品材料的物態(tài)方程和材料動(dòng)力學(xué)問題。為了滿足數(shù)值計(jì)算要求的穩(wěn)定性條件（如SSS 編碼中的Courant 條件），提高計(jì)算精度，并可防止鄰近格點(diǎn)尺度差距過大引起的剛性問題，希望采取一體化的拉格朗日顯式運(yùn)算。力學(xué)方程組可以從給定邊界條件的左端點(diǎn)開始做顯式計(jì)算。把磁擴(kuò)散方程同樣離散為顯式有限差分格式，但其邊界條件不在左、右端點(diǎn)給定，需要通過迭代使得本時(shí)刻的整體磁場(chǎng)分布、磁通量、反電動(dòng)勢(shì)（構(gòu)形的動(dòng)態(tài)電感）與加載或勵(lì)磁電流值達(dá)到高度自洽，才能獲得較高的計(jì)算精度，具體做法見3.7 節(jié)的說明。

避免計(jì)算不穩(wěn)定性和剛性問題的重要措施是編碼中介質(zhì)區(qū)可選擇多種方法進(jìn)行分格。計(jì)算值的振蕩往往發(fā)生于幾何或物理因素造成鄰近格子尺度或物理量急劇變化的區(qū)域，如內(nèi)爆聚心階段的中心范圍、點(diǎn)爆炸和散心爆轟開始階段的核心區(qū)域、高速碰撞的接觸面附近、脈沖電流的集膚層、高功率激光的吸收層和燒蝕層等。如果把這些區(qū)域做小尺度分格，并平滑過渡到周圍正常尺度分格的鄰區(qū)，則有可能避免振蕩，得到物理上合理的計(jì)算結(jié)果。SSS-MHD 編碼中除了通常的等步長(zhǎng)、等質(zhì)量分格方式外，還有在指定分區(qū)中使格子尺度（或質(zhì)量）按一定比例增大或縮小的分格方式。激光吸收層處理和從冷問題開始的等離子體演化等SSS 編碼計(jì)算，就是成功實(shí)例[8-9]。

3.6 SSS-MHD 計(jì)算與驅(qū)動(dòng)電路方程的耦合

負(fù)載構(gòu)形的動(dòng)態(tài)電感變化與驅(qū)動(dòng)器放電過程之間存在重要的相互作用，若要理解這些復(fù)雜過程，首先應(yīng)把MHD 計(jì)算與驅(qū)動(dòng)器電路計(jì)算進(jìn)行耦合，建立準(zhǔn)確預(yù)測(cè)實(shí)際加載電流的能力。與驅(qū)動(dòng)電路的耦合是MHD 計(jì)算能否解決實(shí)際問題的一個(gè)關(guān)鍵，也是SSS-MHD 編碼的重要特色。這個(gè)要求至今仍是對(duì)強(qiáng)電磁驅(qū)動(dòng)實(shí)驗(yàn)計(jì)算模擬能力的重大挑戰(zhàn)。

簡(jiǎn)單實(shí)驗(yàn)裝置中驅(qū)動(dòng)電流的傳輸及成形通常由RLC 電路描述，其計(jì)算是常微分方程組的初值問題。復(fù)雜的大型裝置使用脈沖形成介質(zhì)線技術(shù)，可由一階偏微分方程組描述，并可離散化為線性代數(shù)方程組求解。在與磁流體力學(xué)計(jì)算耦合方面，兩者沒有實(shí)質(zhì)性差別。下面以當(dāng)前使用的集中參數(shù)RLC 電路為例，說明與MHD 計(jì)算的耦合的方法。先把RLC 電路寫為：

式中：第1 個(gè)等號(hào)左邊為電容器組的剩余電壓，Q和C0分別為其剩余電荷量和初始電容值；U0和I分別為充電電壓和負(fù)載電流；L0、LL、R和VS依次為外電路（負(fù)載構(gòu)形）的固定電感、其他可變電感、固定電阻和開關(guān)電壓降（伏安特性），均在數(shù)據(jù)文件中輸入；ε 為負(fù)載構(gòu)形對(duì)于驅(qū)動(dòng)電路的反電動(dòng)勢(shì)，同一構(gòu)形可以有z、θ 等2 個(gè)驅(qū)動(dòng)電路和2 套電路參數(shù)。

驅(qū)動(dòng)電路與負(fù)載構(gòu)形耦合的關(guān)鍵是計(jì)算反電動(dòng)勢(shì)ε，即構(gòu)形路端電壓（電感性）VAK和沿其電流回路的電阻性電壓降VLP之和。根據(jù)積分形式的法拉第定律，VAK為構(gòu)形總磁通量 Φ 的時(shí)間導(dǎo)數(shù)，即VAK=dΦ/dt；根據(jù)歐姆定律，有VLP=ηlLPI（或相應(yīng)的積分式），lLP為構(gòu)形電流回路長(zhǎng)度。任何時(shí)刻，電容器組剩余電壓減去各項(xiàng)電壓降之后應(yīng)等于負(fù)載構(gòu)形的反電動(dòng)勢(shì)。從物理上說，VAK就是構(gòu)形磁通量變化在圍繞其磁場(chǎng)的導(dǎo)體圍線中感應(yīng)的路端電壓。因此，電路耦合不但要求解電路方程，更需要求解磁擴(kuò)散方程組得到構(gòu)形中的磁場(chǎng)分布，進(jìn)行分區(qū)和全構(gòu)形的磁通量計(jì)算。

3.7 SSS-MHD 計(jì)算程序包MHDBLK

力學(xué)計(jì)算之后就可計(jì)算在上一時(shí)刻應(yīng)力、電磁力和各種能量作用下、在已發(fā)生力學(xué)變形并且即時(shí)負(fù)載電流作用下的構(gòu)形中電磁場(chǎng)的變化，也就是求解本時(shí)刻的電路方程和磁擴(kuò)散方程組，這就是SSSMHD 編碼中體現(xiàn)MHD 多場(chǎng)耦合的程序包MHDBLK 的工作（圖5）。它的功能有3 項(xiàng)：(1)電路方程子程序（DLRK），依據(jù)上一時(shí)刻（j）已迭代n次的反電動(dòng)勢(shì)ε 和負(fù)載電流I值解算電路方程（調(diào)用子程CURRNT 可設(shè)定該方程形式），得到負(fù)載電流n+1 次迭代值In+1；(2)磁擴(kuò)散方程組子程序（MAGDFU），依據(jù)算出的In+1，求得磁場(chǎng)定標(biāo)點(diǎn)處值，計(jì)算全部構(gòu)形范圍的磁場(chǎng)分布；(3)依據(jù)得到的磁場(chǎng)分布，計(jì)算負(fù)載構(gòu)形n+1 次迭代的總磁通量及反電動(dòng)勢(shì)ε 的迭代值。若ε 迭代收斂不符合要求，則把n+1 次I和ε 的值輸入DLRK 子程序，進(jìn)行下一次迭代計(jì)算，直到第n次和第n+1 次的反電動(dòng)勢(shì)差值小于千分之一，才輸出本時(shí)刻j+1 的負(fù)載電流和磁場(chǎng)分布的確定值。

圖5 SSS-MHD 編碼中計(jì)算程序包MHDBLK 的組成Fig. 5 Flow chart for the routine package MHDBLK in the SSS-MHD code

若無法獲得驅(qū)動(dòng)器電路方程的資料，加載電流也可選擇指定電流饋入方式，如實(shí)驗(yàn)測(cè)量的實(shí)際負(fù)載電流或者另行估算的負(fù)載電流。此時(shí)，指定電流的數(shù)據(jù)以電流函數(shù)（CURFOM）形式饋入，直接把本時(shí)刻j+1 的指定電流值輸入MAGDFU，算得磁場(chǎng)后直接輸出，與電路無關(guān)，不進(jìn)行迭代。SSS-MHD 編碼整體為顯式運(yùn)算，但通過同一時(shí)刻的負(fù)載電流與磁腔磁場(chǎng)的自洽需通過迭代計(jì)算，可得到與隱式運(yùn)算相同的效果，而且精度更高。按照這種方法，把上時(shí)刻j到本時(shí)刻j+1 的迭代擴(kuò)大到整個(gè)力學(xué)計(jì)算范圍也是可行的，但增加如此多的計(jì)算量對(duì)提高精度不一定有意義。

4 SSS-MHD 編碼應(yīng)用于平面構(gòu)形磁驅(qū)動(dòng)實(shí)驗(yàn)的典型算例

SSS-MHD 編碼針對(duì)各類實(shí)驗(yàn)研究以及國(guó)外高水平實(shí)驗(yàn)結(jié)果，進(jìn)行了大量驗(yàn)算和改進(jìn)工作。本節(jié)列舉的平面構(gòu)形磁驅(qū)動(dòng)實(shí)驗(yàn)（ICE）算例分為2 類：材料樣品準(zhǔn)等熵壓縮和高速金屬飛片發(fā)射。

圖6 是上述2 類磁驅(qū)動(dòng)實(shí)驗(yàn)構(gòu)形的原理示意圖，圖6(a)～(b)都是單條片（strip-line）構(gòu)形；圖6(c)是方柱筒構(gòu)形，其芯柱（陰極）和側(cè)面板（陽電極）相當(dāng)于4 個(gè)單條片，它們并聯(lián)構(gòu)成了四邊筒形立體結(jié)構(gòu)。樣品自由面或窗口界面速度歷史的測(cè)量采用激光干涉技術(shù)，如VISAR（適用于任意反射面的速度干涉儀）、PDV（光子多普勒測(cè)速儀）等。當(dāng)材料樣品受到很強(qiáng)壓縮脈沖時(shí)，其自由表面會(huì)發(fā)生微噴射現(xiàn)象，干擾激光干涉測(cè)速過程，這種情況必須采用力學(xué)阻抗適當(dāng)?shù)墓鈱W(xué)單晶體片（LiF、藍(lán)寶石等）作為窗口，嚴(yán)密覆蓋樣品表面以抑制噴射，測(cè)量的數(shù)據(jù)則是樣品/窗口界面的速度歷史。

圖6 磁驅(qū)動(dòng)準(zhǔn)等熵壓縮實(shí)驗(yàn)的典型負(fù)載構(gòu)形Fig. 6 Typical loading configurations for magnetically-driven isentropic compression experiments (ICE)

4.1 磁驅(qū)動(dòng)平面材料樣品準(zhǔn)等熵壓縮實(shí)驗(yàn)（算例1～4）

算例1 是對(duì)美國(guó)Sandia 實(shí)驗(yàn)室Z 裝置實(shí)驗(yàn)Z-1220[33]的驗(yàn)算，該實(shí)驗(yàn)采用類似圖6(c)的6061-T6 鋁合金4 層方筒構(gòu)形，將一個(gè)側(cè)面（陽電極板）加工為4 段樣品，厚度h分別為0.616、0.868、1.118 和1.371 mm，相當(dāng)于一個(gè)4 階的臺(tái)階靶。每段樣品外表面粘有鍍?cè)鐾改さ腖iF 光學(xué)窗口。樣品與窗口的界面速度用VISAR 測(cè)量。實(shí)驗(yàn)加載電流峰值約18 MA，上升沿約300 ns。美國(guó)Z-1220 實(shí)驗(yàn)數(shù)據(jù)以及ALEGRA-1D 編碼與SSS-MHD 編碼的計(jì)算模擬都示于圖7。

圖7 對(duì)于美國(guó)Sandia 實(shí)驗(yàn)室磁驅(qū)動(dòng)等熵壓縮實(shí)驗(yàn)Z-1220[33]的SSS-MHD 模擬（算例1）Fig. 7 SSS-MHD simulations for theisentropic compression experiment Z-1220[33]at the Sandia National Laboratories, USA (example 1)

SSS-MHD 計(jì)算是依據(jù)圖7(a)的加載電流實(shí)驗(yàn)波形進(jìn)行的，此例的矩形環(huán)腔周長(zhǎng)5.8 cm，換算得知平均電流密度峰值約3.1 MA/cm，等效電流因數(shù)K為1。4 種厚度樣品的2 種數(shù)值模擬速度曲線與實(shí)驗(yàn)數(shù)據(jù)的符合程度都較好。SSS-MHD 計(jì)算給出的材料樣品最高壓力為51 GPa，文獻(xiàn)[33]給出的最高縱向應(yīng)力為54 GPa。2 種模擬的差別在于SSS-MHD 計(jì)算在彈塑性轉(zhuǎn)變階段速度起跳過早，其原因可能是該計(jì)算的樣品材料是純鋁而不是Z 實(shí)驗(yàn)的6061-T6 鋁合金，此外還需要考察計(jì)算中輸入電流波形初始50 ns 段落的細(xì)節(jié)以及鋁樣品中磁場(chǎng)擴(kuò)散的情況，進(jìn)而或可做出相應(yīng)改進(jìn)。對(duì)于模擬速度曲線末期與實(shí)驗(yàn)數(shù)據(jù)的偏差，主要是因?yàn)殇X的彈塑性行為與模擬選用的材料動(dòng)態(tài)本構(gòu)模型之間存在差異。

算例2～4 是對(duì)流體物理研究所CQ 系列脈沖功率裝置上相關(guān)實(shí)驗(yàn)[14,16,29]的驗(yàn)算，材料樣品包括重金屬鉭、輕金屬鋁和以HMX 為基的塑料黏結(jié)炸藥（PBX）JO-9 159，實(shí)驗(yàn)參數(shù)及結(jié)果列于表1。

表1 磁驅(qū)動(dòng)準(zhǔn)等熵壓縮實(shí)驗(yàn)算例2～4 的主要參數(shù)[14,16]Table 1 Parameters of examples 2－4 for magnetically-driven isentropic compression experiments[14,16]

算例2～4 加載電流和樣品表面速度的實(shí)驗(yàn)和計(jì)算結(jié)果示于圖8。由于CQ 系列裝置早期的實(shí)驗(yàn)樣品寬度較窄，算例2（Ta）和算例4（JO-9 159）的等效電流因數(shù)均取為0.6，算例3（Al）的等效電流因數(shù)取為0.4；后續(xù)電流較高的算例中，該因數(shù)大多為0.7 以上。SSS-MHD 編碼具有2 種加載電流輸入方式：電路（CKT）方程組耦合計(jì)算（見圖5）和實(shí)驗(yàn)負(fù)載電流（CUR）數(shù)據(jù)饋入，圖8 將這2 種計(jì)算方式做了比較。早期工作都是使用CKT 計(jì)算的[14-16]，隨著強(qiáng)電流測(cè)量技術(shù)的提高，本文中補(bǔ)充了CUR 計(jì)算，得到了更全面的認(rèn)識(shí)。

圖8 對(duì)于流體物理研究所CQ 裝置磁驅(qū)動(dòng)準(zhǔn)等熵壓縮實(shí)驗(yàn)的SSS-MHD 模擬（算例2～4）Fig. 8 SSS-MHD simulations for magnetically-driven quasi-isentropic compression experiments at IFP (examples 2－4)

將圖8 中3 種材料樣品臺(tái)階靶速度波形處理后，可獲得它們的準(zhǔn)等熵壓縮線，其高端壓力如表1 右端所示。算例2 中鉭金屬樣品的密度高，同時(shí)該例的加載電流更強(qiáng)，使得鉭的等熵線高端達(dá)80 GPa 以上。根據(jù)圖8(a)可知，如果加載電流密度提高3 倍多，等熵壓縮壓力可提高一個(gè)量級(jí)，能得到重金屬1 TPa以下范圍內(nèi)的等熵線，正如美國(guó)Z 裝置上開展的實(shí)驗(yàn)[34]。當(dāng)然，要把磁驅(qū)動(dòng)加載電流強(qiáng)度提高若干倍，對(duì)驅(qū)動(dòng)器硬件和實(shí)驗(yàn)技術(shù)都極為困難。

表1 中提及了電流計(jì)算（CUR）與電路計(jì)算（CKT）方式的比較。SSS-MHD 編碼CKT 計(jì)算采用的集中參數(shù)RLC 電路可以描述很多實(shí)際驅(qū)動(dòng)電路，但電路中若使用伏安特性復(fù)雜的氣體開關(guān)，則RLC 電路方程就不夠準(zhǔn)確。圖8(d)中CKT 計(jì)算優(yōu)于CUR 方式，原因可能是該算例驅(qū)動(dòng)器CQ-1.5 裝置采用了伏安性能簡(jiǎn)單的爆炸開關(guān)[35]，正適合于CKT 計(jì)算的簡(jiǎn)單RLC 電路描述（早期實(shí)驗(yàn)使用的電流測(cè)量技術(shù)也存在較大誤差）。其余2 例的驅(qū)動(dòng)器CQ-4 電路采用氣體開關(guān)，雖然目前簡(jiǎn)單的CKT 計(jì)算也可達(dá)到偏差低于10%的精度，但明顯不如CUR 方式。本文電路耦合CKT 計(jì)算的原理是正確的，但對(duì)具體電路的描述還有待改進(jìn)。

4.2 磁驅(qū)動(dòng)高速金屬飛片實(shí)驗(yàn)（算例5）

磁驅(qū)動(dòng)高速金屬飛片實(shí)驗(yàn)利用等熵壓縮的磁力斜波，以極高壓力（數(shù)百吉帕以至太帕）驅(qū)動(dòng)作為電極板一部分的金屬飛片高速射出，鋁、銅飛片的速度很容易達(dá)到10 km/s 以上[21]，2011 年美國(guó)Sandia 實(shí)驗(yàn)室已能把亞毫米厚度的鋁飛片發(fā)射到45 km/s[36]，開創(chuàng)了沖擊壓縮實(shí)驗(yàn)的新境界，對(duì)極端物理學(xué)研究具有重要意義。飛片運(yùn)動(dòng)前期由強(qiáng)磁壓驅(qū)動(dòng)，電流衰減的后期則由飛片后側(cè)電弧燒蝕等離子體推進(jìn)，這類電磁熱力強(qiáng)耦合的復(fù)雜實(shí)驗(yàn)必須采取多物理場(chǎng)耦合計(jì)算進(jìn)行設(shè)計(jì)或核算。

算例5 是Lemke 等[36]2011 年所做系列實(shí)驗(yàn)之“11 mm-2 s”的模擬，該系列實(shí)驗(yàn)研究了美國(guó)Sandia實(shí)驗(yàn)室Z 裝置驅(qū)動(dòng)高速金屬飛片的能力，給定了用于直至TPa 范圍的沖擊壓縮標(biāo)準(zhǔn)加載手段。實(shí)驗(yàn)“11 mm-2 s”是文獻(xiàn)[36]中唯一報(bào)道了加載電流波形的具體實(shí)驗(yàn)。該實(shí)驗(yàn)采用單條片構(gòu)形，材料為鋁，電極板寬度為11 mm、長(zhǎng)度為25 mm，構(gòu)形的盲孔為矩形凹槽，槽底即矩形金屬飛片，AK-gap 即磁腔間隙為1 mm。電極板凹槽頂面直接覆蓋7.5 mm 厚的銅靶板，凹槽深度即飛片的飛行距離約為7.8 mm。實(shí)驗(yàn)測(cè)量飛片后自由面速度波形，并提供其剖面參數(shù)及高速碰撞過程的信息。計(jì)算采用的等效電流因數(shù)K為文獻(xiàn)[36]給定值0.76（本文自算值為0.73）。圖9 為SSS-MHD 編碼的模擬結(jié)果及其比較。

圖9 對(duì)于Sandia 實(shí)驗(yàn)室Z 裝置鋁飛片實(shí)驗(yàn)“11 mm-2 s”的SSS-MHD 模擬（算例5）Fig. 9 SSS-MHD simulations for the Al flyer experiment “11 mm-2 s”on the Z machine at the Sandia National Laboratories, USA (example 5)

圖9(a)顯示SSS-MHD 計(jì)算的“11 mm-2 s”飛片速度與原實(shí)驗(yàn)及ALEGRA-2D 編碼模擬相當(dāng)一致，但在3.2 μs 后有些偏高。圖9(b)顯示在磁力推動(dòng)后期（3.05 μs）磁壓力高達(dá)285 GPa。圖9(c)表示飛行后期3.24 μs 時(shí)飛片自由面后固體區(qū)溫度過高（1 000 K 左右），導(dǎo)致密度略低于常值，這個(gè)偏差與圖9(a)后期計(jì)算速度偏高的情況相符。雖然這里沒有進(jìn)行飛片中熔化邊界的計(jì)算，但從圖9(c)中3 條剖面曲線斜率變化的拐點(diǎn)可以判斷，飛片后自由面固態(tài)區(qū)厚度約為0.1 mm，符合Z 裝置飛片實(shí)驗(yàn)通常的數(shù)據(jù)。圖9(b)是與圖9(c)時(shí)間不同的剖面，但與文獻(xiàn)[36]中另一實(shí)驗(yàn)“11 mm-1 s”剖面的數(shù)據(jù)相當(dāng)一致，提供了有依據(jù)的物理圖像。

5 SSS-MHD 編碼應(yīng)用于圓筒構(gòu)形電磁內(nèi)爆和MC-1 發(fā)生器實(shí)驗(yàn)的典型算例

本節(jié)計(jì)算的金屬套筒電磁內(nèi)爆實(shí)驗(yàn)和CJ 型MC-1 發(fā)生器實(shí)驗(yàn)，可使材料樣品中斜波壓力達(dá)到太帕量級(jí)，或可使構(gòu)形軸線處的壓縮磁場(chǎng)達(dá)到700～1 000 T（相當(dāng)于磁壓200～400 GPa），對(duì)于低密度物質(zhì)的壓縮具有特殊意義。

5.1 固體金屬套筒高速電磁內(nèi)爆實(shí)驗(yàn)（算例6）

磁驅(qū)動(dòng)圓筒構(gòu)形的內(nèi)爆實(shí)驗(yàn)可以高效實(shí)現(xiàn)樣品材料的斜波壓縮或套筒的高速度高對(duì)稱性內(nèi)爆運(yùn)動(dòng)，其負(fù)載構(gòu)形即是作為陰、陽電極的內(nèi)外兩層同軸金屬圓筒，其間為圓環(huán)磁腔。這種構(gòu)形的主體—作內(nèi)爆運(yùn)動(dòng)的內(nèi)筒稱為套筒（liner）；外筒稱為回流筒（也可由若干金屬柱組成），做速度不快的外向擴(kuò)展運(yùn)動(dòng)，對(duì)內(nèi)爆影響不大，然而它回流的電路功能不可或缺。由于套筒內(nèi)爆運(yùn)動(dòng)是會(huì)聚的，其電磁驅(qū)動(dòng)力和內(nèi)聚壓力越來越強(qiáng)，電流有效作用時(shí)間和效率都明顯高于平面構(gòu)形。金屬套筒內(nèi)爆速度有可能超過100 km/s，內(nèi)爆壓力可達(dá)到數(shù)太帕，成為高能量密度物理宏觀樣品實(shí)驗(yàn)的主要手段。

算例6 是對(duì)于美國(guó)Z 機(jī)器磁驅(qū)動(dòng)Al/Cu 復(fù)合套筒超高壓內(nèi)爆壓縮實(shí)驗(yàn)[34]的SSS-MHD 數(shù)值模擬。圖10(a)是該實(shí)驗(yàn)裝置的截面示意圖，圖中的陽極—回流筒（厚度0.45 mm、內(nèi)半徑13 mm）和陰極—內(nèi)套筒的驅(qū)動(dòng)層（厚度1 mm、外半徑3.43 mm）由6061-T6 鋁合金制造。外筒的膨脹速度歷史由裝置外圍排列的VISAR 探針測(cè)量。內(nèi)套筒的內(nèi)層—樣品層為銅材料，厚度0.53 mm、內(nèi)半徑1.9 mm（此系列還進(jìn)行了鉭和鋁樣品材料的內(nèi)爆壓縮實(shí)驗(yàn)）。套筒的有效高度為10 mm，內(nèi)窺PDV 探針由6 根PDV 光纖組成，封裝于壁厚0.15 mm、外半徑0.35 mm 的鉑管之中。圖10(b)顯示SSS-MHD 編碼的計(jì)算速度與實(shí)驗(yàn)及ALERGRA-1D 模擬結(jié)果相符。當(dāng)內(nèi)爆中止即銅套筒撞上鉑管時(shí)（約3.016 μs），波形優(yōu)化設(shè)計(jì)的加載電流達(dá)到約16.5 MA，此時(shí)鋁驅(qū)動(dòng)層外半徑約2 mm，電流線密度高達(dá)13 MA/cm，磁壓力超過1 TPa，并且壓縮過程是準(zhǔn)等熵的。圖10(c)顯示銅樣品中最高壓力為1.25 TPa，最高密度達(dá)22.4 g/cm3，可看出SSS-MHD編碼計(jì)算的剖面臺(tái)階與ALEGRA-2D 模擬結(jié)果略有差別，主要原因是介質(zhì)密度接觸間斷的處理存在缺陷，有可能是本文計(jì)算中空間格子尺度偏大等原因造成的，需要改進(jìn)。

圖10 對(duì)于Sandia 實(shí)驗(yàn)室Z 裝置上金屬套筒電磁內(nèi)爆實(shí)驗(yàn)的SSS-MHD 模擬（算例6）Fig. 10 SSS-MHD simulation of magnetically-driven liner implosions on the Z-machine at the Sandia National Laboratories, USA (example 6)

電磁內(nèi)爆壓縮實(shí)驗(yàn)的意義在于獲得太帕等級(jí)的物態(tài)方程數(shù)據(jù)，關(guān)鍵在于如何較準(zhǔn)確地處理實(shí)驗(yàn)數(shù)據(jù)。文獻(xiàn)[34]中經(jīng)過精確的物態(tài)方程計(jì)算，得到圖10(c)中銅樣品的最高壓力、密度標(biāo)定值的偏差均小于1%，是否意味著已具有狀態(tài)方程參考價(jià)值，這些重要問題需要繼續(xù)從實(shí)驗(yàn)和模擬方面進(jìn)行探討。

5.2 圓柱形炸藥內(nèi)爆磁通量壓縮實(shí)驗(yàn)（算例7）

MC-1 發(fā)生器利用圓柱形炸藥驅(qū)動(dòng)金屬套筒內(nèi)爆，把套筒空腔內(nèi)種子磁場(chǎng)的初始磁通量壓縮到軸線周圍小區(qū)域中，從而得到高磁場(chǎng)和高磁壓力[37]。MC-1 發(fā)生器與5.1 節(jié)套筒電磁內(nèi)爆的基本不同點(diǎn)在于MC-1 套筒是圖3 中的圓柱形（N=1）θ 構(gòu)形，而電磁內(nèi)爆套筒屬于z構(gòu)形。MC-1 套筒的環(huán)形“勵(lì)磁”電流物理上是注入套筒腔內(nèi)的種子磁場(chǎng)Bz0在內(nèi)爆套筒中感生的渦電流。從圓柱電磁構(gòu)形的分析可知，z構(gòu)形產(chǎn)生的環(huán)形磁場(chǎng)驅(qū)動(dòng)套筒內(nèi)爆具有較高效率，θ 構(gòu)形壓縮腔內(nèi)軸向磁通量則有較高效率。

CJ-100 型MC-1 裝置是中國(guó)工程物理研究院流體物理所研制的小型MC-1 發(fā)生器[38-39]，圖11 是該裝置的截面示意圖和實(shí)物照片。CJ-100 型裝置的套筒材料為304 不銹鋼，外直徑100 mm、厚度1.5 mm、高度210 mm。炸藥為RHT-901（m(TNT)/m(RDX)=40/60），尺寸為內(nèi)直徑100 mm、外直徑210 mm、高度65 mm，總炸藥量不到3 kg。

圖11 CJ-100 型MC-1 發(fā)生器的示意圖和實(shí)物照片F(xiàn)ig. 11 Schematics and picture of the CJ-100 type MC-1 generator

算例7 是CJ-100 型MC-1 裝置磁場(chǎng)壓縮實(shí)驗(yàn)Shot20150630[39]的SSS-MHD 驗(yàn)算。該實(shí)驗(yàn)套筒內(nèi)部的種子磁場(chǎng)Bz0約為5.5 T。該實(shí)驗(yàn)中軸線附近磁場(chǎng)增長(zhǎng)歷史以及套筒內(nèi)、外半徑變化的模擬計(jì)算示于圖12(a)，在磁探針停止工作之前實(shí)測(cè)磁場(chǎng)與計(jì)算值十分符合。

圖12 CJ-100 型MC-1 發(fā)生器磁通量壓縮實(shí)驗(yàn) 及其性能的SSS-MHD 模擬計(jì)算（算例7）Fig. 12 SSS-MHD simulations for the magnetic flux compression experiment and the performances of the CJ-100 type MC-1 generator (example 7)

計(jì)算的套筒內(nèi)層回彈半徑約2.6 mm，此時(shí)計(jì)算的峰值磁場(chǎng)約為1 450 T。但實(shí)測(cè)磁場(chǎng)峰值為690 T，從圖12(a)看出此時(shí)套筒腔內(nèi)磁場(chǎng)區(qū)的半徑值約3.8 mm，此值相當(dāng)于實(shí)驗(yàn)中磁場(chǎng)測(cè)點(diǎn)的半徑值。圖12(b)是不同種子磁場(chǎng)Bz0之下SSS-MHD 計(jì)算的CJ-100 型發(fā)生器達(dá)到的峰值磁場(chǎng)和套筒回彈半徑，可作為實(shí)驗(yàn)設(shè)計(jì)的參考。圖12(b)中還畫出了該型號(hào)發(fā)生器以前2 發(fā)實(shí)驗(yàn)EX1、EX2 的壓縮磁場(chǎng)，以及對(duì)相應(yīng)磁場(chǎng)測(cè)點(diǎn)半徑值的估計(jì)，數(shù)據(jù)與本例相近。從圖12(b)可知，設(shè)計(jì)適當(dāng)?shù)姆N子磁場(chǎng)值，可得到折中的壓縮磁場(chǎng)值和可利用磁場(chǎng)的空間范圍。

6 結(jié) 束 語

SSS-MHD 編碼是拉格朗日形式一維平面和圓柱面幾何、多介質(zhì)、多組分、多連通區(qū)的彈塑性磁流體力學(xué)多物理場(chǎng)計(jì)算編碼，是一維沖擊爆轟編碼SSS 向磁流體力學(xué)擴(kuò)展而形成，具有25 個(gè)子程序和函數(shù)，共約6 000 條Fortran 語句。SSS-MHD 編碼具有框架式結(jié)構(gòu)，便于靈活增加和擴(kuò)充介質(zhì)類型、物性數(shù)據(jù)及物理功能。目前，除了SSS 編碼原有的功能外，SSS-MHD 編碼主要面向強(qiáng)電磁驅(qū)動(dòng)的高能量密度動(dòng)力學(xué)實(shí)驗(yàn)，特別是為極端材料動(dòng)力學(xué)實(shí)驗(yàn)提供數(shù)值模擬平臺(tái)。

本文中，第1～3 節(jié)闡述了SSS-MHD 編碼的理論基礎(chǔ)、程序結(jié)構(gòu)和特色，第4～5 節(jié)給出了該編碼模擬強(qiáng)電磁驅(qū)動(dòng)實(shí)驗(yàn)的7 個(gè)算例，包括平面構(gòu)形材料磁驅(qū)動(dòng)準(zhǔn)等熵壓縮實(shí)驗(yàn)、磁驅(qū)動(dòng)發(fā)射高速金屬飛片實(shí)驗(yàn)、套筒電磁內(nèi)爆材料壓縮實(shí)驗(yàn)和MC-1 發(fā)生器磁通量壓縮實(shí)驗(yàn)?？傮w看來，SSS-MHD 模擬與實(shí)驗(yàn)及其他方法數(shù)值模擬結(jié)果符合較好，相對(duì)偏差均在5%以下。這些算例的實(shí)驗(yàn)形式廣泛多樣，物理變量范圍寬廣，實(shí)驗(yàn)測(cè)量和模擬計(jì)算精度較高，為考核SSS-MHD 編碼的實(shí)際能力提供了合適的試題，同時(shí)說明該編碼已達(dá)到強(qiáng)電磁驅(qū)動(dòng)實(shí)驗(yàn)一維數(shù)值模擬平臺(tái)的要求，可以推廣應(yīng)用。

通過與該領(lǐng)域先進(jìn)實(shí)驗(yàn)及數(shù)值模擬結(jié)果的核對(duì)，為SSS-MHD 編碼的繼續(xù)改進(jìn)和擴(kuò)展指出了方向。例如，補(bǔ)充流體力學(xué)模型，適應(yīng)內(nèi)爆動(dòng)力學(xué)實(shí)驗(yàn)?zāi)M的需要；在基本方程組層面引入更普遍的非線性連續(xù)介質(zhì)力學(xué)本構(gòu)理論，以適應(yīng)大變形高應(yīng)變率計(jì)算的需要；改進(jìn)單溫等離子體為多溫多流模型，向簡(jiǎn)單輻射磁流體力學(xué)計(jì)算發(fā)展；建立更精細(xì)準(zhǔn)確的驅(qū)動(dòng)器電路和關(guān)鍵器件模型，實(shí)現(xiàn)更準(zhǔn)確的MHD 與驅(qū)動(dòng)電路的耦合計(jì)算；改進(jìn)高功率激光與物質(zhì)耦合計(jì)算模型，建立與SOP (streaked opticapyrometer)等先進(jìn)測(cè)溫技術(shù)配合的功能，進(jìn)一步擴(kuò)大對(duì)高能量密度物理實(shí)驗(yàn)的模擬能力。

本文工作得到了流體物理研究所諸多磁驅(qū)動(dòng)、磁壓縮實(shí)驗(yàn)數(shù)據(jù)的支持，謹(jǐn)向參加有關(guān)實(shí)驗(yàn)的仝延錦、唐小松、宋振飛、程誠(chéng)、李建明、匡學(xué)武、吳剛、稅榮杰、胥超、鄧順益、馬驍?shù)韧局乱灾孕母兄x！