Cite this Article
Zhang Xiu-Lu, Han Yun-Xia, Jia Hong, Qu Nuo, Liu Zhong-Li. Molecular Dynamics Simulations of the Elastic Anisotropy of Pd at Extreme Conditions. Communications in Theoretical Physics, 2018, 69(6): 735  
Permissions
Molecular Dynamics Simulations of the Elastic Anisotropy of Pd at Extreme Conditions
Zhang Xiu-Lu1, Han Yun-Xia2, Jia Hong2, Qu Nuo, Liu Zhong-Li2, †
Laboratory for Extreme Conditions Matter Properties, Southwest University of Science and Technology, 621010 Mianyang, China
College of Physics and Electric Information, Luoyang Normal University, Luoyang 471022, China

 

† Corresponding author. E-mail: zl.liu@163.com

Supported by the National Natural Science Foundation of China (41574076), the Basic Research of Technology Program of China under Grant No. JSHS2014404B002, the Young Core Teacher Scheme of Henan Province under Grant No. 2014GGJS-108, and key project of science and technology research of Henan Provincial Education Department under Grant No. 18A140024

Abstract
Abstract

It is very interesting to discover the elastic properties of engineering material palladium, especially its elastic anisotropy along Hugoniot states. We here investigate the evolution of its high pressure and temperature (PT) elastic ansotropy along Hugoniot using molecular dynamics simulations based on accurate classical interatomic potential. In order to testify the validity of the interatomic potential of Pd in describing the high PT elastic properties, we calculate its isothermal and adiabatic elastic moduli using molecular dynamics method. The obtained data are in good agreement with experimental data. From the isothermal elastic constants, we deduce the Hugoniot acoustic velocities and find that the resulting data are in good agreement with experimental acoustic velocity data. Based on the reliable elastic constants, we further investigate the spacial elastic ansotropy along Hugoniot PT states. It is found that the spacial elastic anisotropy of Pd increases along Hugoniot states.

Keyword:elastic anisotropy;Hugoniot states;extreme conditions
1 Introduction

Shock response characteristics of materials may reflect how they behave under dynamic load with increasing pressure and temperature simultaneously, specifically along a fixed shock orientation. Palladium plays an important role in variety of applications such as dentistry, jewelery, catalysis, metalizing ceramics, hydrogen storage, and so on.[1] In addition, bulk Pd is frequently used as a pressure calibrator in high-pressure experiments because of its stable face-centered-cubic (fcc) structure under high pressure. Theoretical studies of shock response of Pd is of vital importance for it as a pressure calibrator, especially in the regions where experiments cannot reach. The elastic properties of materials at high pressures are necessary for us to understand their underlying static and dynamic responses to pressure, mechanical strength, and the equation of states.

The elastic anisotropy can be conveniently achieved from elastic constants for all crystal systems. Ab initio calculations of elastic constants were often conducted at 0 K without any anharmonic effects, and the inclusion of temperature effects in ab initio elastic constants are always complex and computationally expensive. However, the elastic constants of a crystal decrease by over 20% with increasing temperature.[2–5] The temperature dependence of the elastic constants of a material is crucial for predicting and understanding the mechanical strength and even phase transitions of a crystal.[6] Furthermore, elastic anisotropy resulted from the elastic constants has also many other implications in dislocation dynamics and even geophysical significance of materials.[7] It has been demonstrated that the elastic anisotropy influences the nanoscale precursor textures in shape memory alloys and other materials.[8–9]

It is known that molecular dynamics (MD) simulations naturally includes high-temperature anharmonic effects to any order beyond Debye temperature.[10–11] And they especially apply to high-temperature free energy calculations[12] and thermal fluctuation average, for instance in the calculation of elastic constants. The combination of MD with accurate interatomic potential of crystal systems provides a valid way to yield reliable elastic constants under high PT.

The rest of the paper is organized as follows. In Sec. 2, we describe the details about the used interatomic potential of Pd and the computation method of high PT elastic constants. The results and discussion are addressed in Sec. 3. Section 4 presents a short summary of the conclusions.

2 Methodology

We calculate the high pressure and temperature (PT) elastic constants of Pd using classical MD simulations. In order to accurately describe the interatomic interactions of Pd atoms, we adopt the embedded-atom-method (EAM) potential.[13–14] Our previous study[1] verified the high quality of the EAM potential of Pd in describing the high PT properties. For a system with N atoms, the total potential energy Etot is the sum of the embedding energy F and a pair potential ϕ,

where rij is the distance from atom i to j. Fi(ρi) is the embedding energy, i.e., the energy used to embed atom i into the background electron density ρi of the rest atoms except atom j,

The high-temperature elastic constants of Pd are calculated based on direct deformation of the system by applying negative and positive strain of 5% of the lattice parameters in molecular dynamics simulation. Before applying deformations, the simulation box is fully relaxed in the NPT ensemble to equilibrate at the PT conditions, where N is the number of particles, P is the external pressure, and T is temperature. A Langevin thermostat is used in keeping the deformations of the system, with a damping parameter of 0.01 picoseconds (or 10 timesteps). This thermostat is used in combination with the microcanonical ensemble to perform Brownian dynamics. Deformations are tested to yield converged elastic constants. All the MD simulations are conducted with LAMMPS[15–16] package. The simulation box constructed from the multiplication 11 × 11 × 11 of face-centered-cubic (fcc) conventional unit cell containing 5324 atoms.

3 Results and Discussions
3.1 Isothermal Elastic Constants

The 300 K-isothermal elastic constants as a function of pressure are obtained and illustrated in Fig. 1, in comparison with experimental data at 0 GPa.[17] All the three elastic constants, C11, C12 and C44 of the face-centered-cubic Pd, increase with pressure monotonically. While, it is found that they are actually not the linear functions of pressure. We find they are most fit to the second order polynomial function of pressure,

where C(0) is the elastic constants at 0 GPa, and a and b are fitting parameters.

Fig. 1 (Color online) The 300 K isothermal elastic constants of Pd as a function of pressure. The solid symbols are our calculated data. The lines are the second polynomial fit of the elastic constants. The red symbols are experimental values at 0 GPa.[17]

After fitting the second order polynomial function to the three sets of elastic constants data respectively, we yield the function

for C11, and

for C12, and

for C44.

3.2 Isentropic Elastic Constants

The experimental elastic constants data[18–19] were isentropic data, therefore the calculated isothermal elastic constants should be converted in order to be compared with experiments.[5] For the cubic crystal, the conversion between isothermal and isentropic elastic constants are,[6]

where

where BT = V(∂2F/∂V2)T is the isothermal bulk modulus and β = (∂V/∂T)P/V is the volume thermal expansion coefficient. Such necessary data were achieved directly from our previous work.[12]

The calculated temperature dependence of the isentropic elastic constants are shown in Fig. 2. Good agreement of and with experimental data is achieved overall. However, deviates from experiment slightly, i.e., the most difference is below 35 GPa in the whole temperature range. As shown in Fig. 2, all the calculated elastic constants decrease with increasing temperature. The temperature variation of the isentropic elastic constants can be described accurately by the Varshni equation,[20]

where is the elastic constant at 0 K, s and t are fitting parameters.

Fig. 2 (a) High-temperature isentropic elastic constants of Pd at 0 GPa, in comparison with experimental data. Open symbols are experimental data from Refs. [21], respectively. (b) The calculated isothermal bulk moduli of Pd compared with experiment.[22] All the calculated moduli were fitted to the Varshni equation.[20]
3.3 Sound Velocity Along Hugoniot

The obtained elastic constants and the corresponding Hugoniot PT points are listed in Table 1 and Fig. 3(a), respectively.

Fig. 3 (a) The Hugoniot pressure-temperature relations. (b) The Hugoniot acoustic velocities of Pd compared with the experimental data.[23]
Table 1

The elastic constants along Hugoniot PT states.

.

The phase velocity v and polarization of the three waves along a fixed propagation direction defined by the unit vector ni are given by Cristoffel equation,

where Cijkl is the fourth-rank tensor description of the elastic constants, n is the propagation direction, u the polarization vector, and ρ the density. From the Cristoffel equation, the isotropically averaged aggregate velocities are then derived as,

where vP, vB, and vS are the compressional, bulk, and shear acoustic velocities, respectively. In Eqs. (13) and (15), the average isotropic shear modulus G of polycrystalline according to Voigt-Reuss-Hill approximations from single crystal elastic constants. The density ρ in Eqs. (13), (14), and (15) are derived by ρ = M/V, where M is the molecular mass and V is the corresponding unit cell volume at specific P and T, illustrated in the inset of Fig. 3(a). The equilibrium volume V0 of Pd at 0 GPa and 300 K is 14.86 A3.

The average wave velocity vm is calculated from,

From elastic constant data we can calculate the Debye temperature ΘD from the average acoustic velocity vm via the following equation:

where h is Planck’s constant, k the Boltamann’s constant, NA the Avogadro’s number, n the number of atoms in the unit cell.

The acoustic velocities of Pd as a function of pressure along Hugoniot states are illustrated in Fig. 3(b) in comparison with experimental data. The calculated compressional acoustic velocities are in very good agreement with experiment,[23] suggesting the elastic constants along Hugoniot states are accurate for further elastic anisotropy analyses.

As shown in Fig. 4, the Debye temperature of Pd along the Hugoniot states Debye temperature first has an abrupt jump and then has a slight decrease beyond ∼1000 K. The volume variation of Debye temperature are shown in the inset of Fig. 4. Along Hugoniot states, Debye temperature is also not the linear functiong of volume.

Fig. 4 Debye temperature of Pd at 0 GPa.
3.4 Elastic Anisotropy Along Hugoniot

From Subsecs. 3.1 to 3.3, we see that the calculated elastic constants reflect the high PT elastic behaviors of Pd very well. The mechanical anisotropy has an important implication in crystalline physics and engineering applications,[24] especially in understanding the dynamics response properties of Pd as a typical pressure calibrator. We here investigate the spacial anisotropy of Young’s modulus, shear modulus, and Poisson’s ratio to understand the mechanical anisotropy. Then the spacial orientation variation of the mechanical properties of Pd with PT are analyzed along the Hugoniot states.

The Young’s modulus, shear modulus, and Poisson’s ratio are all dependent on orientation.[24–25] It is specially interesting to uncover the evolution process of these elastic properties along the Hugoniot states. From the elastic constants, we calculate their orientation dependence for Pd along the Hugoniot states. Detailed method can be found in Ref. [25]. In Fig. 5 we plot the evolutional process of their 3D orientation dependence. In the [111] direction and the counterpart directions the maximum of Young’s modulus appears, but the anisotropy of Young’s modulus apparently increases with PT along Hugoniot states. The range of maximum in Young’s modulus becomes slim with increasing PT. Such increase of anisotropy can also be seen for the shear modulus in the second column of Fig. 5, mainly implying in the orientational dependence of its minimum along Hugoniot states (green color). Negative Poisson’s ratio gradually appears and its absolute value increases with increasing PT, also suggesting the increase of the elastic anisotropy of Pd along Hugoniot states.

Fig. 5 (Color online) The spacial orientation variation of the Young’s modulus, shear modulus, and poisson ratio of Pd along the principle Hugoniot. In the multi-colored figures, blue color accounts for maximum, green for minimum positive, red for negative.

In order to understand the variation of elastic anisotropy of Pd along specific directions, the 2D Poisson’s ratio in the (110) plane is calculated and plotted in Fig. 6. From Fig. 6, we see that the maximum, minimum positive, and minimum negative values of Poisson’s ratios are distributed along [110] and [001] directions inhomogeneously. It is noteworthy that along some directions, such as [110], the positive and minimum negative values of Poisson’s ratios occur simultaneously. The simultaneous occurrence of positive and minimum negative Poisson’s ratios in the same direction indicates that the Poisson’s has the possibility of being negative. To understand the negative Poisson’s ratio, we need to recall the meaning of Poisson’s ratio. It is known that Poisson’s ratio is a measure of the extent to which the material contract in the directions perpendicular to the direction of stretching. The minimum negative value can be understood as the expansion in the [001] direction, which perpendicular to the [110] direction of stretching, while the maximum positive value corresponds to the normal contraction in [001] when stretching in [110] direction.

Fig. 6 (Color online) 2D representation of Poissons ratio in the (110) plane for Pd at (a) 35 GPa and 466 K, (b) 95 GPa and 1367 K and (c) 167 GPa and 3477 K; maximum (blue), minimum positive (green) and minimum negative (red).
4 Conclusion

In conclusion, we have investigated the elastic properties of Pd under high PT systematically. Firstly, we calculate the isothermal elastic constants of Pd from molecular dyanmics simulations based on accurate interatomic potential. From the isothermal elastic constants, we deduce the high temperature isentropic elastic moduli of Pd and find they are in good agreement with experimental data at 0 GPa. Secondly, we achieve the acoustic velocities of Pd at high PT along Hugoniot states. The derived compressional wave velocity is also in good agreement with shock wave experiment. Finally, we investigate the elastic anisotropy of Pd along Hugoniot states. It is found that the spacial elastic anisotropy of Pd, reflected by the Young’s modulus, the shear modulus, the acoustic velocity and Poisson ratio, increases with increasing PT along Hugoniot states.

Reference
[1] Liu Z. L. Yang J. H. Cai L. C. et al Phys. Rev. B 83 2011 144113
[2] Featherston F. H. Neighbours J. R. Phys. Rev. 130 1963 1324
[3] Gerlich D. Fisher E. S. J. Phys. Chem. Solids 30 1969 1197
[4] Walker E. Bujard P. Solid State Commun. 34 1980 691
[5] Wang Y. Wang J. J. Zhang H. et al. J. Phys.: Condens. Matter 22 2010 225404
[6] Gülseren O. Cohen R. E. Phys. Rev. B 65 2002 064103
[7] Ledbetter H. Migliori A. J. Appl. Phys. 100 2006 063516
[8] Lloveras P. Castan T. Porta M. et al. Phys. Rev. Lett. 100 2008 165707
[9] Legut D. Friak M. Sob M. Phys. Rev. Lett. 99 2007 016402
[10] Errea I. Calandra M. Mauri F. Phys. Rev. B 89 2014 064302
[11] Car R. Parrinello M. Phys. Rev. Lett. 55 1985 2471
[12] Liu Z. L. Li R. Zhang X. L. et al. Comput. Phys. Commun. 213 2017 122
[13] Daw M. S. Baskes M. I. Phys. Rev. B 29 1984 6443
[14] Foiles S. M. Baskes M. I. Daw M. S. Phys. Rev. B 33 1986 7983
[15] Plimpton S. J. J. Comp. Phys. 117 1995 1
[16] http://lammps.sandia.gov.
[17] Rayne J. A. Phys. Rev. 118 1960 1545
[18] Overton W. C. Gaffney J. Phys. Rev. 98 1955 969
[19] Chang Y. A. Himmel L. J. Appl. Phys. 37 1966 3567
[20] Varshni Y. P. Phys. Rev. B 2 1970 3952
[21] Yoshihara M. Mclellan R. B. Brotzen F. R. Acta Metall. 35 1987 775
[22] Touloukian Y. S. Kirby R. K. Taylor R. E. Lee T. Y. R. Thermophysical Properties of Matter Plenum New York 1970
[23] Duffy T. S. Ahrens T. J. J. Geophys. Res. 97 1992 4503
[24] Liu Z. L. Tao Y. P. Zhang X. L. Cai L. C. Comput. Mater. Sci. 114 2016 72
[25] Marmier A. Lethbridge Z. A. D. Walton R. I. et al. Comput. Phys. Commun. 181 2010 2102