Simultaneous CFD evaluation of wind flow and dust emission in open storage piles

Applied Mathematical Modelling 33 (2009) 3197–3207

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Simultaneous CFD evaluation of wind flow and dust emission in open storage piles I. Diego, A. Pelegry *, S. Torno, J. Toraño, M. Menendez Mining and Civil Works Research Group, Department of Exploitation and Exploration of Mines, School of Mines – Oviedo University, Independencia 13, 33004 Asturias, Spain

a r t i c l e

i n f o

Article history: Received 2 March 2007 Received in revised form 11 October 2008 Accepted 14 October 2008 Available online 7 November 2008

Keywords: Emission factor Coal pile Computational fluid dynamics (CFD) Wind blown particles

a b s t r a c t Dust emission from storage yards is a multivariable problem to be solved not only at any new installation in order to obtain the licenses from the involved authorities but also at existing yards to continue the operation. Engineers have a great variety of methodologies available at the market to estimate such emissions, but in general the process is divided into two independent stages: wind flow analysis and application of emission rates into such wind pattern. This paper summarizes the research developed by this group to link both steps: by using CFX version 10.0, a powerful computational fluid dynamics (CFD) software, the wind flow around the piles is predicted, or even through a complex yard, and at the same time by implementing new subroutines introduced into the standard software, the program is able to give a quantitative evaluation of the total fugitive dust. Ó 2008 Elsevier Inc. All rights reserved.

1. Introduction Several researchers have developed multitude of techniques and methodologies to analyse the wind flow around the piles of open storage yards [1–4] and it is possible to buy standardised software to work out from quite simple to highly complicated simulations: ISC3, Aermod, Fluent, CFX, etc. After knowing the pattern of the wind speeds over the emission surfaces, such parameters are applied into well known standards like EPA, so designers are able to evaluate the potential dust emission of the different alternatives of storage pile distribution. All these calculations need a great effort and coordination between several disciplines. The possibility of obtaining qualitative conclusions is present all along the works but only when solving the final step of total dust generated, quantitative conclusions would be reachable. This limitation of the studies means a time consuming discussion of the different alternatives when designing a new installation. A very accurate and flexible tool for simulating the wind effect around piles is Computational Fluid Dynamic software; Ansys CFX 10.0 is one commercial application. This paper explains how the researchers have developed a further algorithm which is introduced into the standard software to enable the computational fluid dynamics (CFD) code to calculate the dust emission based on the wind vector presence over each portion of the surface and the characteristics of the material released to the atmosphere. In this particular case, the algorithm introduced was based on the emission rates given by the USEPA [5]. Of course, it is obvious that the algorithm is open to be defined based on any emission factor either obtained from on-site experiments or from existing research work.

* Corresponding author. Tel.: +34 985 10 42 54. E-mail address: [email protected] (A. Pelegry). 0307-904X/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2008.10.037

3198

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

2. Mathematical background The commercial general purpose CFX 10.0 software is a second order, pressure/velocity coupled, finite element based control volume method that uses an unstructured grid and a coupled algebraic multi-grid solver. A detailed description of the relevant theory and solution procedure can be found in the Ansys CFX 10.0 user manual [6]. The set of equations which describe the processes of momentum, heat and mass transfer in a moving fluid are known as the Navier–Stokes equations. These partial differential equations were derived in the early 19th century and have no known general analytical solution but can be discretised and solved numerically. These equations are three The continuity equation

oq þ r  ðqUÞ ¼ 0 ot

ð1Þ

The momentum equations

oðqUÞ þ r  ðqU  UÞ ¼ rp þ r  s þ SM ; ot

ð2Þ

where the stress tensor, s, is related to the strain rate by



2 3



s ¼ l rU þ ðrUÞT  dr  U :

ð3Þ

The total energy equation

oðqhtot Þ op  þ r  ðqUhtot Þ ¼ r  ðkrTÞ þ r  ðU  sÞ þ U  SM þ SE ; ot ot

ð4Þ

where htot is the total enthalpy, related to the static enthalpy h(T,p) by

1 htot ¼ h þ U 2 : 2

ð5Þ

The term r  (U  s) represents the work due to viscous stresses and is called the viscous work term. The term USM represents the work due to external momentum sources and is currently neglected. Notation used refers to: q Density r Gradient r divergence operator U Vector of velocity  tensor product p pressure s stress tensor Momentum source SM l Molecular (dynamic) viscosity d Kronecker’s delta function (identity matrix) Total Enthalpy htot k Thermal conductivity T Temperature Energy source SE H Static enthalpy In our case several simplifications can be done, as we will consider an incompressible and Newtonian fluid and small temperature differences. In principle, the Navier–Stokes equations describe both laminar and turbulent flows without the need for additional information. However, turbulent flows at realistic Reynold’s numbers span a large range of turbulent length and time scales, and would generally involve length scales much smaller than the smallest finite volume mesh, which can be practically used in a numerical analysis. The Direct Numerical Simulation (DNS) of these flows would require computing power which is many orders of magnitude higher than available in the foreseeable future. RANS methods (Reynolds Averaged Navier Stokes) show a good compromise between accuracy and calculation effort, being widely used in most of the engineering applications. These are based on describing the unsteady eddies that form the turbulence by their mean effects on the flow, through the Reynold’s stresses. The authors have tested several RANS methods and have selected one of them, the k-epsilon model, according to studies shown in Ref. [3] This CFD program has two approaches to simulate a case where two materials/fluids are present (Multiphase flow): Eulerian–Eulerian multiphase model and the Lagrangian Particle Tracking multiphase model. For the present work, authors have selected the Lagrangian Particle Tracking model so the full particulate phase is modelled by just a sample of individual particles. The tracking is carried out by forming a set of ordinary differential equations in

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

3199

time for each particle. Each particle is injected to obtain an average of all particle tracks. Because each particle is tracked from its injection point to final destination, the tracking procedure is applicable to steady state flow analysis. Several different forces affect the motion of a particle in a fluid. Consider a discrete particle travelling in a continuous fluid medium. In the CFD code, the equation of motion for such a particle was derived as follows (the fluid density is considered much lower than the particles density):

mp

dv p 1 2 ¼ pqf d C D jv f  v p jðv f  v p Þ þ F b ; |{z} 8 dt |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} II

ð6Þ

I

where mp is the particle mass, d is the particle diameter, v is velocity, q is density, CD is the drag coefficient, and Fb is the buoyancy force due to gravity, the subscript f refers to the fluid and the subscript p refers to the particle. The term on the right-hand side of Eq. (6) is a summation of all of the forces acting on the particle expressed in terms of the particle acceleration.  Term I is the drag force acting on the particle.  Term II is the buoyancy force due to gravity. As referred above, medium complexity turbulence model k-epsilon, with a logarithmic wind profile was selected for the present work. Meshes used in these studies are unstructured, composed by tetrahedral and prismatic elements. For all models developed, the close-to-wall area (boundary layer) of up to approximately 3 m above the ground and pile was meshed by using prisms while the rest of the surrounding air was divided into tetrahedral elements. Modelling with hybrid tetrahedral grids with near surface prism layers resulted in a smaller analysis model with better convergence of the solution and better analysis results. A medium resolution mesh was selected thinking on future models of complex geometry. Several options exist for the specification of turbulence quantities at Inlets. Selected one was the medium option, which fixes the turbulence intensity at 0.05 and the turbulence viscosity ratio at 10. For this work convergence is considered acceptable when root mean square normalized value of the residuals is under 105. According to the CFD Solver Manual values at this level indicate tight convergence. For all models, between 100 and 200 iterations were needed to achieve this convergence level, with clean descent curves of the residuals (no oscillatory behaviour observed, which is related to numerical instability), and the time step selection was left in ‘‘Auto” mode, thus selected by the software. The computer used to solve the model was equipped with a Q6800 Intel Core 2 Extreme processor @2.93 GHz with 4 cores and 4 Gb of RAM, over Windows XP 64 bits SP2. Solving times in the meshes finally used in the single pile simulations (around 300,000 elements) were 25 min. Mesh partitioning was tested but no significant time saving was observed in these medium-sized simulations. 3. Usepa explanation As mentioned earlier, the new algorithm implemented within the CFD code is based on the USEPA, therefore a brief explanation of this standard is given below. The emission factor for surface airborne dust subject to disturbances may be obtained in units of grams per square metre (g/m2) per year as follows:

Emission factor ¼ K

N X

Pi ;

ð7Þ

i¼1

where K is the particle size multiplier, N is the number of disturbances per year, and Pi is the erosion potential corresponding to the observed fastest mile of wind for the ith period between disturbances, calculated by Eq. (8a & 8b)

P ¼ 58ðu  ut Þ2 þ 25ðu  ut Þ; P ¼ 0 for u 6 ut

ð8a & 8bÞ

where u* is the friction velocity (m/s), and ut is the threshold friction velocity (m/s). Therefore, the erosion is directly affected by particle size distribution as well as by the frequency of disturbance [7]; because of the nonlinear form of the erosion, each disturbance must be treated separately. 4. Inclusion of emission routine in CFD software The CFD software enables the engineer to select one out of three methods to inject particles: from an inlet, from regions or from a Source Point [6].

3200

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

Optionally the user can define a custom injection region introducing a routine into the program giving the following parameters:     

Particle position (X,Y,Z). Particle Mass Flow Rate. Mean particle diameter. Velocity (u, v, w). Temperature.

When meshing the domain of our simulation, the surface of the pile is divided into hundreds or thousands of faces having all of them a defined dimension and location. The Particle Position where dust is lifted off is assigned to the centre point of each face. The quantity of dust lifted off at each face, named as Particle Mass Flow Rate, is defined by an equation obtained for each material by laboratory tests, past experiences and/or bibliography [8–10]. In our case we are taking USEPA Eqs. (8a) and (8b). As we take the methodology of the USEPA, the parameter used to quantify the dust emission is the friction velocity u* above each subdivision of the pile surface. Therefore, during the modelling process, several cm above the pile a control surface was created to measure the wind speed, so each value was assigned to its corresponding face at the pile [11]. Running the CFD Solver module, the new conditional Fortran’s algorithm, based on Eqs. (3a) and (3b), will analyse if the wind present above each face is enough or not to remove the dust [12], and in case it is removed, it will quantify the emission per disturbance [13]. The implementation can be summarized as follows: First, three pieces of Fortran code are prepared in order to let the CFD solver get the position of the particle injection points and the air velocity field around these areas. A fourth piece of code is also prepared to create the (7) and (8a and 8b) mass injection equations. These four pieces of code are compiled to create a dynamic-link library (DLL), which is then linked to the CFD code by modifying a pre-processor file using CCL code (CFX command language). This is, a previously preprocessing file created with the user-friendly GUI (Graphical User Interface), where the geometry and associated mesh, the complete air phase data and partial particle phase data are defined, is manually modified (adding new CCL code in any ASCII editor) in order to set-up the particularized emission routines. Once this is done the Solver module is launched and results are treated as usual using the postprocessing module. 5. Application Piles with cone or oval shapes are widely investigated [14]. Authors were looking for some other contributions so a semicircular pile shape was selected to be analysed using this model to assess its potential contribution to dust lift off. Fig. 1 shows the computational domain used for this case. Close to 300,000 elements were used to represent the pile and assure a high resolution for the calculations. The wind was defined parallel to the X axis with a logarithmic profile and coal selected as the bulk material to be storaged. A Control Surface was inserted around the pile at 34 cm height from pile surface to calculate the wind over the surface following the equivalent procedure used by Stunder and Arya [14].

Fig. 1. Semicircular pile model.

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

3201

Fig. 2(a). Emission of dust.

Fig. 2(b). Emission of dust.

Based on the process explained above, the CFD code determine the wind value over each face, applies the routine for Mass Injection following the USEPA formula (3a & 3b) or any other factor predefined by the engineer and gives an emission factor per face (named Particle Number Rate) as represented on Fig. 2a and b. The Particle Number Rate is the number of physical particles that a computational particle represents, and it is calculated according to Eq. (9)

Particle number rateðparticles=sÞ ¼

mass flow rateðkg=sÞ : number of positions  the mass of the particleðkgÞ

ð9Þ

As it was expected, higher wind speed at the edges of the pile means greater contribution to the total dust emission compared to the centre part because we have assumed that the material is the same all along the pile (segregation of particles by sizes and other phenomenon occurred during stacking process are not taking into consideration in this paper). It is quite simple to obtain the total erosion thanks to an internal function of the software named Function Calculator which takes the number of faces where emission is ongoing and aggregate the individual quantities.

3202

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

Fig. 3. Pile meshing with 14,113 faces.

Table 1 Variation of emissions vs wind speed and mesh quality. Wind speed

PM 30 emission 15.56m/s@10m Ratio CFD/USEPA Emission values 25m/s@10m Ratio CFD/USEPA Emission values

USEPA (kg)

CFD results based on different quality of meshing (kg) CFD 387 faces

CFD 14113 faces

CFD 39533 faces

66.29

108% 71.708

106% 70.171

114% 75.528

390.42

103% 403.179

96% 373.966

101% 394.933

In order to analyze the sensibility of the results when changing the meshing over the pile surface, the same geometry as in Fig. 1 was simulated but using a different quantity of faces. The case studied with 14,113 faces on the pile surface is shown on Fig. 3. Taking values obtained by applying just USEPA procedures onto our pile, their comparison with the results obtained by our CFD routines for different wind conditions and mesh quality is presented in Table 1. For the cases with 387 and 14,113 faces respectively, the difference on dust emission estimation was in the range of less than 10% while the number of faces was increased by almost 37 times, therefore in most of the cases the improvement in accuracy by increasing the number of faces is neither profitable nor justified because of the higher time for calculation and simulation. To prove this statement, further density of mesh was created, taking up to 39,353 faces. The results were again close to USEPA values [15], with better agreement for the higher wind speeds. It is clear that the difference is originated on the interpolation process done by the CFD when calculating the evolution of the wind values from one face to another. USEPA has not such interpolation, and the entire surface of the piles is divided into few areas of constant wind speed ratio. Another advantage of this system is that the CFD procedure enables the user to define the size of the particles present at the surface of the pile, therefore the new routines are applied in parallel for each particle size. It is feasible therefore to study the differences on the behaviour and trajectories of the particles based on the diameter distribution [16]. It is reasonable to foresee that large particles will be deposited a few metres after leaving the pile while fines will be blown far away creating serious troubles to the environment preservation (Fig. 4). 6. Simulation of a real yard In order to check the real applicability of the CFD simulations modified by these new routines to calculate automatically the dust emission, it is necessary to study a case where the geometry has a higher complexity like a real storage yard. But

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

3203

Fig. 4. Dust concentration downstream.

when talking about bulk terminals there is no limit on the available lay outs, however we could consider that one of the most extended geometry for large piles is two parallel oblong sharp or flat-crested form. A coal and iron ore sea terminal has been selected for study. The operation of the terminal is focused to import both materials, unloading up to panamax type vessels by means of two grab unloaders. The cargo is dispatched to the yard through enclosed galleries and the coal is stocked via one combined stacker reclaimer which runs parallel to the piles along the centre rail track. Due to the desired width of the piles (55 m), the use of payloaders is required to extend the material. Because of this disposition, the yard machine would be not able to reclaim all the pile so the amount of dead coal will be significant. Wheel loaders and mobile hoppers are used to load the trucks or trains. This arrangement is shown on Fig. 5 requiring 775,423 elements (prisms and tetrahedral) for meshing. The total surface exposed to the wind per pile is 38,254 m2. The dominant wind was defined at 45° from the top right. Figs. 6a–c represent the top view of the piles and the wind speed contour at 1, 3 and 7 m height being the total height of the piles selected as 16 m. The normalized wind speed is defined as us/ur, where us is the wind speed near the pile surface and ur is the wind speed at the height of 10 m, in this case 15 m/s. The contour of us/ur around the pile at 22 cm over the surface shows that the downwind pile has lower wind speeds thank to the presence of the other pile [17] (Fig. 7).

Fig. 5. Pile configuration for a coal yard.

3204

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

Fig. 6(a). Wind speed at 1 m height.

Fig. 6(b). Wind speed at 3 m height.

Fig. 6(c). Wind speed at 7 m height.

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

3205

From the total dust emission (approximately 700 kg for each disturbance) pile 1 emits 5% more than pile 2 due to its greater exposure to the wind. Fig. 8 shows in different colours the particles lifted off from each pile and how the trajectories are totally different. It is also interesting to note that in case of a surface treatment (like wet conditions due to water sprays), the particles blown from pile 1 could be entrapped at pile 2. This behaviour could be simulated by the CFD by using a control parameter named restitution which calls the program to keep frozen the particles which fall down to the ground level. Fig. 9 shows how a significant share of the particles removed from pile 1 are not transported far away when a restitution coefficient of zero is used so the particles are considered fixed to the floor when they land on the ground or other pile. The opposite situation, when the particle trajectory is considered under restitution coefficient value equal to one, the particles have much longer trajectories because touching the ground does not mean that they are trapped. Applying the new subroutines following USEPA (8a & 8b), several simulations under different wind speeds and with three alternative values for the threshold friction velocity were performed and results show that for lower wind speeds and greater threshold velocities no dust particles are lifted off from the surface (Table 2).

Fig. 7. Normalized wind contour.

Fig. 8. Trajectories when particles are entrapped.

3206

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207

Fig. 9. Trajectories when particles are not entrapped.

Table 2 Overall yard dust emissions. Wind@10 m

Threshold friction velocity 0.5

1

1.5

5 10 15 20

Total emission (kg) 9.22 860.06 2995.3 6332.77

0 21 693.57 2545.53

0 0 35.34 628.86

7. Conclusions Complexity of the dust emission calculation for medium and high size yards forces the use of multidisciplinary teams of engineers to obtain wind flow pattern for different wind events, emission rates and total dust lift off. Standardized powerful software is easily accessible at the market to predict wind mean flow changes but it is not programmed to give the dust lifted off. The present paper has demonstrated how new subroutines could be implemented to estimate the wind vector over each of the surface areas in which the piles are divided and at a second stage calculate the contribution of each of those faces to the total dust emission following an algorithm predefined by the user. Finally real yard behaviour was studied and the arrangement of two parallel piles confirms how the upwind pile creates a shadow downwind reducing the wind speed over the other pile so there is a different contribution from each pile to the total dust eroded. Therefore, if the piles are placed perpendicular to the dominant wind direction, the total emission will be lower because one pile will help the other one. Acknowledgements The authors would like to thank CFX User’s support for their contribution. Also special thanks to the Ministry of Education and Science of the State of Spain as this work was funded by grant CTM2005-00187/TECBO through a subcontract with the University of Oviedo: ‘‘Particle Atmospheric Contamination: Prediction Models and Prevention Systems in Industrial Environment”. References [1] A. Borges, D. Viegas, Shelter effects on a row of coal piles to prevent wind erosion, J. Wind Eng. Indust. Aerodynam. 29 (1988) 145–154. [2] J.B. Smitham, S.K. Nicol, Physico-chemical principles controlling the emission of dust from coal stockpiles, Powder Technol. 64 (1991) 259–270.

I. Diego et al. / Applied Mathematical Modelling 33 (2009) 3197–3207 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

3207

J.A. Toraño et al, Appl. Math. Modell. (2006). J. Yoon, V. Patel, Numerical model of turbulent flow over a sand dune, J. Hydraulic Eng. 122 (1996) 10–18. US EPA Sections 13.2.4 and 13.2.5, Compilation of air pollutant emission factors. ANSYS CFX – Solver. Release 10.0. J. Xuan, Turbulence factors for threshold velocity and emission rate of atmospheric mineral dust, Atmos. Environ. 38 (2004) 1777–1783. Dale A. Gillette, John Adams, Albert Endo, Dudley Smith, Rolf Kihl, Threshold velocities for input soil particles into the air by desert soils, J. Geophys. Res. 85 (C10) (1980) 5621–5630. P.J. Witt, K. Carey, T. Nguyen, Prediction of dust loss from conveyors using CFD modelling, Appl. Math. Modell. 26 (2) (2002) 297–309. S.T. Parker, R.P. Kinnersley, A computational and wind tunnel study of particle dry deposition in complex topography, Atmos. Environ. 38 (2004) 3867– 3878. J. Toraño, R. Rodriguez, I. Diego, Surface velocity contour analysis in the airborne dust generation due to open storage piles, Eur. Conf. Comput. Fluid Dynam. (ECCOMAS CFD) (2006). K.S. Hayden, K. Park, J.S. Curtis, Effect of particle characteristics on particle pickup velocity, Powder Technol. 131 (2003) 7–14. D. Gillette, A wind tunnel simulation of the erosion of soil: effect of soil texture, sandblasting, wind speed, and soil consolidation on dust production, Atmos. Environ. 12 (1978) 1735–1743. B.J.B. Stunder, S.P.S. Arya, Windbreak effectiveness for storage pile fugitive dust control: a wind tunnel study, J. Air Pollution Control Assoc. 38 (1988) 135–143. G.E. Muleski, Review of Surface Coal Mining Emission Factors, EPA-454/R-95-007. U.S. Environmental Protection Agency, Research Triangle Park, NC, 1991. S. Lee, K. Park, C. Park, Wind tunnel observations about the shelter effect of porous fences on the sand particle movements, Atmos. Environ. 36 (2002) 1453–1463. J. Xuan, A. Robins, The effects of turbulence and complex terrain on dust emissions and depositions from coal stockpiles, Atmos. Environ. 28 (1994) 1951–1960.