Molecular Dynamics Simulation for Structural Evolution of Mixed Ash from Coal and Wheat Straw

: We conducted molecular dynamics (MD) simulations to investigate the structural evolution of molten slag composed of wheat straw (WS) and Shenhua (SH) coal. The content of wheat straw in the slag was varied from 0 to 100 wt%. The MD results indicated a slight reduction in the sharpness of the radial-distribution-function curve of each ion–oxygen pair and a decrease in bonding strength with increasing WS content. WS introduced many metal ions to the ash system, increasing its overall activity. The number of bridging and non-bridging oxygen atoms changed upon straw addition, which affected the stability of the system. There were relatively few highly coordinated Si ions. The number of low-coordination Si was highest for a WS content of 30%, at which the density reached a minimum value. The degree of ash polymerization was analyzed by counting the number (Q) of tetrahedra with the number (n) of the bridging oxygen atoms. With increasing WS content, Q 4 (tetrahedral Si) decreased, whereas Q 3 , Q 2 , Q 1 , and Q 0 increased. Q 4 reached a minimum value for a WS content of 30%, at which point the degree of ion aggregation was the weakest and the degree of disorder was the strongest.


Introduction
Coal gasification is important for the high-efficiency utilization of coal, production of liquid fuel and hydrogen, reduction of iron, and other processes. Co-gasification of biomass and coal increases the energy density of the biomass and decreases the clinkering rate of coal ash owing to its high fluidity. Understanding the network structure for the cogasification ash of biomass and coal is important when designing the gasifier and selecting the operating temperature and fuel ratio [1,2]. Shenhua (SH) coal from the China Shenhua coalfield has good surface activity, high calorific value, and low sulfur and ash content. In particular, the melting temperature of SH coal ash is lower than other coal, which is suitable for a flow bed gasifier with liquid slag discharge. Coal ash is mainly composed of metal oxides and nonmetallic oxides [3,4]. Ions can be divided into three types on the basis of their different structures in slag: network ions, modified ions, and neutral ions. The roles of SiO 2 , Al 2 O 3 , Fe 2 O 3 , and CaO in melt have been studied extensively [5,6]. Because of the high SiO 2 content in coal ash, the structure of silicates is of particular importance for understanding the structure and behavior of slags. Silicate slags contain Si cations surrounded by four tetrahedrally arranged oxygen anions. These SiO 4 4− tetrahedra are joined in chains or rings by bridging oxygen (BO) atoms. Cations, such as Na + , Ca 2+ , Mg 2+ , K + , and Fe 2+ , break these bonded oxygens to form non-bridging oxygen (NBO), O − , and free oxygen, O 2− [7,8].
Compared with coal ash, the content of alkali metal elements such as K, Na, Mg, and Ca in biomass ash is high, and the content of network ions (Si and Al) is low [9]. The addition of biomass ash (such as wheat straw, corn straw, etc.) can reduce the melting point of coal ash [8], changing the flow characteristics of ash at high temperatures [10].
MD is an excellent tool to study the structure with classic dynamics and obtain information such as the chemical bond. This method has already been used in slag [11][12][13][14][15][16]. For example, experimental and MD simulations results indicated that the SiO 4 4− tetrahedron is the most stable unit in the slag, and the bonding stability decreases in the order: Si-O > Al-O > Fe-O > Mg-O > Ca-O [13]. Zhang et al. [14] also point out that alkalis (Na 2 O and K 2 O) have only marginal influence on the structure of molten aluminosilicates, while the Na + ion prefers to locate in the BO/NBO networks, and the K + ion tends to be present in various oxygen tri-clusters. Moreover, Si-O and Al-O networks were observed to depolymerize into a simple structure at high alkalinity, with atomic self-diffusion increasing and the viscosity decreasing with increasing alkalinity [15]. The decrease in viscosity led to an increase in the proportion of NBO in the network. In another study, an increase in the MgO:Al 2 O 3 ratio led to the introduction of more Mg 2+ ions to destroy the network structure, resulting in the conversion of BO to NBO [16]. Moreover, as the basicity increased, more BO was converted to free oxygen.
At present, there are plenty of studies on the structure of coal ash, but few studies on the structure of mixed ash from coal and biomass. In this study, we carried out MD simulations to investigate the effect of wheat straw (WS) content (varied from 0 to 100 wt%) on the structural properties of ash slag, including the radial distribution functions (RDF), mean square displacement (MSD), coordination numbers (CN), distribution of oxygen types, distribution of bond angles, and distribution of Q n .

Ash Samples Preparation
SH coal was supplied by the China Energy Investment Corporation; WS was collected from Bingjiang farm in Zhengjiang city. The industrial analysis of SH coal and WS was carried out using SH coal and WS powder (<100 µm) under the PRC National Standard GB/T 212-2008 [17] and GB/T28731-2012 [18], respectively. The analysis results are shown in Table 1. SH coal ash and WS ash were prepared under the ASTM D3174-12 standard in which the final ashing temperature was set at 575 • C to minimize the loss of volatile elements such as Na [19]. The ash composition was analyzed by XRF (ARL ADVANT'X Intellipower 4200, ThermoFisher Scientific, USA), as shown in Table 2. Simulation of real ash is difficult as as there are some impurities (P 2 O 5 , SO 3 , and Fe 2 O 3 ) in real ash, as shown in Table 2. In this simulation study, six main components in coal ash and WS ash were selected for calculation to simplify the ash system, and the molar compositions are shown in Table 3, of which the content of WS in the slag was varied from 0% to 100 wt%.

Simulation Method
We used the Amorphous Cell Tools module of the Materials Studio 7.1 (MS) software to build an amorphous structural system model. The particles were randomly distributed in the system model to simulate the melt at high temperature, as shown in Figure 1. The density of the system determined the size of the cubic box.
CFF (Compass force filed) and PCFF (polymer consistent force field) in the Amorphous Cell Tools module of Materials Studio are widely used to build irregular cell structure of inorganic compounds. However, an oxygen atom cannot be added in CFF, only in the form of oxide, which is not adjusted for random distribution of atoms. In this work, the PCFF was selected for building the model as O can be added in the form of randomly distributed oxygen atoms. Other atoms were directly added to the force field, namely, Si, Na, Ca, Mg, Al, and K.
Periodic boundary conditions were applied on all sides of the model box to create an infinite system. The time step was set as 1 fs, and the data were saved each 100 steps. In the simulation system, the total number of atoms was set to 4000 and the number of different atoms was decided according to their mole fractions. The molar composition of the different ash samples is shown in Table 3. We selected the atoms with high mole fractions, i.e., Na, Mg, Al, Si, K, Ca, and O. The effective charge of ions was set as follows: Si was +1.89, Al was +1.4175, Ca was +0.945, Mg was +0.945, Na was +0.4725, K was +0.4725, and O was −0.945. Processes 2022, 9, x FOR PEER REVIEW 3 of 13 coal ash and WS ash were selected for calculation to simplify the ash system, and the molar compositions are shown in Table 3, of which the content of WS in the slag was varied from 0% to 100 wt%.

Simulation Method
We used the Amorphous Cell Tools module of the Materials Studio 7.1 (MS) software to build an amorphous structural system model. The particles were randomly distributed in the system model to simulate the melt at high temperature, as shown in Figure 1. The density of the system determined the size of the cubic box.
CFF (Compass force filed) and PCFF (polymer consistent force field) in the Amorphous Cell Tools module of Materials Studio are widely used to build irregular cell structure of inorganic compounds. However, an oxygen atom cannot be added in CFF, only in the form of oxide, which is not adjusted for random distribution of atoms. In this work, the PCFF was selected for building the model as O can be added in the form of randomly distributed oxygen atoms. Other atoms were directly added to the force field, namely, Si, Na, Ca, Mg, Al, and K.
Periodic boundary conditions were applied on all sides of the model box to create an infinite system. The time step was set as 1 fs, and the data were saved each 100 steps. In the simulation system, the total number of atoms was set to 4000 and the number of different atoms was decided according to their mole fractions. The molar composition of the different ash samples is shown in Table 3. We selected the atoms with high mole fractions, i.e., Na, Mg, Al, Si, K, Ca, and O. The effective charge of ions was set as follows: Si was +1.89, Al was +1.4175, Ca was +0.945, Mg was +0.945, Na was +0.4725, K was +0.4725, and O was −0.945. The number and density of each atom and the side length of the cubic model box in the system are listed in Table 4. The main components in the ash slag were feldspar (2.55-2.67 g/cm 3 ) and diopside (3.27-3.38 g/cm 3 ) [20]. Because the addition of biomass increases the diopside content in the ash, we selected a slag density of 3.0 g/cm 3 . Since the sizes of The number and density of each atom and the side length of the cubic model box in the system are listed in Table 4. The main components in the ash slag were feldspar (2.55-2.67 g/cm 3 ) and diopside (3.27-3.38 g/cm 3 ) [20]. Because the addition of biomass increases the diopside content in the ash, we selected a slag density of 3.0 g/cm 3 . Since the sizes of Na, Mg, Al, Si, K, Ca, and O atoms are different, there are differences in the length of the cubic model box to ensure the same total atoms and slag density in different slag systems.

MD Simulation Process
We used potential functions based on the Garofalini potential function [21], which includes a pair potential and three-body potential. The Born-Mayer-Huggins (BMH) function was used as a pair potential function, which is written as where r ij is the interatomic distance between atoms i and j, q i and q j represent the charges of atoms i and j, and ξ represents correction function. A ij , ρ ij , and β ij represent the potential parameters of BMH. The parameters used are listed in Table 5. The form of the three-body potential function can be written as: where R ij is the distance between particles i and j; R jk is the distance between particles j and k; λ ij and γ ij are the potential parameters between particles i and j; λ jk and γ jk are the potential parameters between particles j and k; and ω ijk is a function related to the angle. The parameters of the three potentials are taken from Garofalini's work [22], as shown in Table 6. Simulations were carried out in the NVT ensemble, i.e., keeping the number of particles (N), system volume (V), and temperature (T) fixed. As shown in Figure 2, the system was thermally stabilized for 25 ps at T = 5000 K to ensure full mixing of the particles and was then cooled to 2000 K in 25 ps with a cooling rate of 1.2 × 10 14 K/s. To achieve a state of equilibrium, the system was kept at 2000 K for 25 ps and then cooled to 1623 K in 25 ps with a cooling rate of 1.508 × 10 13 K/s. In the final stage, the system was equilibrated at 2223 K for 25 ps to obtain a uniform ash slag in a molten state.

Radial Distribution Functions (RDFs)
From RDFs, basic structural information, including the degree of material ordering and the degree of correlation of electrons, can be obtained. The formula for the RDF can be written as [23]: where V is the volume of the system; n ij is the number of i atoms included in a spherical shell of thickness; ∆r is located at a distance r; and N is the total number of particles. The distributions of the RDF for each particle in the ash sample are shown in Figure 3. The location of the first peak in the RDF graph represents the average bond length of the atomic pair. For Si-O, the location of first peak in all five samples was unaffected by the addition of WS, which corresponds to a bond length of 1.61 Å. The locations of the first peaks seen in Al-O for SH, SH-WS (10%), SH-WS (30%), and SH-WS (50%) were similar, although the addition of WS shifted the peak slightly, indicating a decrease in bond length from 1.75 to 1.73 Å. The peak values of Si-O and Al-O were consistent with previous reports [21], indicating that the calculated results are reasonable.
The sharpness of the first peak in the RDF curves reflect the bonding ability and interaction between ions [24]. With increasing WS content, the width of the first peak of

Mean Square Displacement
The atoms do not stay in a fixed position in the system; rather, they are in continuous movement. The mean square displacement (MSD) of each atom can be determined from its trajectory as [25]: where ri(0) is the displacement of atom i at time zero, and ri(t) represents the displacement of atom i at time t. The MSD has a linear relationship with the diffusion coefficient (DC), which can be written as [26]:

Mean Square Displacement
The atoms do not stay in a fixed position in the system; rather, they are in continuous movement. The mean square displacement (MSD) of each atom can be determined from its trajectory as [25]: where r i (0) is the displacement of atom i at time zero, and r i (t) represents the displacement of atom i at time t. The MSD has a linear relationship with the diffusion coefficient (DC), which can be written as [26]: The MSDs of Si, Al, Ca, Mg, and K, and the total displacement of all atoms in the system are shown in Figure 4. For all samples, the activity and diffusion coefficient of Si atoms were the lowest, because Si formed the network structure of the system. However, with increasing WS content, the activity of other metal oxides (e.g., K 2 O, MgO, and CaO) in the ash increased. For a WS content of 0-10 wt%, the activity of Ca was the highest among all elements in the ash. Between 30% and 100 wt%, K had the highest activity. The high activity of Ca and K is because these elements do not form the network structure and can move freely in the network. As shown in Figure 4F, the total MSD increased with increasing WS content. A possible explanation for this trend is that the addition of more oxide reduces the complexity of the system. We calculated the diffusion coefficients in the system on the basis of the law of diffusion, as shown in Figure 5. The diffusion coefficients decreased in the following order: D Ca > D Mg > D K > D Al > D Si at WS contents up to 10%. The order of diffusion coefficients was D K > D Mg > D Ca for WS contents between 30% and 100%, and at 100%, D Si > D Al .

Coordination Numbers
The microstructure of the Na2O-MgO-Al2O3-SiO2-K2O-CaO slag system was mainly composed of the silicoaluminate and the network modifier. The coordination numbers (CN) varied with the fraction of WS, as shown in Figure 6. It could be seen that the coordination number of Si was maintained at about four, but decreased with increasing WS content, indicating a reduction in the strength of Si-O bond. We attributed this to the depolymerization of the highly coordinated Si by metal atoms, such as Mg and K. For Al, not only four-coordinate Al-O structures but also a large number of structures with fewer than four coordination bonds formed. With increasing WS content, the fraction of more highly coordinated Al-O structures increased. This indicates that the [AlO4] tetrahedral unit is not as stable as the [SiO4] tetrahedron.

Distribution of Oxygen Types
The oxygen ions in this system can be classified as either BO or NBO on the basis of the cations to which they are connected. As shown in Figure 7a

Coordination Numbers
The microstructure of the Na 2 O-MgO-Al 2 O 3 -SiO 2 -K 2 O-CaO slag system was mainly composed of the silicoaluminate and the network modifier. The coordination numbers (CN) varied with the fraction of WS, as shown in Figure 6. It could be seen that the coordination number of Si was maintained at about four, but decreased with increasing WS content, indicating a reduction in the strength of Si-O bond. We attributed this to the depolymerization of the highly coordinated Si by metal atoms, such as Mg and K. For Al, not only four-coordinate Al-O structures but also a large number of structures with fewer than four coordination bonds formed. With increasing WS content, the fraction of more highly coordinated Al-O structures increased. This indicates that the [AlO 4 ] tetrahedral unit is not as stable as the [SiO 4 ] tetrahedron.

Coordination Numbers
The microstructure of the Na2O-MgO-Al2O3-SiO2-K2O-CaO slag system was mainly composed of the silicoaluminate and the network modifier. The coordination numbers (CN) varied with the fraction of WS, as shown in Figure 6. It could be seen that the coordination number of Si was maintained at about four, but decreased with increasing WS content, indicating a reduction in the strength of Si-O bond. We attributed this to the depolymerization of the highly coordinated Si by metal atoms, such as Mg and K. For Al, not only four-coordinate Al-O structures but also a large number of structures with fewer than four coordination bonds formed. With increasing WS content, the fraction of more highly coordinated Al-O structures increased. This indicates that the [AlO4] tetrahedral unit is not as stable as the [SiO4] tetrahedron.

Distribution of Oxygen Types
The oxygen ions in this system can be classified as either BO or NBO on the basis of the cations to which they are connected. As shown in Figure 7a

Distribution of Oxygen Types
The oxygen ions in this system can be classified as either BO or NBO on the basis of the cations to which they are connected. As shown in Figure 7A Figure 7C. The most abundant NBOs were Al-O-K and Si-O-K, indicating that the tetrahedron is more easily destroyed by K 2 O than other oxides in this network structure. Zhang et al. [23] proposed that the disparity between

Distribution of Bond Angles
The distributions of bond angles for the different samples are shown in Figure 8.

Distribution of Bond Angles
The distributions of bond angles for the different samples are shown in Figure 8

Distribution of Q n
The degree of polymerization is an important parameter for characterizing the integrity of the melt network structure. This can be analyzed by counting the number (Q) of tetrahedra with n bridging oxygen atoms. The influence of the WS fraction on Q n for Si is shown in Figure 9. With increasing WS content up to 30%, Q 4 and Q 5 decreased, whereas Q 0 , Q 1 , Q 2 , and Q 3 increased. We attribute these trends to the depolymerization of the network structure, which results in the ash having good fluidity.

Conclusions
In conclusion, we investigated the structure of ash slags via MD simulations. The content of WS in the slag was varied from 0% to 100wt%. The simulations provided structural information, such as the RDF, MSD, coordination number, bond angle, distribution of oxygen types, and Q n . We draw the following conclusions: (1) By analyzing the RDFs of the ion-oxygen pair, we found that the WS had no obvious influence on the bond length of Si-O. The widths of the first peaks in other ion-oxygen pairs were larger than that of Si-O, indicating a weaker bond strength than Si-O. The coordination number of Al was influenced by the addition of WS, demonstrating that WS influences the stability of the network structure. Moreover, the [AlO4] tetrahedron was found to be less stable than the [SiO4] tetrahedron.

Distribution of Q n
The degree of polymerization is an important parameter for characterizing the integrity of the melt network structure. This can be analyzed by counting the number (Q) of tetrahedra with n bridging oxygen atoms. The influence of the WS fraction on Q n for Si is shown in Figure 9. With increasing WS content up to 30%, Q 4 and Q 5 decreased, whereas Q 0 , Q 1 , Q 2 , and Q 3 increased. We attribute these trends to the depolymerization of the network structure, which results in the ash having good fluidity.

Distribution of Q n
The degree of polymerization is an important parameter for characterizing the integrity of the melt network structure. This can be analyzed by counting the number (Q) of tetrahedra with n bridging oxygen atoms. The influence of the WS fraction on Q n for Si is shown in Figure 9. With increasing WS content up to 30%, Q 4 and Q 5 decreased, whereas Q 0 , Q 1 , Q 2 , and Q 3 increased. We attribute these trends to the depolymerization of the network structure, which results in the ash having good fluidity.

Conclusions
In conclusion, we investigated the structure of ash slags via MD simulations. The content of WS in the slag was varied from 0% to 100wt%. The simulations provided structural information, such as the RDF, MSD, coordination number, bond angle, distribution of oxygen types, and Q n . We draw the following conclusions: (1) By analyzing the RDFs of the ion-oxygen pair, we found that the WS had no obvious influence on the bond length of Si-O. The widths of the first peaks in other ion-oxygen pairs were larger than that of Si-O, indicating a weaker bond strength than Si-O. The coordination number of Al was influenced by the addition of WS, demonstrating that WS influences the stability of the network structure. Moreover, the [AlO4] tetrahedron was found to be less stable than the [SiO4] tetrahedron.

Conclusions
In conclusion, we investigated the structure of ash slags via MD simulations. The content of WS in the slag was varied from 0% to 100 wt%. The simulations provided structural information, such as the RDF, MSD, coordination number, bond angle, distribution of oxygen types, and Q n . We draw the following conclusions: (1) By analyzing the RDFs of the ion-oxygen pair, we found that the WS had no obvious influence on the bond length of Si-O. The widths of the first peaks in other ion-oxygen pairs were larger than that of Si-O, indicating a weaker bond strength than Si-O. The coordination number of Al was influenced by the addition of WS, demonstrating that WS influences the stability of the network structure. Moreover, the [AlO 4 ] tetrahedron was found to be less stable than the [SiO 4 ] tetrahedron. (2) The concentration of BO decreased and that of NBO increased with increasing WS content. This indicates that the network structure can be destroyed by metal ions, resulting in the conversion of BO to NBO.
(3) The O-Si-O bond angle was constant at approximately 107.5 • . The O-Al-O bond angles were wide and varied with increasing WS content. (4) With increasing WS content, the diffusion coefficients of Si 4+ varied only slightly with no obvious trend. The sum of all diffusion coefficients increased with WS content, indicating an increased diffusion capacity of the system. (5) For a WS content of 30%, the Q 4 and Q 5 structural units transformed into Q 3 , Q 2 , Q 1 , and Q 0 structural units to the greatest extent, which corresponded to the lowest degree of polymerization and highest ash fluidity.