Molecular Dynamics Simulation of Viscosity of the CaO, MgO and Al₂O₃ Melts

Abstract

Viscosity of oxide melts is a fundamental physicochemical property that plays a critical role in important technological and natural processes like slag flow, slag/metal separation, volcano eruptions etc. Unfortunately, many oxides melt at extremely high temperatures and are highly corrosive in the liquid state, which makes experimental measurement of viscosity by conventional experimental techniques a difficult, expensive, and sometimes impossible task. In this case it might be helpful to use other methods to determine viscosity, such as modelling or simulation. In the present paper the shear viscosity coefficients of the liquid CaO, MgO, and Al₂O₃ have been simulated via the classical molecular dynamics and compared to the available viscosity data (e.g. experimental viscosities, other model predictions) collected and stored in a databank by one of the authors. Different viscosity calculation techniques (e.g. the Green-Kubo equation, the Einstein relation, and the non-equilibrium molecular dynamics methods) coupled with MD have been employed to ensure a more reliable viscosity calculation. It has been shown that the simulated viscosities of the CaO and MgO melts are close to those calculated by phenomenological viscosity models and represent a plausible estimation for viscosity of these unary systems. It has also been shown that the simulated viscosity of the Al₂O₃ melt is lower than the majority of the available experimental data. In general, it has been demonstrated that the molecular dynamics simulation could provide a reasonable estimation of viscosity, especially if no other viscosity data is available.

1. Introduction

Many pure oxides are characterized by extreme melting temperatures and/or extreme reactivity in the molten state, which precludes experimental investigation of important melt properties such as viscosity, density, etc. Among them the calcium and magnesium oxides are encountered in many industrial slags, melts and refractories, and the aluminum oxide is one of the most common ceramic and glass components. The effects of these oxides on slag viscosity are different: while CaO and MgO are usually considered to be network modifiers, which tend to decrease slag viscosity, Al₂O₃ can act both as a modifier (thereby decreasing viscosity) and as a network former (increasing viscosity). A number of phenomenological models have been developed up to date to calculate viscosity of multicomponent liquid slags (including CaO, MgO and Al₂O₃) as a function of temperature and chemical composition. In the present study four typical viscosity models are selected for comparison with MD results for the CaO and MgO melts: the FactSage,¹⁾ GTT,²⁾ Kondratiev-Jak^3,4) and Starodub⁵⁾ models. It should be mentioned that in order to develop such a viscosity model for a multicomponent oxide system first it is necessary to evaluate viscosities of the pure oxide components, which are used then to model viscosities of binary, ternary and high-order slags. If no experimental data is present for viscosity of a pure oxide, arbitrary viscosity parameters are typically selected for modelling on the basis of extrapolation from corresponding binary or ternary oxide systems. However, such estimation is not always reliable, as the experimental data in binary or ternary systems can be obtained for compositional ranges that are far from the pure oxide. Therefore, there is a need to provide a more reliable estimation for viscosity of those pure oxide melts, for which no experimental data is available. Among other methods, classical molecular dynamics (MD) simulation appears to be a reliable and powerful tool for evaluation of physical properties of molten pure oxides including viscosity.

A number of different techniques have been developed to evaluate the shear viscosity coefficient from MD simulation, which include the equilibrium MD methods (Einstein, Green-Kubo) and non-equilibrium MD methods (periodic perturbation, SLLOD and Müller-Plathe). Compared to phenomenological viscosity models the approach based on MD simulation can provide more reliable and physically rigorous viscosity values.

In this study five different MD methods mentioned above were used to calculate the shear viscosity of selected pure oxide melts (CaO, MgO, and Al₂O₃), and the Born-Mayer-Huggins potential was used to describe interatomic interactions for MD simulations. As no experimental data exists for liquid CaO, the shear viscosity of the pure CaO melt calculated in the present work for the first time is compared to viscosities predicted by the different models. Viscosities of the MgO and Al₂O₃ melts calculated by MD simulation are also compared to the available calculations and experimental data, respectively.

2. Computational Methods

2.1. Molecular Dynamics Simulation

The LAMMPS non-commercial software package⁶⁾ was employed to carry out classical MD simulations with the Born–Mayer–Huggins (BMH) potential, for which the potential energy of interactions between atoms is given by the following equation:

  

$$u\left(r_{ij}\right)=\frac{q_i{q}_j}{r_{ij}}+A_{ij}exp\left(\frac{{\sigma }_{ij}-r_{ij}}{{\rho }_{ij}}\right)-\frac{C_{ij}}{r_{ij}^6}+\frac{D_{ij}}{r_{ij}^8}$$(1)

with i and j being atomic species and r_ij being the interatomic distance between two ions of species i and j. It should be noted that only the van der Waals term was taken into account in the dipolar expansion. Parameters of Eq. (1) with the corresponding references are given in Table 1. Parameters by Alvares et al.⁷⁾ were used for simulation of the CaO melt, while parameters by Arkhipin et al.⁸⁾ were used for simulation of the MgO and Al₂O₃ melts. Parameters from Gutierrez et al.⁹⁾ are given for comparison.

Table 1. Parameters of the Born-Mayer-Huggins potential.

Atom/interaction	q	Aij (eV)	ρij (Å)	σij (Å)	Cij (eV Å6)	Reference
Ca	1.2					7)
Mg	1.2
Al	1.8
O	−1.2
Ca–Ca		0.0035	0.0800	2.3440	20.9856
Ca–O		0.0077	0.1780	2.9935	42.2556
O–O		0.0120	0.2630	3.6430	85.0840
Mg–Mg		38868.5	0.2066	0.0	43.773	8)
Mg–O		12655.4	0.2089	0.0	36.798
Al–Al		0.00290	0.0680	1.5704	14.0498
Al–O		0.0075	0.1689	2.6067	20.7448
O–O		0.01015	0.2140	3.2678	5.0840
Al	1.4175					9)
O	−0.9450
Al−Al		0.00295	0.068	1.5704	14.0524
Al−O		0.00746	0.172	2.6067	34.5814
O–O		0.0120	0.276	3.643	85.1011

Using this empirical potential, classical molecular dynamic simulations were performed with the LAMMPS code. The initial crystalline configuration used for MD simulations for all oxides was a cubic state in the Green-Kubo and Einstein methods and an orthorhombic state (L×L×3L) in the Müller-Plathe method. The initial crystalline configurations were subjected to melting at a temperature of 100 K above the highest temperature used for viscosity calculations, which was carried out in the isobaric–isothermal ensemble (NPT) using a Nosé–Hoover thermostat coupled with an isotropic barostat¹⁰⁾ for 100 ps. The melted configuration was initially equilibrated at a given temperature in the NPT ensemble to determine the mean volume of the system over a period of 100 ps. Subsequently, the system volume was adjusted to the mean value and the system was then equilibrated in the canonical (NVT) ensemble with the coupling constant of 0.1 ps over 100 ps. The initial configuration used in the SLLOD algorithm was constructed by placing molecules randomly within a cubic box. Four distinct configurations were employed in the SLLOD calculations of viscosity in total. All configurations were equilibrated at a given temperature in the NPT or/and NVT ensembles, in accordance with the aforementioned procedure for alternative methods. The MD simulations necessary for the periodic perturbation method were performed only for CaO. The time required for each stage of the simulation depends upon the dimensions of the system in question. The melting stage continued 50, 100 and 200 ps for the small, moderate and large systems, respectively. The equilibration stage in the NPT and NVT ensembles lasted 80, 150 and 300 ps for the small, moderate and large systems, respectively. Periodic boundary conditions were applied in the three directions of space, and the long-range Coulomb interactions were treated with the Ewald summation method. Equations of motion were solved numerically by the Verlet algorithm in the velocity form using the time step of 1 fs. The number of atoms contained within the simulation box was dependent upon the MD method employed for the calculation of viscosity. The cubic box containing 729 Ca or Mg atoms and 729 O atoms for CaO or MgO and 600 Al atoms and 900 O atoms for Al₂O₃ was used for the Green-Kubo, Einstein and SLLOD viscosity calculations. The orthorhombic box (L×L×3L) with 1029 Ca or Mg atoms and 1029 O atoms for CaO or MgO and 768 Al atoms and 1152 O atoms for Al₂O₃ was used for the Müller-Plathe method. The cubic boxes containing 216, 1728 and 5832 Ca and O atoms, respectively, were used for the periodic perturbation method. The temperature ranges were 3200–4000 K, 3100–4000 K, and 2300–3200 K for CaO, MgO and Al₂O₃, respectively. The duration of the production stage is dependent upon a MD method employed for the calculation of viscosity, temperature and oxide. A comprehensive examination of the MD production stage and associated viscosity calculations is provided in Section 3.

2.2. Viscosity Calculation

The shear viscosity can be obtained from an equilibrium MD simulation via the stress tensor autocorrelation function by the Green-Kubo formula:¹¹⁾

  

$$\eta =\frac{1}{3}\frac{V}{k_{B}T}{\sum }_{\alpha ,\beta }{\int }_0^{\infty }⟨P_{\alpha \beta }\left(t\right)P_{\alpha \beta }\left(0\right)⟩dt$$(2)

where V is the system volume, k_B is Boltzmann’s constant, T is the absolute temperature, t is the time, P_αβ is the off-diagonal pressure tensor component, α, β refer to x, y, z and α ≠ β.

Following Nevins and Spera,¹²⁾ two independent stress tensor components (P_xx – P_yy and P_yy – P_zz) were included to estimate the shear viscosity. The standard deviation of the calculated shear viscosity ‘components’ represents the error of this estimation:

  

$$\begin{array}{l}\mathrm{\Delta }\eta =\\ \sqrt{{\displaystyle \frac{{\left({\eta }_{xy}-{\eta }_{av}\right)}^2+{\left({\eta }_{xz}-{\eta }_{av}\right)}^2+{\left({\eta }_{yz}-{\eta }_{av}\right)}^2+{\left({\eta }_{xxyy}-{\eta }_{av}\right)}^2+{\left({\eta }_{yyzz}-{\eta }_{av}\right)}^2}{5}}}\end{array}$$(3)

where η_av is the arithmetic mean of the five independently determined shear viscosity components η_xy, η_xz, η_yz, η_xxyy and η_yyzz.

The shear viscosity can also be determined via the following Einstein relation:¹³⁾

  

$$\eta =\underset{t\to \infty }{lim}\frac{1}{6}\frac{V}{k_{B}T}{\sum }_{\alpha ,\beta }\frac{d}{dt}{⟨{\left({\int }_{t_0}^{t_0+t}{P}_{\alpha \beta }\left(t^{\prime }\right)d{t}^{\prime }\right)}^2⟩}_{t_0}$$(4)

where V is the system volume, k_B is Boltzmann’s constant, T is the absolute temperature, t is the time, P_αβ is the off-diagonal pressure tensor component, α, β refer to x, y, z and α ≠ β.

Alternatively, the shear viscosity can be obtained from a non-equilibrium MD simulation by applying an external perturbation to the system. The shear viscosity relates the maximum shear rate s (equal to the gradient of x-component of the fluid velocity regarding z-direction) with the flux of transverse linear momentum j_z(p_x):¹⁴⁾

  

$$\eta =\frac{-1}{s}{j}_z\left(p_x\right)$$(5)

where $s=\left(\frac{\partial v_x}{\partial z}\right)$. It is important to use a high shear rate to obtain reliable statistics, but not so high that the system moves too far from its equilibrium. The system will be able to relax when the inverse shear rate is longer than the relaxation time.

The SLLOD method¹³⁾ is used to simulate the planar Couette flow by planar movement of the box wall. The momentum flux can be considered as an off-diagonal (eg xz) component of the stress tensor:

  

$$j_z\left(p_x\right)=⟨P_{xz}⟩$$(6)

where P_xz is the off-diagonal pressure tensor component.

The periodic simulation box is divided into slabs along the z coordinate in the Müller-Plathe method.¹⁵⁾ A momentum flux (Eq. (7)) and a fluid velocity gradient are generated by momentum swaps between slabs.

  

$$j_z\left(p_x\right)=\frac{p_x}{2tA}$$(7)

where p_x is the total momentum transferred in a simulation, t is the simulation time and A is the box face area.

The error bar for the shear viscosity in the Müller-Plathe method is calculated as follows:¹⁶⁾

  

$$\mathrm{\Delta }\eta \le \eta \left(\frac{\mathrm{\Delta }J}{J}+\frac{\mathrm{\Delta }G}{G}\right)$$(8)

where J is the flux equal to the slope of the exchanged momentum with the simulation time, determined by a least squares fit, G is the slope of the linear velocity profile, also obtained by a least-squares fit, and the error bars ΔJ and ΔG were both obtained from a least squares fit.

Equations (9) and (10) can be used to obtain the shear rate and momentum flux, respectively, when a periodic external force creates a periodic force field in the periodic perturbation method:¹³⁾

  

$$j_z\left(p_x\right)=\rho \frac{a}{k}$$(9)

where ρ is the density, a is the external acceleration, and:

  

$$s=-vk$$(10)

where v is the velocity profile amplitude, k = 2π/l_z, l_z is the box length.

3. Results and Discussion

3.1. General MD Results

The production stage of MD simulation was performed for a duration of 1 ns for all oxides, and temperatures and the off-diagonal components of the stress tensor were written into a file every five timesteps in the Green-Kubo and Einstein methods. The off-diagonal components of the stress tensor were further used to calculate i) the stress autocorrelation functions (SACFs), which are employed in the Green-Kubo calculation, and ii) the mean square displacement (MSD), which is used in the Einstein method of viscosity calculation.

Figure 1 demonstrates the normalised SACFs obtained for CaO, MgO and Al₂O₃ at 3500, 3600 and 2800 K, respectively. The autocorrelation functions for CaO and MgO decay significantly after 0.5 ps, while for Al₂O₃ it decays after 1 ps. The autocorrelation functions obtained at other temperatures exhibit the same behaviour. The shear viscosity calculated by integrating the SACF using Eq. (2) and the trapezoidal integration method is shown in Figs. 2(a)–2(c) as a function of time (marked as “average”) for CaO, MgO and Al₂O₃ at the aforementioned temperatures. The shear viscosities calculated for each off-diagonal stress tensor component (P_xy, P_xz, P_yz) and two independent differences between the diagonal components (P_xx – P_yy and P_yy – P_zz) are also shown in Figs. 2(a)–2(c). The mean square displacement calculated using the off-diagonal components of the stress tensor (see the expression in the angle brackets in Eq. (4)) is given in Figs. 3(a)–3(c) as a function of time for CaO, MgO and Al₂O₃ at 3500, 3600 and 2800 K, respectively. In accordance with the recommendation by Hess,¹³⁾ the time derivative of MSD was calculated for the initial linear interval of the curve (approximately 5 ps) rather than in the limit of infinite time. The error of the Einstein method of viscosity calculations was estimated in a similar manner to that used for the Green-Kubo method, as expressed by Eq. (3). Figures 4(a)–4(c) depict the shear viscosity calculated by the Einstein relation (4) as a function of time for CaO, MgO and Al₂O₃ at the aforementioned temperatures.

Fig. 1. Normalised stress autocorrelation functions for: CaO, MgO, and Al₂O₃ at 3500, 3600 and 2800 K, respectively.

Fig. 2. Shear viscosity calculated via Green-Kubo equation as a function of time for: (a) CaO at 3500 K, (b) MgO at 3600 K, (c) Al₂O₃ at 2800 K.

Fig. 3. Mean square displacement of off-diagonal components of stress tensor as a function of time for: (a) CaO at 3500 K, (b) MgO at 3600 K, (c) Al₂O₃ at 2800 K.

Fig. 4. Shear viscosity calculated via Einstein relation as a function of time for: (a) CaO at 3500 K, (b) MgO at 3600 K, (c) Al₂O₃ at 2800 K.

During the MD simulations with the Müller-Plathe algorithm, the orthorhombic simulation box was divided into 20 slabs in the z direction (the direction of the momentum flux exchange). During the production stage, the velocity, temperature, density profiles (see Fig. 5) and the transferred momentum were written into a file every 500 timesteps. Initially, simulations were carried out with the momentum exchange period (w) varying at intervals of 1, 10, 20 and 40 timesteps at 3500 K for CaO and MgO, and 2700 K for Al₂O₃, and with a single atom being exchanged between layers. Simulations at the remaining temperatures were carried out with the value of w corresponding to the linear velocity profile and minimal deviations in the temperature and the density profiles. The production stage lasted 1 ns for CaO at 3200–3400 K and MgO at 3100–3300 K, and 500 ps for the rest of temperatures for both oxides. The duration of the production stage in the Al₂O₃ simulation was 1.5 ns at 2400 K, 1.25 ns at 2500 K, 1 ns at 2600 K, and 750 ps for the other temperatures. The moment of time, when the system reached the steady state, was controlled by a dependence of the cumulative shear viscosity upon time as shown in Fig. 6. The viscosity calculation was based solely on the MD data obtained in the steady state.

Fig. 5. MD results for MgO at 3500 K: (a) velocity, (b) temperature, (c) density profiles.

Fig. 6. Viscosity calculated via Müller-Plathe algorithm as a function of time for CaO and MgO at 3500 K, and Al₂O₃ at 2700 K. η* denotes shear-rate-dependent viscosity.

Four distinct trajectories were obtained during the production stage in the MD simulations with the SLLOD method, which lasted 700 ps for CaO at 3200–3400 K, 700 ps for MgO at 3100–3300 K, and 500 ps for the remaining temperatures for both oxides. The duration of the production stage in the Al₂O₃ simulation was 1 ns at 2400–2600 K, and 750 ps for the other temperatures for all four trajectories. During the MD simulations with the SLLOD method, the xz and yz components of the stress tensor, which correspond to the xy direction of shear strain, were written into the file every 100 timesteps. The moment of time, when the system reached the steady state, was controlled by a dependence of the cumulative shear viscosity upon time (see Figs. 7(a)–7(c)). The resulting viscosity was calculated by averaging all trajectories, and the standard deviation was obtained as the error of estimation. Although one trajectory in the Al₂O₃ simulations displays a systematic tendency to be higher than the others at all temperatures (Fig. 7(c)), this is not a sufficient reason to exclude this trajectory from the calculation.

Fig. 7. Viscosity calculated via SLLOD algorithm as a function of time for: (a) CaO at 3500 K, (b) MgO at 3500 K, (c) Al₂O₃ at 2700 K. η* denotes shear-rate-dependent viscosity.

It is a common practice to extrapolate the results of non-equilibrium MD simulations obtained at various shear rates to a zero shear rate. However, this approach was not employed in the Müller-Plathe and SLLOD methods to reduce computational time.

MD simulations with the periodic perturbations method (Eqs. (9)–(10)) were carried out with three distinct values of k and three distinct values of acceleration a for each k at 3300, 3500, 3700 and 3900 K. During the production stage, the value of the velocity amplitude was written into the file every 10 timesteps. In order to calculate viscosity, the mean velocity amplitude value for a steady state (with minor fluctuations) was employed. The production stage was completed in 20 ps for all accelerations and temperatures in the small system simulations. In the moderate-size system simulations, the production stage was completed in 40 ps at 3300 K for all accelerations and 20 ps for all accelerations and remaining temperatures. And, in the large system simulations, the production stage was completed in 40 ps for all accelerations and temperatures. The viscosity values calculated for each a and k were subsequently extrapolated to zero acceleration at a given k using a linear law. Finally, these values were extrapolated to zero k (that corresponds to the infinite system size) using Eq. (11).¹³⁾ The standard deviation of the y-intercept in the last extrapolation was accepted as error for the viscosity calculation.

  

$$\eta \left(k\right)=\eta \left(0\right)\left(1-a{k}^2\right)+O\left(k^4\right)$$(11)

3.2. Calculation Results for Liquid CaO, MgO and Al₂O₃

Figure 8 demonstrates viscosity of the CaO melt calculated by five different MD (Green-Kubo, Einstein, periodic perturbation, SLLOD and Muller-Plathe) methods at different temperatures. Additionally shown a viscosity value at 3200 K calculated via the MD simulated oxygen diffusion coefficient and Stokes-Einstein equation.⁷⁾ It can be seen that all the MD methods used in the present study provide very close results, although two methods used equilibrium MD simulation and three methods employed non-equilibrium MD simulation. It can also be seen that the viscosity value calculated via the oxygen diffusion coefficient is lower than that calculated by the Green-Kubo equation by an order of magnitude. The viscosity of the pure CaO liquid was evaluated in the present work for the first time, and it should be noted that viscosity values obtained by MD simulation in oxide systems can typically be taken as lower estimates of experimental viscosities. Tables 2, 3, 4 provide the actual temperature and viscosity values along with the corresponding uncertainties calculated via the five different MD methods discussed above. The error bars shown in Fig. 8 for the Green-Kubo and Müller-Plathe viscosities were calculated using Eqs. (3) and (8), respectively. The calculation uncertainty of temperature was evaluated statistically from the temperature values recorded during MD simulations.

Fig. 8. Viscosity of CaO melt calculated in the present work by five different MD methods.

Table 2. Viscosity of CaO melt obtained via Green-Kubo and Einstein relations.

		Green-Kubo equation		Einstein relation
Temperature (K)	±ΔT (K)	η (mPa∙s)	±Δη (mPa∙s)	η (mPa∙s)	±Δη (mPa∙s)
3203.37	63.05	2.074	0.065	2.087	0.099
3299.91	65.53	2.004	0.064	1.957	0.078
3399.99	66.52	1.833	0.073	1.79	0.080
3497.96	73.84	1.700	0.049	1.717	0.064
3604.68	74.56	1.546	0.085	1.583	0.144
3701.55	75.36	1.391	0.059	1.353	0.114
3793.29	74.82	1.295	0.044	1.318	0.089
3900.50	79.06	1.065	0.031	1.067	0.060

Table 3. Viscosity of CaO melt obtained via Müller-Plathe and SLLOD methods.

Method	Temperature (K)	±ΔT (K)	Shear rate (ps−1)	η (mPa∙s)	±Δη (mPa∙s)
Müller-Plathe	3199.04	57.03	0.16	2.166	0.009
	3300.59	59.81	0.17	1.953	0.010
	3399.49	61.81	0.18	1.809	0.007
	3500.65	64.40	0.20	1.686	0.012
	3600.69	63.70	0.21	1.546	0.010
	3698.06	66.79	0.23	1.421	0.005
	3800.58	68.30	0.24	1.316	0.004
	3902.08	69.41	0.26	1.226	0.006
	3997.16	73.78	0.28	1.110	0.006
SLLOD	3200.02	68.45	0.15	2.212	0.017
	3299.57	71.06	0.17	2.007	0.024
	3399.80	72.10	0.18	1.865	0.017
	3499.91	74.97	0.19	1.706	0.026
	3600.18	75.59	0.20	1.591	0.025
	3700.23	79.55	0.21	1.484	0.016
	3800.06	81.67	0.22	1.388	0.019
	3900.82	84.42	0.25	1.295	0.019
	3999.84	87.86	0.26	1.207	0.015

Table 4. Viscosity of CaO melt obtained via periodic perturbation method.

Temperature (K) ±ΔT (K)	k, A−1	a, A∙ps−2	η (mPa∙s)
3300.34±24.25	0.312	100	0.585
		60	1.039
		15	1.531
	0.156	13	1.494
		7	1.758
		1.6	1.795
	0.104	3	1.850
		1.6	1.985
		0.3	1.827
3500.32±26.8	0.308	90	0.300
		60	0.944
		15	1.710
	0.153	13.5	1.197
		6.5	1.474
		1.35	1.731
	0.103	3.5	1.599
		0.35	1.515
3698.95±28.39	0.304	90	0.586
		60	0.802
		15	1.256
	0.152	13	1.104
		7	1.297
		1.6	1.380
	0.101	3	1.377
		1.6	1.440
		0.3	1.430
3899.81±29.23	0.300	80	0.378
		50	0.843
		15	1.089
	0.150	13	0.949
		7	1.136
		1.6	1.184
	0.100	3	1.153
		1.6	1.188
		0.3	1.098

Figure 9 shows viscosities of the MgO melt calculated in this work by two equilibrium MD methods (Eqs. (1) and (3)) and two non- equilibrium MD methods (Eqs. (5), (6), (7)) and obtained in two similar studies: i) Karki et al.¹⁷⁾ by the ab initio MD simulation based on the density functional theory, and ii) Leu et al.¹⁸⁾ using the theory of the significant liquid structures developed by Eyring and co-workers.¹⁹⁾ Karki et al.¹⁷⁾ employed the density functional theory to obtain the interaction force field, which was then used to carry out MD simulation and further to calculate viscosity using the Green-Kubo equation (Eq. (2)). It can be seen that the MD results obtained in the present work are very close to the MD results by Karki et al.¹⁷⁾ Tables 5–6 provide the actual temperature and viscosity values along with the corresponding uncertainties calculated via the four different MD methods discussed above.

Fig. 9. Viscosity of MgO melt calculated in the present work and values by Karki et al.¹⁷⁾ and by Leu et al.¹⁸⁾

Table 5. Viscosity of MgO melt obtained via Green-Kubo and Einstein relations.

		Green-Kubo equation		Einstein relation
Temperature (K)	±ΔT (K)	η (mPa∙s)	±Δη (mPa∙s)	η (mPa∙s)	±Δη (mPa∙s)
3101.04	64.92	2.437	0.08	2.463	0.08
3197.66	64.18	2.264	0.114	2.283	0.254
3293.56	66.93	1.982	0.155	2.015	0.181
3402.35	70.12	1.767	0.084	1.793	0.072
3497.36	73.09	1.74	0.035	1.753	0.115
3601.25	69.56	1.531	0.117	1.456	0.097
3696.71	72.88	1.452	0.091	1.481	0.145
3799.76	73.77	1.337	0.048	1.341	0.17
3897.37	80.97	1.371	0.048	1.325	0.042
3997.51	79.85	1.249	0.035	1.291	0.185

Table 6. Viscosity of MgO melt obtained via Müller-Plathe and SLLOD methods.

Method	Temperature (K)	±ΔT (K)	Shear rate (ps−1)	η (mPa∙s)	±Δη (mPa∙s)
Müller-Plathe	3100.99	55.81	0.18	2.376	0.011
	3200.48	57.33	0.19	2.209	0.01
	3299.85	58.67	0.22	1.974	0.007
	3396.07	61.43	0.23	1.831	0.012
	3501.55	61.21	0.24	1.723	0.008
	3599.44	65.21	0.24	1.583	0.007
	3701.53	66.82	0.28	1.467	0.006
	3802.19	68.62	0.30	1.383	0.004
	3898.55	67.84	0.31	1.301	0.004
	3998.24	70.83	0.33	1.199	0.003
SLLOD	3099.93	66.47	0.17	2.396	0.018
	3200.21	68.46	0.18	2.223	0.033
	3300	70.99	0.20	2.046	0.022
	3399.77	72.41	0.22	1.873	0.024
	3499.89	75.22	0.23	1.756	0.02
	3599.77	77.52	0.25	1.641	0.027
	3700.22	79.44	0.27	1.528	0.03
	3800.23	80.61	0.28	1.443	0.011
	3899.73	84.17	0.30	1.377	0.02
	3999.52	86.51	0.32	1.275	0.015

Viscosity of the Al₂O₃ melt was calculated by four different methods: Green-Kubo and Einstein (equilibrium MD) and Müller-Plathe and SLLOD (non-equilibrium MD). Figure 10 shows the MD results and the viscosity error bars evaluated by Eqs. (3) and (8). It can be seen that viscosities calculated by equilibrium MD methods (the Green-Kubo and Einstein equations) are slightly higher than viscosities calculated by non-equilibrium MD (the Muller-Plathe and SLLOD methods). Tables 7–8 provide the actual temperature and viscosity values along with the corresponding uncertainties calculated via the four different MD methods discussed above.

Fig. 10. Viscosity of Al₂O₃ melt calculated in the present work.

Table 7. Viscosity of Al₂O₃ melt obtained via Green-Kubo and Einstein relations.

		Green-Kubo equation		Einstein relation
Temperature (K)	±ΔT (K)	η (mPa∙s)	±Δη (mPa∙s)	η (mPa∙s)	±Δη (mPa∙s)
2401.54	49.88	11.637	0.646	11.571	1.074
2497.68	51.46	11.377	0.965	10.7	0.625
2600.67	50.26	9.256	0.418	8.718	0.364
2697.43	51.72	7.945	0.26	7.734	0.338
2799.98	55.17	7.454	0.359	7.514	0.495
2901.24	60.39	6.128	0.403	5.964	0.625
3003.07	58.23	5.752	0.648	5.764	0.59
3099.1	60.35	5.17	0.456	5.148	0.767
3196.16	62.94	5.131	0.537	5.207	0.615
3297.53	64.16	4.46	0.302	4.452	0.44

Table 8. Viscosity of Al₂O₃ melt obtained via Müller-Plathe and SLLOD methods.

Method	Temperature (K)	±ΔT (K)	Shear rate (ps−1)	η (mPa∙s)	±Δη (mPa∙s)
Müller-Plathe	2401.68	45.16	0.17	10.124	0.102
	2499.92	46.01	0.20	8.808	0.05
	2597.21	47.15	0.22	7.777	0.062
	2700.29	49.86	0.25	6.924	0.032
	2801.1	52.27	0.27	6.225	0.034
	2901.54	53.13	0.30	5.611	0.022
	3000.4	56.27	0.32	5.218	0.019
	3100.16	58.84	0.35	4.721	0.014
	3198.83	59.89	0.38	4.349	0.013
	3298.48	60.41	0.41	4.023	0.016
SLLOD	2400.1	51.82	0.16	9.341	0.328
	2500.07	54.7	0.19	8.076	0.284
	2600.23	57.26	0.21	7.268	0.237
	2700.36	58.14	0.23	6.554	0.248
	2799.73	60.08	0.25	5.98	0.201
	2900	62.39	0.27	5.434	0.181
	2999.98	64.33	0.29	5.012	0.156
	3099.93	66.35	0.31	4.613	0.147
	3199.83	68.42	0.33	4.289	0.145
	3299.68	71.79	0.35	3.986	0.125

Since the Müller-Plathe and SLLOD algorithms are non-equilibrium molecular dynamics methods, viscosity calculated by these methods depends on the shear rate, exhibiting a decrease with an increase in shear rate.²⁰⁾ This could explain the discrepancy between the MD results obtained for Al₂O₃ by equilibrium methods (Green-Kubo and Einstein) and those derived from non-equilibrium methods (Müller-Plathe and SLLOD). It can be observed that the discrepancy is more pronounced at lower temperatures (higher viscosities) than at higher temperatures (lower viscosities). Viscosities of CaO and MgO were calculated by the Müller-Plathe and SLLOD methods at significantly higher temperatures, so this discrepancy may not be as noticeable as in the case of Al₂O₃.

3.3. Comparison and Discussion

For comparison with the available viscosity data (experiments or model calculations) for the CaO, MgO and Al₂O₃ melts it may be helpful to fit MD results into an Arrhenius-type viscosity-temperature expression, which parameters can be used for further analysis and discussion. In the present study the classical Arrhenius equation (11) was employed:

  

$$\eta =Aexp\left(\frac{E_a}{RT}\right)$$(11)

where A is the pre-exponential term, T is the absolute temperature, R is the universal gas constant and E_a is the activation energy of viscous flow. The Arrhenius parameters obtained by fitting MD results into Eq. (11) for MgO, CaO and Al₂O₃ are given in Table 9. Since the MD results for CaO and MgO are close to each other for all the calculation methods, the values of A and E_a for CaO and MgO were averaged over five and four different MD methods, respectively, while the Arrhenius parameters for Al₂O₃ were calculated for each of the four MD methods used.

Table 9. Parameters of the Arrhenius equation for the CaO, MgO and Al₂O₃ viscosities obtained by different MD methods.

Oxide	Method	logA (Pa·s)	Ea (kJ·mol−1)
CaO	Average	−4.10	88.07
MgO	Average	−3.90	76.39
Al2O3	Green-Kubo	–3.52	73.71
	Einstein	–3.46	70.07
	Müller-Plathe	–3.46	67.20
	SLLOD	–3.37	61.46

As no experimental data is available for the CaO melt, Fig. 11 shows comparison between viscosities (as the decimal logarithm vs reciprocal temperature) calculated by four different (FactSage, Kondratiev-Jak, Starodub and GTT) viscosity models and the average Arrhenius dependence obtained in the present work. It can be seen that the slope of the Arrhenius dependence, which is proportional to the activation energy of viscous flow (Eq. (11)), is similar to the values obtained from the above viscosity models. It can also be seen from Fig. 11 that the FactSage, GTT and Kondratiev-Jak models underestimate viscosity of the pure CaO melt, as these model curves lie substantially below the MD average curve. It may indicate that extrapolation from the corresponding binary and ternary systems does not always provide a correct estimation of viscosity of a pure oxide.

Fig. 11. Comparison of CaO viscosity calculated via MD simulation (average Arrhenius) and by different viscosity models.

Figure 12 shows comparison between the average Arrhenius viscosity and calculations by the aforementioned viscosity models for the MgO melt in the logarithmic format as a function of the reciprocal temperature. It can be seen that the GTT model might underestimate viscosity of pure MgO, while the viscosity curve by the Kondratiev-Jak model provides a more realistic estimation of the MgO melt viscosity. It can be seen from Fig. 12 that the activation energy obtained in the present study is close to that of the Kondratiev-Jak model, and lower than those of the other models.

Fig. 12. Comparison of MgO viscosity calculated via MD simulation (average Arrhenius) and by different viscosity models.

The available experimental data²¹⁾ for Al₂O₃ were roughly divided into two large groups dependent on a method used to determine viscosity: those obtained via conventional viscometry methods (rotational cylinder techniques, oscillation crucible technique) and those measured indirectly by levitation methods. The rotational cylinder techniques involve mainly the traditional rotating crucible method,²²⁾ while the oscillation crucible technique represents the so-called Shvidkovskii method of dumping torsional oscillations of a cylindrical crucible.²³⁾ Both methods are suitable for viscosity measurement of liquid Al₂O₃ from its melting point up to 2300–2400°C. The levitation methods involve indirect determination of the melt viscosity (altogether with the surface tension) from damping oscillations induced on a levitated drop²⁴⁾ and are characterized a bigger experimental uncertainty (reproducibility) than traditional viscosity measurements.

Figure 13 demonstrates agreement of the MD simulation results obtained in this work for Al₂O₃ (presented as the four Arrhenius dependencies, which parameters are given in Table 9) with the available experimental data²¹⁾ and other MD simulation by Jahn and Madden²⁵⁾ and Jahn,²⁶⁾ who employed the aspherical ion model with more complex interatomic interactions (including ion shape deformations and polarisation) and the Green-Kubo equation. It can be seen that viscosities calculated by equilibrium MD methods (the Green-Kubo and Einstein equations) are slightly closer to the experimental values obtained by levitation than viscosities calculated by non-equilibrium MD (the Muller-Plathe and SLLOD methods). The activation energies obtained in the present study vary from 61.5 to 73.7 kJ·mol⁻¹, which is generally lower than those reported by conventional viscosity measurements^22,23) (95–128 kJ·mol⁻¹) but higher than those obtained by levitation methods²⁴⁾ (43.2 kJ·mol⁻¹).

Fig. 13. Comparison of calculated (MD) and experimental viscosity for Al₂O₃ melt.

It may be interesting to compare the activation energies of viscous flow obtained by MD for Ca and Mg oxides with those obtained for other compounds of Ca and Mg. For example, a recent study by Takeda and co-authors²⁷⁾ of fluorides viscosities reported the activation energies for CaF₂ and MgF₂ to be 47.3 and 32.4 kJ·mol⁻¹, respectively. These authors suggested that, since the molten fluorides are ionic liquids, the cohesive Coulomb force plays a significant role in viscosity: the stronger Coulomb interactions, the higher viscosity. The Coulomb force is directly related to the radius and valence of ions, which was supported by a decrease of viscosity of alkali fluorides from Li to Cs. However, viscosity of MgF₂ was reported in the same work to be lower than that of CaF₂, which does not support this assumption. Additionally, the viscosities of alkaline earth chlorides presented in the same paper²⁷⁾ demonstrate a reverse trend: a decrease in viscosity in the order BaCl₂ > SrCl₂ > CaCl₂ > MgCl₂. Moreover, the activation energies of viscous flow of alkaline earth fluorides are similar to each other, which implies that the electrostatic interaction and ion radius both have a little effect on this parameter. Nevertheless, the activation energies for CaO and MgO obtained in the present work are higher than those of the corresponding fluorides and also chlorides,²⁶⁾ which may indicate that the interatomic interactions in these oxides are stronger than that in corresponding salts.

It is not quite clear what part of the force field is responsible for higher activation energies of oxides. Nevertheless, several authors^24,25,28,29) have achieved a notable enhancement in calculation of certain properties for molten oxide systems, including viscosity, by incorporating polarisation effects and modification of short-range interactions in the force field. The latter interactions are particularly important in calculation of the pressure tensor that is used to compute viscosity.

As the present study aimed to compute and compare MD results obtained by various viscosity calculation methods for three oxides across a wide temperature range, a simple empirical potential (like BMH) was used to ensure faster simulations and lower computational cost than more complex potentials that take into account ion-dipole and dipole-dipole interactions. The detailed investigation of the influence of various force fields on oxide properties will be the subject of subsequent work.

In total, the viscosity calculation results obtained in the present work indicate that MD simulation even with a simple interatomic potential (such as the BMH potential) can provide a reasonable estimation of viscosity. However, MD results obtained with the BMH potential should be considered as a lower estimate. For a more reliable estimation of viscosity a more complex potential should be used.

4. Conclusions

Equilibrium and non-equilibrium molecular dynamics simulations coupled with various viscosity calculation methods (e.g. Green-Kubo, Müller-Plathe, SLLOD) were carried out for liquid CaO, MgO, and Al₂O₃ using the well-known Born-Mayer-Huggins interatomic potential, which parameters were extracted from the literature.

Viscosity of the pure CaO melt was evaluated by MD simulation for the first time using two equilibrium and three non-equilibrium MD methods of viscosity calculation. The MD results showed very close agreement for each other. Since no experimental data is available for CaO, its viscosity was compared to predictions of the selected viscosity models. It was demonstrated that the MD results in the form of the Arrhenius dependence are higher than most of the model predictions and can now be used as a more reasonable estimation for CaO viscosity.

Viscosity of the pure MgO melt evaluated by MD simulation was compared to two theoretical calculations and showed a reasonable agreement with the other MD simulation that employed the density functional theory to obtain the interatomic force field. Again, different MD viscosity calculation methods employed provide similar results for viscosity of the MgO melt. Comparison of the MD results (as the average Arrhenius dependence) with other viscosity model predictions showed that it is agreed reasonably well with the predictions, and the Arrhenius activation energy obtained by MD is close to those values calculated by the viscosity models.

Viscosity of the pure Al₂O₃ melt was evaluated by MD simulation and was then compared to the available experimental data. The comparison demonstrated that MD simulation with a simple interatomic potential (like BMH) can provide a lower estimation of viscosity in relation to the experimental points. It should be mentioned that equilibrium MD methods provide slightly higher viscosity values for the Al₂O₃ melt than non-equilibrium MD methods. The activation energy values obtained for Al₂O₃ were lower than those calculated from viscosities obtained by conventional viscometry methods and higher than those obtained from levitation methods.

In general, the activation energies of CaO and MgO obtained by MD calculations were higher than the values obtained from viscosity data of corresponding fluorides and chlorides, which may indicate that the interatomic interactions in these oxides are stronger than that in corresponding salts. For a more reliable estimation of viscosity a more complex potential should be used, which will be the scope of further work.

Statement for Conflict of Interest

The authors declare no conflicts of interest related to the conduct of this research.

Acknowledgments

This work was financially supported within the framework of the State Program ‘Chemical Thermodynamics and Theoretical Materials Science’ at the Lomonosov Moscow State University (N◦121031300039-1). The authors would also like to express their gratitude to S. Kuzovchikov (Chemical Thermodynamics laboratory, Lomonosov Moscow State University) for his assistance in MD data processing.

References

• 1)   A. N.  Grundy,  H.  Liu,  I. H.  Jung,  S. A.  Decterov and  A. D.  Pelton: Int. J. Mater. Res., 99 (2008), 1185. https://doi.org/10.3139/146.101752
• 2)   G.  Wu,  E.  Yazhenskikh,  K.  Hack,  E.  Wosch and  M.  Müller: Fuel Process. Technol., 137 (2015), 93. https://doi.org/10.1016/j.fuproc.2015.03.025
• 3)   A.  Kondratiev and  E.  Jak: Fuel, 80 (2001), 1989. https://doi.org/10.1016/S0016-2361(01)00083-7
• 4)   L.  Forsbacka,  L.  Holappa,  A.  Kondratiev and  E.  Jak: Steel Res. Int., 78 (2007), 676. https://doi.org/10.1002/srin.200706269
• 5)   K.  Starodub,  G.  Wu,  E.  Yazhenskikh,  M.  Müller,  A.  Khvan and  A.  Kondratiev: Ceram Int., 45 (2019), 12169. https://doi.org/10.1016/j.ceramint.2019.03.121
• 6)   A. P.  Thompson,  H. M.  Aktulga,  R.  Berger,  D. S.  Bolintineanu,  W. M.  Brown,  P. S.  Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton: Comp. Phys. Comm., 271 (2022), 10817. https://doi.org/10.1016/j.cpc.2021.108171
• 7)   C. M. S.  Alvares,  G.  Deffrennes,  A.  Pisch and  N.  Jakse: J. Chem. Phys., 152 (2020), 084503. https://doi.org/10.1063/1.5141841
• 8)   A. S.  Arkhipin,  A.  Pisch,  I. A.  Uspenskaya and  N.  Jakse: Ceramics, 7 (2024), 1187. https://doi.org/10.3390/ceramics7030078
• 9)   G.  Gutierrez,  A. B.  Belonoshko,  R.  Ahuja and  B.  Johansson: Phys. Rev. E, 61 (2000), 2723. https://doi.org/10.1103/PhysRevE.61.2723
• 10)   G. J.  Martyna,  D. J.  Tobias and  M. L.  Klein: J. Chem. Phys., 101 (1994), 4177. https://doi.org/10.1063/1.467468
• 11)   M. P.  Allen and  D. J.  Tildsley: Computer Simulation of Liquids, 2nd ed., Oxford University Press, New York, (2017), 641. ISBN 978–0–19–880320–1
• 12)   D.  Nevins and  F. J.  Spera: Mol. Sim., 33 (2007), 1261. https://doi.org/10.1080/08927020701675622
• 13)   B.  Hess: J. Chem. Phys., 116 (2002), 209. https://doi.org/10.1063/1.1421362
• 14)   B. D.  Todd and  P. J.  Daivis: Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications, Cambridge University Press, Cambridge, (2017), 359. ISBN 978-0-521-19009-1
• 15)   F.  Müller-Plathe: Phys. Rev. E, 59 (1999), 4894. https://doi.org/10.1103/PhysRevE.59.4894
• 16)   P.  Bordat and  F.  Müller-Plathe: J. Chem. Phys., 116 (2002), 3362. https://doi.org/10.1063/1.1436124
• 17)   B. B.  Karki,  J.  Zhang and  L.  Stixrude: Geophys. Res. Lett., 40 (2013), 94. https://doi.org/10.1029/2012GL054372
• 18)   A.-L.  Leu,  S.-M.  Ma and  H.  Eyring: Proc. Nat. Acad. Sci. USA, 72 (1975), 1026. https://doi.org/10.1073/pnas.72.3.1026
• 19)   H.  Eyring and  M. S.  Jhon: Significant Liquid Structures, John Wiley & Sons, New York, (1969), 149. ISBN 471-24982-3
• 20)   P. T.  Cummings and  G. P.  Morriss: J. Phys. F: Metal Phys., 17 (1987), 593. https://doi.org/10.1088/0305-4608/17/3/007
• 21)  The oxide viscosity database OxiVis, 2004-2024, A. Kondratiev, (private database, available on request: [e-mail]).
• 22)   G.  Urbain: Rev. Int. Haut. Temp. Refract., 19 (1982), 55.
• 23)   V. P.  Elyutin,  V. I.  Kostikov,  B. C.  Mitin and  Y. A.  Nagibin: Zh. Fiz. Khim. (Rus. J. Phys. Chem.), 43 (1969), 579.
• 24)   P.-F.  Paradis and  T.  Ishikawa: Jpn. J. Appl. Phys., 44 (2005), 5082. https://doi.org/10.1143/JJAP.44.5082
• 25)   S.  Jahn and  P. A.  Madden: J. Non-Cryst. Solids, 353 (2007), 3500. https://doi.org/10.1016/j.jnoncrysol.2007.05.104
• 26)   S.  Jahn: Am. Miner., 93 (2008), 1486. https://doi.org/10.2138/am.2008.2873
• 27)   O.  Takeda,  Y.  Hoshino,  Y.  Anbo,  K.-I.  Yanagase,  M.  Aono and  Y.  Sato: Int. J. Thermophys., 36 (2015), 648. https://doi.org/10.1007/s10765-014-1828-5
• 28)   T.  Dufils,  N.  Folliet,  B.  Mantisi,  N.  Sator and  B.  Guillot: Chem. Geol., 461 (2017), 34. https://doi.org/10.1016/j.chemgeo.2016.06.030
• 29)   Y.  Ishii,  M.  Salanne,  T.  Charpentier,  K.  Shiraki,  K.  Kasahara and  N.  Ohtori: J. Phys. Chem. C, 120 (2016), 24370. https://doi.org/10.1021/acs.jpcc.6b08052