-
准等熵压缩实验对于材料物态方程、高能量密度物理以及武器物理的研究等具有重要意义。目前, 国外公开发表的金属材料准等熵压缩实验数据中, 大多数工作得到的样品应力-应变关系曲线都包含了大量信息[1-2]:准等熵过程产生的熵增导致实测温度高于等熵压缩温度, 因此, 相同比容下的流体静水压多出一个由熵增引起的增量; 此外, 准等熵过程温升较低, 一般不会出现冲击软化现象, 且连续加载使应力-应变状态更易处于材料的上屈服面, 材料的动态屈服强度和偏应力较冲击加载情况更高, 偏应力占总应力的比重也更高, 故相同比容下, 应力-应变曲线中有偏应力的贡献。为此, 要得到基于静水压的热力学状态方程参考线, 就必须考虑流动应力和熵增问题, 流动应力的问题我们已初步讨论过[3], 本研究重点关注熵增问题。
在开展准等熵压缩实验时, 人们经常思考的一个问题就是准等熵和等熵的区别, 准等熵过程中的熵增究竟有多少?熵增一般与高应变率和大变形相关, 可以分为力学耗散和热耗散两部分。力学耗散主要由各种非弹性变形机理引起, 如相变和变形中的率效应等。热耗散则主要由不可逆的热传导过程引起。力学耗散的大小直接决定大幅度应力波的演化及一些典型特征, 还可用于区分力学做功和热传导的贡献。力、热贡献的不同反过来又可用于确定温度剖面。所以, 热力学耗散的研究对理解高压下的材料状态至关重要。近年来, Ding等人[4]已经开始在准等熵压缩实验数据分析中, 考虑应变率因素及黏性耗散问题, 以分析材料内部不同拉格朗日位置处的原位粒子速度及材料声速。
本研究基于激光加载和磁驱动加载准等熵压缩金属铝实验, 结合反积分-正积分迭代计算方法, 讨论由不同加载应变率带来的材料热力学耗散及熵增问题, 并进一步分析应变率相关的黏性耗散及热传导对材料应力-应变曲线和拉格朗日声速的影响。
-
以往对准等熵压缩样品进行分析计算时, 大多采用双曲型方程描述应力波的传播, 数值计算较易进行。但考虑样品材料的剪切强度、温升、应变率率相关性等热动力学响应关系后, 控制方程已不再是单纯的双曲型方程; 尤其在考虑能量方程及热传导后, 控制方程变为抛物型方程, 由样品后表面(样品/窗口界面或自由面)的数据出发进行反向积分计算较为困难。为此, 需要开展反积分-正积分迭代计算。
首先进行反积分初始条件的预处理:由反积分得到材料加载面信息后, 通过流体动力学方程计算(这里称其为正积分), 将反积分不能考虑的热传导等信息考虑在内, 得到样品后表面的速度、应力和温度等信息。然后, 以样品后表面处的速度、应力和温度等作为初始条件, 向样品内部反向计算, 并针对材料的不同动力学行为, 确定其力学响应关系(计算中需要对材料在准等熵压缩实验中的力学响应关系进行理论分析, 根据冲击动力学相关原理及冲击实验结果, 给出初步关系式, 进而参与反向及正向计算以确定其参数值); 以台阶靶共同的加载历史(后表面反射应力波到达加载面之前)以及后表面处计算结果和测量结果的比较作为判据, 迭代计算台阶靶样品的力学信息, 确定部分力学响应关系式的参数初值。
材料在等熵加载下, 可用宏观流体动力学的连续介质守恒方程进行描述。按照应力波研究惯例, 使用拉格朗日坐标下的质量守恒方程、动量守恒方程和材料本构关系(物性方程)作为控制方程
$\partial v / \partial t=v_{0}(\partial u / \partial x)$ $p=F(v)$ $\sigma=p+\sigma^{\prime}+Q$ $\partial u / \partial t=-v_{0}(\partial \sigma / \partial x)$ $\partial e / \partial t=-\sigma(\partial v / \partial t)+\kappa\left(\partial^{2} T / \partial x^{2}\right)$ 式中:u、v、p、e、σ、σ′和Q分别表示粒子速度、比容、压力、比内能、应力、偏应力和黏性力, 下标“0”表示初始状态; x、t分别为空间和时间变量; p=F(v)根据所采用状态方程的不同, 有相应的形式, 一般选Hugoniot线、等熵线或Grüneisen状态方程; 偏应力σ′与屈服强度Y及剪切模量G关联密切, 材料在屈服之前, 即σ′≤2Y/3时, 满足本构关系
$\frac{\partial \sigma^{\prime}}{\partial t}=2 G\left(-\frac{\partial u}{\partial x}+\frac{1}{3 v} \frac{\partial v}{\partial t}\right)$ 屈服之后满足σ′=2Y/3。
在准等熵压缩实验中, 由于材料的加载历史较难确定, 一般用激光干涉测速仪(VISAR)对材料自由面或样品/窗口界面进行测量, 得到自由面或界面的粒子速度历史。由此, (1)式~(5)式的积分计算步骤需要重新安排为反向积分步骤, 控制方程如下
$\partial_{\sigma} / \partial x=-\left(1 / v_{0}\right)(\partial u / \partial t)$ $p=\sigma-\sigma^{\prime}-Q$ $v=F^{-1}(p)$ $\partial u / \partial x=\partial\left(v / v_{0}\right) / \partial t=\partial \varepsilon / \partial t$ $\partial e / \partial t=-\sigma(\partial v / \partial t)$ 式中:ε为应变。需要说明的是, 在反积分计算中无法考虑热传导问题, 原因之一是考虑热传导后的控制方程经时空转换后不再是双曲型方程, 解的唯一性无法得到保证。
-
力学耗散中的黏性耗散可直接表述为应变率相关的黏性应力。需要注意的是, 由于强冲击波阵面和测量得到的速度波剖面结构很难分辨, 因此, 这里的黏性耗散讨论仅适用于斜波/弱冲击加载, 对强冲击加载不适用。
不考虑弹塑性问题, 应力可分为压力和黏性力两部分之和, 即
$\sigma=p+Q$ 压力需要由状态方程进行计算, 这里选用Grüneisen状态方程
$p-p_{\mathrm{H}}=(\gamma / v)\left(e-e_{\mathrm{H}}\right)$ 式中:γ是Grüneisen系数; 下标“H”表示Hugoniot状态, 为参考状态
$ p_{\mathrm{H}}=p_{0}+\frac{c_{0}^{2}\left(v_{0}-v\right)}{\left[v_{0}-\lambda\left(1-v / v_{0}\right)\right]^{2}} $ $e_{\mathrm{H}}=e_{0}+\left(p_{\mathrm{H}}+p_{0}\right)\left(v_{0}-v\right) / 2$ 式中:c0、λ为Hugoniot系数。对于金属铝有c0=5.33 km/s, λ=1.338。黏性力选用应变率相关形式
$ Q=\eta \sqrt{\dot{\varepsilon}} $ 式中:
为应变率, η为黏性系数。对铝样品, 一般取η =0.1 MPa·s1/2, 进一步可观察材料的非等熵性。
使用比熵s和应变率
作为状态变量, 比内能e等物理量都可以进行如下表述
$ \dot{e}=(\partial e / \partial \varepsilon)_{s} \dot{\varepsilon}+(\partial e / \partial s)_{\varepsilon} \dot{s} $ 由热力学公式有
$(\partial e / \partial \varepsilon)_{s}=p / \rho_{0},(\partial e / \partial s)_{\varepsilon}=T$ 式中:ρ0=1/v0, 为材料密度; T为温度。故比内能的变化可写为
$ \rho_{0} \dot{e}=p \dot{\varepsilon}+\rho_{0} T \dot{s} $ 由于能量守恒方程还具有形式
$ \rho_{0} \dot{e}=p \dot{\varepsilon}+Q \dot{\varepsilon} $ 因此, 可得由黏性力和应变率给出的熵增变化情况
$ \rho_{0} T \dot{s}=Q \dot{\varepsilon} $ 同样地, 还可以得到由应变率给出的压力和温度的变化公式
$ \begin{aligned} & \dot{p}=B \dot{\varepsilon}+\rho T \gamma \dot{s}\end{aligned} $ $ \begin{aligned} & \dot{T}=T \gamma_{0} \dot{\varepsilon}+T \dot{s} / c_{n} \end{aligned} $ 式中:B=-v0(∂p/∂v)s, 是体积模量; cv为比热容, 由下式计算
$c_{v}=3 N_{\mathrm{A}} k+a\left(v / v_{0}\right)^{b} T$ 式中:NA、k分别为Avogadro常数和Boltzmann常数; 对于铝样品, 系数a=1.38 mJ/(mol·K2), b=1.8。
下面以Smith等人的激光驱动台阶Al靶实验数据[5]为例, 分析黏性耗散的影响。台阶靶厚度分别为10、20和30 μm, 后面贴有LiF透明窗口, 加载应力波在10 ns内上升到近100 GPa, 应变率峰值达到108 s-1。图 1(a)所示为不同靶厚条件下, Smith等人在Al/LiF界面处测量得到的速度剖面。为了对黏性力有直观的认识, 图 1(b)给出了10 μm厚的Al靶内3个不同位置(x=0, 5, 10 μm)处计算得到的黏性力剖面, 剖面峰值对应于应变率最高的位置, 约1 GPa的黏性力占该处总应力的5%~10%。本节计算均由反积分预处理得到的加载面应力作为加载条件, 进行正积分计算得到结果。
图 1(b) 计算所得10 μm Al靶中不同位置处的黏性力
Figure 1(b). Calculated viscous stress at 3 different positions of the 10 μm Al target in the laser driving ICE
黏性力对各物理量计算结果的影响如图 2所示。考虑黏性力后, 速度和应力剖面都比不考虑黏性力时低, 拉格朗日声速也较低, 尤其在应变率最高处(约26 ns), 两种声速值的差异达到2%, 从图 2中可以清楚地看到黏性耗散带来的影响。速度剖面中, 在27.6 ns时出现两条曲线相交的情况, 这是由于该处的应变率曲线呈平台状, 黏性力引起的应力增幅达到一个峰值, 从而导致速度增幅的较大跨越。在应力-应变曲线中, 考虑黏性力后, 应力-应变曲线较高, 同一应变下的应力最大变化为7%, 和实验测量应力-应变曲线之间的距离缩短了55%, 黏性耗散带来的熵增在此有所体现。应变率相关的黏性力对认识材料波形也非常重要, 尤其在应变率较高的地方; 随着压力的升高, 黏性力占总应力的比值不断下降, 其影响也越来越弱。总体来说, 和实验直接获取拉格朗日声速的不确定度0.3%(其中:测速不确定度2%, 时间不确定度0.3%, 样品厚度不确定度10%)[6]相比, 应变率相关的黏性带来的影响不容忽略。
图 2 激光驱动加载10 μm Al靶实验中,考虑黏性力前、后Al/LiF界面的应力、速度、拉格朗日声速及应力-应变曲线
Figure 2. Profiles of stress, particle velocity, Lagrangian sound speed and stress-strain relation at the Al/LiF interface in the laser driving 10 μm Al ICE with or without viscous force
下面分析黏性力大小的影响。根据20 μm Al靶的实验数据, 对不同位置处的物理量进行计算和分析。图 3(a)给出了计算得到的4个不同拉格朗日位置(x=5, 10, 15, 20 μm)处的速度剖面, 峰值速度约为3 km/s。对应的压缩比和应变率的情况在图 3(b)和图 3(c)中给出。由于应力波的追赶, 位置越靠近Al/LiF界面, 速度剖面越陡, 上升沿时间越短, 应变率越大, 相应的黏性力也越大, 达到相同速度值的压缩比越大。图 3(d)是应力-应变曲线、Hugoniot线及等熵线的比较, 可以看出, 上升沿时间越短, 应力越高。在应变率峰值阶段, 应力-应变曲线高于Hugoniot线和等熵线。等熵线始终处于应力-应变曲线下方, 这也是熵增的一种体现。图 3(e)为温度随压缩比的变化情况, 上升沿时间越短, 对应的温升越高, 相比于等熵情况, 温升高出近180 K, 但远小于冲击加载情况(3 000 K)。图 3(f)为熵增情况, 同样地, 上升沿时间越短, 黏性力越大, 熵增也越大, 为102 J/(kg·K)量级, 远小于冲击时的103 J/(kg·K)。
图 3 20 μm Al靶内不同位置处的计算结果
Figure 3. Calculated results at 4 different Lagrangian positions of the 20 μm Al target in the laser driving ICE
在x为5和20 μm两个拉格朗日位置处, 最大应变率差异为51 μs-1, 反映到黏性力中, 差异为715 MPa, 占当地总应力的5%;反映到温度变化(13 K)中, 则占温度变化总值的10%;反映到熵增中, 差异为16.65 J/(kg·K), 占变化总量的18%。由此可见, 在107 s-1量级的应变率范围内, 应变率的微小差异都能带来温度和熵增的显著变化。当然, 随着压力的升高, 应变率减小, 该效应也逐渐减弱。
上面的激光加载实验应变率约为108 s-1, 而磁压驱动等熵加载实验的应变率通常只有105~106 s-1, 因此黏性力也会大幅度减小。下面考虑一发CQ-4装置[7]上自由面Al靶的等熵加载实验。Al靶厚1.012 mm, 在700 ns内被加载到17 GPa。图 4所示为实验测量的样品自由面速度剖面以及计算得到的加载压力。图 5给出了5个不同拉格朗日位置(x=0, 0.2, 0.4, 0.6, 0.8 mm)处应变率和熵增的情况。由于存在自由面卸载, 离自由面越近, 应变率越低, 压缩比越小, 熵增自然也越小; 加载面峰值处的熵增低至8 J/(kg·K)以下, 而对应的冲击熵增为140 J/(kg·K)。
图 4 CQ-4等熵加载1.012 mm Al靶的实验和计算结果
Figure 4. Measured and calculated results of the 1.012 mm Al on CQ-4
-
在动高压加载、尤其是冲击加载实验中, 由于较高的应变率和温升, 热传导引起的耗散在总应力中占有一定比重, 直接影响材料压力-比容曲线的判定。为了体现正-反积分方法对抛物型方程的计算能力, 也为了避免热传导引起的耗散被力学耗散湮灭(热传导时间相比于斜波加载时间足够长), 不考虑黏性影响, 仅计算热传导项。对于考虑热传导后的控制方程, 如果不由应变率出发计算温度, 还可以在控制方程中直接加入温度计算公式。根据热力学方程有
$\mathrm{d} e=c_{v} \mathrm{d} T+\left(T \gamma_{\rho c_{v}}-p\right) \mathrm{d} v$ 积分后可得到温度计算公式
$ T=T_{i} \exp \left[\frac{p}{T_{i} c_{v}}\left(v-v_{i}\right)+\frac{e-e_{i}}{T_{i} c_{v}}+\gamma_{0}\left(\frac{v_{i}-v}{v_{0}}\right)\right] $ 式中:下标i表示上一个已知点的信息。相比于等熵温度计算公式Ts=T0exp[γ0(1-v/v0)], 增加的两项主要体现在非等熵效应上。
为了观察到较明显的温度差异, 仍选用应变率较高的Smith激光加载Al台阶靶数据[5]为例, 计算分析考虑热传导(导热系数为204 W/(m·K))后温度变化的影响, 并重点分析在偏应力计算中温度对本构模型的影响。采用修正的率相关SCG模型参与计算
$Y=\left[Y_{T}\left(\dot{\varepsilon}_{p}, T\right)+Y_{\mathrm{A}}\left(1+\beta \varepsilon_{p}\right)^{n}\right] G(p, T) / G_{0}$ $G(p, T)=G_{0}\left[1+\frac{G_{p}^{\prime}}{G_{0}} \cdot \frac{p}{\sqrt[3]{v_{0} / v}}-\frac{G_{T}^{\prime}}{G_{0}}(T-300)\right]$ $\dot{\varepsilon}_{p}=\left\{\frac{1}{C_{1}} \exp \left[\frac{2 U_{k}}{k T}\left(1-\frac{Y_{T}}{Y_{p}}\right)^{2}\right]+\frac{C_{2}}{Y_{T}}\right\}^{-1}$ 式中:
表示塑性应变率; 初始屈服强度和初始剪切模量分别为Y0=0.29 GPa, G0=27.1 GPa; 加工硬化函数(1+βε)n中,
表示屈服强度热激发部分的贡献, Peierls应力Yp=1.0 GPa, 二者通过塑性应变率((29)式)相关联; 2Uk=0.62 keV; k是Boltzmann常数; C1=0.71 μs-1, C2=0.012 MPa·s。
经反积分预处理得到加载情况, 将其作为正积分计算的边界条件; 在正积分计算中考虑热传导后, 得到的Al/LiF界面处的速度剖面, 计算结果和实验测量结果的比较在图 6中给出。
图 6 激光加载实验中考虑热传导后计算所得Al/LiF界面速度与测量结果[5]的对比
Figure 6. Comparison of calculated and measured velocity profiles at the Al/LiF interface in the laser driving ICE
下面分析考虑热传导前、后, 温度、屈服强度及剪切模量的差别(见图 7~图 9)。以10 μm Al靶的计算结果为例, 考虑热传导后, 温度峰值比未考虑热传导时降低了约40 K; 该温升差别也反应在SCG模型中, 屈服强度和剪切模量在其峰值处分别升高了0.014和0.8 GPa, 各占总屈服强度和总剪切模量的比值约为0.5%。
图 7 激光加载实验中考虑热传导前、后,计算所得Al/LiF界面处的温度
Figure 7. Comparison of temperature at the Al/LiF interface with or without heat conduction in the laser driving ICE
图 8 激光加载实验中考虑热传导前、后,计算所得材料的屈服强度
Figure 8. Comparison of yield strength with or without heat conduction in the laser driving ICE
图 9 激光加载实验中考虑热传导前、后,计算所得材料的剪切模量
Figure 9. Comparison of shear modulus with or without heat conduction in the laser driving ICE
图 10~图 12给出了Smith等人采用Lagrange分析方法得到的应力-应变曲线, 与之相比可以看出, 当选取Y0=0.29 GPa时, 我们计算得到的应力-应变曲线略低。在没有更好的途径能反映高应变率影响的情况下, 我们考虑通过直接增加Y0值来弥补应变率的影响。增加Y0至1.0 GPa后, 应力-应变曲线的确更靠近Smith等人的结果。图 11和图 12给出了相应的拉格朗日声速、屈服强度和偏应力的情况, 可以看出, 当Y0为0.29和1.0 GPa时, 对应的弹性屈服极限分别为1.0和4.2 GPa。这些结果说明目前采用的本构模型[8-9](10-4~106 s-1)还不能很好地表述高应变率(108 s-1)的情况。此外, 初始屈服强度的变化虽然可以直接影响动态屈服强度, 并通过2Y/3影响偏应力的计算, 但它仍是一种相对简单的形式, 还需要进一步探寻应变率与屈服强度以及应变率与偏应力的关系。
图 10 采用不同初始屈服强度的率相关SCG模型计算得到的应力-应变曲线
Figure 10. Stress-strain curve by SCG model with different initial yield strengths
-
通过反向积分和正向积分相结合的计算, 对准等熵压缩实验中广泛存在的力学耗散和熵增问题进行了讨论, 通过设置应变率相关的黏性力, 实现了高应变率加载下材料的力学耗散计算。对激光加载金属铝实验, 加载应力波在10 ns内上升到近100 GPa, 应变率峰值为108 s-1, 应变率引起的温升差异约为180 K, 熵增差异约为250 J/(kg·K); 对磁压驱动加载金属铝实验, 加载应力波在700 ns内上升到17 GPa, 较低的应变率(约105 s-1)带来的熵增低至8 J/(kg·K)以下。并对准等熵压缩实验中材料热传导引起的热耗散进行了分析, 激光加载金属铝实验中, 热传导导致温度降低了约40 K, 该温度差反映在率相关的SCG本构模型中, 导致屈服强度和剪切模量的峰值分别增加约0.014和0.8 GPa, 占总屈服强度和总剪切模量的0.5%左右。
准等熵压缩实验中金属铝的热力学耗散分析
Mechanical and Thermal Dissipation Analysis of Aluminum in Quasi-Isentropic Compression
-
摘要: 为获取高压下材料的纯热力学压力-比容参考线和完全物态方程,减去应力-应变曲线中的其它信息,对准等熵压缩实验中由加载应变率引起的黏性耗散和热传导引起的热耗散做了分析讨论。基于反积分计算和流体动力学积分计算相结合的方法,根据激光加载(约108 s-1)和磁驱动准等熵压缩(约105 s-1)的实验数据,对材料声速、应力-应变曲线、温度和熵增等物理量进行计算,分析了不同应变率与该物理量的关系;还对热传导和SCG本构模型进行了计算,分析了热传导引起的温度变化对材料屈服强度、剪切模量和拉格朗日声速的影响。结果表明:激光加载实验中,应变率引起的温升差异约为180 K, 熵增差异约为250 J/(kg·K),热传导导致温度下降40 K; 磁驱动准等熵压缩应变率较低,引起的熵增变化小于8 J/(kg·K)。Abstract: To study the pressure-specific volume (p-v) reference line and equation of state from the stress-strain curve of material at high pressure, the viscous dissipation due to loading strain rate and thermal dissipation due to irreversible heat conduction in quasi-isentropic compression experiment (ICE) were discussed and analyzed.A backward integration and forward integration method was used to analyze the data in laser driving and magnetic pressure driving ICE with different strain rates.For viscous dissipation, the sound speed, stress-strain curve, temperature and entropy production during loading were obtained, and the relations between high strain rate and these physical quantities were discussed.For thermal dissipation, through the calculation of the thermal conduction and SCG constitutive model, the variation of temperature and the corresponding yield strength, shear module and sound speed were presented.The results show that:in the laser driving ICE with a high strain rate of 108 s-1), the temperature variation caused by high strain rate is about 180 K, and the entropy production due to heat conduction is about 250 J/(kg·K); and in the magnetic pressure driving with a relatively low strain rate of 105 s-1, the entropy production is less than 8 J/(kg·K).
-
Key words:
- backward integration method /
- viscous force /
- strain rate /
- thermal conduction
-
图 6 激光加载实验中考虑热传导后计算所得Al/LiF界面速度与测量结果[5]的对比
Figure 6. Comparison of calculated and measured velocity profiles at the Al/LiF interface in the laser driving ICE
-
[1] Seagle C T, Davis J P, Martin M R, et al. Shock-ramp compression: Ramp compression of shock-melted tin[J]. Appl Phys Lett, 2013, 102(24): 244104. doi: 10.1063/1.4811745 [2] Davis J P. Experimental measurement of the principal isentrope for aluminum 6061-T6 to 240 GPa[J]. J Appl Phys, 2006, 99: 103512. doi: 10.1063/1.2196110 [3] 张红平, 王桂吉, 李牧, 等.准等熵压缩下金属钽的屈服强度分析[J].高压物理学报, 2011, 25(4): 321-326.
Zhang H P, Wang G J, Li M, et al. Yield strength analysis of tantalum in quasi-isentropic compression[J]. Chinese Journal of High Pressure Physics, 2011, 25(4): 321-326. (in Chinese)[4] Ding J L, Asay J R. Material characterization with ramp wave experiments[J]. J Appl Phys, 2007, 101(7): 073517. doi: 10.1063/1.2709878 [5] Smith R F, Eggert J H, Jankowski A, et al. Stiff response of aluminum under ultrafast shockless compression to 110 GPa[J]. Phys Rev Lett, 2007, 98(6): 065701. doi: 10.1103/PhysRevLett.98.065701 [6] Davis J P. Charice 1.0: An IDL application for characteristics-based inverse analysis of isentropic compression experiments, SAND 2007-4984[R]. Albuquerque: Sandia National Laboratories, 2007. [7] Wang G J, Zhao J H, Zhang H P, et al. Advances in quasi-isentropic compression experiments at institute of fluid physics of CAEP[J]. Eur Phys J Special Topics, 2012, 206(1): 163-172. doi: 10.1140/epjst/e2012-01597-y [8] Steinberg D J, Lund C M. A constitutive model for strain rates from 10-4 to 106 s-1[J]. J Appl Phys, 1989, 65(4): 1528-1533. doi: 10.1063/1.342968 [9] Steinberg D J. A rate-dependent constitutive model for molybdenum[J]. J Appl Phys, 1993, 74(6): 3827-3831. doi: 10.1063/1.355316 -