Technical Comments

Shock Wave-Induced Melting in Argon by Atomistic Simulation

See allHide authors and affiliations

Science  21 Nov 1997:
Vol. 278, Issue 5342, pp. 1474-1476
DOI: 10.1126/science.278.5342.1474

A. B. Belonoshko (1) considers an interesting and important problem of atomistic simulation of shock wave–induced melting. On the basis of molecular dynamics (MD) simulations with empirical Buckingham potential (2, 3) of the argon Hugoniot, Belonoshko (1) found two discontinuities “that may bracket a mixed-phase region of solid and liquid” and concludes that “this is an intrinsic feature of the Hugoniot crossing the Ar melting curve and does not require the addition of any solid-solid phase transition.” This conclusion is not new [see, for example, figure 9 in (4) for Al and compare with figure 2a in (1)], but the applied method (1), the generalizability of the results, and the implications for the analysis of shock-wave experiments on iron are noteworthy. There are, however, several problems with this report (1).

First, the result is based on the incorrect statement (1, p. 955) that “argon transforms to liquid directly from face-centered-cubic (fcc) phase without any high-pressure [ P ] and high-temperature [ T ] solid-solid phase transformations …”; and Belonoshko makes a reference toexperimental work (5) in this regard. But modeling with empirical interatomic potentials (IP) requires careful analysis of the possible stable configurations of atoms. To understand it, I have plotted the energies of fcc (three-layered structure in terms of closed-packed structures), dhcp (double-hexagonal, closed-packed, four-layered structure), and bcc (body-centered-cubic) structures as function of volume (Fig. 1). While the energy of the bcc structure is always higher than the energy of the closed-packed structures, the energy of one closed-pack structure cannot differ from that of another (6). This is easy to understand, because Buckingham potential (2) is a short-range potential, and all closed-packed structures have exactly the same first two coordination shells (7). I made calculations for all possible closed-packed structures up to 10 layers and got the same results. Consequently, solid “Ar” [now we are talking about modeling a system with IP (2) that can reproduce some properties of real Ar] at given P and T can have many different polymorphs. Even if the starting configuration was fcc, shock wave can produce other polymorphs, and these solid-solid phase transitions can be responsible for discontinuities on simulated Hugoniot. Belonoshko's report (1) does not contain an analysis of such possibilities (8).

Figure 1

Total energy of argon at T = 0 K as function of volume calculated with Buckingham potential (2) for fcc, dhcp, and bcc structures. Lines for closed-packed fcc and dhcp structures are not distinguishable.

Belonoshko (1) states that the calculated shock wave velocity values are somewhat larger than the experimental values because of the initial higher density of Ar in his simulations. Then he finds that “the calculatedU p-U s data agree with experiment” (9), where U p is piston velocity. But overestimation of the difference (U p-U s) in the report (1) gives overestimation of bulk modulus of solid argon by 25 to 35%. This is an unsatisfactory result, especially because the Buckingham IP (2) allows one to reproduce both static and dynamic compression data accurately (2). The source of error is an unrealistic configuration of the computational cell (1). In shock-wave experiments, the samples contain a large amount of crystallites in random orientation. In contrast, the initial configuration of atoms contains an ideal fcc Ar crystal (1) with the [100] normal to propagation direction of the shock wave. As a result, Belonoshko (1) studied elastic properties of Ar in the [100] direction only. Elastic constants of Ar at different P can be calculated (Table 1). C11 is ∼30% bigger than the bulk modulus KT.

Table 1

Elastic constants of Ar at 1000 K and different pressures.

View this table:

An unrealistic computational cell configuration (1) can produce effects of recrystallization under stress conditions. It was shown (11), for metals with fcc structure by a combination of MD and Monte Carlo (MC) calculations, that under stress applied normally to [100], atoms change their positions, and the whole computational cell recrystallizes in a configuration with [111] normal to the loading direction. Following the procedure described by Selinger et al. (11), I made calculations for Ar with the Buckingham IP (2) and found that it also recrystallizes under uniaxial stress. This result may be obtained even without numeric modeling, because when stress is applied along the [111] direction of a cubic crystal, the energy of deformation is less than that for the same stress applied in the [100] direction. Recrystallization should affect the calculated Hugoniot (1) and, in particular, mimic real melting and produce mixed region on Hugoniot, because “a small increase in U p does not provide sufficient energy to permit the sample to jump to the liquid branch of the Hugoniot” (1) because energy is spent on recrystallization.

The T of liquid Hugoniot was underestimated (1) by at least 2000 K in comparison with the experimental data or with the previous calculations (2). This error may be the result of a wrong application of the periodic boundary conditions (1). Indeed, shock wave disturbs the boundary between periodically repeating elementary cells in the x and y directions; in the report by Belonoshko (1), these have sizes from 5×5×60 to 10×10×120 unit cells (or, approximately, from 25×25×300 Å to 50×50×600 Å at 10 GPa). As a result, “computational samples” in the report (1) are systems of small particles. It has been shown both experimentally and theoretically that materials with small (less than 500 Å mean size) particles have different phase relations than bulk phases and, in particular, T of solid-solid and solid-liquid phase transitions may change by several hundred degrees (12). My calculations show that even simple dislocations, with Burgers vector ½<110> in fcc Ar lattice repeated each 50 Å at 150 kbar, decrease melting T by ∼900 K. Belonoshko's conclusion (1) that “melting occurs without any so-called “overshooting” (when a substance metastably remains solid even when its T is higher than the melting T) also may be wrong, because melting was observed in a nonrealistic supermicrocrystalline phase of Ar.

Even though inert gas Ar and metal iron (Fe) are different in nature (13), Belonoshko (1) applied his analysis of shock wave–induced melting of Ar for iron (Fe). But more significantly, there is no evidence (1) of two discontinuities on the curve Vp( P ). Only one point in figure 3 in the report (1) represents “second discontinuities,” and this could be a result of statistical fluctuations in the numerical calculations [the point at ∼50 kbar on figure 3 in (1) is also offset from the smooth line, for example] or a result of artifactual melting or recrystallization.

Thus, I conclude that the suggested method of atomistic simulation of shock wave–induced melting (1) is not correct and cannot be applied to both Ar and Fe.


Shock Wave-Induced Melting in Argon by Atomistic Simulation

Response: I am happy that L. S. Dubrovinsky expresses interest in my report (1). He raises a number of issues with regard to it: (i) there could be other close-packed polymorphs responsible for the discontinuities on the simulated Hugoniot because the potential (1) is short-ranged; (ii) the computational cell in my numerical experiment has unrealistic configuration because it is a single crystal; (iii) whether melting or recrystallization is not adequately addressed; (iv) the cell was too small for a simulation of a shock wave, and simulation of melting requires much larger cells; (v) the liquid Hugoniot is underestimated by at least 2000 K; and (vi) the results of the simulation on Ar cannot be used in the case of Fe. My responses to these issues follow.

1) Because the energies (E) of all close-packed modifications (fcc, hcp, dhcp, and so on) as a function of volume are within the error of calculations [figure 1 in (1)], there is no reason to expect solid-solid structural transformations. The Rankine-Hugoniot relations (equations 1 and 2 in my report) represent dependence between P, V,U s and U p. BecauseEmbedded Imagewhere T is a constant temperature, the transformation from fcc to dhcp and similar transitions will not produce any discontinuities on the Hugoniot if the E(V) function is the same for fcc and dhcp polymorphs [as it is for the short-range potential (1)]. I chose fcc as the starting structure because Ar is stable in the fcc phase according toexperiment (2). The potential reproduces properties of fcc phase accurately (3). I monitored the structural information during all simulation runs, and no hcp or dhcp configurations was observed. The molecular dynamic simulation provides complete information on the structure, including the coordinates of all atoms at any time during a simulation run; more complete information than any experiment. The statement that an analysis of the radial distribution function (RDF) does not help in such cases (reference 8 in the comment) seems equivalent to rejection of any attempt to obtain an x-ray pattern of shocked material, because a structure factor is a Fourier transformed RDF. Because the energies of all closed packed structures are essentially the same (figure 1 in the comment), there is no reason for appearance of other polymorphs from initial fcc structure.

2) Simulations of a sample that contains a number of particles comparable with those in a real sample are not possible at this time for technical reasons. Therefore, caution should be taken in interpreting simulated data. Because the simulated sample is an oriented crystal and there is no shock-wave data on Ar single crystals, one has to take into account the possible influence of orientation. There are experimental data on single and polycrystalline Ni samples, which show that U S in a single crystal can be lower by 10% (4) and higher by 17%, depending on the orientation of the sample to the propagating shock wave, as compared with those in polycrystalline samples. The simulated data are in agreement with experiment, if one takes into account the effect of orientation as assessed from data on Ni.

3) When material is subjected to high strain beyond the yielding limit, a liquid-like structure appears, as observed in experiments (5) and in simulations (6). The appearance of such a structure is responsible for the first discontinuity on the simulated Hugoniot (1). With respect to the following analysis, it is not essential what exactly occurred—re-crystallization or melting. Structures in both cases can be indistinguishable (5,6). The important point is that this is definitely not a solid-solid transition, and there is a discontinuity on the Hugoniot.

4) In my report (1), I state that calculations have been checked by doubling the size of the cross section, and no significant changes of results were observed. Because the dependence of the results on the number of particles (in this case, on the number of unit cells in the x and y directions) is asymptotic (the results do not change after exceeding some sufficiently large number), it is clear that if I did not see any changes of the results when doubling the size from five to ten unit cells, then the results would not change with further increase in size. The Lindemann (7) criteria for melting is a first approximation, and for Ar it is a good one (8). Mean square displacement of atoms [this is within the Lindemann criteria; see (9) for details] depends on the number of unit cells in a computational cell of cubic shape (Fig. 11). Ten cells are more than enough for obtaining reliable results. Dubrovinsky's reference (reference 12 in the comment) to the data, without providing details, is not convincing because properties of small particles could be different from bulk properties for several reasons (for example, the structure of adsorbate is that of the adsorbent).

Figure 1

Mean square displacement (MSD) of Ar atoms with the initial configuration as n × n × n fcc unit cells at P = 50 kbar and T = 300 K. Statistical errors are rather large at smaller n (from 2 to 5; about 5%) and quite small at larger n (at 7 and 30, the errors are less than 1%). Minimum at about n = 4 is in accordance with theory (14). The MSD asymptotically approaches the limit, and the difference between the limit and MSD atn > 7 is small (less than 3%).

5) The calculated melting temperature from the Hugoniot (that is, one of the points on the liquid Hugoniot) is in agreement with experiments and calculations, and Dubrovinsky is not correct on this point. Ross (10, table 5) gives T = 4235 K at P = 22.53 GPa as a point on the liquid branch of the Ar Hugoniot. This is the highest pressure within the range of my calculations. Figure 2A in my report gives T = 3700 K at the same P. Therefore, the difference is about 500 to 600 K, which is definitely less than “at least 2000 K” and can be accounted for by different potentials and methods used in the calculations.

6) Argon and iron are, of course, different substances. In my report, I reasoned that if two discontinuities on the Ar Hugoniot (1) do not require an assumption about solid-solid transition, it is not necessary to presume that the two discontinuities on the Fe Hugoniot (11) are the result of a solid-solid transition. Two discontinuities can exist without a solid-solid structural transition of the material, in accord with the theory of shock waves in solids. In addition, my simulations (1) showed the approximate size of discontinuities one could expect in Ar, which was comparable with the size of discontinuities in Fe (11). The corresponding figures showing discontinuities in curves of velocities of rarefaction waves in Ar (1) and in Fe (11) are similar (12). This comparison supports my original conclusion (1).


Navigate This Article