Selection of turbulence model in ventilation modeling for blind stopes

The authors present the results of mathematical modeling of working area ventilation in a blond stope using the methods of computational aerodynamic and different turbulence models. The analytical model is verified by comparing the model and in situ data.


Introduction
Underground excavations in modern mines can stretch for tens and hundreds kilometers summarily. The working areas being the main consumers of fresh air are located very far from access and air supply openings. It becomes increasingly more difficult to feed fresh air in the working areas as fresh air escapes with return ventilation current, and the total air leakage grows with longer distances to the air consumers. Thus, the goal is to make the best use of air in the working areas [1,2]. The first step on the way to the goal is the field research of the dynamics of air flows in a working area, selection of an appropriate mathematical model of ventilation in this areas, and verification of the model by comparing the numerical data and in situ measurements. Later on, the model can be used in calculation of harmful emission sin blasting or due to operation of combustion engines.
In recent years, the dramatic advance in the computing science has enabled aerodynamic investigations using mathematical modeling software products. The modeling stages include: -geometrization of a test object; -discretization of c computational domain into finite elements (mesh generation); -selection of a turbulence model; -formulation of the boundary and initial conditions; -solver setting. By now there is a variety of models to calculate turbulent flows. The classification of the existing approaches to the numerical modeling of turbulent flows is based on the comprehensiveness of resolution of turbulent fluctuations and their energy spectrum. In this regard, three basic approaches are: direct numerical simulation (DNS), large eddy simulation (LES) and Reynolds averaged Navier-Stokes (RANS) [3].
The DNS approach requires numerical resolution of scales of space and time, starting from large to the smallest (Kolmogorov microscales). For this reason, DNS of flows, which is of practical significance in engineering, seems to remain unachievable for a long time, first, due to limited capabilities of modern computers and, second, because of excessiveness of information provided by DNS.
The LES method is in-between the models which use the RANS equations and DNS. The idea to calculate explicitly large scales of turbulence to model smaller eddies using the so-called subgrid-scale  [3,4].
This study used DNS and LES owing to limited computational power and time.

Mathematical model of ventilation in the working area and numerical results
The Navier-Stokes equations make it possible to obtain the transfer equations for average values which arbitrarily fluctuate in the turbulent flow. In particular, the equations for the average velocity components were for the first time derived by Reynolds in 1894 and were named after the researcher. The Reynolds equations are derived using decomposition of the instantaneous velocity value into the sum of the average vale and the fluctuation component.
The Reynolds decomposition and the averaging procedure help derive a system of equations for an average field of velocities. This system is named the Reynolds averaged Navier-Stokes (RANS). The RANS system is open, and it is necessary to find the turbulence stresses for the system closure. This problem is unsolvable with only RANS equations and Reynolds decomposition. Actually, using initial and averaged Navier-Stokes equations, the transfer equations can be derived for the turbulence stresses. At the same time, these equations contain higher-order correlations which remain unknown (and, thus, need to be closed as well).
The closure problem requires a physically justified mathematical model to calculate turbulence stresses. This is achieved using a turbulence model [3].
Semi-empirical turbulence models in the modern software products can be grouped into: -the algebraic models; -the models with one differential equation of transfer of turbulence characteristics; -the models with two differential transfer equations (two-parameter models); -the models with many equations. This study discusses only two-parameter turbulence models: standard k-ε, k-ω, RNG k-ε and SST.
The standard k-ε model allows calculating some turbulent flows at a sufficient accuracy for engineering and at the acceptable convergence. This model uses the transfer equation of the turbulent kinetic energy k and its dissipation rate ε. The model sufficiently describes properties of free shear flows. The advantages of the standard k-ε model are low accuracy in modeling flows with detachment from surfaces and the necessity of of special approaches to calculating wall flows [5].
The shortcomings of the standard k-ε models are absent in the turbulence model k-ω [6]. It uses two transfer equations: one for the turbulent kinetic energy k, and the other for the turbulent fluctuation frequency ω. This model yields the best solution in the wall layer. The disadvantage of the k-ω model, as against the k-ε model, is a strong dependence of the calculations on the values of ω set in the inlet section.
The RNG k-ε model is obtained from the Renormalization Group theory and is free from some shortages of the initial model. The transfer equations for the kinetic turbulent energy k and its dissipation rate ε have the analogous structure, while the added analytical dependence for the turbulence viscosity improves calculation accuracy at small Reynolds numbers. The experience of using the RNG k-ε model, as compared with the standard k-ε model) shows better agreement of the calculation and experiment results for some flows, in particular, at relatively small Reynolds numbers, high curvature of flow lines and in the ranges of high deformations of the velocity field [3,4].
The Shear Stress Transport model (SST) combines advantages of the k-ε and k-ω models. This is a superposition of the k-ε and k-ω models multiplied by the weight functions F1 and (1 − F1), respectively. The function F1 is constructed so that to be equal to one at the upper limit of the boundary layer and to vanish while approaching the wall. Thus, the SST model uses the modified k-ω model, meant for the description of the large-scale coherent structure, in the interior domain and the standard k-ε model, intended for the small-scale turbulence resolution, in the exterior domain [4]. The SST model is proved to be efficient in calculations of detached flows with small detachment zone [7,8].  3 We tested applicability of these two-parameter turbulence models in calculation of air distribution within a working area of a blind stope after expansion of an opening cut in retreat mining. Figure 1 shows the working area ventilation scheme and the computational domain geometry. The initial data for the model were obtained during field research in a working area of NorNickel's Oktyabrsky Mine. Air flow rates were measured in the opening cut, in the check sections, in intake and return air currents. Fresh air to ventilate the working area in the opening cut at the boundary with the blind stope is boosted by a local ventilation fan via the flow pipeline. At the outlet of the pipeline having diameter of 0.8 m, the average air flow velocity is set as 10 m/s. Return leaves the computational domain along the surface at which the zero static pressure is set. On the sidewalls, adherence is assigned. The boundary layer calculation takes into account high roughness of the sidewalls, which is equal to 0.03 m. the opening cut has a section of 30 m 2 and a length of 26 m; the stope has a section of 150 m 2 and a length of 25 m.
The computational domain was split into grid of tetrahedral elements. First, some calculations were performed on grids of different density toward independence of the solution from the cell size of the grid. Finally, we selected the grid with the cell size of 0.1 m, with increased density nearby the walls.
At the wall a prismatic boundary layer was set using the formula which allows analytical determination of height of the first wall element depending on the conditions imposed on a dimensionless parameter [9]: where is the density of air; is the dynamic viscosity of air; is the shear velocity. The value of should range from 30 to 300 for the standard k-ε model and standard wall functions [10] and should not exceed 3 for the k-ω model.
The calculations were performed using the finite volume method SIMPLE [9,10]. The iterative procedure was run until all relative closure error of the unknown parameters of air were less than 10 -4 .
The total number of elements in each model was around 7 millions, and the number of nodes was round 3 millions. The time step was selected automatically using CFL-Based so that to satisfy the