Finite Difference Analysis of Transient Heat Transfer in Surrounding Rock Mass of High Geothermal Roadway

Based on finite difference method, a mathematical model and a numerical model written by Fortran language were established in the paper. Then a series of experiments were conducted to figure out the evolution law of temperature field in high geothermal roadway. Research results indicate that temperature disturbance range increases gradually as the unsteady heat conduction goes on and it presents power function relationship with dimensionless time. Based on the case analysis, there is no distinct expansion of temperature disturbance range after four years of ventilation, when the temperature disturbance range R = 13.6.


Introduction
At present, solid mineral resource exploitation can reach 4 km underground [1]. With the rapid development of China's economy, the production and consumption of solid mineral resources increase year by year, resulting in the gradual exhaustion of shallow resources and escalation of mining depth. Taking coal as an example, in the past decade, the annual increasing rate of coal production was 8.6%, and the yield reached 366.9 billion t at 2013 [2]. The shallow minerals were exhausted at the majority of mining districts because of long-time, high-intensity, and large-scale mining activities, and the mining deepens downward 8 to 12 m every year [3]. Up till 2013, there were about 200 state-owned coal mines exploiting 800 m beneath the ground, among which 47 coal mines conducted their mining activities under −1000 m. Most of these coal mines are located in mideastern China such as Shandong, Jiangsu, and Anhui province. The average mining depth is −1086 m and maximum value is −1501 m in Suncun coal mine of Xinwen Mining Group, Shandong province [4]. It is predicted that in east, middle, north, and northeast China, most coal mines will conduct their deep exploiting activities from −1000 m to −1500 m by 2020.
With increasing mining depth, heat hazard usually takes place and becomes more and more prominent. High geothermal temperature is a decisive cause of heat hazard in deep mining. Generally, the geothermal temperature gradient of China's coal seams is 2.5∼3 ∘ C/hm, but at abnormal geothermal areas it can reach 4∼5 ∘ C/hm. It is thus estimated that, with the depth and temperature of constant temperature zone of 30 m and 15 ∘ C, respectively, the original rock temperature can reach 40∼45 ∘ C at −1000 m [5]. Actually, the rock temperature at −1000 m of many heat hazard mines is close to or exceeds the upper limit. For example, in Sanhejian mine of Xuzhou Mining Group the original rock temperature can reach 46.8 ∘ C, in Suncun mine of Xinwen Mining Group it is 45 ∘ C, and in Zhaolou mine of Yanzhou Mining Group it is 43 ∘ C [6]. Heat hazard is prone to occur under the combined influence of high geothermal temperature, heat of air compression, and heat dissipation of machine. According to incomplete statistics, mining activities in about 150 coal mines are affected by heat hazard in China, and their workface temperature is usually in the range of 30∼40 ∘ C. Heat has already become the sixth mine hazard after roof falling, gas, water inrush, fire, and mine dust. Beside deteriorating the air environment of mines, high geothermal temperature 2 Mathematical Problems in Engineering can change the physical and mechanical properties of rocks and further affect the stability of surrounding rock and heat conduction process [6][7][8][9].
Temperature that is not disturbed by human engineering activities is named original rock temperature [10,11]. Once the roadway excavation of high geothermal mine was conducted, convective heat transfers between surrounding rock and low-temperature airflow and instantly breaks the balanced temperature distribution in strata; the geothermal energy constantly transferring to airflow may result in underground air temperature increase. Meanwhile, the surrounding rock temperature decreases from surface to deeper area by and by, and the cooling scale of surrounding rock continuously expands, which is considered as temperature disturbance range (TDR). The rock inside TDR is called heat regulating circle (HRC). Because of the heat resistance in the HRC, the temperature decreasing amplitude declines with the distance going deeper until reaching a temperature boundary which is relatively unchanging, and its value is similarly equal to the original rock temperature. Instead of losing huge amount of geothermal heat quickly in a short time, the HRC can release the heat into the airflow little by little so as to buffer the heat dissipation. In this perspective, research on temperature distribution of HRC has important theoretical and practical significance.
Temperature distribution of the surrounding rock can be measured by thermocouples in holes drilling, and this shortcut method has been successfully applied in Huaibei and Pingdingshan Mining Groups in China. The field measure process indicates that, with the ventilation moves on, the surrounding rock temperature drops constantly with a smaller and smaller decreasing amplitude [12,13]. Other scholars measured the original rock temperature by drilling holes deep or shallow. Besides, with the development of computer technology, numerical method has also been applied in temperature field studies of surrounding rocks; for instance, by numerical computation, Zhang concluded that, under the coupling effect of airflow and seepage, the symmetrical distribution characteristics of temperature field are influenced by seepage field rather than the airflow [14]. Meanwhile, Zhang et al. suggested that the temperature distribution fits exponential function in radical direction; the temperature gradient and heat flux at roadway surface are maximum; the relation between surrounding rock temperature and ventilation time similarly meets Hill function, and the fitting formula of the HRC was acquired [15]. According to Fan and Zhang, the scale of HRC was calculated by mathematical model and its asynchronous semiexplicit difference form; however, because of the specific research purpose, general rules of the TDR were not summarized [16]. The research of [17] points out that thermal diffusivity greatly affects both the TDR and temperature distribution of surrounding rock, while roadway geometry size and convective heat transfer coefficient of surrounding rock and airflow have little impact on the TDR. Though the characteristics and evolution laws of roadway temperature distribution have been fully reported in existing literatures, topics concerned about TDR, especially the general rules of the HRC expansion, still demand further studies.
In this paper, the authors established a mathematical model to reveal the expansion laws of the HRC based on the FDM and summarized the general rules of the TDR.

Mathematical Model of Roadway Temperature Field
To facilitate the evolution laws study of roadway temperature distribution, the following assumptions were proposed at first [18]: (1) The roadway cross section is round ( = 0 ) and its surrounding rock is a hollow cylinder with infinite radius ( → ∞).
(2) The surrounding rock is continuous, homogeneous, and isotropous, and its thermophysical parameters are constants.
(3) Initially when = 0, the surrounding rock temperature at each point is uniform and equal to its original rock temperature ( 0 ).
(4) The airflow temperature only changes along the axial direction and remains uniform on the same radical cross section. The heat transfer conditions are identical at the circumferential direction.
Meanwhile, it is assumed that there is an ideal heat insulation layer with infinitesimal thickness, infinitely great external diameter, and an inner diameter equal to roadway diameter; it also has the properties of homogeneity, continuity, and isotropy and heat insulation. In this case, if a roadway with length of was divided into infinite segments by the above layer and assured that inner bore of ideal heat insulation layer was overlapped to roadway cross section, then each segment of the roadway surrounding rock would be infinitely thin heat conduction slice without heat transfer between adjacent slices. Consequently, the heat would only transfer along radical direction and the conduction would become one-dimensional. These assumptions neglect heat transfer in axial direction in the surrounding rock mass. In fact, the heat flux in axial direction is negligibly small, especially when the rock mass temperature has not been disturbed on a large scale.
Based on the above assumptions, the law of conservation of energy, and Fourier's law, taking a microunit from the assumed surrounding rock slice, differential equations (mathematical model) of heat conduction under onedimensional cylindrical coordinate system were established, as shown in Figure 1 and (1). One has where = / stands for the thermal diffusivity of roadway surrounding rock; stands for the surrounding rock temperature; stands for time; stands for the radical coordinate of surrounding rock.
Derived from law of conservation of energy and Fourier's law, (1)   describe the heat conduction phenomenon of surrounding rock, the following monodromy conditions should be limited: (1) Geometrical condition: it is semi-infinite plate with inner radius of 0 .
(2) Physical condition: rock density , specific heat capacity , and heat conductivity are constant and free from the influence of temperature variation; no internal heat source exists in the isotropous surrounding rock.
(4) Boundary condition: the boundary condition of surrounding rock heat conduction reflects the interaction between surrounding rock and its ambience during the process of heat conduction; it reflects the characteristics of heat conduction on surrounding rock boundary. As a semi-infinite hollow round slice of infinitely great diameter, the temperature at the area infinitely far away from the circle center tends to approach the original rock temperature unlimitedly; that is The other boundary of roadway surrounding rock is the surface, where = 0 . It is difficult to acquire the temperature function (the first boundary condition) due to complex boundary conditions and complicated process of convective heat transfer between surrounding rock and airflow at roadway surface. Many scholars assumed that the rock temperature is constant and unchanging at roadway surface, which was obviously appropriate. If the heat flux of roadway surface was given (the secondary boundary condition), all problems would be settled easily, but either rock temperature field or surface heat flux coefficient function should be known as a premise. The third boundary condition is most extensively utilized because of the continuous convective heat conduction between surrounding rock surface and airflow. For convective heat transfer under steady airflow condition and temperature, the surface convective heat transfer coefficient ℎ can be seen as constant and the function of airflow temperature with time and space is considered as known; then where is the exterior normal of roadway surface. The surface rock temperature and temperature gradient ( / ) are unknown. The following mathematical model for surrounding rock temperature field is established as a consequence:

FDM of Surrounding Rock Temperature
where temperature of regional nodes ( Δ , Δ ) is denoted as .

Mathematical Problems in Engineering
Unsteady unit / in (4) applied backward difference scheme while / and 2 / 2 adopted central difference scheme; then where (Δ ) which represents the bottom step of the remainder is one order and can be omitted. Introducing (6) into (5) gives which can also be written as Introducing , , and , then Introducing (9) into (8) gives The matrix form is ] .
Equation (12) is a tridiagonal matrix, which can be solved by Thomas Algorithm method. Discrete initial conditions give In terms of boundary conditions, the left boundary meets the third boundary condition while the right boundary meets the first boundary condition. When = 0 , in order to ensure the discrete precision uniformity with extensive function, one order central difference is applied to the third boundary condition ( / ) = ℎ( − ) and it gives Combining (9) and (14) gives By using the same analysis method, the first boundary condition of the left boundary = 0 will be discretized as

Stability of Difference Scheme.
Because of marching method is involved to solve the FDM of unsteady heat conduction, if deviation was introduced from step , then deviation in step +1 would be inevitably produced and even increased. If this deviation could be controlled, disappears gradually, or has a boundary, then the established FDM is steady. In this paper, Fourier series is adopted to assure the stability of the difference scheme. Setting +1 = V +1 exp ( ℎ) and introducing it to (9) give Excessive divisor is acquired as Setting 1 = Δ /2 Δ and 2 = Δ /Δ 2 and combining (18) give To all 1 and 2 , | ( , )| ≤ 1, which meets Von Neumann condition, the difference scheme equation (8) is steady absolutely.
Based on the above difference scheme, a numerical program for unsteady temperature field of roadway surrounding rock is written by Fortran language.

Numerical Test Program
According to the correspondent mathematical model, simulated roadway has a round cross section, with radius 0 of 2 m which is determined by the conversion of hydraulic diameter. Generally, the thickness of HRC would not exceed 40 m, so the external boundary is set as 60 m after proper expansion. Besides, as the thermal environment of simulated roadway is incapable of resulting in large variation of the thermophysical properties of surrounding rock, rock density , specific heat capacity , and heat conductivity can be regarded as constants, where = 1.2 W/m⋅K and = 0.8 × 10 −6 m 2 /s ( = / ).
From the assumptions of surrounding rock heat conduction models, the inner boundary matches the third boundary condition, while the external boundary belongs to the first boundary condition. The original rock temperature is set as 37 ∘ C. In general, the convective heat transfer coefficient ℎ of turbulent airflow condition is 20-100 W/m 2 ⋅K, and it is rarely affected by the airflow temperature because of its slight temperature variation in roadways. The value of ℎ can be calculated by following equations: where is the airflow heat conductivity coefficient, 0.026 W/m⋅K; is the surrounding rock heat conductivity coefficient, 1.2 W/m⋅K; is the coefficient of thermal diffusivity, 0.8 × 10 −6 m 2 /s; ] is the airflow kinematic viscosity coefficient, 14.4 × 10 −6 m 2 /s; and is the hydraulic diameter, 4 m.
In the simulation process, the airflow temperature is constantly 20 ∘ C and the external boundary temperature of roadway surrounding rock is equal to the original rock temperature; that is, = 0 . Based on the above assumptions, a one-dimensional hollow cylinder slice model with inner radius of 2 m and external radius of 62 m was built up. To shorten the simulation test time, the model is divided into 600 grids by equal step length of 0.1 m. The simulation time is 3 years long.
The test program is illustrated in Table 1. Figure 3 presents the cloud picture of roadway temperature field in 30000 h. It is obvious that the TDR extends continuously with the progress of ventilation and reaches its maximum at about 25 m, but detailed characteristics cannot be observed from the figure. To detect the evolution laws of rock temperature field in detail, the  experiment data was rehandled. To obtain the universal law between TDR, space, and time, dimensionless method was applied to the above factors, where the dimensionless temperature Θ = ( − )/( 0 − ), the dimensionless radius = / 0 , and the dimensionless time = / 2 0 . The relationship between TDR and its corresponding dimensionless time is illustrated in Figure 4. Figure 4 elaborates that, with the progress of unsteady heat conductivity, TDR extends gradually with a declining rate (the range boundary is lower than 1% of original rock temperature). It is clear that the TDR presents exponential function with dimensionless time = 2.91 0.46 + 0.98.

Numerical Test Results.
Certain deviation exists in (22) because theoretically when is equal to 0, should be 1, and the temperature disturbance depth is 0. The deviation may be produced by calculation precision error. For example, if the step length is set as 0.1 m when dividing the spatial grids, experimental data only exists at the nodes integral times of 0.1 m, which would cause a deviation of 0.1 m in TDR statistics.  In order to acquire a more practical formula for TDR calculation, the fit function was defined as follows: Introducing (23) to fit the data in Figure 4 which gives = 2.75 can get a goodness of 0.9988 and a standard error of 0.02, which can meet the needs of engineering application. Introducing = / 2 0 to (23) gives where ( − 1) 0 is the temperature disturbance depth, m; is the thermal diffusivity coefficient of surrounding rock, m 2 /s; is the heat conduction time, s. To calculate the maximum range of temperature disturbance, the paper set 6 months or so as an interval ( = 3) and plotted the relation curve between dimensionless temperature and depth of surrounding rock in Figure 5.
It is clear that the rock temperature changes distinctly during the first 3 years, then this variation gradually diminishes. Four years later when the TDR is about 13.6 (25.2 m), there is no apparent expansion. The value is similar to the theoretical value calculated by (24) (14.47), which provides evidence for the steady temperature field calculation of roadway surrounding rock and drilling detection of original rock temperature.

Discussions on HRC.
HRC is a key concept in studying the roadway temperature field. Initially, HRC was proposed in view of the effect of periodic variation of surface air temperature on surrounding rock temperature distribution when the surface air was supplied into the underground. In opposite, low-temperature airflow would be heated by roadway surrounding rock temperature, and high-temperature airflow would be cooled. In a narrow sense, HRC is the surrounding rock affected by periodical changing air temperature. Afterwards, this concept was also introduced to the research of cooling range of high geothermal roadway. At present, HRC has no uniform definition. Some scholars hold the opinion that the HRC is a static concept, which is the maximum scale of surrounding rock whose temperature has an apparent variation after the heat conduction of roadway surrounding rock is basically steady, while others think that the HRC is a dynamic concept, and it would extend constantly with the continuous progress of heat transfer. In this paper, the authors adopt the latter opinion. Strictly, the basis of the HRC is a temperature disturbance scale affected by ventilation. Therefore, this paper defines the TDR as a scale with temperature difference more than 1% of original rock temperature in the progress of unsteady heat transfer.
Studies of the roadway temperature disturbance range have realistic significance. For example, the depth of drilling hole for measuring rock temperature should be larger than the TDR when conducting in situ measurement of the temperature field. In this case, the external boundary is the maximum value of the TDR because of small variation of rock temperature after long-time ventilation. So the rock temperature field at this time can be simplified as steady. The calculation results in this paper present exponential function between TDR and ventilation time. Seriously, there is no fixed boundary of the temperature disturbance district because the TDR will increase with the progress of ventilation. However, as the order of the power function between the TDR and ventilation is less than 1, with longer ventilation time, the increasing rate of the TDR declines continuously. In that way, there is a relatively steady area of the temperature disturbance district and it can be called "maximum temperature disturbance range" after long-time ventilation. In fact, in mining practice, the judgment standard for the TDR is to seek a point that its temperature keeps on being less than 1% of the original rock temperature in an allowed long interval, rather than to find the external boundary when the surrounding rock temperature is steady.

Conclusions
In this paper, a Fortran language source procedure for surrounding rock temperature field calculation and the corresponding evolution laws was written on the basis of onedimensional heat conductivity mathematical model and the FDM. The results show that the TDR of roadway surrounding rock extends constantly with the progress of ventilation. Meanwhile, it presents power function with dimensionless time of unsteady heat conduction. The findings will provide evidence for original rock temperature measurement in underground mines.