RDX(210)單晶準(zhǔn)等熵壓縮實驗的數(shù)值分析*

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

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

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

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

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

1 數(shù)學(xué)模型

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

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

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

1.1 控制方程組

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

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

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

式中:k 為粘性系數(shù)。

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

進行計算，其中

1.2 數(shù)值算法

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

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

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

計算中，自由面處的壓力為零，密度為環(huán)境密度;而窗口界面處的壓力和比容需要計算得出。本文中對窗口界面處的應(yīng)力波做簡單波假設(shè)，根據(jù)

可逐步得到界面處的p(t)。式中c 為當(dāng)?shù)芈曀佟Ｐ枰⒁獾氖?上式中的v 是窗口材料的比容。樣品比容還需要進一步通過樣品材料的準(zhǔn)等熵參考線式(4)計算得到。這樣就避免了傳統(tǒng)的積分計算，不對窗口材料內(nèi)其他拉格朗日位置點進行計算，而僅計算界面信息，節(jié)省了計算時間，提高了計算效率。

2 數(shù)值結(jié)果與分析

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

圖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 加載歷史

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

2.2 力學(xué)響應(yīng)關(guān)系式

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

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

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

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

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

3 結(jié) 論

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

［1］ Vandersall K S，Reisman D B，F(xiàn)orbes 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，F(xiàn)orbes 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，F(xiàn)urnish M D，Chau R，et al.Shock Compression of Condensed Matter.New York:American Institute of Physics，2007:1165-1168.

［6］ 蔡進濤，王桂吉，趙劍衡，等.固體炸藥的磁驅(qū)動準(zhǔn)等熵壓縮實驗研究［J］.高壓物理學(xué)報，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.