# Study on the Radial and Axial Temperature Fields of Coal Roadway under Seasonal Fluctuation Boundary

Authors: Huaitao Song; Yueping Qin; Jiansong Wu; Gongzhong Wang; Sen Li; Cong Niu; Zeqi Wu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The underground ventilation in coal mine comes from the ground air so that the seasonal fluctuation of atmospheric temperature inevitably has influence upon its climate. The paper obtains a simplified mathematical calculation method of the radial and axial temperature fields to furnish a supplement to air temperature prediction of coal roadway with the emphasis on the influence of seasonal characteristic of inlet airflow. The results reveal that the propagation of radial and axial temperature wave has the similar characteristic of amplitude decay and phase delay. In the radial direction, the main factor is the dimensionless annual cycle time; and in the axial direction, it is influenced by the location-related dimensionless time, the Nusselt number, the Biot number, and the dimensionless annual cycle time. In addition, for the wave amplitude, they present the negative exponential trend as the location increases, the overall shape of which resembles a folded page; for the delayed phase, they present the linear variation trend as the location increases, the overall shape of which resembles a parallelogram. In the paper, the evolution mechanism of temperature field of coal roadway under seasonal fluctuation boundary is profoundly comprehended, which has important guiding significance for the field engineer.


---

## Body

## 1. Introduction

Coal is the biggest energy sources in China. By statistics, it accounted for about 60.4 percent in the primary energy consumption in 2017 [1]. However, as the shallow coal resource is gradually reduced and even exhausted, the mining depth is growing year by year; and China's third resource prediction showed that the coal seam buried depth greater than 1000 m accounted for 59.5 percent of the total forecast. Therefore, it inevitably will suffer from more and more severe heat hazard menace in the not-too-distant future, and is directly linked to endanger the worker’s health, reduce the labor productivity, and increase the frequency of accidents and so on [2–4]. The situation has seriously restricted its safe and efficient exploitation so that government puts more efforts on the solution of the thermal damage problem.The thermal damage governance of high-temperature mine mainly is based on two aspects [5]: the prediction of airflow temperature of coal mine and the implementation of reasonable cooling technical measures. Meanwhile, the accurate prediction for mine climate is a crucial factor that determines the optimal thermal control measures to evaluate the labor conditions and improve the mine environment. However, the variation rules of airflow temperature in underground mine are very complex because of the influence of many factors, principally including the heat dissipation from surrounding rock and the inlet initial temperature from atmosphere [6]. In this perspective, the research on the temperature characteristic under seasonal fluctuation boundary, which includes radial and axial temperature distributions of surrounding rock and airflow, has profound significance to perfect the basic theory of wind temperature prediction in high-temperature mine.There have been many studies focused on the strata heat release and temperature field of underground coal mine roadway. With regard to the radial heat dissipation from surrounding rock, Щсрбаиъ. АП firstly [7] proposed a dimensionless parameter to quantify the heat dissipating capacity. Then an analytical solution for the parameter was established by Cen et al. [8] and he pointed that it was composed of Bessel function and Kelvin function. Sun [9] proposed a simplified analytic expression to calculate the parameter. Other studies also identified that it was an important physical quantity which described the heat dissipation of roadway surrounding rock [10, 11]. After that, Qin et al. [12–14] presented many relation curves between the parameter and Fourier modulus and analyzed the strata heat release using numerical simulation methods. Zhang et al. [15] established a theoretical model of thermal diffusion using finite different method, obtaining the temperature distribution of surrounding rock under the coupling effect of airflow and seepage. At the same time, in the application of business software, Bai et al. [16] analyzed the heat transfer between surrounding rock and airflow based on Fluent software in the case of dry roadway. Danko et al. [17] investigated the transfer phenomena of heat and mass on the wall of roadway based on CLIMSIM software. In the experimental aspect, Zhang et al. [18–20] deduced that the similarity criterion of heat conduction of roadway surrounding rock was the Fo (Fourier number), the Bi (Biot number), and the R (Radius of roadway) in the heat transport process and the Nu (Nusselt number), the Re (Reynold number), and the Pr (Prandtl number) in the heat convection process between surrounding rock and airflow and then certified the feasibility of modeling experiment method. Moreover, the influence of specific factors to the characteristic of temperature field also was analyzed, including the virgin rock temperature [21], the moisture of roadway wall [22–25], the section shape of roadway [26], heterogeneous rocks and anisotropy [27], and ventilation time and ventilation volume [28]. These above factors reflected their influences to heat release of surrounding rock from different sides, which involved the disturbance range of airflow temperature, heat dissipating capacity of surrounding rock, and the distribution rule of temperature field.In addition, for the axial temperature field, Roy et al. [29] developed a computer program to predict the time-dependent climatic change in an underground airway. After that, Liu [30] presented a method of approximate calculation of airflow temperature in mine based on the heat radiation from country rock, the depth of mine and the humidity of the airflow. Gao et al. [31] analyzed the influence of water evaporation place distribution on temperature and humidity of airflow. Xie [32] discussed some influence factors through the accurate observation and analysis of the thermodynamic parameters. Liu [33] deduced the basic law of temperature change in an intake shaft and pointed that it changed with the temperature of airflow at the inlet and the shaft depth. Roy [34] presented a simplified computational approach for prediction of transient climatic conditions in underground mine airways.It was well known that there were two aspects of the periodic change of atmospheric temperature: diurnal variation in a day and annual variation in a year. The ventilation of coal mine came from the ground air so that the regular fluctuation of atmospheric temperature predicatively had influence on its climate. The variation time of diurnal periodicity was relatively short compared with the annual periodicity, and its temperature difference was also small, so it had little influence on the airflow temperature in mine; by contrast, the field measurement indicates that the mining face had obvious annual periodic fluctuation characteristics [35]. However, as mentioned above, most of research only took into account the properties of surrounding rock media and took the temperature of airflow as a constant but ignored its seasonal fluctuation factor. Yi et al. [36] proposed that cooling measures of high-temperature mine should take different implementation on the basis of the heat hazard severity of the season, but without revealing when to implement the optimal cooling methods.In this paper, a technique of variable separation method is introduced to decompose the air inlet temperature of coal mine into two portions: fluctuation boundary and steady boundary. Here, as to the analysis above, we only discuss the influence of seasonal fluctuation of inlet airflow temperature to the radial and axial temperature fields of coal roadway. Through comparing with the analytic solution of a half-infinite plane on even ground, we previously [37] analyzed the mathematical solution of radial temperature field of surrounding rock in detail, and further take a supplement to it, the correctness of which is verified by the experiment in the paper. On this basis, the mathematical solution of axial temperature distribution can be obtained from the law of conservation of energy. Ultimately, all the factors that affect the wave amplitude and phase angle of inlet airflow temperature are synthetically evaluated from the whole space-time perspective when to implement the suitable cooling measures.


## 2. Mathematical Model

### 2.1. Fundamental Assumptions

In the present work, the following basic assumptions have been made during the implementation of heat dissipation model of surrounding rock.(i) As shown in Figure1, the cross-section of roadway is simplified as the shape of circular and the heat transmission within the surrounding rock is only along radial direction. Γ1 and Γ2 represent the inner boundary and the outer boundary, respectively.Figure 1
Influence domain of air temperature.(ii) The surrounding rockmass is homogeneous and isotropic, the original temperature field of which is a constant.(iii) The effect of moisture content and thermal radiation of surrounding rock wall is ignored in this study.(iv) The airflow is laminar and incompressible, and its pressure does not vary as temperature; the air inlet temperature of coal mine is considered as a trigonometric function, which is shown in Figure2.Figure 2
Schematic diagram of seasonal airflow and surrounding rock.


### 2.2. Governing Equation

The governing equation of transient heat conduction of roadway surrounding rock can be established based on the law of conservation of energy and the Fourier's law in the cylindrical coordinate.(1)λrρrCr∂Tr∂t=∂2Tr∂r2+1r∂Tr∂rThe boundary conditions are(2)q=-λr∂Tr∂nΓ1=hTrΓ1-TfTΓ2=TguThe initial conditions are(3)Tt=0=TguTf=T¯f+Afcos⁡2πtτ+φwhereTr is the temperature field of surrounding rockmass, °C; λr is the thermal conductivity of rock, W/(m•°C); ρr is the density of rock, Kg/m3; Cr is the specific heat capacity of rock, J/(Kg∙°C); Tgu is the initial temperature of surrounding rockmass, °C; r is the radial depth, m; h is the convective heat transfer coefficient, W/(m2•°C); q is the heat flux, W/m2; Tf is the seasonal airflow temperature, °C; T-f is the average annual temperature, °C; t is the time of heat exchange between airflow and surrounding rockmass, s; τ is the time of one year (365×24×3600 s). Af is the amplitude of airflow temperature; φ is the phase.


### 2.3. Dimensionless Analysis

For simplicity, the following dimensionless numbers can be introduced to make it more widely applied:(4)θ=Tr-T¯fTgu-T¯f;R=rr0;Bi=hr0λr;Fo=artr02;FoT=arτr02;θAf=AfTgu-T¯fwhereθ is the dimensionless excess temperature; R is the dimensionless radius; r0 is the radius of roadway, m; ar is the thermal diffusivity of surrounding rock, m2/s; Bi is the Biot number; Fo is the Fourier number; FoT is the dimensionless annual cycle time; θAf is the dimensionless fluctuation amplitude.Therefore, (1) can be simplified as(5)∂θ∂Fo=∂2θ∂R2+1R∂θ∂RThe boundary and initial condition also can be similarly obtained based on the above analysis, which can be expressed as(6)-∂θ∂RR=1=BiθR=1-θfθR=∞=1θFo=0=1θf=θAfcos⁡2πFoFoT+φ


### 2.4. Discretized the Governing Equation

Obviously, it is very important that the property of boundary condition is properly analyzed to obtain the accurate result. It is well known that there are three common boundary conditions for heat conduction problem, the first-, second-, and third-type boundary conditions (i.e., Dirichlet, Neumann and Robin, respectively). Compared with the last one, the first two is easily intelligible, which can be taken as its particular case. Unfortunately, the seasonal fluctuation boundary belongs to a more special case. The heat dissipation problem under this condition becomes so complex that it is difficult to implement the numerical simulation straightway. On this ground, we introduce a technique of variable separation method to decompose the complicated fluctuant boundary into two simply parts: one for the average annual temperature of airflowTf- and the other for the periodic temperature with characteristic of cosine wave AfTf′. In addition, this cosine wave AfTf′ can also be divided into two aspects: one for the wave amplitude Af and the other for the cosine function Tf′. Note that the physical quantity Tf′ is a nondimension parameter, which can be defined as(7)Tf′=cos⁡2πtτ+φSimilarly, (1) also consists of two parts: one Ta- represents the influence of the average annual temperature Tf- to the temperature filed of surrounding rock; the other AfTb′ represents the influence of the periodic temperature AfTf′. By that analogy, (5) can also be divided into two parts: θa and θAfθb. (8)Tf=T-f+AfTf′Tr=T-a+AfTb′θ=θa+θAfθbwhere(9)θa=Ta¯-T¯fTgu-T¯fθb=Tb′Therefore, the corresponding definite condition ofθa and θb can be defined as (10)∂θa∂Fo=∂2θa∂R2+1R∂θa∂R∂θa∂RR=1=BiθaR=1θaFo=0=1θaR=∞=1(11)∂θb∂Fo=∂2θb∂R2+1R∂θb∂R∂θb∂RR=1=BiθbR=1-cos⁡2πFoFoT+φθbFo=0=0θbR=∞=0So the problem for (1) is converted to solve the formulation of (10) and (11), respectively. To the best knowledge of the authors, most of the research have studied the impact of the average annual temperature Tf- but neglecting the periodic fluctuation factors Tf′. For this reason, we only discuss the variation regularity of temperature field of roadway surrounding rock under the influence of seasonal fluctuation factors in this article.


## 2.1. Fundamental Assumptions

In the present work, the following basic assumptions have been made during the implementation of heat dissipation model of surrounding rock.(i) As shown in Figure1, the cross-section of roadway is simplified as the shape of circular and the heat transmission within the surrounding rock is only along radial direction. Γ1 and Γ2 represent the inner boundary and the outer boundary, respectively.Figure 1
Influence domain of air temperature.(ii) The surrounding rockmass is homogeneous and isotropic, the original temperature field of which is a constant.(iii) The effect of moisture content and thermal radiation of surrounding rock wall is ignored in this study.(iv) The airflow is laminar and incompressible, and its pressure does not vary as temperature; the air inlet temperature of coal mine is considered as a trigonometric function, which is shown in Figure2.Figure 2
Schematic diagram of seasonal airflow and surrounding rock.


## 2.2. Governing Equation

The governing equation of transient heat conduction of roadway surrounding rock can be established based on the law of conservation of energy and the Fourier's law in the cylindrical coordinate.(1)λrρrCr∂Tr∂t=∂2Tr∂r2+1r∂Tr∂rThe boundary conditions are(2)q=-λr∂Tr∂nΓ1=hTrΓ1-TfTΓ2=TguThe initial conditions are(3)Tt=0=TguTf=T¯f+Afcos⁡2πtτ+φwhereTr is the temperature field of surrounding rockmass, °C; λr is the thermal conductivity of rock, W/(m•°C); ρr is the density of rock, Kg/m3; Cr is the specific heat capacity of rock, J/(Kg∙°C); Tgu is the initial temperature of surrounding rockmass, °C; r is the radial depth, m; h is the convective heat transfer coefficient, W/(m2•°C); q is the heat flux, W/m2; Tf is the seasonal airflow temperature, °C; T-f is the average annual temperature, °C; t is the time of heat exchange between airflow and surrounding rockmass, s; τ is the time of one year (365×24×3600 s). Af is the amplitude of airflow temperature; φ is the phase.


## 2.3. Dimensionless Analysis

For simplicity, the following dimensionless numbers can be introduced to make it more widely applied:(4)θ=Tr-T¯fTgu-T¯f;R=rr0;Bi=hr0λr;Fo=artr02;FoT=arτr02;θAf=AfTgu-T¯fwhereθ is the dimensionless excess temperature; R is the dimensionless radius; r0 is the radius of roadway, m; ar is the thermal diffusivity of surrounding rock, m2/s; Bi is the Biot number; Fo is the Fourier number; FoT is the dimensionless annual cycle time; θAf is the dimensionless fluctuation amplitude.Therefore, (1) can be simplified as(5)∂θ∂Fo=∂2θ∂R2+1R∂θ∂RThe boundary and initial condition also can be similarly obtained based on the above analysis, which can be expressed as(6)-∂θ∂RR=1=BiθR=1-θfθR=∞=1θFo=0=1θf=θAfcos⁡2πFoFoT+φ


## 2.4. Discretized the Governing Equation

Obviously, it is very important that the property of boundary condition is properly analyzed to obtain the accurate result. It is well known that there are three common boundary conditions for heat conduction problem, the first-, second-, and third-type boundary conditions (i.e., Dirichlet, Neumann and Robin, respectively). Compared with the last one, the first two is easily intelligible, which can be taken as its particular case. Unfortunately, the seasonal fluctuation boundary belongs to a more special case. The heat dissipation problem under this condition becomes so complex that it is difficult to implement the numerical simulation straightway. On this ground, we introduce a technique of variable separation method to decompose the complicated fluctuant boundary into two simply parts: one for the average annual temperature of airflowTf- and the other for the periodic temperature with characteristic of cosine wave AfTf′. In addition, this cosine wave AfTf′ can also be divided into two aspects: one for the wave amplitude Af and the other for the cosine function Tf′. Note that the physical quantity Tf′ is a nondimension parameter, which can be defined as(7)Tf′=cos⁡2πtτ+φSimilarly, (1) also consists of two parts: one Ta- represents the influence of the average annual temperature Tf- to the temperature filed of surrounding rock; the other AfTb′ represents the influence of the periodic temperature AfTf′. By that analogy, (5) can also be divided into two parts: θa and θAfθb. (8)Tf=T-f+AfTf′Tr=T-a+AfTb′θ=θa+θAfθbwhere(9)θa=Ta¯-T¯fTgu-T¯fθb=Tb′Therefore, the corresponding definite condition ofθa and θb can be defined as (10)∂θa∂Fo=∂2θa∂R2+1R∂θa∂R∂θa∂RR=1=BiθaR=1θaFo=0=1θaR=∞=1(11)∂θb∂Fo=∂2θb∂R2+1R∂θb∂R∂θb∂RR=1=BiθbR=1-cos⁡2πFoFoT+φθbFo=0=0θbR=∞=0So the problem for (1) is converted to solve the formulation of (10) and (11), respectively. To the best knowledge of the authors, most of the research have studied the impact of the average annual temperature Tf- but neglecting the periodic fluctuation factors Tf′. For this reason, we only discuss the variation regularity of temperature field of roadway surrounding rock under the influence of seasonal fluctuation factors in this article.


## 3. Numerical Analysis

### 3.1. Physical Model Discretization

The given physical model of temperature filed is simplified as a circular roadway, as shown in Figure1, the problem of which has been transformed to solve the axisymmetric heat conduction. Therefore, the discretization of mesh for the physical model can be shown in Figure 3, and the calculated domain from 0 to N depicts from the wall of roadway to a certain depth of surrounding rock. The bold line near the centre denotes the coal roadway, and each thin line represents every cutoff point of the gridding and is divided into N sections (i.e., N=1,2,3,4…). The rule of the gap distance between two adjacent gridding nodes obeys the geometric progression due to the intense heat exchange near the roadway boundary. At the same time, the imaginary line in Figure 3 is located in the middle of two adjacent gridding nodes, and any two adjacent imaginary lines represent one integral domain, which can be considered as a complete cylinder ring. But it can be seen that the integral domain model of the boundary for the node 0 and the node N, comparing with the internal, should be defined as half of the cylinder ring since it only contains one imaginary line.Figure 3
Schematic diagram of mesh generation of surrounding rock.


### 3.2. Mathematical Model Discretization

The thermal balance equations of each cylinder ring within the surrounding rock can be established based on the law of energy conservation, which can be expressed as follows:(12)Ri+Ri-1Ri-Ri-1θi-1k+1+Ri+Ri+1Ri+1-Riθi+1k+1-Ri+Ri-1Ri-Ri-1+Ri+Ri+1Ri+1-Ri+Ri+1-Ri-1Ri+1+Ri-1+2Ri2ΔFoθik+1=-Ri+Ri-1Ri-Ri-1θi-1k-Ri+Ri+1Ri+1-Riθi+1k+Ri+Ri-1Ri-Ri-1+Ri+Ri+1Ri+1-Ri-Ri+1-Ri-1Ri+1+Ri-1+2Ri2ΔFoθikIn the same way, the node 0 and the node N can be displayed as well:(13)2Bi+1+R1R1-1θ0k+1-1+R1R1-1θ1k+1+R1-1R1+32ΔFoθ0k+1=-2Bi+1+R1R1-1θ0k+1+R1R1-1θ1k+R1-1R1+32ΔFoθ0k+4Bicos⁡2πFoFoT+φ(14)θN=0where the subscripts i-1, i, and i+1 denote the node number of the corresponding position; the superscripts k and k-1 denote the time nodes; ΔFo denotes the time step at thekth.Therefore, a complete set of heat conduction equation of surrounding rock has been established by three equations, (12), (13), and (14).The numerical calculation program is developed using VB software based on the aforementioned analysis. In order to reduce the error, the independence of grid and time step has been verified. At last, the radial grid based on the equal proportion is divided by 50 nodes, and the proportionality coefficient is 1.1. The time step is divided as the equal interval, the spacing of which is 1000−1 of the dimensionless annual cycle.


## 3.1. Physical Model Discretization

The given physical model of temperature filed is simplified as a circular roadway, as shown in Figure1, the problem of which has been transformed to solve the axisymmetric heat conduction. Therefore, the discretization of mesh for the physical model can be shown in Figure 3, and the calculated domain from 0 to N depicts from the wall of roadway to a certain depth of surrounding rock. The bold line near the centre denotes the coal roadway, and each thin line represents every cutoff point of the gridding and is divided into N sections (i.e., N=1,2,3,4…). The rule of the gap distance between two adjacent gridding nodes obeys the geometric progression due to the intense heat exchange near the roadway boundary. At the same time, the imaginary line in Figure 3 is located in the middle of two adjacent gridding nodes, and any two adjacent imaginary lines represent one integral domain, which can be considered as a complete cylinder ring. But it can be seen that the integral domain model of the boundary for the node 0 and the node N, comparing with the internal, should be defined as half of the cylinder ring since it only contains one imaginary line.Figure 3
Schematic diagram of mesh generation of surrounding rock.


## 3.2. Mathematical Model Discretization

The thermal balance equations of each cylinder ring within the surrounding rock can be established based on the law of energy conservation, which can be expressed as follows:(12)Ri+Ri-1Ri-Ri-1θi-1k+1+Ri+Ri+1Ri+1-Riθi+1k+1-Ri+Ri-1Ri-Ri-1+Ri+Ri+1Ri+1-Ri+Ri+1-Ri-1Ri+1+Ri-1+2Ri2ΔFoθik+1=-Ri+Ri-1Ri-Ri-1θi-1k-Ri+Ri+1Ri+1-Riθi+1k+Ri+Ri-1Ri-Ri-1+Ri+Ri+1Ri+1-Ri-Ri+1-Ri-1Ri+1+Ri-1+2Ri2ΔFoθikIn the same way, the node 0 and the node N can be displayed as well:(13)2Bi+1+R1R1-1θ0k+1-1+R1R1-1θ1k+1+R1-1R1+32ΔFoθ0k+1=-2Bi+1+R1R1-1θ0k+1+R1R1-1θ1k+R1-1R1+32ΔFoθ0k+4Bicos⁡2πFoFoT+φ(14)θN=0where the subscripts i-1, i, and i+1 denote the node number of the corresponding position; the superscripts k and k-1 denote the time nodes; ΔFo denotes the time step at thekth.Therefore, a complete set of heat conduction equation of surrounding rock has been established by three equations, (12), (13), and (14).The numerical calculation program is developed using VB software based on the aforementioned analysis. In order to reduce the error, the independence of grid and time step has been verified. At last, the radial grid based on the equal proportion is divided by 50 nodes, and the proportionality coefficient is 1.1. The time step is divided as the equal interval, the spacing of which is 1000−1 of the dimensionless annual cycle.


## 4. Result

### 4.1. Mathematical Model of Temperature Distribution along the Radial Direction

To study the regularity of periodic temperature fluctuation of roadway surrounding rock mass, we refer to the temperature distribution function of a semi-infinite medium on even ground [38].(15)Tb′X,Fo=αexp⁡-βXcos⁡2πFoFoT-γX-ηwhere(16)α=1+2Bi-1πFoT+2Bi-2πBi-1/2β=πFoTγ=πFoTη=arctg1+BiFoTπ-1From the analysis of (15), the amplitude α and phaseη is closely related to Bi and Fo. Moreover, the attenuation rate of amplitude α follows the law of negative exponent as the position X increases and the retardation rate of phaseηis linearly related to the position X. The correlation coefficients between them are the function of wave period FoT. It can be seen that the temperature field of surrounding rock under periodic fluctuation has the similar characteristics to the semi-infinite medium’s. Therefore, assume the nondimensional temperature field of surrounding rock under periodic boundary has the similar mathematical expressions, which can be defined as (17)θbX,Fo=Aexp⁡-BXcos⁡2πFoFoT-CX-Dwhere X represents the depth from the surface of surrounding rock to the interior and X=R-1.Based on (17), the wall temperature wave Tb0′(X=0) can be obtained as (18)Tb0′Fo=Acos⁡2πFoFoT-DBy analyzing and comparing with the numerical result of (11), we obtain the similar functional relations between A, B, C, D, and FoT, Bi for (17), which can be expressed using (19). The errors of them relative to the theoretical results are all less than 1%. The detailed derivation process can be referred to in the literature [37].(19)A=2.1325Bi-1FoT-0.3015+1-1B=1.6703FoT-0.3942C=1.7129FoT-0.4803D=arctan⁡0.6122FoT0.4353Bi+0.2784FoT0.3639+1-1


### 4.2. Experimental Verification

The field test in coal mine has certain difficulties to verify the correctness of numerical calculation results due to the limitation of objective conditions like expensive cost, long period, hostile site environment, and so on. Therefore, a method of similarity experiment has been implemented. In addition, because the axial temperature field is the result of theoretical derivation from the radial temperature distribution based on the above analysis, here, we will only verify the temperature distribution rule of the radial direction.


#### 4.2.1. Experimental System

As shown in Figures4 and 5, the experimental platform which is designed according to the similarity criterion is composed by three systems, experimental mimic system, air temperature control system, and temperature field monitoring system, respectively. The geometric parameter of mimic system welded with stainless steel plates is 1000 mm ×1000 mm ×500 mm, and the radius of analogous roadway is 50 mm; the outer surface of test model is covered by the thermal insulation material. The design parameter of control system for inlet air temperature involves the control range -20~150°C, deviation ≤ ±1°C, accuracy ± 0.1°C and air velocity 0~6 m/s. The temperature monitoring system includes the temperature sensor and the data acquisition unit Datataker800, which can be implemented by collecting temperature signal in real time. Figure 6 shows the measuring points layout of thermal resistance in the experimental mimic system, including one in the center of roadway to control and monitor the inlet air temperature, one at the surface of roadway to measure the wall temperature, and the rest within the surrounding rock, the overall outline of which likes a spiral in order to reduce the interplay between sensors.Figure 4
Picture of experiment platform.Figure 5
Schematic diagram of experiment platform.Figure 6
Layout of measurement points in the cross-section.In the process of experiment, the experimental installation cannot keep the inlet air temperature change continuity due to the limitation of its own performance. So in order to reduce its work loads and make the inlet air temperature accord with the nature of the trigonometrical functions, the implement of temperature control in each cycle is averagely divided into 24 pieces based on (3), which then is inputted into the temperature control system.In addition, since the experiment period is relatively short and the temperature of indoor is less affected by that of the outside, we assume that the temperature within surrounding rock model is equal to that of the ambient in initial moment, which can be approximated as a constant 25°C. Thus, the average temperature of inlet air is also set as 25°C and its impact can be ignored. In this paper, the plaster is selected as the main material of similar simulation experiment, in which the powder of iron oxide is added to reinforce its heat conductivity, which is shown in Figure7. The wave amplitude is set as 20°C.Figure 7
Simulation material.
(a)
Gypsum (b)
Iron oxideThen, Table1 lists the thermal physical parameters of the test model, and the experimental schemes can be shown in Table 2. Moreover, the convective heat transfer coefficient between air and surrounding rock model can be calculated through (20), which is based on the Dittus-Boelter formula [39].(20)h=0.0115λfR2RVf0.8vfar-0.4Table 1
Thermal physical parameters of the test model.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (m2/s) Experimental Material 952 0.673 1060 6.67E-07 —— —— Air 1.29 0.025 1010 1.92E-05 1.48E-05 0Table 2
Experimental schemes.
Experimental Schemes Velocity Amplitude Cycle Period Convective Heat Transfer Coefficient Biot Number Dimensionless Period (m/s) (°C) (h) (W/(m2•°C)) (Bi) (FoT) A1 3 20 4 55.42 4.12 3.84 A2 3 20 8 55.42 4.12 7.68 A3 3 20 12 55.42 4.12 11.52 A4 3 20 16 55.42 4.12 15.37 B1 2 20 12 40.06 2.98 11.52 B2 3 20 12 55.42 4.12 11.52 B3 4 20 12 69.76 5.18 11.52 B4 5 20 12 83.39 6.20 11.52where Vf is the velocity of airflow, m/s; vf is the Kinematic viscosity of airflow, m2/s.


#### 4.2.2. Mathematical Model Validation

The wave amplitude and delay phase of numerical calculation are compared with the experimental model. It can be observed from Figures8–11 that the computational results of mathematical model match with the experimental results very well, which shows that this calculation method has high computational accuracy in this paper.Figure 8
Curves of the amplitude of temperature wave with depth at different cycle times.Figure 9
Curves of the amplitude of temperature wave with depth at different wind speeds.Figure 10
Curves of the phase of temperature wave with depth at different cycle times.Figure 11
Curves of the phase of temperature wave with depth at different wind speeds.It can be seen from Figures8 and 9 that the amplitudes of temperature wave show the law of negative exponent as the radial depth increases, and from Figures 10 and 11 that the phases show the linear variation.Furthermore, Figures8 and 10 present that the amplitude and phase of the wall temperature wave (R=0) are very little influenced by the cycle; but in the interior (R>0), the shorter the period, the greater the rate of decay and delay. It is well-know that the mechanism of heat transfers that the high-temperature object transmits its internal energy to the low-temperature object. Hence, when the period of temperature fluctuation has been shortened, the temperature difference per unit time has been increased so that the rate of heat transfer also increases, the characteristic of which exists in the wall and interior. So it not only affects the wall temperature wave of surrounding rock but also the internal.In addition, Figures9 and 11 present, for the wall temperature wave(R=0), that the greater the wind speed, the larger the amplitude but the smaller the phase. Meanwhile, within the surrounding rock (R>0), the gradient of the decay or delay under the condition of different wind speeds is exactly the same. Therefore, the wind speed only affects the wall temperature of surrounding rock but not the internal. The reason is that the heat transport mode between the wall of surrounding rock and air is through convection, which has greatly related to the speed of air current from (20); that is, the parameter A has a close relationship with FoT, Bi; but in the interior, it is through the conduction, which has only relationship with the thermal diffusivity of surrounding rock, as well the parameter C is only a function of FoT.


### 4.3. Mathematical Model of Temperature Distribution along the Axial Direction

According to the above-mentioned analysis, the interior of surrounding rockmass is going to form a temperature field that fluctuates periodically after the seasonal airflow ventilating into the underground mine, and in turn the heat release from the surrounding rock is bound to affect periodically the air. Therefore, the airflow temperature distribution along the axial direction is the result of the interaction between the inlet airflow and the heat release from rock, which makes the air temperature prediction complicated.Based on the theoretical analysis and field observation, we can assume the axial temperature distribution as follows:(21)Taf=Pcos⁡2πtτ+Qwhere P is the amplitude; Q is the phase. It is obvious that they both are associated with the concrete position in the roadway and not related to the time when the inlet air temperature is determined. At the same time, the wall temperature of the same location can be obtained from (18), which can be expressed as follows:(22)TW=PAcos⁡2πtτ+Q-DIn the absence of moisture evaporation, all of the heat lost from the surrounding rockmass is gained by the wind. Therefore, as shown in Figure12, a heat-balance equation in the unit length of dL can be established. (23)GCfdTaf=qUdLFigure 12
Schematic diagram of temperature distribution along the axial direction.where G is the mass flow of air, m3/s; Cf is the specific heat of air, J/(Kg•°C); U is the section perimeter of roadway, m; dTaf is the temperature difference in the distance of dL; L is the length from the ventilation inlet to a site used the airflow, m.(24)G=SVfρfwhere S is the section area of roadway, m2; ρf is the density of air, Kg/m3.Substituting (21), (22), and (2) into (23), we can obtain that (25)cos⁡2πtτ+QdP-Psin⁡2πtτ+QdQ=PWAcos⁡2πtτ+Q-D-cos⁡2πtτ+QdLwhere(26)W=hUGCfSimplify the above equation(27)cos⁡2πtτ+QdP-Psin⁡2πtτ+QdQ=PWAcos⁡D-1cos⁡2πtτ+QdL+APWsin⁡Dsin⁡2πtτ+QdLIt should be noticed that (27) is strict equality in both side at any time and any location. At the same time, there is no correlation between the amplitude P and the phase Q. Therefore, the following formula can be obtained:(28)cos⁡2πtτ+QdP=PWAcos⁡D-1cos⁡2πtτ+QdL-Psin⁡2πtτ+QdQ=APWsin⁡Dsin⁡2πtτ+QdLFurther simplification, we can obtain that(29)P=P0exp⁡WLAcos⁡D-1Q=Q0-AWLsin⁡Dwhere P0 and Q0 are the initial amplitude and the phase in the inlet of roadway, respectively. From (3), that is to say,(30)P0=AfQ0=φFrom the discussion above, the calculation formula of amplitude and phase of the air temperature along the axial direction has been deduced. So the conclusion can be obtained that the amplitude shows the law of negative exponent along the axial direction of the roadway, and the phase changes linearly.Substituting (29) and (30) into (21), the axial temperature distribution of roadway can be obtained as follows:(31)Taf=Afexp⁡WLAcos⁡D-1cos⁡2πtτ+φ-AWLsin⁡Dbecause of(32)WL=hULGCf=21r2LVfλfρfCfhrλf=2FofNufwhereNuf represents the Nusselt number; Fof represents the dimensionless cycle time that related location. Therefore, in the axial direction, there also is a dimensionless cycle time Fof′ on the synchronous occurrence of the maximum value of temperature, and it can be calculated as follows based on (31) and (32):(33)Fof′=πNufAsin⁡DAs a result, we can obtain the dimensionless expression of temperature distribution along the roadway.(34)θFf=θfexp⁡2FofNufAcos⁡D-1cos⁡2πFoFoT+φ-2FofNufAsin⁡D


### 4.4. Project Application of Mathematical Model of Temperature Field

#### 4.4.1. Calculation Scheme

This paper takes a specific mine as an example in order to systematically analyze the influencing factors of temperature field. The relevant parameters are shown in Table3 [40].Table 3
Thermal physical parameters of surrounding rock and air.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Wave Amplitude Cycle Time Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (°C) (s) (m2/s) —— Surrounding rock 2500 3.25 1000 1.30E-06 —— —— —— —— Air 1.29 0.025 1010 1.92E-05 16 365×24×3600 1.48E-05 0However, for any specific coal mine, the main factors that can be controlled by human are generally only the ventilation velocity and the section area of roadway, the influence of which on the amplitude and phase is analyzed based on the coal mine safety regulations. It stipulates that the wind speed in the main intake and return airways should not exceed 8 m/s, and the net section of roadway must meet the needs of pedestrian, transportation, ventilation, safety facilities, equipment installation, overhaul, and construction. Therefore, the following calculation schemes can be listed in Table4.Table 4
Calculation schemes.
Calculation Schemes Velocity Radius Convective Heat Transfer Coefficient Biot Number Dimensionless Annual Period Nusselt Number (m/s) (m) (W/(m2•°C)) (Bi) (FoT) (Nuf) A1 0.5 3 4.46 4.12 4.56 535.44 A2 1 3 7.77 7.17 4.56 932.26 A3 3 3 18.71 17.27 4.56 2245.09 A4 5 3 28.15 25.99 4.56 3378.41 A5 7 3 36.85 34.02 4.56 4421.96 B1 3 1 23.31 7.17 41.00 932.26 B2 3 2 20.29 12.49 10.25 1623.16 B3 3 3 18.71 17.27 4.56 2245.09 B4 3 4 17.66 21.74 2.56 2826.08 B5 3 5 16.89 25.99 1.64 3378.41


#### 4.4.2. Result Analysis

It can be seen from Figures13 and 14 that the farther the entrance of roadway along the axial direction, the larger the amplitude of temperature wave, and it shows the law of negative exponent as the axial depth increases. In addition, the cumulative distance from the ventilation inlet to the working face is generally approximately 6000 m, and the variation range of the amplitude is around 2~6°C, which is consistent with that of field observation at the working face [41].Figure 13
Curves of the amplitude of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 14
Curves of the amplitude of temperature wave at different roadway radius under the condition of the wind speed 3m/s.It can be seen from Figures15 and 16 that the farther the entrance of roadway along the axial direction, the more delayed the phase of temperature wave, and it shows the linear variation rule as the axial depth increases. In addition, when the cumulative distance is about 6000m, the variation range of the delayed phase is around 0.5~2, namely, the temperature wave of the working face lag behind that of the ventilation inlet about 30~116 days. At the same time, the field observations at the working face shows that the maximum value of ground temperature in a year appears in the mid-to-late July and that of working face in early September [38], which agrees with the calculation.Figure 15
Curves of the phase of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 16
Curves of the phase of temperature wave at different roadway radius under the condition of the wind speed 3m/s.Figures13 and 15 present that the bigger the wind speed, the smaller the attenuation rate of wave amplitude and the slower the delay rate of the phase. The reason is that the heat exchange time between air and surrounding rock is shortened when the wind speed increased, so the attenuation rate and delay rate at unit time is decreased.Figures14 and 16 present that the bigger the radius of roadway, the smaller the attenuation of wave amplitude and the slower the delay rate of the phase. Due to the fact that the wind speed is a constant, the rate of heat exchange is also assumed to be constant. Moreover, the ventilation quantity in unit time inevitably increases as the roadway radius enlarges, which lead to a decreased attenuation rate and delay rate at unit length.On the incorporation of the radial and axial temperature distribution, the synthesis analysis of temperature field in three-dimensional space can be obtained when the wind speed is 3m/s and the radius is 3m.Figure17 shows that the surrounding rock and air have the attenuation effect on the amplitude of temperature wave, and both present the law of negative exponent as the corresponding location increases, the overall shape of which resembles a folded page. Figure 18 shows that they both have the postponed effect on the phase and present the linear variation trend as the corresponding location increases, the overall shape of which resembles a parallelogram.Figure 17
Variation law of wave amplitude in three-dimensional space.Figure 18
Variation law of delayed phase in three-dimensional space.Moreover, compared with the air, the decay rate or delay rate of temperature wave within surrounding rock is relatively fast, which is affected by the dimensionless annual cycle time (thermal diffusivity of rock, radius of roadway, cycle time, and radial depth). However, in the axial direction, they are influenced by the location-related dimensionless time (thermal diffusivity of air, radius of roadway, and air velocity), the Nusselt number, the Biot number, and the dimensionless annual cycle time. Therefore, whether the decay of amplitude or the delay of phase, the effect of the rock is greater than that of the airflow.Figure19 shows that the ventilation distance has great influence on the amplitude and the phase of air temperature wave. The location-related cycle time along the axial direction can be calculated based on (33) and has close connection with the ventilation distance and velocity. In this instance, it is about 4.3 hours or the distance is 46600m. As mentioned previously, the ventilation rate in mine is relatively small and the maximum is 8m/s. Therefore, it is difficult to achieve a complete cycle fluctuation along the axial direction at one time.Figure 19
Variation law of axial temperature field in a time cycle.From Figure20, the temperature of air current may be greater than that of surrounding rock in the hottest season, and then the air could dissipate constantly heat to the surrounding rock. The other way round, the coldest season could present the opposite characteristic. That how to make good use of these features to govern the mine heat damage needs further discussion.Figure 20
Variation law of axial temperature field with the time.


## 4.1. Mathematical Model of Temperature Distribution along the Radial Direction

To study the regularity of periodic temperature fluctuation of roadway surrounding rock mass, we refer to the temperature distribution function of a semi-infinite medium on even ground [38].(15)Tb′X,Fo=αexp⁡-βXcos⁡2πFoFoT-γX-ηwhere(16)α=1+2Bi-1πFoT+2Bi-2πBi-1/2β=πFoTγ=πFoTη=arctg1+BiFoTπ-1From the analysis of (15), the amplitude α and phaseη is closely related to Bi and Fo. Moreover, the attenuation rate of amplitude α follows the law of negative exponent as the position X increases and the retardation rate of phaseηis linearly related to the position X. The correlation coefficients between them are the function of wave period FoT. It can be seen that the temperature field of surrounding rock under periodic fluctuation has the similar characteristics to the semi-infinite medium’s. Therefore, assume the nondimensional temperature field of surrounding rock under periodic boundary has the similar mathematical expressions, which can be defined as (17)θbX,Fo=Aexp⁡-BXcos⁡2πFoFoT-CX-Dwhere X represents the depth from the surface of surrounding rock to the interior and X=R-1.Based on (17), the wall temperature wave Tb0′(X=0) can be obtained as (18)Tb0′Fo=Acos⁡2πFoFoT-DBy analyzing and comparing with the numerical result of (11), we obtain the similar functional relations between A, B, C, D, and FoT, Bi for (17), which can be expressed using (19). The errors of them relative to the theoretical results are all less than 1%. The detailed derivation process can be referred to in the literature [37].(19)A=2.1325Bi-1FoT-0.3015+1-1B=1.6703FoT-0.3942C=1.7129FoT-0.4803D=arctan⁡0.6122FoT0.4353Bi+0.2784FoT0.3639+1-1


## 4.2. Experimental Verification

The field test in coal mine has certain difficulties to verify the correctness of numerical calculation results due to the limitation of objective conditions like expensive cost, long period, hostile site environment, and so on. Therefore, a method of similarity experiment has been implemented. In addition, because the axial temperature field is the result of theoretical derivation from the radial temperature distribution based on the above analysis, here, we will only verify the temperature distribution rule of the radial direction.


### 4.2.1. Experimental System

As shown in Figures4 and 5, the experimental platform which is designed according to the similarity criterion is composed by three systems, experimental mimic system, air temperature control system, and temperature field monitoring system, respectively. The geometric parameter of mimic system welded with stainless steel plates is 1000 mm ×1000 mm ×500 mm, and the radius of analogous roadway is 50 mm; the outer surface of test model is covered by the thermal insulation material. The design parameter of control system for inlet air temperature involves the control range -20~150°C, deviation ≤ ±1°C, accuracy ± 0.1°C and air velocity 0~6 m/s. The temperature monitoring system includes the temperature sensor and the data acquisition unit Datataker800, which can be implemented by collecting temperature signal in real time. Figure 6 shows the measuring points layout of thermal resistance in the experimental mimic system, including one in the center of roadway to control and monitor the inlet air temperature, one at the surface of roadway to measure the wall temperature, and the rest within the surrounding rock, the overall outline of which likes a spiral in order to reduce the interplay between sensors.Figure 4
Picture of experiment platform.Figure 5
Schematic diagram of experiment platform.Figure 6
Layout of measurement points in the cross-section.In the process of experiment, the experimental installation cannot keep the inlet air temperature change continuity due to the limitation of its own performance. So in order to reduce its work loads and make the inlet air temperature accord with the nature of the trigonometrical functions, the implement of temperature control in each cycle is averagely divided into 24 pieces based on (3), which then is inputted into the temperature control system.In addition, since the experiment period is relatively short and the temperature of indoor is less affected by that of the outside, we assume that the temperature within surrounding rock model is equal to that of the ambient in initial moment, which can be approximated as a constant 25°C. Thus, the average temperature of inlet air is also set as 25°C and its impact can be ignored. In this paper, the plaster is selected as the main material of similar simulation experiment, in which the powder of iron oxide is added to reinforce its heat conductivity, which is shown in Figure7. The wave amplitude is set as 20°C.Figure 7
Simulation material.
(a)
Gypsum (b)
Iron oxideThen, Table1 lists the thermal physical parameters of the test model, and the experimental schemes can be shown in Table 2. Moreover, the convective heat transfer coefficient between air and surrounding rock model can be calculated through (20), which is based on the Dittus-Boelter formula [39].(20)h=0.0115λfR2RVf0.8vfar-0.4Table 1
Thermal physical parameters of the test model.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (m2/s) Experimental Material 952 0.673 1060 6.67E-07 —— —— Air 1.29 0.025 1010 1.92E-05 1.48E-05 0Table 2
Experimental schemes.
Experimental Schemes Velocity Amplitude Cycle Period Convective Heat Transfer Coefficient Biot Number Dimensionless Period (m/s) (°C) (h) (W/(m2•°C)) (Bi) (FoT) A1 3 20 4 55.42 4.12 3.84 A2 3 20 8 55.42 4.12 7.68 A3 3 20 12 55.42 4.12 11.52 A4 3 20 16 55.42 4.12 15.37 B1 2 20 12 40.06 2.98 11.52 B2 3 20 12 55.42 4.12 11.52 B3 4 20 12 69.76 5.18 11.52 B4 5 20 12 83.39 6.20 11.52where Vf is the velocity of airflow, m/s; vf is the Kinematic viscosity of airflow, m2/s.


### 4.2.2. Mathematical Model Validation

The wave amplitude and delay phase of numerical calculation are compared with the experimental model. It can be observed from Figures8–11 that the computational results of mathematical model match with the experimental results very well, which shows that this calculation method has high computational accuracy in this paper.Figure 8
Curves of the amplitude of temperature wave with depth at different cycle times.Figure 9
Curves of the amplitude of temperature wave with depth at different wind speeds.Figure 10
Curves of the phase of temperature wave with depth at different cycle times.Figure 11
Curves of the phase of temperature wave with depth at different wind speeds.It can be seen from Figures8 and 9 that the amplitudes of temperature wave show the law of negative exponent as the radial depth increases, and from Figures 10 and 11 that the phases show the linear variation.Furthermore, Figures8 and 10 present that the amplitude and phase of the wall temperature wave (R=0) are very little influenced by the cycle; but in the interior (R>0), the shorter the period, the greater the rate of decay and delay. It is well-know that the mechanism of heat transfers that the high-temperature object transmits its internal energy to the low-temperature object. Hence, when the period of temperature fluctuation has been shortened, the temperature difference per unit time has been increased so that the rate of heat transfer also increases, the characteristic of which exists in the wall and interior. So it not only affects the wall temperature wave of surrounding rock but also the internal.In addition, Figures9 and 11 present, for the wall temperature wave(R=0), that the greater the wind speed, the larger the amplitude but the smaller the phase. Meanwhile, within the surrounding rock (R>0), the gradient of the decay or delay under the condition of different wind speeds is exactly the same. Therefore, the wind speed only affects the wall temperature of surrounding rock but not the internal. The reason is that the heat transport mode between the wall of surrounding rock and air is through convection, which has greatly related to the speed of air current from (20); that is, the parameter A has a close relationship with FoT, Bi; but in the interior, it is through the conduction, which has only relationship with the thermal diffusivity of surrounding rock, as well the parameter C is only a function of FoT.


## 4.2.1. Experimental System

As shown in Figures4 and 5, the experimental platform which is designed according to the similarity criterion is composed by three systems, experimental mimic system, air temperature control system, and temperature field monitoring system, respectively. The geometric parameter of mimic system welded with stainless steel plates is 1000 mm ×1000 mm ×500 mm, and the radius of analogous roadway is 50 mm; the outer surface of test model is covered by the thermal insulation material. The design parameter of control system for inlet air temperature involves the control range -20~150°C, deviation ≤ ±1°C, accuracy ± 0.1°C and air velocity 0~6 m/s. The temperature monitoring system includes the temperature sensor and the data acquisition unit Datataker800, which can be implemented by collecting temperature signal in real time. Figure 6 shows the measuring points layout of thermal resistance in the experimental mimic system, including one in the center of roadway to control and monitor the inlet air temperature, one at the surface of roadway to measure the wall temperature, and the rest within the surrounding rock, the overall outline of which likes a spiral in order to reduce the interplay between sensors.Figure 4
Picture of experiment platform.Figure 5
Schematic diagram of experiment platform.Figure 6
Layout of measurement points in the cross-section.In the process of experiment, the experimental installation cannot keep the inlet air temperature change continuity due to the limitation of its own performance. So in order to reduce its work loads and make the inlet air temperature accord with the nature of the trigonometrical functions, the implement of temperature control in each cycle is averagely divided into 24 pieces based on (3), which then is inputted into the temperature control system.In addition, since the experiment period is relatively short and the temperature of indoor is less affected by that of the outside, we assume that the temperature within surrounding rock model is equal to that of the ambient in initial moment, which can be approximated as a constant 25°C. Thus, the average temperature of inlet air is also set as 25°C and its impact can be ignored. In this paper, the plaster is selected as the main material of similar simulation experiment, in which the powder of iron oxide is added to reinforce its heat conductivity, which is shown in Figure7. The wave amplitude is set as 20°C.Figure 7
Simulation material.
(a)
Gypsum (b)
Iron oxideThen, Table1 lists the thermal physical parameters of the test model, and the experimental schemes can be shown in Table 2. Moreover, the convective heat transfer coefficient between air and surrounding rock model can be calculated through (20), which is based on the Dittus-Boelter formula [39].(20)h=0.0115λfR2RVf0.8vfar-0.4Table 1
Thermal physical parameters of the test model.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (m2/s) Experimental Material 952 0.673 1060 6.67E-07 —— —— Air 1.29 0.025 1010 1.92E-05 1.48E-05 0Table 2
Experimental schemes.
Experimental Schemes Velocity Amplitude Cycle Period Convective Heat Transfer Coefficient Biot Number Dimensionless Period (m/s) (°C) (h) (W/(m2•°C)) (Bi) (FoT) A1 3 20 4 55.42 4.12 3.84 A2 3 20 8 55.42 4.12 7.68 A3 3 20 12 55.42 4.12 11.52 A4 3 20 16 55.42 4.12 15.37 B1 2 20 12 40.06 2.98 11.52 B2 3 20 12 55.42 4.12 11.52 B3 4 20 12 69.76 5.18 11.52 B4 5 20 12 83.39 6.20 11.52where Vf is the velocity of airflow, m/s; vf is the Kinematic viscosity of airflow, m2/s.


## 4.2.2. Mathematical Model Validation

The wave amplitude and delay phase of numerical calculation are compared with the experimental model. It can be observed from Figures8–11 that the computational results of mathematical model match with the experimental results very well, which shows that this calculation method has high computational accuracy in this paper.Figure 8
Curves of the amplitude of temperature wave with depth at different cycle times.Figure 9
Curves of the amplitude of temperature wave with depth at different wind speeds.Figure 10
Curves of the phase of temperature wave with depth at different cycle times.Figure 11
Curves of the phase of temperature wave with depth at different wind speeds.It can be seen from Figures8 and 9 that the amplitudes of temperature wave show the law of negative exponent as the radial depth increases, and from Figures 10 and 11 that the phases show the linear variation.Furthermore, Figures8 and 10 present that the amplitude and phase of the wall temperature wave (R=0) are very little influenced by the cycle; but in the interior (R>0), the shorter the period, the greater the rate of decay and delay. It is well-know that the mechanism of heat transfers that the high-temperature object transmits its internal energy to the low-temperature object. Hence, when the period of temperature fluctuation has been shortened, the temperature difference per unit time has been increased so that the rate of heat transfer also increases, the characteristic of which exists in the wall and interior. So it not only affects the wall temperature wave of surrounding rock but also the internal.In addition, Figures9 and 11 present, for the wall temperature wave(R=0), that the greater the wind speed, the larger the amplitude but the smaller the phase. Meanwhile, within the surrounding rock (R>0), the gradient of the decay or delay under the condition of different wind speeds is exactly the same. Therefore, the wind speed only affects the wall temperature of surrounding rock but not the internal. The reason is that the heat transport mode between the wall of surrounding rock and air is through convection, which has greatly related to the speed of air current from (20); that is, the parameter A has a close relationship with FoT, Bi; but in the interior, it is through the conduction, which has only relationship with the thermal diffusivity of surrounding rock, as well the parameter C is only a function of FoT.


## 4.3. Mathematical Model of Temperature Distribution along the Axial Direction

According to the above-mentioned analysis, the interior of surrounding rockmass is going to form a temperature field that fluctuates periodically after the seasonal airflow ventilating into the underground mine, and in turn the heat release from the surrounding rock is bound to affect periodically the air. Therefore, the airflow temperature distribution along the axial direction is the result of the interaction between the inlet airflow and the heat release from rock, which makes the air temperature prediction complicated.Based on the theoretical analysis and field observation, we can assume the axial temperature distribution as follows:(21)Taf=Pcos⁡2πtτ+Qwhere P is the amplitude; Q is the phase. It is obvious that they both are associated with the concrete position in the roadway and not related to the time when the inlet air temperature is determined. At the same time, the wall temperature of the same location can be obtained from (18), which can be expressed as follows:(22)TW=PAcos⁡2πtτ+Q-DIn the absence of moisture evaporation, all of the heat lost from the surrounding rockmass is gained by the wind. Therefore, as shown in Figure12, a heat-balance equation in the unit length of dL can be established. (23)GCfdTaf=qUdLFigure 12
Schematic diagram of temperature distribution along the axial direction.where G is the mass flow of air, m3/s; Cf is the specific heat of air, J/(Kg•°C); U is the section perimeter of roadway, m; dTaf is the temperature difference in the distance of dL; L is the length from the ventilation inlet to a site used the airflow, m.(24)G=SVfρfwhere S is the section area of roadway, m2; ρf is the density of air, Kg/m3.Substituting (21), (22), and (2) into (23), we can obtain that (25)cos⁡2πtτ+QdP-Psin⁡2πtτ+QdQ=PWAcos⁡2πtτ+Q-D-cos⁡2πtτ+QdLwhere(26)W=hUGCfSimplify the above equation(27)cos⁡2πtτ+QdP-Psin⁡2πtτ+QdQ=PWAcos⁡D-1cos⁡2πtτ+QdL+APWsin⁡Dsin⁡2πtτ+QdLIt should be noticed that (27) is strict equality in both side at any time and any location. At the same time, there is no correlation between the amplitude P and the phase Q. Therefore, the following formula can be obtained:(28)cos⁡2πtτ+QdP=PWAcos⁡D-1cos⁡2πtτ+QdL-Psin⁡2πtτ+QdQ=APWsin⁡Dsin⁡2πtτ+QdLFurther simplification, we can obtain that(29)P=P0exp⁡WLAcos⁡D-1Q=Q0-AWLsin⁡Dwhere P0 and Q0 are the initial amplitude and the phase in the inlet of roadway, respectively. From (3), that is to say,(30)P0=AfQ0=φFrom the discussion above, the calculation formula of amplitude and phase of the air temperature along the axial direction has been deduced. So the conclusion can be obtained that the amplitude shows the law of negative exponent along the axial direction of the roadway, and the phase changes linearly.Substituting (29) and (30) into (21), the axial temperature distribution of roadway can be obtained as follows:(31)Taf=Afexp⁡WLAcos⁡D-1cos⁡2πtτ+φ-AWLsin⁡Dbecause of(32)WL=hULGCf=21r2LVfλfρfCfhrλf=2FofNufwhereNuf represents the Nusselt number; Fof represents the dimensionless cycle time that related location. Therefore, in the axial direction, there also is a dimensionless cycle time Fof′ on the synchronous occurrence of the maximum value of temperature, and it can be calculated as follows based on (31) and (32):(33)Fof′=πNufAsin⁡DAs a result, we can obtain the dimensionless expression of temperature distribution along the roadway.(34)θFf=θfexp⁡2FofNufAcos⁡D-1cos⁡2πFoFoT+φ-2FofNufAsin⁡D


## 4.4. Project Application of Mathematical Model of Temperature Field

### 4.4.1. Calculation Scheme

This paper takes a specific mine as an example in order to systematically analyze the influencing factors of temperature field. The relevant parameters are shown in Table3 [40].Table 3
Thermal physical parameters of surrounding rock and air.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Wave Amplitude Cycle Time Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (°C) (s) (m2/s) —— Surrounding rock 2500 3.25 1000 1.30E-06 —— —— —— —— Air 1.29 0.025 1010 1.92E-05 16 365×24×3600 1.48E-05 0However, for any specific coal mine, the main factors that can be controlled by human are generally only the ventilation velocity and the section area of roadway, the influence of which on the amplitude and phase is analyzed based on the coal mine safety regulations. It stipulates that the wind speed in the main intake and return airways should not exceed 8 m/s, and the net section of roadway must meet the needs of pedestrian, transportation, ventilation, safety facilities, equipment installation, overhaul, and construction. Therefore, the following calculation schemes can be listed in Table4.Table 4
Calculation schemes.
Calculation Schemes Velocity Radius Convective Heat Transfer Coefficient Biot Number Dimensionless Annual Period Nusselt Number (m/s) (m) (W/(m2•°C)) (Bi) (FoT) (Nuf) A1 0.5 3 4.46 4.12 4.56 535.44 A2 1 3 7.77 7.17 4.56 932.26 A3 3 3 18.71 17.27 4.56 2245.09 A4 5 3 28.15 25.99 4.56 3378.41 A5 7 3 36.85 34.02 4.56 4421.96 B1 3 1 23.31 7.17 41.00 932.26 B2 3 2 20.29 12.49 10.25 1623.16 B3 3 3 18.71 17.27 4.56 2245.09 B4 3 4 17.66 21.74 2.56 2826.08 B5 3 5 16.89 25.99 1.64 3378.41


### 4.4.2. Result Analysis

It can be seen from Figures13 and 14 that the farther the entrance of roadway along the axial direction, the larger the amplitude of temperature wave, and it shows the law of negative exponent as the axial depth increases. In addition, the cumulative distance from the ventilation inlet to the working face is generally approximately 6000 m, and the variation range of the amplitude is around 2~6°C, which is consistent with that of field observation at the working face [41].Figure 13
Curves of the amplitude of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 14
Curves of the amplitude of temperature wave at different roadway radius under the condition of the wind speed 3m/s.It can be seen from Figures15 and 16 that the farther the entrance of roadway along the axial direction, the more delayed the phase of temperature wave, and it shows the linear variation rule as the axial depth increases. In addition, when the cumulative distance is about 6000m, the variation range of the delayed phase is around 0.5~2, namely, the temperature wave of the working face lag behind that of the ventilation inlet about 30~116 days. At the same time, the field observations at the working face shows that the maximum value of ground temperature in a year appears in the mid-to-late July and that of working face in early September [38], which agrees with the calculation.Figure 15
Curves of the phase of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 16
Curves of the phase of temperature wave at different roadway radius under the condition of the wind speed 3m/s.Figures13 and 15 present that the bigger the wind speed, the smaller the attenuation rate of wave amplitude and the slower the delay rate of the phase. The reason is that the heat exchange time between air and surrounding rock is shortened when the wind speed increased, so the attenuation rate and delay rate at unit time is decreased.Figures14 and 16 present that the bigger the radius of roadway, the smaller the attenuation of wave amplitude and the slower the delay rate of the phase. Due to the fact that the wind speed is a constant, the rate of heat exchange is also assumed to be constant. Moreover, the ventilation quantity in unit time inevitably increases as the roadway radius enlarges, which lead to a decreased attenuation rate and delay rate at unit length.On the incorporation of the radial and axial temperature distribution, the synthesis analysis of temperature field in three-dimensional space can be obtained when the wind speed is 3m/s and the radius is 3m.Figure17 shows that the surrounding rock and air have the attenuation effect on the amplitude of temperature wave, and both present the law of negative exponent as the corresponding location increases, the overall shape of which resembles a folded page. Figure 18 shows that they both have the postponed effect on the phase and present the linear variation trend as the corresponding location increases, the overall shape of which resembles a parallelogram.Figure 17
Variation law of wave amplitude in three-dimensional space.Figure 18
Variation law of delayed phase in three-dimensional space.Moreover, compared with the air, the decay rate or delay rate of temperature wave within surrounding rock is relatively fast, which is affected by the dimensionless annual cycle time (thermal diffusivity of rock, radius of roadway, cycle time, and radial depth). However, in the axial direction, they are influenced by the location-related dimensionless time (thermal diffusivity of air, radius of roadway, and air velocity), the Nusselt number, the Biot number, and the dimensionless annual cycle time. Therefore, whether the decay of amplitude or the delay of phase, the effect of the rock is greater than that of the airflow.Figure19 shows that the ventilation distance has great influence on the amplitude and the phase of air temperature wave. The location-related cycle time along the axial direction can be calculated based on (33) and has close connection with the ventilation distance and velocity. In this instance, it is about 4.3 hours or the distance is 46600m. As mentioned previously, the ventilation rate in mine is relatively small and the maximum is 8m/s. Therefore, it is difficult to achieve a complete cycle fluctuation along the axial direction at one time.Figure 19
Variation law of axial temperature field in a time cycle.From Figure20, the temperature of air current may be greater than that of surrounding rock in the hottest season, and then the air could dissipate constantly heat to the surrounding rock. The other way round, the coldest season could present the opposite characteristic. That how to make good use of these features to govern the mine heat damage needs further discussion.Figure 20
Variation law of axial temperature field with the time.


## 4.4.1. Calculation Scheme

This paper takes a specific mine as an example in order to systematically analyze the influencing factors of temperature field. The relevant parameters are shown in Table3 [40].Table 3
Thermal physical parameters of surrounding rock and air.
Medium Density Thermal Conductivity Specific Heat Capacity Thermal Diffusivity Wave Amplitude Cycle Time Kinematic Viscosity Initial Phase (Kg/m3) W/(m•°C) J/(Kg•°C) (m2/s) (°C) (s) (m2/s) —— Surrounding rock 2500 3.25 1000 1.30E-06 —— —— —— —— Air 1.29 0.025 1010 1.92E-05 16 365×24×3600 1.48E-05 0However, for any specific coal mine, the main factors that can be controlled by human are generally only the ventilation velocity and the section area of roadway, the influence of which on the amplitude and phase is analyzed based on the coal mine safety regulations. It stipulates that the wind speed in the main intake and return airways should not exceed 8 m/s, and the net section of roadway must meet the needs of pedestrian, transportation, ventilation, safety facilities, equipment installation, overhaul, and construction. Therefore, the following calculation schemes can be listed in Table4.Table 4
Calculation schemes.
Calculation Schemes Velocity Radius Convective Heat Transfer Coefficient Biot Number Dimensionless Annual Period Nusselt Number (m/s) (m) (W/(m2•°C)) (Bi) (FoT) (Nuf) A1 0.5 3 4.46 4.12 4.56 535.44 A2 1 3 7.77 7.17 4.56 932.26 A3 3 3 18.71 17.27 4.56 2245.09 A4 5 3 28.15 25.99 4.56 3378.41 A5 7 3 36.85 34.02 4.56 4421.96 B1 3 1 23.31 7.17 41.00 932.26 B2 3 2 20.29 12.49 10.25 1623.16 B3 3 3 18.71 17.27 4.56 2245.09 B4 3 4 17.66 21.74 2.56 2826.08 B5 3 5 16.89 25.99 1.64 3378.41


## 4.4.2. Result Analysis

It can be seen from Figures13 and 14 that the farther the entrance of roadway along the axial direction, the larger the amplitude of temperature wave, and it shows the law of negative exponent as the axial depth increases. In addition, the cumulative distance from the ventilation inlet to the working face is generally approximately 6000 m, and the variation range of the amplitude is around 2~6°C, which is consistent with that of field observation at the working face [41].Figure 13
Curves of the amplitude of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 14
Curves of the amplitude of temperature wave at different roadway radius under the condition of the wind speed 3m/s.It can be seen from Figures15 and 16 that the farther the entrance of roadway along the axial direction, the more delayed the phase of temperature wave, and it shows the linear variation rule as the axial depth increases. In addition, when the cumulative distance is about 6000m, the variation range of the delayed phase is around 0.5~2, namely, the temperature wave of the working face lag behind that of the ventilation inlet about 30~116 days. At the same time, the field observations at the working face shows that the maximum value of ground temperature in a year appears in the mid-to-late July and that of working face in early September [38], which agrees with the calculation.Figure 15
Curves of the phase of temperature wave at different wind speeds under the condition of the roadway radius 3m.Figure 16
Curves of the phase of temperature wave at different roadway radius under the condition of the wind speed 3m/s.Figures13 and 15 present that the bigger the wind speed, the smaller the attenuation rate of wave amplitude and the slower the delay rate of the phase. The reason is that the heat exchange time between air and surrounding rock is shortened when the wind speed increased, so the attenuation rate and delay rate at unit time is decreased.Figures14 and 16 present that the bigger the radius of roadway, the smaller the attenuation of wave amplitude and the slower the delay rate of the phase. Due to the fact that the wind speed is a constant, the rate of heat exchange is also assumed to be constant. Moreover, the ventilation quantity in unit time inevitably increases as the roadway radius enlarges, which lead to a decreased attenuation rate and delay rate at unit length.On the incorporation of the radial and axial temperature distribution, the synthesis analysis of temperature field in three-dimensional space can be obtained when the wind speed is 3m/s and the radius is 3m.Figure17 shows that the surrounding rock and air have the attenuation effect on the amplitude of temperature wave, and both present the law of negative exponent as the corresponding location increases, the overall shape of which resembles a folded page. Figure 18 shows that they both have the postponed effect on the phase and present the linear variation trend as the corresponding location increases, the overall shape of which resembles a parallelogram.Figure 17
Variation law of wave amplitude in three-dimensional space.Figure 18
Variation law of delayed phase in three-dimensional space.Moreover, compared with the air, the decay rate or delay rate of temperature wave within surrounding rock is relatively fast, which is affected by the dimensionless annual cycle time (thermal diffusivity of rock, radius of roadway, cycle time, and radial depth). However, in the axial direction, they are influenced by the location-related dimensionless time (thermal diffusivity of air, radius of roadway, and air velocity), the Nusselt number, the Biot number, and the dimensionless annual cycle time. Therefore, whether the decay of amplitude or the delay of phase, the effect of the rock is greater than that of the airflow.Figure19 shows that the ventilation distance has great influence on the amplitude and the phase of air temperature wave. The location-related cycle time along the axial direction can be calculated based on (33) and has close connection with the ventilation distance and velocity. In this instance, it is about 4.3 hours or the distance is 46600m. As mentioned previously, the ventilation rate in mine is relatively small and the maximum is 8m/s. Therefore, it is difficult to achieve a complete cycle fluctuation along the axial direction at one time.Figure 19
Variation law of axial temperature field in a time cycle.From Figure20, the temperature of air current may be greater than that of surrounding rock in the hottest season, and then the air could dissipate constantly heat to the surrounding rock. The other way round, the coldest season could present the opposite characteristic. That how to make good use of these features to govern the mine heat damage needs further discussion.Figure 20
Variation law of axial temperature field with the time.


## 5. Conclusions

In this work, the evolution mechanism of temperature field of coal roadway under seasonal fluctuation boundary is elaborated in detail, involving the radial and axial direction, to better implement the air temperature prediction in the underground high-temperature mine. The main conclusions can be drawn.The seasonal characteristics of the inlet air temperature can be regarded as a simple harmonic wave, and one dimensional unsteady heat conduction differential equation and its periodic boundary condition is established. The heat-balance equation is derived by finite volume method.The correctness of the temperature distribution function of surrounding rock of roadway is verified by the laboratory experiments. In addition, for the wall temperature wave, the greater wind speed, the larger wave amplitude but the smaller delayed phase; and they are almost not influenced by the cycle period. But in the interior, they appear to be the same attenuation or delay trend under the conditions of different wind speed; and at the same time the shorter the period, the greater the rate of decay and delay.On this basis, the mathematical model of temperature distribution along the axial direction is obtained from the heat-balance equation. The wind speed and roadway radius have great influence on the amplitude and the phase of temperature wave of air. The bigger the wind speed, the smaller the attenuation rate of wave amplitude and the slower the delay rate of phase. The bigger the roadway radius, the smaller the attenuation rate of wave amplitude and the slower the delay rate of phase.The radial and axial propagations of temperature wave have the similar characteristics of the amplitude decay and phase delay. For the wave amplitude, they both present the law of negative exponent as the location increases, the overall shape of which resembles a folded page. For the phase angle, they both present the linear variation trend as the location increases, the overall shape of which resembles a parallelogram.Compared with the air, the decay rate or delay rate of temperature wave within surrounding rock is relatively fast, which is affected by the dimensionless annual cycle time (thermal diffusivity of rock, radius of roadway, cycle time, and radial depth). However, in the axial direction, they are influenced by the location-related dimensionless time (thermal diffusivity of air, radius of roadway and air velocity), the Nusselt number, the Biot number, and the dimensionless annual cycle time. Therefore, whether the decay of amplitude or the delay of phase, the effect of the rock is greater than that of the airflow.In the paper, the seasonal variation of the airflow temperature in mine causes the problem of seasonal heat disaster; the principle analyzed scientifically can be guided engineering practice.



---

*Source: 2120618-2019-08-08.xml*