RDX(210)單晶準等熵壓縮實驗的數值分析*

張紅平，李 牧，趙劍衡，孫承緯，祝文軍

(中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理國防科技重點實驗室，四川 綿陽621900)

研究含能材料的物態性能及熱力學參數，在國防工業中具有重要意義。但含能材料不同于普通材料，采用傳統的沖擊實驗物態測量方法，在高溫高壓狀態下樣品材料容易發生反應，給測量帶來困難。近年來迅速發展的準等熵壓縮實驗測量方法則可以解決該困難，它在獲取高壓狀態的同時，仍可保持樣品材料內較低的溫升，含能材料不易發生反應。目前已有一些實驗結果［1-5］和相關實驗工作［6］，但對該類實驗數據的處理和分析還不完善。

依靠樣品內嵌量計和多臺階樣品實驗提供數據的拉格朗日方法，只能在有限個離散點上進行積分，無法得到樣品內部任意位置的信息。而傳統的積分方法雖然可以得到樣品內部信息，但仍存在不足:電磁加載實驗中一般用測到的電流波形計算得到樣品加載面上的信息，進而向樣品內部進行積分計算，這種方法通常需要考慮磁擴散效應和焦耳加熱導致金屬極板的電離情況，而目前這方面的數學模型建立還不完善，計算也比較復雜，誤差較大。而用反積分方法［7］向樣品內部反向計算，不僅可以避免磁擴散和電極板電離等情況的計算，還可以得到樣品內部任意位置處的信息。

RDX 作為多種含能復合材料的基本成分，對其物態參數和力學性能進行研究具有重要意義。目前的準等熵壓縮實驗中固體單晶RDX(210)材料［5］的溫升較小，沒有發生固液相變，因此其強度效應不可忽略，不能用傳統的流體模型近似處理。鑒于此，本文中將在RDX(210)實驗數據的反積分解讀和分析中加入流體線彈塑性模型，對RDX(210)的強度效應進行描述。基于所測的樣品后界面的速度歷史，向樣品內部反演計算加載面信息。通過臺階靶樣品實驗，以共同的加載歷史為判據，調整樣品力學響應關系式的參數，最終獲取樣品材料的物態參數，并分析其力學性能。

1 數學模型

含能材料的準等熵壓縮屬于一維平面壓縮實驗，可用單軸應變條件下的力學性能進行數學描述。樣品材料的壓縮狀態可用Lagrange 形式的質量方程?v/?t = ?u/(ρ0?x)和動量方程?u/?t =－?σ/(ρ0?x)描述。其中σ、u、v 分別為軸向應力、速度和比容，ρ0為初始密度。在單軸應變下，總應力為σ=p+4τ/3+q，其中p 為靜水壓力，τ 為剪應力，q 為粘性力。

上述質量和動量方程的封閉需要補充一個力學響應關系式才可進行計算。在準等熵壓縮實驗中，可選取壓縮到εS時的準等熵參考線［8］

式中:c0為初始聲速，p0為初始壓力，Γ0為Grüneisen 系數，ε=1－v/v0為體應變。

1.1 控制方程組

在電磁加載實驗中，為了避免磁擴散效應和焦耳加熱效應的數學模型建立和計算，Z 機器上的實驗結果也曾實現過測量得到樣品前界面的速度剖面［4-5］，以此為邊界條件，向樣品內部積分計算，直至得到樣品后界面的速度、壓力和比容，進而通過和后界面實驗測量速度剖面的比較，驗證其積分方法的有效性。但這種做法需要相同加載條件下至少4 個不同厚度臺階樣品的動態壓縮實驗，目前國內的實驗條件還無法做到。因此為了建立適用于更多準等熵壓縮實驗的數據處理技術，這里利用實驗測到的樣品后界面的速度歷史作為“初始條件”，用反積分方法向樣品內部反向計算，對應的反積分控制方程為

式中:p=σ－4τ/3－q，樣品材料在屈服之前τ ≤Y/2(Y 為屈服強度)，滿足線彈性本構關系［7］

式中:ν 為泊松比。屈服之后τ=Y/2，粘性力取PIC 粘性關系［9］

式中:k 為粘性系數。

在v=F(p)中，可利用

進行計算，其中

1.2 數值算法

為保證解的二階精度，上述反積分控制方程組對應的差分格式為

式中:p(x－dx，t)=σ(x－dx，t)－4τ(x－dx，t)/3－q(x－dx，t)，dx 和dt 分別表示空間步長和時間步長。

樣品后自由面/窗口界面的速度歷史作為輸入情況，由式(6)求出材料內部某位置節點處的應力歷史，由式(7)求得比容歷史，由式(8)求得該位置的速度歷史;逐次向內可得到加載面的應力、比容和速度歷史。其中比容歷史的計算采用二階Runge-Kutta 法(PEC 預估校正格式)對式(4)進行計算。

計算中，自由面處的壓力為零，密度為環境密度;而窗口界面處的壓力和比容需要計算得出。本文中對窗口界面處的應力波做簡單波假設，根據

可逐步得到界面處的p(t)。式中c 為當地聲速。需要注意的是:上式中的v 是窗口材料的比容。樣品比容還需要進一步通過樣品材料的準等熵參考線式(4)計算得到。這樣就避免了傳統的積分計算，不對窗口材料內其他拉格朗日位置點進行計算，而僅計算界面信息，節省了計算時間，提高了計算效率。

2 數值結果與分析

為了驗證考慮流體線彈塑性模型后的反積分方法的有效性，選取Sandia 實驗室Z 機器上的1489 號實驗的部分RDX 數據［5］進行處理和分析。Z1489 號實驗中對單晶RDX(210)的準等熵壓縮測量數據表明，RDX(210)未發生反應，且表現出明顯的彈塑性能。圖1 給出了M.R.Baer 等［4-5］測量的2 種厚度h 分別為0.592 和0.724 mm 的帶LiF 窗口的RDX(210)樣品的界面速度剖面。在350 ns 內后界面粒子速度上升到了300 m/s，從圖中可以清晰地看到應力波在上升過程中表現為雙波結構，上升沿前70 ns表現為彈性波，經過一個應力平臺轉變為塑性波。為了向樣品內部進行反演計算，邊界信息僅有速度歷史是不夠的。圖2 中給出了用簡單波假設計算的樣品后界面上的壓力歷史，和進一步通過準等熵參考線計算得到的比容歷史，壓力和比容歷史都有明顯的雙波結構。

圖1 不同厚度RDX(210)的后界面速度歷史［5］Fig.1 Particle velocity profiles at the interfaces between LiF windows and RDX(210)crystals with different thicknesses［5］

圖2 不同厚度RDX(210)的后界面比容歷史和壓力歷史Fig.2 Specific volume and pressure profiles measured at the interfaces between LiF windows and RDX(210)crystals with different thicknesses

2.1 加載歷史

經過反積分計算以后，加載面的速度歷史和壓力歷史如圖3 所示。圖3 中還給出了本文計算結果和文獻［4-5］中的速度加載歷史的比較，可以看出在上升階段，速度結果吻合較好。不同厚度樣品的壓力加載歷史在500 ns 之前也僅有細微差別。本文數值模擬結果在加載面上的有效計算時間是根據不同厚度樣品材料的共同加載歷史確定的。即在較薄樣品的加載應力波到達窗口界面，經窗口反射后回到加載面之前為有效計算時間。在該時間點之后，薄厚樣品加載面上的信息發生變化，不再具有相同的加載歷史。基于此，可以估算出樣品材料在加載面上的有效時間段終點約為500 ns。

2.2 力學響應關系式

反積分方法以不同厚度樣品的相同加載歷史為判據，不斷調整力學響應關系式中的參數，以此獲取數學上更準確的力學響應關系式。在本文所選的等熵參考線方程中，共有4 個參數Γ0、c0、λ 和Y。采用下山單純形優化方法進行參數搜索。初始值Γ0=1.29，c0=2.17 km/s，λ=2.78，Y=439 MPa，經過參數搜索后得到Γ0=1.615，c0=2.72 km/s，λ=3.48，Y=600 MPa。按照參數在沖擊加載下的物理表述，在該RDX(210)的準等熵壓縮實驗中:Grüneisen 系數增大，即定容條件下壓力對單位體積內能的變化率增大;零壓下聲速變大，和材料壓縮程度相關的參數λ 變大;屈服強度增大，即準等熵壓縮下樣品的彈性屈服極限較Hugoniot 彈性極限高。但這些表述用于準等熵壓縮實驗的物理描述還有待推敲，這是由于力學響應關系式及準等熵參考線中的參數在準等熵壓縮實驗中，其物理意義已經發生變化，原先的物理解釋只適用于沖擊加載情況，這里還需要根據力學響應關系式重新定義。另外這里的參數搜索過程是一個數學過程，得到的參數是一個純數學解，真實的物理意義還有待進一步考證。

圖3 不同厚度RDX(210)的加載面速度和壓力Fig.3 Particle velocity and pressure at loading face in RDX(210)crystal，Baer’s driving velocity［5］is here for comparison

2.3 拉格朗日聲速

圖4 (a)給出了樣品內部不同位置處的粒子速度剖面，從圖中可以看出，由于后界面LiF 窗口的阻抗遠高于RDX 阻抗，所以在窗口界面會反射一個壓縮波，該反射波與入射應力波疊加后導致樣品材料在緊靠窗口位置處的壓力和密度升高，而粒子速度取決于該應力下窗口材料的速度。圖4(b)給出的是RDX 內部距離樣品加載面不同長度l 的位置處拉格朗日縱波聲速cL隨粒子速度的變化情況，其中=?σ/?ρ，可以看出cL首先經過彈性區，由于塑性變形而下降，隨后完全進入塑性區。不同位置處樣品發生塑性變形后，cL隨粒子速度的變化也不相同，離樣品加載面越遠，越靠近高阻抗窗口，cL的斜率越大，即縱波聲速cL相對粒子速度更快，這正是壓縮應力波在帶高阻抗窗口的樣品內部的傳播情況。

圖4 RDX(210)內部不同位置處的速度歷史和拉格朗日聲速Fig.4 Particle velocity profiles and Lagrangian sound velocities at different positions in RDX(210)crystal

圖5 RDX(210)的準等熵p-v 參考線和應力應變曲線Fig.5 Stress-strain curves，quasi-isentrope and Lagrangian sound velocity for quasiisentropically loaded RDX(210)crystal

與粒子速度和應力等變量相比，在一維應變狀態下，比容作為基本變量可以更直觀地表述材料內部的力學參量變化規律。圖5 中以比容為橫坐標，給出了樣品中壓力p、應力σ 和拉格朗日縱波聲速cL的變化規律。可以看出，在屈服極限點之前，材料處于線彈性階段，σ 和cL以相同的斜率隨比容的減小而線性增長。對于線彈塑性模型，可以得到應力屈服極限σIEL=Y(1－ν)/(1－2ν)，圖5 中給出了該屈服點v=0.542 cm3/g，σIEL=0.77 GPa。另外圖中還給出了該RDX(210)準等熵壓縮實驗的準等熵參考p-v 線，且p-v 線位于σ-v 線的下方，2 條參考線的差即為偏應力的貢獻。

3 結 論

針對RDX(210)的準等熵壓縮實驗數據，建立了考慮線彈塑性效應的反積分數學模型，通過比較數值模擬結果和文獻結果，說明本文的數學模型可以恰當地描述該準等熵壓縮實驗。基于反積分方法可以獲取樣品內部信息的優勢，給出了計算得到的RDX(210)內部不同位置的拉格朗日聲速和粒子速度，這對分析應力波在樣品內部的傳播過程有積極意義。基于數學搜索技術，用反積分方法得到了RDX(210)的準等熵參考線。但這種材料力學響應關系式的數學獲取過程還需要進一步分析其物理含義。

［1］ Vandersall K S，Reisman D B，Forbes J W，et al.Isentropic compression experiments performed by LLNL on energetic material samples using the Z accelerator［R］.UCRL-TR-236063，LLNL，2007.

［2］ Hare D E，Forbes J W，Reisman D B，et al.Isentropic compression loading of octahydro-1，3，5，7-tetranitro-1，3，5，7-tetrazocine(HMX)and the pressure-induced phase transition at 27 GPa［J］.Applied Physics Letters，2004，85(6):949-951.

［3］ Hooks D E，Hayes D B，Hare D E，et al.Isentropic compression of cyclotetramethylene tetranitramine(HMX)single crystals to 50 GPa［J］.Journal of Applied Physics，2006，99(12):124901.

［4］ Baer M R，Hall C A，Gustavsen R L，et al.Isentropic loading experiments of a plastic bonded explosive and constituents［J］.Journal of Applied Physics，2007，101(3):034906.

［5］ Baer M R，Hobbs M L，Hall C A，et al.Isentropic compression studies of energetic composite constituents［C］∥Elert M，Furnish M D，Chau R，et al.Shock Compression of Condensed Matter.New York:American Institute of Physics，2007:1165-1168.

［6］ 蔡進濤，王桂吉，趙劍衡，等.固體炸藥的磁驅動準等熵壓縮實驗研究［J］.高壓物理學報，2010，24(6):401-408.CAI Jin-tao，WANG Gui-ji，ZHAO Jian-heng，et al.Magnetically driven quasi-isentropic compression experiment of solid explosion［J］.Chinese Journal of High Pressure Physics，2010，24(6):401-408.

［7］ Hayes D.Backward integration of the equations of motion to correct for free surface perturbations［R］.SAND2001-1440.Sandia National Laboratories，2001.

［8］ 莫建軍.金屬等熵壓縮線計算及爆轟加載等熵壓縮實驗研究［D］.綿陽:中國工程物理研究院，2006:46-49.

［9］ 孫承緯.一維沖擊波和爆轟波計算程序SSS［J］.計算物理，1986，3(2):142-154.SUN Cheng-wei.SSS:A code for computing one dimensional shock and detonation wave propagation［J］.Chinese Journal of Computational Physics，1986，3(2):142-154.