A large amount of dust is produced in the process of coal mining and transportation, which seriously affects the safe production of enterprises and the health of employees. In recent years, the problem of dust management has been paid more and more attention. (Luo et al. 2017; Reed et al. 2018, 2019; Shahan and Reed 2019; Yao et al. 2020; Tan et al. 2020). Dust reduction by wetting as the most commonly used method is far more effective than ventilation. Coal seam water injection is a fundamental and efficient method to prevent and control dust from its source (Cheng et al. 2012). In terms of mining and tunneling, a variety of hydrating dust collectors have been popularized (Wang et al. 2017, 2014).

The wettability of coal is closely related to its metamorphic degree, chemical structure, industrial components, bedding structure, pore characteristics and liquid properties (Meng et al. 2019; Wang et al. 2019; Xu et al. 2017). A lot of research results for the wettability of coal dust are characterized by the contact angle—a macroscopic property, leaving the microscopic property somewhat undiscovered. Thus, experimental results always can be inaccurate.

Molecular dynamics (MD) is a simulation method that enables atoms or molecules to interact with each other in a system with given conditions (time, temperature, pressure and other macro thermal conditions). The dynamic evolution process can be presented and various properties of the system can be calculated by MD. MD can reveal the law of material evolution from the micro system of atomic level, which can be used as a micro explanation (Sharma 2019) and it is also a tool to verify theories. The reliability of modern molecular dynamics simulation methods, combined with the First-Principle (or ab initio calculations) and other quantum methods can be considered noteworthy (Glukhova 2017).

Yet, molecular dynamics simulation, an effective research method on micro systems, has some difficulties in the practice of coal systems with complex chemical structure and an irregular arrangement. At present, it is a common practice to choose simulation researches of other materials and apply their parameters to coal related systems, so there is no recognized optimal set of simulation methods and parameters for coal. In the past, H₂O molecular model with a fixed boundary and a droplet shape with molecules of magnitude of only a few thousand have been established to study the simulation and contact angle calculations for solid materials. Qiang et al. and Zhang et al. (2018, 2020) doubt the practical significance and rationality of this approach of the H₂O system’s macroscopic shape at such a small micro scale. Still, this approach is a common practice in MD simulations to investigate the wettability of solid materials. As a consequence, there is barely research using the microscopic molecular distribution of H₂O on the surface of coal dust to directly reflect the wettability of coal dust, which is considered more reasonable.

In this paper, through the molecular dynamics simulation of coal dust systems, at different scales wetted under different simulation conditions and modeling methods, the distribution of micro scale wettability on the coal surface is revealed. Furthermore, the difference and reliability of results under various simulation conditions are revealed. Ultimately, more reliable simulation and modeling methods, procedures and parameters are obtained, which is a preliminary attempt to change the status quo of over simplifying the coal system and empiricism of the molecular dynamics simulation methods of other materials.

2.1 Modeling and optimization

In this paper, the Australia’s Yalluourn coal structure studied by Kumagai is selected as the basic model of coal (Xia 2017). Firstly, a single coal and H₂O molecular model were sketched in Materials Studio. The structure of coal is shown in Fig. 1.

Structure of Yalluourn brown coal molecule

In order to improve the accuracy and efficiency of the simulation, geometric optimization, namely energy minimization, is needed for the established coal and H₂O molecular models. For the initially established molecular model with an imperfect initial structure, the CASTEP module was used to run a geometry optimization task to minimize energy (Imai et al. 2001). The molecular configuration after the geometric optimization is shown in Fig. 2 and these changes include the atomic position, bond length and bonding patterns of the coal molecular.

Molecular structure after CASTEP geometric optimization respectively a Lignite; b H₂O

After the optimization of a single coal and H₂O molecule, the molecular layers were established by Amorphous Cell module. The coal molecular number was temporarily set as 10 and the target density as 0.7 g/cm³ which is an approximate general density of lignite. The geometric optimization was carried out by DMol3 to obtain an ideal initial coal molecular layer structure (Xu et al. 2007). The layer of coal after establishment is shown in Fig. 3a. The energy minimized layer is shown in Fig. 3b. Compared with the same position of the coal molecular layer in Fig. 3a, the optimized layer in Fig. 3b is more compact, which means that the optimization result is well.

Building and optimizing lignite system(small) a Lignite system(small) built by Amorphous Cell; b Lignite system(small) after DMol3 geometric optimization

Then, since the research object is solid–liquid, the position and velocity of coal as a solid should not change with the simulation. This process can be completed by setting Constrain command for coal molecule. The validity of this approach will be verified in Sect. 3.2. At the same time, the Amorphous Cell module is adopted to establish a H₂O molecule layer with a number of 100, as shown in Fig. 4a. Also, H₂O is liquid phase, whose target density is set as 1.0 g/cm³. Geometric optimization was carried out to get the ideal initial H₂O molecular layer structure. Subsequently, the energy of the molecular layer is minimized by the CASTEP module too, which is shown in Fig. 4b. It can be seen that the optimized molecular layer becomes sparser and has lower energy by analyzing this figures.

Building and optimizing H₂O system(small) a H₂O molecular system (small) built by Amorphous Cell; b H₂O molecular system (small) after CASTEP geometric optimization

Establishing molecular layer respectively and their preliminary structure optimization, this paper adopted Build Layers command for combining the two systems into one. Both options of layer building were attempted, that is, 2D periodic surface structure and 3D periodic crystal structure as provided by Build Layers, shown in Fig. 5a and b. For the 2D periodic surface structures, the areas were 452.098 Ų, 750.690 Ų and 955.567 Ų for small, medium and large-scale Coal-H₂O Molecular Systems respectively; and for 3D periodic crystal structures, the volume of the simulation box was 14,378.358 ų. The results of the two modeling methods are discussed in Sect. 3.1. The established layers and surface structures were also minimized with DMol3, as shown in Fig. 5c. It shows that partially discrete coal molecular converge and H₂O molecular become uniform and loose distribution after optimization. For coal, it is a solid substance whose density is higher, so the coal molecular tends to shrink and becomes denser. In addition, H₂O is a fluid and low-density substance, so it becomes loose. Therefore, the optimization result is reasonable.

The Coal-H₂O Molecular System(small) built by Build Layers and optimization a 2D periodic surface system; b 3D periodic lattice system; c The optimized 3D periodic lattice system

2.2 Force field and ensemble selection

The MD simulation of the Coal dust-H₂O molecular system studied in this paper was based on the Forcite module in Materials Studio. In this paper, the ab initio force field compass was used to calculate the forces of other substances with the same element for replacing the experiment-based potential energy function (Savin and Mazo 2020). All the calculations were carried out under the COMPASS force field, which is suitable for organic or some inorganic molecules or generally used to calculate the properties of materials (Liu et al. 2020). The Coal dust-H₂O molecular system studied in this paper was an open system with a constant temperature 285 K and the number of particles was assumed to be constant. The canonical ensemble NVT was therefore selected for the molecular dynamics simulation in this paper.

2.3 Parameters setting

For the surface structure, in the dynamics settings of Forcite, the initial velocity was set to random allocation. The Anderson method was used to provide a constant temperature. Main parameters are shown in Table 1. In addition to, the charge was automatically distributed; for the charge addition method, the Ewald method was selected to calculate the Coulomb interaction (Abdulnour Y. Toukmaji 1996); whilst the atom-based method was selected to calculate the van der Waals interaction.

2.4 Pretreatment of simulation

Prior to running the molecular dynamics simulation, anneal was applied to fully relax the system (Yan and Luo 2014). Specifically, five annealing cycles were run under the NVT ensemble, with the initial temperature of 298 K, heated to 500 K and then cooled to 298 K. The number of steps was 5000; the time step was 1 fs; and the simulation time was 5 ps. Through annealing, the total energy of the system was even smaller and reached a more stable state—this was a more appropriate structure for any molecular dynamics to start with. The output configuration is shown in Fig. 6. After the step of annealing, compression-decompression treatments were carried out on crystal Coal dust-H₂O Molecular System in NPT by Forcite dynamics, with the pressure set as 1, 0.8, 0.6, 0.4, 0.2 and 0.1 MPa respectively by the same simulation parameters as annealing. The influence of this treatment as a variable on crystal system simulation will be discussed in Sect. 3.1. The output configurations of each compression-decompression are shown in Fig. 7 (the above treatment cannot be carried out on the surface systems). Not only can this treatment gradually reach the pressure value of approximately real atmospheric pressure, but it also helped the system reach an even more stable state on the basis of annealing, which made the dynamic simulation results more reliable. After annealing and compression-decompression treatments, the molecular dynamics simulation of the relatively stable crystal Coal dust-H₂O Molecular System was carried out with the same basic dynamics parameters as mentioned above.

Coal-H₂O Molecular System after geometric optimization and anneal (small; crystal)

Compression-decompression treatment on 3D periodic lattice system(small). a 1 MPa; b 0.8 MPa; c 0.6 MPa; d 0.4 MPa; e 0.2 MPa; f 0.1 MPa

3.1 Small-scale coal dust-H₂O molecular system under various conditions

The configurations are shown in Fig. 8a-d. At the same time, the files generated by Frame can be used to view the displacement of the H₂O molecule system throughout the simulation process. It can be seen from the relative concentration analysis graph (Fig. 8e) that in the selected direction (0, 1, 0), and H₂O molecules were distributed in the range of 17 Å to 48 Å from the bottom to the top, while the lignite molecules’ distribution ranged from 4 to 36 Å. The spatial distribution of the coal was irregular, which was determined by the model construction and geometric optimization. Its molecules were set in fixed positions and did not move with the dynamics, so there was no further analysis. The molecular layer of lignite and H₂O molecules began to coincide at approximately the 36 Å plane which was the initial position of lignite on its + Z side and continued to 17 Å in the interior of lignite system, accounting for 59.375% of the Z axis length of the coal dust in this system.

Forcite MD simulation results of small-scale 2D periodic surface Coal-H₂O Molecular System at 285 K and NVT. a Step 25,000; b Step 50,000; c Step 75,000; d Step 100,000; e Relative concentration by Concentration Profile Analysis

In the same system, the temperature rises to 295 K and other conditions remain unchanged. The H₂O and lignite molecules were respectively distributed in the range of 17 Å to 50 Å and 2.8 Å to 35.3 Å and the H₂O molecules were concentrated in the range of 25 Å to 44 Å. The molecular layer of lignite and H₂O molecules began to coincide at 35.3 Å plane and continued to 16.8 Å in the lignite system. According to Fig. 9, little difference was shown for the internal permeability and surface distribution on lignite molecules compared with 285 K, so the simulation of other different temperatures was discontinued.

Forcite MD simulation results of small-scale 2D periodic surface Coal-H₂O Molecular System at 293 K and NVT. a Dynamics output configuration; b Relative concentration by concentration profile analysis

Although in the initial selection of a reasonable simulation ensemble, NVE was excluded considering the practical environment, its irrationality has not been proven yet. Therefore, the author has carried out the NVE dynamics simulation again under the same other conditions. Figure 10 shows that the diffusion and infiltration of H₂O molecules in this ensemble are stronger than those under the NVT ensemble with a trend of diffusion around the coal dust and the accumulation on the contact surface between the coal and H₂O molecular layers is significantly less than that of the NVT simulation.

MD simulation results of small-scale 2D periodic surface Coal-H₂O molecular system (285 K, NVE). a Dynamics output configuration; b Relative concentration analysis

With other conditions remaining constant, the simulation output configuration of the three-dimensional periodic cell system without pressurization is shown in Fig. 11a and the H₂O molecules are disordered and highly diffused. The simulation output configuration after annealing and NPT pressurization of the same model and its relative concentration analysis are shown in Fig.11b, c. From Fig.11a, b, it is evident that the H₂O molecules are more constrained as compared to those without annealing pressurization. According to Fig. 10c, the H₂O and lignite molecules were distributed in the range of 20 Å to 65.5 Å and 3.5 Å to 30.5 Å. The H₂O molecules were relatively concentrated in the range of 28 Å to 55 Å. The two layers began to coincide at the plane of about 30.5 Å and continued to 20 Å in the interior of the lignite system, accounting for about 38.89% of the Z axis length of coal dust in this system. From 30.5 Å to 55 Å on the surface of lignite molecular system, the number of H₂O molecules reached its peak part, accounting for 90.741% of the Z axis length of the coal dust system. The distribution of H₂O molecules in the area above 55 Å gradually decreased to zero. It was apparent that the infiltration of H₂O molecules into the interior of lignite was weaker than that of the surface system and the Z axis distribution range of lignite surface increased, however, there is no accumulation sign with droplet flow. Consequently, 3D periodic cell systems are not appropriate for simulating coal dust wetting of the Coal dust-H₂O Molecular Systems in this paper, though NPT pretreatment can indeed be an optimization approach to some extent.

Forcite MD simulation results of small-scale 3D periodic crystal Coal-H₂O molecular system (285 K, NVT). a Without compression-decompression treatment; b With compression-decompression treatment; c Relative concentration analysis (compression-decompression treated)

This paper considers that the possible reason for this difference is that the three-dimensional cell system adopts the periodic boundary conditions in X, Y and Z directions. After the establishment of the Coal dust-H₂O Molecular System as a 3D periodic cell, a vacuum layer with a certain thickness was attached between the molecular layers of lignite and H₂O by default. This was operated as "lignite molecular layer—vacuum layer—H₂O molecular layer—vacuum layer" by Forcite in the Z axis direction. This type of perpetually periodic structure is technically known as images of the lignite and H₂O molecular layer above and below the simulation box and when H₂O molecules run out of the box, there are corresponding images that run into the box to maintain the constant particle number. In order to study the interaction of the two layers of materials, the cut-off radius cannot be set too small, which can be a solution in other circumstances. Consequently, this influence cannot be avoided under 3D periodic boundary conditions. It does not have a significant influence on the XOY plane in this study; on the contrary, it can reflect the properties of a huge macro system into the established small simulation box, which greatly saves in terms of computing resources. However, it can interfere with the interaction of H₂O molecules and coal dust in the Z axis direction. The surface model with 2D periodic boundary has no periodicity in the Z axis direction and the practice of this paper has verified that the interference of mirror image is avoided as a result, of which the comparison of Fig. 8a and Fig. 11a is an example.

In this paper, small systems are simulated prior to try various simulation conditions and gradually explore the superior simulation conditions such as ensemble, thermodynamic parameters, modeling methods which are suitable for the studied substance, systems and processes. The results show that 2D periodic surface structure perform well and it is necessary to take sufficient NPT pretreatment into consideration if a 3D periodic crystal structure is adapted. Besides, temperature variation was proven to be irrelevant under the simulation conditions of this paper. Thus, the most reliable results of medium and large systems can be obtained in the shortest time based on the research of small-scale coal—H₂O molecular system under various conditions.

3.2 Medium-scale coal dust-H₂O molecular system

The medium-scale Coal dust -H₂O Molecular System in this paper is consisted of 200 H₂O molecules and 20 lignite molecules. Firstly, the model was established according to the aforementioned modeling method and the fixed constraints of coal molecules were removed to carry out a full geometric optimization. Then, the NVE annealing, which was different from the previous NVT annealing, was carried out along with geometric optimization as shown in Fig. 12a. Finally, the first NVT dynamic simulation with the same parameters was carried out. Before the second NVT simulation, the original NVT annealing and fixing constraints of the coal layer were applied as before, with the same kinetic parameters as shown in Fig. 13a.

MD simulation results of medium-scale 2D periodic surface Coal-H₂O molecular system (285 K, NVT dynamics, annealed in NVE, coal layer unfixed). a After NVE anneal; b Dynamics output configuration; c Relative concentration analysis

MD simulation results of medium-scale 2D periodic surface Coal-H₂O Molecular System (285 K, NVT dynamics, annealed in NVT, coal layer fixed). a After NVT anneal; b Dynamics output configuration; c Relative concentration analysis

Figure 12b, c show the output structure image and relative concentration analysis after the simulations respectively. It is evident that the coal molecular layer displacement in the first simulation was extremely large. Due to the unconstrained Cartesian coordinates and fractional coordinates, the coal molecules were fully aggregated, forming a structure closer to the actual coal molecular state. Yet the simulation result image shows that the H₂O molecular diffusion was intense and did not show obvious contact characteristics with coal. It is speculated that the whole and local positions of the coal layer were constantly changing, making it difficult for H₂O molecular layer to fully interact with it. In addition, the annealing was carried out under NVE, an isolated system and the intense diffusion of H₂O molecules in consistency with the results of previous section, also affected the dynamic simulation that followed. As a result, NVT annealing and fixed coal molecular layer simulation are suitable to the focus of this research. The second simulation result, adopting the above methods, showed that with the expansion of the simulation system, the H₂O molecules' penetration into the coal dust and the outward diffusion effect of H₂O molecules throughout the system were both weakened. At the end of the simulation, as shown in Fig. 13b, the H₂O molecular layer was more concentrated on the coal layer surface of the larger system. One possible reason for this is the stronger attraction of the larger coal dust system to H₂O molecules. This result may better reflect the contact characteristics of H₂O molecules as a whole with coal dust. The concentration analysis (Fig. 13c) shows that the coal molecular layer was distributed in the range of 4.8–41.3 Å from the bottom to the top of the crystal cell, while the H₂O molecule was distributed in the range of 31.9–53.2 Å. The furthest H₂O molecule only penetrated to approximately 9 Å into the coal dust, accounting for 24.66% of the distribution range of the Z axis of the coal molecular layer, with a low relative concentration value. However, the H₂O molecules were concentrated in a small range of coal dust surface with a high relative concentration value, accounting for 32.60% of the Z axis distribution range of the coal molecular layer, demonstrating clear characteristics of the wetting process and concentration mutation which is a sign of a boundary.

3.3 Large-scale coal dust-H₂O molecular system

The large-scale Coal dust-H₂O Molecular System in this paper is consisted of a micro droplet system composed of 500 H₂O molecules and a coal dust system composed of 50 lignite molecules and is built as a 2D periodic surface model by the Build Layers command. Due to its large size, the lignite molecular layer was annealed and geometrically optimized with higher precision than the previously small and medium-sized systems before simulation. The temperature range was enlarged to 300–1000 K and the number of cycles was increased to 50. The output configuration is shown in Fig. 14a.

MD simulation results of large-scale 2D periodic surface Coal-H₂O Molecular System. a After NVT anneal; b Dynamics output configuration; c Relative concentration analysis

The configuration and relative concentration analysis were obtained as shown in Fig.14b, c after dynamic simulation. Instead of wide-scale diffusion of the H₂O layer, its stable displacement towards the coal molecular layer was obvious as expected. Thus, under this condition, the H₂O molecules had sharper aggregation and horizontal extension tendencies on the surface of the coal dust.

All of the Molecular Dynamic simulations in this paper have been completed. Typical relative concentration analysis results were compiled and sorted as shown Table 2. The differences in simulation results under different conditions are clearly presented.

Employing molecular dynamics simulations under various conditions for investigation, this paper establishes suitable model and determines dynamics simulation method in the wetting process. the results were compiled as followed:

1. (1)

   The simulated results obtained by annealing and calculation under NVT ensemble are superior to those obtained by NVE ensemble. The diffusion and infiltration of H₂O molecules in the NVT ensemble are weaker than that in the NVE ensemble, so it is more suitable for simulating the process of coal dust wetting.

2. (2)

   In modeling, the 2D periodic surface structure is more reasonable than the 3D periodic crystal structure. This model avoids the periodic mirror-image interference in the Z axis direction of the study and the results can better reflect the research focus of wettability. After the simulation, 3D particle distribution can also be done to analyze quantitatively.

3. (3)

   The wetting behavior of H₂O molecules can be better reflected in a large scale system of coal dust and H₂O molecules. The results of relative concentration analysis showed that the relative concentration range of H₂O molecules in the optimal simulation results for each system scale decreased with the increase of the system scale, which are 56.25%, 32.60% and 25%, respectively, indicating that H₂O molecules tended to adhere to the surface of coal dust in the large system.