ICE自由面臺(tái)階靶數(shù)據(jù)處理的直接計(jì)算方法*

金云聲，孫承緯，趙劍衡，羅斌強(qiáng)，王桂吉，譚福利

（中國(guó)工程物理研究院流體物理研究所，四川 綿陽(yáng) 621999）

對(duì)固體材料動(dòng)態(tài)壓縮和狀態(tài)方程的研究是行星地質(zhì)物理、材料合成、凝聚態(tài)物理等學(xué)科的重要內(nèi)容[1]。普遍使用的動(dòng)態(tài)壓縮形式有沖擊壓縮和斜波壓縮。沖擊壓縮通?？赏ㄟ^(guò)高速飛片、強(qiáng)脈沖激光等加載來(lái)獲得，提供絕熱壓縮下的狀態(tài)信息。而斜波壓縮可通過(guò)爆轟產(chǎn)生膨脹氣體[2]、梯度飛片[3]、磁驅(qū)動(dòng)[4-7]、激光燒蝕樣品材料或獨(dú)立氣庫(kù)膜材料產(chǎn)生蒸汽等加載來(lái)獲得[8-9]，可提供近似等熵壓縮路徑上的物質(zhì)狀態(tài)信息。斜波壓縮雖然相比沖擊壓縮發(fā)展較晚，但已成為極端壓力加載下材料動(dòng)力學(xué)特性研究中重要的加載手段，并吸引研究人員付出大量努力去解讀相關(guān)實(shí)驗(yàn)數(shù)據(jù)和提取材料壓縮特性信息[10-17]。

斜波壓縮實(shí)驗(yàn)中可測(cè)量物理量主要是界面粒子速度歷史，并根據(jù)是否使用窗口可分為窗口界面速度歷史和自由面速度歷史。對(duì)于窗口界面速度歷史，通常采用的數(shù)據(jù)處理方法有增量式阻抗匹配法、迭代求解方法、轉(zhuǎn)換函數(shù)法等；而對(duì)于自由面速度歷史，通常采用的數(shù)據(jù)處理方法有二分之一近似法、迭代求解方法、直接計(jì)算方法等。將這些方法按接近物理理論模型的程度來(lái)進(jìn)行劃分時(shí)，依次為偏離最大的增量式阻抗匹配法和二分之一近似法、處在中間的迭代求解方法以及偏離最小的直接計(jì)算方法。因轉(zhuǎn)換函數(shù)法采用了數(shù)學(xué)上的攝動(dòng)分析方法來(lái)實(shí)現(xiàn)[13-14]，對(duì)于是否具有合理物理含義的問(wèn)題還需要進(jìn)一步研究，所以未將其列在上述劃分之中。

增量式阻抗匹配法和二分之一近似法是近似程度最大的方法，但是因具有使用簡(jiǎn)便、不增加額外不確定模型等特點(diǎn)，長(zhǎng)期應(yīng)用于斜波壓縮實(shí)驗(yàn)研究的初期。

迭代求解方法是在不能直接給出材料壓縮參數(shù)與界面速度歷史之間關(guān)聯(lián)的條件下將兩界面速度歷史作為匹配關(guān)系采用猜想和迭代的方式逼近隱含理論值的方法。雖然該方法存在引入無(wú)關(guān)人為因素、迭代收斂失效等風(fēng)險(xiǎn)，但因具有普適性，是目前國(guó)內(nèi)外使用較多的方法[15,20]。

直接計(jì)算方法是基于材料壓縮參數(shù)與界面速度歷史之間物理聯(lián)系所建立的近似程度較小的一類方法，是相對(duì)于迭代求解方法而取名于本文的一種直接求解方法。2010年Ockenden等[16]首先采用曲線擬合的方式建立了自由面速度歷史曲線與Lagrange聲速間的關(guān)系，但這對(duì)速度歷史曲線的光滑性具有一定的限制。隨后Hinch[17]采用數(shù)值的方式建立了該聯(lián)系，彌補(bǔ)了該方法在曲線光滑性上的不足。然而由于上述2種方法均只在理想流體場(chǎng)景中進(jìn)行求解并且缺少必要的檢驗(yàn)過(guò)程，到目前為止在斜波壓縮領(lǐng)域知之甚少，很少被列入主流數(shù)據(jù)處理方法之中。

本文中通過(guò)深入分析斜波壓縮流場(chǎng)，獨(dú)立于直接計(jì)算方法中上述2種形式關(guān)聯(lián)的建立過(guò)程，在Euler坐標(biāo)系下采用幾何方法建立了Euler聲速與自由面速度歷史之間的關(guān)系，實(shí)現(xiàn)了直接計(jì)算過(guò)程，并通過(guò)數(shù)值實(shí)驗(yàn)和實(shí)驗(yàn)數(shù)據(jù)的檢驗(yàn)，分析了該方法分別在有、無(wú)強(qiáng)度效應(yīng)數(shù)據(jù)處理中的可靠性。

1 未知EOS下斜波壓縮流場(chǎng)的直接計(jì)算過(guò)程

很多情況下，材料介質(zhì)的EOS未知或已知程度不太高。此刻如何利用實(shí)驗(yàn)可測(cè)量數(shù)據(jù)求解整個(gè)流場(chǎng)并獲取材料壓縮特性參數(shù)的問(wèn)題，在材料壓縮科學(xué)研究中具有重要的實(shí)際意義。

圖1 未知EOS條件下求解流場(chǎng)示意圖Fig.1 Schematic of calculating the flow field without EOS

經(jīng)實(shí)驗(yàn)測(cè)量可獲得相同壓力歷史加載下2條不同厚度樣品的自由面速度歷史。首先對(duì)這2條自由面速度歷史曲線進(jìn)行數(shù)值離散，其過(guò)程中要求曲線上速度離散點(diǎn)對(duì)應(yīng)數(shù)值和速度間隔相同，獲得分別對(duì)應(yīng)于2個(gè)厚度 和 的自由面速度歷史h1h2和以及自由面跡線和然后以獲得2條自由面跡線為基準(zhǔn)線，用特征線的形式在一幅圖中畫(huà)出整個(gè)流場(chǎng)，示意圖如圖1所示，其中i和j分別為右行和左行特征線的角標(biāo)，曲線(ξi,τi)為壓力活塞跡線，h1和h2分別為階梯靶實(shí)驗(yàn)中2個(gè)樣品的厚度，曲線和分別為2個(gè)樣品的自由面跡線，為2個(gè)樣品自由面發(fā)出左行特征線的黎曼不變量，和分別為2個(gè)樣品內(nèi)坐標(biāo)。具體步驟可參考文獻(xiàn)[18]。

由特征線理論，左、右行特征線的特征不變量可用式（1）確定：

式中：為初始?jí)毫Γǔ海?，p0分別為2個(gè)樣品自由面的壓力，ρ為密度，c為Euler聲速，p為壓力。

再由實(shí)驗(yàn)測(cè)量自由面速度歷史曲線的離散要求，有：

并由自由面常壓條件，可知

將式（2）和式（3）代入式（1），可確定整個(gè)流場(chǎng)內(nèi)左、右行特征線的黎曼不變量：

根據(jù)特征關(guān)系，流場(chǎng)內(nèi)格點(diǎn)的粒子速度u(i,j)可直接求出，如式（5）所示。但由于EOS未知，聲速c(i,j)還無(wú)法直接計(jì)算?，F(xiàn)在需要采取一定的措施，使流場(chǎng)獲得簡(jiǎn)化。

式中：f1為與未知EOS相關(guān)的函數(shù)。

式（5）可簡(jiǎn)化為：

式中：g1為可由表達(dá)式確定的函數(shù)，g2(x)為f1(g1(x))。

由式（6）可知，Euler聲速c(i,j)與兩角標(biāo)之差（k=i-j）相關(guān)。曲線l1和格點(diǎn)上變量k的取值均為0，所以該曲線格點(diǎn)處的聲速c(i,j)具有相同的值，記為c1（由自由面常壓條件c1=c0）；曲線l2和格點(diǎn)上變量k的取值均為1，對(duì)應(yīng)曲線格點(diǎn)聲速c(i,j)取相同值c2；曲線l3和格點(diǎn)上變量k的取值均為2，對(duì)應(yīng)曲線格點(diǎn)聲速c(i,j)取相同值c3；以此類推，其余聲速取值表達(dá)式為經(jīng)過(guò)以上分析和處理，流場(chǎng)內(nèi)未知聲速變量數(shù)目獲得大量減少。

在x-t空間內(nèi)，在經(jīng)歷由自由面發(fā)出的卸載波之前，右行特征線為直線，其斜率與粒子速度、Euler聲速相關(guān)。流場(chǎng)內(nèi)格點(diǎn)的粒子速度已由自由面速度歷史曲線完全確定，此時(shí)只要知道該簡(jiǎn)單波區(qū)右行特征線的直線斜率，即可求解對(duì)應(yīng)Euler聲速。

直線的斜率可由落在直線內(nèi)不重合的2點(diǎn)來(lái)確定。從圖（2）中容易看出， αi對(duì)應(yīng)右行特征線直線內(nèi)2點(diǎn)，可選其右行特征直線分別與 β1和對(duì)應(yīng)2條左行特征曲線的交點(diǎn)，即點(diǎn)經(jīng)簡(jiǎn)化后，Euler聲速ck可寫(xiě)成：

為了獲取以上2個(gè)交點(diǎn)，分別將自由面粒子跡線和設(shè)置為基準(zhǔn)線，利用該跡線格點(diǎn)上已知的粒子速度和Euler聲速c1分別計(jì)算曲線l2和內(nèi)格點(diǎn)坐標(biāo)，其中包含與 α2特征直線對(duì)應(yīng)2個(gè)交點(diǎn)，即點(diǎn)和點(diǎn)再利用式（7）可將c2求出，進(jìn)而又可計(jì)算曲線l3和的格點(diǎn)坐標(biāo)。以此類推，可依次求解ck和對(duì)應(yīng)lk+1曲線上的格點(diǎn)坐標(biāo)。

最終，我們獲得了 β1和對(duì)應(yīng)2條左行特征曲線格點(diǎn)的坐標(biāo)、粒子速度、Euler聲速和Lagrange聲速。

2 數(shù)值計(jì)算

利用特征線方法進(jìn)行正向計(jì)算時(shí)，只要知道加載面速度歷史以及Euler聲速與原位粒子速度關(guān)系曲線，整個(gè)流場(chǎng)即可確定。同時(shí)，該Euler聲速與原位粒子速度關(guān)系曲線又可利用直接計(jì)算方法由2個(gè)不同厚度樣品自由面速度歷史曲線獲得。通過(guò)比對(duì)正向計(jì)算中人為設(shè)定的Euler聲速與原位粒子速度關(guān)系曲線和利用直接計(jì)算方法獲得相應(yīng)關(guān)系曲線，即可實(shí)現(xiàn)對(duì)直接計(jì)算方法的檢驗(yàn)過(guò)程。

人為設(shè)定的加載面速度歷史曲線以及Euler聲速與原位粒子速度關(guān)系曲線（該Euler聲速與原位粒子速度關(guān)系曲線不具有真實(shí)物理含義，僅用來(lái)檢驗(yàn)本文中方法的可靠性），分別如下式所示：

式中：初始聲速C0不妨取為5 328.0 m/s。

在以上人為設(shè)置的參數(shù)下，通過(guò)計(jì)算獲得了2種厚度（2.7 mm和3.0 mm）樣品的自由面速度歷史，如圖2所示。

利用直接計(jì)算的特征線方法由2種厚度自由面速度歷史又計(jì)算了整個(gè)流場(chǎng)，并將2個(gè)流場(chǎng)疊加在一幅圖中顯示，如圖3所示，其中綠色表示簡(jiǎn)單波區(qū)，藍(lán)色表示復(fù)雜波區(qū)，黑色實(shí)線為粒子跡線，2條紅色實(shí)線為自由面發(fā)出的第一條左行特征線。

圖2 加載面速度歷史和兩個(gè)樣品自由面速度歷史Fig.2 The input data of loading surface velocity history and the free surface velocity history by calculation

圖3 特征線方法計(jì)算獲得流場(chǎng)Fig.3 The calculated flow field by characteristic method

將人為設(shè)置和計(jì)算獲得的2條Euler聲速與原位粒子速度關(guān)系曲線進(jìn)行了比較，如圖4和圖5所示。比對(duì)結(jié)果顯示，2條關(guān)系曲線符合的非常好。

圖4 聲速與原位粒子速度關(guān)系比對(duì)（利用自由面速度上升段速度歷史數(shù)據(jù)）Fig.4 Comparison of sound speeds between different methods(from the rising section of the free surface velocity history

圖5 聲速與原位粒子速度關(guān)系比對(duì)（利用自由面速度歷史下降段歷史數(shù)據(jù)）Fig.5 Comparison of sound speeds between different methods(from the falling section of the free surface velocity history

本節(jié)中的檢驗(yàn)結(jié)果說(shuō)明，直接計(jì)算方法在無(wú)強(qiáng)度效應(yīng)的數(shù)據(jù)處理中可準(zhǔn)確逼近設(shè)定的理論值，具有較好的可靠性。

3 強(qiáng)度效應(yīng)實(shí)驗(yàn)數(shù)據(jù)的處理

當(dāng)實(shí)驗(yàn)數(shù)據(jù)含有明顯的強(qiáng)度信息時(shí)，數(shù)據(jù)處理的難度會(huì)增加很多。最簡(jiǎn)便的處理方法有針對(duì)自由面速度歷史的二分之一近似方法和針對(duì)窗口界面速度歷史的增量式阻抗匹配法，但因近似程度過(guò)大而常常導(dǎo)致計(jì)算值與真值的偏差較大。經(jīng)進(jìn)一步的改進(jìn)和發(fā)展，又有了迭代方法和轉(zhuǎn)換函數(shù)法。但迭代方法中需要引入預(yù)先設(shè)定的強(qiáng)度模型，增加了不確定人為因素；而轉(zhuǎn)換函數(shù)法作為一種偏數(shù)學(xué)的方法，將轉(zhuǎn)換關(guān)系建立在曲線之間的外形上，尚缺少對(duì)其處理過(guò)程的物理解釋，有待進(jìn)一步發(fā)展。針對(duì)強(qiáng)度效應(yīng)問(wèn)題目前還沒(méi)有對(duì)應(yīng)理論比較完備的數(shù)據(jù)處理方法。

本文的直接計(jì)算方法作為一種新增方法，可處理含有強(qiáng)度效應(yīng)的實(shí)驗(yàn)數(shù)據(jù)。下面通過(guò)兩部分內(nèi)容的檢驗(yàn)，分析其可靠性。

首先，利用LSDYNA計(jì)算軟件在Johnson-Cook強(qiáng)度模型下對(duì)不同厚度的材料模型（為保持軸向的一維性，設(shè)置的半徑遠(yuǎn)大于厚度）施加相同壓力歷史，計(jì)算了不同厚度自由面速度歷史，如圖6所示。

采用直接計(jì)算方法和二分之一近似方法，由速度歷史計(jì)算了聲速與原位粒子速度關(guān)系曲線，并與LSDYNA計(jì)算中設(shè)定的相應(yīng)理論關(guān)系曲線進(jìn)行了比對(duì)，比對(duì)結(jié)果如圖7所示。

從圖7的比對(duì)結(jié)果中可以看出，直接計(jì)算方法獲得的結(jié)果在具有強(qiáng)度效應(yīng)場(chǎng)景中遠(yuǎn)優(yōu)于相比二分之一近似法，并且加載路徑上的數(shù)據(jù)處理結(jié)果幾乎與預(yù)設(shè)的理論曲線重合。

其次，利用文獻(xiàn)中真實(shí)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了處理和比較。有學(xué)者在Z磁驅(qū)動(dòng)裝置上[19]分別對(duì)厚度為2 594 μ m 和2 799 μ m的平面金屬銅樣品進(jìn)行了相同加載壓力歷史下的階梯靶實(shí)驗(yàn)，測(cè)量了自由面速度歷史，如圖8所示，并由此獲得了斜波壓縮數(shù)據(jù)。本文中利用該自由面速度歷史實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理，獲得計(jì)算結(jié)果與文獻(xiàn)中給出的結(jié)果進(jìn)行了對(duì)比，比對(duì)結(jié)果如圖9所示。

圖6 LSDYNA計(jì)算獲得速度歷史（考慮了強(qiáng)度）Fig.6 Free surface velocity history with strength effect calculated by LSDYNA

圖7 聲速與原位粒子速度關(guān)系曲線的對(duì)比結(jié)果Fig.7 Comparison of the relations between Lagrange sound speed and in-situ particle velocity by different methods

圖8 階梯靶實(shí)驗(yàn)中金屬銅的自由面速度歷史Fig.8 The free surface velocity curves of Cu in step target experiment

圖9 聲速與自由面速度關(guān)系曲線Fig.9 Relation between Lagrange sound velocity and free surface velocity

從Lagrange聲速結(jié)果比對(duì)中可以看出，直接計(jì)算方法給出的結(jié)果與文獻(xiàn)中的結(jié)果符合較好，其計(jì)算結(jié)果好于二分之一近似方法和反積分方法。

從以上兩部分檢驗(yàn)中可以看出，含有強(qiáng)度效應(yīng)場(chǎng)景中直接計(jì)算方法計(jì)算結(jié)果不僅與預(yù)設(shè)理論值具有較好的逼近，而且與文獻(xiàn)獲得結(jié)果符合較好。

4 結(jié) 論

本文針對(duì)含有自由面邊界的斜波壓縮流場(chǎng)進(jìn)行研究，實(shí)現(xiàn)了未知EOS下斜波壓縮流場(chǎng)的直接計(jì)算過(guò)程，并通過(guò)進(jìn)行數(shù)值實(shí)驗(yàn)和實(shí)驗(yàn)數(shù)據(jù)檢驗(yàn)發(fā)現(xiàn)，該方法分別在有、無(wú)強(qiáng)度效應(yīng)數(shù)據(jù)處理中均有較好的可靠性。這不僅將其方法推向斜波壓縮領(lǐng)域成為主流數(shù)據(jù)處理方法，而且為探索具有較完整理論的強(qiáng)度效應(yīng)數(shù)據(jù)處理方法提供了可靠的新途徑。

感謝中國(guó)工程物理研究院流體物理研究所張紅平博士在利用反積分方法進(jìn)行數(shù)據(jù)處理過(guò)程中提供的幫助。