Abstract
Mining activities have increased owing to the rise in the social demand for minerals. Thermal hazards have become a major health and safety consideration in mines. The thermal environment of a working face is related to the air parameters at the bottom of shaft. The objective of this study is to accurately predict the air temperature at the bottom of a shaft in a mine with the ventilation time over 3 years. For this purpose, a mathematical model of the heat and mass exchange between the surrounding rock of the shaft and the downcast air is established by utilizing the finite volume method. The C++ languages are used for numerical calculations. The results are in agreement with the measured data. The effects of the relative humidity of air at the inlet of the shaft, the surface moisture coefficient of the shaft surface, and the physical parameters of the rock on the air parameters at the shaft bottom are studied in detail. Equations for calculating the enthalpy increase of air per 100 m in shaft with the depth of 1300 m were established by using cluster analysis. This equation provides a theoretical basis for predicting the air parameters along the shaft with the ventilation time over 3 years.
1. Introduction
The number of high-temperature mines is increasing with mining depth. According to statistics [1], there were almost 100 mines with a depth of more than 1000 m in China by 2015, among which approximately 47 were coal mines with an average mining depth of 1086 m. In the next 10 to 15 years, 50% of iron ore resources, 33% of nonferrous metal resources, and 53% of coal resources will be mined below 1000 m. The geothermal gradient distribution in China is high in the east and south and low in the west and north. The geothermal gradient within a depth of 1000 m (1.8–4.3°C/100 m) is linear and varies negligibly. The surface temperature at a depth of 1000 m is generally 40–45°C in northeast China and the southeast coast and 35–45°C in the central Sichuan Basin. The geothermal temperatures at a depth of 2000 m typically reach 70–80°C [2]. Therefore, thermal hazards are becoming an increasing concern in deep mines. It is crucial to effectively judge the thermal hazards of the working face of a mine and accurately predict the temperature and relative humidity of air.
The calculation of the heat and mass exchange between the surrounding rock of an air supply shaft and the airflow in the shaft can provide a basis for the prediction of the temperature and relative humidity of air on the working face.
Numerous studies have investigated the heat and mass transfer between airflow and excavation roadways and working surfaces. As early as 1919, South Africa began to study the thermodynamic law of mine air flow, and Jappe [3] puts forward the basic idea of the wind temperature calculation. Lambrechts [4] analyzed a large number of measured temperature changes in shafts with an average depth of 1258.8 m. They studied the relationship between the wet-bulb temperature of airflow in a roadway and natural rock temperature, initial wet-bulb temperature, distance, air velocity, and wall moisture content. Starfield [5] used the finite difference method to calculate the wall temperature of a roadway, and a formula was established using mathematical fitting to rapidly predict the dry-bulb and wet-bulb temperatures of airflow in the roadway. The formula can be used in practice. Ross [6] predicted the temperature, humidity, and other related climate parameters by numerical simulation. Lowndes and Hargreaves [7, 8] developed the ventilation and climate model of the fast tunneling roadway. Onder [9] studied the change of air supply volume and its influence on ventilation in the long-distance tunneling roadway. Kertikov [10] calculated the air flow temperature and humidity of the local ventilation drivage roadway. Sasmito [11] established a three-dimensional unsteady local thermal nonequilibrium model to evaluate thermal storage and the heat transfer between ventilation air and rock pit. In 2020, based on the similarity theory, Zhu [12] designed a facility to investigate effects of airflow velocity, thermal diffusion coefficient, moisture content, and thermal insulation material on the temperature field distribution in the surrounding rock. In order to save the energy and control mine environment, many scholars began to study how to use surrounding rocks or abandoned mines for thermal storage. Ghoreishi–Madiseh [13] and Zhu [14] studied the heat storage and heat regulation of surrounding rocks; Loredo [15] studied the heat storage ability of the isolated surveyed mine voids.
Most of the researches focus on the underground space, and there are few researches on the heat exchange process between the downcast air and shaft. In China, the method of differences and regression analysis, which were first proposed by Cen and Hou [16], was used to predict the air temperature at a shaft in 1989. At present, the simple difference iterative method [16], measured statistical method [4, 17], and PSO-BP neural network prediction method [18] are primarily used to calculate the air temperature at a shaft. However, these methods are too simple or geographically limited to estimate the air temperature accurately. And they neglect the latent heat effect of water vapor. In this work, according to the characteristics of coupled heat and mass exchange in the fluid–solid systems, the finite volume method is adopted to establish a mathematical model of the heat and water vapor exchange between the surrounding rock of an air supply shaft and the downcast air in the shaft. The formula that can accurately calculation the air parameters of the shaft bottom in different geological conditions and inlet air condition is given by using the C++ languages.
2. Mathematical Model
The bottom of the vertical shaft is an important node of the entire mine ventilation network. The temperature and relative humidity of the airflow in the shaft are related to the thermal environment of the entire mine. The heat and mass exchange between the shaft and the downcast air in the mine are different from those in a horizontal tunnel, and they have the following characteristics: (1) the airflow produces autocompression heat as it moves to the bottom of the shaft because of height difference. (2) The temperature of the shaft wall surface changes because of the geothermal gradient. (3) The shaft may pass through the underground aquifers. In this case, there is infiltration water on the shaft wall, and this water participates in the heat and mass exchange. (4) There may be heat dissipation by equipment and drainage pipes in the shaft.
The following assumptions are proposed to simplify the problem to be studied:(1)The air supply shaft is a hollow cylinder. The surrounding rock is isotropic and homogeneous, and the original temperature of the surrounding rock changes linearly with the depth(2)Only the heat and mass exchange between the downcast air and surrounding rock is considered, and the influence of infiltration water and other heat dissipation is neglected(3)The temperature and humidity of air are evenly distributed on the cross section of the shaft, and the change in the kinetic energy of the airflow is not considered
Based on these assumptions, the heat and mass exchange between the shaft and the airflow are shown in Figure 1. The equations of the heat and mass balance between the airflow and surrounding rock can be established as follows:
Equations (1) and (2) consider the surface of the shaft wall to be completely wet, but this is not the case. Therefore, the above equations must be revised. Japanese scholars first put forward the moisture coefficient of the wet wall rock [19], and most Chinese scholars have adopted the moisture coefficient to calculate the moisture exchange of the wet wall rock [20]. The value ranges from 0.05 to 0.4, in which the complete wet wall is 0.4 and the wet wall in one quarter (i.e., 0.1). After introducing the moisture coefficient of the roadway surface, equation (2) can be written as
2.1. Finite Volume Meshing
As shown in Figure 2, based on the assumption that the temperature and humidity of air are uniformly distributed on the cross section of the shaft in the axial direction, it can be considered that airflow parameters change only in the axial direction. The axial direction from the entrance of the shaft is set as in the meshing space.
As the temperature field of the surrounding rock is axisymmetric, the determination of the thermal conductivity of the surrounding rock can be simplified to a two-dimensional problem in the radial and axial directions. The finite difference grid is shown in Figure 2, where the meshing space is set as ∆r and ∆z.
2.2. Mathematical Model
2.2.1. Equations of Airflow Parameters
Air velocity and density are assumed to be 1 m/s and 1.2 kg/m3, respectively. The relative strength of convection and diffusion at adjacent nodes can be measured based on the grid Peclet number (Pe), which is defined as
The thermal conduct of air between adjacent nodes is negligible. According to the principle of energy and mass balance, the change in air temperature in the shaft is equal to the convective heat transfer between the downcast air and the surface of the surrounding rock and the autocompression. Additionally, the change in the humidity ratio of air is equal to the mass exchange between the airflow and the surface of the surrounding rock. Therefore, the equations of air parameters at different nodes are as follows:
The evaporation coefficient is based on the Lewis relation.
The equation can be applied to gases or liquids with ranging from 0.6–2500 and ranging from 0.6–100. can be approximately considered as when .
2.2.2. Equations for the Surrounding Rock
The differential equation of heat conduction in a two-dimensional cylindrical coordinate system is given by equation (8).
2.2.3. Boundary Conditions for the Surrounding Rock
The boundary conditions in the radial direction are as follows: the inner interface is the convective heat transfer interface between the surrounding rock surface and downcast air. This is regarded as the fourth type boundary condition. The outer surface is the surface of the cooling zone surface. It is a given temperature interface which is regarded as the first type boundary condition. The temperature is the original rock temperature.
The boundary conditions in the axial direction are as follows: the entrance of the shaft in the axial direction is set as the soil layer of constant temperature. Thus, the temperature at the section at can be regarded as a temperature in the soil layer of constant temperature. Based on literature [16], the thickness of the stable cooling zone is considered as approximately 10–20 m, and the rock temperature within the cooling zone is a power function distribution. The radial temperature gradient close to the shaft surface is considerably larger than the axial temperature gradient. The surrounding rock temperature rapidly approaches the original rock temperature as the radius increases. Therefore, the axial heat transfer can be neglected compared to the radial heat transfer at the shaft bottom. The bottom of the interface between the shaft wall and rock can be considered as an adiabatic surface. The mathematical model of the surrounding rock temperature is as follows:
The finite volume method is used to discretize differential equation (9). According to the energy conservation equation, it is assumed that the temperature at a node represents the temperature of the control volume for internal node (, ), which can be written as
This can also be written as
Its stability condition is .
2.2.4. Temperature of the Shaft Surface
The temperature of the shaft surface can be obtained using equation (5). Before solving equation (5), it is necessary to determine the density of saturated water vapor on the shaft surface, which is a single-value function of temperature. There are numerous empirical formulas for calculating the partial pressure of saturated water vapor; most of them are complex functional relations. To simplify the calculation, the relationship between the density of saturated water vapor and temperature between 10–45°C is expressed as the following linear equation:
The relative error between the calculated value calculated by equation (12) and the exact value is within 2.3%. Then, the surface temperature of the shaft can be calculated as
The convective heat transfer coefficient, , is calculated using the formula proposed by Sherbani [21].
3. Results and Discussion
3.1. Grid Independence Verification and Ventilation Time Selection
The geological conditions of a gold mine in Laizhou in the Shandong province in China are considered. According to the geological survey results of the Fourth Institute of Geology and Mineral Exploration of Shandong Province, the main composition of the surrounding rock is metagabbro when the depth is less than 1300 m. The thermal conductivity of metagabbro is 2.2 W/m°C [22].The temperature in the soil layer of constant temperature is 15°C, and the geothermal gradient is 2.2°C/100 m [23]. The radius of the air supply shaft is 2.5 m, and the air supply rate is 30 kg/s.
In order to select the appropriate space step size, mesh independence verification is carried out. As shown in Figures 3 and 4, grids are selected to calculate the air flow temperature in the shaft, when the step sizes are considered as and , and the air flow curve almost coincides. The variation of the air temperature in the inlet shaft with the ventilation time is also calculated. As shown in Figure 4, when the inlet air temperature is 15°C, the air temperature along the shaft gradually decreases with the increase of the ventilation time and tends to be consistent. The temperature curves of the ventilation time at 2 and 3 years are almost coincident. At the beginning of ventilation, the shaft wall temperature is higher, and the heat releases into the air are higher. As the ventilation time increases, the shaft wall is gradually cooled, and the heat release to the air is reduced, so the airflow temperature along the shaft decreases with time. And the radius of the heat adjusting circle expands continuously, but the radius expands more slowly with the time increases and approaches to the steady state gradually. So, the temperature curves of the ventilation time at 2 and 3 years are almost coincident. In the following calculation, the step sizes are considered as and , and the ventilation time is considered 3 years.
3.2. Validation of the Model
The comparison of the calculation results and the measured data obtained from the literature [24] is shown in Figure 5. The calculation condition in literature [24] is as follows: the constant temperature layer with a value of 15.7°C; the mean value of the geothermal gradient is 2.91°C/100 m. The shaft radius is 3.5 m, and the speed of the airflow is 3 m/s. The physical properties of the wall rock are as follows: the underground 0-99 m is a mudstone layer with a thermal conductivity of 2.1 W/(m·°C); the sandstone layer is distributed at 100-500 m underground, and its thermal conductivity is 4.1 W/(m·°C); the limestone layer is distributed at 501-666 m underground, and its thermal conductivity is 3.1 W/(m·°C). The calculated values exhibit strong agreement with the measured values. Therefore, the results obtained using the model are reliable. Thus, the model is used to calculate the air parameters at the shaft bottom under different inlet air parameters and different thermophysical characteristics of the rock. The results are described below.
3.3. Effect of Temperature and Relative Humidity of Inlet Air on Temperature and Enthalpy of Air along the Intake Shaft
As in literature [24], in this work, it is assumed that is 0.05, the inlet air temperature is 10–25°C, and relative humidity is 50–80%. As shown in Figures 6 and 7, the temperature and relative humidity of air in the shaft are evidently correlated with the inlet air temperature. The variation in the air temperature in the shaft depends on the heat dissipation of the surrounding rock and the autocompression.
The heat dissipation of the surrounding rock is composed of the sensible heat transfer and latent heat transfer between the air and surrounding rock surface. The sensible heat transfer may be negative when the inlet air temperature is higher than the temperature of the surrounding rock surface. The air dry-bulb temperature decreases when the absolute value of the sensible heat transfer is more than the autocompression. The air dry-bulb temperature starts increasing when the value of the sensible heat transfer is positive or the absolute value of the sensible heat transfer is less than the autocompression. However, the enthalpy of the airflow along the intake shaft increases because of air self-compression and wall heat dissipation.
3.4. Effect of Temperature and Relative Humidity of Inlet Air on the Air Temperature at the Shaft Bottom
The working efficiency and thermal comfort of workers depend on the air parameters on the working face, which are related to the air parameters at the bottom of the shaft. The influence of the temperature and relative humidity of inlet air on the air temperature at the shaft bottom is analyzed when , as shown in Figure 8. The air temperature at the shaft bottom increases with the temperature and relative humidity of inlet air. As relative humidity increases, the evaporation of water on the intake shaft wall decreases and the latent heat decreases. This increases the air temperature at the shaft bottom. Figure 6 shows that the relative humidity of inlet air significantly influences the airflow temperature along the shaft. However, as the depth reaches 600 m, the air temperature in the shaft increases linearly with the shaft depth.
3.5. Effect of the Surface Moisture Coefficient on Air Parameters at the Shaft Bottom
In addition to the relative humidity of inlet air, the surface moisture coefficient of the wall surface is related to the amount of water vapor that enters the airflow from the wall. Therefore, the surface moisture coefficient of the wall surface strongly affects the air temperature at the shaft bottom. The inlet air at 15°C is considered as an example to calculate the change in the air temperature at the shaft bottom with varying from 0–0.3, as shown in Figure 9. There is no water evaporation when the surface moisture coefficient is 0, and the heat transferred from the surrounding rock is converted into the sensible heat of air. Under the influence of rock heat dissipation and autocompression, the air temperature at the shaft bottom is up to 28.73°C, and it is not affected by the relative humidity of inlet air. As the surface moisture coefficient increases, the air temperature at the shaft bottom decreases. The air temperature at the shaft bottom does not change significantly after the surface moisture coefficient becomes more than 0.2.
The effect of the surface moisture coefficient on the air enthalpy at the shaft bottom is analyzed, as shown in Figure 9. When the inlet air temperature is 15°C, the air enthalpy at the shaft bottom increases slightly from the dry wall surface to the wet wall surface. However, the change is not as evident as the increase in the surface moisture coefficient. The air enthalpy at the shaft bottom changes negligibly after the surface moisture coefficient becomes greater than or equal to 0.1.
Heat conduction and transfer are expressed in cylinder coordinates by the following equations:
is in a range of 1.36–2.34 when the inlet air temperature is 15°C. However, it has been reported in literature [19] that the evaporation heat transfer of the actual underground roadway is approximately 0.2 to 1.57 times the sensible heat of the completely dry roadway surface. This implies that the range of is 1.2–2.57. Thus, the results of our study agree with literature. The effect of on heat transfer is not evident in the case of . The main reason for the increase in heat transfer is that the air temperature along the shaft decreased as the latent heat of evaporation increases. Thus, the change in the air enthalpy at the shaft bottom is negligibly affected by the surface moisture coefficient.
3.6. Effect of Air Velocity on the Air Temperature at the Shaft Bottom
The relative humidity of mine supply air is generally high. Therefore, the effect of the air velocity on the air temperature at the shaft bottom is calculated at a relative humidity of 80% and different temperatures when , as shown in Figure 10. The increase in the air velocity significantly affects the air temperature at the shaft bottom when the air velocity is small. This effect is stronger when the inlet air temperature is low and vice versa. When the inlet air temperature is 25°C, the air velocity negligibly affects the air temperature at the shaft bottom. When the inlet air velocity increases to 1.70 m/s (air volume: 40 kg/s), the air temperature at the shaft bottom does not change significantly with the air velocity at all temperatures. Therefore, when the inlet air temperature is low, the air temperature at the shaft bottom decreases considerably as the air velocity increases. However, when the inlet air temperature is high, air velocity does not affect the air temperature at the shaft bottom.
3.7. Effect of Rock Thermophysical Properties on the Air Temperature at the Shaft Bottom
The temperature at the same depth of the surrounding rock increases with the geothermal gradient. Additionally, the thermal conductivity of the rock influences rock heat dissipation. Inlet air temperature, relative humidity, and the surface moisture coefficient are set as 15°C, 80%, and, and the air supply rate is 30 kg/s to calculate. The results are shown in Figure 11. When the geothermal gradient increases from 2°C/100 m to 4°C/100 m, the heat transfer from the surface of the shaft increases by about 2.28 times. When the thermal conductivity of the rock increases from 1 W/m·°C to 4 W/m·°C, the heat transfer from the surface of the shaft increases by about 3.6 times. Moreover, the effect of the thermal conductivity of the rock is stronger in rocks with a larger geothermal gradient. As the thermal conductivity of the rock increases from 1 W/m·°C to 4 W/m·°C, the heat transfer from the surface of the shaft increases by 75.8KW and 173.3KW when the geothermal gradient is 2°C/100 m and 4°C/100 m, respectively. The geothermal gradient and thermal conductivity of the rock have similar effects on the air temperature at the shaft bottom.
3.8. Comprehensive Effect and Equation Fitting
The effect of the abovementioned factors on the air parameters at the shaft bottom is further examined. The previous analysis results show that when the wind speed exceeds 1.7 m/s, the airflow parameter changes very little, while the wind speed in the inlet shaft often exceeds 3 m/s [4, 24], so the wind speed is set as 3 m/s in the calculation. According to the literature [3], the main geothermal gradient of China is 1.8-4.38°C/100 m, and the rock thermal conductivity is 0.2-4.2 w/m·°C [16], so the calculated parameters are set as shown in Table 1. And the radius of the shaft is set as 2.5 m. C++ languages is used to calculate the temperature and enthalpy of air along the shaft under the conditions listed in Table 1. As shown in Table 2, the equations for calculating the difference enthalpy per 100 m of air along the shaft of the depth less than 1300 m are obtained using the clustering analysis method.
According to Figure 9, although the surface moisture coefficient can affect the air temperature, it has little influence on air enthalpy. Therefore, the equations can be combined into a unified equation when the surface moisture coefficients are different. It can also be seen from the equation that in the vertical shaft excluding other heat sources, the influence of autocompression is dominant, and the heat transfer of the surrounding rock only has significant influence when the geothermal gradient and rock thermal conductivity are both high. With the increase of the depth, the original rock temperature of the rock increases so in the deep mine, the heat dissipated in the surrounding rock of the shaft will also increase gradually.
The equations in Table 2 show the difference enthalpy per 100 m along shaft. When the shaft radius is 2.5 m, the error of the enthalpy difference calculated by the fitting formula and the program under the dry condition is within 2%, and the error under the wet condition is within 8%.The enthalpy of the flow at the bottom of the inlet parameter in literature [24] is also calculated by formula in Table 2. The calculated results are in good agreement with the measured value in literature [24] as shown in Figure 12.The results of the fitting equations can provide a reference for designers to estimate the air parameters at the bottom of the intake shaft in different.
4. Conclusions
A mathematical model is built according to the heat and the mass exchange principle between a shaft and the downcast air with the ventilation time of 3 years by employing the finite volume method. C++ languages is utilized to calculate the heat and mass exchange. The calculated results are in agreement with the measured data in literature [24]. Equations for calculating the enthalpy increase of air per 100 m in shaft of the depth less than 1300 m were established by using cluster analysis. This equation provides a theoretical basis for predicting the air parameters at the shaft bottom with the ventilation time over 3 years.
The main conclusions of this study are as follows:(1)The parameters of inlet air affect the parameters of downcast air along the intake shaft. The influence of relative humidity on the air temperature along the intake shaft is stronger at a high inlet air temperature(2)There is an evident positive relationship between the air temperature at the shaft bottom and the temperature and relative humidity of inlet air. The surface moisture coefficient affects the air temperature at the shaft bottom. However, the air temperature at the shaft bottom does not changes significantly when the surface moisture coefficient becomes greater than or equal to 0.2. The surface moisture coefficient negligibly affects the air enthalpy at the shaft bottom(3)The increase in air velocity evidently affects the air temperature at the shaft bottom when the inlet air temperature is low. However, when air velocity increases to 1.70 m/s, the air temperature at the shaft bottom does not change significantly as air velocity increases(4)The air temperature at the shaft bottom and the heat released from the surface of the shaft increase with the geothermal gradient and thermal conductivity of the rock
Nomenclature
| : | area, m2 |
| : | specific heat of air, J/kg·K |
| : | air moisture, g/kg(a) |
| : | moisture coefficient |
| : | area of the roadway section, m2 |
| : | gravitational acceleration, m/s2 |
| : | air mass flow of the shaft, kg/s |
| : | convective heat transfer coefficient, w/m2·K |
| : | Evaporation coefficient, kg/m2·s |
| : | the enthalpy of air, kJ/kg(a) |
| : | the enthalpy of saturated air at original rock temperature, kJ/kg(a) |
| : | Lewis number |
| : | Peclet number |
| : | Prandtl number |
| : | total heat exchange, w |
| : | significant heat exchange, w |
| : | latent heat exchange, w |
| : | latent heat of vaporization of water, kJ/kg |
| : | shaft radius, m |
| : | cooled zone radius, m |
| : | Schmidt number |
| : | air temperature, °C |
| : | rock temperature, °C |
| : | original rock temperature, °C |
| : | temperature in soil layer of constant temperature, °C |
| : | geothermal gradient, °C/100 m |
| : | the difference enthalpy of air per 100 m, kJ/100 m |
| : | circumference of the roadway, m |
| : | air speed, m/s |
| : | depth of the shaft, m. |
| : | surrounding rock thermal diffusivity, m2/s |
| : | thermal conductivity of rock, W/m·°C |
| : | thermal conductivity of air, W/m·°C |
| : | roughness coefficient of roadway |
| : | relative humidity, % |
| : | air density, kg/m3 |
| : | moisture separation coefficient |
| : | time, s. |
| in: | inlet air |
| out: | outlet air |
| : | the surrounding rock surface. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (No. 51134005 and No. 52174098) and Human Engineering Research Centre of Energy Efficiency and Material Technology in green and low-carbon building.