Abstract
In the process of coal-gas outburst, the gas-solid two-phase flow of pulverized coal and gas will induce large-scale damage to underground mining. In this work, in order to study the dynamic evolution law of gas-solid two-phase flow of coal-gas outburst, and clarify the short-time destructive effect of outburst dynamic phenomenon, an EDEM-FLUENT coupled model is constructed to realize the numerical simulation of coal-gas outburst two-phase flow. The simulation results show that the flow velocity of pulverized coal in a short time (tens of milliseconds) approaches to the maximum velocity rapidly, which can reach 60 m/s. The velocity of pulverized coal is inversely proportional to the particle size. The initial acceleration of particles is slow, and then the acceleration increases rapidly. The accumulation angle of pulverized coal in the outburst hole is about 21° and is far less than the natural accumulation angle of pulverized coal, which is consistent with the general phenomenon of coal-gas outbursts. The results also show that the shock wave overpressure of the coal-gas outburst of numerical simulation is similar to that of the physical simulation test, and the difference is less than 10%. And the pressure change in the roadway is completed in short time. In terms of spatial relationships, the pressure at the position closest to the outburst mouth is not the maximum, and the maximum pressure exists at a certain position away from the outburst mouth.
Introduction
Coal-gas outburst is an extremely complex coal gas dynamic phenomenon (Du et al., 2020; Si et al. 2021a; Wang and Du, 2020). When coal-gas outburst occurs, a large number of pulverized coal and gas flow into the roadway or stope instantaneously, forming the gas-solid two-phase flow of pulverized coal and gas in the roadway or stope (Alexeev et al., 2004; Black, 2019; Huang et al., 2020; Wang et al., 2021). This kind of special gas-solid two-phase flow has the nature of a storm, which can go against the wind flow and fill thousands of meters long roadway. Moreover, it could destroy underground facilities and lead to casualties. When it is serious, it can cause gas explosion and coal dust explosion, leading to larger-scale damage (Lama and Bodziony, 1998; Liu and Harpalani, 2013; Guo et al., 2020; Sheng and Derek, 2016; Wang et al., 2013, 2020; Wang and Guo, 2019). In recent years, with the increase of coal mining depth in China, the outburst risk of coal seams is increasing. Many shallow non-outburst coal seams gradually transform into outburst coal seams after the mining activities enter the deep level (Du and Wang, 2019; Liu et al., 2021; Si et al., 2021b). Coal-gas outburst accidents have become the most important factor seriously threatening the safety production of coal resource (An et al., 2019; Du et al., 2018; Wang and Du, 2019; Xin et al., 2021; Zhao et al., 2018).
The movement process of outburst coal gas two-phase flow involves complex and high dynamic solid-gas coupled effect. Therefore, for the high-pressure pulverized coal gas two-phase flow with impact effect formed after outburst, the combination of physical experiment and numerical simulation is one of the important means to study coal-gas outburst two-phase flow. Litwiniszyn (1985, 1995) and Topolnicki (1995) first simulated the process of coal-gas outburst through similar experiments. They measured the change of gas pressure in the process of outburst, and analyzed the energy transformation and energy balance in the process of outburst. Zhu (2010) and Jin (2017) have also studied the morphology, formation mechanism and migration law of coal-gas outburst two-phase flow through the outburst simulated experimental device, and analyzed the propagation attenuation law and influencing factors in the roadway. In terms of numerical simulation, Otuonye and Sheng (1994) used MAC Cormack explicit finite difference method to study the shock wave and gas migration law in the roadway after coal-gas outburst. The simulation results show that the propagation speed of shock wave is ahead of gas flow. Li et al. (1991) established the mathematical equation of the movement of coal and rock mass thrown out in roadway space during coal-gas outburst. They solved it with basic language, and obtained the throw distance of coal samples with different particle sizes under the conditions of different coal seam dip angle, different outburst hole inclination angle and different gas/coal particle initial velocity. Kang (2011) studied in detail the variation law of gas state parameters in the process of gas combustion, explosion and discharge from two aspects of experiment and numerical simulation. Sun et al. (2012) assumed that the total pressure in the hole during the outburst process was constant. In their work, the Mixture multi-phase flow model was used to study the flow field distribution and axial velocity of outburst fluid after the coal gas two-phase flow is ejected. It is found that the instantaneous expansion of a large number of incomplete expansion gas flow at the outburst mouth in the roadway space can produce huge dynamic effects, which seriously damages the underground production equipment and ventilation facilities. Wang et al. (2011, 2012) and Zhou et al. (2014, 2020) built an experimental system to study the propagation law of coal-gas outburst shock wave in vertical roadway, and analyzed the causes of outburst shock wave. On this basis, the geometric model of straight roadway is established, and the damage effect of shock wave is numerically simulated. Meanwhile, the distribution laws of shock pressure, velocity and gas relative mass concentration in roadway were obtained. Taking the coal-gas outburst accident in Daping Coal Mine as an example, Yang (2009) demonstrated that the variation law of gas concentration after outburst conforms to the law of power exponent attenuation. Based on this, they carried out numerical simulation research and obtained the change law of gas outburst and abnormal emission. Xue et al. (2011, 2014) simulated the distribution of pore pressure, principal stress and yield zone in the initial stage of gas outburst in coal roadway driving. They analyzed the influence of gas pressure, buried depth and coal mechanical strength on outburst, and studied the migration and propagation characteristics of gas shock wave after outburst. Xu et al. (2006) and An et al. (2013) have also enriched the technical means of coal-gas outburst research from different aspects using computer numerical simulation method.
It can be seen that most of the current literature on numerical simulation of coal-gas outburst mainly focus on the change of gas outburst shock wave, or the change of stress field and strain field of coal and rock in outburst area. However, there are few reports on numerical simulation of outburst coal gas two-phase flow. Therefore, in this work, the EDEM-FLUENT coupled numerical model of outburst coal gas two-phase flow was established, and the dynamic characteristics of outburst coal-gas two-phase flow were simulated and analyzed, which has certain guiding significance for the mechanism of coal-gas outburst.
Coupled principle of EDEM-FLUENT
Coupled simulation process of EDEM-FLUENT
In order to study the dynamic characteristics of coal gas two-phase flow in coal-gas outburst, the coupled numerical modeling of discrete element software EDEM and finite element software FLUENT is selected to analyze the movement laws of outburst shock wave and pulverized coal. In the coupled simulation, the software FLUENT is mainly responsible for calculating the law of outburst gas and the pressure propagation change of the shock gas flow, while the EDEM software is mainly responsible for calculating the law of movement of the outburst pulverized coal in the roadway. Meanwhile, the data interaction between the two software is realized through the coupled interface.
Based on the operating principle of DEM discrete element software and CFD finite element software, it can be known that neither of them can complete the simulation and analysis of dense phase solid-gas two-phase flow alone. It is very difficult to develop a software that has the performance of two types of software at the same time, but these two types of software can be linked by establishing the DEM-CFD coupled framework. And the numerical simulation of solid-gas two-phase flow could be completed based on their respective computational advantages. The typical EDEM-FLUENT coupled solution process is shown in Figure 1.

Simulation of solid-gas two-phase flow by EDEM-FLUENT coupled.
According to Figure 1, the operating logic of the EDEM-FLUENT coupled operation is that after the simulation starts, a flow field calculation of a time step is first performed by FLUENT, and then a time step calculation is performed by EDEM according to the initial conditions. Then, the flow field data of this step are obtained from FLUENT, and the effect of fluid on particles is calculated according to the flow field data. Meanwhile, the motion law of particles is solved. Moreover, the position of the particles from EDEM is obtained by FLUENT and the effect of the particles on the fluid is solved. Then the cycle calculation of the next time step is performed. It can be found that the calculation of EDEM and FLUENT is not at the same time, and there is a time step misalignment in operation. However, if the time step is set short enough, the time misalignment effect of a time step can be ignored.
There are two calculation methods for gas-solid two-phase flow, namely Euler-Lagrangian approach and Euler-Euler approach. The application scenario of Euler-Lagrangian approach is generally adopted for two-phase flow with less solid phase, which is characterized by fast calculation speed, but large calculation error for dense solid phase. As for the Euler-Euler approach, in the dense solid two-phase flow model, more mathematical methods of solid-phase and solid-phase, solid-phase and gas-phase are added in the Euler-Euler approach. Euler-Euler approach is more accurate, but it takes up a lot of computer resources. In view of the two-phase flow characteristics of coal-gas outburst, the Euler-Euler approach is selected in this simulation work.
Particle contact model
The interaction and feedback analysis between the solid phases in the numerical simulation of dense solid-gas two-phase flow are the core of the discrete element software calculation. Various particle contact models are actually collision analysis between solid phases, and the solid phase forces and moments calculated by different contact models are different. According to the different application scenarios, different contact models need to be adopted. The commonly used contact models in EDEM software include Hertz-Mindlin contact model, Hertz-Mindlin bonded contact model, Linear Cohesion contact model, Moving Plane contact model, etc. Among them, the Hertz-Mindlin contact model is the default model of the software. This contact model can also be called “elastic-damping-friction contact mechanics model” (Gu and Yang, 2019; Wu and Feng, 2008), as shown in Figure 2. This model has been successfully applied in a variety of dense-phase solid-gas two-phase flow EDEM simulation analysis. Therefore, Hertz-Mindlin contact model is used to analyze the movement law of pulverized coal group in this paper.

Contact model of Hertz-Mindlin.
The positive contact force Fcni between two pulverized coal particles is expressed as:
In the formula, E′ represents the equivalent elastic modulus of the particle. R′ represents the equivalent particle radius, and α represents the degree of collision. The expression of elastic modulus and equivalent particle radius is:
In the formulas, E, v, R indicates the elastic modulus, Poisson's ratio and radius of pulverized coal, respectively. Subscripts 1 and 2 represent the codes of two pulverized coal particles.
The expression of the normal damping force Fdni between two pulverized coal particles is:
In the formula, m' represents the equivalent mass of pulverized coal. β is the coefficient. Sn represents the normal stiffness of the coal. vnrel represents the normal velocity of the relative movement of two pulverized coal particles. The expressions are:
In the formulas, m1 and m2 represent the mass of the two pulverized coal particles respectively. e represents the coefficient of restitution during the collision. Sn represents the normal overlap amount. v1 and v2 respectively represent the respective speeds of the pulverized coal before the collision. n represents the relative speed direction at the time of collision.
In the formula, r1 and r2 represent the center position vector of the coal particles.
Particle tangential contact force Fcti is as follows,
In the formula, St represents the tangential stiffness of the coal powder. δ represents the tangential overlap of the two particles at the moment of collision.
The tangential stiffness St is expressed as:
In the formula, G′ represents the equivalent shear modulus of coal, and its expression is:
Where G represents the shear modulus of coal.
The tangential damping force when pulverized coal contacts is Fdti, which could be expressed as,
In the formula, vtrel indicates the relative tangential velocity when coal pulverized coal collides, which is calculated as:
The setting of the rolling friction coefficient of pulverized coal in the numerical simulation is also very important, which can be expressed by the torque of the friction surface:
In the formula, μr represents the rolling friction coefficient of the pulverized coal on the simulated roadway. Ri represents the distance between the center of the pulverized coal and the collision point. wi represents the unit angular velocity vector of the pulverized coal at the collision point.
The establishment of numerical model of outburst pulverized coal-gas two-phase flow
The geometric dimensions of the numerical model are based on the physical simulation test system developed in previous work (Wang and Wu, 2017; Wang K and Wang L, 2019; Wu et al., 2020), and the full-scale model of the physical simulation test system is carried out to have good comparative data. The total length of the model is 9.41 m, and the outburst cavity is a cylinder with a diameter of 0.2 m and a length of 0.3 m. In order to avoid the small structure existing in the modeling to increase the burden on the mesh generation and numerical solution, the pressure sensors in the simulated roadway are not created, as shown in Figures 3 and 4. During the movement of coal-gas two-phase flow, parameters such as pulverized coal velocity and shock wave pressure can be obtained by setting numerical monitoring points in the model. The geometric model is established and the mesh is generated using ANSYS ICEM. The grid files are respectively substituted into the software of FLUENT and EDEM. Then the simulation parameters are set, and the coupled solver is used to solve the problem.

Overall schematic diagram of the system.

Schematic diagram near the outburst cavity.
In the EDEM software, the particle parameters should be as real as possible, but the particle size of the pulverized coal used in the physical simulation test in our previous work is 10∼35 mesh, 35∼80 mesh, 80∼200 mesh, 200 mesh or smaller, which is too small for numerical simulation (Wang, 2020). Taking the largest particle size group 10∼35 mesh as an example, the corresponding size is about 1 mm. If the maximum accumulation theory of pulverized coal is used for accumulation and filling in the outburst cavity, the number of pulverized coal in the outburst cavity will reach ten million levels. The consumption of computer resources by the number of tens of millions of particles in EDEM far exceeds the computing and storage limits of personal computers or ordinary workstations. Therefore, in this numerical simulation, the pulverized coal is enlarged into three types of size, with equivalent diameters of 4 mm, 5 mm, and 6 mm. These three diameters were enlarged by 10 times correspond to 40 mesh, 32 mesh, and 28 mesh, respectively. In order to ensure as much as possible that the pulverized coal after the size enlargement is similar to the corresponding real-size pulverized coal movement, the pulverized coal density in the software is modified based on the similarity theory. According to the force analysis of coal particles, there are mainly two kinds of forces, namely the horizontal drag force and the vertical gravity. The relationship between gravity and particle size is cubic, and the relationship between drag force and particle size is square. As the particle size increases ten times, in order to ensure that the ratio of gravity to drag force in the simulation remains unchanged, the density of pulverized coal should be set to one-tenth of that of real pulverized coal. So, the density of coal in the numerical simulation is set as 145 kg/m3.
According to our previous work (Wang et al., 2019), pulverized coal rarely undergoes further crushing during the outburst process. Therefore, the force crushing effect of the pulverized coal during the outburst process is not considered in the numerical simulation.
The porosity of the pulverized coal used in the physical test after natural compaction in the cavity is about 50%, that is, the volume of pulverized coal and the volume of air are each half of the cavity (Wang K, Wang L et al., 2019). The parameter setting of software refers to this porosity. The total volume of the cavity is 10,52,40,000 mm3, and the volumes of pulverized coal with diameters of 4 mm, 5 mm, 6 mm are 268.1 mm3, 523.6 mm3, and 904.8 mm3 respectively. Meanwhile, the total volume of these three particle diameters of coal is required to be equal in the assumption, so the numbers of these three particles in the model are 6543, 3350, and 1939 respectively, with a total of 11,832 granular coal powders, as shown in Figure 5. In the numerical simulation analysis of this work, the initial gas pressure is set to 0.75 MPa. Once the simulation is started, the pulverized coal particles move rapidly until the end when the energy is depleted.

Pulverized coal setting in outburst cavity.
Simulation result analysis
Movement of pulverized coal particles
The move process of pulverized coal in the outburst roadway system at different times after the outburst occurs is shown in Figure 6. The color of the particles in the figure represents the movement speed of the pulverized coal. Due to the randomness of the velocity parameters of a single particle, the foremost 10 particles of each particle size in the pulverized coal flow are selected, and their average velocity is calculated, as shown in Figure 6.

Movement of pulverized coal at different times. (a) 0 s. (b) 0.02 s. (c) 0.04 s. (d) 0.06 s. (e) 0.08 s. (f) 0.1 s.
As shown in Figure 7, the average speeds of pulverized coal with diameters of 4 mm, 5 mm, and 6 mm were rapidly accelerated in the first 0.04 seconds and approached the highest speed, reaching 41.2 m/s, 48.4 m/s and 60.9 m/s, respectively. As the figure shows the movement of the front-most pulverized coal particles, they are all in the front end, and the solid phase is very sparse, so there is no collision between them to form mutual obstruction or acceleration. As a result, the velocity and particle size show a good inverse proportional relationship.

Comparison of movement velocity of three kinds of particle size coal.
The numerical simulation migration speed of coal particles with the diameter of 5 mm is compared with the results calculated in previous theoretical calculation work (Wang et al., 2020), as shown in Figure 8. The theoretical calculation and numerical simulation results are not exactly the same of necessity. In order to better compare the movement trends of the particles, the two curves in Figure 8 are placed at the same position at the beginning and last point (The actual operation is to multiply the theoretical calculation result by a coefficient of 0.924). The difference between theoretical calculation and numerical simulation is represented by the difference in curve trend. It can be seen from the figure that in the initial stage, the speed of the pulverized coal particles calculated theoretically is faster, and then its acceleration slows down, and the speed is caught up by the numerical simulation. The reason is that there is an acceleration process in the gas flow in the numerical simulation, and the drag force provided to the pulverized coal particles is relatively small when the numerical simulation starts.

Comparison of pulverized coal velocity in theoretical calculation and numerical simulation.
Distribution law of outburst pulverized coal in roadway
If the pulverized coal moves freely, their moving distance will be very large, but the laboratory space is limited. Therefore, a collecting box is set at the end of the simulated roadway. The collecting box can block the further movement of coal powder and collect the coal powder moving to the distance. Figure 9 shows how the pulverized coal moves and accumulates in the collecting box during the movement of the outburst coal-gas two-phase flow (The color in the figure indicates the speed of the pulverized coal particles). And Figure 10 shows the accumulation of pulverized coal near the outburst cavity after the numerical simulation is completed (The color in the figure indicates the particle size of the pulverized coal particles: red is the diameter of 6 mm, green is the diameter of 5 mm, and blue is the diameter of 4 mm). The pulverized coal accumulation angle in the outburst cavity in Figure 10 is about 21°, which is much smaller than the natural pulverized coal accumulation angle, conforming to the general law of coal-gas outbursts.

Moving coal in the collecting box.

Accumulation of pulverized coal after the outburst.
The accumulation of pulverized coal in every area of the roadway could be easily counted using the built-in function of the software EDEM. As shown in Figure 11, the numerical simulation results are compared with the physical simulation results of previous work (Wang et al., 2020). According to the subsection of physical simulation, the results of numerical simulation are divided into different areas. From left to right, there are No.1 to No.4 subsections of simulated roadway and the collecting box. It can be seen from the figure that the accumulation characteristics of pulverized coal in the roadway in the numerical simulation and the physical simulation are similar. The ratio of the total outburst coal amount of the numerical simulation and the physical simulation to the total coal load is 61.5% and 55.4%, respectively. In view of the fact that the numerical simulation setting involves the debugging of a large number of parameters in the two software of EDEM and FLUENT, and considering the complexity of the solid-gas interaction itself, the fact that numerical simulation can realize the analysis of the particle motion law is encouraging. Thus, the difference between 61.5% and 55.4% is not big under the circumstances.

Statistical comparison between physical simulation and numerical simulation of pulverized coal accumulation.
The particle size of the pulverized coal at the farthest end of the model is used to analyze the outburst sorting characteristics of pulverized coal. Figure 12 shows the mass distribution of the three particle sizes in the collecting box. It can be seen from the figure that the pulverized coal has almost no sorting characteristics under the numerical simulation. This is because that the numerical simulation parameter settings are similar to the experimental conditions. The shock wave lasts for a short time, so there is no period of uniform motion. Taking pulverized coal with the diameter of 4 mm as an example, the deceleration distance after the speed reaches 30 m/s is 17.4 m, which is far beyond the total length of the simulated roadway of 8 m. In this case, it is difficult to form a sorting phenomenon.

Mass distribution of pulverized coal with different particle size in collecting box.
Outbuest shock wave pressure
Figure 13 shows the pressure monitoring points set to collect data at 2.27 m, 4.27 m, 6.27 m and 8.27 m, which are in the same position as the pressure sensor of the physical simulation (Wang K, Wang L et al., 2019). The 4 curves in the figure represent the change process of the outburst pressure data of the four monitoring points. It should be noted that the collected value is static pressure, and the pressure sensor in this physical simulation is the total pressure (stagnation pressure). These two values can be converted into each other through aerodynamic theory. The calculation results show that there is little difference between the numerical simulation and the experimental data, and the difference is within 10%.

Pressure change in roadway.
In the numerical simulation, the pressure change in the roadway is completed in an instant. The pressure at the 1# monitoring point increased slowly after the sudden change. The pressure at 2# sensor position is the maximum value, and then gradually decreases, indicating that the shock wave pressure in the numerical simulation has a tendency of increasing first and then decreasing.
Conclusion
In order to study the dynamic evolution law of gas-solid two-phase flow of coal-gas outburst, and clarify the short-time destructive effect of outburst dynamic phenomenon, an EDEM-FLUENT coupled model is constructed to realize the numerical simulation of coal-gas outburst two-phase flow. The main conclusions are as follows.
The numerical simulation of coal-gas outburst two-phase flow was carried out, and the movement law of pulverized coal was analyzed. The velocity of pulverized coal can reach 60 m/s in a short time (tens of milliseconds). The velocity of pulverized coal is inversely proportional to the particle size. Comparing the difference between the theoretical calculation and the numerical simulation, it can be found that the theoretical calculation of particle acceleration is faster in the initial stage, and then the acceleration slows down. And the velocity is overtaken by the velocity result of numerical simulation. However, the overall movement trend is consistent. The accumulation angle of pulverized coal in the outburst hole is about 21°, and far less than the natural accumulation angle of pulverized coal, conforming to the general law of coal-gas outburst. The pressure change of shock wave in the simulated roadway is completed in an instant. The pressure at the position nearest to the outburst mouth is not the maximum. The maximum pressure appears at a certain position from the outburst mouth and has the characteristics of a sudden increase.
Article highlights
An EDEM-FLUENT coupled model was constructed for coal-gas outburst two-phase flow. The dynamic evolution law of gas-solid two-phase flow of coal-gas outburst was simulated and analyzed. The short-time destructive effect of outburst dynamic phenomenon was studied and clarified.
Footnotes
Author contributions
Feng Du: Formal analysis, Writing-original draft. Yangyang Guo: Investigation, Writing-review & editing, Supervision. Liang Wang: Methodology, Investigation, Writing-review & editing, Supervision. Chao Xu: Formal analysis. Aitao Zhou: Formal analysis. Gongda Wang: Conceptualization, Formal analysis.
Data availability statement
The data used to support the findings of this study are available from the corresponding author upon request.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was funded by the National Natural Science Foundation of China (52004291, 51874314, 51974161), the State Key Laboratory Cultivation Base for Gas Geology and Gas Control (Henan Polytechnic University) (WS2020A04), the Key Laboratory of Mine Disaster Prevention and Control (Shandong University of Science and Technology) (JMDPC202104), the Fundamental Research Funds for the Central Universities (2021YJSAQ25).
