ISIJ International
Online ISSN : 1347-5460
Print ISSN : 0915-1559
ISSN-L : 0915-1559
Molecular Dynamics Simulations of Iron/Graphite Interfacial Behaviors: Influence of Oxygen
Yajun YinWen LiHujun ShenJianxin ZhouHai NanMingsen DengXu ShenZhixin Tu
Author information

2018 Volume 58 Issue 6 Pages 1022-1027


Graphite morphology and the liquid flow behaviour in lower part of blast furnace show a close relationship to the interfacial energy between iron and graphite. In this study, molecular dynamics simulations were employed to explore the position preference of an anti-spheroidizing element (such as oxygen atom) in melt cast iron as well as its influence on the iron/graphite interfacial energy. First, our results showed that oxygen atoms were more favorably adsorbed on the prism plane of graphite than on the basal plane. Second, MD study revealed that the adsorbed oxygen atoms would influence the wettability between molten iron and graphite prism plane, leading to a significant change in interfacial energy. For instance, the interfacial energy decreased with increasing number of adsorbed oxygen atoms.

1. Introduction

Graphite morphology is essential to the performance in the process of iron casting. It has been acknowledged that spheroidal graphite is a traditional graphite morphology widely existing in the pure Fe–C–Si alloys.1,2) Some chemical elements, including sulfur(S), oxygen(O), magnesium(Mg), and cerium(Ce), are considered to be important factors influencing the morphology of graphite, in particular oxygen and sulfur are well-known anti-spheroidizing elements.3,4,5,6) Thus, controlling the chemical compositions appears to be critical in the industry to obtain the spheroidal graphite in casting. A variety of experiments have been conducted to study the effects of altering chemical compositions on the graphite morphology. In 1970s, some researchers found that some anti-spheroidizing elements (such as S and O) tend to be adsorbed at the iron/graphite interface, directly affecting the morphology of graphite. More specifically, graphite tends to be flake-like when these anti-spheroidizing elements are adsorbed at the iron/graphite interface, and otherwise it would become spheroidal.7,8,9) Although substantial efforts have been made to reveal the mechanism of spheroidal graphite formation, our understanding of spheroidal graphi te is still limited, leading to various debatable theories. For instance, one leading theory is that the graphite morphology relies on the preferred direction of graphite growth and the orientation of graphite growth would be greatly affected by the anti-spheroidizing elements adsorbed at the iron/graphite interface.10,11) In another theory, some researchers believe that the spheroidal graphite can be formed by curving stacked flake-like graphite12,13,14,15) and the interfacial energy can be lessened by reducing the bending capacity of flake-like graphite owing to the presence of the anti-spheroidizing elements adsorbed at the iron/graphite interface. Although two debating theories exist, both of them agree with the importance of interfacial energy (or interfacial tension). To explore the role of interfacial energy and its influencing factors, lots of experiments have been undertaken, demonstrating that surface tension of melt iron can be influenced by impurities (such as C, O, S),16,17,18) and graphite tends to be flake-like as the iron/graphite interfacial energy decreases in the presence of anti-spheroidizing elements,19,20,21) implying the close relationship among influencing chemical elements, iron/graphite interfacial energy and graphite morphology. Moreover, wetting between liquid iron and coke influences liquid flow in the lower part of a blast furnace, which strongly affects the operation of the furnace.22,23)

Two hundred years ago, Thomas Young proposed a method by treating the contact angle of a liquid lying on a solid surface as the result of the mechanical equilibrium of a droplet resting on the solid plane.24) Although a sessile drop method can be used to measure interfacial energy,25,26) it is still challenging for this method to reveal the influence of interfacial energy on the graphite morphology under various surface conditions. Fortunately, Molecular dynamics (MD) simulation, as an important complementary tool to experiment, has been proved to be practical, since it’s convenient to control the surface condition in MD simulation.27,28) In particular, MD simulations have been extensively used for the wettability research since 1995,29) providing comparable results to experiment as well as molecular kinetic theory.30) Since then, numerous works by MD approach have been done to disclose the connection between liquid/solid interaction energy and wettability,31,32) the dewetting properties of Cu drop on nanopillared graphene,33) the structure and dynamic behaviors of two immiscible liquid metal on graphene,34) and as well as to calculate the solid-liquid contact angle with a 3D configuration method.35) In addition, the coupling of experimental and MD methods were recently used to investigate the wettability of water on different substrates by Rishi Raj et al.36) and Javad Rafiee et al.37) Unfortunately, few studies have been done to explore the effects of anti-spheroidizing elements on the iron/graphite interfacial energy from a microscopic level, such that there is limited understanding of the physics behind the iron/graphite interfacial behaviors. In this work, MD simulations have been carried out to advance our understanding of the role of oxygen as the anti-spheroidizing element influencing the iron/graphite interfacial energy and our studies will benefit to the research improving the cast iron performance.

2. Computational Details

To obtain accurate results from MD simulations, reliable force fields are required. Currently, various force fields have been parameterized and improved to describe a variety of systems.38,39,40) These force fields, relying on the employed potential functions, have their advantages and weakness, such that it has become attractive to combine different potential functions to maximize the performance of MD simulations. embedded-atom method (EAM) is a model which is particularly appropriate for metallic system. Reactive force field (ReaxFF) potential uses distance-dependent bond-order functions to represent the contributions of chemical bonding to the potential energy, which is suitable for the C–O interaction. In this work, the EAM potential41) and ReaxFFpotential42,43) were used for describing the interactions between Fe atoms in the metal cluster and for the C–O interactions respectively. Since the states of oxygen (molecular oxygen and atomic oxygen) in the molten iron are complicated and it was found that molecular oxygen would have negative impact on graphite morphology,44) we assume that oxygen exists as a free atom. In melt iron, the Fe–C and Fe–O interactions are van der Waals forces, and can be well described by Lennard-Jones (LJ) potential. The LJ potential has been extensively used for the Fe–O interactions in earlier works.45,46,47) Thus, in this work the LJ potential was employed for the Fe–O and Fe–C interactions. The LJ potential is given by   

E=4ε[ ( σ r ) 12 - ( σ r ) 6 ]  r< r c (1)
Where ε is the depth of the potential well, σ is the finite distance at which the inter-particle potential is zero, and rc is the cutoff. All LJ parameters are given in Table 1.
Table 1. Parameters of LJ (12-6) Potential.45,46,53)
ε (eV)σ (Å)

All the MD simulations were carried out with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).48) For each case, energy minimization was performed and followed by heating systems from 0.1 K to 1900 K under NVT condition, and a subsequent NVT production run was carried out. During the NVT simulation, the Nose-Hoover method49,50) was employed to govern the constant temperature at 1900 K. The integration time step of 0.2 fs was adopted for all MD simulations and a cutoff value of 12 Å was set for van der Waals (vdW) interactions. Visual Molecular Dynamics (VMD)51) and the Open Visualization Tool (OVITO)52) programs were used for visualization and analysis.

3. Result and Discussion

3.1. Adsorption of Oxygen

To investigate the position preference of oxygen atoms adsorbed on graphite surfaces, we built two types of graphite surface: the prism plane [0001] and basal plane [1010]. The Fe cluster, containing 2691 Fe atoms (with body-centered cubic(BCC) structure) and 36 oxygen atoms, was placed above the two graphite surfaces, as shown in Fig. 1. In the Fe cluster, oxygen atoms were evenly or randomly placed above the prism and basal planes, as respectively shown in Figs. 1(a) and 1(b). Then, the size of the simulation box was set to 2.7 nm*3.3 nm*5.3 nm for all MD simulations. Oxygen atoms were fixed during the heating process and then were released in the NVT production run for 120 ps while the graphite was fixed in the entire simulations.

Fig. 1.

Initial structure of oxygen adsorption MD simulation model with and without Fe atoms. (Online version in color.)

Snapshots were taken from two independent MD simulations, showing the adsorption process of oxygen atoms on the two different graphite surfaces, presented in Fig. 2(a). From the top panel of Fig. 2(a), it is seen that oxygen atoms, evenly placed above the planes in the beginning of MD simulation, have a preference of adsorption on the prism plane over the basal plane. Similarly, from the bottom panel of Fig. 2, one can see that oxygen atoms, randomly placed above the planes in the beginning of MD simulation, have the same preference of adsorption on the prism plane over the basal plane. Thus, MD results suggest that the prism plane can arrest oxygen atoms more easily than the basal plane (shown in Fig. 2(b)), in consistence with the observation from experiment.7,8,9) In fact, this is reasonable because the prism plane provides more dangling bonds for oxygen adsorption than the basal plane.54,55) Therefore, when oxygen concentration is low (and the solubility of oxygen in melt iron is low), there are much more oxygen on the graphite prism plane. In addition, MD simulations reveal two dominant locations of oxygen atoms preferentially adsorbed on the prism plane: one (type I) is located above the concave of graphite surface and the other (type II) is located above the edge of two carbon atoms, shown in Figs. 3(a)–3(b). The distances between adsorbed oxygen atoms and their closest carbon atoms were measured around 1.4 Å (C–O single bond is about 1.36 Å), suggesting that this is the chemical adsorption.

Fig. 2.

The oxygen motion snapshots of two MD simulations (a) and the number of oxygen atoms on different planes (1 for oxygen evenly placed and 2 for randomly) (b) at the temperature of 1900 K. (Online version in color.)

Fig. 3.

(a) Schematic of oxygen adsorption (b) number of adsorbed oxygen atoms. (Online version in color.)

3.2. Iron/graphite Interfacial Energy and the Influence of Oxygen

Figure 4 shows two types of interface for liquid droplet on a solid substrate.24,56) The solid will not dissolve into the liquid and the interface will remain flat when the liquid is saturated under thermodynamic equilibrium conditions (Fig. 4(a)). In this case, the Young’s equation allows the calculation of interfacial energy as follows:   

γ SV = γ LS + γ LV cosθ (2)
γSV, γSL and γLV means substrate-vapour, liquid-vapour and liquid-substrate interfacial energy respectively. If the liquid is not saturated, then the equilibrium shape involves a finite curvature for the solid-liquid interface (Fig. 4(b)) due to the dissolution of the substrate, and in this case Young’s equation is no longer applicable. The three interfacial energy contributions are then related to the dihedral angles formed by the planes tangent to the interfaces through the sine rule:   
γ LV sin Ø S = γ SV sin Ø L = γ LS sin Ø V (3)
In this work, the graphite solid substrate was treated as a fixed plane so that the Young’s equation is applicable.
Fig. 4.

Schematic drawing showing the equilibrium shape of liquid sessile drop on a solid substrate (redrawn after Saiz et al.57)). (a) The liquid does not dissolve the substrate (b) The liquid partially dissolves the substrate.

According to the MD results given above, it is known that the prism plane can arrest oxygen atoms more easily than the basal plane. This leads us to investigate the influence of anti-spheroidizing elements on the prism plane. In the wetting simulation, 1380 oxygen (13.52 atoms/nm2) atoms were placed above the graphite prism plane with the proportion of 1:1 between oxygen and carbon atoms, as shown in Fig. 5. The simulation box, containing 8633 Fe atoms and 15720 carbon atoms, has the dimension of 10.4 nm*10.0 nm*3.9 nm. During the MD simulation, the graphite and oxygen atoms were fixed, the time step was set to be 1 fs, and the temperature of iron was set to 2000 K under the NVT ensemble.

Fig. 5.

Initial structure for the Fe cluster on the graphite substrate covered with oxygen atoms. (Online version in color.)

Figure 6 shows six sequential snapshots for the first 50 ps simulation in front and top views. These snapshots show the wetting process of molten iron on the substrate: first, the iron cluster starts with the contraction along the x-y plane and the protrusion in z direction; then, the contraction along the x-y plane continues, creating the iron cluster with a half-moon shape within 20 ps; after 20 ps, the iron cluster continuously spreads out along the x-y plane until an equilibrium is reached, and the spreading ability is a representation of wettability. The cross-sectional views of velocity in z direction after the relaxation is presented in Fig. 7, showing that the top atoms of iron cluster have higher velocity than the bottom atoms. This suggests that the wetting behavior of iron droplet should be dominantly influenced by these top atoms at the very beginning.

Fig. 6.

Snapshots of front views and top views for the first 50 ps, 13.52 oxygen atoms/nm2. (Online version in color.)

Fig. 7.

Cross-sectional views in x-z direction, the distribution of velocity in z direction after the relaxation. (Online version in color.)

In order to explore the influence of covering density of oxygen atoms on iron/graphite interfacial energy, three systems were built with different number (690, 345 and 0) of oxygen atoms covering on the graphite prism plane, as shown in Fig. 8. Thus, in these wetting simulations the covering densities of oxygen atoms are 6.76 atoms/nm2, 3.38 atoms/nm2 and 0 atom/nm2 respectively. The potential energy between iron droplets and substrates (ED–S) is defined as the sum of Fe–C interaction energy (EFe–C) and Fe–O interaction energy (EFe–O), such as ED–S=EFe–C+EFe–O, the EFe–C and EFe–O is given by Eq. (1). Figure 9 illustrates ED–S as a function of time, showing that the four simulated systems reached equilibrium within 60 ps. From Fig. 9, one can see that oxygen atoms strengthen the interactions between iron droplets and substrates, resulting in a decrease in the height of iron droplet (shown in Fig. 10). For instance, the height of iron droplet is about 4.5 nm without oxygen atoms and it decreases to 4.3 nm, 4.0 nm, and 3.6 nm when 345, 690, and 1380 oxygen atoms were added respectively. The decrease in the height of iron droplet causes an increase in the contact area (characterized with droplet base radius rB) between droplet and substrate (shown in Fig. 11), which reflects an increase in the spreading ability of iron droplet. These results indicate that oxygen atoms increase the wettability of iron cluster on the graphite substrate.

Fig. 8.

Simulation models with (a) 690 O atoms (b) 345 O atoms (c) no O atom on the graphite substrates. (Online version in color.)

Fig. 9.

Interaction energy between droplets and substrates, as a function of time.

Fig. 10.

Height of droplets as a function of time.

Fig. 11.

rB of equilibrium states in the simulations, obtained by 8 measurements.

Meanwhile, wettability is correlated with the contact angle (CA) between iron droplet and graphite substrate. The macroscopic CA θ is related to the microscopic CA through the modified Young’s equation.58) The relationship of CA θ, line energy τ, droplet base radius rB and interfacial energies of relevant phases can be described as:   

γ SV = γ SL + γ LV cosθ+ τ r B (4)
Young’s equation is recovered for macroscopic droplets, for the condition of 1/rB→∞. The macroscopic contact angle is defined as cosθ=(γSVγSL)/γLV, and the equation above can be rewritten as:   
cosθ=cos θ - τ γ LV 1 r B (5)

The cosθ is, thus, linearly related to the droplet base curvature 1/rB. Therefore, it is feasible to calculate the macroscopic CA through MD method. Figure 12 shows the calculated CAs of four simulated systems from MD simulations. From Fig. 12, it is seen that the CA is about 112.86° (±2.9°) without oxygen atoms and it decreases to 102.76° (±1.43°), 91.9° (±1.41°) and 78.65° (±1.73°) with 345 oxygen atoms, 690 oxygen atoms and 1380 oxygen atoms respectively placed on the substrate surface, suggesting that the oxygen atoms adsorbed on the graphite surface increase the wettability of iron droplet on graphite prism plane, in consistence with the results given above.

Fig. 12.

CAs of iron droplets on graphite prism planes covered with different number of oxygen atoms, obtained by 8 measurements.

According to the Eq. (3), different from the macroscopic interfacial energies of relevant phases, the microscopic interfacial energies are not only related to the contact angle, but also related to line energy τ and droplet base radius rB. According to the linear relation between cosθ and 1/rB, cosθ=cosθ when 1/rB=0. From MD simulations, the values of cosθ were calculated under four different surface conditions (the covering densities of oxygen atoms are 13.52 atoms/nm2, 6.76 atoms/nm2, 3.38 atoms/nm2 and 0 atom/nm2 respectively), given in Fig. 13. The values of cosθ can be obtained by the interpolation at 1/rB=0, showing that cosθ decrease as the number of oxygen atoms increases. As a result, the interfacial energy decreases with the increase of oxygen atoms adsorbed on the graphite prism plane, in agreement with other studies.44,59)

Fig. 13.

cosθ as a function of 1/rB with linear fit (dash lines).

4. Conclusions

In this work, molecular dynamics simulations were used to investigate the influence of oxygen atom on the iron/graphite interfacial behaviors. Based on these MD results, the following conclusions can be drawn:

(1) In comparison with the graphite basal plane, oxygen atoms were more favorably to be arrested by the prism plane, due to its different surface condition, in consistence with experiment. This suggests the reliability of our MD simulations. From the MD simulations, we discover two dominant locations (type I and type II) of oxygen atoms adsorbed on the prism plane.

(2) The wetting process of molten iron on the substrate was illustrated by the MD simulations, showing that the wetting behavior of iron droplet should be dominantly influenced by the top Fe atoms in the iron cluster.

(3) The wetting simulations under different surface conditions reveal that oxygen atoms adsorbed on the graphite substrate can strengthen the interactions between iron droplet and substrate (owing to strong Fe–O interactions), leading to an increase in the wettability of iron droplet, which are characterized with a decrease in the height of iron droplet, an increase in the contact area, and a decrease in the contact angle.

(4) Thus, based on these MD results, it becomes clear about the fundamental physics behind the influence of oxygen atoms on the iron/graphite interfacial behaviors: oxygen atoms adsorbed on graphite substrate strengthen the interactions between iron liquid and graphite substrate, reducing the dynamics of the bottom Fe atoms in iron cluster on the one hand and facilitating the spreading (or lateral) motion of the top Fe atoms in iron cluster toward the graphite substrate on the other hand.


This research was financially supported by a project supported by the National Nature Science Fund Projects of China (No. 51305149), China Postdoctoral Science Foundation (2014M552035), and Special Funding of China Postdoctoral Science Foundation (2015T80795), and Natural Science Foundation of Guizhou Province (QKH-JZ[2014]2003).

© 2018 by The Iron and Steel Institute of Japan