Nonequilibrium Molecular Dynamics Simulations of Coal Ash

Both molecular dynamics (MD) and nonequilibrium molecular dynamics (NEMD) simulations were performed to simulate coal ashes using the Guillot-Sator model in this work. The structural and transport properties of coal ashes at high temperatures have been obtained. Superheating of coal ash system with anorthite crystal structure initial configuration has been observed for MD simulation which explains the discrepancy between previous MD simulation results and FactSage thermochemical calculations. The fluxing effects of both calcium oxide and sodium oxide have been investigated systematically through MD and NEMD simulations. Moreover, the viscosities of coal ash systems have been computed by two methods: (1) Stokes-Einstein equation; (2) NEMD simulations. Estimations of viscosities for various coal ash systems based on Stokes-Einstein equation exhibit a strong temperature dependence of viscosity, which agrees with previous experimental results. On the other hand, NEMD simulation results that showed a strong shear-thinning feature, failed to reproduce this strong temperature dependence of viscosity, possibly due to the short simulation time. Nevertheless, NEMD simulations not only provide us detailed information about atoms dynamics under shear, but also allow us to model the coal ash system far from equilibrium which cannot be accessed by thermodynamics calculation using software like FactSage.


Introduction
The mineral formation (e.g., reaction among silica melt, CaO and Al 2 O 3 ) during coal combustion took place due to heterogeneous nature of feeding coals [1][2][3][4][5][6][7][8]. These reactions could cause slagging and fouling problems in coal fired power plants which caused reduction in combustion efficiency [9][10][11][12][13]. Therefore, to facilitate the long-term stable operation of entrained flow gasifier, which is one of the core parts in IGCC (integrated gasification combined-cycle technology) with high power generation efficiency and low greenhouse gas emission, a good understanding of coal ash fusibility as well as flow properties at high temperature is required [14][15][16]. Extensive experimental studies, e.g., slag viscosity measurement [17,18], XRD [19], FTIR [20], SEM [20], and thermomechanical analysis [21], have recently been performed to investigate the compositions and fusibility of coal ashes. As the desired slag tapping temperatures for different gasifiers typically range from 1573 to 1673 K [17,18], the addition of fluxing agent for those coals with high ash melting temperature is usually necessary [22]. Hence, the effects of various fluxing agents, such as Na 2 O, K 2 O, MgO, TiO 2 , on lowering the melting temperatures of different coal ashes (e.g., highcalcium coal ash [23], brown coal ash [24]) have been carefully examined recently. Based on those experimental results, various empirical models [17,18,25] have been proposed to relate the viscosities and fusibilities of various coal ashes to their compositions at high temperatures. Moreover, with the development and wide applications of thermochemical databases and software such as FactSage [26], the predictions of coal ashes fusibilities as well as their structures according to their compositions become extremely convenient and accurate on many occasions.
However, since those empirical models and thermochemical software can only address thermodynamics of coal ashes in equilibrium, molecular dynamics simulation, which can track the motions of individual atoms, has a big advantage in studying complex phase behaviors of coal ash dynamically. B. Guillot et al. [27,28] developed a simple interionic potential (Guillot-Sator model) to describe a nine component system (K 2 O-Na 2 O-CaO-MgO-FeO-Fe 2 O 3 -Al 2 O 3 -TiO 2 -SiO 2 ), which successfully reproduced many structural, thermodynamic, and transport properties of a number of natural silicate melts. And this model was adopted by Dai et al. [19,29,30] to investigate the fusion properties of various coal ashes using molecular dynamics simulations. By combining with experiments and FactSage thermochemical calculations, Dai et al. [29] believed they were able to reproduce the melting temperatures of coal ashes and elucidate the effects of fluxing agents like calcium oxide [29] on coal ash fusion. Nevertheless, determination of liquid-glass transition temperature T g is not easy [27,31], and depends on the applied cooling rate as well as on the size of the simulation system [32]. On the other hand, the viscosity or even the rheological properties of coal ash at slag tapping temperature is of great interest, as it is critical for gasifier operation. Therefore, the application of Guillot-Sator model in combination with methods like the Green-Kubo integrand [33,34] or nonequilibrium molecular dynamics (NEMD) simulation could be more informative compared to conventional MD simulations. Although the simple silicate melts shear viscosities have been studied using molecular dynamics simulations decades ago [35], a systematic NEMD study of multiple components coal ash melts equipped with Guillot-Sator model can still shed some light on coal ash slagging behavior under gasification conditions. In addition, recent quantum chemistry calculations of aluminosilicate [36] might help improving the simple model like Guillot-Sator potentials, as the lacking of many-body interactions of that model obstructs its reproduction of crystal structure, which can be alleviated by combining quantum chemistry results.
Using the Guillot-Sator model [27], we seek to study systematically transport and structural properties of aluminosilicate with calcium oxide and sodium oxide at varying temperatures. Our nonequilibrium molecular dynamics simulation results which can bridge the gap between coal ash microstructure and its flow properties is not only of interest in its own right, but is a good starting point for developing a new atomistic model incorporating many-body interactions among different kinds of atoms. With the help of molecular simulations of coal ashes, the mineral formation at high temperatures during coal combustion/gasification which is difficult to access by experimental studies can be investigated numerically even when the simulation systems is nonequilibrium.

Molecular Model
As discussed above, the natural silicate melts model developed by B. Guillot et al. [27,28] was employed to simulate our coal ash systems, where the pairwise interaction energy (between the i-th and j-th atoms) is modeled by simple Born-Mayer-Huggins (BMH) function as follows.
where r ij is the distance between atoms i and j; q i is the effective charge associated with the ith atom; A ij and ρ ij are energy parameters for the repulsive forces between atoms i and j; and C ij is the energy parameter for dispersive forces. These potential coefficients for various types of atoms are shown in the Table 1 below. Note that the cation-cation interactions in this model are described only by electrostatics repulsion (the first term on the right hand side of Equation (1)) [27].

Simulation Method
In this work, all simulations were performed using the LAMMPS software package [37]. Visualization of atoms ( Figure 1) as well as calculation of various radial pair distribution functions were carried out by VMD [38]. And the initial configurations of all our simulations except where otherwise stated were prepared by randomly positioning different numbers of those atoms discussed above in the simulation box, according to their concentrations. For the molecular dynamics simulations, simulation system was first equilibrated by performing a 500 ps isothermal-isobaric (NPT) simulation, as the volume of the system was stabilized. And then, another 500 ps NPT simulation was performed to track the atoms motions in our simulation systems with a 0.25 femtosecond ( f s) timestep, the same as in the equilibration stage. The long-distance electrostatic interaction was computed by Ewald method [39] with a 8.5 Å cutoff distance. And the cutoff distance for all other interactions is 12 Å. The apparent diffusion coefficient D of atoms in molten coal ashes was determined using the following equation during molecular dynamics simulations: where MSD is the average mean square displacement of all atoms in the simulation box; r i (t) is the position of atom i at time t; and r i (0) is the initial position of atom i. The shear viscosities of coal ashes with different compositions were evaluated in the constant volume nonequilibrium molecular dynamics (NEMD) simulations with applying shear flows to our systems. And the shear flows were generated by combining SLLOD method [40,41] and deforming simulation box, as shown in Figure 1. Due to the external perturbation caused by shear flow, the timestep for all our NEMD simulations was set to 0.05 f s. In addition, to minimize the effects of the structural variation of coal ash as a result of shear stress, 100 ps NEMD simulations were performed to achieve steady state before the final 250 ps NEMD simulations. And the shear rate for all our NEMD simulations except where otherwise stated was set to 3 ns −1 .

Initial Configuration Effects
Although many methods have been developed for the computation of melting point of a compond, predicting melting points by molecular dynamics simulations is still a challenging problem [31]. Due to the existence of superheating during coal combustion, or hysteresis phenomenon, calculation of melting points by simply heating a perfect lattice until the lattice breaks down when the system is assumed to be at melting temperature might result in a significant overestimation of melting point [7,42]. Hence, the so-called "direct" methods involving the direct simulation of melting process, such as hysteresis method [43], the voids method [44], and solid-liquid interface-based methods [45], have been developed to avoid this overestimation. Since the strict thermodynamics definition of melting point is the temperature at which the free energy of crystal phase is equivalent to the free energy of liquid phase, a bunch of methods involving free energy calculation, e.g., Hoover and Ree's single-occupancy cell method [46], Frenkel and Ladd's Einstein crystal method [47], and λ-integration method developed by Grochola [48] and extended by Maginn et al. [49], have also been developed which can predict the melting temperature even more accurately as those methods allow researchers to minimize the effects of hysteresis phenomenon. And this hysteresis phenomenon explaines the overestimation of melting point of anorthite by MD simulations against FactSage thermochemical calculations which has been observed by Dai et al. [29], as demonstrated below.
Because the Guillot-Sator model [27,28] was developed for the simulation of natural silicate melts, it is hardly possible to reproduce crystal structures of natural silicate by this specific model. However, this Guillot-Sator model still allows us to evaluate the superheating of a simulation system, as the hysteresis phenomenon exists even for simple monatomic molecules [42]. Hence, before shifting our focus onto molten aluminosilicates, we seek to elucidate the initial configuration effects on molecular dynamics simulations.
A comparison of the isothermal-isobaric MD simulations results with two different initial simulation configurations is shown in Figure 2a,b. One initial configuration of atoms in simulation box is the crystal structure of anorthite under high pressure [50] as shown in Figure 2a, which has also been explored by Dai et al. [29] with MD simulations recently. For the other initial configuration, all atoms in the simulation box were randomly positioned. Note that the half blue-half white calcium atoms in Figure 2a indicate that the occupation of calcium atoms is only 0.5. And an example of the simulation box of 3 × 3 × 3 periodic supercell is shown in Figure 2b. A 500 ps NPT MD simulation was carried out for both the anorthite crystal structure initial configuration and the randomly positioned atoms initial configuration. And the densities of this CaO·Al 2 O 3 ·2SiO 2 system as functions of time with these two initial configurations are plotted in Figure 2c. Then the second NPT MD simulation was performed in which the average mean square displacement of all atoms as a function of time was evaluated for both initial configurations. The MSDs are plotted against time in Figure 2d. As shown in Figure 2d, the average MSD of all atoms with randomly positioned atoms as the initial configuration increases rapidly with time, since the MSD after 200 ps simulation reaches 100 Å 2 which indicates that the simulation system is in liquid phase. On the contrary, the average MSD of all atoms with anorthite crystal structure as the initial configuration remains a small number around 1 Å 2 which is about 1 atom size after 500 ps MD simulation. This small MSD suggests that the CaO·Al 2 O 3 ·2SiO 2 system with crystal structure initial configuration remains a crystal after 1 ns MD simulation (two simulations combined), which agrees with the 2050 K melting point obtained by Dai et al. [29] with MD simulations. However, Dai et al. [29] obtained the ≈ 1800 K melting point result using FactSage software which can be considered as the accurate number. The gap between these two methods can be bridged by the superheating phenomenon observed in our simulations, as the atoms have been confined around their equilibrium locations in the crystal structure. The thermal fluctuation itself is not strong enough as a perturbation to drive the system away from crystal structure which is a metastable state in the free energy landscape. This can be confirmed by Figure 2c, as the density of CaO·Al 2 O 3 ·2SiO 2 system with crystal structure decreases instantly once the MD simulation begins, and then remains a constant for the rest of the simulation (red curve). This is because the average interatomic distance within anorthite crystal structure under high pressure will increase instantly after applying low pressure (1 atm) during NPT simulation, and will not change any further since all atoms are only vibrating around their equilibrium positions (MSD≈ 1 Å 2 ). For the initial configuration with randomly positioned atoms in a large simulation box to avoid the overlapping, the density increases gradually and then reaches a plateau (blue curve). The density of crystal structure is slightly below the density of anorthite melt (Figure 2c), for the liquid phase can be more compact in this condition.

Effects of Calcium Oxide on Aluminosilicate Melts
The MSD as a function of time is plotted in Figure 3a at different temperatures for anorthite system with randomly positioned atoms as its initial configuration. As shown in Figure 3a, the MSD increases almost linearly with time in the log-log plot, which indicates that the motions of atoms are just normal diffusion. Note that the MSD as a funtion of time becomes rough and deviates from staight line at low temperatures, due to the poor sampling as the travelling distances of atoms decreases dramatically at low temperatures. The average apparent diffusion coefficient as a funtion of inverse of temperature which can be calculated by Equation (2) is plotted in Figure 2b. And the relationship between diffusion coefficient D and the inverse of temperature obeys the Arrhenius equation, which agrees with results of B. Guillot et al. [27]. In addition, it has been observed that the volume of our simulation system increases linearly with temperature for isothermal-isobaric simulations, as shown in Figure 3c. This observation agrees with previous studies [27,28], and implies that no phase transition occurs even when temperatures decreases below 1400 K. Nevertheless, this phenomenon suggests a supercooling rather than a low melting temperature (<1400 K), for the crystal nucleation is a rare event and extremely unlikely to be observed during a molecular dynamics simulation. Although it is not easy to compute the melting points of aluminosilicate, how the calcium oxide affects the phase transition of aluminosilicate can be investigated by MD simulation in an indirect manner. The diffusion coefficients D of atoms within CaO-Al 2 O 3 -SiO 2 systems with various CaO weight percentages are plotted against the inverse of temperature 1000/T in Figure 4. And the ratio between the number of Al atoms and the number of Si atoms is fixed to be 1. Again, the relationship between the diffusion coefficient D and inverse of temperature 1000/T obeys the Arrhenius equation for CaO-Al 2 O 3 -SiO 2 with different CaO contents. Besides, the diffusion coefficients increases as CaO weight percentage increases from 5 wt% to 40 wt%, which implies that the melting temperature of CaO-Al 2 O 3 -SiO 2 decreases continuously as CaO content increases, which is different from FactSage predictions [29] that there is a local peak of CaO-Al 2 O 3 -SiO 2 melting temperature around CaO weight percentage =20%. This discrepancy of melting temperature between MD simulation predictions and FactSage thermochemical calculations can be attributed to two aspects: (1) the supercooling of MD simulation systems; (2) the lack of many-body interactions among different atoms of Guillot-Sator model [27,28]. For simplicity, B. Guillot et al. only considered the pair interactions which is described by Equation (1). Therefore, it is unlikely to reproduce every crystal structures by this simple model. As a result, the FactSage calculation which considers the phase transitions between different crystal structures based on the chemical compositions of the modeling system can provide much more accurate predictions on melting point. To fully reproduce the melting process of aluminosilicate with different fluxing agents like calcium oxide, or sodium oxide, people should not only use a more accurate as well as sophisticated melting point calculation method such as λ-integration method [48], but also should develop more accurate molecular models incorporating many-body interactions. Nevertheless, simulations results in this work demonstrate that even this simple Guillot-Sator model can predict general trend of the relationship between melting points and calcium oxide contents.
The microstructure of CaO-Al 2 O 3 -SiO 2 system can be evaluated through the radial pair distribution function g(r) between atoms of type A and B, which is defined by the following equation: where ρ B local is the averaged atom density of type B around atoms A; N A and N B are total number of atom A and B, respectively; and r ij is the distance between the ith atom (type A) and jth atom (type B).
The radial pair distribution function (RDF) g(r) curves of different atomic pairs for the CaO-Al 2 O 3 -SiO 2 system with 15 wt% CaO at various temperatures are plotted in Figure 5. Since all RDF curves are very smooth, the simulation systems have been fully equilibrated and properly sampled in these conditions. As shown in Figure 5, the equilibrium interatomic distances, i.e., the first peak locations in Figure 5 [29]. Besides, the RDF curves stretch out slightly as the temperature increases from 1800 K to 2200 K. This explains the previous results that system volume increases as temperature increases. In addition, since the average separation distance between atoms increases at high temperatures, atoms become more mobile, thus, resulting in a larger diffusion coefficient, as demonstrated in Figures 3b and 4. How the CaO content affects the radial pair distribution funtions g(r) of Si-O and Al-O pairs is plotted in Figure 6. And the change of RDF curves caused by varying CaO contents is negligible. By comparing results shown in Figures 4 and 6, it can be concluded that that the main reason of the increase of average diffusion coefficient of all atoms in CaO-Al 2 O 3 -SiO 2 system as CaO content increases is simply due to the more mobile Ca atoms, for the microstructure consist of other atoms remains almost unchanged. However, the adding of Ca atoms might cause the phase transition from mullite to anorthite which can hardly been observed by this simple model. Therefore, the simple pair interaction described by Equation (1) is not sufficient to elucidate the complex phase behaviors of aluminosilicates.

Effects of Sodium Oxide on Aluminosilicate Melts
How the addition of sodium oxide affects the aluminosilicate melts has also been investigated, since a lot of brown coals contain large amount of sodium element [24]. For all the simulations of Na 2 O-CaO-Al 2 O 3 -SiO 2 systems, the ratios among Ca, Al, Si atom numbers are fixed to be 1:2:2. As the Na 2 O content increases from 0 wt% to 10 wt%, the average diffusion coefficient of all atoms in Na 2 O-CaO-Al 2 O 3 -SiO 2 system increases continuously, as demonstrated in Figure 7. Therefore, the Na 2 O works as a fluxing agent just like CaO does. The relationship between diffusion coefficient D and 1000/T still follows the Arrhenius equation. Note that the slight deviation of D as a function of 1000/T from straight line with 4 wt% & 6 wt% of Na 2 O might be a result of poor sampling instead of inherent properties of the simulation systems.  Figure 8, as the temperature increases from 1500 K to 2000 K. Similarly, the increase of temperature enlarges the average separation distances between different atoms, for the thermal fluctuations are stronger at high temperatures. The increase of kinetic energy of atoms at higher temperature can drive those atoms further away from their equilibrium locations, which causes a larger separation between cation atoms and oxygen atoms. As the atoms are moving far away from their neighboring atoms, there are fewer neighboring atoms obstructing their diffusions, which explains the increase of diffusion coefficient with temperature shown in Figure 7 [24].
The radial pair distribution functions of different atom pairs with changing Na 2 O concentration at 2000 K are plotted in Figure 9. Like Ca atoms, the effects of addition of sodium atom on g(r) can be negligible, which means the increase of D with Na 2 O content is again due to the higher mobility of sodium atoms. Still, the diffusion coefficient D increases due to the increase of Na 2 O weight percentage. This could suggest that the addition of sodium oxide decreases the melting temperature of aluminosilicate which agrees with the experimental observations [24]. Although some previous studies reported that the change of RDFs implies phase transition [29], similar phenomena have not been observed in this study as all the simulation systems remain liquid phase even for very low temperature agreeing with B. Guillot's results [27] except the MD simulation with crystal structure initial configuration.

Non-Equilibrium Molecular Dynamics Simulations
Viscosities of molten slags composed mainly of silicates must be evaluated for gassifier operations. And silicates viscosities exhibit a strong temperature denpendence [51], as the viscosities of glass can decrease by 10 orders of magnitude when temperature increases from 773 to 1273 K [52]. There are multiple methods for estimating shear viscosity by molecular dynamics simulations [35]. One of them is the Green-Kubo treatment [33] which is to integrate the stress autocorrelation function. Theoretically, the Green-Kubo treatment allows people to access all the information about the viscoelasticity of a system. However, the obtaining of a smooth stress autocorrelation function is usually very computational expensive. Thus, here in this work, the nonequilibrium molecular dynamics (NEMD) simulation is adopted to compute the viscosities of both CaO-Al 2 O 3 -SiO 2 and Na 2 O-CaO-Al 2 O 3 -SiO 2 systems. The steady-state shear flow in non-equilibrium molecular dynamics simulations can be introduced by applying the classical Lees-Edwards boundary conditions [53]. For our NEMD simulations, the shear flow was generated by a combination of deforming simulation box and applying SLLOD algorithms [41]. Because the size of our simulation box is extremely small compared to the experiments, the shear rates applied to the simulation systems are very large, as shown in Figures 10-12.
Shear viscosity as a function of shear rate for Na 2 O-CaO-Al 2 O 3 -SiO 2 system with 2 wt% Na 2 O at 2000 K is plotted in Figure 10. The viscosity is given by following equation: whereγ is the shear rate; P xy is the shear stress. A strong shear-thinning feature of aluminosilicate melts have been observed as demonstrated in Figure 10, as the simulated viscosity of Na 2 O-CaO-Al 2 O 3 -SiO 2 system decreases by 50-fold when the shear rateγ increases from 1 ns −1 to 1 ps −1 . This strong shear-thinning phenomenon originates from the strong perturbation introduced by shear flow. Since the shear rateγ is larger than the experimental values by several orders of magnitude, the oxygen bridges which help the formation of long distance correlated structures break down instantly after introducing shear flow. In other words, the aluminosilicate melts will behave more like a simple Lennard-Jones liquid under high shear.  The effects of CaO contents as well as temperature on the viscosity of CaO-Al 2 O 3 -SiO 2 system is demonstrated in Figure 11. For all different temperatures, the viscosity of aluminosilicate melts decreases by about 40% as CaO weight percentage increases from 5% to 40%. This result agrees with our previous results on diffusion coefficients, as shown in Figure 4. The reduction of system viscosity as a result of addition of CaO confirms that CaO can work as fluxing agent in aluminosilicate system. At 1250 K (orange curve in Figure 11), viscosity first decreases when CaO content increases from 5 wt% to 15 wt%, and then it increases with CaO content ranging from 15 wt% to 20 wt%. Finally, the viscosity decreases as CaO content increases, again. These results at 1250 K seems to suggest that the melting temperature of aluminosilicate follows the similar trend of the relationship between viscosity and CaO content, which happens to be the same as FactSage predictions [29]. However, this agreement between orange curve in Figure 11 and the FactSage calculation resutls from previous research [29] is merely a coincidence, as this agreement has not been observed by comparing diffusion coefficients of aluminosilicates with various CaO percentages demonstrated in Figure 4. Moreover, the strong temperature dependence of viscosity has not been observed through our NEMD simulaitons, although the simulated viscosity at 1500 K is of the same order as FactSage predictions [30]. This failure of predicting temperature dependence of viscosity is attributed to the limited size of simulation box and, more importantly, the very limited simulation time [35]. Besides, the extremely high shear rateγ leaves atoms no time to respond to external perturbations, which also results in an underestimation of viscosity especially for low temperatures. Although for nonequilibrium molecular dynamics simulation a lot of issues might arise due to the lack of sufficient sampling, NEMD sitll has its own unique advantages, for it can provide detailed structural information under shear directly which is difficult to access through Green-Kubo treatment.
The simulated shear viscosity computed by NEMD as a function of Na 2 O content for the anorthite systems with sodium oxide is plotted in Figure 12. The viscosity of Na 2 O-CaO-Al 2 O 3 -SiO 2 system decreases almost linearly by about 40% as sodium oxide weight percentage increases from 0 wt% to 10 wt%, which confirms the fluxing effect of sodium oxide. NEMD simulations results for 4 different temperatures are shown in Figure 12. They suggest that viscosity decreases as temperature increases, which agrees with experimental observations. However, the small difference between viscosities for 2000 K and that for 1250 K (∼20%) does not show the huge effects of temperature on viscosity. Again, this big discrepancy between NEMD predictions and experimental results [51] can be attributed to the poor sampling and short simulation time. For the same reason, NEMD simulations cannot be expected to give valid structual information and show how those different atoms respond to this fast shear correctly.
Note that the viscosity of silicate melts can also be estimated by Stokes-Einstein equation [35] as follows.
where a is the radius of diffusing spheres; and D is the diffusion coefficient. For simplicity, a is assumed to be 2 Å, which can be approximated by the radial pair distribution functions g(r) shown in Figures 8 and 9. Due to the strong attractive force between cations and oxygen anions, the diffusions of cations and oxygen anions are not uncorrelated. However, for simplicity, all atoms in the simulation system are assumed to be independent diffusing spheres. As shown in Figures 4 and 7, the diffusion coefficients of atoms for various aluminosilicate melts are provided. Let's take the CaO-Al 2 O 3 -SiO 2 system with 5 wt% CaO at 2000 K as an example. The average diffusion coefficient of all atoms within this specific system is D =4.25×10 −2 Å 2 /ps. Therefore, based on Equation (5), the viscosity of this system can be obtained which is 0.0172 Pa · s. Calculation of viscosity based on Stokes-Einstein equaion is merely a rough approximation, due to the correlation between cation diffusion and anion diffusion. Nevertheless, this method demonstrates the strong temperature dependence of silicate melts. For example, in Figure 7   From a physics point of view, as the aluminosilicate temperature approaches the melting point, more and more electrons are leaving their ground states, which will alter the wavefunctions of the system. As a result, an empirical force field parameterized at one temperature is supposed to fail for another temperature. Moreover, the simple model like Guillot-Sator model [27] without many-body interactions can never reproduce the various crystal structures observed in experiments and tabulated in FactSage. Nevertheless, software such as FactSage which is based on thermodyamics calculation of the system cannot give accurate predictions of systems far away from equilibrium. The coal ash/slag at high temperature, due to the complex phase behaviors as well as chemical reactions, it is not surprising that coal ash system is far away from thermodynamic equilibrium. Therefore, methods like NEMD simulation which allow people to track the dynamics of such non-equilibrium systems have their own unique advantages. As shown in this work, computation of viscosities using NEMD is still very challenging for coal ash systems. But it is a good starting point for further development of coal ash models which can predict the complicated phase behaviors at different temperatures. Hopefully, with the help of faster computational facilities, nonequilibrium molecular dynamics simulations of large systems comparable to experiments with sufficient long simulation time would become possible.

Conclusions
Using the Guillot-Sator model [27], the relationship between the diffusion coefficients of coal ash systems and calcium oxide/sodium oxide weight percentages can be obtained. It has been found that the diffusion coefficient increases continuously as fluxing agent (CaO/Na 2 O) content increases, which implies that the melting point of coal ash decreases continuously with fluxing agent content. And this result disagrees with FactSage predictions [29] which suggests that many-body interactions as well as much longer NEMD simulaitions are necessary to predict the complicated phase behaviors of nonequilibrium coal ash systems in the gasifiers. A potential solution to this problem is ReaxFF (reactive force field) MD simulation [54] which has been widely used in simulation of coal combustion/gasification [55]. Since this ReaxFF which allows chemical reactions to take place during MD simulations is originally parameterized based on DFT calculations of many hydrocarbons chemical reactions [54], it can be parameterized for coal ashes systems in the same manner, as many-body interaction is already included. A combination of reactive force field for organic molecules [54,55] and reactive force field for mineral formation reactions is a complete model for coal combustion/gasification. However, for such complete model, initial configuration effects and superheating still exist during molecular dynamics simulations, due to thermal fluctuation is not strong enough to drive the simulation system out of metasatble state, as discussed above. Therefore, the simulation results in this work can provide some clues for further force field development. And such molecular dynamics simulations should be cross-check with chemistry and mineralogical compositions of slags in the gaisifiers to minimize the inaccuracy of reactive force field.