電磁驅動高能量密度動力學實驗的一維磁流體力學多物理場數值模擬平臺：SSS-MHD*

孫承緯，陸 禹，趙繼波，羅斌強，谷卓偉，王桂吉，張旭平，陳學秒，周中玉，李 牧，袁 紅，張紅平，王剛華，孫奇志，文尚剛，譚福利，趙劍衡，莫建軍，蔡進濤，金云聲，賀 佳，種 濤，趙小明，劉倉理

（1. 中國工程物理研究院流體物理研究所，四川 綿陽 621999；2. 上海激光等離子體研究所，上海 201800；3. 深圳技術大學先進材料測試技術研究中心，廣東 深圳 518118；4. 深圳技術大學大數據與互聯網學院，廣東 深圳 518118；5. 中國工程物理研究院化工材料研究所，四川 綿陽 621999；6. 中國工程物理研究院應用電子學研究所，四川 綿陽 621999；7. 中國工程物理研究院，四川 綿陽 621999）

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

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

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

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

圖1 適合SSS-MHD 編碼模擬的各種類型電磁驅動高能量密度物理實驗Fig. 1 Kinds of the magnetically-driven high-energy-density physics experiments suitably simulated with the SSS-MHD code

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

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

1 彈塑性磁流體力學基本方程組

連續介質力學運動的基本規律為質量、動量和能量守恒定律，從一般的連續介質動力學方程組可以方便地導出磁流體力學的基本方程組[26]。首先，在動量方程和能量方程中分別增加與電磁力和電阻產生熱量（焦耳熱）相關的項，把向量形式的連續介質動力學方程組改寫如下，依次為連續性方程（質量守恒）、運動方程（動量守恒）和能量方程（能量守恒）：

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

式(4)是SSS-MHD 編碼實際使用的能量方程，其優點是簡化了張量運算。后文中，能量方程中無外體力F項，外源項Q寫成WL即單位質量介質吸收的激光輻射功率。

在沖擊動力學或爆炸力學中，對連續介質通常采用流體彈塑性的力學模型加以描述。將應力張量Σ分解為 Σ=-pI+T，寫成分量式：

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

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

法拉第定律：

安培定律：

再增加廣義歐姆定律：

前面方程組已含有粒子速度u為未知函數，式(7)～(9)提供了求解E、J和B等3 個電磁量的方程，其中 μ 為介質的磁導率。由于安培定律，式(1)～(3)中僅出現磁場B，在該意義下系統被稱為磁流體力學方程組。進一步化簡，可從式(7)～(9)中消去E和J，得到決定B的磁擴散方程：

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

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

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

運動方程：

能量方程：

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

磁擴散方程組：

式中：Cθ=rNBθ稱為磁感應強度的廣義 θ 分量，在柱對稱情形中有重要意義。

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

為了加以區分，將拉格朗日形式方程組中的r改寫為R，粒子速度u改寫為U，徑（縱）向應力 σrr簡寫為 σ 。每個拉格朗日格點的質量保持不變體現了質量守恒，也省去了連續性方程。為此，需要添加一個速度方程，即U等于R的時間偏導數。

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

經過上述推導和轉換，得到SSS-MHD 編碼使用的拉格朗日形式一維彈塑性磁流體力學方程組如下。

連續性方程和速度方程：

運動方程：

能量方程：

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

磁擴散方程組：

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

3 SSS-MHD 編碼的結構和功能

3.1 計算編碼的結構

SSS-MHD 編碼的結構立足于SSS 程序原有的框架[7]，保持彈塑性流體動力學計算的各個子程序模塊構架，進行適應MHD 計算的改造。首先，在導程序中增加MHD 建模的內容，在總卡（GENERL）和區卡（COMPNT）中引入與MHD 計算相關的數組、變量和指數等，然后，在專門的MHD 建模子程序（MHDZON）里再對有關區卡賦以MHD 負載構形的類型指數、構形參數以及初始條件等，完成MHD 計算建模。接著，在相關子程內部和計算循環節點處加入控制MHD 計算的語句，如在計算循環開始處設定邊界條件狀態及參數的同時，檢查和調整空腔和磁腔的狀態及參數。

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

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

原編碼SSS 是一個框架式的通用程序，具有多種物理場耦合計算的功能，可以包容或擴展各種物態方程和物性數據庫，加入有關局部性物理過程的計算（例如激光與物質的熱和沖量耦合、平衡電離的Saha 方程、動態斷裂的NAG 模型和含能材料的燃燒爆炸反應等），適合平、柱、球對稱幾何的動力學問題，設有多種邊界條件和初始條件可供選擇，能方便地實現沖擊或爆轟的加載模擬[7]，這些特點已全部為SSS-MHD 編碼沿用。下面各小節將具體闡述SSS-MHD 計算編碼的關鍵內容和特色。

3.2 MHD 負載構形的統一建模

在磁流體力學計算模型的建模中，由于力學與電磁量具有不同方向的坐標分量，首先需要設定力學構形和相應的加載磁場類型，然后賦予各個MHD 計算區域及格點處的電磁學物性參數和相應初始值，必要時在整體模型邊界處補充電磁學邊界條件。其中的實質性工作是負載構形及磁場的分類和表述。

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

圖3 SSS-MHD 編碼中實驗負載構形的種類（紅線表示回流導體）[29]Fig. 3 Types of experimental configurations in the SSS-MHD code(red curves standing for return conductors)[29]

圖1 中各類實驗的構形和磁場均可抽象為圖3 的類型。負載構形中的磁場可以是加載電流感生的（持續勵磁），也可以是初始時由外界給定的（瞬態種子磁場或定常勵磁）；推動構形運動的可以是加載電流自身形成的洛倫茲力（磁驅動），也可以利用炸藥爆轟或沖擊的力量（磁壓縮）。

3.3 空腔和磁腔

SSS-MHD 編碼沿用了SSS 編碼的空腔功能，并擴展到帶磁場的磁腔情形。計算模型通常為順序排列的介質分區所組成的單連通區，如果隊列中存在真空間隔并不妨礙歐拉形式計算，但對于拉格朗日計算則發生“零質量格點”的困難，通常需采取分區分階段計算合成或實行一維ALE 策略，這樣可能使得編碼復雜、機時增多，計算精度也受影響。

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

圖4 SSS 編碼中空腔區及鄰區中格點的設置Fig. 4 Meshes in the cavity and nearby in the code SSS

磁流體力學計算中磁空腔較常見，其力學運動可沿用SSS 編碼的處理方法，但需增加電磁學計算，并保持兩岸點磁場及磁力對稱相等的處理。經過改進，SSS-MHD 編碼實現了一維磁流體力學多連通區的拉格朗日一體化計算，具備了對于較復雜實驗及裝置進行數值模擬的能力。

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

3.4 磁腔磁場的計算

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

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

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

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

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

本節中論述的平面磁腔磁場近似計算方法對于SSS-MHD 編碼的磁擴散計算具有重要意義，因為磁腔中的磁場Bz或Cθ都是常值，根據式(24)～(26)可由構形加載總電流的即時值直接確定。該總電流（時間函數）可以作為輸入數據給定，也可從電路方程組計算獲得。因而，這些磁腔的磁場值可以用作為磁擴散方程計算中可靠的定標值，避免了隱式計算的麻煩和誤差。該問題在理論上有更嚴謹的方法，如Lemke 等[32]利用Poisson 方程計算了含有四邊形環形磁腔的多連通區平面二維磁流體力學問題，但并不能得到有普遍意義的結果。

3.5 顯式計算格式

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

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

3.6 SSS-MHD 計算與驅動電路方程的耦合

負載構形的動態電感變化與驅動器放電過程之間存在重要的相互作用，若要理解這些復雜過程，首先應把MHD 計算與驅動器電路計算進行耦合，建立準確預測實際加載電流的能力。與驅動電路的耦合是MHD 計算能否解決實際問題的一個關鍵，也是SSS-MHD 編碼的重要特色。這個要求至今仍是對強電磁驅動實驗計算模擬能力的重大挑戰。

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

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

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

3.7 SSS-MHD 計算程序包MHDBLK

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

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

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

4 SSS-MHD 編碼應用于平面構形磁驅動實驗的典型算例

SSS-MHD 編碼針對各類實驗研究以及國外高水平實驗結果，進行了大量驗算和改進工作。本節列舉的平面構形磁驅動實驗（ICE）算例分為2 類：材料樣品準等熵壓縮和高速金屬飛片發射。

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

圖6 磁驅動準等熵壓縮實驗的典型負載構形Fig. 6 Typical loading configurations for magnetically-driven isentropic compression experiments (ICE)

4.1 磁驅動平面材料樣品準等熵壓縮實驗（算例1～4）

算例1 是對美國Sandia 實驗室Z 裝置實驗Z-1220[33]的驗算，該實驗采用類似圖6(c)的6061-T6 鋁合金4 層方筒構形，將一個側面（陽電極板）加工為4 段樣品，厚度h分別為0.616、0.868、1.118 和1.371 mm，相當于一個4 階的臺階靶。每段樣品外表面粘有鍍增透膜的LiF 光學窗口。樣品與窗口的界面速度用VISAR 測量。實驗加載電流峰值約18 MA，上升沿約300 ns。美國Z-1220 實驗數據以及ALEGRA-1D 編碼與SSS-MHD 編碼的計算模擬都示于圖7。

圖7 對于美國Sandia 實驗室磁驅動等熵壓縮實驗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 計算是依據圖7(a)的加載電流實驗波形進行的，此例的矩形環腔周長5.8 cm，換算得知平均電流密度峰值約3.1 MA/cm，等效電流因數K為1。4 種厚度樣品的2 種數值模擬速度曲線與實驗數據的符合程度都較好。SSS-MHD 計算給出的材料樣品最高壓力為51 GPa，文獻[33]給出的最高縱向應力為54 GPa。2 種模擬的差別在于SSS-MHD 計算在彈塑性轉變階段速度起跳過早，其原因可能是該計算的樣品材料是純鋁而不是Z 實驗的6061-T6 鋁合金，此外還需要考察計算中輸入電流波形初始50 ns 段落的細節以及鋁樣品中磁場擴散的情況，進而或可做出相應改進。對于模擬速度曲線末期與實驗數據的偏差，主要是因為鋁的彈塑性行為與模擬選用的材料動態本構模型之間存在差異。

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

表1 磁驅動準等熵壓縮實驗算例2～4 的主要參數[14,16]Table 1 Parameters of examples 2－4 for magnetically-driven isentropic compression experiments[14,16]

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

圖8 對于流體物理研究所CQ 裝置磁驅動準等熵壓縮實驗的SSS-MHD 模擬（算例2～4）Fig. 8 SSS-MHD simulations for magnetically-driven quasi-isentropic compression experiments at IFP (examples 2－4)

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

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

4.2 磁驅動高速金屬飛片實驗（算例5）

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

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

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

5 SSS-MHD 編碼應用于圓筒構形電磁內爆和MC-1 發生器實驗的典型算例

本節計算的金屬套筒電磁內爆實驗和CJ 型MC-1 發生器實驗，可使材料樣品中斜波壓力達到太帕量級，或可使構形軸線處的壓縮磁場達到700～1 000 T（相當于磁壓200～400 GPa），對于低密度物質的壓縮具有特殊意義。

5.1 固體金屬套筒高速電磁內爆實驗（算例6）

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

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

圖10 對于Sandia 實驗室Z 裝置上金屬套筒電磁內爆實驗的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)

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

5.2 圓柱形炸藥內爆磁通量壓縮實驗（算例7）

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

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

圖11 CJ-100 型MC-1 發生器的示意圖和實物照片Fig. 11 Schematics and picture of the CJ-100 type MC-1 generator

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

圖12 CJ-100 型MC-1 發生器磁通量壓縮實驗 及其性能的SSS-MHD 模擬計算（算例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)

計算的套筒內層回彈半徑約2.6 mm，此時計算的峰值磁場約為1 450 T。但實測磁場峰值為690 T，從圖12(a)看出此時套筒腔內磁場區的半徑值約3.8 mm，此值相當于實驗中磁場測點的半徑值。圖12(b)是不同種子磁場Bz0之下SSS-MHD 計算的CJ-100 型發生器達到的峰值磁場和套筒回彈半徑，可作為實驗設計的參考。圖12(b)中還畫出了該型號發生器以前2 發實驗EX1、EX2 的壓縮磁場，以及對相應磁場測點半徑值的估計，數據與本例相近。從圖12(b)可知，設計適當的種子磁場值，可得到折中的壓縮磁場值和可利用磁場的空間范圍。

6 結 束 語

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

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

通過與該領域先進實驗及數值模擬結果的核對，為SSS-MHD 編碼的繼續改進和擴展指出了方向。例如，補充流體力學模型，適應內爆動力學實驗模擬的需要；在基本方程組層面引入更普遍的非線性連續介質力學本構理論，以適應大變形高應變率計算的需要；改進單溫等離子體為多溫多流模型，向簡單輻射磁流體力學計算發展；建立更精細準確的驅動器電路和關鍵器件模型，實現更準確的MHD 與驅動電路的耦合計算；改進高功率激光與物質耦合計算模型，建立與SOP (streaked opticapyrometer)等先進測溫技術配合的功能，進一步擴大對高能量密度物理實驗的模擬能力。

本文工作得到了流體物理研究所諸多磁驅動、磁壓縮實驗數據的支持，謹向參加有關實驗的仝延錦、唐小松、宋振飛、程誠、李建明、匡學武、吳剛、稅榮杰、胥超、鄧順益、馬驍等同志致以衷心感謝！