電爆炸金屬橋箔早期過程中電磁-熱-力多物理場耦合建模與分析

2024-02-24 09:54:44王剛華唐久棚歐海彬闞明先段書超

高壓物理學報 2024年1期

王剛華，謝龍，肖波，王強，唐久棚，歐海彬，闞明先，段書超

（中國工程物理研究院流體物理研究所, 四川綿陽 621999）

利用大電流流經金屬箔/金屬絲時的歐姆加熱作用，可以將金屬箔/絲汽化電離，產生急速膨脹的高能量密度等離子體，進而可以作為高壓驅動源。該類高壓驅動源應用廣泛，例如：基于該原理發展的電炮可以用于開展超高速飛片發射和材料動力學特性研究[1-4]；利用金屬箔電爆炸發射小尺寸飛片來制造高安全性、低感度的沖擊起爆雷管[5-10]等。

電爆炸過程中，導體被快速加熱，短時間內經歷固-液-氣/等離子體多個狀態，歐姆電阻可能增長幾個數量級。因此，采用數值模擬方法研究這一過程并非易事。早期的研究主要基于電路方程和電導率模型以實現對實驗電流的模擬和溫度的估算[11-12]。1981 年，董玉斌等[12]將金屬箔電爆炸過程分為沸點前和沸點后兩個階段，分階段計算了放電電流和橋箔電阻的變化以及金屬箔電爆炸驅動飛片的加速過程，實現了一維流體動力學方程、放電回路方程以及沙哈方程的耦合求解。1989 年，Osher 等[13]報道了用于計算金屬箔電爆炸全過程的一維磁流體動力學程序，該程序與2010 年Saxena 等[14]報道的基于Burgess 電導率模型的一維磁流體程序SEEFA 一樣，均對金屬箔區域做零維處理。后來，羅斌強[15]、賀佳等[16]和Luo 等[17]結合大量實驗完善了電箔爆炸一維磁流體程序，該程序能較好地模擬爆炸箔的放電波形，同時給出與實驗基本一致的飛片速度波形。

然而，對于沖擊片雷管而言，提升雷管工作的可靠性需要對橋箔電爆炸驅動飛片的過程有更為深入的認識，一維簡化的磁流體模型雖然能夠預估飛片的速度，但不能給出飛片的形狀、溫度、速度的空間分布這些可靠起爆炸藥的關鍵信息，也無法給出爆炸箔上電流密度、磁場的空間分布，因而無法為具體的改進設計提供參考。因此，發展更全面的物理模型和更高維度的模擬程序是深入理解沖擊片雷管爆炸物理過程的基礎。2017 年，英國原子武器研究中心（Atomic Weapons Establishment，AWE）和美國圣地亞實驗室（Sandia National Laboratories，SNL）報道了他們利用三維磁流體模擬程序ALEGRA MHD 開展的爆炸箔起爆的數值仿真工作[18]，首次獲得了與實驗照相結果基本相似的圖像，是目前報道的最先進的沖擊片雷管模擬能力。

近年來，我國的爆炸箔數值仿真工作取得了不少進展，在工程應用中發揮了較大作用[19-23]。然而，目前國內相關數值模擬工作還存在物理模型不全、空間維度降階等問題，對于最為關鍵的橋箔三維結構和電磁場分布造成的等離子體空間分布等問題缺乏研究模型和研究工具。例如，北京理工大學與中國工程物理研究院化工材料研究所合作開展了三維沖擊片雷管的數值模擬研究，對橋箔進行了簡化處理，使用的是等效電阻模型，無法分析電磁場和電流的三維空間分布。雖然基于磁流體的數值模型在物理上較為全面，但是現有的報道還是一維的，中國工程物理研究院流體物理研究所在二維磁流體方面的研究正在進行，發展二維、三維磁流體數值模擬能力尚需時日。

盡管橋箔電爆炸的全過程模擬極其困難，但是其后續復雜的物理過程是在早期能量沉積和分布的基礎上演化而來的，因此，在大量等離子體形成和爆炸膨脹之前，弄清楚橋箔上的電流分布和溫度分布演化過程極其重要，有利于開展橋箔性能改進設計。本研究綜合考慮電磁場在三維結構上的擴散、歐姆加熱效應、材料受熱膨脹效應，擬建立三維電爆炸橋箔的電磁-熱-固體力學耦合模型，開展橋箔早期電磁能量沉積過程模擬。這一模型充分考慮爆炸橋箔發生劇烈汽化膨脹前的各種物理效應，是一種宏觀尺度上的全三維、全物理模擬方法，對清晰理解爆炸橋箔早期過程和指導工程設計具有重要意義。

1 電磁-熱-固體力學耦合模型

橋箔的早期能量沉積過程取決于電流密度和電導率，而電導率又依賴于溫度場分布。由于早期溫度相對較低，橋箔膨脹較小，可以使用固體力學來描述橋箔的膨脹過程。為此，建立了電磁-熱-固體力學耦合的物理模型以模擬橋箔早期的能量沉積過程。采用電磁模塊進行橋箔的電磁和電流演化分析，采用固體傳熱模塊考慮橋箔的溫度變化，采用固體力學模塊模擬橋箔的早期受熱膨脹情況。

電磁模塊求解如下控制方程。

對于橋箔區域，使用A形式的矢量勢公式（安培定律）描述

對于空氣區域，使用標量勢公式（磁通量守恒）描述

式中：H為磁場強度，A/m；B為磁通密度，T；E為電場，V/m；A為磁矢勢，Wb/m； σ為電導率，S/m； μ0為真空磁導率，H/m； μr為相對磁導率；Vm為磁標勢，A。

固體傳熱模塊求解如下控制方程

式中： ρ為密度，kg/m3；cp為比定壓熱容，J/(kg·K)；T為溫度，K；u為速度，m/s；k為熱導率，W/(m·K)；Qe為熱源，W/m3，耦合計算中熱源為歐姆加熱，即Qe=J·E。

固體力學模塊求解如下控制方程

式中：s為位移，m；S為應力張量，Pa；Fv為外力，N。

耦合計算中調用熱膨脹模塊，即考慮由熱效應引起的熱應變

式中： α(T) 為隨溫度變化的熱膨脹系數，K-1； εth為彈性應變；C為彈性模量張量，Pa；Tref為參考溫度，K。

計算初始時刻，電場、磁場強度、電流密度等均為零，溫度為常溫。電磁邊界條件主要包括：對稱面上指定磁標勢V0，區域邊界上法相磁場為零，空氣與橋箔交界面上切向磁場連續。固體傳熱模塊邊界條件：將所有面均設定為絕熱條件。力學邊界條件主要包括：對稱面上位移為零，其余面為自由面。計算中將圖1 所示的實驗電流作為激勵電流，對應的回路電阻為280 mΩ，回路電感為310 nH，電容為0.25 μF，充電電壓為5.1 kV。橋箔材料為銅，考慮了銅材料的熱導率（k）隨溫度（T）和密度（ ρ）的變化以及電導率隨溫度和密度的變化情況（如圖2 和圖3 所示）。

圖1 實驗獲得的電流歷史曲線Fig.1 Experimental current history curve

圖2 銅的熱導率Fig.2 Thermal conductivity of copper

圖3 銅的電導率Fig.3 Conductivity of copper

具體求解過程中，3 個模塊的耦合關系如下：電磁求解模塊通過電流密度J將電流分布傳遞給熱傳導模塊，傳熱模塊通過溫度T與固體力學模塊進行耦合；同時，材料溫度的改變通過電導率模型傳遞給電磁模塊，用于計算實時的電磁場分布；固體力學計算產生的幾何變形，也會耦合傳遞到電磁模塊和傳熱模塊。

為了減少計算量，采用如圖4 所示的1/2 幾何模型進行分析，并通過鏡像還原整體情況（圖4 中未顯示外圍空氣區）。由于物理模型只適用于小變形，因此計算時使用拉氏網格，網格劃分采用自適應四面體網格，在所關心的橋箔區采用密網格，空氣區使用較粗網格，總網格數約為2.084×107，如圖5 所示。圖4 中，由于橋箔的厚度（z方向）很薄，僅為5 μm（為了便于顯示，圖4 中將厚度放大了10 倍），遠小于x、y方向上的尺寸（毫米量級），為了減小計算量，劃分網格時，銅箔區在z方向只劃分了2 層網格。為便于后續沿標記線進行物理量分布分析，在如圖4 所示的橋面上標記兩條線L1、L2，其中：L1 為橋面上對稱軸線，L2 為邊緣輪廓線，弧長為l。

圖4 爆炸箔幾何1/2 模型示意圖（單位：mm）Fig.4 Schematic diagram of the 1/2 model of explosion foil (Unit: mm)

圖5 計算使用的網格Fig.5 Grid diagram used for calculation

2 計算分析

2.1 基于常數電導率的簡要分析

采用常溫電導率模型分析金屬橋箔的磁場和電流密度分布，不考慮材料電導率、熱導率隨溫度和密度的變化。由于常溫電導率模型只能用于分析放電早期的電磁場，為了更接近實際情況，需要進一步考慮歐姆加熱的影響，因此，只對常數電導率模型的計算結果作簡要介紹。

鑒于材料性質不發生變化，可采用歸一化電流進行分析。圖6 給出了橋箔區附近的磁場強度分布。受幾何構型的影響，磁場強度在橋區和過渡區的交界處表現出明顯的非均勻性，特別是邊緣附近非均勻性更為顯著。因此，邊緣附近的電流密度也很大，如圖7 所示。

圖6 表面磁場強度分布Fig.6 Magnetic field intensity distribution of upper and lower surfaces

圖7 電流密度與方向分布Fig.7 Current density and directional distribution

依據圖7 的電流密度分布，繪制電流密度的高度圖，如圖8 所示。圖9 給出了沿L1、L2 的電流密度隨弧長l的分布，弧長起點為橋箔左端。從圖8 和圖9 可以非常直觀地看出：隨著過渡區變窄，電流密度逐步增加，在橋區和過渡區交界處的4 個角點上電流最為集中；在橋區主體部分雖然還存在一定的邊界效應，但橋面主體部分的電流分布較為均勻。

圖8 電流密度高度Fig.8 Current density height

圖9 線上電流密度分布Fig.9 Current density distribution on lines

爆炸箔上的電流分布決定了各處產生的電功率密度Pe（圖10），也隨即決定了各處的溫升和早期等離子體的分布，因為在電功率高的地方，材料被率先加熱，進而汽化形成等離子體。

圖10 電功率密度高度Fig.10 Electric power density height

本節的所有分析是基于常數電導率模型開展的，即爆炸箔的電導率在空間上均勻分布，因此，電流和磁場只與爆炸箔的幾何尺寸相關，各個時刻獲得的磁場強度、電流密度的空間分布類似，不同時刻電流的大小僅影響其絕對值。然而，實際上，金屬的電導率與其自身的溫度和密度密切相關，需要考慮這種因素的影響。

2.2 基于Burgess 電導率的計算分析

常數電導率模型只能分析放電早期的電磁場，在中后期，為了更接近實際情況，需要進一步考慮歐姆加熱的影響。事實上，隨著電流加熱形成等離子體，材料的電參數會發生復雜的變化，為此，在電導率模型中引入隨溫度變化的Burgess 電導率模型，并引入固體傳熱模塊和固體力學模塊，以分析橋箔放電一段時間內的電磁場分布情況。

圖11 給出了引入Burgess 電導率模型、固體傳熱模塊和固體力學模塊后計算得到的不同時刻爆炸箔上的磁場強度分布。可以看出，與常數電導率模型類似，改進后的模型計算結果顯示，磁場先集中分布在爆炸箔4 個角點和爆炸箔邊緣，這些地方的磁場梯度較大，磁場分布極不均勻。在橋區，由于寬度變窄，磁場相對較大，邊緣效應更加明顯，特別是橋區與過渡區銜接的4 個角點的磁場強度很大。與常數電導率模型計算結果不同的是，隨著時間的增長，磁場向爆炸箔中間擴散，中間區域的磁場逐步增大。磁場擴散呈現兩種擴散特征：一是從4 個角點出發的點擴散，二是從邊緣向里的線擴散。

圖11 不同時刻的磁場強度分布Fig.11 Magnetic field intensity distribution at different instants

電流分布主要由兩個因素決定：構型和電導率。前述的磁場角點和邊緣效應是構型引起的。由于電流密度對應磁場強度的旋度，因此，磁場不均勻的地方電流密度更大。圖12 給出了由改進后的模型計算獲得的不同時刻的電流密度分布。可以看出：電流一開始集中在邊緣，然后向中間擴散，此時構型是電流分布的主導因素，在70 ns 時，橋區主體的電流分布基本均勻；此后，電導率因素占據主導，邊緣由于溫升快，其電導率低于中心處的電導率，使得電流向中心集中；后續電流分布在兩種因素的競爭中不斷演化。圖13 給出了90 ns 時刻沿線L1、L2 的電流密度分布，可以看出，其大致特征與常數電導模型率給出的結果一致，只是橋區的電流分布更均勻。

圖12 不同時刻的電流密度分布Fig.12 Current density distribution at different instants

圖13 90 ns 時刻沿線的電流密度分布Fig.13 Current density distribution on lines at 90 ns

爆炸箔的主要溫升機制是歐姆加熱，溫度分布取決于電流密度、電導率和時間。圖14 給出了溫度分布的演化。可以看出，早期邊緣加熱效應明顯，4 個角點的溫升相比于其他地方更快，溫度也更高。但是，隨著磁場和電流的擴散作用，電流向溫度較低的中間區域集中，橋區主體部分的溫度分布逐漸變得勻滑。圖15 給出了沿線L1、L2 的溫度分布。兩條曲線的平臺區對應橋區，兩條曲線的平臺有錯位，這是L1、L2 的弧長不一致導致的。兩條曲線的平臺有高度差，反映了橋區邊緣與橋區中心有溫度差。92 ns 時，橋區大部分區域的溫度均在1 245～1 345 K之間，溫度差在100 K 以內，這意味著在該結構設計下爆炸箔的溫度分布已十分均勻。通常情況下，在設計沖擊片雷管時，將橋箔爆炸時刻設計在電流峰值附近，雖然計算時刻與電流峰值時刻相距時間較長，尚不能確定后續溫度分布是否更加均勻，但是從目前的計算結果來看，溫度均勻性已經較好。此外，爆炸箔緊鄰橋區的過渡區中有相當多的區域溫升較為明顯，意味著該處沉積了較多能量。在沖擊片雷管設計中，可以充分利用該處產生的等離子體能量來驅動飛片，使飛片獲得更多動能。因此，從提高能量利用率的角度來講，設計沖擊片雷管時可考慮使炮筒口徑略大于橋區。

圖14 不同時刻的溫度分布Fig.14 Temperature distribution at different times

圖15 92 ns 時刻線上的溫度分布Fig.15 Temperature distribution on lines at 92 ns

圖16 給出了90 ns 時爆炸箔沿z軸的膨脹速度分布，4 個角點附近的膨脹速度相對較快，橋區邊緣的膨脹速度比橋區中心慢，其原因是邊緣處金屬箔除了z向膨脹以外，還向寬度方向膨脹。

圖16 90 ns 時刻沿z 軸的速度分布Fig.16 Velocity distribution alone z-axis at 90 ns

圖17 顯示了采用改進模型計算的溫度分布與文獻[16]報道結果的對比，其中：右圖為3 個時刻的溫度分布灰度圖，左圖為文獻[16]中ALEGRA MHD3D 的模擬結果，中間圖為文獻[16]給出的實驗照片。可以看出，文獻[16]的模擬結果沒有抓拍到實驗觀測到的邊緣加熱效應，而本研究的模

圖17 模擬獲得的溫度分布與文獻[16]的對比Fig.17 Comparison between simulated temperature distribution and Ref.[16]

擬結果較為清晰地顯示了前期的邊緣加熱效應和持續時間較長的角加熱效應，4 個角點附近的擴散現象也與實驗結果定性相符。

以上分析只適應于電爆炸箔的初步設計及其工作過程的早期情況。在角點出現明顯的受熱膨脹以后，建立的模型受制于固體力學模型所能描述的范圍，無法進行后續過程的模擬。事實上，隨著電流加熱形成等離子體，材料的電參數會發生復雜的變化，其動力學行為將是復雜的磁流體力學過程，需要開展更為深入的磁流體力學建模研究。

3 結論

建立了三維電磁-熱-固體力學耦合求解模型，模擬了爆炸箔在大電流加載下的早期受熱膨脹過程。基于常數電導率模型，給出了爆炸箔的磁場、電流密度、電功率密度的空間分布。基于Burgess 電導率模型，分析了早期膨脹過程中橋箔的磁場、電流密度、溫度的演化情況，指出影響其分布的決定因素為構型和橋箔電導率分布。模擬了爆炸箔的膨脹速度，計算結果較好地捕捉到了實驗中觀測到的邊緣加熱效應和角加熱效應。