Molecular dynamics modeling of the Hugoniot states of aluminum

In this study, molecular dynamics (MD) simulations coupled with multi-scale shock technique (MSST) are used to predict the Hugoniot curve P H , Gr¨uneisen coefﬁcient γ and melting temperature T m of single crystal (SC) and nanocrystalline (NC) aluminum (Al) with grain sizes of 6 and 60 nm at dynamic high pressure. The linear relation between the shock wave velocity and particle velocity is reproduced, and the results indicate that there is nearly no difference for the Hugoniot of SC and NC Al, which could be explained by the fact that the grain size effect on P H can be negligible at high pressure. Some empirical models are used to predict γ and T m , which exhibit an opposite behavior. In addition, it is found that the melting pressure and temperature are 107.5 GPa, 3063 K for SC Al, while they are 109.5 GPa, 3082 K for NC Al, which have a reasonable agreement with the published work. © 2018 Author(s). All article content, except where otherwise noted, is licensed a Creative Commons


I. INTRODUCTION
In high-pressure physics, knowledge of the shock Hugoniot, Grüneisen coefficient (γ) and melting temperature (T m ) is essential to investigate the equation of state (EOS) of metallic materials. 1-3 The Hugoniot relation has been an increasing concern, because it could be used to characterize the Grüneisen EOS and shock-induced melting of metals. Comprehensive theoretical, 4,5 experimental, 6,7 and numerical [8][9][10] investigations have been performed on the Hugoniot relation, Grüneisen EOS, and T m . As a metal with low melting point, aluminum (Al) is one of the most important and widely studied materials under extreme pressure. Mitchell et al. 11 measured the Hugoniot curve of Al with a twostage-light gas gun at the pressure of 30 -430 GPa. Wang et al. 12 investigated the Hugoniot curve of Al using the classical mean-field approach when the pressure was beyond 1000 GPa. Choudhuri and Gupta 13,14 used the laser interferometry to measure the Hugoniot relation for 1050 Al under 70 GPa. They examined the crystalline anisotropy effect of single crystals, and found that the influence of crystal orientation on the Hugoniot behavior was significant. Ju et al. 15 also investigated the shock melting of single crystal Al by multi-scale shock technique (MSST) with the help of LAMMPS, and in turn the Hugoniot relation and the Hugoniot melting curve of single crystal Al was obtained availably.
In the Grüneisen EOS, γ is an important parameter, because it could help understand the melting essence of solid at high pressure. Barton et al. 16 calculated γ at high pressure through molecular dynamics (MD), and reported that the volume dependence of γ(V ) is highly sensitive to the EOS of a Correspondence should be addressed to: xiangguozeng@scu.edu.cn, wangfang cq1978@163.com a material. Various models regarding to the volume dependence of γ have been presented. 17,18 For example, Burakovsky et al. 19 investigated the volume dependence of γ of Al at low-pressure using a concrete equation, which was based on the assumption that γ is an analytic function of specific volume. Jacobs et al. 20 also presented another expression to describe the volume dependence of γ. Nie et al. 21 compared these two functions, and found the latter was more appropriate to Al. Furthermore, Peng et al. 22 compared three expressions that described the pressure or volume dependence of γ, which were reported by Nie et al., 17 Fang et al., 23 and Cui et al., 24 respectively. The results showed that the method presented by the former was better at high pressure compared with two other ones. T m related to γ 0 is of importance to investigate the nature of material melting under intensive shock loading. Hänström et al. 25 determined the melting curve of Al up to 50 GPa using Lindemann and Simon melting equations, which showed a reasonable agreement reported by previous study. 26 Wang et al. 27 proposed a model to reveal the melting of Al at extreme pressure, and found that it could better reproduce the high pressure melting process compared to the previous empirical functions. Vočadlo et al. 28 determined the melting curve of Al up to 150 GPa using ab initio MD, which exhibited an agreement with shock data. 25 Zou et al. 29 studied T m of Al under a wide range of pressure based on the Lindemann's law and the Debye model, and the results showed a good agreement with experimental data. 25 Nie et al. 21 also studied T m of Al on the basis of γ presented by Burakovsky et al. 19 and Jacobs et al., 20 respectively, and concluded that these two equations of T m were almost the same melting curve. Shao et al. 13 investigated the melting behavior in single crystal Al under the shock loading and unloading through MD, and found that the melting pressure was about 125 -155 GPa under the shock loading and 72 GPa for unloading.
In fact, P H (V ), γ(V ) and T m (V ) had been successfully investigated by many numerical simulations, including a non-equilibrium molecular dynamics (NEMD) method, 16 an equilibrium molecular dynamics (EMD) method 30 and the first principle. 28 However, the time scale and model size required by NEMD are limited, because more computational time is needed when the model size increases. In addition, it becomes impossible that the first principle method with poor model size is used to simulate the shock propagation. A method coupled with MSST proposed by Reed et al. 31 attempted to solve the above problems. Compared with the NEMD, EMD and the first principle method, MSST can efficiently save the calculation time with a small model of MD system, and guarantee the simulation accuracy. Because the calculation time is approximate linear to the size of MD model, the MSST is employed to simulate shock compression in condensed matters in the present study.
Regarding to the relation between P H and EOS at high pressure, there are two different issues involved. One is to determine P H based on the EOS, and the other is to calculate the EOS on the basis of P H . It is widely accepted that the shock Hugoniot data is used to calculate γ and T m . The detail investigation path of the present paper is firstly utilize MSST to obtain a linear equation between the shock wave velocity (u s ) and the particle velocity (u p ) for single crystal (SC) and nanocrystalline (NC) Al. The Hugoniot curve is determined based on the Hugoniot relation. In addition, γ and T m are calculated to obtain the melting characteristic of Al at high pressure.

A. The Hugoniot curve (P H )
The Hugoniot relation is governed by three conversation equations: mass, momentum and energy, which are expressed by ρ u s − u p = ρ 0 u s (1a) where ρ, u p , u s , P, and E are the density, the particle velocity, the shock wave velocity, the pressure, and the internal energy per units mass, respectively. V = 1/ρ represents the specific volume of the shocked material. The subscript "0" denotes that these qualities are in the initial state without shock. P H is characterized by the pressure in the Hugoniot state.
Under a wide range of pressure, the relation between u s and u p is linear where C 0 is the volume speed at zero pressure, and λ is a material constant. These quantities are obtained through fitting u s and u p . According to Eqs. (1a), (1b), and (2), the Hugoniot curve is obtained when P 0 = 0 (3)

B. Grüneisen coefficient (γ)
The coefficient γ is expressed by 5 where P c (V ) is the cold pressure. t is used to represent various models for characterizing γ. When t = 0, Eq. (4) represents the Slate model γ s (V ). When t = 1, it denotes the Dugdale-MacDonald model γ DM (V ). When t = 2, Eq. (4) it becomes the Free-volume model γ f (V ). The isentrope P s (V ) is approximately satisfied with Eq. (4) at the point with V = V 0 and P = 0, and consequently Eq. (4) can be written When t = 0, 1 and 2, we have It is assumed that the first-and second-order derivatives of P s (V ) and P H (V ) are identical when V = V 0 and P = 0. While P H (V 0 ) and P H (V 0 ) can be directly derived from P H given from experiments and MD simulations, the following equation is obtained when Eq. (2) is known The substitution of Eq. (7) into Eq. (6) produces The dependence of γ on volume is commonly expressed by some simpler empirical formulas, which are as follows: where n and A are material constants, and they approximately equal to 1.18, and 1.8, respectively. In previous work, 20 the equation of γ for Al was expressed as Accordingly, Eqs. (9) and (10) can be used to describe the dependence of γ on volume for Al. However, it needs to choose a more appropriate one under the given γ 0 .

C. Melting temperature (T m )
T m is an important material parameter under extreme shock compression. The famous Lindemann melting criterion is written as According to Eqs. (9) and (10), Eq. (11) is also deduced to where T m0 is the melting temperature at zero pressure, and it is 933.5 K for Al. 19

III. MD COMPUTATIONAL DETAILS
It is widely accepted that the embedded-atom-method (EAM) potential presented by Winey et al. 32 can be used to describe the interactions between atoms in Al. Such a potential has also been proved to be capable of providing an accurate analysis of the Hugoniot curve of Al at extreme pressure by Liao et al. 33 Thus, this potential will be used in our subsequent simulations. Three models are shaped to investigate the differences of P H, γ, and T m of SC and NC Al under intensive shock. Figure 1(a) represents the configuration of SC Al, and its size is about 7.28 × 7.28 × 7.28 nm 3 with 23,328 atoms. Figure 1 Before applying the shock loading, three models are firstly equilibrated at 300 K for 10 ps in NPT ensemble so that they reach a state that neglects the residual stress in the system. The periodical boundary conditions are applied along three directions. After that, the model is applied by a shock wave of 9 -14 km/s through MSST along the x-axis [1 0 0] direction. It needs to be mentioned that the computation time should be enough to make sure that the shock wave reaches a steady state after the shock. In other words, the particle velocity, pressure, temperature, and other physical quantities are in a steady state. Consequently, the statistic results are obtained from these steady physical quantities.

IV. RESULTS AND DISCUSSION
According to previous work, 15 u s ranged from 9 to 14 km/s. In order to assess the validity of model, u s is set in the range between 9 and 14 km/s in the present study. As shown in Figure 2, the relation of u p and u s for SC Al is obtained through MSST. The square represents this work and the red line shows a good linear relation by fitting these squares. The red fitting line is given by u s = 5.280 + 1.295×u p . The blue line is characterized by u s = 5.533 + 1.213×u p . The errors of C 0 and λ are only 4.79 % and 6.33 %, respectively, which proves the accuracy.
On the other hand, the experimental data was compared with this work. The magenta line, orange line, and cyan line were determined from experiments by Chijioke et al., 34 Choudhuri et al., 14 and Marsh, 35 respectively. The magenta line was given by u s = 5.386 + 1.339×u p at up to 300 GPa. The orange line was given by u s = 5.361 + 1.290×u p at 40 -70 GPa, and the cyan line was given by u s = 5.350 + 1.340×u p at 10 -50 GPa. It is found that the orange line is better close to the red line compared to the magenta line and the cyan line. It is also found that the red line deviates from the magenta line, the orange line and the cyan line as u p increases, because there is a bigger slope λ compared with the red line.
To investigate the influence of GS on the relation between u p and u s for NC Al, u p is also calculated when u s = 9 -14 km/s, which is shown in Figure 3. It is evident that u s is very close to each other in different GSs, which suggests that the GS has an insignificant effect on the linear relation. These discrete data is fitted to be a line given by u s = 5.275 + 1.310×u p for u s = 9 -14 km/s. Compared  with SC Al, it is found that C 0 and λ are close with an error of less than 2 %. This suggests that there is almost no difference between SC and NC Al, which could be explained by the fact that the GS effect could be negligible at the extreme pressure. The comparative details of this work with the previous studies are listed in Table I. The maximum relative error of C 0 is 4.79 % and 4.89 % for SC and NC Al, respectively, and correspondingly the maximum relative error of λ is 6.33 % and 7.40 %, respectively, demonstrating the rationality of this work. Figure 4 shows the Hugoniot curve determined by Eq.
(2) is known. The square, circle and triangle represent the results from Mitchell et al., 36 Nellis et al., 37 and Wang et al., 38 respectively. The magenta curve represents the result from Ju et al., 15  curve represent this work for SC and NC Al, respectively. When V /V 0 is larger than 0.65, the curves calculated by the MSST are well consistent with the results of Mitchell et al. 36 and Wang et al. 38 When V /V 0 is less than 0.65, the curves show an insignificant difference with discrete datum from the three researchers. In fact, all curves are slightly over-estimated comparing with result of Nellis et al. 37 Due to the negligible error of C 0 and λ for SC and NC Al, it shows that there has no noticeable difference between these two curves. In addition, it is found that the curves are nearly overlapped with P H (V ) by Ju et al., 15 suggesting that this work is in accordance with other's results. Thus, it is concluded that the EAM potential is appropriate for describing the dynamic behavior, such as γ and T m of Al at high pressure. While λ of SC and NC Al is equal approximately, λ in NC Al is employed to calculate γ 0 of the Slate, Dugdale-Macdonald and Free-volume models through Eqs. (8a)-(8c), and then substitute γ 0 into Eqs. (9) and (10) to obtain γ. To validate the accuracy of this work, experimental data, 35 and MD data 15 are used to give a comparison between Eq. (9) and Eq. (10) for NC Al, as seen in Figure 5. The empirical Eqs. (9) and (10) provide a way to predict γ(V ) at the high pressure. It is shown that γ decreases quickly when V 0 /V is lower. As V 0 /V increases, the decreasing ratio degrades gradually and the difference of γ becomes smaller for the same color line. In addition, the departure point of γ 0 is as follows: γ f (V 0 ) < γ DM (V 0 ) < γ S (V 0 ) when V 0 /V =1. Since λ (1.310) in this work is smaller than that (1.340) by Choudhuri et al., 35 the solid line is above the dashed line in Figure 5. A similar case is found in Ju et al. 15 and this work.

and the blue curve and cyan
As shown in Figure 5, it is worthy noting that the black line is under other lines. In other words, γ(V ) determined by Eq. 9(a) exhibits a more apparent decreasing trend with the increase of V 0 /V, compared with those determined by Eqs. (9b)-(9d) and (10). This suggests that γ characterized by black line has the maximum change rate, which is dependent upon V 0 /V. In fact, Eq. (10) (cyan line) has been proved to be more suitable for describing γ of Al. 21 Contrarily, Eq. (9a) (black line) is inadequate for describing γ of Al at high pressure. In Figure 5(a), the green line is closer to the cyan line compared with other color lines. Thus, it indicates that Eq. (9c) (green line) is more suited than other equations when describing γ. A similar phenomenon is further found in the Dugdale-MacDonald model, as shown in Figure 5 The difference error of γ for the same model reaches the minimum value at V 0 /V =2.5, especially for the cyan line. In the Slate model, the cyan line is close to the green line ( Figure 5(a)). In the Dugdale-MacDonald model, the cyan line exceeds the green line and is close to the red line ( Figure 5(b)).
In the Free-volume model, the cyan line exceeds the green, red, and blue lines ( Figure 5(c)), which shows that γ described by the cyan line has the minimal variance rate with V 0 /V. It is evident that the red and green lines can approximately characterize γ(V ) at higher pressure when the cyan line is taken as a criterion. The calculation of γ 0 is based on the slate model, namely Eq. (8a). The substitution γ 0 into Eqs. (12a)-(12d) produces T m . Due to an adiabatic process under shock compression, the temperature rises sharply to the melting point. As seen in Figure 6, T m shows an increasing trend. It is justified that Eq. (12d) characterized by the blue line could successfully describe T m . 21 From the figure, it is evident that the red line is above the blue curve, and the black and green lines are below the blue line, except for the initial compressive stage of Al. In details, the red line is closer to the blue line when V 0 /V <1.3, and deviates from the blue line when V 0 /V >1.3. It indicates that Eq. (12b) is not suitable for describing the melting temperature of Al at the higher compression. Meanwhile, it is found that the effect of γ 0 on T m is insignificant when the compression is low. Thus, Eqs. (12a) and (12b) can be used to approximately predict T m . When V 0 /V < 1.6, the black line is closer to the blue one. However, the green line is never close to the blue one during the compressive process. The phenomenon indicates that Eq. (12a) is more suitable to describe the melting temperature of Al at the compression stage of V 0 /V < 1.6. Therefore, the red line represents a higher T m , and the black and green lines exhibit a lower T m than the blue line. Also, Eq. (12b) can approximately describe T m when 1.0 < V 0 /V <1.6. Importantly, the effect of γ 0 on T m becomes more significant when V 0 /V increases, because a smaller γ 0 produces a smaller T m . This presents an opposite trend compared to γ. Figure 7 demonstrates the relationship between the pressure and temperature of SC and NC Al with u s of 8.5 -14 km/s. Compared to other previous work, the results of this study show a better agreement with the experimental data 39 and MD results 15 at the pressure within 0 -110 GPa. When the pressure is beyond 110 GPa, the dispersity among the MD simulation results becomes prominent. By increasing the pressure, the temperature predicted by this work is larger than those by the EAM potential (MD1) and the corrected EAM potential (MD2). In addition, there exists a mismatch between the green line and the purple line, suggesting that the phase transition of SC Al from the solid state to liquid one occurs. The green and purple lines represent the solid and liquid Hugoniot curve, respectively. It is also found that the melting behavior occurs when the pressure and temperature are about 107.5 GPa, 3063 K for SC Al, and 109.5 GPa, 3082 K for NC Al, respectively, which agrees with the melting pressure of 100 -130 GPa and the melting temperature of 3100 K provided by others. 40 However, it was reported that the beginning melting pressure was 93.6 GPa for SC Al, 15   not. The difference between this work and previous work might be attributed to the number and size of defects, 15 including dislocations, melting voids, and interactive forces between solid and liquid atoms. The results demonstrate that there is no noticeable difference between SC and NC Al through the same potential. However, there are some differences among three potentials when applying to describe the interaction force between atoms at extreme pressure. Thus, the effect of other potentials will be considered in future work.

V. CONCLUSIONS
In this work, the Hugoniot curve (P H ), the Grüneisen coefficient (γ) and the melting temperature (T m ) of aluminum at high pressure are predicted through the MSST method. The linear relation between the shock wave velocity and the particle velocity as well as the Hugoniot curves have been reproduced in SC and NC Al with grain sizes of 6 and 60 nm which agrees well with the previous studies. It is found that the linear relation and the Hugoniot curve are almost identical for SC and NC Al. This is attributed to because the fact that the grain size cannot nearly affect the dynamic properties of Al at high pressure. γ is found to have a decreasing trend as V 0 /V increases, and contrarily T m exhibits an increasing trend. There is a negative effect between γ 0 and γ, and there is a positive effect between γ 0 and T m when V 0 /V is higher. The melting pressure and temperature are found to be 107.5 GPa, 3063 K for SC Al and 109.5 GPa, 3082 K for NC Al, respectively, highlighting the accuracy of the analytical model. However, it should be mentioned that the findings of this work will be applied to experimental polycrystalline systems, but this would have to be considered in further study. Also, the EOS for different phases of Al needs to be taken into account to check the reliability of such a potential in high pressure and temperature regions.