Experimental and Numerical Study of Air Flow Reversal Induced by Fire in an Inclined Mine Working

: Effective fire prevention in mine workings and tunnels requires a thorough theoretical analysis of the heat and mass transfer processes within these structures. This involves using established models to calculate non-isothermal air flow dynamics in long tunnels and mine workings. While the ventilation of tunnels has been extensively studied, significant challenges persist regarding mine ventilation systems, particularly due to their complex and branched topology. This study aimed to address these challenges and gaps in mine ventilation. We designed a laboratory bench to simulate an inclined mine working with a heat source (fire) and validated a mathematical model of heat and mass transfer in such settings. Using experimental measurements, we verified the model’s accuracy. It is important to note that our experimental and theoretical analyses focused solely on the thermal effects of a fire, without considering the release of harmful impurities. Using the validated model, we conducted multiparameter simulations to identify the conditions leading to the formation of a thermal slug in an inclined mine working and the subsequent reversal of air flow. The simulation data enabled us to determine the dependency of the critical heat release rate on the aerodynamic parameters of the mine working. Additionally, we evaluated the changes in average air density within a mine working at the critical heat release rate. These findings are crucial for the further development of a network-based method to analyze air flow stability in mine ventilation networks during fires.


Introduction
An important aspect of the safe operation of any underground structures for various purposes is the prevention of fires and effective fire management [1].A particularly high level of danger is posed by fires that occur in long tunnels and branched mining systems, where the objects of fire (transport, conveyor, ventilation equipment, etc.) are capable of releasing large amounts of harmful substances into the environment.In addition to this, high gas concentrations and temperatures near the fire object can lead to strong thermal convection, which not only complicates the ventilation of underground spaces but also makes it unstable [2,3].
Fires in tunnels and mine workings are fairly common phenomena, and have been repeatedly noted in the scientific literature [4,5].As Brake [6] rightly noted, even if all reasonable precautions are taken, there will always be at least a few plausible scenarios that could lead to a major underground fire in any mine.Therefore, important areas of research are not only minimizing the risks of underground fires, but also the development of approaches and methods for localizing and eliminating the consequences of fires [7,8].
The development of new approaches and methods for localizing and eliminating fires in tunnels and mine workings requires an understanding of the patterns of development of underground fires, as well as their influence on the distribution of heat and mass transfer parameters of air flows in tunnel systems or mine workings.The analysis of such patterns has, to date, been largely studied for subway tunnels and other transport systems.The works presented in the literature describe both experimental studies of the dynamics of air flow parameters in the presence of powerful heat sources [9][10][11] and theoretical calculations using modern models and numerical methods.When implementing the latter, three-dimensional models of heat and mass transfer processes in individual sections of tunnels are usually utilized [12,13].When applied to mines, researchers also use onedimensional models of mine ventilation networks [14,15].This is important for assessing the potential gas contamination of the mine as a whole as a result of a long-term fire, as well as for developing measures for the evacuation of miners [16].
One of the most important issues when analyzing fires in tunnels and mine workings is to understand the direction of movement of the hot gas-air mixture.Zhou and Smith [17] identified the possible stages of loss of stability of the gas-air flow in tunnels and the occurrence of reverse flows due to gas or thermal convection: (1) Unidirectional flow; (2) Partial return flows; (3) Completely reverse flows.
The authors noted that the occurrence of one or another stage depends on the mean air velocity in the cross-section of the mine workings.The concept of critical air velocity, which is key in determining the stages of development of return air flows during a fire, is discussed in this study.
Another important parameter is the critical heat release rate (HRR) of the fire, which is usually associated with the critical air velocity by one or another empirical relationship.For example, Oka and Atkinson [9] proposed the following dependence, obtained experimentally from a model bench simulating a tunnel: V c = 0.35q * 0.3 , q * ≤ 0.12 0.35, q * > 0.12 (1) where q * = Q/ ρ 0 C p T 0 g 1/2 H 5/2 is the dimensionless HRR of the fire; Q is the dimensional HRR of the fire (W); ρ 0 is the inlet air density (kg/m 3 ); T 0 is the inlet air temperature ( • C); C p is the specific heat capacity of air (J/(kg• • C)); H is the tunnel height (m); and g is the gravitational acceleration (m/s 2 ).Subsequently, many studies in the literature have been devoted to clarifying this empirical relationship: Atkinson and Wu [18] took into account the influence of the excavation inclination angle, Wu and Bakar [19] proposed determining the dimensionless HRR of a fire through the equivalent diameter of the excavation, Lee and Ryou [20] studied the propagation of smoke in tunnels at different ratios of their transverse dimensions, and Vauquelin [21] studied the effect of channel width on the critical air velocity.The maximum temperature of the gas-air mixture at the roof of the tunnel [22,23] and the length of the return flows of the gas-air mixture [24,25] were also studied.Vauquelin and Telle [26] introduced a new term: confinement velocity.This is the minimum longitudinal velocity induced by extraction that is necessary to prevent smoke layer development after the last exhaust vent has been activated.
Tang et al. [27] studied the influence of the tunnel aspect ratio on the temperature distribution in tunnels.Zhang et al. [28] investigated the thermal effect of the latent heat of vaporization in concrete linings.Zhou et al. [29] utilized computer-vision-based deep learning methods to analyze smoke diffusion characteristics and achieved real-time predictions of smoke temperature fields.Several researchers examined the dynamics of the temperature field and the harmful effects of fire when implementing natural ventilation, which is often used in tunnels [30][31][32].
Most modern work is aimed at studying the characteristics of the spread of heat and smoke during fires in subway and road tunnels, which have their own characteristics in terms of organizing ventilation [32,33].At the same time, mines have their own specifics of ventilation, inclination angles of mine workings, and branched systems of mine workings ventilated by one or more fans [34,35].When calculating the mine ventilation processes in the event of fires, it is important to take into account not only the local effects of thermal convection at the site of the fire but also their subsequent impact on the entire ventilation network as a whole.For example, fires can lead to a significant redistribution of air flows in mine workings due to the effect of thermal draught [36,37].Today, several approaches are being developed to account for the complex processes of return air flows in mine ventilation systems: one-dimensional models of ventilation networks [38], modified one-dimensional models representing individual workings as a network of branches in a ventilation graph [39], as well as CFD simulations [40].
All these approaches have not been sufficiently developed to date and have certain problems (difficulty in parameterizing models, computational complexity, failure to take into account some effects, such as 'thermal slug' of the air flow, etc.).One of the promising and currently undeveloped areas of research is the improvement of operational methods for calculating air distribution during the development of underground fires in mines by introducing effective empirical dependencies into one-dimensional models of ventilation networks to account for the thermal effect of a fire on air flow parameters.Moreover, such empirical dependencies can be obtained based on multiparameter modeling of heat and mass transfer in a relatively small area of the mine near the fire source in a threedimensional setting.Obtaining such dependencies is the global goal of this study.At the same time, in this study, we solved the primary issue related to the validation of the mathematical model by comparing the data of numerical calculations with the data of experimental measurements of the parameters of the air stream on a laboratory bench that we developed, simulating an inclined mine working.

Experimental Study
As part of this study, an aerodynamic test bench was designed and built, which is a physical model of an inclined working (bremsberg or slope) connected to ventilation and haulage drifts (see Figure 1b).Figure 1a  The test bench consists of a steel channel with a rectangular cross-section measuring 300 mm by 500 mm along almost its entire length, representing a scaled-down model of mine working.In contrast, the cross-sectional area of mine workings in modern mines i usually significantly larger.A constant cross-section is maintained everywhere except fo  The test bench consists of a steel channel with a rectangular cross-section measuring 300 mm by 500 mm along almost its entire length, representing a scaled-down model of a mine working.In contrast, the cross-sectional area of mine workings in modern mines is usually significantly larger.A constant cross-section is maintained everywhere except for a small area near the heating source (in the central area of the bench) and near the fan (in the upper part of the bench), where the cross-section transitions from a circular profile at the fan outlet to a rectangular profile of 300 mm by 500 mm.Photos of the assembled test bench are shown in Figure 2. The test bench consists of a steel channel with a rectangular cross-section measurin 300 mm by 500 mm along almost its entire length, representing a scaled-down model of mine working.In contrast, the cross-sectional area of mine workings in modern mines usually significantly larger.A constant cross-section is maintained everywhere except fo a small area near the heating source (in the central area of the bench) and near the fan (i the upper part of the bench), where the cross-section transitions from a circular profile the fan outlet to a rectangular profile of 300 mm by 500 mm.Photos of the assembled te bench are shown in Figure 2. Air movement in the channel is ensured by a duct fan, which allows for the regulatio of performance.In this case, the average air flow velocity in the channel can vary from 0 to 3 m/s.The fan power is 187 W. The head-flow characteristic of the fan is shown in Figu 3. Air movement in the channel is ensured by a duct fan, which allows for the regulation of performance.In this case, the average air flow velocity in the channel can vary from 0.5 to 3 m/s.The fan power is 187 W. The head-flow characteristic of the fan is shown in Figure 3. Air movement in the channel is ensured by a duct fan, which allows for the regulation of performance.In this case, the average air flow velocity in the channel can vary from 0.5 to 3 m/s.The fan power is 187 W. The head-flow characteristic of the fan is shown in Figure 3.The bench also provides the ability to change the angle of inclination of the channel in the range from −15° to 15°, which approximately corresponds to the possible range of inclination angles of real mine workings, with some margin.This paper, however, describes the results of experimental measurements and comparative analysis with model data only at angles of 0° and −13°.The bench also provides the ability to change the angle of inclination of the channel in the range from −15 • to 15 • , which approximately corresponds to the possible range of inclination angles of real mine workings, with some margin.This paper, however, describes the results of experimental measurements and comparative analysis with model data only at angles of 0 • and −13 • .
A heating source was installed in the central part of the channel.It is a ducted electric heater with a maximum HRR of 22 kW.The HRR from the heating source allows for stepwise regulation with a minimum step of 2.4 kW.At the location where the heating source is installed, the cross-section of the channel changes locally due to the design features of the heaters.The release of harmful impurities in the experimental bench has not yet been considered.Thus, the subject of the study was only the thermal effect of the fire.
The bench has two ventilation openings, one located in the upper part and the other in the lower part.The openings are equipped with valves that allow for manual adjustment.For the purposes of the study described in this paper, the valve of opening No. 2 was hermetically sealed, while the valve of opening No. 1 was examined in both the closed and open positions.In addition to the openings, the channel had two more aerodynamic connections with the external atmosphere: the inlet section in its upper part and the outlet section in its lower part.The outlet section of the bench is also equipped with a valve, which can be in an open or partially closed state.
The range of air flow velocities and HRR values at the bench were determined to comply with the conditions of similarity of heat and mass transfer processes between the test bench and real conditions in mines.The scaling of air velocity and heat flow was carried out using the Froude scaling method described in [9]: (2) Appl.Sci.2024, 14, 6840 where V is the air flow velocity (m/s); W is the HRR from the heat source (W); and L is the characteristic size of the channel (m).The indices "b" and "m" indicate the bench and mine working, respectively.If expression (2) follows directly from the condition of conservation of the Froude number, then expression (3) additionally includes the condition of conservation of temperatures [41].An important point is that when choosing the main parameters of the bench based on similarity, the Reynolds number was not taken into account.It is impossible to simultaneously satisfy the Froude and Reynolds criteria without changing the type of medium being studied in the test bench.As noted in [42], when the turbulent flow is fully developed, the viscous stress force becomes less significant for the formation of the large-scale flow structure, i.e., the true value of the Reynolds number becomes less significant.The Froude criterion, which is the ratio of buoyancy forces to flow inertia forces, is key when analyzing large-scale flow structures formed under conditions of thermal convection or an ignition source.However, within the framework of this work, we selected the parameters of the bench to maintain the range of Reynolds numbers at which the regime of fully developed turbulence is realized in the channel (Re > 10, 000).
For the same reason, associated with the low importance of viscous forces, researchers do not use the Rayleigh criterion, which characterizes the ratio of buoyancy forces to viscous forces, when developing smaller benches for analyzing fires in tunnels and mine workings.The Rayleigh criterion, however, is important when considering problems in media moving in laminar or transient flow regimes [43], as well as in problems with larger temperature differences [44].
It is important to note that the HRR from local ignition sources in mine workings significantly depends on time [45].The short duration of peak HRR values in practice leads to an equivalent HRR often being used, which is the average HRR value over a given time period of the fire [46].The equivalent value of HRR from fires involving various mining equipment (loading and hauling machines, drilling rigs, and dump trucks) reaches 8 MW.The equivalent value of HRR from the ignition of extended heat sources, such as conveyor belts, is slightly higher [6], but they were not the subject of this study.For this reason, considering (3), it is sufficient to reproduce the HRR of the heating source as 30.8 kW in a laboratory bench, corresponding to the smallest equivalent hydraulic diameter of mine workings in potash and polymetallic mines in Russia (3.6 m).
This value is slightly higher than the maximum HRR currently achievable in laboratory conditions (22 kW).This indicates that the developed test bench is primarily applicable for analyzing the patterns of heat and mass transfer processes in underground workings with a characteristic cross-sectional diameter of more than 4 m.However, in our opinion, it can also be useful in analyzing air exchange in mine workings with a smaller cross-section, allowing one to qualitatively capture the basic patterns of heat and mass transfer processes.
The following instruments were used to carry out measurements at the bench: 1.

5.
A Leica Disto D2 rangefinder (Leica geosystems, Aarau, Switzerland) with an absolute measurement error of ±0.0015 m.The locations on the test bench where measurements were taken are shown in Figure 4.
The locations on the test bench where measurements were taken are shown in Figure 4.  Measurements of air flow velocity and temperature in each section under study were carried out in three vertical planes, with a vertical measurement step of 5 cm.The distance between the vertical measurement planes in each section was 20 cm.A schematic diagram of the points for measuring air flow parameters is also shown in Figure 4.

Theoretical Study
An experimental study on a laboratory bench was used to validate the mathematical model of heat and mass transfer in an inclined channel (a dimensionally reduced model of a mine working).This model was further needed to perform a series of calculations of the dynamics of air flows when varying the parameters of a mine working, as well as when considering more complex systems of several mine workings.The mathematical model was based on the Reynolds-averaged continuity equation, the Navier-Stokes equations, and the energy equation, written for a three-dimensional flow domain: where ρ is the air density (kg/m 3 ); u is the velocity vector (m/s); t is the time (s); p is the static pressure (Pa); µ is the molecular viscosity (m 2 /s); µ t is the turbulent viscosity (m 2 /s); e is the specific air energy (J/(kg• • C)); T is the temperature ( • C); and λ and λ t are, respectively, the molecular and turbulent thermal conductivities (W/(m• • C)).
To close this system of equations, the Boussinesq hypothesis was used, which assumes that turbulent viscosity can be represented as a scalar quantity.In this case, turbulent thermal conductivity was assumed to be proportional to turbulent viscosity: where k is the turbulent kinetic energy (m 2 /s 2 ); ε is the rate of dissipation of turbulent kinetic energy (m 2 /s 3 ); A 0 and A s are the parameters of the realizable k-epsilon turbulence model [46]; S is the flow deformation rate tensor (1/s); and Sc t = 0.7 is the turbulent Schmidt number.
To determine the parameters k and ε, we used the realizable k-epsilon turbulence model, which includes the following equations: where σ k and σ ε are turbulent analogues of the Prandtl number; S = √ 2S : S is the invariant of the tensor S of flow deformation rates; C 1 , C 1ε , and C 2 are empirical coefficients of the model; G k is the generation of the turbulent kinetic energy due to velocity field inhomogeneities (kg/m/s 3 ); and G b is the generation of the turbulent kinetic energy due to buoyancy forces (kg/m/s 3 ): where β is the coefficient of thermal expansion of air (1/ • C) and Pr t = 0.85 is the turbulent Prandtl number.
The system of Equations ( 4)-( 11) was supplemented with boundary conditions.At the inlet and outlet areas, a fixed pressure condition was set; on the channel walls, a noslip condition and convective heat transfer condition were applied.A fixed heat flux was imposed on the wall of the heat source.Inside the three-dimensional computational domain, a cross-section was selected where a boundary condition of the 'fan' type was specified (see Figure 5).The air flow passing through this cross-section experienced a prescribed pressure jump ∆P(V), determined based on the head-flow characteristic (Figure 3).However, the head-flow characteristic ∆P(V) in the ANSYS model did not coincide with the curve in Figure 3 due to neglecting shock losses at the inlet and outlet of the flow from the bench, as well as due to several other factors.The polynomial function ∆P(V) was chosen so that, at an HRR of zero and a given angle of rotation of the valve at the outlet, the air flow rate in the inclined part of the channel matched that in the experiment.
It should be noted that Figure 5 shows openings No. 1 and No. 2 on opposite sides of the channel.The same orientation was used on the laboratory bench shown earlier in Figure 2. In the laboratory bench, this orientation of openings No. 1 and No. 2 was more convenient from the point of view of ease of access and for performing measurements.
An important aspect during the initial parameterization of the model was the choice of the heat transfer coefficient through the channel wall.Part of the heat imparted to the flow from the heat source is removed through the walls of the channel into the surrounding atmosphere.Moreover, this thermal factor is more pronounced for the section of the channel below the heat source.Therefore, to account for this thermal factor, the heat transfer coefficient was selected to achieve the best agreement between the model and measured values of average air temperatures in sections below the heat source (see Figure 4).The effective value of the heat transfer coefficient was 150 W/(m²• • C).
It should be noted that Figure 5   An important aspect during the initial parameterization of the model was the choice of the heat transfer coefficient through the channel wall.Part of the heat imparted to the flow from the heat source is removed through the walls of the channel into the surrounding atmosphere.Moreover, this thermal factor is more pronounced for the section of the channel below the heat source.Therefore, to account for this thermal factor, the heat transfer coefficient was selected to achieve the best agreement between the model and measured values of average air temperatures in sections below the heat source (see Figure 4).The effective value of the heat transfer coefficient was 150 W/(m²•°C).
The numerical solution of the system of Equations ( 4)-( 11) with specified boundary conditions was carried out in the Fluent module of the ANSYS software package (version The numerical solution of the system of Equations ( 4)- (11) with specified boundary conditions was carried out in the Fluent module of the ANSYS software package (version 2021 R2) using the finite volume method.The velocity and pressure fields were linked using the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm [47,48].Spatial discretization was carried out using schemes of second-order accuracy (for pressure and velocity equations) and first-order accuracy (for equations of transfer of turbulent characteristics).Time discretization was performed using an implicit first-order scheme.Gradients of scalar quantities were calculated using the least squares method.
To speed up the numerical calculations, a parallelization procedure was used on central processor cores (24 cores).The process of finding a solution at each time iteration was also carried out iteratively using the Multigrid Cycles algorithm [49].The subiterative procedure at each time iteration continued until relative residuals of 10 −4 were achieved for all variables.One calculation took from 3 to 16 h, depending on the initial and boundary conditions.
The mesh was generated in the Meshing module of the ANSYS software package.An irregular tetrahedral finite-volume mesh with a prismatic boundary layer near solid walls was constructed.The mesh parameters were selected to ensure the solution's independence from the mesh.In this work, mesh convergence was studied by controlling the pressure drop in the inclined section of the domain for various initial values of the mass flow entering the computational domain.The results of the comparison of the controlled parameters are given in Table 1.
Table 1 shows that when the cell size changes from 0.3 m to 0.2 m, the change in pressure drop in the study area was no more than 3% for the entire range of considered flow rates.Thus, for further calculations, a mesh with a cell size of 0.3 m on the free surface was selected.The final number of cells was 1,586,296.

Model Validation
At the first stage of the study, experimental measurements of the distribution of the air temperature and air velocity in the cross-sections of the test bench were carried out.These data were used to validate the numerical model.The inclined (−13 • ) and horizontal positions of the bench were considered.
Figure 6 shows theoretical curves and experimental points for the longitudinal component of the air flow velocity, while Figure 7 shows the theoretical curves and experimental points for the air flow temperature.The scenario considered had an HRR value equal to 14 kW and the angle of inclination of the channel was equal to 0 • .The measured parameter is marked along the X-axis, and the immersion depth of the sensor, measured from the upper wall of the flow area of the bench, is marked along the Y-axis.The temperature and flow velocity sensors were immersed down through holes in the upper wall of the channel, so the values along the Y-axis are given with a minus sign.The profiles were built for section No. 3 (near the exit from the channel) and section No. 2 (after the heat source).
Considering that during the intense heating of the air flow, the measured values periodically changed over time (i.e., the distribution characteristics of the air flow were non-stationary); Figures 6 and 7 display not only the average values of the measured velocities and temperatures of the air flow, but also the minimum and maximum observed data values for a sample of measurement data at individual points.These are depicted using elements of the box-and-whisker plots.
In addition to a zero-inclination angle, the case of an inclination angle of −13 • (downward air flow in the channel) was studied, and the fan performance was also varied.Figure 8 shows diagrams of the air velocity and temperature in section No. 3 (near the exit from the channel) for experimental points and numerically calculated curves.The HRR of the source was also set to 14 kW.
A comparison of experimental and theoretical distributions of the air flow velocity and temperature was also carried out using other HRR values, but they were not included in this paper.Overall, the obtained data indicate good agreement between the numerical simulation and laboratory experiments at the considered channel inclination angles.Minor quantitative discrepancies between individual experimental points and numerical simulation data were observed, but overall, at a qualitative level, the model curves and experimental points demonstrated similar patterns of changes in air flow velocities and temperatures across the cross-section.The small discrepancies (less than 5%) were likely associated with errors in the measuring equipment, as well as inaccuracies in the selected mathematical model.These errors can arise due to simplifications in geometry (such as when creating boundary surfaces for the heater) and in physics (for example, the use of the RANS approach).Considering that during the intense heating of the air flow, the measured values periodically changed over time (i.e., the distribution characteristics of the air flow were nonstationary); Figures 6 and 7 display not only the average values of the measured velocities and temperatures of the air flow, but also the minimum and maximum observed data values for a sample of measurement data at individual points.These are depicted using elements of the box-and-whisker plots.In addition to a zero-inclination angle, the case of an inclination angle of -13° (downward air flow in the channel) was studied, and the fan performance was also varied.Figure 8 shows diagrams of the air velocity and temperature in section No. 3 (near the exit from the channel) for experimental points and numerically calculated curves.The HRR of the source was also set to 14 kW.
A comparison of experimental and theoretical distributions of the air flow velocity  In both considered sections, there was good agreement between the model and experimental data for the average values of velocity and temperature.This agreement was largely ensured by the preliminary calibration of the mathematical model (selection of the fan operating point and the heat transfer coefficient through the wall of the channel).The similarity in the variation of the studied flow parameters across the channel's cross-sections near the corresponding average sections further validated the accuracy of the selected mathematical model.
The conclusion drawn demonstrates the suitability of RANS turbulence models for analyzing the patterns of steady non-isothermal flows in inclined channels with flow stratification induced by intense heat sources.This conclusion is justified as we analyzed timeaveraged flow characteristics at a distance from the heat source, which are consistent with those in the existing literature [50].However, it is important to note that RANS models may not accurately capture some small-scale structural features of convective flows in inclined channels and can lead to overestimated air temperatures in the vicinity of heat sources [51].In contrast, LES can offer higher accuracy by resolving a broader range of turbulent scales [52,53].Nevertheless, we consider the accuracy provided by RANS models to be sufficient for the purposes of our study.

Numerical Simulation of Air Flow in an Inclined Channel
The validated mathematical model was used to conduct a series of numerical calculations of heat and mass transfer in the channel under several HRR values and fan operating modes.The fan modes were selected to maintain a specified average velocity in the channel, considering an HRR of zero from the heat source and with the valve fully open at the outlet.Six velocity points were chosen for the study.The pressure jump at the fan interface was determined for each velocity point (see Table 2).This information allowed for an assessment of the air resistance of the flow domain.In both considered sections, there was good agreement between the model and experimental data for the average values of velocity and temperature.This agreement was largely ensured by the preliminary calibration of the mathematical model (selection of the fan operating point and the heat transfer coefficient through the wall of the channel).The similarity in the variation of the studied flow parameters across the channel's cross-sections near the corresponding average sections further validated the accuracy of the selected mathematical model.
The conclusion drawn demonstrates the suitability of RANS turbulence models for analyzing the patterns of steady non-isothermal flows in inclined channels with flow stratification induced by intense heat sources.This conclusion is justified as we analyzed time-averaged flow characteristics at a distance from the heat source, which are consistent with those in the existing literature [50].However, it is important to note that RANS models may not accurately capture some small-scale structural features of convective flows in inclined channels and can lead to overestimated air temperatures in the vicinity of heat sources [51].In contrast, LES can offer higher accuracy by resolving a broader range of turbulent scales [52,53].Nevertheless, we consider the accuracy provided by RANS models to be sufficient for the purposes of our study.

Numerical Simulation of Air Flow in an Inclined Channel
The validated mathematical model was used to conduct a series of numerical calculations of heat and mass transfer in the channel under several HRR values and fan operating modes.The fan modes were selected to maintain a specified average velocity in the channel, considering an HRR of zero from the heat source and with the valve fully open at the outlet.Six velocity points were chosen for the study.The pressure jump at the fan interface was determined for each velocity point (see Table 2).This information allowed for an assessment of the air resistance of the flow domain.It should be noted that the pressure differentials presented here are significantly lower than the characteristic fan pressures in Figure 3.There is no contradiction here for two main reasons: (1) the actual operating point of the fan during the experimental measurements was considerably shifted to the right, into a region of lower pressures, which also resulted in lower pressures during the experiments; (2) additionally, the model did not account for shock losses due to the narrowing of the flow entering the fan from the atmosphere, nor for the expansion of the flow leaving the channel through its outlet.
The numerical simulation enabled the identification of the principal flow regimes, as illustrated in Figure 9, for a thermal power of 22 kW and a channel inclination angle of −13 • .The figure depicts the velocity vector component aligned with the channel axis.The color scheme in this figure is designed such that red corresponds to downward air flows driven primarily by the fan pressure, while blue corresponds to return air flows influenced by thermal convection.At this stage of the study, our primary interest lay in determining the point of air flow reversal in an inclined channel under varying heat source powers, different degrees of ventilation valve opening near the channel exit, and various fan operation modes.During numerical calculations, it was assumed that valve at opening No. 1 (see Figure 5) was open.From a practical perspective, the system under consideration resembles parallel mine workings, with one potentially affected by a fire.By adjusting the ventilation valve opening, we varied the air resistance proportions between the two parallel workings, At the maximum pressure drop (21.09Pa) at the outlet section of the channel, a stable unidirectional flow of air was observed from left to right and downwards.As the pressure decreased, return flows began to form in the outlet section, becoming comparable to the forward flow.At a pressure drop of 8.45 Pa, the return flow in the upper part of the channel reached opening No. 1 and exited the channel through this boundary.In this scenario, the return flow of hot air did not propagate further up the channel towards the fan itself.Near the lower section of the channel, the average velocity in both directions over the cross-section approached zero.This indicates a situation resembling a thermal plug, where the inertial force of the flow is, on average, balanced by the buoyancy force.
However, was is also a situation where the return flows reached the fan and completely occupied the cross-section of the domain.An example of this is illustrated by the calculated velocity field for a minimum pressure of 1.06 Pa (see Figure 9c).These three examples demonstrate that the fan pressure, which provides a given air velocity in an inclined channel at W = 0, is an important parameter that significantly influences the flow regime at a non-zero value of W (e.g., direct flow, reversed flow, and thermal slug).Since the prototype of our model is an inclined mine ventilated by the main fan of the mine, such an analysis is very important to us.Of course, in our case, the fan with the parameters from Table 2 should be interpreted as an equivalent, whose effect is equal to the effect of the mine fan minus the losses due to the friction resistance of the remaining network.
At this stage of the study, our primary interest lay in determining the point of air flow reversal in an inclined channel under varying heat source powers, different degrees of ventilation valve opening near the channel exit, and various fan operation modes.During numerical calculations, it was assumed that valve at opening No. 1 (see Figure 5) was open.From a practical perspective, the system under consideration resembles parallel mine workings, with one potentially affected by a fire.By adjusting the ventilation valve opening, we varied the air resistance proportions between the two parallel workings, thereby investigating its influence on the likelihood of air flow reversal in the branch with the fire.
Figure 10 depicts the dependencies of the air mass flow at the outlet on the fan operation mode, ventilation valve opening degree at the outlet section, and HRR from the heat source, which were derived from multiparameter simulations.Positive mass flow values indicate reversed air flow in the inclined channel (from bottom to top), while negative values correspond to air flow in the same direction as at an HRR of zero (from top to bottom).
Figure 10 shows that at a minimum initial pressure drop, the air flow was reversed for all considered values of the HRR.As the initial pressure drop increased, the reversed air flow rate decreased in absolute value and eventually reached zero.The lower the HRR of the heat source, the lower the pressure at which the zero-flow point was reached.With a further increase in the initial pressure drop, the mass flow rate became negative, indicating that within this parameter range, the air flow did not globally change its direction after the heating source was turned on.
An interesting feature of the dependencies presented in Figure 10 is the approximate coincidence of the intersection points of the curves characterizing the degree of opening of the valve of opening No. 3 with the ordinate axis at a fixed HRR from the source.This holds true for all considered HRR values (22,33,44, and 55 kW).Thus, there was a weak dependence of the moment of reversing the air flow on the air resistance of the subsequent ventilation network behind the outlet.That is, the reversal of the air flow depends, to a greater extent, on the power of the fire and the initial pressure drop in the inclined part of the channel.
Based on the data from the multiparameter numerical simulations, we built a linear regression model to predict the mass air flow in the outlet section of the channel using three variable parameters: the initial pressure drop (∆P), the degree of opening of the ventilation valve at the outlet (φ), and the HRR (W).The form of the linear regression model with dimensional coefficients is presented below: Figure 10 shows that at a minimum initial pressure drop, the air flow was reversed for all considered values of the HRR.As the initial pressure drop increased, the reversed air flow rate decreased in absolute value and eventually reached zero.The lower the HRR of the heat source, the lower the pressure at which the zero-flow point was reached.With a further increase in the initial pressure drop, the mass flow rate became negative, indicating that within this parameter range, the air flow did not globally change its direction after the heating source was turned on.
An interesting feature of the dependencies presented in Figure 10 is the approximate coincidence of the intersection points of the curves characterizing the degree of opening of the valve of opening No. 3 with the ordinate axis at a fixed HRR from the source.This holds true for all considered HRR values (22,33,44, and 55 kW).Thus, there was a weak dependence of the moment of reversing the air flow on the air resistance of the subsequent ventilation network behind the outlet.That is, the reversal of the air flow depends, to a greater extent, on the power of the fire and the initial pressure drop in the inclined part of the channel.This model, despite its simplicity, has two advantages.Firstly, it allows for the analysis and comparison of the degree of influence of each model factor on the predicted value . Secondly, linear regression addresses our main interest-studying the conditions for the reversal of an air flow in an inclined channel.
In general, Figure 10 shows that the flow rate was a nonlinear function of the parameters ∆P, φ, and W.However, the proposed linear regression had an acceptable error and was applicable for quantitative analysis in a small region near the intersection of the curves with the y-axis.Considering the goodness-of-fit of the model (12), it should be noted that the linear regression (12) had an R 2 value of 0.8 for the entire range of data points, which is not very high.
The goodness-of-fit could be improved by considering nonlinear regression.However, determining a nonlinear regression that accurately predicts the air mass flow across a wide range of the parameters ∆P, φ, and W is a challenging task.Moreover, it is not of primary interest from a practical standpoint, as our primary goal was to quantitatively capture the moment when the mass flow Q (out) m transitions from positive to negative values.We associate this transition with the reversal of the mass flow in the inclined channel.For this purpose, linear regression was sufficient.In this case, the mass flow rate was calculated integrally over the entire cross-section of the channel.The moment of the onset of the reversal of the mass air flow, therefore, will be earlier than the moment of the appearance of partial return flows, which we did not study here.
The regressors ∆P, φ, and W take only positive values, suggesting that another option is to consider a log-linear regression model of the form: where ∆P 0 = 1 Pa and W 0 =1 kW are reference values for the pressure drop and HRR, introduced to make the logarithmic expression dimensionless.The log-linear regression (13) had a lower R 2 value (0.67) compared to model (12).Therefore, we continued working with model (12).Given the condition Q (out) m = 0, we can express the W value from expression (12).This will be the critical HRR of the fire at which the air flow will reverse; it can be written as a function of the parameters ∆P and φ: From this dependence, it is clear that the degree of valve opening at the outlet affects the critical HRR less than the initial pressure drop in the channel.For example, with an initial pressure drop of 8 Pa, varying the valve opening degree from 10% to 60% changed the critical HRR from 22.3 kW to 18.7 kW.With an initial pressure drop of 13 Pa, varying the valve opening degree in the same range changed the critical HRR from 59.2 kW to 55.5 kW.Based on these quantitative estimates, we also observed good agreement between the predicted data of model (12) and the calculated data in Figure 10a,d.A similar linear regression was built for the air mass flow rate at the inlet zone: Given the condition Q (in) m = 0, we can express the value of W from Equation (15).This represents the HRR of the heat source at which the air flow on the fan will reverse: In this case, the value W * cr is significantly higher than W cr .At an initial pressure drop of 8 Pa, W * cr ranged from 129 to 141 kW as the valve opening degree varied from 10% to 60%.This confirms the intuitive fact that it is much more difficult to upset the air flow from a fan than in one of the branches of a parallel connection located downstream of the fan.
An important consideration is that the obtained dependencies ( 14) and ( 16) were derived under test bench conditions and need to be scaled for real mine workings.According to dependencies (2) and (3), when changing the linear size of the object under study from L m to L e , the HRR should be scaled by a factor of L e L m 2.5 . Given that the pressure drop ∆P is a quadratic function of the mean air velocity in the cross-section, ∆P should be scaled by a factor of L e L m .The parameter φ is relative, and no scaling is needed for it.As a result, dependencies ( 14) and ( 16) take the form: These expressions can be used to estimate the critical HRR leading to the reversal of an air flow in a mine working with a given inclination angle of −13 • as a result of the action of a powerful heat source.
A numerical simulation of air flow reversal in inclined channels with an increased characteristic channel size L confirmed the validity of formulas ( 17) and (18).For example, with a 10-fold increase in the characteristic channel size, formulas ( 17) and ( 18) predicted pressure drop values at a fixed HRR that differed from the numerical model by only 7.8%.
The adaptation of this relationship in the calculation of air distribution in one-dimensional mine ventilation networks involves incorporating critical HRR values into the mass balance equation for each independent network circuit, especially those affected by fire.For instance, in the context of a simplified ventilation network comprising two circuits, one of which features a significant heat source (refer to Figure 11), the following system of equations should be formulated: where Q is the volumetric air flow (m 3 /s); R is the air resistance (N s 2 /m 8 ); H is the fan pressure in branch No. 1 (Pa); h = z I I − z I I I > 0 denotes the height difference between nodes II and III (m); and z is the vertical coordinate (m).Here, the indices for Q, R and ρ indicate the branch number (see Figure 11).The orientation of the branches was chosen such that at an HRR of zero, the air flow rates in all three branches are positive.Nodes I and III are considered atmospheric, meaning that the outgoing flow maintains a specified atmospheric temperature  and density  , regardless of the thermal characteristics of the inflowing air.
As the HRR in branch No. 2 increases, the air density within this branch decreases due to heating.According to a simple advective transport model, the change in air temperature caused by the HRR from a source at the midpoint of the branch can be estimated as follows: where  is the volumetric air flow rate in branch No. 2, calculated upstream of the heating source (m 3 /s).As the flow approaches the point of thermal slug, preceding the reversal of the air flow, the parameter  will approach  .When  is reached, the air flow in branch No. 2 should become zero, indicating that The orientation of the branches was chosen such that at an HRR of zero, the air flow rates in all three branches are positive.Nodes I and III are considered atmospheric, meaning that the outgoing flow maintains a specified atmospheric temperature T 0 and density ρ 0 , regardless of the thermal characteristics of the inflowing air.
As the HRR in branch No. 2 increases, the air density within this branch decreases due to heating.According to a simple advective transport model, the change in air temperature caused by the HRR from a source at the midpoint of the branch can be estimated as follows: where Q 2 is the volumetric air flow rate in branch No. 2, calculated upstream of the heating source (m 3 /s).As the flow approaches the point of thermal slug, preceding the reversal of the air flow, the parameter W will approach W cr .When W cr is reached, the air flow in branch No. 2 should become zero, indicating that In Equation ( 24), we included a factor of ½ because only the lower half of branch No. 2, below the heat source, experiences an increase in air density (see Equation ( 22)).
In our simple advective transport model, the second term on the right-hand side of Equation ( 24) tends towards zero as the denominator approaches zero.However, near the heating source, the air velocity can be significantly non-zero due to convection effects (see Figure 9b).Therefore, the accurate prediction of ρ 2 requires three-dimensional numerical simulation.In this study, the average air density in an inclined channel was computed under various HRR values, fan operational modes, and outlet valve positions (see Figure 12).A pressure drop of 13.62 Pa corresponded most closely to the point where the air flow in the channel began to overturn at an HRR value of 55 kW, as observed in Figure 10.Even with this pressure drop, the air density did not approach zero but remained above 0.8 kg/m³ across all thermal power values.It is noteworthy that as the outlet valve near the outlet zone opened wider, the minimum air density asymptotically approached a value A pressure drop of 13.62 Pa corresponded most closely to the point where the air flow in the channel began to overturn at an HRR value of 55 kW, as observed in Figure 10.Even with this pressure drop, the air density did not approach zero but remained above 0.8 kg/m 3 across all thermal power values.It is noteworthy that as the outlet valve near the outlet zone opened wider, the minimum air density asymptotically approached a value dependent on the HRR according to a nonlinear relationship.The exact nature of this relationship could not be fully determined from the data in Figure 12, as these points did not precisely coincide with the condition of air flow overturning (the intersection of the curves with the ordinate axis in Figure 10).The closest matches occurred at HRR values of 22 kW and 55 kW, where the minimum air density values were approximately 0.96 and 0.8 kg/m 3 , respectively.Considering that the air density was 1.2 kg/m 3 at an HRR of zero, we observed that the minimum average air density decreased with increasing critical HRR at a rate of approximately 0.007 kg/m 3 /kW.This information can be utilized in estimating ρ 2 instead of using the relationship in Equation (24).This provides a lower-bound estimate of the resulting thermal draught in the ventilation network caused by the fire, since we used the maximum angle of inclination of the mine workings and the minimum average air density, which occurs under conditions of minimal air exchange between the fire-affected workings and adjacent areas.

Conclusions
The main results of the study are as follows: 1.A test bench was designed and assembled to study the patterns of heat and mass transfer in inclined mine workings with intense heat sources.2. Experimental measurements of temperature and air flow velocity distributions in the test bench at two angles of inclination, different heat release powers, and varying air flow resistances validated the mathematical model of non-isothermal turbulent air flow in an inclined channel.3. Multi-parameter modeling of air flow in an inclined channel was conducted at different heat release rates of the fire source, fan operating modes (pressure drops), and air resistances.The results were processed and generalized to build a linear regression model, which can determine the mass flow of air in an inclined channel under various parameters.This model was used to determine the dependence of the critical heat release rate on the aerodynamic parameters of real mine workings.4. A weak dependence of the moment of reversing the air flow on the air resistance of the ventilation system surrounding the considered inclined channel when a heat source was revealed.The overturning of the air flow largely depended on the heat release rate of the fire and the initial pressure drop in the inclined part of the test bench.
The main practical value of the results for mine ventilation is that the developed dependencies can be implemented in one-dimensional models of heat and mass transfer in mine ventilation networks and calculation methods.This will enable a more accurate consideration of the thermal effects of a fire on air flow parameters in complex mining systems without significantly increasing the computational complexity.However, further research is needed to analyze the influence of the channel inclination angle on the critical thermal power and to compare the experimental and theoretical results with known data from the literature.

2 Figure 1 .
Figure 1.Schematic diagrams of the mine (a) and aerodynamic test bench (b).

Figure 1 .
Figure 1.Schematic diagrams of the mine (a) and aerodynamic test bench (b).

Figure 1 .
Figure 1.Schematic diagrams of the mine (a) and aerodynamic test bench (b).

Figure 2 .
Figure 2. Photographs of the test bench: general view of the bench tilted at the maximum angle (a local expansion of the flow area near the heater (b), general view of the bench in a horizontal pos tion (c), and the area near the opening no. 1 (d).

Figure 2 .
Figure 2. Photographs of the test bench: general view of the bench tilted at the maximum angle (a), local expansion of the flow area near the heater (b), general view of the bench in a horizontal position (c), and the area near the opening no. 1 (d).

Figure 2 .
Figure 2. Photographs of the test bench: general view of the bench tilted at the maximum angle (a), local expansion of the flow area near the heater (b), general view of the bench in a horizontal position (c), and the area near the opening no. 1 (d).

Figure 3 .
Figure 3.The head-flow characteristic of the fan.

Figure 3 .
Figure 3.The head-flow characteristic of the fan.

Figure 4 .
Figure 4.The locations of measurement points.Measurements of air flow velocity and temperature in each section under study were carried out in three vertical planes, with a vertical measurement step of 5 cm.The distance

Figure 4 .
Figure 4.The locations of measurement points.
shows openings No. 1 and No. 2 on opposite sides of the channel.The same orientation was used on the laboratory bench shown earlier in Figure 2.In the laboratory bench, this orientation of openings No. 1 and No. 2 was more convenient from the point of view of ease of access and for performing measurements.

Figure 5 .
Figure 5. Geometric model created in ANSYS Design Modeler.

Figure 5 .
Figure 5. Geometric model created in ANSYS Design Modeler.
Appl.Sci.2024, 14, x FOR PEER REVIEW 11 of 23 temperature and flow velocity sensors were immersed down through holes in the upper wall of the channel, so the values along the Y-axis are given with a minus sign.The profiles were built for section No. 3 (near the exit from the channel) and section No. 2 (after the heat source).

Figure 6 .
Figure 6.Comparative analysis of calculated profiles and experimentally measured values of the longitudinal component of the air velocity in cross-sections No. 2 (a) and No. 3 (b).

Figure 6 .
Figure 6.Comparative analysis of calculated profiles and experimentally measured values of the longitudinal component of the air velocity in cross-sections No. 2 (a) and No. 3 (b).

23 Figure 7 .
Figure 7. Comparative analysis of calculated profiles and experimentally measured values of the air temperature in cross-sections No. 2 (a) and No. 3 (b).

Figure 7 .
Figure 7. Comparative analysis of calculated profiles and experimentally measured values of the air temperature in cross-sections No. 2 (a) and No. 3 (b).

Figure 8 .
Figure 8. Comparative analysis of calculated profiles and experimentally measured values of the air velocity (a) and temperature (b) in cross-section No. 3 at an inclination angle equal to -13°.

Figure 8 .
Figure 8. Comparative analysis of calculated profiles and experimentally measured values of the air velocity (a) and temperature (b) in cross-section No. 3 at an inclination angle equal to −13 • .

Figure 9 .
Figure 9. Distributions of the longitudinal component of the velocity vector in an inclined channel under different fan operating modes: high performance mode, Δ = 21.09Pa (a), intermediate mode, Δ = 8.45 Pa (b), and low performance mode Δ = 1.06 Pa (c).

Figure 9 .
Figure 9. Distributions of the longitudinal component of the velocity vector in an inclined channel under different fan operating modes: high performance mode, ∆P = 21.09Pa (a), intermediate mode, ∆P = 8.45 Pa (b), and low performance mode ∆P = 1.06 Pa (c).

Figure 10 .
Figure 10.Dependence of the air mass flow at the outlet on the fan operating mode, the valve opening degree at the outlet, and the HRR: W = 22 kW (a), W = 33 kW (b), W = 44 kW (c), and W = 55 kW (d).

Figure 10 .
Figure 10.Dependence of the air mass flow at the outlet on the fan operating mode, the valve opening degree at the outlet, and the HRR: W = 22 kW (a), W = 33 kW (b), W = 44 kW (c), and W = 55 kW (d).

23 Figure 11 .
Figure 11.One-dimensional model of the ventilation network corresponding to the test bench under consideration.

Figure 11 .
Figure 11.One-dimensional model of the ventilation network corresponding to the test bench under consideration.

23 Figure 12 .
Figure 12.Dependencies of average air density on the fan operating mode, the valve opening degree at the outlet, and the HRR: W = 22 kW (a), W = 33 kW (b), W = 44 kW (c), and W = 55 kW (d).

Figure 12 .
Figure 12.Dependencies of average air density on the fan operating mode, the valve opening degree at the outlet, and the HRR: W = 22 kW (a), W = 33 kW (b), W = 44 kW (c), and W = 55 kW (d).

Table 1 .
Analysis of solution independence from the mesh.

Table 2 .
Relationship between air velocity in the channel and pressure jump at the fan interface.