# Multifield Coupled Dynamic Simulation of Coal Oxidation and Self-Heating in Longwall Coal Mine Gob

Authors: Shuanglin Song; Shugang Wang; Shuang Jiang; Yuntao Liang; Peiyu Hu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The spontaneous combustion of residual coal in coal mine gob has long been a problem and poses a threat to the safe production of coal. Therefore, it is of great significance to conduct an in-depth study of the oxidation and self-heating progress of residual coal in the gob. Considering that the geometric dimensions and physical characteristics of the gob will change during the advance of the working face, the purpose of the present paper is to determine how the coal self-heating develops during and after coal mining. A fully coupled transient model including gas flow, gas species transport, and heat transfer in the gob and the butt entries, as well as heat transfer in the surrounding strata, is developed to quantify the evolution of coal self-heating in gob during and after mining. The model was solved by COMSOL Multiphysics package and then verified by comparing the field data with the simulated data. On this basis, parametric studies including the influences of the surrounding strata temperature, airflow temperature, coal-rock particle size, and advance rate of the working face on coal oxidation and self-heating in the gob were conducted. The results show that a tailing phenomenon of the high-temperature area is formed on the air inlet side of the gob during mining, and the temperature in the high-temperature zone decreases gradually due to the accumulation and compaction of the gob and heat dissipation to the surrounding strata. Also, although the temperature in gob increases gradually after the stopping of mining, the high-temperature area migrates towards the working face. Moreover, when the temperature of the surrounding strata is consistent, different ventilation temperatures have no obvious effect on the maximum temperature of the gob at the initial mining stage, whereas the higher ventilation temperature results in the higher self-heating temperature after several days of mining. Finally, the smaller average particle size or faster advance rate results in a lower maximum temperature of gob.


---

## Body

## 1. Introduction

Spontaneous combustion of residual coal in gob is the most common type of fire in coal mines, accounting for 60% of underground fire accidents in China [1]. Spontaneous combustion of coal is a phenomenon that commonly occurs in coal mines, which results from the self-ignition of coal following its self-heating in a mine heading or its immediate surrounding [2]. Crushed coal left in the gob and air flowing through the gob with a certain speed, as well as high concentration of oxygen, are the primary determinants of the beginning of the self-heating process, which can begin the spontaneous combustion of coal [2, 3]. The spontaneous combustion of residual coal in gob is a dynamic development process, which is the result of the comprehensive action of many factors and has strong complexity, variability, and contingency [4, 5]. In the process of studying the flow, heat transfer, and mass transfer in gob, scholars analyzed the mechanism of the coal-oxygen combination theoretically and studied the distribution law of the airflow field, temperature field, and concentration field in the gob by using the numerical simulation method.In the research of mathematical models, Zhang and Wang [6, 7] established a mathematical model of gas seepage in a gob based on 2D porous medium steady flow, and the heat exchange with surrounding rock was ignored. Ding et al. [8] established a 2D nonlinear seepage equation for the gob and derived a solution for the nonlinear seepage equation. Zhang et al. [9] established a mathematical model of the unsteady temperature field in the gob based on the theory of convective heat and mass transfer. Jiang and Zhang [10] established a 3D flow field mathematical model based on the classic linear Darcy law. Du and Yang [11] built a 3D mathematical model for the unsteady temperature field in the gob using theories of fluid mechanics, heat transfer, and mass transfer. Li et al. [12–14] established a nonsteady mathematical model of coal spontaneous combustion to disclose the distribution laws of air leakage, oxygen concentration, and temperature in the gob based on the air leakage continuity equation, oxygen concentration field equation, and temperature field equation. At the same time, based on mine pressure observation, Li et al. [12] established the finite element seepage mathematical model of gob permeability coefficient distribution by describing the expansion coefficient of collapsed rock. Xu [15] established an unsteady model of the temperature field in the gob to calculate the relationship between the temperature distribution and the advance rate of the working face. Deng et al. [16] established an oxygen concentration distribution model based on heat transfer and conservation of composition mass. Shan et al. [17] established a mathematical model of seepage in the gob based on dispersion dynamics and seepage theories and conducted a related analysis with the finite element algorithm. Hua [18] established a gas fuzzy seepage mathematical equation of the gob based on the irregularity of the gob and the fuzzy mathematical theory. In recent years, with the development of numerical simulation technology, related scholars have established mathematical models for coal spontaneous combustion in the gob. Tan et al. [19] established a one-dimensional model for spontaneous combustion of coal in the gob under moving coordinates through introducing a moving coordinate system to process variations in physical boundaries, which achieved a steady-state transformation of the temperature field. Xia et al. [20] established a thermo-fluid-solid-chemical coupling model to simulate the deformation of coal in the process of coal oxidation and the relationship between permeability and coal temperature. Liu [21], Zhu et al. [22] and Zhang [5] established a coupling model for spontaneous combustion of coal in the gob, which not only realized the dynamic simulation of coal oxidation and self-heating, but also simplified the calculation process. Taraba and Michalec [23] developed a three-dimensional, single-phase model with a continuously advancing longwall face and studied the effect of the longwall face advancing speed on the oxidation heat production as well as the evolution of the gases in the gob area.In the research of factors affecting the spontaneous combustion of coal in the gob, Brodny and Tutak [24] determined the impact of ventilation airflow volume on the location of the special hazard zone with endogenous fire in the gob and pointed out that the ventilation rate of the air had a significant impact on the location and the range of the special hazard spontaneous fire zone. Wen [25] thought that coal spontaneous combustion in the gob was mainly determined by the thickness of residual coal, air leakage intensity, oxygen concentration, advance rate of working face, and spontaneous combustion period. Then, the risk of spontaneous combustion of coal in the gob under a different advance rate of working face and air leakage intensity was predicted [26]. In the model established by Ejlali et al. [27] and Lohrer et al. [28], the effect of moisture on the heat transfer process was considered. Yang et al. [29] theoretically analyzed the influence of changes in air volume and pressure on the spontaneous combustion status of coal in the gob and divided the range of the oxidation zone in the gob by numerical simulation under different conditions. Wang et al. [30] established a dynamic model of the porosity in the gob with the variation of position. Li et al. [31] explained the abrupt change of temperature rise and variability of factors of coal spontaneous combustion based on the carving coal-rock “O” ring theory of the gob and the nonuniform distribution characteristics of coal-rock and oxygen consumption. Tutak and Brodny [3] studied the influence of the tensile strength of roof rocks on the extent of the zone with a particularly high risk of spontaneous coal combustion in caving goaves of the longwalls ventilated with the Y-type system. Yuan and Smith [32] studied the effects of coal properties (coal surface area and heat of reaction) on spontaneous heating in longwall gob (mined-out) areas by using FLUENT software.Since spontaneous combustion of coal in the gob is a multifield coupling process involving airflow, stress, fracture, oxygen concentration, and temperature, some scholars have studied the evolution of coal spontaneous combustion in high-temperature areas with the multifield coupling research method [21, 33–36]. Meanwhile, there are a variety of unsteady factors that can affect the process of spontaneous combustion of coal, such as the air pressure disturbance caused by the variation of the ventilation system, the enlarged physical dimension of the seepage field near working face, and the changes of physical properties of porous skeletons affected by the overlying strata and compaction in the gob, which requires that the physical model established in numerical simulation research be as close to the site conditions as possible. Therefore, the dynamic grid technology [37] of CFD software was employed by some scholars to obtain the heterogeneous distribution of porosity in time and space, and the simulation result was more in line with the actual situation. Besides, some scholars combined COMSOL Multiphysics® with MATLAB to generate a new mine gob grid and took the field information of the last time step as the initial value of the current time step. Thus, the continuity of simulated physical quantity was achieved in time and space, which reflects the heating process of coal spontaneous combustion in gob to a large extent [4, 38, 39].However, the above studies did not consider the influence of factors such as the surrounding strata temperature and ventilation temperature on the spontaneous combustion of residual coal in the gob. Besides, the evolution of the high-temperature area in the mined-out gob still needs to be discussed when stopping mining, for which the self-heating process has lasted for several days since the start of the mining activity.Considering that the gob environment is always dynamically changed in the mining process, the openings of the butt entries as well as the surrounding strata were taken as the boundary to jointly solve the physical fields of these different areas in this study. In this way, the dynamic evolution of coal oxidation and self-heating in the gob was studied during and after the coal mining, which can not only improve the prediction level of coal spontaneous combustion but also provide a new technical means for delineating the location of hidden fire sources in the gob.


## 2. Problem Description

Longwall mining is an exploitation method used in flat-lying, relatively thin tabular coal seams in which a long face is established to extract the coal. The length of longwall working face in China is generally 80 m–250 m, and there are 300 m long working faces in recent years. The simplest roadway system of longwall working face includes open-off cut, the butt entries (including working face, section headentry, and section tailentry), and the gob. The coal mining activity starts from the open-off cut. One side of the open-off cut is the coal wall to be mined, and coal is dropped, loaded, and transported along the coal wall of the working face. With the development of coal mining, gob is gradually formed. With the continuous advance of the working face, the size of the mined-out area continues to expand.The layout of a typical longwall section with “U type” ventilation system is shown in Figure1. The working face advances with a rate of uwork, resulting in a decreasing area of the entity coal and an increasing area of the gob. The airflow rate in the entries can be controlled by using bleeders and regulators so that the methane accumulation is kept away from coal mining. Because of the void in the gob, the air flows through the gob with a certain speed and carries high concentration of oxygen into the gob as well, resulting in the coal oxidation process, which generates heat in this reaction. The heat can be continuously produced if air supply is not eliminated. The generated heat is transported out of the gob by conduction through coal mass, convection through flow, and radiation. If the rate of heat generation is greater than that dissipated to the external environment, the excessive heat will result in temperature rise in the gob.Figure 1
2D schematic diagram of longwall section with “U type” ventilation system.Under the dual effects of gravity and mining stress, the overlying strata continue to fall into the mined-out area, forming a porous medium filling space with certain distribution laws. According to the distribution of mining stress, the gob can be divided into three parts: coal pillar supporting influence zone, the stress transition area, and the recompaction area [40]. Different areas show different stacking patterns, with significant heterogeneity. In the area close to the working face and butt entries, because of the support of coal pillar, the gob is not pressed by the cantilever beam and overlying rock. The collapsed coal and rock mass are in the state of free accumulation. With good permeability and large air leakage flow, the basic conditions of coal spontaneous combustion are easy to meet. Therefore, these areas are also high-risk areas of coal oxidation and even spontaneous combustion. Because of the large mining intensity of the longwall mining method, the gob formed in the underground is also large, and the residual coal in the gob is easy to ignite naturally.


## 3. Multifield Coupled Dynamic Model

### 3.1. Assumptions

(1)
The length of the strike and dip of the gob is much larger than the carving height of the overlying strata caused by coal seam mining, so the variations of the physical parameters in the height of the gob are ignored; i.e., a 2D model is assumed and implemented(2)
Physical parameters of the air and coal gangue are not varied with temperature(3)
The effect of water vapor in the process of coal oxidation and self-heating is ignored


### 3.2. Dynamic Evolution of Porosity and Permeability in the Gob

The dynamic evolution of the longwall section during the mining process is shown in Figure1. With the advance of the working face, rock and residual coal in the free carving zone are gradually compacted under the effect of gravity and mining stress. The bulking factor is continuously decreased in the compaction process, and its value can be expressed by [31](1)Kpx,y=Kp,min+Kp,max−Kp,minexp−axdx1−exp−ξaydy,ξ1<1,where Kp (x, y) is the bulking factor (distribution function) of the gob, dimensionless; Kp, max is the initial caving and crushing expansion coefficient, dimensionless; Kp, min is the crushing expansion coefficient after caving compaction, dimensionless; αy and αx are the attenuation rates at the strike and dip direction of the crushing expansion coefficient in the gob, respectively (m−1), which are related to the physical and mechanical properties of the caving coal and rock, the geological characteristics of strata, the load of mining stress, and the caving duration (the advance of the working face); dy and dx are the distances from any point (x, y) in the gob to a fixed surrounding solid boundary and the working face boundary, respectively (m); ξ is the adjustment coefficient for controlling the distribution morphology of the “O” ring model (the larger the ξ is, the more likely the quadrate will be presented).If the advance rate of the working face isuwork, the distance from any point in the gob (x, y) to a fixed surrounding solid boundary of the gob and the working face over time can be expressed as(2)dx=x−x0+uworkt,dy=0.5L−y−0.5L.The height of the caving zone and the porosity can be expressed as [41](3)H=MKpx,yKp,max−1,ε=1−1Kpx,y,where H is the caving height (m); ε is the porosity; and M is the mining height of the working face (m).According to Ergun’s experimental research [42], the gob permeability can be calculated by the porosity and volume average particle size:(4)k=ε3dp21501−ε2,where k is the permeability of the gob (m2); and dp is the particle diameter of the porous medium of the gob (m).


### 3.3. Control Equation of Coal Oxidation and Self-Heating in the Gob

Assuming that the gas in the gob is incompressible, the continuity equation of gases and fluids can be written as(5)ρ∇⋅u=Qm,where u is the Darcy velocity field in the gob (m·s−1); and Qm is the source (sink) term of the gas (kg·m−3·s−1).The permeability in the gob has obvious heterogeneous characteristics. The airflow is in a turbulent state when its velocity is large in the gob near the working face, while the airflow is in a laminar or peristaltic state when its velocity is low inside the gob. Therefore, the Darcy–Brinkman equation was used to describe the airflow in the gob [43].(6)ρεu⋅∇uε=∇⋅−pI+Κ−μk−1+βερu+Qmε2u,where K=μ/ε∇u+∇uT−2/3∇⋅uI; p is the gas pressure in the gob (Pa); μ is the dynamic viscosity of the mixed gas in the gob (Pa⋅s); u is the magnitude of the velocity (m·s−1), for the 2D model u=u2+v2; ρg is the density of the mixed gas in the gob (kg·m−3); β is defined as(7)β=cFk,where cF is the gas drag coefficient and defined as(8)cF=1.75150ε3.The gas concentration equation in the gob can be expressed as [44](9)ε∂Ci∂t+u⋅∇Ci=∇⋅εDe,i∇Ci+εri,where Ci is the substance concentration (mol⋅m−3); De, i is the effective diffusion coefficient (m2·s−1), defined as ε4/3Di; and ri is the reaction rate expression for the species (mol·m−3·s−1). The oxygen concentration is the most important component among a variety of gas species in the gob, which determines the rate of coal oxidation reaction. The chemical reaction between coal and oxygen can be simplified as [32](10)Coal+O2⟶CO2+0.1CO+heatIn this paper, only oxygen (O2) and carbon monoxide (CO) are calculated in the model. The reaction rate of CO can be determined from the stoichiometric relationship between O2 and CO in reaction (10), and the O2 consumption rate is related to the oxygen concentration and coal oxidation parameters, which can be expressed by the Arrhenius equation:(11)rO2=−H1HCO2nAe−E/RT,where rO2 is the consumption rate of O2 (mol·m−3·s−1); CO2 is the oxygen concentration (mol·m−3); A is the preexponential factor (s−1); E is the activation energy (J·mol−1); R is the gas constant (J·mol−1·K−1); n is the reaction order, which is unity in this study; and H1 is the thickness of the residual coal in gob (m).Based on the conservation of energy, energy equations of air and coal-rock can be written as(12)ερcpf∂Tf∂t+ρcpfu⋅∇Tf=∇⋅ελf∇Tf+61−εhsfdpTs−Tf,(13)1−ερcs∂Ts∂t=∇⋅1−ελs∇Ts+1−εqs−q0,u+∇F−61−εhsfdpTs−Tf,where the subscripts s and f represent solid and fluid, respectively; T is the temperature (K); c and cp are solid specific heat and fluid constant pressure specific heat, respectively (J·kg−1·K−1); t is time (s); h is convective heat transfer coefficient (W·m2·K−1); λ is the thermal conductivity (W·m−1·K−1); F is the radiant flux vector (W·m−2); qs is the solid heat generation term (W·m−3). As it is related to the coal-oxygen reaction process, qs can be calculated with [45](14)qs=−rO2ΔH,where ΔH is the amount of thermal generation during the oxidation reaction of coal (J·mol−1), which can be obtained through coal sample experiments.The radiant fluxF should not be ignored when the local temperature of the coal-rock is high. It can be calculated in various methods. A simple approximation method can be used for the coal-rock; that is, radiative heat transfer is approximately processed into a diffusion process that is similar to heat conduction [45]. It can be calculated as follows:(15)F=βRT3∇T=λr∇T,where βR is constant (W·m−1·K−4); and λr is the radiation thermal conductivity (W·m−1·K−1).Considering that enough heat exchange can be conducted between air and coal-rock with the low seepage velocity in the coal-rock, the fluid temperature is consistent with the temperature of the coal-rock mass (orTf = Ts = T). By adding (12) and (13), the heat transfer equation can be expressed as [44](16)ρcm∂T∂t+ρcpfu⋅∇T=∇⋅λm∇T+1−εqs−q0+∇F,where the subscript m is the average volume, and ρcm=ερcpf+1−ερcs, λm=ελf+1−ελs.During the advance of the working face, the surface heat fluxq0 at the boundary of the top and bottom plates can be calculated by [13](17)q0=−uworkλsρcsπd2T−TW.


### 3.4. Control Equations of the Butt Entries

Since the airflow in the working face, headentry, and tailentry are generally turbulent in the coal mine, a turbulence model is required for simulation. Two-equation models, algebraic turbulence models, and large eddy simulation models are common turbulence models. The LVEL turbulence model [46] used in this paper is based on Prandtl’s mixed length theory, which calculates turbulent viscosity stably and efficiently based on local speed and distance from the nearest wall and is suitable for calculating internal flow. The gas composition and heat transfer equation are shown as follows:(18)∂Cb,i∂t+u⋅∇Cb,i=∇⋅Di∇Cb,i,ρcpf∂Tb∂t+ρcpu⋅∇Tb+∇⋅λf∇Tb=0,where Cb, i is the substance concentration (mol·m−3); Di is the gas diffusion coefficient (m2·s−1) in free flow; Tb is the air temperature in the butt entries (K). The gas species of O2 and CO are considered in equation (18), and no reactions are considered beyond the gob area.


### 3.5. Solid Heat Transfer Equation for Surrounding Strata

The heat transfer of the surrounding strata is primarily dominated by heat conduction, which can be expressed as(19)ρcs∂Ts∂t+∇⋅λs∇Ts=0.


### 3.6. Cross-Coupling Relations and the Simulation Environment

The established model has a coupling relationship in three aspects. First, the gob energy equation and the gas composition equation were coupled mutually inside the gob and were related to the velocity field solved by the gob flow equation. Secondly, the gas heat transfer equation and gas concentration equation of the butt entries were related to the internal velocity field inside them. Finally, the values of temperature, airflow velocity, pressure, and gas concentration at the interface between different computational domains were equal and continuous. The coupling relationship is shown in Figure2.Figure 2
Multiphysics coupling relationship of the model.The established model consists of several nonlinear partial differential equations in space and time and cannot be analytically solved. Therefore, a multiphysics modeling environment, the COMSOLMultiphysics®, is implemented. COMSOLMultiphysics employs the Finite Element Method and is capable of handling PDEs and DAEs. A segregated solution approach was chosen for the simulation, which splits the solution process into substeps.Besides, with the advance of the working face, the area of the gob will be continuously increased, and the physical size of the calculation area will be constantly varied (as shown in Figure3). Therefore, the calculation grid will be also adjusted accordingly. The deformation geometry model established by COMSOLMultiphysics realized the free deformation of the mesh by setting the moving speed of the calculation boundary, thereby simulating the dynamic growth of the gob.Figure 3
Change of computing domain during the dynamic advancement of the gob.


## 3.1. Assumptions

(1)
The length of the strike and dip of the gob is much larger than the carving height of the overlying strata caused by coal seam mining, so the variations of the physical parameters in the height of the gob are ignored; i.e., a 2D model is assumed and implemented(2)
Physical parameters of the air and coal gangue are not varied with temperature(3)
The effect of water vapor in the process of coal oxidation and self-heating is ignored


## 3.2. Dynamic Evolution of Porosity and Permeability in the Gob

The dynamic evolution of the longwall section during the mining process is shown in Figure1. With the advance of the working face, rock and residual coal in the free carving zone are gradually compacted under the effect of gravity and mining stress. The bulking factor is continuously decreased in the compaction process, and its value can be expressed by [31](1)Kpx,y=Kp,min+Kp,max−Kp,minexp−axdx1−exp−ξaydy,ξ1<1,where Kp (x, y) is the bulking factor (distribution function) of the gob, dimensionless; Kp, max is the initial caving and crushing expansion coefficient, dimensionless; Kp, min is the crushing expansion coefficient after caving compaction, dimensionless; αy and αx are the attenuation rates at the strike and dip direction of the crushing expansion coefficient in the gob, respectively (m−1), which are related to the physical and mechanical properties of the caving coal and rock, the geological characteristics of strata, the load of mining stress, and the caving duration (the advance of the working face); dy and dx are the distances from any point (x, y) in the gob to a fixed surrounding solid boundary and the working face boundary, respectively (m); ξ is the adjustment coefficient for controlling the distribution morphology of the “O” ring model (the larger the ξ is, the more likely the quadrate will be presented).If the advance rate of the working face isuwork, the distance from any point in the gob (x, y) to a fixed surrounding solid boundary of the gob and the working face over time can be expressed as(2)dx=x−x0+uworkt,dy=0.5L−y−0.5L.The height of the caving zone and the porosity can be expressed as [41](3)H=MKpx,yKp,max−1,ε=1−1Kpx,y,where H is the caving height (m); ε is the porosity; and M is the mining height of the working face (m).According to Ergun’s experimental research [42], the gob permeability can be calculated by the porosity and volume average particle size:(4)k=ε3dp21501−ε2,where k is the permeability of the gob (m2); and dp is the particle diameter of the porous medium of the gob (m).


## 3.3. Control Equation of Coal Oxidation and Self-Heating in the Gob

Assuming that the gas in the gob is incompressible, the continuity equation of gases and fluids can be written as(5)ρ∇⋅u=Qm,where u is the Darcy velocity field in the gob (m·s−1); and Qm is the source (sink) term of the gas (kg·m−3·s−1).The permeability in the gob has obvious heterogeneous characteristics. The airflow is in a turbulent state when its velocity is large in the gob near the working face, while the airflow is in a laminar or peristaltic state when its velocity is low inside the gob. Therefore, the Darcy–Brinkman equation was used to describe the airflow in the gob [43].(6)ρεu⋅∇uε=∇⋅−pI+Κ−μk−1+βερu+Qmε2u,where K=μ/ε∇u+∇uT−2/3∇⋅uI; p is the gas pressure in the gob (Pa); μ is the dynamic viscosity of the mixed gas in the gob (Pa⋅s); u is the magnitude of the velocity (m·s−1), for the 2D model u=u2+v2; ρg is the density of the mixed gas in the gob (kg·m−3); β is defined as(7)β=cFk,where cF is the gas drag coefficient and defined as(8)cF=1.75150ε3.The gas concentration equation in the gob can be expressed as [44](9)ε∂Ci∂t+u⋅∇Ci=∇⋅εDe,i∇Ci+εri,where Ci is the substance concentration (mol⋅m−3); De, i is the effective diffusion coefficient (m2·s−1), defined as ε4/3Di; and ri is the reaction rate expression for the species (mol·m−3·s−1). The oxygen concentration is the most important component among a variety of gas species in the gob, which determines the rate of coal oxidation reaction. The chemical reaction between coal and oxygen can be simplified as [32](10)Coal+O2⟶CO2+0.1CO+heatIn this paper, only oxygen (O2) and carbon monoxide (CO) are calculated in the model. The reaction rate of CO can be determined from the stoichiometric relationship between O2 and CO in reaction (10), and the O2 consumption rate is related to the oxygen concentration and coal oxidation parameters, which can be expressed by the Arrhenius equation:(11)rO2=−H1HCO2nAe−E/RT,where rO2 is the consumption rate of O2 (mol·m−3·s−1); CO2 is the oxygen concentration (mol·m−3); A is the preexponential factor (s−1); E is the activation energy (J·mol−1); R is the gas constant (J·mol−1·K−1); n is the reaction order, which is unity in this study; and H1 is the thickness of the residual coal in gob (m).Based on the conservation of energy, energy equations of air and coal-rock can be written as(12)ερcpf∂Tf∂t+ρcpfu⋅∇Tf=∇⋅ελf∇Tf+61−εhsfdpTs−Tf,(13)1−ερcs∂Ts∂t=∇⋅1−ελs∇Ts+1−εqs−q0,u+∇F−61−εhsfdpTs−Tf,where the subscripts s and f represent solid and fluid, respectively; T is the temperature (K); c and cp are solid specific heat and fluid constant pressure specific heat, respectively (J·kg−1·K−1); t is time (s); h is convective heat transfer coefficient (W·m2·K−1); λ is the thermal conductivity (W·m−1·K−1); F is the radiant flux vector (W·m−2); qs is the solid heat generation term (W·m−3). As it is related to the coal-oxygen reaction process, qs can be calculated with [45](14)qs=−rO2ΔH,where ΔH is the amount of thermal generation during the oxidation reaction of coal (J·mol−1), which can be obtained through coal sample experiments.The radiant fluxF should not be ignored when the local temperature of the coal-rock is high. It can be calculated in various methods. A simple approximation method can be used for the coal-rock; that is, radiative heat transfer is approximately processed into a diffusion process that is similar to heat conduction [45]. It can be calculated as follows:(15)F=βRT3∇T=λr∇T,where βR is constant (W·m−1·K−4); and λr is the radiation thermal conductivity (W·m−1·K−1).Considering that enough heat exchange can be conducted between air and coal-rock with the low seepage velocity in the coal-rock, the fluid temperature is consistent with the temperature of the coal-rock mass (orTf = Ts = T). By adding (12) and (13), the heat transfer equation can be expressed as [44](16)ρcm∂T∂t+ρcpfu⋅∇T=∇⋅λm∇T+1−εqs−q0+∇F,where the subscript m is the average volume, and ρcm=ερcpf+1−ερcs, λm=ελf+1−ελs.During the advance of the working face, the surface heat fluxq0 at the boundary of the top and bottom plates can be calculated by [13](17)q0=−uworkλsρcsπd2T−TW.


## 3.4. Control Equations of the Butt Entries

Since the airflow in the working face, headentry, and tailentry are generally turbulent in the coal mine, a turbulence model is required for simulation. Two-equation models, algebraic turbulence models, and large eddy simulation models are common turbulence models. The LVEL turbulence model [46] used in this paper is based on Prandtl’s mixed length theory, which calculates turbulent viscosity stably and efficiently based on local speed and distance from the nearest wall and is suitable for calculating internal flow. The gas composition and heat transfer equation are shown as follows:(18)∂Cb,i∂t+u⋅∇Cb,i=∇⋅Di∇Cb,i,ρcpf∂Tb∂t+ρcpu⋅∇Tb+∇⋅λf∇Tb=0,where Cb, i is the substance concentration (mol·m−3); Di is the gas diffusion coefficient (m2·s−1) in free flow; Tb is the air temperature in the butt entries (K). The gas species of O2 and CO are considered in equation (18), and no reactions are considered beyond the gob area.


## 3.5. Solid Heat Transfer Equation for Surrounding Strata

The heat transfer of the surrounding strata is primarily dominated by heat conduction, which can be expressed as(19)ρcs∂Ts∂t+∇⋅λs∇Ts=0.


## 3.6. Cross-Coupling Relations and the Simulation Environment

The established model has a coupling relationship in three aspects. First, the gob energy equation and the gas composition equation were coupled mutually inside the gob and were related to the velocity field solved by the gob flow equation. Secondly, the gas heat transfer equation and gas concentration equation of the butt entries were related to the internal velocity field inside them. Finally, the values of temperature, airflow velocity, pressure, and gas concentration at the interface between different computational domains were equal and continuous. The coupling relationship is shown in Figure2.Figure 2
Multiphysics coupling relationship of the model.The established model consists of several nonlinear partial differential equations in space and time and cannot be analytically solved. Therefore, a multiphysics modeling environment, the COMSOLMultiphysics®, is implemented. COMSOLMultiphysics employs the Finite Element Method and is capable of handling PDEs and DAEs. A segregated solution approach was chosen for the simulation, which splits the solution process into substeps.Besides, with the advance of the working face, the area of the gob will be continuously increased, and the physical size of the calculation area will be constantly varied (as shown in Figure3). Therefore, the calculation grid will be also adjusted accordingly. The deformation geometry model established by COMSOLMultiphysics realized the free deformation of the mesh by setting the moving speed of the calculation boundary, thereby simulating the dynamic growth of the gob.Figure 3
Change of computing domain during the dynamic advancement of the gob.


## 4. Model Verification

To verify the reliability of the multifield coupling model established in this paper, test data acquired from the E12604 working face of Cuijiazhai Coal Mine of Kailuan Group [47] were compared with the calculation results obtained in this paper. The advancing distance of the E12604 working face is 570 m, and the average distance in the dip direction is 135 m. The coal mining method is longwall mining on strike, the coal seam thickness is 3 m–3.2 m, and the mining height is 2.8 m–3 m. The onsite data monitoring point was located 2.5 m above the floor of tailentry pillar close to the gob and 120 m away from the open-off cut of the E12604 working face [48]. According to the actual situation, the boundary and initial conditions of the mathematical model are shown in Tables 1 and 2, and the numerical calculation parameters are shown in Table 3.Table 1
Boundary conditions.
BoundaryAB/CD/EF/FG/GH/IJ/AJDEBCBI/HI/EH/CM/MN/DNBoundary conditionsn⋅λ∇T=0v=vin, cO2=c0, cCO=0, T=Tinp=0outlet, n⋅λ∇T=0, n⋅Di∇ci=0No slip wallTable 2
Initial conditions.
GobSurrounding strataButt entryInitial conditionsT=TW, CO2=C0, CCO=0Ts=TWTb=Tin, Cb,O2=C0, Cb,CO=0Table 3
Numerical simulation parameters.
ParameterValueInlet wind velocityvin1.7 m·s−1Reaction ordern1Preexponential factorA0.424 s−1Activation energyE21.9 kJ·mol−1Initial temperature of the surrounding rockTW17°CAir densityρf1.1 kg·m−3Heat conduction coefficient of airλf0.026 W·m−1·K−1Specific heat of aircp1012 J·kg−1·K−1Air inlet temperatureTin17°COxygen concentration in the inlet of airflowC09.375 mol·m−3Heat of oxidation reactionΔH400 kJ·mol−1Oxygen diffusion coefficientD2 × 10−5 m2·s−1Density of coal gangueρs1300 kg·m−3Heat conduction coefficient of coal gangueλs0.2 W·m−1·K−1Specific heat of coal ganguec1003 J·kg−1·K−1Dynamic viscosityμ1.84 × 105 Pa·sComparison results between the numerical simulation and the measured data are shown in Figure4. Figures 4(a) and 4(b) show the changes in gas temperature and oxygen concentration at the monitoring points in front of the working face, respectively. It can be seen from Figure 4(a) that the simulation and measured results of the temperature of the gob at the early stage of mining were very close. However, the simulation temperature at the measured point was gradually higher than the measured temperature with the increase of advance distance. The main reason for the above phenomenon may be that there will be gas inflow from the adjacent gob in the actual mine. The simulation oxygen concentration in Figure 4(b) was consistent with the measured value. The initial oxygen concentration at the working face was high. However, with the increase of the advance distance, as the coal oxidation gradually consumed the oxygen in the gob, the oxygen concentration decreased continuously. It can be generally considered that the dynamic variation of the coal oxidation process in the gob with the advance of the working face could be well described by the model established in this paper.Figure 4
Comparison between the numerical results and the field data.
(a)(b)The simulation result of the variation of porosity in the gob during the advance of the working face is shown in Figure5. It can be seen from this figure that the distribution of porosity in the gob is not uniform. Overall, the porosity of the area near the working face and the two butt entries is large due to the supporting effect of the hydraulic support and the coal pillars, while the porosity of the area far away from the working face and the double roadway is small, which is mainly due to the good compaction under the action of greater ground stress. Therefore, the distribution of porosity in the gob is characterized by typical O-ring. With the continuous advance of the working face, the scope of the gob area is expanding, and the length of the gob on the strike is increasing. Although the size of the “O” ring increased continuously, the distribution law of porosity in the gob remained unchanged. The minimum porosity in the gob is 0.15 at the initial stage and 0.13 at the later stage, while the porosity near the support of the working face is always maintained at about 0.33, indicating that the gob is gradually compacted with the advance of the working face.Figure 5
Dynamic variation of porosity in the gob during the advance of the working face. (a) 10 days after the advance of the working face. (b) 20 days after the advance of the working face. (c) 40 days after the advance of the working face. (d) 60 days after the advance of the working face.
(a)(b)(c)(d)The simulation result of the variation of the temperature field during the advance of the working face is shown in Figure6. As can be seen from Figures 6(a)–6(d), there is a temperature rise zone within a certain range of the gob during mining. This area first appears near the headentry near the working face. With the advance of the working face, the coverage of the temperature rise zone in the gob gradually expands to the air return side and moves upward and forward. And as the distance in the advancing direction increases, its relative distance from the working face remains unchanged. When the working face advances for a certain distance, the range of the temperature rise zone in the gob moves forward instead of increasing with the advance of the working face, showing a shape of “wing section” on the chart. It can be seen from Figure 6(e) that the high-temperature zone in the gob can still be seen after the mining is stopped and is approaching the terminal line. The terminal line is about 100 m away from the center of the high-temperature zone with the highest temperature of 64.2°C.Figure 6
Dynamic change of temperature field during the advance of the working face. (a) 10 days after the advance of the working face. (b) 20 days after the advance of the working face. (c) 40 days after the advance of the working face. (d) 60 days after the advance of the working face. (e) 10 days after the end of the mining.
(a)(b)(c)(d)(e)The simulation result of the variation of the oxygen concentration field in the gob during the advance of the working face is shown in Figure7. From Figures 7(a)–7(d), it can be seen that the oxygen concentration in the gob is asymmetrically distributed. The high oxygen concentration distribution area in the gob on the air inlet side is deeper than that on the return side. The isoline distribution of oxygen concentration and the relative distance between the working face is fixed. It can be seen from Figure 7(e) that the shape of the oxygen distribution area in the gob at 10 days after the stopping of mining has little change compared with the normal mining period, which is still asymmetrically distributed and the width at the air inlet side is larger than that of the air return side. However, a decline in the total width is observed. The reason is that the path for transporting oxygen in the gob is cut off after the closure of the working face. As a result, the concentration of oxygen in the gob decreases gradually with the oxidation process.Figure 7
Dynamic variation of the oxygen concentration field. (a) 10 days after the advance of the working face. (b) 20 days after the advance of the working face. (c) 40 days after the advance of the working face. (d) 60 days after the advance of the working face. (e) 10 days after the end of the mining.
(a)(b)(c)(d)(e)The simulation result of the variation of the CO concentration field in the gob during the advance of the working face is shown in Figure8. It can be seen from Figures 8(a)–8(e) that the shape of the CO concentration field is similar to that of the oxygen concentration field. But the difference is that the area with a relatively high concentration of CO has a low concentration of oxygen. The CO concentration in the gob gradually increases after the mining starts, extending from the open-off cut to the air return corner of the working face. High concentration of CO is distributed throughout the gob in-depth with the advance of the working face. It is mainly caused by continuous coal oxidation. According to the (16), CO is produced constantly while oxygen is consumed. Due to the weak airflow in the deep gob, CO, which cannot be diluted, accumulates continuously. The temperature of coal-rock next to the working face is relatively low with weak coal-oxygen reaction. As a result, the CO concentration is low under the dilution effect of airflow.Figure 8
Dynamic variation of the CO concentration field. (a) 10 days after the advance of the working face. (b) 20 days after the advance of the working face. (c) 40 days after the advance of the working face. (d) 60 days after the advance of the working face. (e) 10 days after the end of the mining.
(a)(b)(c)(d)(e)


## 5. Sensitivity Analysis

### 5.1. Effect of Ventilation Temperature

The simulation result of the maximum temperature in the gob under varied ventilation temperatures with the temperature of the surrounding strata of 20°C is shown in Figure9. As can be seen from the figure, the maximum temperature in the gob increases rapidly and then decreases slowly after reaching the maximum temperature over time under varied ventilation temperatures. The reason is that the advance distance of the working face is short at the beginning, and the compaction of the carving rock in the gob is low with small air leakage resistance and large air leakage. Under good oxygen supply conditions, residual coal in an oxygen-rich environment leads to an intense coal-oxygen reaction with a large amount of heat release. Meanwhile, few heats are required for heating the environment of the gob due to a small range of the gob, causing a quick temperature rise in the gob. With the continuous advance of the working face, the strike length of the gob increases with the gob range. However, with more and more low-temperature coal and rock entering the gob, the total amount of oxidation heat of the residual coal no longer increases, and the temperature of the gob decreases slightly. Moreover, although the higher the ventilation temperature, the greater the maximum temperature of the gob, the increase of the maximum temperature is smaller than that of the ventilation temperature. Besides, the maximum temperatures of the gob with different ventilation temperatures are similar at the beginning of the advance of the working face. For example, when the ventilation temperature is 10°C, the maximum temperature on the 30th day is 55.2°C; and when the ventilation temperature is 30°C, the maximum temperature on the 30th day is 59.4°C. The ventilation temperature increases by 10°C, while the increase in the maximum temperature is less than 5°C. The reason is that, before entering the gob, the air conducts heat exchange with the walls of the roadway so that the temperature is reduced.Figure 9
The maximum temperature in the gob under varied ventilation temperatures.


### 5.2. Effect of Temperature of the Surrounding Strata

The simulation result of the maximum temperature in the gob under varied temperatures of the surrounding strata is shown in Figure10. As can be seen from the figure, the maximum ambient temperature in the gob is mainly affected by the temperature of the surrounding strata at the initial mining period, and the temperature in the gob is always higher than that of the surrounding strata temperature. For example, when the surrounding strata temperature is 20°C, the maximum ambient temperature of the gob is about 24°C. With the continuous advance of the working face, the maximum ambient temperature of the gob is gradually consistent under different surrounding strata temperatures. The reason is that the oxidation of residual coal gradually becomes dominant. In this way, the gob is gradually heated by the heat produced by the oxidation of the residual coal, so that the ambient temperature in the gob gradually tends to be consistent. Furthermore, the tendency of the maximum temperature in the gob increases rapidly and then declines slowly, which is similar to that in Figure 9.Figure 10
The maximum temperature in the gob under varied temperatures of the surrounding strata.


### 5.3. Effect of Coal-Rock Particle Size

The simulation result of the maximum temperature in the gob under varied particle sizes of the coal-rock is shown in Figure11. As can be seen from the figure, the larger the average particle size, the larger the value of maximum temperature in the gob. For example, the maximum temperature with a particle size of 0.04 m on the 30th day is 56.3°C, and the maximum temperature with a particle size of 0.12 m is 104.2°C. The reason is that the permeability is proportional to the square of the particle size (see (4)). The increase of particle size can improve gob permeability, enhance airflow in the gob, ensure sufficient oxygen supply, enhance the degree of coal-oxygen reaction, and increase the temperature significantly.Figure 11
The maximum temperature in the gob under varied particle sizes.


### 5.4. Effect of the Advance Rate

The simulation result of the maximum temperature in the gob under varied advance rates is shown in Figure12. As can be seen from the figure, the higher the advance rate of the working face, the lower the maximum temperature in the gob. For example, when the advance rate is 2 m/d, the maximum temperature in the gob is close to 194°C, and when the advance rate is 8 m/d, the maximum temperature in the gob is less than 47°C. The reason is that the faster the working face advances, the faster the residual coal in the gob enters the asphyxiation zone. The coal-oxygen reaction is suppressed with less heat release; thus, the temperature in the gob rises slowly.Figure 12
The maximum temperature in the gob under varied advance rates.


## 5.1. Effect of Ventilation Temperature

The simulation result of the maximum temperature in the gob under varied ventilation temperatures with the temperature of the surrounding strata of 20°C is shown in Figure9. As can be seen from the figure, the maximum temperature in the gob increases rapidly and then decreases slowly after reaching the maximum temperature over time under varied ventilation temperatures. The reason is that the advance distance of the working face is short at the beginning, and the compaction of the carving rock in the gob is low with small air leakage resistance and large air leakage. Under good oxygen supply conditions, residual coal in an oxygen-rich environment leads to an intense coal-oxygen reaction with a large amount of heat release. Meanwhile, few heats are required for heating the environment of the gob due to a small range of the gob, causing a quick temperature rise in the gob. With the continuous advance of the working face, the strike length of the gob increases with the gob range. However, with more and more low-temperature coal and rock entering the gob, the total amount of oxidation heat of the residual coal no longer increases, and the temperature of the gob decreases slightly. Moreover, although the higher the ventilation temperature, the greater the maximum temperature of the gob, the increase of the maximum temperature is smaller than that of the ventilation temperature. Besides, the maximum temperatures of the gob with different ventilation temperatures are similar at the beginning of the advance of the working face. For example, when the ventilation temperature is 10°C, the maximum temperature on the 30th day is 55.2°C; and when the ventilation temperature is 30°C, the maximum temperature on the 30th day is 59.4°C. The ventilation temperature increases by 10°C, while the increase in the maximum temperature is less than 5°C. The reason is that, before entering the gob, the air conducts heat exchange with the walls of the roadway so that the temperature is reduced.Figure 9
The maximum temperature in the gob under varied ventilation temperatures.


## 5.2. Effect of Temperature of the Surrounding Strata

The simulation result of the maximum temperature in the gob under varied temperatures of the surrounding strata is shown in Figure10. As can be seen from the figure, the maximum ambient temperature in the gob is mainly affected by the temperature of the surrounding strata at the initial mining period, and the temperature in the gob is always higher than that of the surrounding strata temperature. For example, when the surrounding strata temperature is 20°C, the maximum ambient temperature of the gob is about 24°C. With the continuous advance of the working face, the maximum ambient temperature of the gob is gradually consistent under different surrounding strata temperatures. The reason is that the oxidation of residual coal gradually becomes dominant. In this way, the gob is gradually heated by the heat produced by the oxidation of the residual coal, so that the ambient temperature in the gob gradually tends to be consistent. Furthermore, the tendency of the maximum temperature in the gob increases rapidly and then declines slowly, which is similar to that in Figure 9.Figure 10
The maximum temperature in the gob under varied temperatures of the surrounding strata.


## 5.3. Effect of Coal-Rock Particle Size

The simulation result of the maximum temperature in the gob under varied particle sizes of the coal-rock is shown in Figure11. As can be seen from the figure, the larger the average particle size, the larger the value of maximum temperature in the gob. For example, the maximum temperature with a particle size of 0.04 m on the 30th day is 56.3°C, and the maximum temperature with a particle size of 0.12 m is 104.2°C. The reason is that the permeability is proportional to the square of the particle size (see (4)). The increase of particle size can improve gob permeability, enhance airflow in the gob, ensure sufficient oxygen supply, enhance the degree of coal-oxygen reaction, and increase the temperature significantly.Figure 11
The maximum temperature in the gob under varied particle sizes.


## 5.4. Effect of the Advance Rate

The simulation result of the maximum temperature in the gob under varied advance rates is shown in Figure12. As can be seen from the figure, the higher the advance rate of the working face, the lower the maximum temperature in the gob. For example, when the advance rate is 2 m/d, the maximum temperature in the gob is close to 194°C, and when the advance rate is 8 m/d, the maximum temperature in the gob is less than 47°C. The reason is that the faster the working face advances, the faster the residual coal in the gob enters the asphyxiation zone. The coal-oxygen reaction is suppressed with less heat release; thus, the temperature in the gob rises slowly.Figure 12
The maximum temperature in the gob under varied advance rates.


## 6. Conclusion

A fully coupled transient multiphysical model including gas flow, gas species transport, and heat transfer in the gob and the butt entries, as well as heat transfer in the surrounding strata, is developed to quantify the evolution of coal self-heating in gob during and after mining activities. The rationality of the model was verified by comparing it with the field monitoring data. By analyzing the dynamic variation of the temperature field and oxygen concentration field in the gob during the advance of the working face, some useful conclusions were drawn.(1)
A tailing phenomenon of the high-temperature area in the temperature field is formed along the air inlet side of the gob. The temperature in the high-temperature area formed during the advance of the working face decreases gradually under the accumulation and compaction of the gob and the heat dissipation of the surrounding strata.(2)
The maximum temperature in the gob increases rapidly and then decreases slowly in the advance process of the working face. When mining stops, the temperature in the gob increases continuously, while the high-temperature zone moves towards the working face.(3)
At the same surrounding strata temperature, different ventilation temperatures nearly have no effect on the maximum temperature of the gob at the initial stage of mining. However, the higher the ventilation temperature, the larger the value of maximum temperature in the gob at the later stage of mining.(4)
The smaller the average particle size of gob material or the higher the advance rate of the working face, the lower the maximum temperature in the gob.



---

*Source: 9075657-2020-07-22.xml*