Superconductivity of Solid Hydrogen under Extreme Pressure
-
摘要: 氢元素在常压下具有最简单的晶体结构和物理性质。随着压强增加,氢单质发生相变,由绝缘体转变为金属,被称为金属氢。数值模拟表明,金属氢具有高温超导电性,因此,金属氢研究也被称为高压物理领域的“圣杯”课题。利用基于密度泛函理论的第一性原理计算方法,对固体氢在极端高压(0.5~5.0 TPa)下的结构和超导电性开展了系统研究。研究结果表明:固体氢的高压相变序列为I41/amd→oC12→cI16;对于同一种结构,随着压强增加,电声耦合系数减小,费米面处电子态密度减小,特征振动频率增加,超导转变温度发生小幅变化;在2.0 TPa压强下,固体氢的超导转变温度高达418 K(库伦赝势经验值μ*=0.10)。研究工作将为金属氢及其超导电性的后续理论和实验研究提供参考。Abstract: Hydrogen has the simplest crystal structure and physical properties at ambient pressure. As the pressure increases, hydrogen undergoes phase transition from insulator to metal, which being called metallic hydrogen. The numerical calculations also indicate that metallic hydrogen has high-temperature super-conductivity, thus the metal hydrogen is also known as the holy grail of physics subject. In this paper, the structural phase transition and superconducting transition temperature (Tc) of solid hydrogen under extreme high pressure 0.5–5.0 TPa were studied by first principles based on density functional theory, which may provide knowledge reserved for subsequent theoretical and experimental studies of metallic hydrogen and its superconductivity. The results show that the phase transition sequence of solid hydrogen under extreme high pressure is: I41/amd→oC12→cI16. For the same structure, with the increase of pressure, the electron-phonon-induced interaction decreases, the density of electronic states at the Fermi surface decreases, the vibration frequency increases, and the superconducting transition temperature changes. When the pressure is 2.0 TPa, the oC12 structure of solid hydrogen can obtain the highest Tc of 418 K (coulomb pseudopotential parameter μ*=0.10). This work provides a reference for further theoretical and experimental research on metallic hydrogen and its superconductivity.
-
It is very important for mining and civil construction to predict the morphology distribution of cracks induced by blasting. Hence, many researchers have paid their attention to dynamic fracture behavior of rocks due to drilling and blasting operations[1-3]. A number of experiments and numerical simulations have been conducted to investigate the blasting-induced fractures in the near borehole zone as well as in the far field[4-5]. In order to gain high fidelity in simulating the complex responses of rocks subjected to blasting loading, a realistic constitutive model is required. In the last 20 years, various macro-scale material models have been proposed, from relatively simple ones to more sophisticated, and their capabilities in describing actual nonlinear behavior of material under different loading conditions have been evaluated[6].
During blasting operation, chemical reactions of explosive in borehole occur rapidly, and instantaneously a shock/stress wave applies to borehole wall. Initially, a crushed zone around the borehole is developed by the shock/stress wave. Then, a radial shock/stress wave propagates away from borehole, and its magnitude decreases. Once the radial shock/stress drops below the local dynamic compressive strength, no shear damage occurs. At the same time, a tensile tangential stress with enough strength can be developed behind the radial compressive stress wave, which results in an extension of the existing flaws or a creation of new radial cracks. If there is a nearby free boundary, the incident compressive stress wave changes to a tensile stress wave upon reflection, and reflects back into the rock. In this case, if the dynamic tensile strength of rock is exceeded, spall cracks appear close to the free boundary.
The purpose of this paper is to conduct a numerical study on borehole blasting-induced fractures in rocks. First, a dynamic constitutive model for rocks based on the previous work of concrete[7] is briefly described, and the values of various parameters in the model for granite are estimated. The model is then employed to simulate the borehole blasting-induced fractures in granitic rocks. Comparisons between the numerical results and the experimental observations are made, and a discussion is given.
1. Dynamic Constitutive Model for Rocks
A number of models for concrete-like materials, such as TCK model[8], HJC model[9], RHT model[10], K&C model[11], have been developed. The sophisticated numerical models are increasingly used as they are capable of describing the material behavior under high strain rate loading. However, these models have been found to have some serious flaws, and cannot predict the experimentally observed crack patterns or exhibit improper behavior under certain loading conditions[7, 12-15].
In the following, a dynamic constitutive model for rocks is briefly described according to equation of state (EOS) and strength model, based on the previous work on concrete[7].
1.1 EOS
A typical form of EOS is the so-called p-
α relation, which is proved to be capable of representing brittle material’s response behavior at high pressures, and it allows for a reasonably detailed description of the compaction behavior at low pressure ranges as well, as shown schematically in Fig.1.pcrush corresponds to the pore collapse pressure beyond which plastic compaction occurs, andplock is the pressure when porosityα reaches 1,ftt is the tensile strength,ρ0 is the initial density,ρs0 refers to the density of the initial solid.Figure 1. Schematic diagram of EOS[7]The EOS for compression (p≥0) is given by
p=K1ˉμ+K2ˉμ2+K3ˉμ3 (1) where p denotes pressure, K1, K2, K3 are constants, and
ˉμ is defined byˉμ=ραρ0α0−1=αα0(1+μ)−1 (2) where
ρ is the current density,μ=ρ/ρ0−1 specifies the volumetric strain,α0 =ρ s0/ρ 0 andα=ρs/ρ represent the initial porosity and the current porosity, respectively.ρs refers to the density of fully compacted solid. Physically,α is a function of the hydro-static pressure p, and is expressed asα=1+(α0−1)(plock−pplock−pcrush)n (3) where n is the compaction exponent.
When material withstands hydro-static tension, the EOS for tension (p<0) is given by
p=K1ˉμ (4) α=α0(1+p/K1)/(1+μ) (5) 1.2 Strength Model
The strength model takes into account various effects, such as pressure hardening, damage softening, third stress invariant (Lode angle) and strain rate. The strength surface Y, shown schematically in Fig.2, can be written as[7]
Y={3(p+ftt)R(θ,e) p<0[3ftt + 3p(fcc−3ftt)/fcc]R(θ,e) 0⩽p⩽fcc/3{fcc+Bfc′[p/fc′−fcc/(3fc′)]N}R(θ,e) p>fcc/3 (6) where p is the hydro-static pressure, parameters B and N are constants, R(θ,e) is a function of the Lode angle θ and the tensile-to-compressive meridian ratio e,
fc′ is the static uni-axial compressive strength, the compressive strengthfcc and the tensile strengthftt are defined byfcc=fc′Dm_tηc (7) ftt=ftDtηt (8) where
ft is the static uni-axial tensile strength.Dm_t is the compression dynamic increase factor due to strain rate effect only, and can be expressed as[7, 16]Dm_t=(Dt−1)ft/fc′+1 (9) where
Dt is the tension dynamic increase factor determined byDt={tanh[(lg˙ε˙ε0−Wx)S][FmWy−1]+1}Wy (10) where
Fm ,Wx ,Wy andS are experimental constants,˙ε is the strain rate, and˙ε0 is the reference strain rate, usually taken˙ε0=1.0 s−1 .ηc is the damage function for compression, which can be expressed asηc={l+(1−l)η(λ)λ⩽λmr+(1−r)η(λ)λ>λm (11) where l and r are constants[7],
λm is the value of shear damage (λ ) when strength reaches its maximum value under compression.η(λ) is defined asη(λ)=aλ(λ−1)exp(−bλ) (12) in which a and b can be determined by setting
η(λ) = 1 and∂η∂λ=0 whenλ =λ m.ηt is the damage function for tension which can be written asηt=[1+(c1εtεfrac)3]exp(−c2εtεfrac)−εtεfrac(1+c31)exp(−c2) (13) where c1 and c2 are constants[7],
εt denotes the tensile strain andεfrac is the fracture strain.The residual strength (
rfc′ ) surface for rocks, shown schematically in Fig.2, can be obtained from Eq.(6) by settingftt=0 andfcc=rfc′ , vizY={3pR(θ,e) 0<p⩽rfc′/3{rfc′+Bfc′[p/fc′−rfc′/(3fc′)]N}R(θ,e) p>rfc′/3 (14) 2. Numerical Simulations
Granite is selected for investigating the dynamic fractures which result from borehole blast loading.
2.1 Evaluation of Various Parameters in the Model
Table 1 lists the values of the various parameters used in the dynamic constitutive model for granite. As to how to determine the values of the various parameters in the model, more details are presented in [7, 15-17].
p-α relation ρ0/(kg∙m−3) pcrush/MPa plock/GPa n K1/GPa K2/TPa K3/TPa 2660 50.5 3 3 25.7 −3 150 Strength surface Strain rate effect fc′/MPa ft/MPa B N G/GPa Fm Wx 161.5 7.3 2.59 0.66 21.9 10 1.6 Strain rate effect Shear damage Wy S ˙ε0/s−1 λs λm l r 5.5 0.8 1.0 4.6 0.3 0.45 0.3 Lode effect Tensile damage e1 e2 e3 c1 c2 εfrac 0.65 0.01 5 3 6.93 0.007 Fig.3 shows the comparison of the strength surface between Eq.(6) (with B=2.59, N=0.66) and the triaxial test data for granite[17]. It can be seen from Fig.3 that a good agreement is obtained. Similarly, Fig.4 shows the tensile strengths/dynamic increase factor obtained by Eq.(10) and the test results of various rocks at different strain rates[18-23]. It is clear from Fig.4 that a good agreement is achieved.
Figure 3. Comparison of the strength surface between Eq.(6) (with B=2.59, N=0.66) and the triaxial test data for granite[17]2.2 Numerical Results
In the following, numerical simulations are carried out for the response of the granite targets subjected to borehole blasting loading. The dynamic fracture behavior of two kinds of granite samples are studied, namely, cylindrical sample as reported in the literature and square sample as examined in our own laboratory.
2.2.1 Cylindrical Rock Sample
In consideration of the sizes of the cylindrical granite samples prepared for laboratory-scale blasting experiments by Dehghan Banadaki and Mohanty[17] (with a diameter of 144 mm, a height of 150 mm and a borehole diameter of 6.45 mm), a circular plane strain model with an outer diameter of 144 mm is made in our simulation, as shown in Fig.5, being a scaled close-up view of the borehole region. Multi-material Euler solver is used for modeling PETN explosive, polyethylene and air. Lagrangian descriptions are used for modeling the copper tube and granite.
The material model and the properties of PETN explosive, polyethylene, air and copper tube used in the simulation are given in Ref.[17]. The values of various parameters in the constitutive model for granite are listed in Table 1.
Fig.6 shows the comparison of the peak pressures between our simulation results of the present model, the numerical results[17], and the experimental results by Dehghan Banadaki and Mohanty[24]. It can be seen from Fig.6 that a good agreement is obtained.
In order to characterize the damping behavior of stress in granite, the peak pressure p in granite is expressed in an exponential form as
pp0=(dd0)−γ (15) where p0 is the peak pressure on the borehole wall, d0 is the initial radius of the borehole, d is the distance from the center point of the borehole,
γ is an index. It is evident from Fig.6 that Eq.(15) withγ =1.6 correlates well with the experimental results.Fig.7 shows the comparison between the crack patterns predicted numerically based on the present model and the one observed experimentally in the cylindrical granite sample[17]. It is clear that a relatively good agreement on the crack pattern is obtained. It is also clear that the stress waves produce three distinct crack regions in the cylindrical granite sample: densely populated smaller cracks around the borehole, a few large radial cracks propagating towards the outer boundary, and circumferential cracks close to the sample boundary which are due to the reflected tension stress.
Figure 7. Comparison of the crack patterns between the numerical prediction and the experiment of the cylindrical granite sample[17]In order to make an assessment of the contributions of both the compression/shear stress and the tensile stress to the crack patterns, Fig.8 shows the numerically predicted crack pattern which results from the tension stress only. Fom Fig.8 and Fig.7(a), it is apparent that there are virtually few changes in crack patterns, both having the large radial and circumferential cracks caused by the same tensile stress, however, the smaller cracks around the borehole in Fig.8 are much less than those in Fig.7(a). In another word, crack patterns are mainly caused by tensile stress, and smaller cracks around borehole are created largely by compression/shear stress.
2.2.2 Square Granite Sample
The laboratory-scale single-hole blasting tests are also carried out in order to validate further the accuracy and the reliability of the present model. Two square granite samples with a side length of 400 mm and a height of 100 mm are employed in the experiments. The borehole diameter is 4 mm and a series of concentric rings is drawn on the top surface of the samples (see Fig.9), so that the damage regions induced by detonation can be assessed visually.
A cylindrical RDX explosive enclosed by an aluminum sheath (see Fig.10) is tightly installed in the borehole of the No.1 sample, while an unwrapped RDX explosive is inserted into the borehole of the No.2 sample. The density of the RDX is 1700 kg/m3, and the material model and properties of the RDX explosive and the aluminum sheath used in the experiments are given in ref.[25]. The values of the various parameters in the constitutive model for granite are listed in Table 1.
Fig.11 and Fig.12 show the comparisons of the crack patterns between the numerical predictions from the present model and the ones observed experimentally in the square granite samples. It can be seen from Fig.11 and Fig.12 that good agreements are obtained. It should be mentioned here that No.1 sample receives less damage due to less RDX explosive used in the test, and that No.2 sample is broken up into four major pieces due to more RDX explosive employed in the experiment. Severe damages and small cracks are induced in the vicinity of the boreholes of both samples, as can be seen clearly from Fig.11(b) and Fig.12(b).
3. Conclusions
A numerical study on the borehole blasting-induced fractures in rocks is conducted in this paper, using a dynamic constitutive model developed previously for concrete. Two kinds of granite rocks are simulated numerically, one in the cylindrical form and the other in the square form. The numerical results are compared with the corresponding experiments. Main conclusions can be drawn as follows.
(1) The crack patterns predicted numerically from the present model are found to be in good agreement with the experimental observations, both in cylindrical and square granite samples subjected to borehole blasting loading.
(2) The peak pressures predicted numerically based on the present model are found to be in good agreement with the test data.
(3) Crack pattern observed experimentally in the rock sample is mainly caused by the tensile stress, while the smaller cracks in the vicinity of the borehole are created largely by compression/shear stress.
(4) The consistency between the numerical results and the experimental observations demonstrates the accuracy and reliability of the present model. Thus the model can be used in the numerical simulations of the response and the failure of rocks under blasting loading.
-
图 1 (a) 极端高压下各结构相对于I41/amd相固体氢的焓差随压强的变化,(b) 考虑零点振动能后极端高压下各结构相对于oC12相固体氢的焓差随压强的变化
Figure 1. (a) Calculated enthalpies of various structures as a function of pressure with respect to I41/amd structure; (b) calculated enthalpies of various structures with the inclusion of ZPE as a function of pressure with respect to the oC12 structure
图 2 3.0 TPa下I41/amd相(a)、oC12相(b)(不同颜色的原子代表不同层)、cI16相(c)固体氢的结构和oC12相电子局域函数(d)(均使用VESTA绘制[28])
Figure 2. I41/amd structure (a), oC12 structure (b) (spheres of different colors indicate different layers), and cI16 structure (c), and ELF profile of the oC12 structure (d) at 3.0 TPa (All structures and the ELF profile were drawn by using VESTA[28])
表 1 优化后的结构参数
Table 1. Optimized structural parameters
Structure p/TPa Lattice parameters Atomic coordinates (fractional) I41/amd 1.0 a=b=1.067 9 Å, c=2.993 9 Å,
α=β=γ=90°H1(4a) 0 1 0 oC12(Cmcm) 3.0 a=0.947 9 Å, b=2.812 3 Å, c=2.319 0 Å,
α=β=γ=90°H1(8f) 0.500 0 0.140 5 0.580 7 H2(4c) 0.500 0 0.595 6 1.250 0 cI16(I¯43d) 4.0 a=b=c=1.929 3 Å, α=β=γ=90° H1(16c) 0.955 5 0.455 5 0.044 5 表 2 不同压强下各结构的超导转变温度及超导计算中的重要参数
Table 2. Calculated superconducting parameters and the superconducting transition temperature for I41/amd, oC12, cI16 under different pressures
Structure p/TPa λ ωlog/K NEF/
(states·Ry−1/atom)Tc/K μ*=0.10 μ*=0.13 μ*=0.17 I41/amd 0.5 1.94 1 784.5 0.23 339 307 280 1.0 1.87 2 039.0 0.20 382 354 320 1.5 1.53 2 539.5 0.19 385 348 310 oC12 2.0 2.65 1 580.2 0.18 418 391 356 3.0 1.77 2 324.7 0.16 398 356 318 4.0 1.48 2 822.3 0.14 388 353 312 cI16 3.0 1.73 2 303.0 0.16 356 345 307 4.0 1.42 2 782.6 0.14 392 362 322 5.0 1.27 3 159.0 0.13 295 258 217 -
[1] WIGNER E, HUNTINGTON H B. On the possibility of a metallic modification of hydrogen [J]. The Journal of Chemical Physics, 1935, 3(12): 764–770. doi: 10.1063/1.1749590 [2] PICKARD C J, NEEDS R J. Structure of phase Ⅲ of solid hydrogen [J]. Nature Physics, 2007, 3(7): 473–476. doi: 10.1038/nphys625 [3] MCMAHON J M, CEPERLEY D M. Ground-state structures of atomic metallic hydrogen [J]. Physical Review Letters, 2011, 106(16): 165302. doi: 10.1103/PhysRevLett.106.165302 [4] EREMETS M I, DROZDOV A P, KONG P P, et al. Semimetallic molecular hydrogen at pressure above 350 GPa [J]. Nature Physics, 2019, 15(12): 1246–1249. doi: 10.1038/s41567-019-0646-x [5] DIAS R P, SILVERA I F. Observation of the Wigner-Huntington transition to metallic hydrogen [J]. Science, 2017, 355(6326): 715–718. doi: 10.1126/science.aal1579 [6] LOUBEYRE P, OCCELLI F, DUMAS P. Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen [J]. Nature, 2020, 577(7792): 631–635. doi: 10.1038/s41586-019-1927-3 [7] ASHCROFT N W. Metallic hydrogen: a high-temperature superconductor? [J]. Physical Review Letters, 1968, 21(26): 1748–1749. doi: 10.1103/PhysRevLett.21.1748 [8] CUDAZZO P, PROFETA G, SANNA A, et al. Ab initio description of high-temperature superconductivity in dense molecular hydrogen [J]. Physical Review Letters, 2008, 100(25): 257001. doi: 10.1103/PhysRevLett.100.257001 [9] YAN Y, GONG J, LIU Y H. Ab initio studies of superconductivity in monatomic metallic hydrogen under high pressure [J]. Physics Letters A, 2011, 375(9): 1264–1268. doi: 10.1016/j.physleta.2011.01.045 [10] MAKSIMOV E G, SAVRASOV D Y. Lattice stability and superconductivity of the metallic hydrogen at high pressure [J]. Solid State Communications, 2001, 119(10/11): 569–572. [11] SZCZȨS̀NIAK R, JAROSIK M W. The superconducting state in metallic hydrogen under pressure at 2 000 GPa [J]. Solid State Communications, 2009, 149(45/46): 2053–2057. [12] MCMAHON J M, CEPERLEY D M. High-temperature superconductivity in atomic metallic hydrogen [J]. Physical Review B, 2011, 84(14): 144515. doi: 10.1103/PhysRevB.84.144515 [13] LIU H Y, WANG H, MA Y M. Quasi-molecular and atomic phases of dense solid hydrogen [J]. The Journal of Physical Chemistry C, 2012, 116(16): 9221–9226. doi: 10.1021/jp301596v [14] KRESSE G, FURTHMÜLLER J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set [J]. Physical Review B, 1996, 54(16): 11169–11186. doi: 10.1103/PhysRevB.54.11169 [15] BLÖCHL P E. Projector augmented-wave method [J]. Physical Review B, 1994, 50(24): 17953–17979. doi: 10.1103/PhysRevB.50.17953 [16] PERDEW J P, BURKE K, ERNZERHOF M. Generalized gradient approximation made simple [J]. Physical Review Letters, 1996, 77(18): 3865–3868. doi: 10.1103/PhysRevLett.77.3865 [17] MONKHORST H J, PACK J D. Special points for Brillouin-zone integrations [J]. Physical Review B, 1976, 13(12): 5188–5192. doi: 10.1103/PhysRevB.13.5188 [18] TOGO A, CHAPUT L, TADANO T, et al. Implementation strategies in phonopy and phono3py [J]. Journal of Physics: Condensed Matter, 2023, 35(35): 353001. doi: 10.1088/1361-648X/acd831 [19] TOGO A. First-principles phonon calculations with phonopy and phono3py [J]. Journal of the Physical Society of Japan, 2023, 92(1): 012001. doi: 10.7566/JPSJ.92.012001 [20] GIANNOZZI P, BARONI S, BONINI N, et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials [J]. Journal of Physics: Condensed Matter, 2009, 21(39): 395502. doi: 10.1088/0953-8984/21/39/395502 [21] MIGDAL A B. Interaction between electrons and lattice vibrations in a normal metal [J]. Soviet Physics JETP, 1958, 34(7): 996–1001. [22] ÉLIASHBERG G M. Interactions between electrons and lattice vibrations in a superconductor [J]. Soviet Physics JETP, 1960, 11(3): 696–702. [23] SCALAPINO D J, SCHRIEFFER J R, WILKINS J W. Strong-coupling superconductivity.Ⅰ [J]. Physical Review, 1966, 148(1): 263–279. doi: 10.1103/PhysRev.148.263 [24] VIDBERG H J, SERENE J W. Solving the Eliashberg equations by means of N-point Padé approximants [J]. Journal of Low Temperature Physics, 1977, 29(3/4): 179–192. [25] HANFLAND M, SYASSEN K, CHRISTENSEN N E, et al. New high-pressure phases of lithium [J]. Nature, 2000, 408(6809): 174–178. doi: 10.1038/35041515 [26] GREGORYANZ E, LUNDEGAARD L F, MCMAHON M I, et al. Structural diversity of sodium [J]. Science, 2008, 320(5879): 1054–1057. doi: 10.1126/science.1155715 [27] MCMAHON M I, GREGORYANZ E, LUNDEGAARD L F, et al. Structure of sodium above 100 GPa by single-crystal X-ray diffraction [J]. Proceedings of the National Academy of Sciences of the United States of America, 2007, 104(44): 17297–17299. [28] MOMMA K, IZUMI F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data [J]. Journal of Applied Crystallography, 2011, 44(6): 1272–1276. doi: 10.1107/S0021889811038970 [29] ALLEN P B, DYNES R C. Transition temperature of strong-coupled superconductors reanalyzed [J]. Physical Review B, 1975, 12(3): 905–922. doi: 10.1103/PhysRevB.12.905 [30] DROZDOV A P, EREMETS M I, TROYAN I A, et al. Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system [J]. Nature, 2015, 525(7567): 73–76. doi: 10.1038/nature14964 [31] MA L, WANG K, XIE Y, et al. High-temperature superconducting phase in clathrate calcium hydride CaH6 up to 215 K at a pressure of 172 GPa [J]. Physical Review Letters, 2022, 128(16): 167001. doi: 10.1103/PhysRevLett.128.167001 [32] KONG P P, MINKOV V S, KUZOVNIKOV M A, et al. Superconductivity up to 243 K in the yttrium-hydrogen system under high pressure [J]. Nature Communications, 2021, 12(1): 5075. doi: 10.1038/s41467-021-25372-2 [33] DROZDOV A P, KONG P P, MINKOV V S, et al. Superconductivity at 250 K in lanthanum hydride under high pressures [J]. Nature, 2019, 569(7757): 528–531. doi: 10.1038/s41586-019-1201-8 [34] SOMAYAZULU M, AHART M, MISHRA A K, et al. Evidence for superconductivity above 260 K in lanthanum superhydride at megabar pressures [J]. Physical Review Letters, 2019, 122(2): 027001. doi: 10.1103/PhysRevLett.122.027001 -