液氘单次冲击压缩的QMC模拟研究

李名锐; 周刚; 李志康; 耿宝刚; 范如玉; 赵南

doi:10.11858/gywlxb.2015.01.001

液氘单次冲击压缩的QMC模拟研究

doi: 10.11858/gywlxb.2015.01.001

李名锐^{1, 2,},
周刚^2,*,,
李志康²,
耿宝刚²,
范如玉^{1, 2},
赵南²

1.
清华大学工程物理系，北京 100084
2.
西北核技术研究所，陕西西安 710024

详细信息

作者简介:
李名锐(1983—)，男，博士研究生，主要从事稠密液氢模拟研究.E-mail limrui_nint@163.com

通讯作者:
周刚(1964-), 男, 博士, 研究员, 主要从事爆炸力学与冲击力学研究.E-mail:gzhou@nint.ac.cn

中图分类号: O521.2
计量
- 文章访问数: 6411
- HTML全文浏览量: 2248
- PDF下载量: 238
出版历程
- 收稿日期: 2013-03-20
- 修回日期: 2014-05-07

Single Shock Compression of Fluid Deuterium by QMC Simulation

LI Ming-Rui^{1, 2
,},
ZHOU Gang^{2,*
,},
LI Zhi-Kang²,
GENG Bao-Gang²,
FAN Ru-Yu^{1, 2},
ZHAO Nan²

1.
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
2.
Northwest Institute of Nuclear Technology, Xi'an 710024, China

摘要

摘要: 液氢在较宽压力与温度范围内具有复杂的物理特性，须从分子层面构建精确模型开展探索研究。利用电子-核耦合的CEIMC法模拟液氘单次冲击实验，分析了液氘冲击特性，当压力达50.3 GPa时液氘具有最大压缩率4.48，在110 GPa冲击压力附近未发现有压缩率急剧增大的迹象。选用合适的Al基板材料模型，建立了液氘单次压缩状态与实验条件间的关系，总结了单次冲击实验规律。得到的状态方程与现有动高压实验结果一致，也与经修正后的100 GPa以上压力的Omega激光实验值吻合，说明采用基于共振价键理论的波函数后，CEIMC法可应用于液氘的冲击模拟。
- 单次冲击压缩 /
- 液氘 /
- Hugoniot曲线 /
- QMC法
Abstract: Since liquid hydrogen has some complex physical behaviors in a relatively wide range of pressure and temperature, one more accurate and appropriate ab-initio simulation method is urgently required.We used the Coupled Electron-Ion Monte Carlo (CEIMC) method to simulate the single shock compression experiment of deuterium, and systematically studied its thermodynamic properties under high pressure.It is shown that the liquid deuterium principle Hugoniot has a maximum compression ratio of 4.48 at about 50.3 GPa, but the compression ratio does not increase gradually around 110 GPa.By utilizing an appropriate aluminum impedance match model, we have established the correlations between the states of compressed deuterium and the shock or particle velocities in aluminum taken before the shock wave crossing the Al-D₂ interface, and contemporarily summarized the physical processes of liquid deuterium under single shock compression in detail.The calculation results are not only in very good agreement with the existing experiments via different high-pressure generating technologies, but also consistent with the data above 100 GPa corrected from the Omega Laser driven experiment.Therefore, the CEIMC method combined with the wave function based on the resonant valence bond theory is very suitable for the simulations of shocked liquid deuterium.
- single shock compression /
- liquid deuterium /
- Hugoniot curve /
- QMC method

HTML全文

1. 引言

现代动高压冲击压缩技术是实验室内获取物质高温高压极端状态的一种常用手段。研究液氢及其同位素的冲击特性在现代天体物理、惯性约束核聚变(ICF)以及氢金属化理论等方面具有重要的科学意义与应用价值, 长期以来一直是实验与理论研究需解决的关键问题之一^[1]。1995年, Nellis等^[2]利用二级轻气炮多次冲击压缩法, 在140 GPa、3 000 K条件下首次证实了液态金属氢的存在。之后, 人们采取强激光^[3-6]、Z装置^[7-8]、球形内聚爆轰^[9-11]等多种动高压设备研究了液氘的单次冲击压缩特性, 但不同加载手段的结果间仍有明显差异, 如Nova激光单次压缩液氘的最大压缩率为6.0^[3-5], Z装置最大压缩率为4.3^[8], Omega激光发现液氘压缩率在110 GPa左右会急剧增大到5.0^[6]。目前100 GPa以上压力区内的实验数据仍具争议且较为稀缺, 许多实验现象及规律还需从理论上系统分析获得。

理论研究方法主要有化学模型法和ab-initio法两大类。尽管化学模型法如Sesame^[12]、Ross^[13]、Saumon-Chabrie^[14]、FVT^[15]等计算速度快, 且在低压区内能获得与气炮实验吻合较好的冲击结果, 但由于无法描述出电子的束缚态和激发态, 且高温高压下预先假定的粒子间有效势可能会失效, 化学模型法在100 GPa以上压力范围内的适用性较差。基于第一性原理的ab-initio法, 如DFT^[16]、QMD^[17-18]、PIMC^[19]等是研究高温稠密物质状态方程较为精确的方法, 但目前大部分研究集中在对100 GPa以下压力区的模拟, 对更高压力范围的冲击实验研究较少, 不同方法得到的液氘最大压缩率也存在差异, 且许多ab-initio法因在低压区内不能预估出与轻气炮、Z装置、强激光等实验一致的结果, 使精度受到了质疑, 如赝势及近似交换-关联函数等限制了DFT法在极端高压条件下的应用, 正比于温度倒数的时间步长降低了PIMC法在10 000 K以下低温区的模拟精度等。为考虑核的剧烈运动效应, 提高计算精度, Ceperley等^[1]提出了电子-核耦合的MC方法(Coupled Electron-Ion Monte Carlo Method, CEIMC), 并成功预测了液氢的电导率、自由能、分子离解等, 发现CEIMC法比较适用于高温高压液氢系统, 认为该方法是当前ab-initio MC模拟中基于能量计算最为精确的方法, 可作为PIMC法在低温区的应用互补, 但至今仍未应用于冲击实验的模拟中。为此, 本研究将利用CEIMC法来模拟液氘单次冲击实验, 在经受住现有实验结果考核的前提下, 预估液氘的最大压缩率, 探寻100 GPa以上压力范围内是否存在新相致使液氘压缩率突变至5.0, 分析离解及电离等效应对冲击特性的影响。选用合适的Al基板材料模型, 建立液氘单次压缩状态与实验条件间的关系, 总结单次冲击实验规律。

2. ab-initio模拟方法

为获得液氘在较宽压力温度范围内的精确状态数据, 可从两个方面入手:(1)采用精度较高的ab-initio模拟算法, (2)构建能较好描述高压液氢特性的电子波函数。本研究采用西北核技术研究所提供的程序进行模拟计算, 其模拟算法为基于能量计算最为精确的且考虑核运动的CEIMC法^[1]。波函数是在共振价键理论基础上, 由Jastrow项与反对称孪生函数乘积AGP项(Jastrow Antisymmetrized Geminal Product, Jastrow+AGP)构成的JAGP波函数, 该波函数能描述多种结构的共振键, 考虑了更多电子间关联能, 且满足大小一致性, 具有较高的计算精度, 其一般形式可参见文献[20]。

原则上讲, 波函数应包含所有可能的原子轨道, 但过于复杂计算任务也十分庞大, 需在效率与精度间进行平衡。本研究将针对每个原子仅选用一个轨道形式函数作为基组来构造JAGP波函数, 以便能减少优化参数数目, 节省计算机时, 而这个轨道形式函数由电子可能占据主要轨道的线性组合表示出。对AGP项, 采取Staler型函数轨道(STO)基组1s、2s、2p、3s轨道的线性组合

$\begin{aligned} \varphi_{s}(r) &=C_{1 s}\left[\mathrm{e}^{-z_{1} r}+C_{2 s} r \mathrm{e}^{-z_{2} r}+C_{3 s} r^{2} \mathrm{e}^{-z_{3} r}+C_{2 p}(x+y+z) \mathrm{e}^{-z_{4} r}\right] \\ &=C_{1 s}\left[\mathrm{e}^{-z_{1} r}+p_{1} r \mathrm{e}^{-z_{2} r}+p_{2} r^{2} \mathrm{e}^{-z_{3} r}+p_{3}(x+y+z) \mathrm{e}^{-z_{4} r}\right] \end{aligned}$

(1)

式中:C_1s为归一化因子, z₁、z₂、z₃、z₄、p₁、p₂、p₃为待优化参数。

对三体Jastrow项, 采取Gaussian型轨道函数(GTO)基组

$\varphi_{g}(r)=\mathrm{e}^{-x r^{2}}$

(2)

为改善波函数的结点面, 当满足r_s < 1.75或T > 15 000 K时在构建的JAGP波函数中引入Backflow项^[21]。采取正则系综, 模拟盒内氢原子数固定N=64, 采取经典理论处理核的运动, H₂-H₂分子间相互作用选取Silvera-Goldman势^[22], 氢分子内H-H作用势选用Kolos-Wolniewicz势^[23]。由于JAGP波函数参数量级与选取的轨道基函数N_orb成N_orb(N_orb-1)/2关系, 优化计算量十分庞大, 为兼顾计算精度与效率, 采用仅需求解一阶导数的基于随机游走的GBRW法来优化能量及其方差组合, 即min{ξ〈E〉+(1-ξ)〈σ²〉}, 取ξ=0.95以混入少许比例的能量方差值。采用周期性边界条件(PBC)和扭转平均边界条件(TABC)^[1], 其中根据不同温度及密度条件等分80~300份扭转角。前期运行2 000~3 000步用于平衡计算, 之后依据密度和温度进行20 000~30 000步的统计计算, 且每隔40步将数据输出用于分析。每运行1步核计算步将随机处理(移动、平动或转动)4个H₂分子, 每运行1步电子计算步将重复10个参数优化步, 其中以处理500个参数为一循环, 直至依次循环优化完所有参数方构成1个参数优化步。

在1.75≤r_s≤2.20、2 000 K < T < 70 000 K范围内计算出148个液氘状态点。在此基础上, 先固定液氘密度, 依据多项式p(E)=α₀+α₁T+α₂T²+α₃T³分段拟合出压力p、能量E随温度T的变化关系, 然后以密度作为变量进行插值, 得到任意(T, ρ)点的状态值(p, E)。由于CEIMC法计算选取的密度点较多, 密度间隔不大, 且在讨论的密度、温度范围内压力与能量随密度均未出现突变的现象, 从后续分析可以证实该处理方法可行。

3. 液氘单次冲击特性

冲击作用下, 液氘在波阵面上满足动量、质量及能量守恒的Rankine-Hugoniot关系

$p-p_{0}=\rho_{0}\left(D-u_{\mathrm{p}}\right)\left(u_{\mathrm{p}}-u_{0}\right)$

(3)

$V=V_{0}\left[1-\left(u_{\mathrm{p}}-u_{0}\right) /\left(D-u_{\mathrm{p}}\right)\right]$

(4)

$E-E_{0}=\left(p+p_{0}\right)\left(V_{0}-V\right) / 2$

(5)

式中:p₀、ρ₀、E₀、V₀、u₀分别代表冲击前液氘的压力、密度、能量、体积及粒子速度; p、ρ、E、V、u_p、D分别代表单次冲击后液氘的压力、密度、能量、体积、粒子速度及冲击波速度; V₀=1/ρ₀, V=1/ρ。利用CEIMC法计算出液氘状态点(p, ρ, E)后, 结合上述波阵面守恒关系, 便可以分析液氘的冲击特性。

图 1是不同实验及计算方法获得的液D₂单次冲击压缩结果, 液氘的初始条件^[19]为ρ₀= 0.171 g/cm³, T₀=20 K, E₀=-15.886 eV, p₀=0。在100 GPa以下低压区内, 本研究计算出的冲击Hugoniot曲线位于Desjarlais^[17]及Collins^[24]的QMD结果之间, 在50.3 GPa压力处具有最大压缩率η_max=4.48, 与Z装置等现代动高压装置测得的最大压缩率η_max=4.25~4.50相符^[7-8-25]。尽管在100~120 GPa压力范围内, 本研究曲线落在了球形内聚爆轰实验误差棒内^[9-11], 但在100 GPa以上的高压区域内, 未发现液氘压缩率如同Hicks^[6]激光实验结果一样出现突然增加的迹象, 相反液氘压缩率会随压力增加而逐渐减小。当p=100 GPa时, 压缩率η=4.37, p=240 GPa时, 压缩率η=4.22, 最终将收敛至理想气体的压缩率η=4.0。Knudson等^[26]针对石英(Quartz)进行了压力至1.6 TPa的冲击压缩实验, 重新校准获得了精度更高的石英压力标准模型, 利用该模型修正Hicks激光实验后发现, Hicks实验结果将整体左移变硬, 且大部分测量数据的误差棒都落在了本研究预估曲线上, 在100~120 GPa压力范围内液氘压缩率突然增加的现象消失。整体而言, 本研究给出的液氘Hugoniot曲线位于多种化学模型及ab-initio结果的中间某位置, 且基本通过了所有加载类型实验的误差棒, 特别是通过了目前仅有的单次冲击压力高于100 GPa的实验误差棒(即修正后的Hicks激光实验结果)。因此, 本研究方法具有较高的可信度, 较其他ab-initio方法而言, 的确能改善冲击特性的模拟精度, 基于共振价键理论的JAGP波函数, 可将CEIMC法的应用范围从Ceperley的10 000 K以下低温区^[1]推广至70 000 K的高温区内, 可应用于动高压冲击实验的模拟研究, 同时也说明在100 GPa以上高压区内亟需开展Z装置或球形内聚装置实验以验证模型的正确性。

图 1 液氘冲击Hugoniot曲线的实验及计算结果

Figure 1. Experimental and calculation results for deuterium Hugoniot curve

下载: 全尺寸图片幻灯片

在低于40 GPa的压力区内, 液氘单次冲击Hugonoit曲线变化缓慢, 诸多曲线及数据点压缩在一起不便于分析。为考核本方法在该区域内的适应性, 在轻气炮实验^[25]范围内重新获得了液氘单次冲击压力及温度曲线, 如图 2所示。

图 2 轻气炮范围内液氘压力、温度与密度关系

Figure 2. Pressure and temperature versus density for deuterium in the gas gun range

下载: 全尺寸图片幻灯片

由图 2可以看出, 本研究结果较Lenosky等^[16]的DFT结果更接近实验值, 与实验相比, 压力最大偏差约2.4 GPa, 温度值最大偏差约600 K, 考虑到轻气炮实验温度测量不确定度为100~200 K, 体积与压力不确定度约为1%~2%, 本研究结果与实验值吻合较好。进一步分析计算值与实验值之间的差别, 认为在40 GPa以下低压区内获得的冲击温度(或能量)略低于实验值的现象, 可能是本方法未考虑液氘零点振动能的缘故。

在p≈20 GPa、ρ≈0.581 7g/cm³、T≈4 000 K处为冲击Hugoniot曲线的转折点, 之后压力及温度向下弯曲变软, 这是液氘在20 GPa时发生离解吸收能量的缘故, 与Bonev等^[27]的QMD结果(ρ≈0.587 g/cm³、T≈4 500 K)及其他研究者^[17]的结果相一致。随着温度升高, 液氘的离解度变大, 数目增多的自由粒子会使系统压力变大, 在p≈30 GPa、ρ≈0.726 3 g/cm³、T≈5 440 K后液氘冲击压力及温度曲线开始向上弯曲, 冲击压力随温度近似成线性关系(如图 3所示), 说明此时粒子热动能占主导因素, 高温液氘逐渐表现出理想气体的特性。在稍过最大压缩率η_max=4.48(对应p≈50.3 GPa、T≈10 000K)的状态点p≈53 GPa、ρ≈0.765 g/cm³、T≈11 000 K后, 冲击压力、温度及能量将随压缩率成近直线上升, 如图 4所示, Desjarlais^[17]利用QMD法认为此时有90%的液氘分子离解。

图 3 液氘压力、压缩率与冲击温度关系

Figure 3. Pressure and compression ratio versus temperature for D₂

下载: 全尺寸图片幻灯片

图 4 液氘温度、压力及能量与压缩率关系

Figure 4. Temperature, pressure and energy versus compression ratio for D₂

下载: 全尺寸图片幻灯片

图 5是单次压缩下液氘冲击波速度及压力随粒子速度的变化关系。可以看出, Ross模型^[13]预估的结果偏软, Sesame 72模型预估的结果偏硬^[12], 而本研究曲线与实验最为接近, 特别是在粒子速度低于约21 km/s或压力低于约105 GPa的范围内, 研究结果与实验值吻合得非常好。当粒子速度大于约21 km/s或压力大于约105 GPa时, 本研究结果较Hicks实验值^[6]略显偏硬, 但计算曲线仍过Hicks实验的大部分误差棒, 偏差不超过5%。

图 5 液氘冲击波速度及压力与粒子速度关系

Figure 5. Shock velocity and pressure versus particle velocity in D₂

下载: 全尺寸图片幻灯片

从图 1分析可知, 利用更为精确的石英模型^[26]修正Hicks实验后, Hicks激光实验结果会相应变硬左移而靠近本研究预估的Hugoniot曲线, 由此推测, 经修正后Hicks实验的冲击波速度及压力也会相应左移变硬而落在本研究预估曲线上。依据计算结果可拟合出液氘冲击波速度D(km/s)与粒子速度u_p(km/s)间的分段线性关系

$D=\left\{\begin{array}{ll} (3.1243 \pm 0.0041)+(1.0743 \pm 0.0005) \times u_{\mathrm{p}} & 7.2151 \leqslant u_{\mathrm{p}} \lt 9.1949 \\ (4.9769 \pm 0.0049)+(0.8725 \pm 0.0005) \times u_{\mathrm{p}} & 9.1949 \leqslant u_{\mathrm{p}} \lt 11.4769 \\ (-0.6035 \pm 0.0022)+(1.3273 \pm 0.0001) \times u_{\mathrm{p}} & 11.4769 \leqslant u_{\mathrm{p}} \lt 35.0126 \end{array}\right.$

(6)

图 6是计算得到的单次压缩下液氘粒子速度及冲击波速度随冲击温度的变化结果。当冲击温度低于70 000 K时, 液氘粒子速度和冲击波速度在T≈5 440 K、p≈30 GPa之前均随温度增大而迅速变大, 但在T≈11 000 K、p≈53 GPa之后, 增长率du_p/dT、dD/dT逐渐减小, 所以离解反应初期对液氘冲击波速度及粒子速度的影响较大, 而当离解完全后其影响将变得很小。倘若液氘在低于70 000 K冲击温度范围内先离解而后电离, 可以猜测D原子电离生成两自由粒子D⁺与e对冲击波速度及粒子速度的影响不及D₂离解生成两个D原子的大。

图 6 液氘粒子速度和冲击波速度与温度关系

Figure 6. Particle velocity and shock velocity versus temperature for D₂

下载: 全尺寸图片幻灯片

4. 液氘冲击状态与实验条件关系

尽管动高压的加载手段各异, 但当冲击波进入基板后其传播规律及对液氘样品的压缩规律基本相近, 计算原理也基本相同。在保证计算结果能与实验进行比较的前提下, 本研究选用金属Al做样品盒的基板材料。当基板冲击波初次到达基板与液氘界面时, 基板将会向液氘侧等熵膨胀, 向基板内反射一中心稀疏波对基板形成部分卸载, 同时透入液氘样品内一冲击波对液氘形成单次冲击压缩。由冲击波基础知识可知, 界面两侧的压力与粒子速度应保持连续, 界面两侧的状态点即为Al基板等熵卸载曲线与液氘单次冲击压缩曲线的交点, 而液氘的冲击Hugoniot曲线已求得, 因此, 该问题的核心便归结为求解Al基板的卸载曲线。

为提高模拟精度, 采用Hicks的建议, 即将Sesame 3717, 3715, 3719, Kerley’s3700及qEOS等Al材料标准模型方程进行平均作为计算所需的等熵卸载线^[6]。Al基板的初始密度2.740 g/cm³, 冲击波速度u_{s, Al}与粒子速度u_{p, Al}关系选用Hicks给出的最优拟合结果

$u_{\mathrm{s}, \mathrm{Al}}=\left\{\begin{array}{lr} (9.449 \pm 0.020)+(1.324 \pm 0.016) \times\left(u_{\mathrm{p}, \mathrm{Al}}-3.0220\right) & u_{\mathrm{p}, \mathrm{Al}} \lt 6.763 \\ (17.992 \pm 0.078)+(1.167 \pm 0.026) \times\left(u_{\mathrm{p}, \mathrm{Al}}-9.8381\right) & 6.763 \leqslant u_{\mathrm{p}, \mathrm{Al}} \lt 30 \end{array}\right.$

(7)

Al基板的Grüneisen参数γ_Al(V)可根据文献数据^[28]拟合成有关体积V的多项式

$\gamma_{\mathrm{Al}}(V) / V=\Gamma(V) / V \quad V_{1} \leqslant V \leqslant V_{2}$

(8)

式中:Γ(V)=0.657 76+1.173 18V+1.365 19V²+15.786 56V³, V₁=0.162 24 cm³/g, V₂=0.374 13 cm³/g, 在体积V₁、V₂两端区域将应用γ₀/V₀=γ/V近似。不失一般性, 将Al基板的粒子速度作为冲击实验的输入量, 而不纠结于加载手段是飞片撞击还是激光驱动。

图 7是单次压缩下液氘冲击波速度及压力随Al基板入射冲击波速度的变化结果。与图 5相类似, 在100 GPa以下的压力范围内, 本研究预估曲线与实验结果吻合得很好, 而当压力高于100 GPa时本研究结果较Hicks等^[6]激光实验值略显偏硬, 偏差也不超5%, 但未发现液氘冲击曲线有由硬变软的迹象。同样推测, 若利用更为精确的石英模型^[26]修正Hicks实验后, Hicks激光实验值也会左移变硬而落在本研究预估曲线上, 说明本研究选用的Al基板材料模型合适。

图 7 液氘冲击波速度及压力与Al基板入射冲击波速度关系

Figure 7. Shock velocity and pressure in D₂ versus incident shock velocity in Al

下载: 全尺寸图片幻灯片

拟合图 7(a)的计算结果可以得到单次压缩下液氘冲击波速度D(km/s)与基板入射冲击波速度u_{s, Al}(km/s)的关系

$D=\left\{\begin{array}{ll} (-4.8182 \pm 0.0087)+(1.4847 \pm 0.0008) \times u_{\mathrm{s}, \mathrm{Al}} & 10.7439 \leqslant u_{\mathrm{s}, \mathrm{Al}} \lt 11.8753 \\ (-1.5910 \pm 0.0078)+(1.2135 \pm 0.0006) \times u_{\mathrm{s}, \mathrm{Al}} & 11.8753 \leqslant u_{\mathrm{s}, \mathrm{Al}} \lt 13.4789 \\ (-10.4974 \pm 0.0018)+(1.8468 \pm 0.0001) \times u_{\mathrm{s}, \mathrm{Al}} & 13.4789 \leqslant u_{\mathrm{s}, \mathrm{Al}} \lt 30.4717 \end{array}\right.$

(9)

图 8是计算获得的单次冲击下液氘温度、能量、粒子速度及压缩率随Al基板入射粒子速度的变化结果。

图 8 液氘温度、能量、粒子速度及压缩率和Al基板粒子速度关系

Figure 8. Temperature, energy, particle velocity and compression ratio in D₂ versus initial particle velocity in Al

下载: 全尺寸图片幻灯片

单次冲击后液氘的粒子速度u_{p, D₂}(km/s)与Al基板的初始粒子速度u_{p, Al}(km/s)成线性关系

$u_{\mathrm{p}, \mathrm{D}_{2}}=(0.65646 \pm 0.01659)+(1.70713 \pm 0.00142) u_{\mathrm{p}, \mathrm{Al}}$

(10)

当Al基板初始粒子速度u_{p, Al}=8.4 km/s时, 液氘达最大压缩率η_max=4.48。依据图 7(b)结果及(7)式, 可以得到液氘经单次冲击后的压力p₁(GPa)随Al基板初始粒子速度u_{p, Al}(km/s)的变化关系

$p_{1}=7.09759-1.45465 \times u_{\mathrm{p}, \mathrm{Al}}+0.82798 \times u_{\mathrm{p}, \mathrm{Al}}^{2}-0.00473 \times u_{\mathrm{p}, \mathrm{Al}}^{3}$

(11)

式中:4.0≤u_{p, Al}≤18.0。液氘的冲击温度T及能量E随Al基板初始粒子速度成上凹曲线, 其增长率dT/u_{p, Al}、dE/u_{p, Al}在u_{p, Al}≈5 km/s(对应p≈20 GPa)与u_{p, Al}≈8.7 km/s(对应p≈53 GPa)之间即离解反应区内变化最快, 而在两端区域即液氘处于双原子分子态或单原子分子态时的变化比较平缓。

5. 结论

采用CEIMC法系统研究了液氘的单次冲击特性, 获得的冲击Hugoniot曲线基本通过了所有动高压加载方式实验的误差棒。当压力为20 GPa时Hugoniot曲线变软, 当压力超过30 GPa后认为液氘的热动能占主导, 当压力达到50.3 GPa时液氘具有最大压缩率4.48, 与实验值4.25~4.50相吻合, 在110 GPa冲击压力附近未发现有压缩率急剧增大的迹象。通过分析液氘单次冲击下的冲击波速度与粒子速度, 拟合出两者的函数关系式, 认为离解反应初期对液氘冲击特性的影响较大。撇开各类动高压加载手段, 选择合适的Al基板模型并以基板入射粒子速度为实验输入条件, 建立了液氘单次冲击状态与实验条件间的关系, 模拟结果在100 GPa以下压力区内与实验相吻合, 在100 GPa以上区域内较Hicks激光实验略显偏硬, 但相差不超过5%, 推测与经精确石英模型修正Hicks实验后的数据相比差距会更小。说明基于共振价键理论的JAGP波函数可将CEIMC法的应用范围从Ceperley的10 000 K以下低温区推广至70 000 K的高温区内, 可应用于动高压冲击实验的模拟研究, 但仍需100 GPa以上不同加载方式实验数据或其他精确理论数据的考核与验证。

图 1 液氘冲击Hugoniot曲线的实验及计算结果

Figure 1. Experimental and calculation results for deuterium Hugoniot curve

下载: 全尺寸图片幻灯片

图 2 轻气炮范围内液氘压力、温度与密度关系

Figure 2. Pressure and temperature versus density for deuterium in the gas gun range

下载: 全尺寸图片幻灯片

图 3 液氘压力、压缩率与冲击温度关系

Figure 3. Pressure and compression ratio versus temperature for D₂

下载: 全尺寸图片幻灯片

图 4 液氘温度、压力及能量与压缩率关系

Figure 4. Temperature, pressure and energy versus compression ratio for D₂

下载: 全尺寸图片幻灯片

图 5 液氘冲击波速度及压力与粒子速度关系

Figure 5. Shock velocity and pressure versus particle velocity in D₂

下载: 全尺寸图片幻灯片

图 6 液氘粒子速度和冲击波速度与温度关系

Figure 6. Particle velocity and shock velocity versus temperature for D₂

下载: 全尺寸图片幻灯片

图 7 液氘冲击波速度及压力与Al基板入射冲击波速度关系

Figure 7. Shock velocity and pressure in D₂ versus incident shock velocity in Al

下载: 全尺寸图片幻灯片

图 8 液氘温度、能量、粒子速度及压缩率和Al基板粒子速度关系

Figure 8. Temperature, energy, particle velocity and compression ratio in D₂ versus initial particle velocity in Al

下载: 全尺寸图片幻灯片

参考文献(28)

[1]	Morales M A, Pierleoni C, Ceperley M D. Equation of state of metallic hydrogen from coupled electron-ion Monte Carlo simulations[J]. Phys Rev E, 2010, 81(2): 021202. http://www.ncbi.nlm.nih.gov/pubmed/20365556
[2]	Nellis W J, Weir S T, Mitchell A C. Minimum metallic conductivity of fluid hydrogen at 140 GPa(1.4 Mbar)[J]. Phys Rev B, 1998, 59(5): 3434-3449.
[3]	da Silva L B, Celliers P M, Collins G W, et al. Absolute equation of state measurements on shocked liquid deuterium up to 200 GPa(2 Mbar)[J]. Phys Rev Lett, 1997, 78(3): 483-486. doi: 10.1103/PhysRevLett.78.483
[4]	Collins G W, da Silva L B, Celliers P, et al. Measurements of the equation of state of deuterium at the fluid insulator-metal transition[J]. Science, 1998, 281(1): 1178-1181. http://www.ncbi.nlm.nih.gov/pubmed/9712579
[5]	Collins G W, Celliers P M, da Silva L B, et al. Equation of state measurements of hydrogen isotopes on Nova[J]. Phys Plasmas, 1998, 5(5): 1864-1869. doi: 10.1063/1.872857
[6]	Hicks D G, Boehly T R, Celliers P M, et al. Laser-driven single shock compression of fluid deuterium from 45 to 220 GPa[J]. Phys Rev B, 2009, 79(1): 014112. doi: 10.1103/PhysRevB.79.014112
[7]	Knudson M D, Hanson D L, Bailey J E, et al. Equation of state measurements in liquid deuterium to 70 GPa[J]. Phys Rev Lett, 2001, 87(22): 225501. doi: 10.1103/PhysRevLett.87.225501
[8]	Knudson M D, Hanson D L, Bailey J E, et al. Principal hugoniot, reverberating wave, and mechanical reshock measurements of liquid deuterium to 400 GPa using plate impact techniques[J]. Phys Rev B, 2004, 69(14): 144209. doi: 10.1103/PhysRevB.69.144209
[9]	Belov S I, Boriskov G V, Bykov A I, et al. Shock compression of solid deuterium[J]. JETP Lett, 2002, 76(7): 433-435. doi: 10.1134/1.1528696
[10]	Boriskov G V, Bykov A I, Ilkaev R I, et al. Shock compression of liquid deuterium up to 109 GPa[J]. Phys Rev B, 2005, 71(9): 092104. doi: 10.1103/PhysRevB.71.092104
[11]	Grishechkin S K, Gruzdev S K, Gryaznov V K, et al. Experimental measurements of the compressibility, temperature, and light absorption in dense shock-compressed gaseous deuterium[J]. JETP Lett, 2004, 80(6): 398-404. doi: 10.1134/1.1830656
[12]	Kerley G I. Equation of state and phase diagram of dense hydrogen[J]. Phys Earth Planet Inter, 1972, 6(1): 78-82. http://www.sciencedirect.com/science/article/pii/0031920172900362
[13]	Ross M. Linear-mixing model for shock-compressed liquid deuterium[J]. Phys Rev B, 1997, 58(2): 669-677. http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=PRBMDO000058000002000669000001&idtype=cvips&gifs=Yes
[14]	Saumon D, Chabrier G. Fluid hydrogen at high density: Pressure ionization[J]. Phys Rev A, 1992, 46(4): 2084-2100. doi: 10.1103/PhysRevA.46.2084
[15]	田春玲, 经福谦, 顾云军, 等.高温高密度氢(氘)的物态方程——离解效应研究[J].高压物理学报, 2007, 21(1): 8-14. http://d.wanfangdata.com.cn/Periodical/gywlxb200701002 Tian C L, Jing F Q, Gu Y J, et al. Equation of state for fluid hydrogen and deuterium: Dissociation effects studies[J]. Chinese Journal of High Pressure Physics, 2007, 21(1): 8-14. (in Chinese) http://d.wanfangdata.com.cn/Periodical/gywlxb200701002
[16]	Lenosky T J, Bickham S R, Kress J D, et al. Density-functional calculation of the Hugoniot of shocked liquid deuterium[J]. Phys Rev B, 2000, 61(1): 1-4. http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=PRBMDO000061000001000001000001&idtype=cvips&gifs=Yes
[17]	Desjarlais M P. Density-functional calculations of the liquid deuterium Hugoniot, reshock, and reverberation timing[J]. Phys Rev B, 2003, 68(6): 064204. doi: 10.1103/PhysRevB.68.064204
[18]	Lenosky T J, Kress J D, Collins L A. Molecular-dynamics modeling of the hugoniot of shocked liquid deuterium[J]. Phys Rev B, 1997, 56(9): 5164-5169. doi: 10.1103/PhysRevB.56.5164
[19]	Militzer B, Ceperley D M. Path integral monte carlo calculation of the deuterium hugoniot[J]. Phys Rev Lett, 2000, 85(9): 1890-1893. doi: 10.1103/PhysRevLett.85.1890
[20]	Casula M, Attaccalite C, Sorella S. Correlated geminal wave function for molecules: An efficient resonating valence bond approach[J]. J Chem Phys, 2004, 121(15): 7110-7126. doi: 10.1063/1.1794632
[21]	Holzmann M, Ceperley D M, Pierleoni C, et al. Backflow correlations for the electron gas and metallic hydrogen[J]. Phys Rev E, 2003, 68(4): 046707. doi: 10.1103/PhysRevE.68.046707
[22]	Silvera I F, Goldman V V. The isotropic intermolecular potential for H₂ and D₂ in the solid and gas phase[J]. J Chem Phys, 1987, 69(9): 4209-4213. doi: 10.1063/1.437103
[23]	Kolos W, Wolniewicz L. Potential-energy curves for the, and C¹Π_u states of the hydrogen molecule[J]. J Chem Phys, 1965, 43(7): 2429-2441. doi: 10.1063/1.1697142
[24]	Collins L A, Bickham S R, Kress J D, et al. Dynamical and optical properties of warm dense hydrogen[J]. Phys Rev B, 2001, 63(18): 184110. doi: 10.1103/PhysRevB.63.184110
[25]	Nellis W J, Mitchell A C, Thiel V M, et al. Equation of state data for molecular hydrogen and deuterium at shock pressures in the range 2-76 GPa(20-760 kbar)[J]. J Chem Phys, 1983, 79(3): 1480-1486. doi: 10.1063/1.445938
[26]	Knudson M D, Desjarlais M P. Shock compression of quartz to 1.6 TPa: Redefining a pressure standard[J]. Phys Rev Lett, 2009, 130(22): 225501. http://www.ncbi.nlm.nih.gov/pubmed/20366104
[27]	Bonev S A, Militzer B, Galli G. Ab initio simulations of dense liquid deuterium: Comparison with gas gun shock wave experiments[J]. Phys Rev B, 2004, 69(1): 014101. http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=PRBMDO000069000001014101000001&idtype=cvips&gifs=Yes
[28]	Knudson M D, Lemke R W, Hayes D B, et al. Near-absolute Hugoniot measurements in aluminum to 500 GPa using a magnetically accelerated flyer plate technique[J]. J Appl Phys, 2003, 94(7): 4420-4431. doi: 10.1063/1.1604967