Numerical study on calculation model of airflow temperature and humidity prediction in humid airway

In deep coal mining, the high temperature and humidity environment are one of the main disasters threatening workers and equipment. To predict the temperature and humidity of airflow in the airway accurately, this study establishes a numerical calculation model of the surrounding rock temperature, airflow temperature, and humidity under the conditions of considering the moisture evaporation and heat absorption of the airway wall, the condensation and heat release of supersaturated water vapor in the airflow, and the influence of air compression or expansion. Combined with the two‐dimensional (2D) heat conduction model (considering the influence of the axial surrounding rock heat conduction of the airway), the axial strip partial moisture calculation model of the airway is proposed, and then the uniform moist wall model and the axial strip partial moist wall model are established. The influence of different parameters on the temperature and humidity of airflow is discussed by analyzing the results of numerical simulation of the uniform wet surface model (UWSM) and partly wet surface model (PWSM) of the airway wall. The results show that the surrounding rock temperature, especially the temperature distribution of the airway wall, calculated by using the 2D heat conduction model under the local humid airway conditions, is more realistic. The increment of airflow temperature decreases with the increase of wet area ratio, airflow temperature at the entrance and airway radius, and increases with relative humidity at the entrance. The calculation accuracy of temperature and humidity prediction for large ventilation networks can be further improved by using the equivalent humidity of the UWSM to replace the correction coefficient of the wet area ratio of the PWSM.


| INTRODUCTION
Coal resources are the main component of China's energy structure. With industry development, shallow coal resources are increasingly exhausted, and the normalization of deep mining is inevitable. With the increase of mining depth, the influence of surrounding rock on coal mining also increases. Generally (not affected by special geothermal environment), the average temperature of surrounding rock will increase by 3°C for every 100 m increase in the depth below the constant temperature zone. The high temperature and high humidity working environment not only affects the working efficiency of operators but also poses a serious threat to their physical and mental health and life safety. 1 Therefore, in the use of deep underground space or mining, to achieve the purpose of controlling the environment, it is necessary to study the change law of underground airflow temperature and humidity and its influencing factors, and accurately predict them, to carry out reasonable ventilation planning, achieve the control of airflow temperature and humidity, which is of great significance to solve the mine heat disaster.
The research on the heat conduction of surrounding rock and the heat & moisture exchange between surrounding rock and airflow started earlier at home and abroad. McPherson 2 summarized the time change of heat budget in mine airways and the basic theory of airway wall water evaporation. Danko and Bahrami 3 used software (Multiflux and CLIMSIM) to simulate the heat and moisture transfer between roadway surrounding rock and airflow. Cygankiewicz and Knechtel 4 used the computer simulation method to study the influence of surrounding rock on airflow temperature and enthalpy and the influence of original rock temperature and rock mass temperature on heat dissipation of surrounding rock. Gao and Zhang 5 proposed a quadratic curve fitting formula for the relationship between saturated air moisture content and temperature. Zhang et al. 6 extended the starfield prediction formula and applied it to predict airflow temperature in dry airways. Zhang et al. 7 established the mathematical model of unsteady heat dissipation of surrounding rock based on the twodimensional (2D) heat conduction model and explored the factors affecting the heat conduction of surrounding rock. Then he concluded that the lithology of surrounding rock, ventilation speed, ventilation time and equivalent diameter of roadway are important factors affecting the unsteady heat dissipation of surrounding rock.
Based on the numerical model of surrounding rock and airflow established, some scholars adopted different methods to process and analyze the numerical calculation model, and obtained some simplified formulas and dimensionless numerical models. Song et al. 8 established the unsteady 2D heat conduction model of surrounding rock, analyzed the governing equations of heat conduction of surrounding rock in the cylindrical coordinate system and rectangular coordinate system using dimensionless criteria, and deduced the calculation method of heat dissipation of surrounding rock. Qin et al. 9 established a dimensionless 2D heat conduction model of surrounding rock by introducing dimensionless parameters such as Biot number and Prandtl number, and discretized the dimensionless numerical model by using the finite volume method. Yi et al. 10 revealed the seasonal air temperature affected the airflow temperatures and observed the slight temperature variations on the return route, then get the radius of the heat-adjusting layer was between 28 and 33 m. Zhang and Yang 11 analyzed the sensible heat, latent heat and total heat loss by the influencing of ventilation time and found that the wetness factor has a strong effect on surrounding rock heat release.
However, in the actual application process, the commonly used calculation programs have large errors in the prediction of mine climate. 12 Especially for some partly wet airways, there is a lack of a perfect calculation model. Some studies put forward two basic concepts of humidity and wet area ratio in the study of the airway. [13][14][15] In the numerical simulation, humidity (relative to the effect of a completely wet airway wall) is a common factor, 16 which is defined as the ratio of the influence effect on the airway wall relative to the complete wet. This definition is a qualitative concept and cannot be given according to the physical parameters of the airway. However, the wet area ratio is defined as the ratio of the wet area of the airway. Humidity and wet area ratio are different concepts, but they have not yet been defined as quantitative differences. In the research on the wall model of the partially wet airway, the most representative is the research on the influence of Starfield and Uchino on the zonal wet part of the circumferential wall of the vertical z-axis airway section and put forward the theory of equivalent humidity in different parts of the zonal wet airway. However, there are few reports on the study of the equivalent humidity of some wet walls in the axial direction (airflow direction).
In this study, the axial surrounding rock heat flow component of the airway is incorporated into the calculation model. Considering the initial surrounding rock temperature, air inlet temperature and humidity, roadway inclination, roadway wall humidity, ventilation time and mutual conversion of water vapor, the finite difference method is employed to simulate and calculate the temperature and humidity of the surrounding rock airflow. The quantitative relationship between the humidity and the wet area ratio under the humid environment (axial strip humid wall) is studied, and the numerical calculation model is applied to the climate prediction of some humid environments to improve the accuracy of predicting the temperature and humidity of the humid roadway airflow.

| Heat conduction equation of surrounding rock
The thermal characteristics of the surrounding rock are affected by the direction, stress, local groundwater movement and rock types and change slowly with the change of temperature. In the study of the underground coal mine climate, it is assumed that the thermal properties of the surrounding rocks involved remain unchanged, and the airway section is circular. The relationship between the airway and surrounding rock and the schematic diagram of the numerical difference model in surrounding rock is shown in Figure 1. The coordinates of surrounding rock location are defined as z r γ ( , , ), the z-axis is the axial length of the roadway (m), r is the radial distance from the z-axis (m), γ is the circumference angle (rad), R is the roadway radius (m), g (m/s 2 ) is the gravitational acceleration, Ѱ IN is the relative humidity at the entrance of airway, Θ IN is the airflow temperature at the entrance of airway (°C), Θ j and χ j are the airflow temperature (°C) and moisture content (kg/kg) in section j, ΔΘ and Δχ are the airflow temperature increment (°C) and moisture content increment (kg/kg) in section j, Q (m 3 /s) is airflow, θ ij (°C) is the surrounding rock temperature of node (i, j), ΔZ and ΔR is the step length along the axial and radial direction of the airway, m. Based on the airway and surrounding rock model established, the unsteady 2D heat conduction equation of surrounding rock in a cylindrical coordinate system can be derived, as shown below: (1) where θ is the surrounding rock temperature (°C); t is ventilation time (s); α is the thermal diffusion coefficient of surrounding rock (m 2 /s), α = λ/(ρC p ), λ is the thermal conductivity of surrounding rock (W/m°C); ρ is the density of surrounding rock (kg/m 3 ); C p is the specific heat of surrounding rock (J/kg°C); z is the axial length of the roadway (m); r is the radial distance from the z-axis (m); γ is the circumference angle (rad).
According to the assumption that the surrounding rock is homogeneous and isotropic, the heat conduction in the circumferential direction cannot be considered in the actual calculation process; that is, ∂θ/∂γ = 0 and ∂ 2 θ/ ∂γ 2 = 0, the unsteady heat conduction equation of surrounding rock is: Equation (2) is a 2D heat conduction model equation, and the equation neglecting the ∂ 2 θ/∂z 2 term is a 1D heat conduction model. To simplify the calculation process in practical application, the heat flux of surrounding rock in the axial direction of the airway is ignored. However, for the narrow strip of the wet airway wall, the axial temperature gradient is relatively large. To ensure the applicability and accuracy of the model, the influence of axial heat conduction is considered in the modeling process.
Equation (2) simultaneously considers the heat flow of surrounding rock along the radial and axial directions of the airway. When the heat flow of surrounding rock along the airway axis is ignored (∂ 2 θ/∂z 2 = 0), the 1D heat conduction equation can be obtained:

| Airflow temperature and humidity equation
Two heat transfer phenomena are considered when the airway wall is wet. 17 The first is sensible heat, which is F I G U R E 1 Schematic figure of airway and strata model.
defined as q C (W/m 2 ); The second is the latent heat from the state transfer of water and steam on the wet wall, which is defined as q L (W/m 2 ). When the latent heat is converted into sensible heat, the surrounding rock temperature θ unsteady state changes. The high-temperature rock makes the airway wall temperature higher through heat conduction. When the airflow is through the airway wall, the airflow temperature is lower than the airway wall temperature. A certain temperature difference exists, transferring heat between the roadway wall and the airflow. Therefore, in the process of underground ventilation, when the airflow passes through the wet airway wall, the water vapor molecules in the boundary layer will enter the unsaturated airflow, causing the water vapor molecular concentration in the boundary layer to decrease and return to the unsaturated state. At this time, the water evaporates, making the water vapor in the boundary layer reach the saturation state again. This is the convective mass transfer process between the airway wall and the airflow. It is assumed that the surrounding rock temperature is equal to the original rock temperature when the airway excavation is completed, and the surrounding rock temperature does not change with time at the location far enough from the airway. At the same time, it is defined that the temperature and humidity of airflow are equal on the section perpendicular to the airway axis, and the effect of surrounding rock thermal radiation on the airflow is not considered. According to Newton's cooling law, when there is a temperature difference between the surface of an object and its surroundings, the heat lost per unit area per unit time is proportional to the temperature difference, especially in the process of forced convection. In this case, the sensible heat in airflow is calculated as follows: where α is convective heat transfer coefficient (J/m 2°C ); θ W is the tunnel wall temperature (°C); Θ is the average temperature of airflow (°C). Airflow temperature along z-direction Θ calculated as follows: where U is the perimeter of the roadway (m); ρ a is the air density of airflow (kg/m 3 ); C′ p is the constant pressure specific heat of airflow, C′ p = 1004 J/kg°C; Q is air volume (m 3 /s). The partial pressure of water vapor on the completely wet wall is saturated, but the actual airway wall is partially wet. The mass concentration of water vapor in airflow shall be determined according to the partial pressure of steam. When the water vapor concentration on the wet wall is higher than that of the airflow, the water evaporation on the wall continues. The water evaporation rate w (kg/sm 2 ) is calculated from the following equation: where ϕ is the humidity, representing the ratio of mass transfer coefficient of a partial wet wall to that of a complete wet wall 6 ; β is the mass transfer coefficient, according to Lewis law, β = α/(ρ a C′ p ) (m/s); m W is the wall water vapor concentration; m is the concentration of water vapor in the airflow.
If m W > m, it will promote the continuous evaporation of wall water. Set up χ (kg/kg) is the moisture content, χ S is the saturated moisture content (kg/kg), then m and m W can be expressed as: Therefore, the humidity increase equation along the airflow direction can be expressed as: where p′ a is dry air density, p′ a = 1.2 kg/m 3 . At the same time, relative humidity Ѱ can be expressed as: When χ > χ S , the moisture in the airflow reaches the supersaturated state. The numerical algorithm based on heat and mass balance correction χ and Θ proposed by Sasaki et al. 18 is used to estimate the excess water vapor condensed into the fog in the wind flow.

| Boundary and initial condition setting
In the actual production of coal mine, due to the influence of water sprinkling, gushing, and seepage, the airway wall is often in a damp state. Under the humid condition of the airway wall, in addition to considering the convection heat transfer between the airway wall and the airflow, the convection mass transfer between the airway wall and the airflow should also be considered, that is, the transition of the state of water and steam from the wet wall. During the transition of water and steam, the moisture on the wet wall needs to absorb heat from the outside, that is, the latent heat required for water evaporation, expressed in q L (W/m 2 ). This paper assumes that this part of heat (latent heat) is provided by the airway wall.
Assuming that for a circular airway with radius R, the total heat transferred from the wet wall to the airflow q W is equal to the heat generated by the temperature gradient of the surrounding rock perpendicular to the airway wall at the boundary (r = R), then: where h l is the latent heat of water evaporation, h l = 2.50 × 10 6 -2370θ W (J/kg). Equation (11) is the boundary condition at the wall to determine the wall temperature in the numerical calculation θ W . Another temperature boundary condition in the calculation of surrounding rock heat conduction is that the surrounding rock temperature at the outer boundary far away from the radial direction is assumed to be constant (θ = θ 0 ). θ 0 is the original rock temperature, that is, the initial temperature of the surrounding rock at the beginning of ventilation (t = 0). The airflow temperature and relative humidity at the entrance of the airway (z = 0) are, respectively, expressed as Θ IN (°C) and Ѱ IN . Other conditions are set as follows: θ 0 = 30°C, a = 1.10 × 10 −6 m 2 /s, V (V = Q/(πR 2 )) = 4.0 m/s, α = 13.0 W/m 2°C , Bi is the Biot number (Bi = α R/λ) = 5.2 R.

| Airway wall wet surface models
In this paper, two models of the uniform wet surface model (UWSM) and partly wet surface model (PWSM) are discussed, as shown in Figure 2. 1. UWSM Set model wet cells to be sufficiently tiny so that the dry and wet parts cannot be distinguished, and it is regarded as a uniform humidity state, as shown in Figure 2A. ϕ indicates the degree of wetting of the mode, ϕ = 1 means a completely wet wall, and ϕ = 0 refers to the drywall surface. Assumption the parameter ϕ is constant in the roadway. The model can be applied to the ventilation simulation of a large ventilation network.

PWSM
The wall moist state of the model is composed of regular axial spaced moist strips, and the axial airway wall and airflow have nonuniform heat exchange. Define wet part ϕ = 1, dry part ϕ = 0. Suppose ΔL W is the width of the wet part, L S is the span between adjacent wet parts, as shown in Figure 2B, then the humidity rate ϕ a can be expressed as: The PWSM established in this paper is totally different from the model described in Starfield 19 and Uchino. 15 The equations proposed by Starfield and Gibson 20 used the 2D heat conduction model for the prediction algorithm of θ W to avoid complex heat conduction calculations. On the basis of this research, this paper establishes a 3D heat conduction model of the effect of axial heat conduction, proposes an axial strip moisture calculation model, and applies it to the calculation and analysis of ventilation networks.

| Numeral calculations
For the difference calculation of wet airway, the dimensionless analysis method is used to deal with the differential equations of surrounding rock temperature and airflow temperature and humidity. First, the where θ*, θ * W , Θ*, Θ* IN , r*, z*, and τ are the dimensionless quantities of variables θ, θ W , Θ, Θ IN , r, z, and t mentioned above. τ is Fourier number.
Substitute the relevant unitless variables in Equation (13) into the unsteady heat conduction Equation (1) of surrounding rock, and replace the variables in the formula to obtain the dimensionless unsteady heat conduction differential equation of surrounding rock: When τ = 0, set the airflow temperature of each point in the roadway at this moment equal to the inlet airflow temperature, and the temperature of each point in the surrounding rock is the original rock temperature.
When τ > 0, the boundary condition of the tunnel entrance (z* = 0) can be expressed as: Θ*= Θ* = 0 . At the outer boundary of the surrounding rock (r* = ∞), θ θ *= * = 1 0 . At the roadway wall (r* = 1): where Bi is the Biot number, Bi αR λ = / . The dimensionless calculation formula of dimensionless latent heat generated by water evaporation of the airway wall can be obtained, as shown in Equation (16): According to Taylor expansion and Figure 1, the dimensionless temperature at node i + 1 and i − 1 in the surrounding rock at time l is expressed by the dimensionless temperature at node i, which can be expressed as: For the studied surrounding rock area, assuming that Δr* is very small, its third power and higher power are infinitesimal, so the terms of the third power and higher power can be ignored. Subtract Equations (17) and (18) to get: The first-order partial derivative solution of the dimensionless temperature of surrounding rock along the radial direction can be obtained from Equation (19): Add Equations (17) and (18) to get: According to Equation (21), the solution of the second-order partial derivative of the dimensionless temperature of surrounding rock along the radial direction of the roadway is as follows: Similarly, the solutions of the first-order partial derivative of the dimensionless temperature of the surrounding rock with respect to time and the secondorder partial derivative of the dimensionless temperature of the surrounding rock along the axis of the roadway can be obtained as follows: In the formula, Δτ is the time step, Δr* is the space step of dividing the surrounding rock area r* along the radial direction of the roadway, and Δz* is the space step of dividing the distance z* along the axis direction of the roadway.
The dimensionless surrounding rock temperature analytic formula of nodes (i, j) at time l + 1 can be obtained, as shown in Equation (26): Therefore, when the original rock temperature, ventilation time, and other parameters are known, the dimensionless temperature of surrounding rock at any time and any location can be calculated by using this analytic formula.
According to Taylor expansion, the dimensionless temperature of surrounding rock at node i = 2 and i = 1 is expressed by the dimensionless temperature of surrounding rock at node i = 0, which can be expressed as: The finite difference method is used to calculate the surrounding rock temperature and airflow temperature and humidity. The model is shown schematically in Figure 1. In the ventilation network analysis, L = 100 m is a reasonable airway length interval. 5,18 The calculation area of the airway surrounding rock selected in this paper is 0 ≤ z ≤ 100 m (L = 100 m) and 1 ≤ r/R ≤ 46 (far boundary). The calculation area is divided axially into 500 (or 1000) equal intervals (Δz = 0.2 or 0.1 m), divided radially into 80 intervals (minimum Δr = 0.1R), realizing each boundary (wall) element ϕ all within 0 ≤ ϕ ≤ 1 is given separately. The calculation model consists of two heat conduction models and two wet wall models:  Table 1. It can be seen that the calculation results of Model I are basically consistent with those of Danko and CLIMSIM models (developed by McPherson); Model II is also very close to the TUNNEL calculation results developed by Hemp. It shows that the numerical model and calculation program used in this study are reasonable. In addition, the models listed in the table only consider the influence of radial heat conduction of surrounding rock.

| Numerical calculation results analysis
On the basis of comparison and verification with other commonly used 1D numerical models, the 2D numerical model and the 1D numerical model are used to simulate the partly wet airway, and analyze the difference between the results of Model II and Model IV on the calculation of surrounding rock temperature. And the influence of parameter ϕ, ϕ a , ΔL W , L S , Θ IN , Ѱ IN , R on θ, Θ, χ is discussed. Then the relationship between ϕ and ϕ a is studied.
When 2R = 4 m, θ 0 = 30°C, Θ IN = 10°C, Ѱ IN = 0.6, ΔL W = 5.0 m, L S = 50 m, t = 62, the temperature distribution of surrounding rock corresponding to Models II and IV is shown in Figure 3. Compared with the 1D heat conduction model, the calculated results under the two conditions are significantly different. In the 2D heat conduction model, the influence of axial heat conduction is considered, and the influence of the axial heat conduction of surrounding rock along the z-axis direction on the surrounding rock temperature decreases with the increase of its distance from the wall.  contour curvature of PWSM is affected by its axial width ΔL W obviously. When 2R = 4 m, θ 0 = 30°C, Θ IN = 10°C, Ѱ IN = 0.6, ϕ a = 0.2, t = 2 years, the axial distribution of tunnel wall temperature (θ W ) calculated by Models II, III, and IV are shown in Figure 5. θ W value of PWSM fluctuates under the influence of strip wet part. However, θ W value of UWSM axial variation is very small, and the value is between the maximum and the minimum value calculated by the PWSM.  Obtained from Model I, θ W numerical calculation program can be directly replaced by the dimensionless temperature gradient equation proposed by starfield or Gibson. 20 Therefore, in the complex humid environment described in Model IV, it is only necessary to introduce equivalent ϕ, the temperature and humidity of airflow in a large-scale mine ventilation network can be predicted.
In this study, the correction coefficients C W and C′ W of temperature and humidity are introduced to the correct the equivalent ϕ in ΔΘ and Δχ of the Model I. After correction, equivalent ϕ can be expressed as, In different humid environments, the calculation results of C W and C′ W are shown in Figure 10. When ΔL W ≥ 2R, ϕ a ≤ 0.3, the values of C W and C′ W are about 0.5, in this calculation, C W and C′ W are independent with Θ IN , Ѱ IN , R, and t. Indicates that ΔΘ and Δχ of PWSM can be calculated through ϕ = 0.5ϕ a in UWSM under the conditions of ΔL W ≥ 2R.
ϕ is a theoretical value and cannot be assigned by visual roadway parameters; ϕ a is easier to measure in the actual airway. Using Equation (13), the equivalent value of ϕ can be given according to the actual humidity of the airway under certain conditions. Therefore, a simple uniform moisture model is used to simulate the temperature and humidity of airflow.

| Comparison and analysis between actual measurement and simulation
In this paper, the working face, track transportation airway, strap-laneway and drift of the Xinzhi coal mine in Shanxi Province of China are selected for field measurement and simulation comparison. Parameter   information of surrounding rock, airflow and airway, [22][23][24] and so forth, as shown in Table 2. Figure 11 shows that the simulation results are basically consistent with the actual measurement results. In the whole airway, the measured values of airflow temperature and relative humidity fluctuate with the change of airway length. It is considered that this is due to the complexity of the actual airway, such as the impact of various environmental conditions in the measurement process.
From Figures 12-14, it can be found that the airflow temperature first decreases and then gradually increases with the increase of the airway length, while the relative humidity first increases rapidly and then tends to be stable with the increase of the airway length. Affected by the characteristics of the dry airflow, it will quickly absorb moisture in the wet airway first. When the relative humidity reaches about 99%, the moisture content in the airflow tends to be saturated, and its increase with the change of airway length is not obvious.
From Figures 11-15, it can be seen that the calculated results of inlet and outlet airflow of working face, left-wing track transportation airway, left-wing straplaneway and 540 adit primarily in keeping with the measured results. In addition to surrounding rock heat dissipation and air self-compression, the heat sources in the actual airway also include heat generated by hot water and the operation of electromechanical equipment, and so forth. The paper only considers the main reasons for the high-temperature environment in the mine, ignoring the impact of other secondary heat sources on the temperature and humidity of the airflow in the airway. In addition, the underground conditions are complex, and there are coupling effects of many factors, which affect the change of airflow temperature and humidity, leading to the deviation of some measuring points.  The comparison between the measured and simulated data in the Xinzhi coal mine shows that the prediction model can meet the needs of mine climate prediction in the actual production process and proves the applicability and accuracy of the prediction model of airflow temperature and humidity in wet airway proposed in this paper.

| CONCLUSIONS
Based on the theory of 2D surrounding rock temperature field, this paper establishes a numerical calculation model of surrounding rock temperature, together with radial heat conduction and axial heat conduction. Supported by this numerical calculation model, a numerical calculation model of tunnel air temperature and humidity under the 2D heat conduction model is established to calculate and simulate the surrounding rock temperature distribution under the humid condition of the wall. This paper studies the influence of the wet air ratio of the axial strip humid wall on the airflow temperature and moisture content, quantitatively analyzes the difference of the influence of the humidity and wet air ratio on the airflow climate, and discusses the feasibility of using the equivalent humidity of the uniform humid wall model to replace the correction coefficient of the wet air ratio of some humid wall models. The conclusions are as follows: 1. The 2D heat conduction numerical model is used to calculate the temperature distribution of the airway surrounding rock. The temperature distribution of partly wet surface airway surrounding rock (especially airway wall) is obviously different from that of the uniform wet surface airway model. It is more reasonable to establish the airflow temperature and humidity prediction formula based on the 3D heat conduction model. 2. The temperature increment decreases with the increase of wet air ratio, inlet air temperature and airway radius and increases with the increase of inlet relative humidity; The increment of moisture content increases with the increase of wet air ratio. 3. The equivalent humidity of the UWSM is proposed to replace the correction coefficient of the wet air ratio of the PWSM. When the wet part width of the wall is larger than the airway diameter, the humidity is defined as 50% of the wet air ratio, and the objective airflow temperature and humidity change is obtained using the UWSM. 4. Through the actual measurement of airway data in the Xinzhi coal mine, the difference between the measured and simulated data is compared and analyzed, which verifies the applicability and accuracy of the prediction model proposed in this paper.