Application of Data Envelopment Analysis in the Ventilation and Cooling Efficiency Evaluation of Hot Development Headings

Abstract

The thermal environment of the development headings in deep underground mines tends to be harsh. The design of a reliable auxiliary ventilation system is an important strategy to maintain the thermal environment and air quality in the workplace. Hence, it is of great significance to comprehensively evaluate the effects of different air supply conditions and explore the optimal solution for the auxiliary ventilation system. In this study, cooling efficiency (CE) and ventilation efficiency (VE), which are the two most widely used indices for evaluating the ventilation in public places, were introduced into the field of underground mine ventilation. Numerical simulations of multiple air supply conditions were carried out to obtain the CE and VE, respectively. A data envelopment analysis model was used innovatively to evaluate and rank the efficiency of the different conditions. In the evaluation model, the air supply temperature and air supply velocity were the inputs, whereas CE and VE were the outputs. The preliminary results showed that cases with either the highest VE or highest CE alone may not actually indicate efficient air supply conditions. The optimum efficiency was determined to be at an air supply temperature difference of 6 °C. Thereafter, the air supply rate could be determined using the psychrometric chart based on the cooling and moisture loads. These results can provide theoretical guidance for optimizing air supply conditions and energy-saving control of mining ventilation and cooling systems.

1. Introduction

The high temperature of roadways caused by geothermal energy has become the main factor restricting the safe production of underground mines. The intensive heat transfer between the hot rock surface and airflow, coupled with the heat and exhaust gas produced by diesel equipment during excavation, can create an intolerable hot and humid environment that negatively impacts the performance, overall productivity, and most importantly the safety of the underground workforce. The deterioration in environmental conditions may also lead to thermal discomfort and heat-related illnesses such as thermal stress, heat cramps, heat rash, and heat stroke [1,2]. Hence, mining ventilation and cooling equipment must be used to ensure that the thermal environment of the workplace is within the acceptable limits for miners.

The development headings in the hot underground mines are typically hot workplaces because they are at the end of the ventilation system. A forcing auxiliary ventilation system combined with a heat exchanger is widely used to ventilate and remove heat in these places. Hence, the thermal environment is affected by the air supply conditions of the air duct (i.e., air supply velocity v, air supply temperature t_s, and distance between the outlet of the air duct and the heading face Z_m) and heat produced by the surrounding rock, mining machinery, and mining workers, as shown in Figure 1. Moreover, the air supply conditions of the air duct also determine the energy usage of the ventilation and cooling (VC) systems of mines. Related studies have shown that VC systems are responsible for up to a quarter of the energy consumption of a typical deep-level mine. Hence, reasonable control of the air supply conditions has significant cost-saving potential. Different combinations of v, t_s, and Z_m may provide the same thermal comfort to a mining area but differ in energy consumption. Hence, the optimal combination of v, t_s, and Z_m needs to be determined so that the VC system can run under economic conditions.

To explore the optimal combination of different air supply conditions, computational fluid dynamics (CFD), which is characterized by low costs and easy visualization, has come to play an important role [3]. Scholars have designed case studies and researched the variation laws of methane/dust dispersion, the diesel particulate matter concentration, and the thermal environment by numerically modeling differences in the air supply conditions. Kurnia et al. studied the propagation patterns of methane dispersion with different v [4]. Hasheminasab et al. considered auxiliary ventilation consisting of a brattice and ducted fan and employed different v and Z_m to study variations in the velocity field and methane dispersion [5]. Torano et al. numerically studied the relationship between the air supply conditions and methane/dust behavior in mining areas [6]. Yu et al. and Liu et al. simulated the variation in dust dispersion with different v and Z_m [7,8]. Kurnia et al. proposed various Z_m and v conditions to handle hazardous gases from diesel emissions [9]. Zheng and Tien concluded that increasing v can help reduce the diesel particulate concentrations at the mining face [10]. Sasmito et al. examined different air supply conditions to develop an optimum VC strategy [11]. Wang et al. investigated the influence of different Z_m on the distribution of a roadway’s air temperature and humidity [12]. To sum up, these studies have put forward constructive guidance on the air supply conditions of VC systems. However, it should be noted that, when optimizing these systems, it is necessary to consider the optimal solutions from the perspective of efficiency.

In the field of building and the environment, to evaluate the performance of VC systems, techno-economic analysis has often been used to carry out comparative studies on the impact of different cases on the cost. For example, Mao et al. employed the technique for order preference by similarity to an ideal solution to evaluate the overall performance of the air conditioning system under different air supply conditions and obtained the optimal solution [13]. Esen et al. used the cooling performance coefficients to compare the energy efficiency of a ground-coupled heat pump system and air-coupled heat pump system. The techno-economic analysis clearly shows that ground-coupled heat pump systems are economically preferable to air-coupled heat pump systems for the purpose of space cooling [14]. Habibi et al. used the seasonal coefficient of performance to compare the cooling performance of ground-coupled ejector cooling systems and air-coupled ejector systems. The results revealed that the relative payback period of the two systems can be as low as five years [15]. Ozyogurtcu et al. used the life-cycle cost to evaluate the techno-economic performances of a ventilation system assisted with exhaust air heat recovery, electric heater and solar energy. The results showed that, compared with the conventional ventilation system, using the ventilation system assisted with an electric heater can reduce energy consumption by 86%. The payback period of the system was 5 years and 8 months [16]. Aste et al. used the global cost calculation method to conduct the techno-economic performance of the application of different window features in office buildings under three different representative European climates. The results showed that the use of solar shading devices is the most convenient option in all analyzed cases [17]. In summary, the optimal scheme can be obtained by conducting a techno-economic analysis as long as the relevant data are available.

However, according to the aforementioned research, in the process of the studies on techno-economic analysis, the cost of each case was composed of various types, and appropriate simplification and assumptions must be made in the calculation process, resulting in the uncertainty of evaluation results. Based on the characteristics of the underground VC system, it can be considered that the efficiency of VC systems can be calculated by multiple input–output models composed of multiple indices. The output indices (i.e., VC effects) are composed of ventilation efficiency (VE) and cooling efficiency (CE), while the input indices (i.e., energy consumption) can be represented by v and t_s. Therefore, the optimal air supply conditions are cases that can achieve the optimal input–output balance. However, to the best of our knowledge, the existing methods give weights to these indices to evaluate efficiency, which makes it difficult to ensure objectivity. The data envelopment analysis (DEA) model can objectively weight the input and output indices [18], which can avoid the shortcomings of the aforementioned methods. Hence, it is believed that the DEA model can be used to evaluate the air supply conditions of the VC systems in hot mines.

Therefore, to fill the gaps in the aforementioned studies, the objective of this study is to apply the DEA model to the performance evaluation of ventilation engineering in underground development headings and verify the practicability and accuracy of the evaluation results, so as to realize the intersection of management and engineering. First, the three-dimensional numerical CFD model was established based on the actual situation of a development heading in an underground gold mine. To construct the DEA evaluation model, the VE and CE at different air supply conditions can be obtained by the simulated air velocities, temperatures, and humidity. Finally, the DEA model was used to rank the performances, and the optimal air supply conditions can be proposed based on the evaluation. The findings of this study are expected to provide theoretical support for effectively reducing the energy consumption of VC systems in hot mines.

2. Methodology

2.1. CFD Numerical Simulation

2.1.1. Governing Equations

A commercial CFD code (i.e., Airpak 3.0) was used in this study. To predict the velocity, temperature and moisture fields of the airflow and the air quality in development headings, the governing equations were solved by Airpak 3.0 in a steady state. In the case of underground mining ventilation, the airflow should be considered as the incompressible fluid [19] (i.e., the density of the airflow is constant, which was 1.29 kg/m³ in this study) due to the relatively low velocity. Hence, the governing equations of mass conservation can be written as follows [20]:

$$\nabla ·\overrightarrow{v}=0$$

(1)

where $\overrightarrow{v}$ is the velocity vector (m/s).

Transport of momentum in an inertial reference frame is described by [20]:

$$\frac{\partial }{\partial t}\left(\rho \overrightarrow{v}\right)+\nabla ·\left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-\nabla p+\nabla ·\left(\stackrel{=}{\tau }\right)+\rho \overrightarrow{g}+\overrightarrow{F}$$

(2)

where p is the static pressure (Pa), $\stackrel{·}{\tau }$is the stress tensor (Pa), $\rho \overrightarrow{g}$ is the gravitational body force (Pa), and $\overrightarrow{F}$ includes the pressure, contact force and surface tension in the domain (Pa).

The energy equation for a fluid region can be written as [20]:

$$\frac{\partial }{\partial t}\left(\rho h\right)+\nabla ·\left(\rho h\overrightarrow{v}\right)=\nabla ·\left[\left(k+k_t\right)\nabla T\right]+S_h$$

(3)

where h is sensible enthalpy (J/g), k is the molecular conductivity and k_t is the conductivity caused by turbulent transport, T is the temperature (°C), and S_h includes any heat source defined (J).

The species transport can be expressed as [20]:

$$\frac{\partial }{\partial t}\left(\rho Y_i\right)+\nabla ·\left(\rho \overrightarrow{v}{Y}_i\right)=-\nabla ·{\overrightarrow{j}}_i+S_i$$

(4)

where Y_i is the mass fraction of component i (J/kg), ${\overrightarrow{j}}_i$ is the mass diffusion flux of component i (kg/(m²∙s)), and S_i is the generation rate of other customized mass source terms (kg/(m²∙s)).

In turbulent flows, Airpak computes the mass diffusion in the following form:

$${\overrightarrow{j}}_i=-\left(\rho D_{i,m}+\frac{{\mu }_t}{S{c}_t}\right)\nabla Y_i$$

(5)

where D_i,m is the mass diffusion coefficient for species i in the mixture (m²/s). $S{c}_t$ is the turbulent Schmidt number (𝜇_𝑡/𝜌𝐷_𝑡 where D_t is the turbulent diffusivity). The default $S{c}_t$ is 0.7.

According to the airflow state in the heading, it can be considered that the airflow was fully turbulent. To obtain accurate simulation results, the standard k-ε turbulence model proposed by Launder and Spalding [21] combined with the enhanced wall treatment model were used. The enhanced wall treatment is a near-wall modeling method that combines a two-layer model with reinforcement wall function. The standard k-ε turbulence model is described as:

$$\rho \frac{\partial }{\partial t}k+\rho \frac{\partial }{\partial x_i}\left(ku\right)=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+G_k+G_b-\rho ϵ$$

(6)

$$\rho \frac{\partial }{\partial t}ϵ+\rho \frac{\partial }{\partial x_i}\left(ϵu\right)=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{ϵ}}\right)\frac{\partial ϵ}{\partial x_i}\right]+C_{1ϵ}\frac{ϵ}{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2ϵ}\rho \frac{{ϵ}^2}{k}$$

(7)

where k is turbulent kinetic energy (J), μ is dynamic viscosity of the airflow (Pa·s), μ_t is turbulent viscosity (Pa·s), G_k represents the generation of turbulence kinetic energy due to the mean velocity gradients, and G_b is the turbulent kinetic energy due to buoyancy, ε is the dissipation rate, C_1ε, C_2ε and C_3ε are turbulence model constants, and σ_k and σ_ε are the turbulent Prandtl numbers for k and ε, respectively.

In Airpak, turbulent heat transport is modeled using the following equation:

$$\frac{\partial }{\partial t}\left(\rho E\right)+\frac{\partial }{\partial x_i}\left[u_i\left(\rho E+p\right)\right]=\frac{\partial }{\partial x_i}\left(k_{\mathrm{eff}}\frac{\partial T}{\partial x_i}\right)+S_h$$

(8)

where E is the total energy (J) and k_eff is the effective conductivity.

2.1.2. Geometry Model

The geometry models of the development heading were established with Rhino 6.0 (a CAD tool) based on an actual route heading in Xiadian Gold Mine in China. The model described the scenario in which one piece of load haul dump (LHD) machinery was carrying out the mucking activities in the heading as shown in Figure 2. The section shape of the development heading can be approximately regarded as a rectangle, which had a length of 40 m, a width of 4.6 m, and a height of 4.9 m. The heading face was ventilated by a forcing air duct, which had a section diameter of 0.8 m. The center of the air duct’s section was 4 m away from the floor. The outlet of the air duct was 20 m away from the heading face initially. The LHD had a height of 2.0 m, width of 2.5 m, and length of 8.0 m. During the mucking operation, the front end of the LHD was about 2 m away from the heading face. In the following section, the positive X direction is designated as the pointing from the entrance to the heading face, the positive Y direction is from the left side to the right when facing the entrance, and the positive Z direction is from the floor to the roof.

2.1.3. Mesh Generation and Mesh Independence Tests

The number of elements in a mesh has a great influence on the accuracy of the numerical results and the time-consumption of the numerical simulation. Hence, mesh independence tests are necessary to achieve the balance. First, a coarse mesh with 0.81 × 10⁶ elements (Coarse Mesh) was generated by using the default Hexa unstructured mesher provided by Airpak as shown in Figure 3a. Subsequently, another two finer meshes with 1.13 × 10⁶ elements (Fine Mesh) and 1.56 × 10⁶ (Finer Mesh) elements were generated, respectively, by using the meshing controls as shown in Figure 3a. Airpak also combines the two-layer model with enhanced wall functions. If the near-wall mesh is fine enough to be able to resolve the laminar sublayer (typically y+ ≈ 1), the viscosity-affected near-wall region is completely resolved all the way to the viscous sublayer by the two-layer model. The refined mesh can produce a finer mesh in the boundary layer of the strata and the LHD as shown in Figure 3a. The locations of plumb line (Y = 1 m) and horizontal line (Z = 2.5 m) at the plane of X = 10 m (Figure 2) were selected to show the simulation results as shown in Figure 3b. Figure 3c shows the comparisons of the simulation results of the air velocities with the three different mesh elements. It can be found that the simulation results using the mesh with 1.13 × 10⁶ elements were closer to that of 1.56 × 10⁶ elements. Due to the large deviations from the simulation results obtained with the other two meshes, the accuracy of the simulation using the mesh with 0.81 × 10⁶ elements can be considered poor and the simulated results obtained by the mesh with 1.13 × 10⁶ elements can achieve acceptable balance between accuracy and time-consumption. This study used this meshing scheme. The grid convergence index (GCI) method [22] was also applied to evaluate the impact of this mesh on the results. Using second-order convergence and setting the safety factor and mesh refinement factor to 2.5 and 1.5, respectively, the maximum value of GCI was 2.6%, which shows that the influence of the elements on the results can be negligible for the fine mesh.

2.1.4. Boundary Conditions and Solution Schemes

To collect reliable boundary conditions for the numerical simulations, the data collection and the in-situ measurements were carried out in the development heading. The size of the heading and the production situation conducted in the heading were consistent with the geometry models described in Section 2.1.2. The rock type and the thermophysical parameters are listed in

Table 1. Thermophysical parameters of granite.

Rock Type	Thermal Conductivity (W/(m·K))	Specific Heat Capacity (J/(kg·K))	Thermal Diffusivity ×10−6 (m2/s)	Density (kg/m3)
Granite	2.8	790	1.41	2693

. The virgin rock temperature reached about 45 °C. The dynamic pressure in the forcing duct was maintained at 160 Pa throughout the survey period, which equated to a volumetric flow rate of 6 m³/s and a mass flow rate (M_s) of 7.74 kg/s (the density of the airflow was 1.29 kg/m³). The information was provided by the surveying department of the mining company.

A Testo 440 anemometer was employed to measure the air velocity, and KE-COS-03 data-logging devices were used to record the air temperature and humidity data at 30 s intervals, respectively. The general characteristics of the two devices are shown in

Table 2. General characteristics of the measurement devices.

Equipment	Airway VariableMeasured	Range	Accuracy
KE-COS-03 data-logging device	Air temperature (°C)	−40–80 °C	0.2 °C
	Relative humidity (%)	0–100%	1.5%
Testo 440anemometer	Air velocity (m/s)	0–50 m/s	0.03 m/s + 4% measuredvalue

. As shown in Figure 4, the transducer c1 was installed at the outlet of the air duct to measure t_s. The other six transducers need to be installed across a tunnel cross-section at X = 10 m, where the airflow directions and velocities were relatively stable. Consequently, as shown in Figure 4, these transducers (groups A (a1–a3) and group B (b1–b3)) were mounted 30 m back from the heading face at the plumb lines A and B, respectively. Both lines were approximately 1 m away from the roadway surfaces, as shown in Figure 4.

The measured results of the air temperature variation with respect to time during a mining cycle is illustrated in Figure 5. According to the mining record, the cycle can be divided into the non-operational period and the mucking period. In the non-operational period, since there was no LHD in the heading, the air temperatures were significantly lower than those in the mucking period. This trend can also be clearly observed from Figure 5. In Figure 5, “airway” represents the average measured temperature (t) of the six transducers (i.e., groups A and B), and “outlet of the air duct (t_s)” represents the measurements of transducer c1.

Table 3. Measured results of moisture contents at different periods.

Period of Measuring Time	Measurement Points
	Relative Humidity/%	Moisture Content (g/kg Dry Air)
Before 8:00 a.m.	52	13.95
8:40–9:45 a.m.	53	13.99
9:45–11:15 a.m.	41	14.67
11:15–11:45 a.m.	44	14.62
11:45 a.m.–1:00 p.m.	45	14.35
After 1:00 p.m.	50	14.13

provides the measured moisture contents of the airflow at different periods. It can be observed that the moisture content of the air varied only slightly during the field measurements, which indicates that there was no moisture source in the development heading.

In the development heading, cooled air is passed from the air duct into the heading face. When there is a temperature difference between the continuously exposed rock surface/LHD and airstream, heat transfer will take place. Hence, the main sources of thermal load inside the development heading were the strata heat (Q_s), the heat produced by an LHD (Q_l). Considering the energy conservation law, the thermal balance in the development heading can be expressed as follows:

$$Q_{\mathrm{air}}=Q_s+Q_l$$

(9)

The thermal load of airflow can be determined by:

$$Q_{\mathrm{air}}=M_s\times \left(H_s-H_o\right)$$

(10)

where Q_air is heat taken away by the airflow (kW); Ms is mass flow of the airflow (kg/s), which is 7.74 kg/s; H_s is the enthalpy of the airflow at the outlet of the air duct (kJ/kg), which is 60.78 kJ/kg; and H_o is the enthalpy of the airflow at the measuring points (kJ/kg).

In this study, the airflow density was set as a constant of 1.29 kg/m³. The averaged air temperature and the relative humidity at the outlet of the air duct were about 25 °C and 70%, respectively, based on the field measurements, Hence the H_s was 60.78 kJ/kg. According to Figure 5, before 9:45 a.m., no working activities took place; hence, there were no Q_l. The ventilation system removed only the heat from the strata. The air temperature difference between t_s and t was stable at about 5.4 °C and the difference in the moisture contents was a constant. From Equations (9) and (10), the calculated Q_s was 42.3 kW and the heat transfer area of the strata was 592.5 m². Hence, the heat flux from the strata was 71.4 W/m². After 9:45 a.m., an LHD entered the mining face. At its peak points, LHD performed its heaviest work. The Q_l and Q_s increased the air temperature to more than 36 °C at the peak points. However, the average air temperature difference between t_s and t at the mining period was about 9.4 °C. Hence, the average Q_l can be calculated as 31.3 KW.

The turbulent kinetic energy and dissipation at the outlets can be determined by the turbulence intensity (I) and turbulence length scale (l), respectively, from Equations (12) and (13):

$$Re=\rho \cdot u\cdot \frac{L}{\mu }$$

(11)

$$I=0.16R{e}^{-\frac{1}{8}}$$

(12)

$$𝓵=0.07L$$

(13)

where Re is the flow Reynolds number, L is the hydraulic diameter of the duct or airway (m), and μ is dynamic viscosity (N·s/m²). The factor of 0.07 is based on the maximum value of the mixing length in fully-developed turbulent pipe flow. I and l were 3% and 0.056 m, respectively, at the air duct outlet and 4% and 0.329 m, respectively, at the airway outlet.

For the treatment of the roughness of the rock surface, the average height of the surface texture was set to 0.01 m based on Peltier et al. [23]. Due to the application of the thermal insulating material of rubber plastic sponge on the air duct, the wall of the air duct was adiabatic. The boundary conditions are summarized in

Table 4. Boundary conditions of CFD method.

Boundary	Conditions
Air duct outlet	ts = 25 °C, v = 12 m/s, relative humidity = 70%, I = 3%, l = 0.056 m
Airway outlet	Pressure outlet, I = 4%, l = 0.329 m
Wall of strata	Thermal conductivity 2.8 W/(m·°C), Heat flux 71.4 W/m2, average height of the surface texture 0.01 m
Wall of LHD	Fixed heat flux 31.3 kW
Wall of the air duct	Adiabatic

. The SIMPLE algorithm was used with a second order scheme for the convective terms

2.1.5. Model Validation

The average air temperatures measured at each measuring point of group A and group B at period of 9:45–11:15 a.m. (i.e., the period of mucking operation) were compared with those of the CFD results. The layout of measuring points and measuring equipment are as described above. To avoid the influences of the different temperature scales on the comparison results, the dimensionless temperature of each point should be used. The dimensionless temperatures of each measuring point for the measured (t_dm) and the simulated results (t_ds) can be defined as follows, respectively:

$$t_{dm}=\frac{t_p-t_s}{t_s}$$

(14)

$$t_{ds}=\frac{t_n-t_s}{t_s}$$

(15)

where t_p is the average air temperature measured at each measuring point (°C), t_s is the air supply temperature (°C) and t_n is simulated air temperature at each measuring point (°C).

The deviation between the measured (t_dm) and the simulated results (t_ds) can be calculated as follows:

$$deviation=\frac{\left|t_{dm}-t_{ds}\right|}{t_{dm}}$$

(16)

Table 5. Comparisons between CFD results and the measured data.

Measuring Point	a1	a2	a3	b1	b2	b3
tdm	0.32	0.30	0.29	0.32	0.31	0.33
tds	0.28	0.33	0.33	0.29	0.30	0.31
deviation	12.5%	10.0%	13.7%	9.4%	3.2%	6.1%

shows the raw data of the comparisons. The largest deviation occurred at point a3 with a deviation of 13.7%. However, the average deviation between the measured data and CFD results is less than 15%. The deviation may be caused by the approximate treatment of the geometry of the heading and locations of the measuring points.

2.2. Numerical Studied Cases

The purpose of this study is to evaluate the VC performances brought by different air supply conditions. Therefore, the independent variables of the numerical simulations consist of v, t_s and Z_m. According to the common air supply conditions, the Z_m was set at five different distances, which varied from 5 m to 25 m with 5 m increments. The t_s were varied from 21 °C to 25 °C with 2 °C increments and the v varied from 8 m/s to 12 m/s at 2 m/s intervals based on the related regulations [24].

Table 6. Independent variables and the studied cases.

Air Supply Conditions		Case Number
v (m/s)	ts (°C)	Zm = 5 m	Zm = 10 m	Zm = 15 m	Zm = 20 m	Zm = 25 m
8	21	1.1	2.1	3.1	4.1	5.1
8	23	1.2	2.2	3.2	4.2	5.2
8	25	1.3	2.3	3.3	4.3	5.3
10	21	1.4	2.4	3.4	4.4	5.4
10	23	1.5	2.5	3.5	4.5	5.5
10	25	1.6	2.6	3.6	4.6	5.6
12	21	1.7	2.7	3.7	4.7	5.7
12	23	1.8	2.8	3.8	4.8	5.8
12	25	1.9	2.9	3.9	4.9	5.9

lists the details of all 45 simulated cases. The study aims to employ the DEA model to study the cost–performance relationship between supply air conditions and VC performances for each studied case based on the numerical simulations.

2.3. Evaluation Indices

Choosing the appropriate dependent variables to accurately characterize the VC performance of the auxiliary ventilation system is very important. The purposes of the VC system are to maintain the thermal environment and air quality in the workplace. Hence, the CE and VE were selected as the dependent variables to represent the cooling and ventilation performances, respectively. Note that appropriate modifications should be conducted due to the application scenarios of underground mines.

2.3.1. Index of Cooling Efficiency

In the case of the cooling performance, CE refers to the effective utilization rate of the cold air provided to the workplace to achieve the specified thermal environment. A higher CE indicates that the cooling capacity has been effectively utilized in the occupied zone (i.e., the workplace). The CE can be determined by [25]:

$$CE=\frac{t_e-t_s}{t_m-t_s}$$

(17)

where t_e is the temperature of the airflow at the extract (°C), t_s is the temperature of supply air (°C), and t_m is the mean value of the airflow temperature in the workplace (°C).

However, note that the CE only employs the parameters of dry bulb temperatures. The humidity of the airflow is also one of the important factors affecting the thermal environment. Roghanchi et al. [2] recommended using the “humidex” index to evaluate the thermal environment for hot mines. The advantage of this index is that both the dry bulb temperature and humidity of the airflow are considered. The humidex is calculated as follows:

$$humidex=t+\left[0.5555\times \left(e-10\right)\right]$$

(18)

where t is the dry bulb temperature of air flow (°C), and e is the partial pressure of water vapor of air flow (hPa).

Through numerical simulation, the t and relative humidity of the airflow in a specified zone can be obtained. The psychrometric chart or calculator can be used to determine the e. Thus, for the application scenarios of underground mines, humidex can be used instead of the dry bulb temperature in the original CE equation. The modified CE can be given by:

$$\mathrm{CE}=\frac{\left(humide{x}_{uz}-humide{x}_s\right)}{\left(humide{x}_{oz}-humide{x}_s\right)}$$

(19)

where humidex_uz represents the average humidex of unoccupied zone, humidex_oz represents the average humidex of occupied zone, and humidex_s represents the air supply humidex.

The modified CE of each case can be determined by the difference in the humidex for air flowing in and out of the occupied zone (i.e., the workplace). The higher the value of CE, the better the cooling performance of the VC system. Note that, for the cooling situation, the air supply humidex at the outlet of the air duct (humidex_s) must be lower than the ambient humidex (humidex_uz and/or humidex_oz).

2.3.2. Index of Ventilation Efficiency

Another objective of the ventilation system is maintaining air quality for the workplace. According to the research of Sandberg and Sjoberg [26], the local mean age of air (MAA) can predict the air quality well. The MAA represents the average lifetime of air at a measurement point in the airway relative to the time when it first entered the airway. A higher value of MAA implies that the cool air did not mix properly with the ambient air and that delivering air to an area is difficult. The MAA at a measurement point can be calculated as follows:

$${\tau }_p=\frac{{{\displaystyle \int }}_0^{\infty }{C}_p\left(t\right)dt}{C\left(0\right)}$$

(20)

where τ_p is the local MAA at point p, C_p(t) is the concentration at point p at time = t, and C(0) is the concentration at point p at time = 0.

For a ventilated heading, the lower the value of the MAA, the fresher the airflow in the space. With the assistance of Airpak 3.0, the MAA in a specified zone for each studied case can be obtained. Then, the VE can be calculated by:

$$\mathrm{VE}=\frac{{\tau }_n}{2{\tau }_p}\times 100\%$$

(21)

where τ_n is nominal time constant = airway volume (m³)/air flow rate (m³/s).

Hence, the higher the values of VE, the better the ventilation performance of the VC system. In this study, CE and VE were used to comprehensively evaluate the VC performance under different air supply conditions. The purpose of the VC system is to minimize the humidex of the occupied zone and avoid consuming the limited cooling capacity of the unoccupied zone. In addition, the air quality in the operation area needs to be maintained; hence, MAA needs to be minimized in the occupied zone. According to Equations (8) and (10), the energy consumption needs to be minimized to ensure high values of CE and VE, which are influenced by the air supply conditions, occupied zone, and unoccupied zone, as well as other parameters such as the occupancy pattern and the location and strength of the heat and airflow. In this study, the air temperature, humidity, and MAA of the above-mentioned 45 cases were numerically simulated using CFD, and CE and the VE values were calculated for each case. These values were input into the DEA model for a comprehensive evaluation.

2.4. Evaluation Method

2.4.1. Introduction to the DEA Model

In 1978, Charnes, Cooper, and Rhodes established the DEA model for estimating the efficiency of homogeneous organizational units called decision-making units (DMUs) that use the same inputs to produce the same outputs [18]. In this study, each case listed in Table 6 can be considered as a DMU. The inputs of each case were v and t_s, and the outputs were VE and CE. The efficiency of each case can be quantitatively evaluated according to the ratio of the outputs to the inputs as shown in Figure 6.

Because the evaluation process involved two input indices and two output indices, the efficiency could not be directly calculated. Thus, each index was assigned a certain weight, and the ratio of the weighted output to the weighted input was calculated to generate an index reflecting efficiency. The DEA model was used to objectively weight the input and output indices according to the data rather than subjectively, such as through expert consultation or discussion [18].

Depending on the method of measuring efficiency, the DEA model can be input-oriented, output-oriented, or non-oriented. This study focused on minimizing the inputs by decreasing v and increasing t_s without reducing VE and CE, which is in line with the application scope of the input-oriented method. Figure 7 shows the basic principle of the input-oriented method. Suppose that there are four DMUs with two inputs (x₁, x₂) and one output (y). The x-coordinate is the number of inputs x1 per unit of output (x₁/y). The y-coordinate is the number of inputs x₂ per unit of output (x₂/y). The four DMUs are, respectively, represented by points A–D on the coordinate system. A high efficiency means that a smaller input can produce a greater output. Perpendicular lines can be drawn to the two coordinate axes from the coordinates of the four DMUs. In the area enveloped by the two perpendicular lines and coordinate axes (including the boundary), the coordinate value of any point is less than or equal to that of the evaluated DMU. Among the four DMUs, C and D did not contain any other DMUs in the area enveloped by their perpendicular lines and coordinates, which means that they are at the frontier of technical efficiency. The curve and its extension line formed by the two-point connection are called an efficient frontier. The efficiency of a DMU on the frontier is 1, and the efficiency of a DMU enveloped by the frontier is 0–1. A′ is the projection of point A on the frontier. The efficiency of point A can be expressed as OA′/OA. The efficient part of point A is OA′, while the inefficient part is A′A.

Note that there will be a slack problem in the case of multiple inputs and multiple outputs. As shown in Figure 7, the projection of A on the frontier falls on A′ in the section parallel to the coordinate axis. Thus, DMU A can only achieve weak efficiency, even if the efficiency is improved to A′ on the efficient frontier through technical means. The investment of x₂ needs to be further reduced to improve the efficiency to point C. In the DEA model, the interval parallel to the coordinate axis (e.g., segment A′C) is called the input slack variable. This easily occurs in the basic DEA model with fewer DMUs. As shown in Figure 7, if few DMUs are evaluated, the generated frontier will not be sufficiently refined. Thus, the projection of the evaluated DMU can easily fall in an area parallel to the coordinate axes, which results in the slack problem. Hence, the DEA model needs to be improved to obtain more accurate evaluation results.

2.4.2. CCR and Super-Efficiency DEA Models

There are two types of DEA model: radial and non-radial. Radial models are represented by the Charnes-Cooper-Rhodes (CCR) model, which is the basic DEA model. Assume that a set of observed DMUs {DMU_j; j = 1, ..., n} is associated with m inputs {x_ij; i = 1, …, m} and s outputs {y_rj; r = 1, …, s}. The efficiency of the jth DMU is defined as follows:

$$Eff=\frac{{{\displaystyle \sum }}_r{u}_r{y}_{rj}}{{{\displaystyle \sum }}_i{v}_i{x}_{ij}}$$

(22)

where y_rj is the amount of the rth output from DMU_j; u_r is the weight given to the rth output; x_ij is the amount of the ith input used by DMU_j; v_i is the weight given to the ith input.

Basically, a radial model considers proportional changes in inputs or outputs. Hence, the efficiency of the kth DMU is defined as follows for an input-oriented linear programming model:

$$Eff=\max {{\displaystyle \sum }}_r{u}_r{y}_{rj0}$$

(23)

subject to:

$${{\displaystyle \sum }}_r{u}_r{y}_{rj}-{{\displaystyle \sum }}_r{v}_i{x}_{ij}\le 0$$

$${{\displaystyle \sum }}_i{v}_i{x}_{ij}=1$$

(24)

$$u_r,v_i\ge 0$$

If Eff < 1, the DMUs are inefficient; if Eff = 1, the DMUs are efficient. With the CCR model, multiple DMUs may sometimes be determined to be highly efficient (i.e., Eff = 1), such as C and D in Figure 7. If the number of DMUs is small, most DMUs will have an efficiency value of 1, so the efficiencies of different DMUs cannot be compared.

Another shortcoming of radial models is that they neglect slack when evaluating the efficiency. In many cases, the remaining non-radial slack can be observed. Hence, if such slack is important to evaluating the managerial efficiency, radial approaches may lead to erroneous decisions for evaluating the performance of DMUs. Therefore, many scholars have improved the CCR model and proposed super-efficiency models [27].

Super-efficiency models rank DMUs by assigning an efficiency score greater than 1. A higher efficiency score indicates a more efficient DMU. Tone et al. proposed the non-radial model, which uses a slack-based measure (SBM) to evaluate efficiency [27]. One SBM-type model is the super-efficiency slack-based measure model in input-oriented CCR (Super-SBM-I-C) model, which handles input or output slacks directly and does not assume proportional changes in inputs or outputs. The details of the model are given in the literature [27]. In this study, the CCR and Super-SBM-I-C models were used to evaluate the efficiency of the 45 cases and find the optimal air supply conditions.

3. Results and Analysis

The distributions of the air temperature, humidity, air velocity, and MAA at different Z_m when v = 12 m/s and t_s = 25 °C were present first. Second, the influence of different air supply conditions on CE and VE was considered by changing v and t_s. All cases in Table 6 were simulated, and their corresponding dry bulb temperature, relative humidity, airflow velocity, MAA, CE, and VE results were obtained. The simulation results were analyzed at 1.7 m (i.e., head level for standing occupants on the Z axis) [28], which is denoted as Plane A and the occupied zone was defined as the zone with dimensions of X = 28–38 m, Y = 1–3.6 m and Z = 0–2.2 m as shown in Figure 8.

3.1. Performance Evaluations with Separate Indices

3.1.1. Cooling Efficiency

Figure 9 shows the dry bulb temperature distributions at Plane A for different Z_m (i.e., 5, 10, 15, 20, and 25 m). In the occupied zone, the air temperature gradually increased with Z_m. When Z_m = 5 and 10 m, the average air temperature in the occupied zone was kept within 30 °C. When Z_m > 15 m, the average air temperature in the occupied zone exceeded 30 °C. Because the outlet of the air duct was close to the heading face, the blocking effect of the LHD was avoided and cooled air could be delivered to the occupied zone [12]. No significant temperature difference was observed in the unoccupied zone when Z_m = 5, 10, 15, or 20 m. This indicates that the location of the air duct outlet had little impact on the thermal environment of the unoccupied zone within 20 m of the heading face. This is because the airflow in the unoccupied zone was mainly in the form of reflux and a similar heat transfer occurred between the cool air from the air duct and the LHD/strata, which resulted in little difference in the reflux temperature. When Z_m = 25 m, the airflow temperature in the unoccupied zone differed significantly from that of the other cases. Because the air duct outlet was far from the heading face, the cool air from the air duct formed a non-isothermal jet that exchanged heat with the roof of the strata. This caused the cool air to increase in temperature before entering the heading face and reduced the cooling effect. Afterwards, the airflow carried away the heat produced by the LHD, which resulted in a higher temperature in the unoccupied zone. Thus, the positions of the LHD and air outlet were important factors affecting the air temperature distribution. This is in line with the findings of Kurnia et al. [9]. Figure 10 shows the distributions of the relative humidity at Plane A under different Z_m. Since there is no obvious moisture source in the heading, all the cases show relatively low relative humidity, at about 50–55%. As mentioned above, the decrease in Z_m would reduce the airflow dry bulb temperature in the occupied area. Hence, the relative humidity increased instead, such as in the case of Z_m = 5 m, because of the few variations in the air moisture content. With increase in the airflow temperature, the relative humidity gradually decreased along the axial length of the heading (from the heading face to the outlet of the airway) in all cases.

Figure 10 shows the distributions of the relative humidity at Plane A under different Z_m. Since there is no obvious moisture source in the heading, all the cases show relatively low relative humidity, at about 50–55%. As mentioned above, the decrease in Z_m would reduce the airflow dry bulb temperature in the occupied area. Hence, the relative humidity increased instead, such as in the case of Z_m = 5 m, because of the few variations in the air moisture content. With increase in the airflow temperature, the relative humidity gradually decreased along the axial length of the heading (from the heading face to the outlet of the airway) in all cases.

Figure 11a compares the CE values for different Z_m and v. According to Equation (19), CE of more than 100% means the average humidex of the unoccupied zone is higher than that of the occupied zone. Except for Z_m = 10 m, the other cases showed a trend of CE increasing along with v. However, the growth rate of each case was different. For example, when Z_m = 10 m, CE significantly increased when v was increased from 8 m/s to 10 m/s but not when v was increased from 10 m/s to 12 m/s. CE was maximized when v = 10 m/s and started to decrease at higher v. This indicates that CE does not improve indefinitely with an increasing air supply velocity; a critical air supply velocity exists. When v is higher than the critical air supply velocity, CE may become worse, which is in line with the study of Gan et al. [25]. Thus, the critical air supply velocity needs to be determined when designing a forcing auxiliary ventilation system to reduce the ventilation energy consumption.

Figure 11b compares the CE values for different Z_m and t_s. CE was best at Z_m = 10 m, especially at lower t_s. However, when the temperature increased to 25 °C, CE clearly decreased at Z_m = 10 m. This may indicate that CE decreases when the humidex difference between the unoccupied and occupied zones is reduced. At Z_m = 15 m, three different t_s values had no evident effect on CE. Therefore, to ensure a high CE, v and t_s can be set to reasonable ranges of values for a certain Z_m. In conclusion, considering only CE, the optimum air supply conditions for this case are Z_m = 10 m, v = 10 m/s, and t_s = 21 °C.

3.1.2. Ventilation Efficiency

Figure 12 shows the airflow velocity distributions for different Z_m. All cases showed that the air velocity from the outlet of the air duct decreased with increasing distance and dissipated more after each contact with the rib forced the air to change course. The exception was when the truck blocked the airflow, which resulted in eddies in front of the truck. Eddies were also observed on all four sides of the loader. These eddies resulted from air recirculating from the air duct, as demonstrated by the streamlines. This is consistent with the results of Zheng and Tien [10]. When Z_m = 5, 10, and 15 m, a high-velocity airflow (>1 m/s) appeared around the driver’s cab. According to Roghanchi et al. [29], this airflow velocity is suitable for mining operation. The cool airflow can effectively carry away heat generated by the strata, mining equipment, and miners, while the airflow velocity behind the LHD is relatively low [30]. When Z_m = 20 and 25 m, the zones of high-velocity airflow appeared around the tail of the LHD, whereas the airflow velocity was low at the heading face. This is dangerous because it causes gases and dust to accumulate. At Z_m = 25 m, the LHD blocked the path of fresh air to the heading face, so the airflow could not effectively remove the heat transferred from the heading face.

Figure 13 shows the MAA distributions for different Z_m. In the occupied zone, MAA increased with Z_m. This is because, when the air duct was located close to the heading face, fresh air could flow directly into the occupied zone, which reduced MAA. A small MAA is desirable, especially in the occupied zone, because it indicates that the fresh airflow can effectively reduce the DPM concentration in the space between the rib and tailpipe and the driver’s cab. When Z_m = 5, 10, 15, or 20 m, no significant difference in MAA was observed in the unoccupied zone. At Z_m = 25 m, MAA showed more evident variation in the unoccupied zone. This is because, when the air duct was far from the heading face, most of the fresh airflow went into the unoccupied zone, which is in line with the findings of Xin et al. [19]. Because of the increase in MAA, the occupied zone was filled with air with serious diesel exhaust pollution when Z_m = 25 m.

Figure 14 shows the effect of v on VE at different Z_m. VE decreased slightly with increasing v. This indicates that there is a threshold for the effectiveness of increasing v, which is important for optimizing v in a design. VE was better at Z_m = 5 and 10 m than the other cases. VE was poor at Z_m = 25 m. Figure 14 shows that, if only VE is considered, the optimum air supply conditions are Z_m = 10 m and v = 8 m/s. In summary, if CE and VE are evaluated individually, two different sets of air supply conditions can be obtained. The highest CE was obtained at t_s = 21 °C, which required the largest cooler to provide the lowest t_s values. Therefore, this case does not necessarily provide the most efficient air supply. To explore the most efficient air supply strategy, the DEA model was applied to evaluate the efficiency of 45 ventilation cases.

3.2. Overall Performance Evaluation with the DEA Model

DEA Solver 5.0 was used to evaluate and rank the efficiency of 45 cases with the CCR and Super-SBM-I-C models. The evaluation model is shown in Figure 6, and the evaluation results are presented in

Table 7. CCR efficiency scores and ranks of 45 cases.

DMU (Case)	Zm (m)	Input		Output		CCR-Model
		v (m/s)	ts (°C)	CE/%	VE/%	Score	Rank
1.1	5	8	21	127.1	78.8	0.941	16
1.2	5	8	23	125.4	81.15	0.966	12
1.3	5	8	25	121.8	84.2	0.999	6
1.4	5	10	21	143.8	77.95	0.885	28
1.5	5	10	23	141.1	79.45	0.943	15
1.6	5	10	25	138.3	81.15	0.998	7
1.7	5	12	21	142.7	79.95	0.858	33
1.8	5	12	23	142	80.15	0.925	22
1.9	5	12	25	141.3	80.75	1	1
2.1	10	8	21	158.9	82.25	1	1
2.2	10	8	23	129	81.4	0.971	9
2.3	10	8	25	125.4	84.3	1	1
2.4	10	10	21	144.3	78.2	0.888	27
2.5	10	10	23	141.7	79.7	0.947	14
2.6	10	10	25	138.9	81.25	1	1
2.7	10	12	21	150.6	79.85	0.901	26
2.8	10	12	23	150	80.1	0.970	10
2.9	10	12	25	131.2	80.7	0.979	8
3.1	15	8	21	125	75.75	0.907	24
3.2	15	8	23	121.6	78.65	0.936	18
3.3	15	8	25	118.3	81.7	0.969	11
3.4	15	10	21	128.8	73.3	0.802	39
3.5	15	10	23	128.4	74.7	0.865	31
3.6	15	10	25	125	76.3	0.928	20
3.7	15	12	21	145.2	76.1	0.869	29
3.8	15	12	23	145.1	76.3	0.937	17
3.9	15	12	25	144.4	76.95	1	1
4.1	20	8	21	123.5	75.35	0.901	25
4.2	20	8	23	118.7	78.2	0.929	19
4.3	20	8	25	115.2	81.2	0.963	13
4.4	20	10	21	129.5	72.75	0.803	38
4.5	20	10	23	128.2	74	0.862	32
4.6	20	10	25	126.9	75.4	0.924	23
4.7	20	12	21	134.8	72.75	0.808	36
4.8	20	12	23	134.1	72.95	0.868	30
4.9	20	12	25	132.2	73.4	0.927	21
5.1	25	8	21	83.3	65.6	0.778	40
5.2	25	8	23	79.9	67.75	0.803	37
5.3	25	8	25	77.9	69.4	0.823	34
5.4	25	10	21	99.3	63.4	0.645	45
5.5	25	10	23	98.6	64	0.705	43
5.6	25	10	25	97.2	65.25	0.774	41
5.7	25	12	21	113.3	66.9	0.690	44
5.8	25	12	23	113.5	67.15	0.758	42
5.9	25	12	25	111.8	67.45	0.823	35

. A score of 1 indicates an efficient DMU located on the efficiency frontier. A score of less than 1 indicates an inefficient DMU with room for further optimization. Table 7 also indicates that five DMUs, or 11.1% of the cases, were CCR–DEA-efficient for reducing energy consumption. These were Cases 1.9, 2.1, 2.3, 2.6, and 3.9 at distances Z_m = 5, 10, and 15 m. This is consistent with China’s safety regulations for coal mines, metal mines, and non-metal mines [31,32]. In contrast, setting Z_m = 20 and 25 m yielded CCR–DEA-inefficient solutions with relatively low DEA scores; this indicates an input excess or output shortage, as in other cases of CCR–DEA inefficiency.

DEA can also propose efficiency optimization strategies for inefficient cases. Take Case 5.4 as an example; the improvement strategy for Case 5.4 is given in

Table 8. Optimization strategies for Case 5.4.

DMU (Case) 5.4	Evaluation Results	Improvement Goals	Difference between Goals and Results
v (m/s)	10.0	6.5	−3.5
ts (°C)	21.0	24.3	3.3
CE/%	99.3	99.3	0
VE/%	63.4	63.4	0

. The evaluation results of the CCR-Model show that, to make Case 5.4 the efficient DMU, v should be reduced from 10 m/s to 6.5 m/s and t_s should be raised from 21.0 °C to 24.3 °C under the conditions of keeping CE and VE unchanged. Hence, in Table 8, 6.5 m/s and 24.3 °C were the improvement goals of v and t_s, respectively, for Case 5.4. The differences between 6.5 m/s to 10.0 m/s and 24.3 °C to 21.0 °C are the differences between goals and results of v and t_s, respectively. Therefore, for Case 5.4, it may be necessary to control the cooling load of the heading roadway in order to achieve high efficiency.

As noted previously, the DEA model can avoid the limitations of individual evaluation indices. Evaluating CE or VE separately could lead to the misleading conclusion that Cases 2.6 and 3.9 are inefficient because their performances are not outstanding. However, the DEA model assigns DMUs with the most favourable indices’ higher weight values and thus can make a more objective evaluation. According to the DEA model, Cases 2.6 and 3.9 provide good overall efficiency. Based on the DEA model, Case 2.6 and Case 3.9 were advantageous for control of the air supply temperature.

As noted earlier, the five DMUs with CCR–DEA-efficiency scores of 1 could not be further ranked by the CCR model. Hence, they were further ranked with the Super-SBM-I-C model, and

Table 9. Evaluation ranking of 5 top DMUs using Super-SBM-I-C model.

DMUs	Input		Output		Super-SBM-I-C
	v (m/s)	ts (°C)	CE/%	VE/%	Score	Rank
1.9	12	25	141.3	80.75	1.004	5
2.1	8	21	158.9	82.25	1.166	1
2.3	8	25	125.4	84.3	1.123	2
2.6	10	25	138.9	81.25	1.013	3
3.9	12	25	144.4	76.95	1.012	4

presents the evaluation results. Case 2.1 obtained the highest Super-SBM-I-C score of 1.166. Thus, the cost–performance ratio for CE and VE was optimized when Z_m = 10 m, v = 8 m/s, and t_s = 23 °C. To keep the air temperature in the occupied zone below 30 °C, t_s should be maintained at 24 °C, which is a critical temperature for high efficiency. This is in line with the standards of ventilation and cooling engineering [33]. For an environmental cooling load of about 65 kW, the air supply volume can be set to 240 m³/min, and the air duct should be within 10 m from the heading face. This will effectively provide the required thermal comfort for operators within 6 m of the heading face. These results show that effectively organizing the airflow pattern for the design of auxiliary ventilation is critical for maintaining comfortable conditions in deep underground mining headings. It is also essential to provide fresh and cool airflow to the occupied zone and employ reasonable v and t_s. The cooled air increases the cooling load and encourages heat flow away from the rock. Moreover, although an excessive v leads to an increase in the ventilation energy consumption, it may not necessarily improve the VE. To sum up, through the DEA model, efficient ventilation strategy can be obtained. The optimal layout of the air duct and air t_s are consistent with the current regulations and related research, which fully shows that the DEA model can be used in the field of ventilation optimization in underground development headings.

4. Conclusions

This study evaluated the performance of a VC system installed in a development heading under different air supply conditions. The following conclusions were obtained:

Detailed distributions of the air flow, temperature, relative humidity, and mean age of air inside the development heading were simulated with CFD at varying v and Z_m. Without considering technical efficiency, the ventilation conditions for obtaining high CE and VE were obtained individually. The critical v (8 m/s) and t_s (24 °C) for the VC system were identified. When v was lower than the critical v, it was observed that ventilation could improve the thermal environment by increasing v and reducing t_s. By contrast, when v was higher than the critical v, we found that the VE was worse. Hence, it is necessary to avoid the occurrence of excessive eddy areas or short-circuiting due to an excessively large v, which affects the air quality of development headings.

The DEA model, with v and t_s as inputs and CE and VE as outputs, was used to evaluate mine ventilation. Multiple ventilation schemes were ranked in terms of technical efficiency, and an economical air supply scheme with a reasonable air duct arrangement and air supply temperature difference of 6 °C was determined. For an environmental cooling load of about 65 kW, the air supply volume can be set to 240 m³/min, and the air duct should be within 10 m from the heading face. The other air supply conditions can be determined based on the psychrometric chart. These results are expected to be significantly beneficial for the design of forcing auxiliary ventilation systems with a heat exchanger in development headings.

This study is a preliminary step in the evaluation of ventilation systems for deep underground mines; hence, the study is associated with a few limitations. First, in our assessment, only cost indicators of the VC system terminal equipment were considered, whereas other indicators, such as the investment for the refrigeration plant, chilled water pipelines, and thermal insulation materials, were not included due to a lack of data. Therefore, the assessment is performed merely from a perspective of v and t_s efficiency, rather than a perspective of overall economic benefit. Second, other auxiliary ventilation systems such as overlap systems are also used, and there may be differences in the efficiencies of these systems. These limitations should be addressed in future research to ensure thermal comfort in the headings of deep underground mines.