Atomistic simulations on the interfacial interaction of metallic fuel and structural materials in SFRs - molecular dynamics model for Pu-Fe system

Abstract

Interaction between metallic fuel and steel structures is one of the predominant phenomena in the progress of core disruptive accidents of Sodium-cooled fast reactor. In this study, the atomic diffusion across the interface between Pu and Fe was investigated by using molecular dynamics. The simulation was performed by using Modified Embedded Atom Method (MEAM). The interactions between plutonium and iron atoms were calculated by using the newly developed potential model determined so as to reproduce the material properties of PuFe₂ and Pu₆Fe. The material properties of the compounds predicted with the developed potential were in good agreement with the referenced data. The dissolution or melting at the interface between solid Fe and solid or liquid Pu were simulated by contacting semi-infinite slabs (or liquid layer) of them. Dissolution was observed for all the tested temperature conditions from 800 K up to 1700 K. The melting at the interface was also observed on the interface between solid Fe and PuFe₂ slabs at the temperature approximately 100 K below the melting temperature of PuFe₂ obtained based on the present model.

Keywords:

• molecular dynamics
• eutectic
• modified embedded atom method (MEAM)
• potential model
• melting
• metallic fuel

1. Introduction

Sodium-cooled fast reactor (SFR) is one of the efficient choices of the next-generation nuclear energy supplier. For the analyses of hypothetical core disruptive accidents (CDAs) for such type of reactor, it is of utmost importance to estimate accurately the recriticality induced by the concentration of the molten fuel. When the fuel and structure meet with each other, reaction between metals would occur and induce local complex flow of molten material. Especially the reactions between metals which form the eutectic systems, e.g. U-Fe, Pu-Fe, Fe-Zr, progress at the temperatures even lower than either of the melting temperatures of the pure metals. Thus the evaluation of such metal reactions is a crucial matter to estimate the process of the accidents.

Molecular dynamics (MD), which has been developed especially in these two decades, is one of the powerful techniques for the elucidation of the detailed mechanism of the interaction on the interface between two metals. In the framework of MD, Embedded Atom Method (EAM) [1] and its advanced model, Modified EAM (MEAM) [2], are most commonly used for analyzing the behaviors of metal atoms. Potential parameters for various kinds of metals have been proposed (e.g. [3–7]).

Alloy and eutectic systems have been investigated with EAM or MEAM as well as pure component systems. Webb et al. [8] investigated the melting temperature of Ag-Cu system by using EAM and obtained good consistency with the phase diagram obtained from experiments. They also analyzed intrusion phenomena of liquid Ag into solid Cu plate in their eutectic reaction. Chen et al. [9] performed MD simulation on the diffusion bonding of solid Cu and solid Ag and found that the increase in the system pressure enhances the mutual diffusion at the interface. One of the present authors [10] analyzed the diffusion on the interface between Ag and Cu, and found that the location of liquefaction (i.e. solid-liquid interface), even on microscopic scale, can be estimated by applying the local atomic fraction of the two metals to the phase diagram. He also shows that the surface melting on the interface between the two solid metals remarkably contribute to the mutual diffusion at the interface, which associates with the eutectic melting. Fe-C [11] and Fe-Cu [12] systems have been also examined by several researchers.

However, little attention has been devoted to metallic reaction including nuclear fuel materials. Pu is the only nuclear fuel material for which the potential models based on MEAM have been developed. Baskes [5] proposed the MEAM potential parameters for Pu and investigated the complex behavior of Pu in terms of the distribution of electron density. Baskes et al. [13] calculated the lattice vibration of Pu-Ga alloy by using the potential model developed in the earlier study [14].

In the present work, we focus on the reaction between Pu and Fe metals, the possible reaction to be undergone in CDAs of SFRs. In Pu-Fe system, an eutectic interaction is observed between Pu-metal and Pu₆Fe as shown in Figure 1 [15, 16]. The eutectic temperature is reported to be 686 K, which is lower than the melting temperatures of two terminate materials, such as Pu (913 K) and Fe (1811 K). There are two intermetallic compounds in the Pu-Fe system, PuFe₂ and Pu₆Fe. PuFe₂ melts congruently at 1513 K. Pu₆Fe decomposes peritectically. Between PuFe₂ and Fe-metal, there is another eutectic interaction and the eutectic temperature is reported to be 1440 K.

Figure 1 Phase diagram of Pu-Fe system

The objective of this study is to provide the potential parameters of MEAM for Pu-Fe system and to discuss the mutual diffusion and liquefaction on the Pu-Fe interface. In Section 2, the formalism of MEAM is briefly described. In Section 3, the determination of the potential parameters for the interaction between Pu and Fe atoms is mentioned. The calculational results on the material properties of PuFe₂ and Pu₆Fe are presented in Section 4.1. Section 4.2 introduces the analysis results on the atomic diffusion on the interface between Pu and Fe, and in Section 4.3, the evaluation of the liquefaction at the interface between PuFe₂ and Fe is made. Section 5 gives a summary of this paper.

2. Formalism of MEAM

The analyses in this study are carried out by using the Modified Embedded Atom Method (MEAM) [2]. In the MEAM, the total energy of atoms is described, similar to the original Embedded Atom Method [1], as,

where F(ρ) is the energy to embed an atom i into the background electron density ρ, the background electron density at atom i, Z_i the number of nearest neighbors in what we call the reference structure of a type-i atom, and φ(R_ij ) represents the core-core pair interaction between atoms i and j separated by the distance R_ij . The electron density is obtained through the superposition of angular-dependent ‘partial electron densities’ as,

where ρ ^a(h) _j (R_ij ) (h = 0,1,2,3) is the contribution of atomic electron density of atom j. The parameters α, β and γ represent x, y and z, respectively, and x ^α _ij = R _ij ^α/R_ij . The Equations (2b)–(2d) are equivalent to the power of the cosine dependent of three-body interactions. Equation (2b), for instance, can be rewritten as (ρ⁽¹⁾ _i )²=∑limits _j(≠i)∑limits _k(≠i)ρ ^a(1) _j (R_ij )ρ ^a(1) _k (R_ik )cos θ _jik .

The atomic electron density is given as

where β ^{( h )} is a constant, and R ₀ is the nearest neighbor distance in the reference structure in equilibrium.

The density is related with ρ^(h) _i (h = 0, 1, 2, 3) by the following expressions,

where

The functional form of φ is determined from the total energy per atom in the reference structure, E^u (R), as

where is the background electron density () at an atom in the reference structure, and R means the nearest-neighbor distance. The function E^u (R) is calculated from the universal equation of state,

where E_c is the cohesive energy and

where R ⁰ _i , B_i , Ω _i and E ⁰ _i represent the nearest neighbor distance, bulk modulus, volume per atom and energy in the reference structure, respectively, and δ is an adjustable parameter.

Angular screening is used to limit the interaction range between atoms, as well as the radial cut off often used in MD calculations. With this screening, the interaction between atoms is screened by atoms which are located near the line segment connecting the two concerned atoms. Details of this technique are found in the reference [3].

In this study, Pu was simulated by using the model parameters proposed by Baskes [5]. Plutonium has six allotropes depending on the temperature. This model reproduces the energies and the lattice constants for all the observed solid phases in the experiments. The melting temperature obtained by the model, which is the most important property for the present study, also shows quite good agreement with the experimental value of 913 K.

For Fe, we employed the potential model developed by Lee et al. [6]. This model takes into account the contribution of the second nearest neighbor (2NN) atoms, which is not considered in the original MEAM by Baskes. The form of φ is slightly changed from Equation (6) as

with

where Z₁ and Z₂ are the numbers of the first and the second neighbor atoms, respectively, S is the screening factor, and λ represents the ratio between the second and first neighbor distances. The expression of Equation (10) is identical to Equation (6), for the first neighbor model by Baskes.

The electron density is expressed as,

The energy E^u (R), represented as Equation (7) for Pu in the MEAM, is replaced by the following equation,

Lee et al. [6] showed that this 2NN model reproduces the material properties of Fe quite well such as the lattice constant, melting temperature, elastic constants, etc.

3. Development of interatomic potential parameters

The interatomic potential parameters between unlike (Pu and Fe) atoms were determined so as to reproduce the material properties of the two intermetallic compounds in Pu-Fe systems, i.e., PuFe₂ and Pu₆Fe.

The referred experimental parameters were the formation enthalpy ΔH_f , the lattice constants a ₀ (and c ₀ for Pu₆Fe), and the melting temperature T_m , which are given in Table 1 . The values with asterisks are those referred as the target value for the determination of the parameters in this study.

Table 1 Material properties of PuFe₂ and Pu₆Fe

	−ΔHf (eV/atom)	a 0 (nm)	B (GPa)	Tm (K)
(a) PuFe2 (η cubic)
Reference	0.17^(a)*	0.7189^(c)*	216^(d)*	1512^(f)
	0.094^(b)	0.670^(d)	193^(e)*
		0.681^(e)
Present MEAM	0.19	0.720	207	approx. 1700
Absolute and relative deviation between the	0.02 (11%)	0.01 (1.4%)	9 (4%)	188 (11%)
reference and present model
	−ΔHf (eV/atom)	a 0 (nm)	c 0 (nm)	a 0 /c 0	Tm (K)
(b) Pu6Fe
Reference	0.048^(b)*	1.0405^(d)	0.5349^(d)	1.945^(d)	701^(f)
Present MEAM	0.040	1.082	0.5561	1.945	approx. 800
Absolute and relative deviation between the	0.008 (20%)	0.042 (4%)	0.0212 (4%)	0	99 (12%)
reference and present model

Along with the above parameters, the value of the bulk modulus B is essential for the determination of the potential parameters in MEAM formalism. In this study, the first principle calculations with VASP [20] and WIEN2k [21] codes were performed to obtain the value of B for PuFe₂.

VASP code was used with the spin-polarized generalized gradient approximation (SGGA+U). WIEN2k is based on Full-potential Linearized Augmented Plane Wave Method (FLAPW), and SGGA is also applied in this code. More details for these codes are available in the references shown at the end of this paper.

The basic cell consists of one lattice of PuFe₂ crystal, and 8 × 8×8 and 10 × 10 × 10 of k points were used for VASP and WIEN2k calculations, respectively. The cutoff wave number is set on the basis of k_max = 7/R_mt where R_mt is the smallest atomic sphere radius.

The obtained results of E^U by VASP and WIEN2k are shown in Figure 2 . The results obtained with the two codes show trends similar with each other. The value of B was calculated by using Equation (8) with the constants obtained by fitting Equation (7) to the calculation data. The obtained values of B were 216 GPa for VASP and 193 GPa for WIEN2k.

Figure 2 Relation between total energy E^u and lattice size a ₀ of PuFe₂ obtained with VASP and WIEN2k

The lattice constant and the cohesive energy were also available from the calculation data through similar procedures. The obtained lattice constant (0.670 or 0.681 nm), however, was smaller than the experimental value (0.7189 nm). The calculated cohesive energy (−3.68 or −3.93 eV) results in negative formation enthalpy, which is inconsistent with the experimental observation. In this study, only the value of B is then referred as the target value. The deviations of these calculated values of the lattice constant and the cohesive energy from the experimental ones should indicate the error involved in the estimation of B. However the obtained values of B are close to the experimentally obtained value for UFe₂ (239 GPa) [22], which takes the same lattice structure with PuFe₂ and has a similar lattice constant (0.706 nm) [22] and formation energy (0.080 eV/atom) [18]. These similarities imply the value of B for PuFe₂ to be close to that of UFe₂. It can be then concluded that the obtained values of B would not include essential errors and have enough accuracy for the present study.

Although several attempts were made, no parameter sets have simultaneously reproduced the target properties of PuFe₂ and Pu₆Fe with the conventional model of MEAM for binary system. Furthermore in the case the parameter sets were determined based on the conventional MEAM model, glass-like solidification of Pu-Fe mixture was observed especially at the temperature below 1700 K. To cope with these problems, modification was made on the equation form of the atomic electron density. The pre-exponential factor ρ _0i appearing in Equation (3) was replaced by h-dependent coefficient, namely,

Equation (5) was also modified as

in order to remove the effect of the change in Equation (3′) on the calculations for the pure materials.

As the reference structure, L1₂ structure of the fictitious Pu₃Fe compound was adopted, although the structure of the intermetallic compound is generally taken as the reference structure. In the present case, the test calculation with the reference structure of C15/cF24 structure, which is the structure for PuFe₂, was not succeeded by negative value for the root square calculation appearing in the electron density term. From this result, L1₂ structure, often used for alloy system [12], was then applied. The total energy for one atom in the L1₂ reference structure, is written as follows,

where

In Equation (11) only the first nearest neighbor atoms are taken into account, then the φ_FeFe term does not appear. Note that Equation (12) differs from the original model of MEAM by the change corresponding to Equation (3′), that is, is used for the denominator in the root sign instead of ρ_0i .

The weighting parameter t ^(h) _ave for the interaction with unlike atoms is calculated with the same equation of the original MEAM as

The total energy in Equation (11) is evaluated with the universal equation of state, namely,

where R represents the distance between Pu and Fe atoms.

The obtained values of the potential parameters for Pu-Fe interaction are listed in Table 2 . In the third column of the table, the tangible target parameters or phenomena for the determination of each parameter are also described. The above-described glass-like solidification of Pu-Fe mixture liquid was well suppressed by applying large value (∼3) for ρ⁽³⁾ _0Pu /ρ⁽³⁾ _0Fe . Although this application results in degeneration of the miscibility of Pu and Fe in liquid phase, it can be avoided by lowering ρ⁽⁰⁾ _0Pu /ρ⁽⁰⁾ _0Fe .

Table 2 Potential parameters for Pu-Fe interaction

Parameter	Value	Main contribution
(eV)	3.68	ΔHf of PuFe2
(nm)	0.304	a 0 of PuFe2
α Pu 3Fe	5.6	B of PuFe2
δ Pu 3Fe	0.2	B of PuFe2
ρ^(0) 0Fe /ρ^(0) 0Pu	0.4	Miscibility between Pu and Fe
ρ^(1) 0Fe /ρ^(1) 0Pu	1.5	ΔHf of Pu6Fe
ρ^(2) 0Fe /ρ^(2) 0Pu	0.7	ΔHf of PuFe2 and Pu6Fe
ρ^(3) 0Fe /ρ^(3) 0Pu	3.0	Minimizing solidification of Pu-Fe mixture liquid
Cmax (for all interaction)	2.8	Nominal value
Cmin (Pu-Pu-Fe)	0.8	Minimizing formation of unanticipated structure
Cmin (Pu-Fe-Pu)	0.8	Minimizing formation of unanticipated structure
Cmin (Pu-Fe-Fe)	0.8	Minimizing formation of unanticipated structure
Cmin (Fe-Pu-Fe)	0.36	Assumption (= CminFe )

4. Results

4.1. Material properties of PuFe₂ and Pu₆Fe

The material properties of PuFe₂ and Pu₆Fe predicted with the present potential model are given in Table 1 in comparison with the referenced ones (experiments or first principle calculations). The deviation of the predicted and reference values of ΔH_f , B and a ₀ (and c ₀ for Pu₆Fe) are not more than 11%, except ΔH_f for Pu₆Fe, indicating good agreement between the predicted and reference values. Although the deviation in ΔH_f for Pu₆Fe was 21%, the fact that magnitude of ΔH_f is considerably small as compared to that of PuFe₂ is reproduced with the present model.

The melting temperature T_m was estimated by investigating the advance of melting at the interface between the solid compound and the liquid of the same concentration. The calculation was performed at constant temperatures with zero pressure. The interval of the examined temperature was 100 K. Nosè thermostat [23] was used to keep the temperature at the target value and zero pressure was maintained with Berendsen method [24]. The velocity Verlet algorithm was used with the time step of 1 fs.

The estimated value was approximately 1700 K for PuFe₂ and 800 K for Pu₆Fe, as shown in the table. Both of these values were higher than the experimental ones. One of the reasons for this discrepancy would be that the applied MEAM model for Fe estimates T_m of Fe as 2000 K, which is higher than the experimental value of 1813 K by approximately 200 K.

Figure 3 shows the predicted value of the enthalpy H−H_st where H_st is that at 298.15 K, compared with the experimental one [17]. The present MEAM model reproduces the experimental values for T < 500 K, however, it underestimates for the temperature higher than 500 K. For the temperature between 500 K and 700 K, the experimental result exhibits rapid increase (that is, large specific heat) as compared with those for the other temperatures. This is due to the increase in the specific heat around the Curie temperature of PuFe₂, 633 K, which is not taken into account in the MEAM calculation. This underestimation of the enthalpy in the present MEAM model would be also responsible for the too high melting temperature of PuFe₂ as described above.

Figure 3 Enthalpy of PuFe₂

In the same figure, the calculated enthalpy for η _hex phase of PuFe₂ is also plotted. Although η _hex phase is observed in the temperature range between 1013–1293 K, the calculation was made for the broader temperature range, 700–1500 K. The values of H−H_st were calculated with H_st for the cubic phase. For T ≤ 1100 K the enthalpy for η _hex is slightly lower than that of η _cubic. Since the difference of the enthalpy of the two phases ΔH ≡ H_hex −H_cubic at the transition temperature (1013 K) is associated with the latent heat, ΔH should be positive. It is then known that the present model underestimates H_hex relative to H_cubic at T of 1100 K.

4.2. Simulation of contact of Pu and Fe

In order to explore the atomic diffusion on the interface between Pu and Fe, the simulations were performed by using a system illustrated in Figure 4 , in which pure, semi-infinite Fe slab and Pu slab (or liquid layer) were contacted with each other through an interface parallel to x-y plane.

Figure 4 System geometry for the simulation of contact of two materials. (a) Case 1 and (b) Case 2A–Case 4

The simulation conditions are summarized in Table 3. The simulation temperatures were set at 800 K, 1100 K, 1400 K and 1700 K. The surface orientation of the contacting surface of the solid was (110) or (100). For T = 1100 K (Case 2A and Case 2B), simulations were made for both (100) and (110) orientations of the contact surface of Fe. Periodic boundary condition was imposed on x and y directions for all the conditions.

Table 3 Simulation conditions for the simulation of contact

	Case 1	Case 2A	Case 2B	Case 3	Case 4
Temperature (K)	800	1100	1100	1400	1700
Pu surface	b.c.c. (100)	(Liquid)	(Liquid)	(Liquid)	(Liquid)
Fe surface	b.c.c. (110)	b.c.c. (110)	b.c.c. (100)	f.c.c. (100)	b.c.c. (100)
Lx , Ly (nm)	3.663	2.4616	2.3203	2.5878	2.3396
	2.930	2.0307	2.3203	2.5878	2.3396

Case 1, in which the Pu atoms are also in solid phase, was performed to explore the eutectic point of Pu-Fe alloy observed for X_Pu = 90%. This eutectic reaction could be simulated by contacting Pu₆Fe with Pu; however the contact of pure Pu and Fe is more relevant to the analyses for CDAs then examined. The system temperature in this case was below the melting point of Pu only by 113 K. The bottom atoms of the Pu slab put below the Fe slab (i.e. the atoms on the opposite surface to the Pu-Fe interface) were then constrained to the fixed points in order to avoid unexpected melting of the Pu slab.

For Case 2A–Case 4, the Fe slab was put below the liquid layer of Pu, and the bottom layer atoms of the Fe slab were fixed at constant positions to avoid deformation of the slab.

The cell sizes in x and y directions were kept constant throughout each simulation. The cell sizes for all the cases, except Case 1, were determined based on the equilibrated lattice size of Fe. For Case 1, it was done based on the lattice size of Pu. That is, the cell size L_x and L_y were determined as the lengths corresponding to 10 and 6 lattices of Pu crystal (a ₀ = 0.3663 nm), respectively. The value of L_x does not match the longer side of the equilibrated lattice size of Fe (0.414 nm for (110) surface of body-centered cubic lattice), the Fe lattice was then distributed with reduced length of the longer side to 0.407 nm (98.2% of equilibrium size), which leads to the alignment of 9 lattices of Fe in the x direction of the cell.

The Fe slab and liquid layer (or slab) of Pu were separately equilibrated at the target temperature with zero pressure in the preliminary calculations. The simulations were then initiated by contacting the surfaces of the equilibrated Pu layer (or slab) and Fe slab.

Shown in Figure 5 are the instantaneous distributions of the atoms projected onto x-z plane at t = 2000 ps or 4000 ps. As shown in the figure, the Pu and Fe atoms penetrate into the bulk of the other metal. This indicates that the atomic diffusion with dissolution or melting progressed at the interface.

Figure 5 Distribution of atoms in the simulation of contact. Black and gray particles represent Pu and Fe atoms, respectively. (a) Case 1 (T = 800 K, t = 8000 ps), (b) Case 2A (T = 1100 K, t = 8000 ps), (c) Case 2B (T = 1100 K, t = 2000 ps), (d) Case 3 (T = 1400 K, t = 4000 ps) and (e) Case 4 (T = 1700 K, t = 2000 ps)

For Case 1, in which both the metals are in solid phase, the distributions of the atoms near the interface are disturbed, as seen in surface melting phenomena [25], which is the phenomena wherein the melting occurred in the region molecularly close to the surface (or interface) even below the melting point. In the corresponding temperature (800 K), no surface melting was observed both for the Pu and Fe surfaces on their free surfaces before the contact. It can be thus concluded that the contact enhances the surface melting on the interface between them. This should have served as promotion of the mutual diffusion (namely, mutual dissolution) at the interface as if they are contacted in the liquid phases. Such surface melting is the essential mechanism for the initiation of the eutectic melting in which the two contacted metals melt below either of their melting temperatures [10].

For Case 2A, Case 3 and Case 4, the Fe atoms diffused (dissolved) into the liquid Pu layer more prominently than the Pu atoms did into the solid Fe. The thickness of the mixing region, in which both Pu and Fe atoms coexist, spreads more rapidly than in Case 1. In the highest temperature case, Case 4, almost all the Pu and Fe atoms mixed with each other by 2000 ps.

In Case 2B, the Pu atoms adjacent to the interface formed clear coherent crystalline structure with Fe. Such ‘epitaxial’ crystallization can be seen when the lattice size in the solid state of the corresponding material [26] or the number density of atoms [27] is close to that of the other material. For the present case, epitaxial crystallization was developed for only one layer. Thus by taking into account only the distance between the neighboring atoms on the interface, the distance of the nearest neighbor atoms of Pu in β phase (0.297 nm) is close to the interval of the Fe atoms (equal to the lattice size of 0.290 nm). Then such solidification by the epitaxial crystallization is likely to be developed. This solidification of Pu at the interface prevented the diffusion across the interface. As a result, the Fe atoms diffused into Pu bulk only in early period (∼100 ps) before this epitaxial layer of Pu was developed. Such epitaxial crystallization was seen even when the potential parameters for Pu-Fe interaction were widely varied. No experimental observations, however, have been reported on such crystallization, thus further investigation should be performed for the clarification of this phenomenon.

Figure 6 depicts the profiles of the atomic concentration of Pu, X_Pu , along z axis for all the cases except Case 2B. The distribution of X_Pu was calculated as the number density concentration of Pu atoms, in each layer with a thickness of 0.025 nm parallel to x-y plane over the whole cross section of the simulation cell. The values presented in the figure are those averaged for 10 ps around the corresponding time. The magnitude of the gradient of X_Pu in the mixing region (0 < X_Pu < 1) decreases with time. This indicates that the diffusion progresses with time. The decreasing rate of the gradient apparently depends on the temperature. That is, the gradient of X_Pu in Case 4 (T = 1700 K, the highest temperature case) decays more prominently than those for the lower temperatures.

Figure 6 Profiles of the atomic concentration of Pu (X_Pu ). (a) Case 1, (b) Case 2A, (c) Case 3 and (d) Case 4

The time evolution of the diffusion can be clearly represented by the thickness of the mixing region, η _10-90, where 0.1 < X_Pu < 0.9. The values of η _10-90, as shown in Figure 7 , increased with time. The values of η _10-90 for higher temperatures increase more rapidly than those for the lower temperatures. For all the cases the increasing rate of η _10-90 (i.e. dη _10-90/dt) gradually declined with increase in time; however for t > 2000 ps in Case 3 and t > 1000 ps in Case 4, dη _10-90/dt fell rapidly to zero. These rapid falls should be caused as the result that the concentration of Fe atoms reached the saturation concentration in liquid Pu for Case 3 and that all the Fe atoms dissolved and diffused over the whole region of liquid Pu for Case 4.

Figure 7 The time evolution of η _10-90. Dashed lines represent fitting by Equation (16)

For the other period, the processes which contribute to the variation in η _10-90 with time should be (1) the dissolution process of the Fe atoms into Pu and (2) the diffusive motion of the atoms. By assuming that the diffusive motion is more responsible to η _10-90 than the dissolution process, the evolution of η _10-90 is expected to follow the parabolic rule [28], that is,

where k is a constant which has the same dimension of diffusion coefficient. In Equation (16), η ₀, which does not appear in [28], represents the rapid increase in η _10-90 in the beginning period of the contact of the two materials. This rapid increase should be associated with the disturbance of the atom motions caused by the their contact.

Equation (16) was fitted to all the data of η _10-90 for each case shown in Figure 7, except t > 2000 ps for Case 3 and t > 1000 ps for Case 4. As shown in the figure by dashed lines, the evolution of η _10-90 was well fitted by Equation (16) for all the cases. This indicates that the rates of the dissolution or melting at the solid-liquid interface were dominated by the diffusive motion of the Pu and Fe atoms. This result also implies that the decrease in the increment rate of η _10-90 with time was primarily owing to the decrease in the magnitude of the spatial gradient of X_Pu , that is, ∂X_Pu /∂z.

In Figure 8 , the obtained values of k are plotted versus inverse of the temperature. As shown in the figure, the logarithm of k linearly decreases with T ⁻¹, indicating that k obeys the Arrhenius form. It is worth noting that k for Case 1 (10⁴/T = 12.5 K⁻¹) also meets the fitting line drawn in the figure, though the initial configuration in Case 1 is the coalescence of two solid phases. This indicates that, in Case 1, the melting on the solid-solid interface was also dominated by the diffusion in the liquid phase. For this indication it is required that the region of the surface melting should have expanded before the atoms of the other metal diffused to the point where liquefaction took place. This is consistent with the observation in Figure 5(a) that the surface melting occurred and the diffusion takes place in the melting region. However for this case the evolution rate of surface melting and diffusion is quite low, so further simulation should be needed to estimate the phenomenon in the macroscopic practical system.

Figure 8 Temperature dependence of k in Equation (16). Black circles are MD results, and dashed line is the best fit

4.3. Simulation of contact of PuFe₂ and Fe

Pu-Fe system has the other eutectic temperature around X_Pu ≈ 0.2 as shown in Figure 1. This eutectic point was examined by the simulation of the contact of PuFe₂ compound and pure Fe slabs. The simulation procedures were similar to that described in Section 4.2. The slab of PuFe₂ was put above the Fe slab with the bottom atoms constrained at the fixed positions. The top surface of PuFe₂ slab was free surface. The surface of Fe slab was (100) of face-centered cubic (f.c.c.) lattice.

The simulations were performed for two temperature cases of 1500 K and 1600 K, both of which were lower than the melting temperature of PuFe₂ simulated by the present model (approx 1700 K). The lengths of the cell, L_x and L_y , are 2.9128 nm and 2.9128 nm for T = 1500 K and 2.9232 nm and 2.9232 nm for T = 1600 K, both of them corresponding to 4 lattices of PuFe₂ crystal in x and y directions, respectively. The consequent deviation of the length of the distributed Fe lattice in the corresponding cell sizes from that of equilibrated Fe lattice was below 0.8%. The slabs of Fe and PuFe₂ were separately equilibrated in the preliminary calculations with the positions of the atoms in the top layer of PuFe₂ slab and in the bottom layer of Fe being fixed.

Figure 9 depicts the instantaneous distributions of the atoms projected onto x-z plane at t = 0 ps and 3000 ps. As seen in Figure 9(a) (T = 1600 K), surface melting was observed at the surface of Fe, which occurred before the contact with the PuFe₂ slab. On the other hand the surface of PuFe₂ did not melt at t = 0 ps. The distribution of ρ_N for 1500 K was quite similar to Figure 9(a).

Figure 9 Distribution of atoms projected onto x-z plane in the simulation of contact of PuFe₂ and Fe. Black and gray particles represent Pu and Fe atoms, respectively. (a) 0 ps (T = 1600 K), (b1) 3000 ps (T = 1500 K) and (b2) 3000 ps (T = 1600 K)

At t = 3000 ps, for T = 1500 K (Figure 9(b1)), a few Pu atoms which were contained in PuFe₂ at the initial distribution moved into the Fe slab, however the crystalline structure of Fe slab was still maintained. On the other hand, for T = 1600 K (Figure 9(b2)), the crystalline structure of Fe near the interface became unclear and more Pu atoms diffused into the Fe slabs. This indicates that the present model reproduced the existence of the eutectic melting point between Fe and PuFe₂ around 1600 K.

In Figure 10 the profiles of the number density ρ_N at 100 ps (Figure 10(a) for T = 1600 K) and at 3000 ps (Figures 10(b1) and (b2) for T = 1500 K and 1600 K, respectively) are shown. The initial position of the interface between the Fe and PuFe₂ slabs was z ≈ 2 nm. In all the figures, periodic changes in ρ_N can be seen, which represent the layered distribution of the atoms. This layer structure should be mainly associated with the crystalline structure of Fe and PuFe₂. For t = 100 ps (Figure 10(a)), although surface melting of Fe should advance near the interface, the clear layered distribution of the Fe atoms is maintained.

Figure 10 Distribution of the number density of the atoms, ρ_N . (a)t = 100 ps (T = 1600 K), (b1) t = 3000 ps (T = 1500 K) and (b2) t = 3000 ps (T = 1600 K)

At t = 3000 ps, the distribution was almost unchanged from that at t = 100 ps for T = 1500 K (Figure 10(b1)) except for a few Pu atoms which moved to z ≈ 1.5 nm in the Fe slab. While, for T = 1600 K (Figure 10(b2)), ρ_NPu in 1.1 nm < z < 2.2 nm increases up to 10 nm⁻³, indicating the diffusion of the Pu atoms into the Fe slab. More importantly, for 1.3 nm < z < 3 nm the amplitude of the variation of ρ_NFe in this area is remarkably decreased from that in the other area. This moderation of the layered distribution indicates that the slabs in this area were melted into liquid phase and the crystalline structure was disintegrated.

Such disintegration of the structure can be also quantified with the S factor defined as [29]

where r _i is the x-y position of atom i, k_m is the reciprocal space lattice vector (k_mα = 2πm _α/L _α, with m_α = 1,2, … and α = x or y) in the x-y plane and N is the number of the atoms in the concerned area. When the atoms form a perfect lattice, 〈S_m ²〉 equals to one with m_α = νL_α /a ₀ where ν is an arbitrary positive integer. Shown in Figure 11 is 〈S_m ²〉 for Fe obtained for each x-y layer distribution with m_α = 8 corresponding to L_α /m identical to a ₀ of Fe lattice and a ₀/2 of PuFe₂ lattice.

Figure 11 Distribution of S factor 〈S_m ²〉 for Fe defined in Equation (17). (a) T = 1500 K and (b) T = 1600 K

For T = 1500 K, as shown in Figure 11(a), the distribution of 〈S_m ²〉 is almost unchanged from t = 100 ps until 3000 ps. This means that the crystalline structure was maintained except in the vicinity of the interface (1.7 nm < z < 2.1 nm) where surface melting on the surface of the Fe slab occurred before the their contact. For T = 1600 K, as shown in Figure 11(b), 〈S_m ²〉 in 1.3 nm < z < 3 nm at t = 3000 ps remarkably decreased from that at t = 100 ps, indicating that the lattices in this region disintegrated (melted) and the atoms behaved as liquid phase. This estimated range of the melting region agrees with that in which the layered structure presented in Figure 10 is moderated.

5. Conclusion

Atomic diffusion across the interface between Pu and Fe was investigated in the framework of MEAM. The applied potential parameters for Pu and Fe were those proposed by Baskes [5] and Lee et al. [6], respectively. The potential parameters for Pu-Fe interaction were developed so as to reproduce the material properties of the two compounds, PuFe₂ and Pu₆Fe, measured in the experiments or calculated with the first principle calculations by using VASP and WIEN2k codes. The material properties of the compounds predicted by the developed potential are in good agreement with the experimental data.

The dissolution or melting at the interface between pure Pu and Fe was simulated by contacting semi-infinite slab of Fe with liquid layer of Pu or slab of Pu. Dissolution was observed for all the tested temperature conditions from 800 K up to 1700 K. For the simulation with T = 800 K, in which solid Pu and Fe slabs were contacted with each other, the surface melting occurred at the interface, which enhanced the diffusion at the interface. The time evolution of the thickness of the mixing region, in which 0.1 < X_Pu < 0.9, strongly depended on the temperature, and was found to be dominated by the atomic diffusion in the liquid phase even in the simulation for the contact of the two solids.

The melting at the interface was also observed on the interface between Fe and PuFe₂ in solid phase at the temperature approximately 100 K below the melting temperature of PuFe₂ obtained based on the present model. This indicates that the model reproduces the existence of the eutectic point between Fe and PuFe₂.

Acknowledgements

A part of this study was carried out within the task ‘R&D of the Next Generation Safety Analysis Methods for Fast Reactors with New Computational Science and Technology’ entrusted from the Ministry of Education, Culture, Sports, Science and Technology of Japan.

Notes

b.c.c: body-centered cubic, f.c.c.: face-centered cubic.