ICE自由面台阶靶数据处理的直接计算方法*

金云声，孙承纬，赵剑衡，罗斌强，王桂吉，谭福利

（中国工程物理研究院流体物理研究所，四川 绵阳 621999）

对固体材料动态压缩和状态方程的研究是行星地质物理、材料合成、凝聚态物理等学科的重要内容[1]。普遍使用的动态压缩形式有冲击压缩和斜波压缩。冲击压缩通常可通过高速飞片、强脉冲激光等加载来获得，提供绝热压缩下的状态信息。而斜波压缩可通过爆轰产生膨胀气体[2]、梯度飞片[3]、磁驱动[4-7]、激光烧蚀样品材料或独立气库膜材料产生蒸汽等加载来获得[8-9]，可提供近似等熵压缩路径上的物质状态信息。斜波压缩虽然相比冲击压缩发展较晚，但已成为极端压力加载下材料动力学特性研究中重要的加载手段，并吸引研究人员付出大量努力去解读相关实验数据和提取材料压缩特性信息[10-17]。

斜波压缩实验中可测量物理量主要是界面粒子速度历史，并根据是否使用窗口可分为窗口界面速度历史和自由面速度历史。对于窗口界面速度历史，通常采用的数据处理方法有增量式阻抗匹配法、迭代求解方法、转换函数法等；而对于自由面速度历史，通常采用的数据处理方法有二分之一近似法、迭代求解方法、直接计算方法等。将这些方法按接近物理理论模型的程度来进行划分时，依次为偏离最大的增量式阻抗匹配法和二分之一近似法、处在中间的迭代求解方法以及偏离最小的直接计算方法。因转换函数法采用了数学上的摄动分析方法来实现[13-14]，对于是否具有合理物理含义的问题还需要进一步研究，所以未将其列在上述划分之中。

增量式阻抗匹配法和二分之一近似法是近似程度最大的方法，但是因具有使用简便、不增加额外不确定模型等特点，长期应用于斜波压缩实验研究的初期。

迭代求解方法是在不能直接给出材料压缩参数与界面速度历史之间关联的条件下将两界面速度历史作为匹配关系采用猜想和迭代的方式逼近隐含理论值的方法。虽然该方法存在引入无关人为因素、迭代收敛失效等风险，但因具有普适性，是目前国内外使用较多的方法[15,20]。

直接计算方法是基于材料压缩参数与界面速度历史之间物理联系所建立的近似程度较小的一类方法，是相对于迭代求解方法而取名于本文的一种直接求解方法。2010年Ockenden等[16]首先采用曲线拟合的方式建立了自由面速度历史曲线与Lagrange声速间的关系，但这对速度历史曲线的光滑性具有一定的限制。随后Hinch[17]采用数值的方式建立了该联系，弥补了该方法在曲线光滑性上的不足。然而由于上述2种方法均只在理想流体场景中进行求解并且缺少必要的检验过程，到目前为止在斜波压缩领域知之甚少，很少被列入主流数据处理方法之中。

本文中通过深入分析斜波压缩流场，独立于直接计算方法中上述2种形式关联的建立过程，在Euler坐标系下采用几何方法建立了Euler声速与自由面速度历史之间的关系，实现了直接计算过程，并通过数值实验和实验数据的检验，分析了该方法分别在有、无强度效应数据处理中的可靠性。

1 未知EOS下斜波压缩流场的直接计算过程

很多情况下，材料介质的EOS未知或已知程度不太高。此刻如何利用实验可测量数据求解整个流场并获取材料压缩特性参数的问题，在材料压缩科学研究中具有重要的实际意义。

图1 未知EOS条件下求解流场示意图Fig.1 Schematic of calculating the flow field without EOS

经实验测量可获得相同压力历史加载下2条不同厚度样品的自由面速度历史。首先对这2条自由面速度历史曲线进行数值离散，其过程中要求曲线上速度离散点对应数值和速度间隔相同，获得分别对应于2个厚度 和 的自由面速度历史h1h2和以及自由面迹线和然后以获得2条自由面迹线为基准线，用特征线的形式在一幅图中画出整个流场，示意图如图1所示，其中i和j分别为右行和左行特征线的角标，曲线(ξi,τi)为压力活塞迹线，h1和h2分别为阶梯靶实验中2个样品的厚度，曲线和分别为2个样品的自由面迹线，为2个样品自由面发出左行特征线的黎曼不变量，和分别为2个样品内坐标。具体步骤可参考文献[18]。

由特征线理论，左、右行特征线的特征不变量可用式（1）确定：

式中：为初始压力（常压），p0分别为2个样品自由面的压力，ρ为密度，c为Euler声速，p为压力。

再由实验测量自由面速度历史曲线的离散要求，有：

并由自由面常压条件，可知

将式（2）和式（3）代入式（1），可确定整个流场内左、右行特征线的黎曼不变量：

根据特征关系，流场内格点的粒子速度u(i,j)可直接求出，如式（5）所示。但由于EOS未知，声速c(i,j)还无法直接计算。现在需要采取一定的措施，使流场获得简化。

式中：f1为与未知EOS相关的函数。

式（5）可简化为：

式中：g1为可由表达式确定的函数，g2(x)为f1(g1(x))。

由式（6）可知，Euler声速c(i,j)与两角标之差（k=i-j）相关。曲线l1和格点上变量k的取值均为0，所以该曲线格点处的声速c(i,j)具有相同的值，记为c1（由自由面常压条件c1=c0）；曲线l2和格点上变量k的取值均为1，对应曲线格点声速c(i,j)取相同值c2；曲线l3和格点上变量k的取值均为2，对应曲线格点声速c(i,j)取相同值c3；以此类推，其余声速取值表达式为经过以上分析和处理，流场内未知声速变量数目获得大量减少。

在x-t空间内，在经历由自由面发出的卸载波之前，右行特征线为直线，其斜率与粒子速度、Euler声速相关。流场内格点的粒子速度已由自由面速度历史曲线完全确定，此时只要知道该简单波区右行特征线的直线斜率，即可求解对应Euler声速。

直线的斜率可由落在直线内不重合的2点来确定。从图（2）中容易看出， αi对应右行特征线直线内2点，可选其右行特征直线分别与 β1和对应2条左行特征曲线的交点，即点经简化后，Euler声速ck可写成：

为了获取以上2个交点，分别将自由面粒子迹线和设置为基准线，利用该迹线格点上已知的粒子速度和Euler声速c1分别计算曲线l2和内格点坐标，其中包含与 α2特征直线对应2个交点，即点和点再利用式（7）可将c2求出，进而又可计算曲线l3和的格点坐标。以此类推，可依次求解ck和对应lk+1曲线上的格点坐标。

最终，我们获得了 β1和对应2条左行特征曲线格点的坐标、粒子速度、Euler声速和Lagrange声速。

2 数值计算

利用特征线方法进行正向计算时，只要知道加载面速度历史以及Euler声速与原位粒子速度关系曲线，整个流场即可确定。同时，该Euler声速与原位粒子速度关系曲线又可利用直接计算方法由2个不同厚度样品自由面速度历史曲线获得。通过比对正向计算中人为设定的Euler声速与原位粒子速度关系曲线和利用直接计算方法获得相应关系曲线，即可实现对直接计算方法的检验过程。

人为设定的加载面速度历史曲线以及Euler声速与原位粒子速度关系曲线（该Euler声速与原位粒子速度关系曲线不具有真实物理含义，仅用来检验本文中方法的可靠性），分别如下式所示：

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

在以上人为设置的参数下，通过计算获得了2种厚度（2.7 mm和3.0 mm）样品的自由面速度历史，如图2所示。

利用直接计算的特征线方法由2种厚度自由面速度历史又计算了整个流场，并将2个流场叠加在一幅图中显示，如图3所示，其中绿色表示简单波区，蓝色表示复杂波区，黑色实线为粒子迹线，2条红色实线为自由面发出的第一条左行特征线。

图2 加载面速度历史和两个样品自由面速度历史Fig.2 The input data of loading surface velocity history and the free surface velocity history by calculation

图3 特征线方法计算获得流场Fig.3 The calculated flow field by characteristic method

将人为设置和计算获得的2条Euler声速与原位粒子速度关系曲线进行了比较，如图4和图5所示。比对结果显示，2条关系曲线符合的非常好。

图4 声速与原位粒子速度关系比对（利用自由面速度上升段速度历史数据）Fig.4 Comparison of sound speeds between different methods(from the rising section of the free surface velocity history

图5 声速与原位粒子速度关系比对（利用自由面速度历史下降段历史数据）Fig.5 Comparison of sound speeds between different methods(from the falling section of the free surface velocity history

本节中的检验结果说明，直接计算方法在无强度效应的数据处理中可准确逼近设定的理论值，具有较好的可靠性。

3 强度效应实验数据的处理

当实验数据含有明显的强度信息时，数据处理的难度会增加很多。最简便的处理方法有针对自由面速度历史的二分之一近似方法和针对窗口界面速度历史的增量式阻抗匹配法，但因近似程度过大而常常导致计算值与真值的偏差较大。经进一步的改进和发展，又有了迭代方法和转换函数法。但迭代方法中需要引入预先设定的强度模型，增加了不确定人为因素；而转换函数法作为一种偏数学的方法，将转换关系建立在曲线之间的外形上，尚缺少对其处理过程的物理解释，有待进一步发展。针对强度效应问题目前还没有对应理论比较完备的数据处理方法。

本文的直接计算方法作为一种新增方法，可处理含有强度效应的实验数据。下面通过两部分内容的检验，分析其可靠性。

首先，利用LSDYNA计算软件在Johnson-Cook强度模型下对不同厚度的材料模型（为保持轴向的一维性，设置的半径远大于厚度）施加相同压力历史，计算了不同厚度自由面速度历史，如图6所示。

采用直接计算方法和二分之一近似方法，由速度历史计算了声速与原位粒子速度关系曲线，并与LSDYNA计算中设定的相应理论关系曲线进行了比对，比对结果如图7所示。

从图7的比对结果中可以看出，直接计算方法获得的结果在具有强度效应场景中远优于相比二分之一近似法，并且加载路径上的数据处理结果几乎与预设的理论曲线重合。

其次，利用文献中真实的实验数据进行了处理和比较。有学者在Z磁驱动装置上[19]分别对厚度为2 594 µ m 和2 799 µ m的平面金属铜样品进行了相同加载压力历史下的阶梯靶实验，测量了自由面速度历史，如图8所示，并由此获得了斜波压缩数据。本文中利用该自由面速度历史实验数据进行处理，获得计算结果与文献中给出的结果进行了对比，比对结果如图9所示。

图6 LSDYNA计算获得速度历史（考虑了强度）Fig.6 Free surface velocity history with strength effect calculated by LSDYNA

图7 声速与原位粒子速度关系曲线的对比结果Fig.7 Comparison of the relations between Lagrange sound speed and in-situ particle velocity by different methods

图8 阶梯靶实验中金属铜的自由面速度历史Fig.8 The free surface velocity curves of Cu in step target experiment

图9 声速与自由面速度关系曲线Fig.9 Relation between Lagrange sound velocity and free surface velocity

从Lagrange声速结果比对中可以看出，直接计算方法给出的结果与文献中的结果符合较好，其计算结果好于二分之一近似方法和反积分方法。

从以上两部分检验中可以看出，含有强度效应场景中直接计算方法计算结果不仅与预设理论值具有较好的逼近，而且与文献获得结果符合较好。

4 结 论

本文针对含有自由面边界的斜波压缩流场进行研究，实现了未知EOS下斜波压缩流场的直接计算过程，并通过进行数值实验和实验数据检验发现，该方法分别在有、无强度效应数据处理中均有较好的可靠性。这不仅将其方法推向斜波压缩领域成为主流数据处理方法，而且为探索具有较完整理论的强度效应数据处理方法提供了可靠的新途径。

感谢中国工程物理研究院流体物理研究所张红平博士在利用反积分方法进行数据处理过程中提供的帮助。