- Email: [email protected]
Impact of vegetation on dustiness produced by surface coal mine in North Bohemia
Computers and Mathematics with Applications 78 (2019) 3175–3186
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Impact of vegetation on dustiness produced by surface coal mine in North Bohemia ∗
Hynek Řezníček , Luděk Beneš CTU in Prague, Faculty of Mechanical Engineering, Karlovo namesti 13, 121 35 Prague, Czech Republic
article
info
Article history: Available online 3 June 2019 Keywords: Computational Fluid Dynamics (CFD) Particulate matter (PM) concentration Atmospheric boundary layer (ABL) Vegetative barrier
a b s t r a c t The contribution deals with the numerical simulation of the efficiency of the protective vegetative barriers close to an open coal mine. An assignment has come from the Czech mining company (Severočeské doly a.s.) requesting the assessment of the impact of a mine extension on the environment. The flow field and the concentrations are computed on 2D cuts with a real geometry of Bílina coal mine. System of the RANS equations for viscous incompressible flow with variable density in Boussinesq formulation is used for description of the flows. The two equations turbulence model is used for the closure of this set of equations. The transport equations for different diameters of pollutant are solved. The numerical solver is based on finite volume method and AUSM+ up scheme. Advanced vegetation model which includes both influence of vegetation on the flow field as well as processes inside vegetation is adopted. For estimation of the efficiency of barrier several integral and point criteria are used. © 2019 Published by Elsevier Ltd.
1. Introduction The open coal mines represent one of the biggest sources of dust in Europe. Each extension of the mine or change of the technology create a potential risk factor for the environment and thorough studies of its influence are therefore necessary. The dust spreads to the environment has damaging effect on human health [1]. The Czech mining company plan, due to the shift of the excavation area, to move some technologies (conveyors, loaders) closer to the nearby village. One of the protective measures to reduce impact of these changes are vegetative barriers established on the mine border. One part of the comprehensive influence study of these types of barriers for different wind speed, vegetative height and density is summarized in this paper. From the simulation point of view, there are several problematic issues:
• different terrain scales, typically the dimension of the mine is in order of kilometers but the size of technologies or mine steep are in meters or tenth of meters.
• combination of turbulence sources in Atmospheric Boundary Layer (ABL), the standard turbulence generated by roughness and terrain obstacles has to be combined with additional turbulent sources given by the vegetation. In our model this part is described by the more complex model which takes into account the vegetation structure (similar to turbulence generated in the porous medium). • difference between terrain dimensions (meters to kilometers) in comparison to the size of leaves or needles whose structure is essential in capturing the dust particles. ∗ Corresponding author. E-mail address: [email protected] (H. Řezníček). https://doi.org/10.1016/j.camwa.2019.05.011 0898-1221/© 2019 Published by Elsevier Ltd.
3176
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
The correct capture of the processes associated with vegetation, its influence on the flow field on one hand and processes inside the vegetation affecting the deposition on the other hand constitutes key point of these simulations. Many publications focused on the mathematical modeling of the pollutant deposition on the vegetation are available. Among the most notable the following can be included: reviews on the topic of dry deposition on the vegetation [2,3] and [4] or modeling studies [5,6] and [7] dealing with aerosols and collection efficiency. Studies [8] and [9] are dedicated to numerical simulation of similar problems containing the pollutant dispersion from the coal mine. The first one focus on PM retention in the coal mine under worst wind conditions and the second refers about influence of solid obstacles. Unlike our study, the vegetation is not considered there (or it is modeled very simply as homogeneous block without turbulence generation or PM deposition). In our work the advanced deposition model based on Petroff’s work is used [2]. Petroff et al. in the extensive study derived the deposition model supported by a lot of experimental data which takes in the account aerodynamical condition and aerosol size dependence in the forest canopies. Like in the Petroff’s work, in our model the three effects of the vegetation are considered: effect on the air flow, i.e. slowdown or deflection of the flow, influence on the turbulence levels inside and near the vegetation and filtering of the particles present in the flow. The deposition model reflects four main processes leading to particles deposition on the leaves: Brownian diffusion, interception, impaction and gravitational settling. The Petroff’s model was successfully used in our previous simulations, for more details see [10]. The two different approaches are possible for modeling of the pollution dispersion transport—sectional and moment method [11]. The sectional method is preferred in the case when not large amount of PM fractions is needed. Due to the fact three dust fraction only are relevant for the Coal mine company, the sectional method is used. Based on the most common wind orientations blowing from the mine to the village, two preferred directions were chosen. In each direction two parallel 2D cuts were computed to get better representation. The situation in general 3D was simplified to 2D cuts mainly because of following reasons:
• 72 variants differs by the wind direction, wind speed and tree height was necessary to compute. Such number of variants could not be solved in 3D for time reasons.
• The complete 3D geometry of the surface mine was not available due to the private industrial confidentiality of the mining company. We believe that the neglect of 3D effect in the situation can be justified. The most interesting wind directions are almost perpendicular to the mining slope and the situation mostly recall a 2D flow over a mountain ridge. Also the linear dust sources are almost perpendicular to the wind orientation so they can be easily simplified to the point sources (in 2D). The 3D effects of the 3D flow (in general) are for sure present but they should be less significant due to 2D similarity of the situation [12]. 2. Physical and mathematical models 2.1. Fluid flow Reynolds-averaged Navier–Stokes (RANS) equations for viscous, incompressible, turbulent and stratified flow are used to describe the air flow in Atmospheric boundary layer. This set of equations is simplified by the Boussinesq hypothesis. The pressure p and the potential temperature θ are split into two components—a background part in hydrostatic balance and a perturbation: p = p0 + p′ and θ = θ0 + θ ′ . Changes in density are according to Boussinesq approximation negligible everywhere except in gravity term. Resulting system of equations can be written in the following form:
∇ · u = 0, (1) ( ′) ∂u p + (u · ∇ )u − ∇ · (νE ∇ u) = −∇ + g + S u, (2) ∂t ρ0 (ν ) ∂θ ′ E + ∇ · (θ ′ u) = ∇ ∇θ ′ . (3) ∂t Pr Vector u represents velocity, ρ∗ stands for the reference air density (usually near-ground value), νE = ν +νT is the effective kinematic viscosity which is composed of the laminar (molecular) and turbulent kinematic viscosity. Gravitational term ′ is expressed by g = (0, 0, −g θθ ) where g = 9.81 m/s2 is gravitational constant. Term S u represents momentum sink due 0 to vegetation and Pr = 0.75 is the turbulent Prandtl number. 2.2. Turbulence For closure of previous system of equations, k–ϵ turbulence model in following form is employed.
(( ∂ρ k + ∇ · (ρ ku) = ∇ · µ+ ∂t (( ∂ρϵ + ∇ · (ρϵ u) = ∇ · µ+ ∂t
) ) µT ∇ k + Pk − ρϵ + ρ Sk , σk ) ) µT ϵ ϵ2 ∇ϵ + Cϵ1 Pk − Cϵ2 ρ + ρ Sϵ . σϵ k k
(4) (5)
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3177
The production of turbulent kinetic energy caused by main stream interactions is denoted as Pk and µ stands for the laminar (molecular) dynamic viscosity. The model is completed with a constitutional relation for turbulent dynamic 2 viscosity µT = Cµ ρ kϵ . Source terms Sk and Sϵ of k and ϵ respectively are modified as compared to the standard k −ϵ model and consist of two parts Sk = Skr + Skv resp. Sk = Sϵr + Sϵv : a part expressing road traffic influence and a part expressing vegetation influence. The terms Skr , Sϵr modeling the road traffic sources are adopted from [13]. The sinks and sources due to vegetation influence will be described later. Standard setting of the k–ϵ model with following constants was used: σk = 1.0, σϵ = 1.167, Cϵ1 = 1.44, Cϵ2 = 1.92 and Cµ = 0.09. 2.3. Transport of PM concentration The pollutant density in the air is modeled by the equation for passive scalar transport. The equation for each non dimensional mass fraction cj is as follows:
(ν ) ∂ (ρ cj us ) ∂ρ cj T + ∇ · (ρ cj u) − =∇· ∇ρ cj + ρ Fcj + Scj . ∂t ∂y Sc
(6)
Here us stands for the settling velocity, Fc denotes the pollutant source term and Sc is the vegetation deposition term. The turbulent Schmidt number Sc = 0.72 is used according to the review and the discussion in [14]. The settling velocity us of a spherical particle with the diameter d and density ρp is given by the Stokes’ equation see [15]: us = (d2 ρp gCc )/(18µ), with the correction factor Cc = 1 +
λ d
(2.34 + 1.05 exp(−0.39d/λ)) ,
where λ = 0.066 µm is the mean free path of the particle in the air. 2.4. Model of vegetation Suitable model of the vegetation is a crucial point of this computation. Vegetation deforms flow field, increases level of turbulence and plays significant role in the deposition processes. The vegetation barrier is modeled as a porous block described by a so called Leaf Area Density (LAD ) profile which represents foliated surface area per unit volume. In our computation the horizontally homogeneous forest is assumed. Original LAD profile is multiplied by coefficient representing the vegetation density. The coefficient is considered 1 for the new vegetation, but for the old one which is sparse it is set to 0.3. The LAD profile was adopted from [16] for the broadleave trees and corresponds to the expected vegetation. Three effects of vegetation are considered: the first one is the drag induced by the vegetation. It is modeled as momentum sink inside the vegetation in Eq. (2): S u = −Cd LAD |u|u, where Cd = 0.3 is the drag coefficient [17]. The second effect is the influence on the turbulent quantities. Following [17] the source terms in Eqs. (4) and (5) are written Skv = Cd LAD (βp |u|3 − βd |u|k),
ϵ
Sϵv = Cϵ3 Skv . k
The constants are βp = 1.0, βd = 5.1 and Cϵ3 = 0.9. The particle deposition in the vegetation represents the third process. According to the [2], this effect is given by the term Sc = −LAD ud ρp c in Eq. (6). The term is proportional to the deposition velocity ud which reflects four main processes by which particles depose on the leaves: Brownian diffusion, interception, impaction and gravitational settling. Its value generally depends on wind speed, particle size and vegetation properties. In this study the model from [18] is adopted. 3. Numerical method The system of Eqs. (1)–(6) is solved using artificial compressibility method. Eq. (1) is rewritten as evolution equation for pressure as follows: 1 ∂ p′
β ∂t
+ ∇ · u = 0.
The choice of dimensionless parameter β is discussed e.g. in [19], in our computations it is set to the β = 1000.
(7)
3178
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 1. The situation near the village Braňany and the surface coal mine with 2D cuts (blue lines, presented cut in red), black lines represent new position of technologies and a mineroad (■), planned vegetation is displayed by green stripes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The resulting set of equations is discretized using finite volume method on unstructured grid. The AUSM+up scheme [20], fitting for all speeds, is used for convective terms. Second order accuracy is achieved via the linear reconstruction, where gradients are calculated using least squares approach. To prevent artificial overshooting, Venkatakrishnan limiter [21] is utilized. Gradients needed for diffusive terms are evaluated using Gauss–Green theorem on a dual cell associated with the cell-face (diamond type scheme). This discretization results in a set of ordinary differential equations (in time) solved using implicit BDF2 method. In every time step (outer iteration), first the system of the Navier–Stokes equations (2), (3), (7) is solved, followed by the system of the k − ϵ equations ((4), (5)) and then by the system of the passive scalar equations (6). Values of turbulent viscosity, coupling together turbulence equations with the Navier–Stokes equations, are taken from the previous time step. Each of these non-linear systems is solved by the Newton method. Inner linear systems are solved using GMRES solver. The linear systems are preconditioned by ILU(3) pre-conditioner. Necessary evaluations of the Jacobians are done via finite differences. High computational cost of these operations is reduced by two complementing approaches: via matrix coloring, which profits from the sparseness of the Jacobian matrix, and by calculating the pre-conditioner matrices (as well as the Jacobians) only every 20th time step. Since a steady-state solution is sought, the time step is continuously adapted to accelerate the convergence. The adapting criterion is based on the number of the iterations of the linear solvers in one outer iteration. Time stepping proceeds until a steady-state solution is reached. The in-house solver is written in C++. PETSc library [22] is used for the non-linear system solution. The code has been validated on several cases [23]: The ABL flow solver was tested on the benchmark with rising thermal bubble and the measurements for flow over an isolated 2D hill. The ability to describe correctly behavior of vegetative barrier was validated for the flow in and around a forest canopy and flow around a hedgerow. 4. Computational cases The situation and position of the coal mine and its neighborhood is shown in the orthofotomap in Fig. 1. Braňany village is situated close to the border of the mine. Position of the mining technologies which are major dust sources will be moved to the south, closer to the mine border and village. Therefore the mining company plan to plant trees between the edge of the mine and the village to help the environment inside the village. Effect of these new vegetative barriers was studied. For clarity the new positions of the technologies are denoted by black dashed lines and the new vegetation is displayed by the green stripes inside Fig. 1. The problem is generally 3D, but due to the many variants needed to be modeled (different wind speed and direction, tree heights, PM fractions), only 2D cuts are computed. Based on mutual position of the mine and the village, on the meteorological data and terrain configuration, the two most significant wind directions are chosen. In each direction two characteristic 2D cuts are determined where flow field and PM concentrations are computed. The concentrations of three
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3179
Fig. 2. Scheme of the domain (for chosen cut) with main dust sources (belt conveyors for coal (•) and mine roads (■)), X is a horizontal coordinate along the cut.
Fig. 3. Sketch of the grid. Table 1 Position and intensities of point sources + linear intensities of all uncovered areas inside the mine for all the PM fractions. Type
road 1 conveyor 1 conveyor 2 road 2 conveyor 3
Position [m]
926.18 933.92 976.62 986.35 996.19
Source intensity [g/s] PM 2.5
PM 10
PM 75
0.00075 0.0001 0.0139 0.00015 0.0001
0.00740 0.0029 0.2502 0.00148 0.0029
0.01323 0.0058 0.8173 0.00265 0.0058
1.3 × 10−5
2.7 × 10−5
Source intensity [g/(m s)] linear sources
2 × 10−6
PM fractions (2.5, 10, and 75 µm) and two representative wind speeds (1.7 and 5 m/s) are calculated on each cut. The new vegetative barrier is set with two heights 3 m and 15 m representing young and fully grown trees, old vegetation outside the coal mine is also included with current height (10 m). Together with simulation of the current state it represents 72 cases. The results for a north–south cut marked by the red line in Fig. 1 are presented in this paper. The real terrain geometry of chosen 2D cut is sketched in the Fig. 2. Main PM sources are the conveyor belts and the mine roads but the simulation also considers the open surfaces in the mine as dust sources. The dashed lines represent entrance to the village (x = 0 m) and border of the mine. The computational domain is 1200 × 450 m large. The computational grid, shown in Fig. 3, contains about 50 · 103 cells with the greatest size of 10 × 10 m. The cells are refined close to the ground and in the neighborhood of the forests where the smallest cell is about 0.6 m height. Particles with diameter 2.5 µm, 10 µm and 75 µm and density 1000 kg/m3 are modeled. The sources positions and intensities were obtained from the mining company and are summarized in Table 1. All free surface of the coal mine without vegetation is considered as an area source and it is represented as a line (in 2D). The setting of linear intensities is also listed in Table 1. The ABL is considered as weakly stable stratified (∂ T /∂ z = 0) K/m. As was shown in [10], the stratification plays a minor role in this type of computation. LAD profile plotted in Fig. 4 is adopted from [16]. The relative position of profile maximum (z /h) is set to 0.75 and its value is 1.4 m−1 .
3180
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 4. LAD profile.
4.1. Boundary conditions Following boundary conditions are satisfied on borders of the domain. Inlet: the logarithmic wind profile is prescribed with u = uref at height zref = 10 m (uref ∈ {1.7, 5} m/s). A roughness parameter is set to z0 = 0.1 m. The homogeneous Dirichlet b.c. are satisfied for other quantities v , c, θ ′ and pressure is extrapolated from the domain. Top: velocity is given by Dirichlet b.c. as well as concentration and potential temperature fluctuation; the pressure is extrapolated from the domain. Bottom: no-slip b.c. for velocity components are prescribed, homogeneous Neumann b.c. are given for pressure and potential temperature. No re-suspension of the particles fallen on the ground is allowed. Outlet: homogeneous Neumann b.c. has to be fulfilled for velocity components, concentration and potential temperature. Pressure is prescribed by the barometric formula. For the Turbulent quantities boundary conditions and wall functions from [24] are used. 5. Computational results We have evaluated two types of characteristics—point and integral. The point characteristics, e.g. concentrations in the given point, are more interesting in terms of the impact on village and their inhabitants. The integral characteristics, typically mass flows through cut, are significant in evaluation of influence of the vegetation on mine as the volume source of PM. The computations for all cuts show the same tendencies in dispersion and PM settling, however the specific values of the final quantities are slightly different. The results for the point concentration for the parallel cut corresponding to the same wind direction are presented in the Appendix for the comparison. Here follows the results for the chosen representative cut (red in Fig. 1). 5.1. Point concentration Firstly the point concentration of each PM diameter near the village (x = 0 m) is evaluated, the results are plotted in Fig. 5 for wind speed 1.7 m/s and Fig. 6 for 5 m/s. The percentages in the figures are taken from the case without forest (0 m, blue column = 100%). The significant reduction of the concentration is obtained for PM75, approximately by 83% (the concentration is reduced to 17%). It is not surprising, the flow is decelerated inside the forest and these particles fall down. On the other hand, total concentration of these particles is low. For PM10 the total concentration is highest and therefore they are most interesting for us. In this case the efficiency of the barrier is around 29% for low forest and 49% for high one (corresponding concentrations are approximately 71% and 51%). Similar effect is visible for PM2.5 with efficiency 12% and 23%. The influence of the wind speed is negligible in this case. 5.2. Integral characteristics Figs. 7–10 represent integral characteristics of the flow between first 300 m in z-coordinate whereas the ground has a vertical coordinate zground = 97.5 m in this case, meaning the shown integral characteristics covers approximately first 200 m in height above the ground. Fig. 7 shows the integral of the mass flow through 300 m (in z) at the point x = 0 m. This value represents total mass flow from the mine and corresponds to the source intensity. In comparison to the point concentration, the results are
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3181
Fig. 5. Concentration for each PM diameter at points 1.5 m and 10 m above the ground for different tree height (‘‘0 m’’ symbolize the case without new vegetation), uref = 1.7 m/s.
Fig. 6. Concentration for each PM diameter at points 1.5 m and 10 m above the ground for different tree height (‘‘0 m’’ symbolize the case without new vegetation), uref = 5 m/s.
Fig. 7. Integral mass flow through the first 300 m (in z) for different PM fractions and different tree height, uref = 1.7 m/s.
significantly different. Efficiency of the fully grown barrier for PM75 is 78%, for PM10 32% but it is 2% for PM2.5 only. The efficiency of the barrier for reduction of the source intensity is for the lighter particles negligible. It is given by the different behavior of the particles which is clearly demonstrated in Figs. 8 and 9 where the vertical profiles of the mass flow and concentrations are shown. Flow inside the forest is significantly decelerated and heavier particles sediment on leafs or fall down. For lighter particles different processes play role. Most particles do not enter the forest, significant process is deflection of the flow and increase of the turbulence by vegetation, which leads to the spreading of pollutant to the higher parts of the atmosphere. The shift of the mass flow maxims for the smaller particles to the higher levels of atmosphere documents it well. On the vertical profiles of concentrations can be seen the higher
3182
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 8. Mass flows for different PM fractions in vertical cut x = 0 m, uref = 1.7 m/s.
Fig. 9. Vertical profiles of the concentrations for different PM fractions in vertical cut x = 0 m, uref = 1.7 m/s.
concentration in higher part of ABL for the smaller PM when the vegetation is present. This behavior was studied in details earlier and the results are summarized in [25]. The near ground concentrations are reduced in all cases, see Fig. 9, but for the PM2.5 and partly for PM10 the concentration increase in the higher part of the atmosphere. Similar results we can deduce from Fig. 8. Maxims are reduced in all cases, but their position is shifted up for small particles. Near ground concentrations are always significantly reduced, which is favorable for the village. Influence of the mine as the dust source is significant for heavy and middle particles,
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3183
Fig. 10. Horizontal mass flows for different PM fractions in vertical cut, uref = 1.7 m/s.
for lighter ones it is negligible. Light particles are spread to the higher part of the atmosphere and transported to longer distances. Different point of view is represented by Fig. 10, where the horizontal distribution of the vertical mass flow is shown. Significant dependency on the forest height is clearly visible for the PM10 and the PM75. It is particularly important to notice significant decrease of the number of particles in the forest. The PM2.5 particles are transported independently on the tree height without significant reduction. 6. Conclusions A new method for evaluation of the effects of vegetation barriers on pollution dispersion suitable for industrial use was developed and used for modeling of the particles concentration emitted from the coal mine. Effects of the forest height, particle diameter and wind speed in the case of real topography were studied. Basic processes playing important roles in the deposition and transportation of the particles were identified and shown. Near ground concentrations of all PM fractions are significantly reduced. Concentration of the PM2.5 particles is at 77% of the state without barrier, PM10 52% and PM75 17%, which is favorable for the village inhabitants. Drop in PM10 particles is especially significant, as they form the most important component of the source. The dependency on the wind speed does not seem too significant for common velocities (0 − 5) m/s. Each PM fraction can be dangerous for human health, each in their own way. According to [26] the most dangerous for human health are smallest particles which can fly to pulmonary cells and even infiltrate to the blood. The PM10 fraction can caused asthma because it can penetrate to the bronchi and the larger particles are dangerous in high concentrations when upper respiratory tract cannot filtered them all from the air, because they can damage soft tissue. For the mining company is the most significant reduction of the PM in inhabited areas. The proposed measures have the least impact on the lightest PM2.5 particles. But in this case, efficiency of the barrier is 23%. Efficiency could be increased by a change in the position of the forest or the construction of a protective dust blocking fabrics in front of the village, but these measures are not in the competence of the mining company. The situation is quite different with still dangerous PM10 particles. In this case the proposed vegetative barrier will help to reduce the point concentration in the village almost to the half. Influence of the vegetation on the mine as the volume dust source is significant for heavy and middle particles PM75 and PM10, for lighter one is negligible. Light particles are spreading to the higher part of the atmosphere and are transported for longer distances. Source intensity is reduced in the 15 m forest by 2% for PM2.5, 32% for PM10 and 78% for PM75. Due to the particle mass distribution emitted from the mine, the most important are PM10 particles. Reduction by 17% for young vegetation (3 m) and 31% for fully grown one (15 m) is very important. The limitation of the study lies mainly in 2D simplification of the problem. Although the situation is suitable for 2D because the mine slope is almost perpendicular to the wind direction and the mine is much larger than our domain of interest (as was discussed in the introduction). The 3D simulation would be for sure interesting but it depends mainly if the mining company release the full geometry (in the future). It is necessary to say, the results are very sensitive to exact capture of the fluid flow and parametrization of the vegetation plays significant role.
3184
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 11. Concentration for each PM diameter at points 10 m and 1.5 m above the ground for different tree heights (‘‘0 m’’ is the case without new vegetation), uref = 1.7 m/s (parallel cut).
Fig. 12. Vertical profiles of the concentrations (parallel cut).
Acknowledgments This work was supported by project CAAS, Czech Republic No. CZ.02.1.01/0.0/0.0/16_019/0000778 (from Ministry of Education CR) and grant No. SGS16/206/OHK2/3T/12 of Grant Agency from the CTU, Czech Republic in Prague. Our thanks belong to our colleague Viktor Šíp for helping with the coding and for sharing his great experience.
Appendix Results for the cut parallel to the original (see Fig. 1) are shown in Figs. 11–14. The trends are very similar and the vegetation barrier effectiveness differs only in order of percents. However the concentration of each PM fractions are different due to slightly different geometry of the mining slope.
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3185
Fig. 13. Mass flows for different PM fractions (parallel cut).
Fig. 14. Integral mass flow through the first 300 m (in z) for different PM fractions and different tree heights, uref = 1.7 m/s (parallel cut).
References [1] A.K. Patra, S. Gautam, P. Kumar, Emissions and human health impact of particulate matter from surface mining operation-a review, Environ. Technol. Innov. 5 (2016) 233–249. [2] A. Petroff, et al., Aerosol dry deposition on vegetative canopies, Atmos. Environ. 42 (2008) 3654–3683. [3] S. Janhäll, Review on urban vegetation and particle air pollution - deposition and dispersion, Atmos. Environ. 105 (2015) 130–137. [4] T. Litschke, W. Kuttler, On the reduction of urban particle concentration by vegetation - a review, Meteorol. Zeitsch. 17 (2008) 229–240. [5] A. Tiwary, H. Morvanb, J. Colls, Modelling the size-dependent collection efficiency of hedgerows for ambient aerosols, Aerosol Sci. 37 (2005) 990–1015. [6] L. Matejicek, et al., Spatio-temporal modelling of dust transport over surface mining areas and neighbouring residential zones, Sensors 8 (2008) 3830–3847. [7] J. Steffens, Y. Wang, K. Zhang, Exploration of effects of a vegetation barrier on particle size distributions in a near-road environment, Atmos. Environ. 50 (2012) 120–128. [8] S.A. Silvester, I.S. Lowndes, D.M. Hargreaves, A computational study of particulate emissions from an open pit quarry under neutral atmospheric conditions, Atmos. Environ. 43 (40) (2009) 6415–6424. [9] T. Bodnar, L. Benes, K. Kozel, Numerical simulation of flow and pollution dispersion in the area of opencast coal mine, Numer. Anal. Apl. Math. 52 (2008). [10] V. Šíp, L. Beneš, Modelling the effects of a vegetation barrier on road dust dispersion, Appl. Mech. Mater. 821 (2016) 105–112. [11] V. Šíp, L. Beneš, Dry deposition model for a microscale aerosol dispersion solver based on the moment method, J. Aerosol Sci. 107 (2017) 107–122. [12] P.J. Coelho, J.C.F. Pereira, Finite volume computation of the turbulent flow over a hill employing 2D or 3D non-orthogonal collocated grid systems, Int. J. Numer. Methods Fluids 14 (4) (1992) 423–441. [13] D. Bäumer, B. Vogel, F. Fiedler, A new parameterisation of motorway-induced turbulence and its application in a numerical model, Atmos. Environ. 39 (2005) 5750–5759. [14] Y. Tominaga, T. Stathopoulos, Turbulent Schmidt numbers for CFD analysis with various types of flowfield, Atmos. Environ. 41 (2007) 8091–8099. [15] W.C. Hinds, Aerosol technology: Properties, behavior, and measurement of airborne particles, Wiley, 1999. [16] B. Lalic, D.T. Mihailovic, An empirical relation describing leaf-area density inside the forest for environmental modeling, Appl. Math. Comput. 55 (2004) 641–645. [17] G. Katul, et al., One- and two-equation models for canopy turbulence, Bound.-Layer Meteorol. 113 (2004) 81–109.
3186
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
[18] A. Petroff, L. Zhang, et al., An extended dry deposition model for aerosols onto broadleaf canopies, J. Aerosol Sci. 40 (2009) 218–240. [19] F. Muldoon, S. Acharya, A modification of the artificial compressibility algorithm with improved convergence characteristics, Int. J. Numer. Methods Fluids 55 (2007) 307–345. [20] M.-S. Liou, A sequel to AUSM, part II: AUSM+-up for all speeds, J. Comput. Phys. 214 (2006) 137–170. [21] V. Venkatakrishnan, Convergence to steady state solutions of the Euler equations on unstructured grids with limiters, J. Comput. Phys. 118 (1995) 120–130. [22] PETSc web page, https://www.mcs.anl.gov/petsc/. (Accessed on 30 August 2015). [23] V. Šíp, Numerical Simulations of Microscale Atmospheric Flows and Pollution Dispersion (Ph.D. thesis), FME, Czech Technical Uni in Prague, 2017. [24] P. Richards, R. Hoxey, Appropriate boundary conditions for computational wind engineering models using the k-ϵ turbulence model, J. Wind Eng. Ind. Aerodyn. 46 & 47 (1993) 145–153. [25] L. Beneš, V. Šíp, Numerical optimization of near-road vegetation barriers, in: Proceedings of ECCOMAS Congress 2016 - European Congress on Computational Methods in Applied Sciences and Engineering, Crete Island, 2016, pp. 5–10. [26] C.I. Davidson, R.F. Phalen, P.A. Solomon, Airborne particulate matter and human health: a review, Aerosol Sci. Technol. 39 (8) (2005) 737–749.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Impact of vegetation on dustiness produced by surface coal mine in North Bohemia ∗
Hynek Řezníček , Luděk Beneš CTU in Prague, Faculty of Mechanical Engineering, Karlovo namesti 13, 121 35 Prague, Czech Republic
article
info
Article history: Available online 3 June 2019 Keywords: Computational Fluid Dynamics (CFD) Particulate matter (PM) concentration Atmospheric boundary layer (ABL) Vegetative barrier
a b s t r a c t The contribution deals with the numerical simulation of the efficiency of the protective vegetative barriers close to an open coal mine. An assignment has come from the Czech mining company (Severočeské doly a.s.) requesting the assessment of the impact of a mine extension on the environment. The flow field and the concentrations are computed on 2D cuts with a real geometry of Bílina coal mine. System of the RANS equations for viscous incompressible flow with variable density in Boussinesq formulation is used for description of the flows. The two equations turbulence model is used for the closure of this set of equations. The transport equations for different diameters of pollutant are solved. The numerical solver is based on finite volume method and AUSM+ up scheme. Advanced vegetation model which includes both influence of vegetation on the flow field as well as processes inside vegetation is adopted. For estimation of the efficiency of barrier several integral and point criteria are used. © 2019 Published by Elsevier Ltd.
1. Introduction The open coal mines represent one of the biggest sources of dust in Europe. Each extension of the mine or change of the technology create a potential risk factor for the environment and thorough studies of its influence are therefore necessary. The dust spreads to the environment has damaging effect on human health [1]. The Czech mining company plan, due to the shift of the excavation area, to move some technologies (conveyors, loaders) closer to the nearby village. One of the protective measures to reduce impact of these changes are vegetative barriers established on the mine border. One part of the comprehensive influence study of these types of barriers for different wind speed, vegetative height and density is summarized in this paper. From the simulation point of view, there are several problematic issues:
• different terrain scales, typically the dimension of the mine is in order of kilometers but the size of technologies or mine steep are in meters or tenth of meters.
• combination of turbulence sources in Atmospheric Boundary Layer (ABL), the standard turbulence generated by roughness and terrain obstacles has to be combined with additional turbulent sources given by the vegetation. In our model this part is described by the more complex model which takes into account the vegetation structure (similar to turbulence generated in the porous medium). • difference between terrain dimensions (meters to kilometers) in comparison to the size of leaves or needles whose structure is essential in capturing the dust particles. ∗ Corresponding author. E-mail address: [email protected] (H. Řezníček). https://doi.org/10.1016/j.camwa.2019.05.011 0898-1221/© 2019 Published by Elsevier Ltd.
3176
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
The correct capture of the processes associated with vegetation, its influence on the flow field on one hand and processes inside the vegetation affecting the deposition on the other hand constitutes key point of these simulations. Many publications focused on the mathematical modeling of the pollutant deposition on the vegetation are available. Among the most notable the following can be included: reviews on the topic of dry deposition on the vegetation [2,3] and [4] or modeling studies [5,6] and [7] dealing with aerosols and collection efficiency. Studies [8] and [9] are dedicated to numerical simulation of similar problems containing the pollutant dispersion from the coal mine. The first one focus on PM retention in the coal mine under worst wind conditions and the second refers about influence of solid obstacles. Unlike our study, the vegetation is not considered there (or it is modeled very simply as homogeneous block without turbulence generation or PM deposition). In our work the advanced deposition model based on Petroff’s work is used [2]. Petroff et al. in the extensive study derived the deposition model supported by a lot of experimental data which takes in the account aerodynamical condition and aerosol size dependence in the forest canopies. Like in the Petroff’s work, in our model the three effects of the vegetation are considered: effect on the air flow, i.e. slowdown or deflection of the flow, influence on the turbulence levels inside and near the vegetation and filtering of the particles present in the flow. The deposition model reflects four main processes leading to particles deposition on the leaves: Brownian diffusion, interception, impaction and gravitational settling. The Petroff’s model was successfully used in our previous simulations, for more details see [10]. The two different approaches are possible for modeling of the pollution dispersion transport—sectional and moment method [11]. The sectional method is preferred in the case when not large amount of PM fractions is needed. Due to the fact three dust fraction only are relevant for the Coal mine company, the sectional method is used. Based on the most common wind orientations blowing from the mine to the village, two preferred directions were chosen. In each direction two parallel 2D cuts were computed to get better representation. The situation in general 3D was simplified to 2D cuts mainly because of following reasons:
• 72 variants differs by the wind direction, wind speed and tree height was necessary to compute. Such number of variants could not be solved in 3D for time reasons.
• The complete 3D geometry of the surface mine was not available due to the private industrial confidentiality of the mining company. We believe that the neglect of 3D effect in the situation can be justified. The most interesting wind directions are almost perpendicular to the mining slope and the situation mostly recall a 2D flow over a mountain ridge. Also the linear dust sources are almost perpendicular to the wind orientation so they can be easily simplified to the point sources (in 2D). The 3D effects of the 3D flow (in general) are for sure present but they should be less significant due to 2D similarity of the situation [12]. 2. Physical and mathematical models 2.1. Fluid flow Reynolds-averaged Navier–Stokes (RANS) equations for viscous, incompressible, turbulent and stratified flow are used to describe the air flow in Atmospheric boundary layer. This set of equations is simplified by the Boussinesq hypothesis. The pressure p and the potential temperature θ are split into two components—a background part in hydrostatic balance and a perturbation: p = p0 + p′ and θ = θ0 + θ ′ . Changes in density are according to Boussinesq approximation negligible everywhere except in gravity term. Resulting system of equations can be written in the following form:
∇ · u = 0, (1) ( ′) ∂u p + (u · ∇ )u − ∇ · (νE ∇ u) = −∇ + g + S u, (2) ∂t ρ0 (ν ) ∂θ ′ E + ∇ · (θ ′ u) = ∇ ∇θ ′ . (3) ∂t Pr Vector u represents velocity, ρ∗ stands for the reference air density (usually near-ground value), νE = ν +νT is the effective kinematic viscosity which is composed of the laminar (molecular) and turbulent kinematic viscosity. Gravitational term ′ is expressed by g = (0, 0, −g θθ ) where g = 9.81 m/s2 is gravitational constant. Term S u represents momentum sink due 0 to vegetation and Pr = 0.75 is the turbulent Prandtl number. 2.2. Turbulence For closure of previous system of equations, k–ϵ turbulence model in following form is employed.
(( ∂ρ k + ∇ · (ρ ku) = ∇ · µ+ ∂t (( ∂ρϵ + ∇ · (ρϵ u) = ∇ · µ+ ∂t
) ) µT ∇ k + Pk − ρϵ + ρ Sk , σk ) ) µT ϵ ϵ2 ∇ϵ + Cϵ1 Pk − Cϵ2 ρ + ρ Sϵ . σϵ k k
(4) (5)
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3177
The production of turbulent kinetic energy caused by main stream interactions is denoted as Pk and µ stands for the laminar (molecular) dynamic viscosity. The model is completed with a constitutional relation for turbulent dynamic 2 viscosity µT = Cµ ρ kϵ . Source terms Sk and Sϵ of k and ϵ respectively are modified as compared to the standard k −ϵ model and consist of two parts Sk = Skr + Skv resp. Sk = Sϵr + Sϵv : a part expressing road traffic influence and a part expressing vegetation influence. The terms Skr , Sϵr modeling the road traffic sources are adopted from [13]. The sinks and sources due to vegetation influence will be described later. Standard setting of the k–ϵ model with following constants was used: σk = 1.0, σϵ = 1.167, Cϵ1 = 1.44, Cϵ2 = 1.92 and Cµ = 0.09. 2.3. Transport of PM concentration The pollutant density in the air is modeled by the equation for passive scalar transport. The equation for each non dimensional mass fraction cj is as follows:
(ν ) ∂ (ρ cj us ) ∂ρ cj T + ∇ · (ρ cj u) − =∇· ∇ρ cj + ρ Fcj + Scj . ∂t ∂y Sc
(6)
Here us stands for the settling velocity, Fc denotes the pollutant source term and Sc is the vegetation deposition term. The turbulent Schmidt number Sc = 0.72 is used according to the review and the discussion in [14]. The settling velocity us of a spherical particle with the diameter d and density ρp is given by the Stokes’ equation see [15]: us = (d2 ρp gCc )/(18µ), with the correction factor Cc = 1 +
λ d
(2.34 + 1.05 exp(−0.39d/λ)) ,
where λ = 0.066 µm is the mean free path of the particle in the air. 2.4. Model of vegetation Suitable model of the vegetation is a crucial point of this computation. Vegetation deforms flow field, increases level of turbulence and plays significant role in the deposition processes. The vegetation barrier is modeled as a porous block described by a so called Leaf Area Density (LAD ) profile which represents foliated surface area per unit volume. In our computation the horizontally homogeneous forest is assumed. Original LAD profile is multiplied by coefficient representing the vegetation density. The coefficient is considered 1 for the new vegetation, but for the old one which is sparse it is set to 0.3. The LAD profile was adopted from [16] for the broadleave trees and corresponds to the expected vegetation. Three effects of vegetation are considered: the first one is the drag induced by the vegetation. It is modeled as momentum sink inside the vegetation in Eq. (2): S u = −Cd LAD |u|u, where Cd = 0.3 is the drag coefficient [17]. The second effect is the influence on the turbulent quantities. Following [17] the source terms in Eqs. (4) and (5) are written Skv = Cd LAD (βp |u|3 − βd |u|k),
ϵ
Sϵv = Cϵ3 Skv . k
The constants are βp = 1.0, βd = 5.1 and Cϵ3 = 0.9. The particle deposition in the vegetation represents the third process. According to the [2], this effect is given by the term Sc = −LAD ud ρp c in Eq. (6). The term is proportional to the deposition velocity ud which reflects four main processes by which particles depose on the leaves: Brownian diffusion, interception, impaction and gravitational settling. Its value generally depends on wind speed, particle size and vegetation properties. In this study the model from [18] is adopted. 3. Numerical method The system of Eqs. (1)–(6) is solved using artificial compressibility method. Eq. (1) is rewritten as evolution equation for pressure as follows: 1 ∂ p′
β ∂t
+ ∇ · u = 0.
The choice of dimensionless parameter β is discussed e.g. in [19], in our computations it is set to the β = 1000.
(7)
3178
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 1. The situation near the village Braňany and the surface coal mine with 2D cuts (blue lines, presented cut in red), black lines represent new position of technologies and a mineroad (■), planned vegetation is displayed by green stripes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The resulting set of equations is discretized using finite volume method on unstructured grid. The AUSM+up scheme [20], fitting for all speeds, is used for convective terms. Second order accuracy is achieved via the linear reconstruction, where gradients are calculated using least squares approach. To prevent artificial overshooting, Venkatakrishnan limiter [21] is utilized. Gradients needed for diffusive terms are evaluated using Gauss–Green theorem on a dual cell associated with the cell-face (diamond type scheme). This discretization results in a set of ordinary differential equations (in time) solved using implicit BDF2 method. In every time step (outer iteration), first the system of the Navier–Stokes equations (2), (3), (7) is solved, followed by the system of the k − ϵ equations ((4), (5)) and then by the system of the passive scalar equations (6). Values of turbulent viscosity, coupling together turbulence equations with the Navier–Stokes equations, are taken from the previous time step. Each of these non-linear systems is solved by the Newton method. Inner linear systems are solved using GMRES solver. The linear systems are preconditioned by ILU(3) pre-conditioner. Necessary evaluations of the Jacobians are done via finite differences. High computational cost of these operations is reduced by two complementing approaches: via matrix coloring, which profits from the sparseness of the Jacobian matrix, and by calculating the pre-conditioner matrices (as well as the Jacobians) only every 20th time step. Since a steady-state solution is sought, the time step is continuously adapted to accelerate the convergence. The adapting criterion is based on the number of the iterations of the linear solvers in one outer iteration. Time stepping proceeds until a steady-state solution is reached. The in-house solver is written in C++. PETSc library [22] is used for the non-linear system solution. The code has been validated on several cases [23]: The ABL flow solver was tested on the benchmark with rising thermal bubble and the measurements for flow over an isolated 2D hill. The ability to describe correctly behavior of vegetative barrier was validated for the flow in and around a forest canopy and flow around a hedgerow. 4. Computational cases The situation and position of the coal mine and its neighborhood is shown in the orthofotomap in Fig. 1. Braňany village is situated close to the border of the mine. Position of the mining technologies which are major dust sources will be moved to the south, closer to the mine border and village. Therefore the mining company plan to plant trees between the edge of the mine and the village to help the environment inside the village. Effect of these new vegetative barriers was studied. For clarity the new positions of the technologies are denoted by black dashed lines and the new vegetation is displayed by the green stripes inside Fig. 1. The problem is generally 3D, but due to the many variants needed to be modeled (different wind speed and direction, tree heights, PM fractions), only 2D cuts are computed. Based on mutual position of the mine and the village, on the meteorological data and terrain configuration, the two most significant wind directions are chosen. In each direction two characteristic 2D cuts are determined where flow field and PM concentrations are computed. The concentrations of three
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3179
Fig. 2. Scheme of the domain (for chosen cut) with main dust sources (belt conveyors for coal (•) and mine roads (■)), X is a horizontal coordinate along the cut.
Fig. 3. Sketch of the grid. Table 1 Position and intensities of point sources + linear intensities of all uncovered areas inside the mine for all the PM fractions. Type
road 1 conveyor 1 conveyor 2 road 2 conveyor 3
Position [m]
926.18 933.92 976.62 986.35 996.19
Source intensity [g/s] PM 2.5
PM 10
PM 75
0.00075 0.0001 0.0139 0.00015 0.0001
0.00740 0.0029 0.2502 0.00148 0.0029
0.01323 0.0058 0.8173 0.00265 0.0058
1.3 × 10−5
2.7 × 10−5
Source intensity [g/(m s)] linear sources
2 × 10−6
PM fractions (2.5, 10, and 75 µm) and two representative wind speeds (1.7 and 5 m/s) are calculated on each cut. The new vegetative barrier is set with two heights 3 m and 15 m representing young and fully grown trees, old vegetation outside the coal mine is also included with current height (10 m). Together with simulation of the current state it represents 72 cases. The results for a north–south cut marked by the red line in Fig. 1 are presented in this paper. The real terrain geometry of chosen 2D cut is sketched in the Fig. 2. Main PM sources are the conveyor belts and the mine roads but the simulation also considers the open surfaces in the mine as dust sources. The dashed lines represent entrance to the village (x = 0 m) and border of the mine. The computational domain is 1200 × 450 m large. The computational grid, shown in Fig. 3, contains about 50 · 103 cells with the greatest size of 10 × 10 m. The cells are refined close to the ground and in the neighborhood of the forests where the smallest cell is about 0.6 m height. Particles with diameter 2.5 µm, 10 µm and 75 µm and density 1000 kg/m3 are modeled. The sources positions and intensities were obtained from the mining company and are summarized in Table 1. All free surface of the coal mine without vegetation is considered as an area source and it is represented as a line (in 2D). The setting of linear intensities is also listed in Table 1. The ABL is considered as weakly stable stratified (∂ T /∂ z = 0) K/m. As was shown in [10], the stratification plays a minor role in this type of computation. LAD profile plotted in Fig. 4 is adopted from [16]. The relative position of profile maximum (z /h) is set to 0.75 and its value is 1.4 m−1 .
3180
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 4. LAD profile.
4.1. Boundary conditions Following boundary conditions are satisfied on borders of the domain. Inlet: the logarithmic wind profile is prescribed with u = uref at height zref = 10 m (uref ∈ {1.7, 5} m/s). A roughness parameter is set to z0 = 0.1 m. The homogeneous Dirichlet b.c. are satisfied for other quantities v , c, θ ′ and pressure is extrapolated from the domain. Top: velocity is given by Dirichlet b.c. as well as concentration and potential temperature fluctuation; the pressure is extrapolated from the domain. Bottom: no-slip b.c. for velocity components are prescribed, homogeneous Neumann b.c. are given for pressure and potential temperature. No re-suspension of the particles fallen on the ground is allowed. Outlet: homogeneous Neumann b.c. has to be fulfilled for velocity components, concentration and potential temperature. Pressure is prescribed by the barometric formula. For the Turbulent quantities boundary conditions and wall functions from [24] are used. 5. Computational results We have evaluated two types of characteristics—point and integral. The point characteristics, e.g. concentrations in the given point, are more interesting in terms of the impact on village and their inhabitants. The integral characteristics, typically mass flows through cut, are significant in evaluation of influence of the vegetation on mine as the volume source of PM. The computations for all cuts show the same tendencies in dispersion and PM settling, however the specific values of the final quantities are slightly different. The results for the point concentration for the parallel cut corresponding to the same wind direction are presented in the Appendix for the comparison. Here follows the results for the chosen representative cut (red in Fig. 1). 5.1. Point concentration Firstly the point concentration of each PM diameter near the village (x = 0 m) is evaluated, the results are plotted in Fig. 5 for wind speed 1.7 m/s and Fig. 6 for 5 m/s. The percentages in the figures are taken from the case without forest (0 m, blue column = 100%). The significant reduction of the concentration is obtained for PM75, approximately by 83% (the concentration is reduced to 17%). It is not surprising, the flow is decelerated inside the forest and these particles fall down. On the other hand, total concentration of these particles is low. For PM10 the total concentration is highest and therefore they are most interesting for us. In this case the efficiency of the barrier is around 29% for low forest and 49% for high one (corresponding concentrations are approximately 71% and 51%). Similar effect is visible for PM2.5 with efficiency 12% and 23%. The influence of the wind speed is negligible in this case. 5.2. Integral characteristics Figs. 7–10 represent integral characteristics of the flow between first 300 m in z-coordinate whereas the ground has a vertical coordinate zground = 97.5 m in this case, meaning the shown integral characteristics covers approximately first 200 m in height above the ground. Fig. 7 shows the integral of the mass flow through 300 m (in z) at the point x = 0 m. This value represents total mass flow from the mine and corresponds to the source intensity. In comparison to the point concentration, the results are
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3181
Fig. 5. Concentration for each PM diameter at points 1.5 m and 10 m above the ground for different tree height (‘‘0 m’’ symbolize the case without new vegetation), uref = 1.7 m/s.
Fig. 6. Concentration for each PM diameter at points 1.5 m and 10 m above the ground for different tree height (‘‘0 m’’ symbolize the case without new vegetation), uref = 5 m/s.
Fig. 7. Integral mass flow through the first 300 m (in z) for different PM fractions and different tree height, uref = 1.7 m/s.
significantly different. Efficiency of the fully grown barrier for PM75 is 78%, for PM10 32% but it is 2% for PM2.5 only. The efficiency of the barrier for reduction of the source intensity is for the lighter particles negligible. It is given by the different behavior of the particles which is clearly demonstrated in Figs. 8 and 9 where the vertical profiles of the mass flow and concentrations are shown. Flow inside the forest is significantly decelerated and heavier particles sediment on leafs or fall down. For lighter particles different processes play role. Most particles do not enter the forest, significant process is deflection of the flow and increase of the turbulence by vegetation, which leads to the spreading of pollutant to the higher parts of the atmosphere. The shift of the mass flow maxims for the smaller particles to the higher levels of atmosphere documents it well. On the vertical profiles of concentrations can be seen the higher
3182
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 8. Mass flows for different PM fractions in vertical cut x = 0 m, uref = 1.7 m/s.
Fig. 9. Vertical profiles of the concentrations for different PM fractions in vertical cut x = 0 m, uref = 1.7 m/s.
concentration in higher part of ABL for the smaller PM when the vegetation is present. This behavior was studied in details earlier and the results are summarized in [25]. The near ground concentrations are reduced in all cases, see Fig. 9, but for the PM2.5 and partly for PM10 the concentration increase in the higher part of the atmosphere. Similar results we can deduce from Fig. 8. Maxims are reduced in all cases, but their position is shifted up for small particles. Near ground concentrations are always significantly reduced, which is favorable for the village. Influence of the mine as the dust source is significant for heavy and middle particles,
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3183
Fig. 10. Horizontal mass flows for different PM fractions in vertical cut, uref = 1.7 m/s.
for lighter ones it is negligible. Light particles are spread to the higher part of the atmosphere and transported to longer distances. Different point of view is represented by Fig. 10, where the horizontal distribution of the vertical mass flow is shown. Significant dependency on the forest height is clearly visible for the PM10 and the PM75. It is particularly important to notice significant decrease of the number of particles in the forest. The PM2.5 particles are transported independently on the tree height without significant reduction. 6. Conclusions A new method for evaluation of the effects of vegetation barriers on pollution dispersion suitable for industrial use was developed and used for modeling of the particles concentration emitted from the coal mine. Effects of the forest height, particle diameter and wind speed in the case of real topography were studied. Basic processes playing important roles in the deposition and transportation of the particles were identified and shown. Near ground concentrations of all PM fractions are significantly reduced. Concentration of the PM2.5 particles is at 77% of the state without barrier, PM10 52% and PM75 17%, which is favorable for the village inhabitants. Drop in PM10 particles is especially significant, as they form the most important component of the source. The dependency on the wind speed does not seem too significant for common velocities (0 − 5) m/s. Each PM fraction can be dangerous for human health, each in their own way. According to [26] the most dangerous for human health are smallest particles which can fly to pulmonary cells and even infiltrate to the blood. The PM10 fraction can caused asthma because it can penetrate to the bronchi and the larger particles are dangerous in high concentrations when upper respiratory tract cannot filtered them all from the air, because they can damage soft tissue. For the mining company is the most significant reduction of the PM in inhabited areas. The proposed measures have the least impact on the lightest PM2.5 particles. But in this case, efficiency of the barrier is 23%. Efficiency could be increased by a change in the position of the forest or the construction of a protective dust blocking fabrics in front of the village, but these measures are not in the competence of the mining company. The situation is quite different with still dangerous PM10 particles. In this case the proposed vegetative barrier will help to reduce the point concentration in the village almost to the half. Influence of the vegetation on the mine as the volume dust source is significant for heavy and middle particles PM75 and PM10, for lighter one is negligible. Light particles are spreading to the higher part of the atmosphere and are transported for longer distances. Source intensity is reduced in the 15 m forest by 2% for PM2.5, 32% for PM10 and 78% for PM75. Due to the particle mass distribution emitted from the mine, the most important are PM10 particles. Reduction by 17% for young vegetation (3 m) and 31% for fully grown one (15 m) is very important. The limitation of the study lies mainly in 2D simplification of the problem. Although the situation is suitable for 2D because the mine slope is almost perpendicular to the wind direction and the mine is much larger than our domain of interest (as was discussed in the introduction). The 3D simulation would be for sure interesting but it depends mainly if the mining company release the full geometry (in the future). It is necessary to say, the results are very sensitive to exact capture of the fluid flow and parametrization of the vegetation plays significant role.
3184
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
Fig. 11. Concentration for each PM diameter at points 10 m and 1.5 m above the ground for different tree heights (‘‘0 m’’ is the case without new vegetation), uref = 1.7 m/s (parallel cut).
Fig. 12. Vertical profiles of the concentrations (parallel cut).
Acknowledgments This work was supported by project CAAS, Czech Republic No. CZ.02.1.01/0.0/0.0/16_019/0000778 (from Ministry of Education CR) and grant No. SGS16/206/OHK2/3T/12 of Grant Agency from the CTU, Czech Republic in Prague. Our thanks belong to our colleague Viktor Šíp for helping with the coding and for sharing his great experience.
Appendix Results for the cut parallel to the original (see Fig. 1) are shown in Figs. 11–14. The trends are very similar and the vegetation barrier effectiveness differs only in order of percents. However the concentration of each PM fractions are different due to slightly different geometry of the mining slope.
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
3185
Fig. 13. Mass flows for different PM fractions (parallel cut).
Fig. 14. Integral mass flow through the first 300 m (in z) for different PM fractions and different tree heights, uref = 1.7 m/s (parallel cut).
References [1] A.K. Patra, S. Gautam, P. Kumar, Emissions and human health impact of particulate matter from surface mining operation-a review, Environ. Technol. Innov. 5 (2016) 233–249. [2] A. Petroff, et al., Aerosol dry deposition on vegetative canopies, Atmos. Environ. 42 (2008) 3654–3683. [3] S. Janhäll, Review on urban vegetation and particle air pollution - deposition and dispersion, Atmos. Environ. 105 (2015) 130–137. [4] T. Litschke, W. Kuttler, On the reduction of urban particle concentration by vegetation - a review, Meteorol. Zeitsch. 17 (2008) 229–240. [5] A. Tiwary, H. Morvanb, J. Colls, Modelling the size-dependent collection efficiency of hedgerows for ambient aerosols, Aerosol Sci. 37 (2005) 990–1015. [6] L. Matejicek, et al., Spatio-temporal modelling of dust transport over surface mining areas and neighbouring residential zones, Sensors 8 (2008) 3830–3847. [7] J. Steffens, Y. Wang, K. Zhang, Exploration of effects of a vegetation barrier on particle size distributions in a near-road environment, Atmos. Environ. 50 (2012) 120–128. [8] S.A. Silvester, I.S. Lowndes, D.M. Hargreaves, A computational study of particulate emissions from an open pit quarry under neutral atmospheric conditions, Atmos. Environ. 43 (40) (2009) 6415–6424. [9] T. Bodnar, L. Benes, K. Kozel, Numerical simulation of flow and pollution dispersion in the area of opencast coal mine, Numer. Anal. Apl. Math. 52 (2008). [10] V. Šíp, L. Beneš, Modelling the effects of a vegetation barrier on road dust dispersion, Appl. Mech. Mater. 821 (2016) 105–112. [11] V. Šíp, L. Beneš, Dry deposition model for a microscale aerosol dispersion solver based on the moment method, J. Aerosol Sci. 107 (2017) 107–122. [12] P.J. Coelho, J.C.F. Pereira, Finite volume computation of the turbulent flow over a hill employing 2D or 3D non-orthogonal collocated grid systems, Int. J. Numer. Methods Fluids 14 (4) (1992) 423–441. [13] D. Bäumer, B. Vogel, F. Fiedler, A new parameterisation of motorway-induced turbulence and its application in a numerical model, Atmos. Environ. 39 (2005) 5750–5759. [14] Y. Tominaga, T. Stathopoulos, Turbulent Schmidt numbers for CFD analysis with various types of flowfield, Atmos. Environ. 41 (2007) 8091–8099. [15] W.C. Hinds, Aerosol technology: Properties, behavior, and measurement of airborne particles, Wiley, 1999. [16] B. Lalic, D.T. Mihailovic, An empirical relation describing leaf-area density inside the forest for environmental modeling, Appl. Math. Comput. 55 (2004) 641–645. [17] G. Katul, et al., One- and two-equation models for canopy turbulence, Bound.-Layer Meteorol. 113 (2004) 81–109.
3186
H. Řezníček and L. Beneš / Computers and Mathematics with Applications 78 (2019) 3175–3186
[18] A. Petroff, L. Zhang, et al., An extended dry deposition model for aerosols onto broadleaf canopies, J. Aerosol Sci. 40 (2009) 218–240. [19] F. Muldoon, S. Acharya, A modification of the artificial compressibility algorithm with improved convergence characteristics, Int. J. Numer. Methods Fluids 55 (2007) 307–345. [20] M.-S. Liou, A sequel to AUSM, part II: AUSM+-up for all speeds, J. Comput. Phys. 214 (2006) 137–170. [21] V. Venkatakrishnan, Convergence to steady state solutions of the Euler equations on unstructured grids with limiters, J. Comput. Phys. 118 (1995) 120–130. [22] PETSc web page, https://www.mcs.anl.gov/petsc/. (Accessed on 30 August 2015). [23] V. Šíp, Numerical Simulations of Microscale Atmospheric Flows and Pollution Dispersion (Ph.D. thesis), FME, Czech Technical Uni in Prague, 2017. [24] P. Richards, R. Hoxey, Appropriate boundary conditions for computational wind engineering models using the k-ϵ turbulence model, J. Wind Eng. Ind. Aerodyn. 46 & 47 (1993) 145–153. [25] L. Beneš, V. Šíp, Numerical optimization of near-road vegetation barriers, in: Proceedings of ECCOMAS Congress 2016 - European Congress on Computational Methods in Applied Sciences and Engineering, Crete Island, 2016, pp. 5–10. [26] C.I. Davidson, R.F. Phalen, P.A. Solomon, Airborne particulate matter and human health: a review, Aerosol Sci. Technol. 39 (8) (2005) 737–749.