*4.2. Mathematical Experiment Modeling*

#### 4.2.1. Physical Model and Meshing

According to the orthogonal tests, it is necessary to establish three physical models with different roadway lengths. The spatial relationship model and grid meshing of the ventilation system are shown in Figure 4. The parameters of each test model are shown in Table 2.

**Figure 4.** Physical model and its meshing.



## 4.2.2. Mathematical Model

There are three basic equations of the fluid numerical simulation, which includes the mass conservation equation, the momentum conservation equation, and the energy conservation equation. After the gas outburst accident in the excavation roadway, the airflow should also be satisfied by gas migration. Furthermore, a chemical component exchange exists in this numerical simulation, i.e., CH4, O2, and N2, so the component transport equation should be introduced.

A turbulent flow is a common flow phenomenon in nature, and the fluid is always in the state of turbulence in most engineering problems. The realizable *k-ε* model is adopted because its good performance for solving the adverse pressure gradient and vortex problems.

## 4.2.3. Boundary Conditions

The boundary conditions of the model are shown in Table 3.

From Table 3, the physical model area is the Fluid Type, and the pressure of operation condition is 101.325 kPa. The inlet boundary condition includes Inlet 1 of the intake airway, Inlet 2 of the return airway, outburst position, and auxiliary fan. The inlet boundary type of auxiliary fan inlet is the Velocity-inlet, and the other is the Pressure-inlet. The species of the outburst position is 100% CH4, and the other is 20% oxygen and 80% nitrogen.

The outlet boundary includes Outlet 1 of the intake airway and Outlet 2 of the return airway, the outlet boundary type is also the Pressure-outlet. Furthermore, the species are 20% oxygen and 80% nitrogen either.

The gauge pressure of each boundary condition is calculated as follows: Assuming that both the sectional area of the return airway and excavation roadway is 25 m2, the friction resistance coefficient is 100 × <sup>10</sup>−4, the inlet air quantity of the intake airway and return airway *Q*<sup>10</sup> = *Q*<sup>20</sup> is 75 m3/s, and the air quantity of the auxiliary fan *Qfan* is 4 m3/s, 10 m3/s, and 20 m3/s, respectively. The outlet of the return airway is *Q*<sup>21</sup> = *Q*<sup>20</sup> + *Qfan*, *Q*<sup>12</sup> = *Q*<sup>10</sup> − *Qfan*.

**Table 3.** Boundary conditions of each test.


The relationship between the air resistance and the friction resistance coefficient of the roadway is

$$R = \frac{LIL\alpha}{S^3} \tag{1}$$

where *R* represents the air resistance, kg/m7; *L* represents the length of the roadway, m; *U* represents the perimeter of the roadway section, m; *α* represents coefficient of frictional resistance, kg/m3; *S* represents the sectional area of the roadway, m2.

From Equation (1), the air resistance value of the return airway is *R*<sup>10</sup> = *R*<sup>12</sup> = *R*<sup>20</sup> = *R*<sup>21</sup> = 0.0032 kg/m7, the excavation roadway from the air duct outlet to the return airway is *R*<sup>31</sup> = 0.0256 kg/m7. The air pressure *H*<sup>10</sup> = 16.13 Pa, *H*<sup>12</sup> = *H*<sup>20</sup> = 18 Pa, *H*<sup>21</sup> = 19.97 Pa and the air pressure of the excavation roadway are calculated by the air quantity of the auxiliary fan and the length of the excavation roadway.

Furthermore, the initial gauge total pressure at Inlet 1 is the sum of the pressure difference of the air door and the air pressure of the upper wind side of the intake airway *H*10, and the initial pressure of Outlet 1 is the sum of the pressure difference of air door and *H*12. Suppose the operation pressure is the standard atmospheric pressure (101.325 kPa), the gauge total pressure at Outlet 2 of the return airway is 0 Pa. Then, the air pressure of each inlet and outlet are calculated as shown in Table 3.

#### 4.2.4. Controlling Parameters of Fluid Dynamics

(1) Solver

In the Pressure-based type, the absolute velocity formulation and Pressure-velocity coupling algorithm are usually adopted to solve the problem. Furthermore, using the SIMPLE scheme to calculate the mathematical model, it computes the mass conservation and obtains the pressure field by the mutual correction of pressure and velocity.

(2) Convergence accuracy

A reasonable accuracy is an important parameter to ensure the convergence of the model. The convergence residuals are 10−<sup>6</sup> in this model, and the number of iterations is 500, so the calculation will be finished when the residuals of each variable are less than 10−<sup>6</sup> or the iterations reach to 500 steps. Furthermore, the convergence of the calculation results can be dynamically monitored by checking the iterative residual of each variable.

#### *4.3. Results of Numerical Simulation*

Using Fluent numerical simulation software, the flow field distributions of the outburst gas of the above nine tests are simulated. Figure 5 illustrates the path-line near the air door of each test on the horizontal plane 1 m away from the floor, which is colored by the velocity magnitude.

**Figure 5.** Color map of velocity path-line near the air door.

From Figure 5, the air velocity of L3 is the highest near the air door after the gas burst accident, followed by L2 and L1, and L7 is the smallest. The excavation roadway length of L1, L2, and L3 is 200 m, which has the greatest air velocity magnitude. The larger air velocity has a greater impact on the ventilation system, especially in the air inlet part. Figures 6 and 7 show the air velocity distributions of the two inlet boundaries.

**Figure 6.** Air velocity of Inlet 2 of the return airway.

**Figure 7.** Air velocity of Inlet 1 of the intake airway.

As can be seen from the contour map in Figures 6 and 7, all inlet boundaries are reversed in each test, and the air velocity of Inlet 2 is larger than Inlet 1. Therefore, the influence of the gas outburst accident on the return air system is greater than that on the inlet air system. Extracting the maximum air velocity of the section, the airflow reversal degree is calculated, as shown in Table 4. Meanwhile, the variation curve of air velocity and airflow reversal degree is shown in Figure 8.

**Table 4.** Air velocity and airflow reversal degree of the inlet boundary.


**Figure 8.** Variation of air velocity and airflow reversal degree.

From Table 4 and Figure 8, the airflow reversal degree of Inlet 2 is larger than Inlet 1, that is, because the air pressure distribution of the ventilation system can restrain the outburst gas, and the air pressure variation of the intake airway is smaller than the return airway.

The numerical simulation results show that the airflow reversal degree of 200 m roadway is greater than that of 1000 m roadway and 2000 m roadway. Therefore, the airflow reversal degree decreases with the increase of the roadway length. Furthermore, the outburst air flow pressure of L3, L6, and L9 is 0.3 MPa, and the absolute value of their airflow reversal degree is greater than the outburst pressure of 0.1 MPa and 0.2 MPa in the same length of the excavation roadway. Therefore, the airflow reversal degree increases with the increase of the outburst pressure and the outburst energy.

#### **5. Comprehensive Evaluation of Airflow Disorder**

#### *5.1. Fuzzy Comprehensive Optimization Theory*

According to the main factors that affect the safety and stability of the ventilation system as well as the basic data of orthogonal experiments, the air quantity of the auxiliary fan, the air pressure difference of the air door, and the air flow reversal degree of Inlet 1 and Inlet 2 are used to establish the evaluation model of the airflow disorder on the basis of variable fuzzy theory.

Assuming F = {*Fij*}, *i* = 1,. ... , *m*, *j* = 1,. ... , *n* indicates that the parameter of the *m* test model corresponds to the set of evaluation factors of *n.* Here, *m* = 4. The method of relative membership is adopted to deal with the initial data. There are three situations,

(a) The larger the better

$$Fr\_{ij} = \frac{F\_{ij} - \min\_{1 \le i \le m} \{ F\_{ij} \}}{\max\_{1 \le i \le m} \{ F\_{ij} \} - \min\_{1 \le i \le m} \{ F\_{ij} \}} \tag{2}$$

where *Frij* is the relative membership of the set *Fij* and its range is [0, 1]; max {} indicates the maximum value of the set; min {} indicates the minimum value of the set.

(b) The smaller the better

$$Fr\_{ij} = \frac{\max\_{1 \le i \le m} \{ F\_{ij} \} - F\_{ij}}{\max\_{1 \le i \le m} \{ F\_{ij} \} - \min\_{1 \le i \le m} \{ F\_{ij} \}} \tag{3}$$

(c) The value is equal to 1

$$Fr\_{i\rangle} = 1, \left(F\_{1\rangle} = F\_{2\rangle} = \dots \dots = F\_{mj}\right) \tag{4}$$

The fuzzy partition matrix is defined as follows,

$$\mathcal{U} = \begin{bmatrix} \mu\_{11} & \mu\_{12} & \dots & \mu\_{1m} \\ \mu\_{21} & \mu\_{22} & \dots & \mu\_{2m} \end{bmatrix} = \begin{bmatrix} \mu\_{1i} \\ \mu\_{2i} \end{bmatrix} \tag{5}$$

where *<sup>u</sup>*1*<sup>i</sup>* is subordinate to the superior index and *<sup>u</sup>*2*<sup>i</sup>* is the inferior index, <sup>2</sup> ∑ *k*=1 *uki* = 1, *uki* ∈ [0, 1].

The objective function is the minimal sum of the weight difference squared of the *m* models.

$$\min f(u\_{1i}) = \sum\_{i=1}^{m} \left( \left( u\_{1i} \sqrt{\sum\_{j=1}^{4} \left[ \omega\_j (gr\_j - Fr\_{ij}) \right]^2} \right)^2 + \left( u\_{2i} \sqrt{\sum\_{j=1}^{4} \left[ \omega\_j \left( Fr\_{ij} - br\_j \right) \right]^2} \right)^2 \right) \tag{6}$$

where *ω<sup>j</sup>* is the weight of each factor; *grj* represents the standard superior membership, which is the optimal value of evaluation factor *j*; *brj* represents the standard inferior subordinate, which is the worst value of evaluation factor *j*.

Supposing *d f*(*u*1*i*) *du*1*<sup>i</sup>* = 0, the optimal fuzzy partition matrix is calculated. Then, the influence of different factors is evaluated by *u*1*i*.

#### *5.2. Quantification and Weight Analysis of the Factors*

#### 5.2.1. Factors

The pressure difference of the air door has a strong inhibition of the outburst energy, which belongs to situation (a). The greater the pressure difference of the air door, the greater the pressure energy of the air inlet roadway. The air quantity of the auxiliary fan belongs to situation (a). Because it is opposite to the direction of the outburst gas, the resistance to the migration of the outburst gas increases with the larger air quantity of the auxiliary fan. The airflow reversal degree of Inlet 1 and Inlet 2 also belongs to situation (a), the larger the better.

The above factors consist of the *F* matrix.

#### 5.2.2. Weight

Since the airflow reversal degree is more important than the pressure difference of air door and the air quantity of the auxiliary fan, the weight of the airflow reversal degree is 0.4 and 0.4, respectively, at Inlet 1 and Inlet 2. The weight of the pressure difference of the air door is equal to the air quantity of the auxiliary fan, their value is 0.1, respectively. Then the weight vector *ω<sup>j</sup>* = (0.1, 0.1, 0.4, 0.4)T.

#### *5.3. Results and Discussion*

The fuzzy comprehensive optimization theory is used to evaluate the airflow disorder degree of the mine ventilation network after a gas outburst accident. From Tables 1 and 4, the max {*Fij*} = (1500, 1200, −19.50, −7.22)T, min {*Fij*} = (500, 240, −85.23, −46.77)T. The membership vector of nine groups is calculated by Equations (2)–(6), *u*1*<sup>i</sup>* = (0.4928, 2.2160, 23.5568, 0.1586, 0.0394, 0.3398, 0, 0.0746, 0.1209). The values of each factor and evaluation results are shown in Figure 9.

**Figure 9.** The evaluation results of each test.

From Figure 9, the evaluation result of L3 is the largest, illustrating that this parameter condition has the strongest influence on the ventilation system, and L7 is the smallest. In L3, the length of excavation roadway is 200 m, the outburst pressure is 0.3 MPa, the pressure difference of air door is 1500 Pa, and the air quantity of the auxiliary fan is 600 m3/min. In L7, the length of the excavation roadway is 2000 m, the outburst pressure is 0.1 MPa, the pressure difference of air door is 1500 Pa, and the air quantity of the auxiliary fan is 1200 m3/min.

According to the extreme difference of the range analysis method, the results of the orthogonal experiment are determined. The extreme difference is the difference between the maximum and minimum value of the airflow disorder of each factor. With the increase of the influence range, the disturbance degree of the factor increases. The result is shown in Table 5 and Figure 10. *Ki* represents the sum of levels of factor *i*; *ki* represents the average value of the levels of factor; *R* represents the polar difference.


**Table 5.** Range analysis of the orthogonal experiment.

**Figure 10.** Range analysis of the orthogonal experiment.

From Figure 10, A > B > D > C, the length of the excavation roadway and the outburst pressure are relatively important factors affecting the ventilation system. The airflow reversal degree of the ventilation system increases with the increase of the outburst pressure or the decrease of the length of the excavation roadway.

The energy is continuously attenuated while the outburst gas migrates in the roadway. Adopting the method of reducing the air resistance of the airway or changing the air pressure distribution of the main fans, the initial air pressure in both the roadway and its parallel roadway is easily affected and will be increased in a gas outburst accident. Meanwhile, it will contribute to increase the total air quantity and improve the kinetic energy of the airflow under the condition of the constant capacity of the main fans, so it is harder to reverse the airflow.

#### **6. Conclusions**

This paper analyzes the reversal principle of airflow induced by a coal and gas outburst in the excavation roadway and proposes the indexes of the airflow disorder. Through an orthogonal experiment, it is found that all the inlet boundaries are reversed in each test, and the air velocity of Inlet 2 is larger than Inlet 1, indicating that the influence of the gas outburst accident on the return air system is greater than that on the inlet air system. Moreover, the larger air velocity has a greater impact on the ventilation system, especially the air inlet part.

On the basis of the fuzzy comprehensive optimization theory, the evaluation model of the airflow disorder is established. The results show that the order of effect factors for the stability of the mine ventilation system are the length of excavation roadway > outburst pressure > pressure difference of air door > air quantity of auxiliary fan. The most influential result is that the length of the excavation roadway is 200 m, the outburst pressure is 0.3 MPa, the pressure difference of the air door is 1500 Pa, and the air quantity of the auxiliary fan is 600 m3/min. The airflow reversal degree of the ventilation system increases with the increase of the outburst pressure or decreases with the length of the excavation roadway.

**Author Contributions:** Methodology, J.S.; software, Y.W.; formal analysis, J.C.; investigation, L.L.; resources, J.S., W.H., L.L. and T.L.; data curation, J.S.; writing—original draft preparation, L.L.; writing—review and editing, J.S., J.C. and Z.L.; supervision, J.C. and Y.W.; project administration, J.S.; funding acquisition, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by "the Fundamental Research Funds for the National Natural Science Foundation of China" (51804120, 52074122), "the National Key Research Project of China" (2018YFC08080306), "the Fundamental Research Funds for the Central Universities" (3142018003) and Science and Technology Research Project of Higher Education in Hebei Province (Z2018004).

**Acknowledgments:** The authors would like to thank the anonymous referees for their thoughtful comments and suggested edits that have improved the rigor and presentation of this work.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**

