A numerical study for the assessment of air pollutant dispersion with chemical reactions from a thermal power plant

Abstract

In this work, a computational simulation of the pollutants spread generated during fuel combustion at Ekibastuz SDPP-1 and their chemical reaction in the atmosphere have been presented. Using the example of a real thermal power plant (Ekibastuz SDPP -1), the dispersion NO, NO₂, CO and products NO₂, HNO₃, CO₂ during a chemical reaction with oxygen was modeled. The numerical method was verified by solving three test problems and the received computational solutions were compared with data from the measurement and the computational data of other authors. The purpose of the work was to investigate the pollution level at various ranges from the source. As a result, the mass fraction of the concentration and product during a chemical reaction were determined. According to the obtained data, with increasing distance from the source, the concentration of pollution spreads more widely under the influence of diffusion. The farther the distance from the chimney, the lower the concentration of the substance and products. Thus, the obtained computational values may allow in the future predicting the optimal distance from residential areas for the construction of thermal power plants, at which the concentration of emissions will remain at a safe level.

KEYWORDS:

• Chemical reaction
• air pollutant dispersion
• thermal power plant
• jet in crossflow
• Reynolds averaged Navier-Stokes (RANS)

1. Introduction

Atmospheric air is a vital component of the natural environment. The continuous negative impact on the atmosphere and the unsatisfactory resolution of issues related to its improvement negatively affect the health of the population. Air pollution is created from the activities of industrial enterprises, power plants, cars, which emit hundreds of tons of harmful substances into the atmosphere. Photochemical processes constantly occur in the air basin, leading to the emergence of new compounds, sometimes more harmful than the initial ones. Thermal power plants (TPPs) are one of the main sources of environmental pollution (Figure 1). Thermal power plants produce electricity (up to 75% of the total world electricity generation) and thermal energy. When burning, the bulk of the fuel is converted into waste that enters the environment in the gaseous and solid combustion products form. When burning fuel, a large amount of oxygen is consumed, and a significant amount of combustion products are released.

Figure 1. Ekibastuz SDPP-1, Kazakhstan.

Air pollutants can be divided into couple classes. The first class includes primary air pollutants (Akbarian et al., 2018; Ghalandari et al., 2019; Mayer, 1999; Mou et al., 2017), which are ejected into the environment directly from emission source (usually from fuel combustion) and mainly contain of particulate matter components, carbon monoxide (CO), volatile organic compounds (VOCs) and nitrogen oxides (NOx) – mainly NO (Dunmore et al., 2015). The next class is auxiliary air pollutants that form in the environment, while primary air pollutants undergo physical and chemical reactions (Jacobson, 2005). Chemical reactions lead to the dilution of reagents and the formation of secondary pollutants, which are functions of the pollutants concentration, temperature and chemical kinetics (Fiorina et al., 2015; Wu & Liu, 2018). One of the significant auxiliary air pollutant is Ozone (O₃), which is formed as a result of chemical reactions, containing basically the oxidation of NOx and VOCs in the sunlight presence. Reactive intermediates (hydroxyl (OH) and hydroperoxyl (HO₂) radicals) regulate the chemical decomposition cycle, converting NO to NO₂ and, consequently, the formation of O₃ (Bloss, 2009). In the papers (Li, Li et al., 2018, 2019; Muñoz et al., 2018) the formation mechanisms NO₂ in the chemical composition have been numerically investigated, which consisted of 45 species and 142 reactions. The approved model has been applied to form NO₂ from NO. So nitrogen oxides are greenhouse gases, but these pollutants also react with organic compounds to form smog (ground-level ozone). The sulfur dioxide and nitrogen oxides contribute to the appearance of acid rain. The share of combined heat and power (CHP) plant in anthropogenic emissions of these oxides is 15–45% and 45–65%, respectively (U.S. Environmental Protection Agency). When these pollutants mix with oxygen, water and other chemicals in the air, it forms acids. The oxygen is characterized by high chemical activity, and the air contains 21% oxygen. It should be borne in mind that many substances react with oxygen at room temperature. As a result of chemical interaction of substances with oxygen, new harmful ingredients that are more dangerous for humans can be synthesized. Therefore, the problems of atmospheric air pollution and the implementation of measures to clean it remain relevant today.

However, it should be taken into account that if TPPs are located inside or near a city, then one of the main factors affecting the distribution of concentration will be the surrounding buildings and the location of the building (Li et al., 2008b). Since constructions are man-made barriers to urban atmospheric flow and lead to limited ventilation, particularly for deep street canyons (Salim et al., 2011). The deterioration of urban air quality is due to the combined effects of traffic emissions, dynamics of emissions and chemical composition in the environment (Li et al., 2008b). The study on air pollution in street canyons due to their great practical importance has been examined experimentally, numerically and theoretically in many papers (Ahmad et al., 2005; Bright et al., 2013; Kwak et al., 2013; Li et al., 2012; Vardoulakis et al., 2003; Yazid et al., 2014; Zhong et al., 2016). In paper (Wegner & Huai, 2016) a chemical reactive pollutant was numerically investigated that relates to the dynamics of street canyons. Moreover, the LES model was used to study the dispersion and transfer of reactive O₃, NO₂, NO species in the deep canyon of a city street by connecting with simple NOx-O₃ chemical compounds (Li et al., 2008a; Zhong et al., 2015). Furthermore in works (Kikumoto & Ooka, 2012; Kim et al., 2012; Zhong et al., 2015), using the LES method, simulated the effects of bimolecular chemical reactions on concentration and transport pollutant gases to disperse air pollutants in the atmospheric boundary layer.

There are a several methods for solving like that kind of problems, among which one can distinguish a three dimensional experiment, computational simulation and theoretical investigations. The empirical approach includes conducting a real measurement on actual states, using the required devices and fixtures, making measurements of the required values. Nevertheless, this approach is very complex and very expensive way. Moreover, in densely populated places such experiments are extremely inconvenient, because of the air quality in urban conditions basically depends on the wind flow. Theoretical approaches could be determined by methods such as the semi-empirical model and the Gaussian model, but these are unsatisfactory for accurate pollutant prediction. Thus, for a more profitable and efficient investigation of environmental problems, computational simulation is preferred because of its flexibility in calculations, and so it is not so expensive viable method, which makes it much more convenient.

So from these works, the reactivity of air pollutants was revealed that have a significant effect on the concentration value, thus, the concentration of products in the canyon space was higher than the concentration value of the reactants. This effect is due to a stably circulating stream in the canyon, which promotes the reactive gases mixing, and the reactions take place during the time in which the reactors are held. So (Liu et al., 2017) simulated the formation of synthesis gases, such as CO, CO₂, H₂, CH₄ and C₂H₄. The results of transient processes simulation show that an excess of steam accelerates the gasification of coal and accelerates the reaction; therefore, an optimal steam flow should be chosen as a compromise between the reaction time and the quality of synthesis gas. While in the other works (Mărunţălu et al., 2015; Shi et al., 2006, 2007), gas-phase chemical reactions were simulated. Also in paper (Bulat et al., 2014) used the same strategy to study the production of CO and NOx. So the mechanisms and kinetics of gas-phase reactions O₃ have been the subject of numerical investigations over the past years (Atkinson & Carter, 1984; Atkinson, 1989, 1990; Bahta et al., 1984; Bennett et al., 1987; Gillies et al., 1988; Grosjean, 1990; Hatakeyama &  Akimoto, 1990; Horie & Moortgat, 1991; Izumi et al., 1988; Martinez & Herron, 1988; Niki et al., 1983; Nolting et al., 1988; Paulson et al., 1992; Shorees et al., 1991; Yokouchi & Ambe, 1985; Zozom et al., 1987). Currently, one of the main uncertainties relates to reactions in the atmospheric conditions of biradicals, the degree and nature of decomposition products. There is evidence of product and kinetic studies that reactive radical species are formed during these reactions (Atkinson et al., 1990a; Herron & Huie, 1978; Martinez et al., 1981; Martinez & Herron, 1988; Niki et al., 1987; Paulson et al., 1992), and OH radicals play an important role in this chemical reaction.

However, it should be noted that in order for these elements to enter a chemical reaction, certain atmospheric conditions will be necessary. And another thing to remember is certain conditions during the combustion of fuel. So at the moment, various technologies of oxygen combustion are used in order to reduce carbon dioxide (CO₂) emissions from chimneys. This method is also a promising method for eliminating NOx emissions (Andersson & Johnsson, 2006, 2007; Habib et al., 2017). The combustion process using CFD (Al-Abbas et al., 2011) is a very important process, since the investigation of the combustion reaction mechanisms and the emission model are very important for accurate prediction of the combustion characteristics (Ditaranto & Oppelt, 2011; Khalil & Gupta, 2017a; Khalil & Gupta, 2017b; Kim et al., 2007; Li, Shi et al., 2018; Smith et al., 1999). However, taking into account the entire combustion process in modeling is currently very resource-intensive and a set of simplified chemical mechanisms is usually used to reduce computational efforts (Andersen et al., 2009; Frassoldati et al., 2009). However, it should be noted that carbonyl sulphide (COS), carbon disulphide (CS₂) and carbon monoxide (CO) chemical compounds are formed in reactions involving carbon dioxide (CO₂) and hydrocarbons present in the acid gas stream (Goar, 1989; Goar & Hyne, 1996; Maddox, 1998).

A detailed review of works presented the investigating of the substances flow leaving the pipe in a transverse flow was found in work (Denev et al., 2006). In papers (Andreopoulos & Rodi, 1984; Kelso et al., 1996; Muppidi & Mahesh, 1999; Muppidi & Mahesh, 2008) numerically investigated the velocity field, and in works (Crabb et al., 1981; Majander & Siikonen, 2006; Sherif & Pletcher, 1989) modeled a passive scalar field of concentration. In the article (Yuan et al., 1999; Ziefle & Kleiser, 2009), numerical simulation is performed using the large eddy simulation (LES) method, which gives best values than other methods. The simulation was carried out at two ratios of jet velocity to transverse flow, 2.0 and 3.3, and two different Reynolds numbers (1050, 2100) based on the transverse flow velocity and jet diameter. In paper (Keffer & Baines, 1963), jet trajectories and velocity along it were measured. A complete study was conducted in (Crabb et al., 1981), which measured the mean and fluctuating velocity values using a laser Doppler anemometer near the jet exit. In paper (Wegner & Huai, 2004) turbulent mixing was studied using the LES method. They changed the angle between the stream and the cross flow. The mixing intensified as the angle increased, that is, when the jet was directed against the cross flow. The received numerical data were compared with measurement data (Andreopoulos & Rodi, 1984).

To simulate such real size problems, it is necessary to solve test problems. By these test problems were verified the accuracy of the numerical algorithm. For this objective, in this work, 2D and 3D test problems of the pollutants movement which is exiting from the pipe, perpendicular to the main flow in the channel were investigated. The calculations were performed using the ANSYS Fluent. The obtained data were compared with values from the experiment. Satisfactory agreement with data from measurement was obtained. The values showed a closer agreement with the measurement values than the values of simulations by other authors.

Then, modeling of atmospheric pollution at Ekibastuz SDPP-1 was performed (Figure 1). This power plant is the biggest TPP in the Republic of Kazakhstan Ekibastuz city, Pavlodar region, located on the northern side of the Zhyngyldy Lake, 16 km north of Ekibastuz city. The operating capacity of Ekibastuz SDPP-1 is about 3,500 MW. Construction dimensions: height – 64.0 m, length – 500.0 m and width – 132.0 m. Chimney height: 300.0 m (built in 1980) and 330.0 m (built in 1982). The various heights of chimneys allowed studying the influence of the source height on the pollution dispersion.

2. Mathematical model

CFD modeling of these processes is based on the Navier – Stokes equations, consisting of the continuity and momentum equations (Chung, 2002; Ferziger & Peric, 2013; Issakhov, 2014, 2016; Issakhov et al., 2018; Issakhov & Imanberdiyeva, 2019).

• The continuity equation: (1) $\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_j}=0$(1) • The momentum equation: (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho u_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho u_i{u}_j\right)\\ & =-\frac{\mathrm{\partial }{p}^{\mathrm{\prime }}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left({\mu }_{eff}\left(\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}\right)\right)-\rho \overrightarrow{g}\stackrel{⌢}{n}\end{array}$(2) where $p^{\mathrm{\prime }}$ is the modified pressure and ${\mu }_{eff}$ the effective viscosity. Here $p^{\mathrm{\prime }}=p+\left(2/3\right)\rho k+\left(2/3\right){\mu }_{eff}\left(\mathrm{\partial }{u}_k/\mathrm{\partial }{x}_k\right)$ and ${\mu }_{eff}=\mu +{\mu }_t$, where ${\mu }_t=C_{\mu }\rho \left(k^2/\epsilon \right)$ the turbulence viscosity, $u_i$ velocity components, $\rho$ density, $\overrightarrow{g}$ gravitational force, $\stackrel{⌢}{n}$ normal vector. To carry out this simulation, the temperature dependence was not used, and the temperature was taken as T = 300 K. In order to close the system of the equation a different turbulent models have been applied.

• Chemical reaction:

To account for the chemical reaction, substance B is used, coming out of the pipe, which reacts with substance A of the main flow, as a result substance C is formed. This chemical reaction can be written in this form: $A+B\stackrel{k}{\to }C$where k is the kinetic reaction rate, which is determined by the Arrhenius equation, which is given by $k=A{T}^{b}ex{p}^{\frac{-E_a}{RT}}$where A is the pre-exponential reaction coefficient (${10}^{15}$), T is the temperature in Kelvin, b is the temperature coefficient, R is the universal gas constant, $E_a$ is the activation energy for the reaction (${10}^8$).

• The species equations:

To calculate the species transfer with the chemical reaction, equations for the components $Y_1$ and $Y_2$ taking into account the chemical kinetics of the reaction were used. (3) $\frac{\mathrm{\partial }\rho Y_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho u_i{Y}_i\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(D_{eff}\left(\frac{\mathrm{\partial }{Y}_i}{\mathrm{\partial }{x}_j}\right)\right)+{\omega }_i$(3) where ${\omega }_i=-k_i{Y}_1{Y}_2$.

The resulting product concentration $Y_3$ is calculated according to the Dalton law: (4) $Y_3=1-Y_1-Y_2$(4) The system of equations is used the Cartesian coordinate system $x_i$ (i = 1, 2, 3), in physical space; $Y_i$ the components of species, three components of velocity $u_i,$ p – pressure; μ is the dynamic viscosity, $\vartheta =\mu /\rho$ – is the kinematic viscosity, ρ is the density, μ is the diffusion coefficient, t is the dimensionless time, $D_{eff}=\left(\rho D_{i,m}+\left({\mu }_t/S{c}_t\right)\right)$, where $S{c}_t=\left({\mu }_t/\rho D_t\right)$ – turbulent Schmidt number ($D_t$ – the turbulent diffusivity). The $S{c}_t$ by default is 0.7.

3. Turbulent models

In order to close the system Equations (1)–(4) several turbulent models have been used. It was used different turbulent model:

1. For $k-\epsilon$ turbulent model (k – kinetic turbulent energy and ϵ – turbulent dissipation rate) (Issakhov et al., 2019a, 2019b; Issakhov et al., 2020; Issakhov & Mashenkova, 2019; Issakhov & Omarova, 2020; Launder & Spalding, 1974).

2. For SST $k-\omega$ turbulent model (k- kinetic turbulent energy and ω – the specific dissipation rate) (Issakhov & Mashenkova, 2019; Issakhov & Omarova, 2020; Issakhov & Zhandaulet, 2020a, 2020b, 2020c; Menter, 1993, 1994).

4. Numerical method

The SIMPLE numerical algorithm (semi-explicit method for pressure related equations) by Patankar and Spalding (Patankar, 1980) was used to numerical solve the equations system (1)–(4). The use of this numerical algorithm is described in detail (Issakhov & Baitureyeva, 2019; Issakhov & Mashenkova, 2019; Issakhov & Zhandaulet, 2019; Issakhov et al., 2019a; Muhammad et al., 2020).

5. Test problem for the velocity field

For a more accurate implementation of the selected problem study, the values were first obtained on the basis of the available data from the papers (Ajersch et al., 1995; Keimasi & Taeibi-Rahni, 2001). The test problem area was the 3D channel with a pipe inside it. The area’s parameters as follows: the main channel height 20D, the main channel length 45D, the main channel width 3D, the height of the vertical channel 5D, and the center of the vertical channel is at a distance of 5D from the inlet. The vertical channel was square, with the diameter D = 12.7 mm (Figure 2).

Figure 2. Parameters of the calculated area.

Furthermore, the structured grid was constructed in order to obtain an accurate value: the transverse channel consists of 315 × 140 × 21 nodes, the vertical channel consists of 7 × 35 × 7 nodes, the total number of elements is 926031 (Figure 3).

Figure 3. The computational domain grid.

5.1. Boundary conditions

In paper (Keimasi and Taeibi-Rahni) considered different ratios of jet velocity to transverse flow R (0.5, 1.0, and 1.5). In this study, the ratio R = 0.5 was used; at this ratio, the jet velocity was 5.5 m/s, the cross flow velocity was 11.0 m/s. The air was selected as the medium. The diameter of vertical channel (D) was 12.7 mm. The Reynolds number was determined like that: $Re=\left(\right(\rho V_{jet}D\left)/\mu \right)=4700$. To represent the inlet velocity profile of the transverse flow the 1/7 power law of the wind in the boundary layer was used. $\frac{u}{u_r}={\left(\frac{z}{z_r}\right)}^{\alpha }$At a height more than 2D, the U-velocity was set to a stationary and determined like 11.0 m/s. The V-velocity for the vertical channel was set to a stationary and determined like 5.5 m/s (Figure 4). Here u is the wind speed at the height z, and $u_r$ is the known wind speed at the reference height $z_r$. $\alpha$ _– empirically derived coefficient, which changes depending on the atmosphere stability. Here $\alpha =1/7$ for conditions of neutral stability (Figures 4 and 5).

Figure 4. Boundary conditions of the reactor (OXY plane).

Figure 5. Boundary conditions of the reactor (OYZ plane).

These velocity values were set to repeat the measurement done by (Ajersch et al., 1995). At the entrance boundary of the horizontal flow, conditions were set on the components of the (k-ω)/(k-ϵ) models.

During calculations using the k-ϵ turbulent model, the following initial condition were specified. For turbulent kinetic energy profile was occurred that relies on the initial velocity in order to receive closest agreement with the experiment by (Ajersch et al., 1995). $k=\left\{\begin{array}{ll}\left(0.3+0.3{\left(\frac{0.008-y}{0.0254}\right)}^{\frac{1}{7}}\right)V_{hor}& y<0.008\\ 0.3\cdot V_{hor}& y=0.008\\ 0.3-0.3\cdot {\left(\frac{y-0.008}{0.0254}\right)}^{\frac{1}{7}}{V}_{hor}& \begin{array}{c}0.008<y\\ <0.0254\end{array}\\ 0.0012\cdot V_{hor}& y>0.0254\end{array}\right.$The boundary conditions for dissipate turbulent energy was used like in paper (Versteeg & Malalasekera, 1996): ${\epsilon }_{in}=C_{\mu }^{\frac{3}{4}}\left(k^{\frac{3}{2}}/L_m\right)$where C_μ = 0.09 is an empirical constant.

For the boundary conditions operating as a ‘wall’, the standard conditions were selected, which are offered by the ANSYS Fluent. The periodic boundary conditions for the horizontal channel on the Z axis were established.

The following initial conditions were used for simulation the k-ω turbulent model. For the dissipation rate of turbulent energy, a formula was set that was proposed by (Menter, 1994): $\omega =(1\to 10)\left(V_{hor}/L\right),$where L is the reactor length. The horizontal acceleration at the outlet boundary was determined as 0.

5.2. Comparison of numerical results

The received computational values were compared with the computational values (Keimasi & Taeibi-Rahni, 2001) and values from the measurement studies (Ajersch et al., 1995). In Figure 6, the U-velocity profiles for various turbulent models were matched, the received computational values agree well with the measured values for the transverse section x/D = 0, z/D = 0 (Ajersch et al., 1995), but for the transverse section x/D = 3, z/D = 0 as it can be seen there are several differences in a data. From Figure 7, it could be noted that the V-velocity profiles for various turbulent models characterize the flow behavior well, but, however, deviations from empirical values appear in all turbulence models (Ajersch et al., 1995), except for the SST k-w turbulent model. Among the received computational data using various turbulent models (Figure 7), it is worth highlighting the SST k- w turbulent model, the computational data of this model are closest to the measured values for all cross sections (Ajersch et al., 1995). From Figure 8 it can be noted that the received computational data by the SST k-w model for the V-velocity agree well with the measured values (Ajersch et al., 1995).

Figure 6. Comparison of the profile of U-velocity for various transverse sections.

Figure 7. Comparison of the profile of V-velocity for various transverse sections.

Figure 8. Comparison of the profile of W-velocity for various transverse sections.

During verification of turbulent models (Figures 6–8), the SST k-w model turned out to be the closest and optimal model. In the future, the turbulent SST k-w model was applied for the problem in real sizes, since the LES turbulent model captures low-frequency parameter changes, unlike RANS, when constant average values are obtained, and large computational power is expended.

6. Test problem for the velocity field and concentration components

This section discusses the following test problem. A detailed description of this test problem and experiment is given in (Crabb et al., 1981; Majander & Siikonen, 2006). The test problem area was 3D channel with a pipe inside it. The pipe diameter (jet width) was D = 12.7 mm, which was used as a unit of length. The height of the transverse channel H = 7D, the height of the pipe 2D, the length of the transverse channel L = 13D, the width of the channel W = 8D, and the center of the vertical channel is at a distance $L_1$ = 3D from the entrance (Figure 9).

Figure 9. Map of the calculated area.

Air was chosen as the computational material for the transverse flow, and water vapor for the jet. For the numerical simulation of this problem, the following parameters were taken: dynamic air viscosity was $\mu =1.7894\times {10}^{-5}$ kg·${\text{m}}^2$/sec, density ρ = 1.225 kg/${\text{m}}^3$; the dynamic viscosity of water vapor is $\mu =1.34\times {10}^{-5}$ kg·${\text{m}}^2$/sec, density ρ = 0.5542 kg/${\text{m}}^3$. The initial temperature was selected as 300 K. The Reynolds number was defined as $R{e}_{jet}=\rho V_{jet}D/\mu =46700$ (Crabb et al., 1981; Majander & Siikonen, 2006).

6.1. Construction of the computational grid

An unstructured grid was built for the simulation and was thickened near the wall, around the exit of the pipe and along the ejection path (the size of the computational grid is up to 17 cm). The total number of three-dimensional cells is more than 1054656 elements (Figures 10 and 11).

Figure 10 Computational grid of the computational domain.

Figure 11. (a) horizontal section and (b) vertical section in the central plane, z/D = 0.

6.2. Boundary conditions

To describe the initial velocity profile of the jet, the profile of the substance in the form of an exponent was used. $v_y=v\ast (1-e^{-1510\zeta }),\zeta =min(x;L_x-x)$The initial velocity profile of the transverse flow is given as a constant. The ratio of jet velocity to cross flow velocity was determined like $R=V_{jet}/V{}_{crossflow}$. For this test problem, the ratio is taken as R = 2.3, so the jet velocity was 53.7 m/s, the cross-flow velocity was 23.35 m/s. The remaining boundary conditions are displayed in Figure 12. So the boundary conditions for the jet are taken as the concentration of water vapor is 1, and the concentration of the transverse air flow is taken as 0. And for the cross-flow, the boundary conditions are taken as the concentration of water vapor is taken as 1.

Figure 12. Boundary conditions of the reactor (OYZ plane).

6.3. Comparison of numerical results

Figure 13 compares the mean horizontal velocity profiles with experimental data (Crabb et al. (1981)) and numerical data that were solved using the LES turbulent model (Majander & Siikonen,  2006) with the obtained numerical data in this work for different sections x/D = {−1, 0, 0.5, 0.75, 1, 1.5, 4, 6, 8} and z/D = 0. In Figure 14, the mean concentration profiles were compared with measured values (Crabb et al., 1981) and numerical values that were solved using the LES turbulent model (Majander & Siikonen, 2006) with obtained numerical data in this work for different sections x/D = {1, 2, 4, 6, 8, 10} and z/D = 0. From all the obtained computational values, it can be noted that the received values using the SST k-w turbulent model for the U-velocity agree adequately with the empirical data (Crabb et al., 1981).

Figure 13. Mean flow rates in the central plane, z/D = 0: (Crabb et al., 1981) (o), Majander & Siikonen (2006) (-) LES (-).

Figure 14. The mean proportion of the mixture in the central plane, z/D = 0: (Crabb et al., 1981) (o), Majander & Siikonen (2006) (-) LES (-).

Thus, the obtained numerical results for the concentration component showed much better results than the obtained results using the LES turbulent model. Thus, the obtained values from the computation using the turbulent SST k-w model showed good agreement with the measured values and computational data for both the velocity components and the concentration components.

7. Test problem for a chemical reaction

To implement the test problem for a chemical reaction, a 2D model of a microreactor was constructed based on the data provided by Falconi et al. (2007). A similar work was done by foreign researchers (Schonauer & Adolph, 2005). The microreactor consisted of two channels. The diameter of the pipe (jet width) was D = 1 m, which was used as a unit of length. The height of the transverse channel 4D, the height of the pipe 2D, the length of the transverse channel 16D. The vertical channel was located precisely at the center of the vertical channel along the X-axis and was shifted by a range equal to 2.5D from the origin of the horizontal channel (Figure 15). Here the jet inlets from the pipe, perpendicular to the main flow in the channel. The flow conditions for the jet and the transverse flow were set by two various predetermined velocity profiles. Chemical component B inlets perpendicular to the flow (in 2), component A through the side channel (in 1). Figure 13 shows the configuration of the microreactor under study. Here, it was considered an incompressible fluid with $Re=\left(\right(\rho V_{\mathrm{\infty }}D\left)/\mu \right)=25$, where the chemical components A and B react and form a new substance C. ($A+B\stackrel{k}{\to }C$). The ratio of speeds was expressed through a coefficient $R=U_{jet}/U_{crossflow}=1.5$.

Figure 15. Scheme of the calculated area.

7.1. Initial and boundary conditions

Boundary conditions were determined like: the outlet – ‘Pressure-outlet’, the inlet 1 and inlet 2 – ‘Velocity-inlet’, the walls – ‘Wall’.

The velocity profile in the main channel was described as follows: $u_x=u_{wind}(1-e^{-5\zeta })$, where $\zeta =\text{min}(z,L_z-z)$. The remaining parameters were set constant: w = 0, $Y_A$ = 1, $Y_B$ = 0.

In the inlet pipe: u = 0, $w=2R{u}_{wind}(1-4r^2)$, where r = x, $Y_A$ = 0, $Y_B$ = 1. The remaining boundary conditions are specified as indicated in Table 1 and in Figure 16.

Figure 16. Numbering the outer boundaries of the solution domain.

Table 1. Sequence of variables and equations at the boundaries of the $◯1$–$◯6$ domain.

Var.	$◯1$	$◯2$	$◯3$	$◯4$	$◯5$	$◯6$
u	$u_{wind}(1-e^{-5\zeta })$	u = 0	${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }x}}=0$	u = 0	u = 0	u = 0
W	w = 0	Equation (1)	${\displaystyle \frac{\mathrm{\partial }w}{\mathrm{\partial }x}}=0$	w = 0	w = 0	w = 0
P	Equation (2)	p = $p_{atm}$	Equation (2)	Equation (2)	Equation (2)	Equation (2)
$Y_A$	$Y_A=0$	$Y_A=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_A}{\mathrm{\partial }x}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_A}{\mathrm{\partial }y}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_A}{\mathrm{\partial }x}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_A}{\mathrm{\partial }y}}=0$
$Y_B$	$Y_B=1$	$Y_B=1$	${\displaystyle \frac{\mathrm{\partial }{Y}_B}{\mathrm{\partial }x}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_B}{\mathrm{\partial }y}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_B}{\mathrm{\partial }x}}=0$	${\displaystyle \frac{\mathrm{\partial }{Y}_B}{\mathrm{\partial }y}}=0$

Here $u_{wind}$ varies depending on the selected substance. For substances A and B were indicated as oxygen O₂. For the numerical simulation of this problem, the following parameters were taken: the dynamic viscosity of oxygen was $\mu =1.919\times {10}^{-5}$ kg ${\text{m}}^2$/sec, density ρ = 1.299874 kg/${\text{m}}^3$, speed $u_{wind}=\text{0}\text{.000369074233}\text{m/s}$, hydraulic diameter D = 1 m, the diffusion coefficient was set, as 1.4762674052 m/${\text{s}}^2$. To satisfy the $Sc=1$, the diffusion coefficient was determined like $\text{0}\text{.677}{\text{m}}^{\text{2}}/s$. The rate of chemical reaction was taken as 1. In ANSYS Fluent, all calculations were made in actual sizes, so in this case the actual parameters were indicated.

7.2. A computational grid

To obtain an accurate result, a structured computational grid was constructed. The main part of the channel consists of 641 × 141 computational points, the dimensions of the computational grid for the pipe are 41 × 81 and total the number of nodes is 638,886. The appreciate parameters for computational mesh generation is based on the computational data of Falconi et al. (2007) (Figure 17).

Figure 17. Calculate computational grid.

7.3. Comparison of computational values

The computational values and a comparative analysis for the mass fraction and the velocity field are performed below. Values were matched with numerical values by Schonauer and Adolph (2005) and illustrated a good match. Figure 18 shows a comparison of the velocity streamlines with the computational values received by calculations of Schonauer and Adolph (2005). Figure 19 illustrates a comparative analysis of the velocity contour with the data from computations obtained by calculations of Schonauer and Adolph (2005) for horizontal velocity u and vertical velocity v. Figure 20 presents a comparative analysis of the of the concentrations distribution results obtained of this work and the results obtained by Falconi et al. (2007); Schonauer and Adolph (2005). As can be seen from the obtained data, the contours for the horizontal velocity u and the vertical velocity v showed good values in comparison with the calculations of Schonauer and Adolph (2005). The vortex length behind the pipe for both cases was l = 6.2D. It should also be noted that the distribution of the concentration component and the new concentration component also showed good agreement with numerical calculations (Schonauer & Adolph, 2005) (Figure 20).

Figure 18. The velocity streamlines comparison: the upper result is the calculations of Schonauer and Adolph (2005), the lower result is the obtained results in this work.

Figure 19. Comparative velocity contour analysis: left plots – results from Schonauer and Adolph (2005), the right-hand plots are the values obtained in the course of this work: (a) horizontal speed u (b) vertical speed v.

Figure 20. Comparative analysis of the results of the substance distribution: the left graphs – the values of Schonauer and Adolph (2005), the right-hand plots are the obtained results in this work: (a) concentration of substance A, (b) concentration of substance B, (c) concentration of substance C.

Thus, the obtained values from the computation using the turbulent SST k-w model showed good agreement with the measured values and computational data for both the velocity components, the concentration components and the concentration of product from the chemical reaction.

8. Real size problem

Next, the distribution of pollutants from Ekibastuz SDPP-1 in real sizes was considered (Figure 1). The computational domain is a 3D box with two chimneys inside it (Figure 21). At this Ekibastuz SDPP-1, pollution is produced from 2 chimneys, which have various altitudes of 300.0 and 330.0 m, respectively. The range between the two chimneys is 200.0 m (Figure 22). The diameters of the chimney holes are about 10.0 m. The NO, NO₂, CO, etc. are emitted in large quantities as pollutants from Ekibastuz SDPP-1. These pollutants in the atmosphere react with oxygen and water vapor. This paper discusses 3 major chemical reactions.

1. $2\text{CO}+{\text{O}}_2\to 2\text{C}{\text{O}}_2$

2. $2\text{NO}+{\text{O}}_2\to 2\text{N}{\text{O}}_2$

3. $4\text{N}{\text{O}}_2+2{\text{H}}_2\text{O}+{\text{O}}_2\to 4\text{HN}{\text{O}}_3$

Figure 21. Configuration of computing area for thermal power plant.

Figure 22. Parameters and dimensions of the calculated area.

(1) The chemical reaction of carbon monoxide (CO) with oxygen produces carbon dioxide (CO₂) (carbon dioxide). Because of its toxicity, carbon monoxide is very dangerous for the human body. This danger is exacerbated by the fact that it is odorless and poisoning can occur unnoticed. Carbon monoxide interacts with blood hemoglobin (carboxyhemoglobin is formed), which is 200 times more active than oxygen and reduces the ability of blood to be its carrier. Therefore, even with low concentrations of CO in the air, it has a detrimental effect on health. When a chemical reaction with oxygen forms carbon dioxide. Carbon dioxide easily transmits radiation in the ultraviolet and visible parts of the spectrum, which comes to the earth from the sun and heats it. At the same time, it absorbs infrared radiation emitted by the earth and is one of the greenhouse gases, as a result of which it takes part in the process of global warming.

Besides the greenhouse properties of carbon dioxide, the fact that it is heavier than air is significant. Since the average relative molar mass of air is 28.98 g/mol and the molar mass of CO₂ is 44.01 g/mol, an increase in the proportion of carbon dioxide leads to an increase in air density and, accordingly, to a change in its pressure profile depending on height. Due to the physical nature of the greenhouse effect, such a change in the properties of the atmosphere leads to an increase in the average temperature on the surface (Stamler & Gow, 1998). Since, with an increase in the proportion of this gas in the atmosphere, its large molar mass leads to an increase in density and pressure, at the same temperature an increase in CO₂ concentration leads to an increase in the air capacity and to an increase in the greenhouse effect due to the large amount of water in the atmosphere. Increasing the proportion of water in the air to achieve the same level of relative humidity – due to the low molar mass of water (18 g/mol) – reduces the density of air, which compensates for the increase in density caused by the presence of elevated levels of carbon dioxide in the atmosphere.

(2) As a result of the chemical reaction of nitric oxide (NO) with oxygen, nitrogen dioxide (NO₂) is formed. Nitrogen oxide (II) is odorless, but when inhaled, it can bind to hemoglobin, like carbon monoxide gas converting it into a form (methemoglobin) that cannot transport oxygen (Stamler & Gow, 1998). Combined with oxygen in the air, it forms nitrogen dioxide. The NO₂ impurities paint industrial smokes brown. NO₂ irritates the lungs and can lead to serious health consequences. A concentration of NO₂ above 200 ppm is considered lethal, but already at a concentration above 60 ppm, discomfort and burning in the lungs can occur. NOx emissions are considered one of the main reasons for photochemical smog formation. Elevated concentrations of NOx have a detrimental effect on human health, therefore, standards have been adopted in various countries that limit the maximum allowable concentrations of NOx in the exhaust of power plants boilers, gas turbines, cars, airplanes and other devices. Improving combustion technologies is largely aimed at reducing NOx emissions while improving the energy efficiency of devices.

(3) As a result of the chemical reaction of nitrogen dioxide (NO₂) with oxygen and water vapor, nitric acid (HNO₃) is formed. As mentioned earlier, the resulting products are more harmful than the original. Nitric acid (HNO₃), together with sulfur oxides (SOx), is the cause of the formation of acid rain. Nitric acid is poisonous. According to the degree of impact on the body, it belongs to substances of the 3rd class of danger. It dissolves in fat and can penetrate into the capillaries of the lungs, where it causes inflammation and asthmatic processes, irritates the respiratory tract (Table 2).

Table 2. Parmameters of species.

	Parameters
Species	Density (kg/m^3)	Viscosity (w/m-K)	Cp (Specific Heat) (j/kg-K)	Thermal Conductivity (w/m-K)	Molecular Weight (kg/kmol)	Standard State Entropy (j/kgmol-K)
O2	1.2999	1.919 × 10^−5	piecewise-polynomial	0.0246	31.9988	205026.9
NO	1	1.72 × 10^−5	piecewise-polynomial	0.0454	30.0061	210640
NO2	1	1.72 × 10^−5	piecewise-polynomial	0.0454	46.0055	239909.6
H2O	0.5542	1.34 × 10^−5	piecewise-polynomial	0.0261	18.01534	188696.4
HNO3	1504	0.011	piecewise-polynomial	0.0454	63.01287	266389.2
Air	1.225	1.7894 × 10^−5	1006.43	0.0242	28.966	194336
CO	1.1233	1.75 × 10^−5	piecewise-polynomial	0.025	28.01055	197531.6
CO2	1.7878	1.37 × 10^−5	piecewise-polynomial	0.0145	44.00995	213720.2

8.1. Grid construction

The computational grid was condensed in the area of the pollution trajectory, in which the area where the buildings are located (the size of the computational cell is up to 20 m) and then the concentration of the building (the size of the computational cell is up to 5 m) and the chimneys (the size of the computational cell is up to 2 m) (Figures 23–25). The computational grid consists of 218915 nodes and 1244098 elements (Table 3).

Figure 23. (a) computational grid and (b) bottom view of the grid.

Figure 24. (a) section along the axis OXY and (b) section along the axis OYZ.

Figure 25. View of the chimneys and buildings.

Table 3. Parmameters of geometry.

Geometry parameters (m)	First chimney coordinates (300 m)	Second chimney coordinates (330 m)
X = 10,000	X = 1000	X = 1000
Y = 1800	Y = 0	Y = 0
Z = 4000	Z = 2100	Z = 1900

8.2. Boundary conditions

The boundary conditions were determined like: for the side walls and the top wall – the ‘Symmetry’, for the wind inlet and the chimney hole – the ‘Velocity inlet’, for the outlet – the ‘Pressure Outlet’ and for chimney walls, the ground – the ‘Wall’ (Figure 26).

Figure 26. Boundary conditions of the computational domain.

To represent the inlet velocity profile of the transverse flow the 1/7 power law of the wind in the boundary layer was used, as in Keimasi and Taeibi-Rahni (2001). The jet velocity was 1.5 m/s, the velocity of the transverse flow was 3 m/s. So the ratio of the jet velocity to the velocity of the transverse flow is R = 0.5. The ambient temperature is taken as 300 K. So the boundary conditions for the chimneys are taken as the concentration of different species like 1, and the concentration for different species of the transverse flow is taken as 0. And for all other domain sides are taken as no-flux boundary.

8.3. Numerical results

The numerical modeling was done in the ANSYS Fluent, where all modeling were performed in actual scales. For the computational solution SIMPLE approach was chosen.

Figure 27 shows the concentration distribution of pollution, which is visualized using Volume Rendering. So Figure 27(a-d) displays the distribution of CO and CO₂ concentrations and for different angles of view. From Figure 27(e-h) one can see the spread of concentration for different angles of view. Figure 27(i-l) shows the distribution of concentration NO₂, HNO₃ for different angles of view. So from the obtained results it can be noted that the emitted concentrations and products from a chemical reaction extend very far. But it can also be noted that heavy elements obtained by a chemical reaction quickly precipitate near Ekibastuz SDPP-1. It can also be noted that with increasing distance from the source, a diffusion effect is observed for all kinds of concentration. Figures 28–30 show a comparison of the mass fraction profiles CO, CO₂, NO, NO₂, H₂O, NO₂, HNO₃ for different points: 1500, 2000, 2500, 3000, 10,000 and z = 2000.

Figure 27. Visualization of the spread of concentrations: (a-d) CO and CO₂, (e-h) NO and NO₂, (i-l) NO₂ and HNO₃.

Figure 28. Comparison of profiles of the CO, CO₂ mass fractions in points: x = 1500, 2000, 2500, 3000, 10,000 and z = 2000.

Figure 29. Comparison of profiles of the NO, NO₂ of mass fractions in points: 1500, 2000, 2500, 3000, 10,000 and z = 2000.

Figures 31–33 represent the comparison of the mass fractions of CO, NO, H₂O, NO₂ from two chimneys (300 and 330 m) at different distances at a wind speed of 3 m/s: 1000, 1500, 2500, 3000 and 10,000 m from the origin. As can be seen in the graphs, the mass fraction for a chimney with an altitude of 330 m is higher compared to a chimney of 300 m. Accordingly it can be concluded that the altitude of the chimneys essentially affects the contaminants spreading. Construction of higher chimney is more suitable for environmental safety.

Figure 30. Comparison of profiles of the H₂O, NO₂, HNO₃ mass fraction in points: 1500, 2000, 2500, 3000, 10,000 and z = 2000.

Figure 31. Comparison of profiles of the CO mass fraction at specified points for two various heights of chimneys (300.0 and 330.0 m).

Figure 32. Comparison of profiles of the NO mass fraction at specified points for two various heights of chimneys (300.0 and 330.0 m).

Figure 33. Comparison of profiles of the H₂O mass fraction at the specified points for two various heights of chimneys (300.0 and 330.0 m).

Figures 34–37 show the mass fractions of the products of a chemical reaction CO₂, NO₂, HNO₃ of two chimneys (300 and 330 m) compared at different distances at a wind speed of 3 m/s: 1000, 1500, 2500, 3000 and 10,000 m from the origin. From Figures 34–36 it can be noted that the pollutants do not immediately mix and enter into a chemical reaction, and the reaction process takes place through a certain distance from the chimneys. So active mixing starts from 500 m from the source. However, it should be noted that the concentration of CO and CO₂ practically does not have or has a very small value of concentration on the surface of the earth. This phenomenon is due to the fact that the density of CO is less than the density of air. And the concentration of CO₂, despite the fact that it matters more than the density of the air, could be explained by the fact that the impurity value CO₂ has a very small value. And also an important factor is that the formation of concentration CO₂ occurs around the distribution of CO. So from the data it can be concluded that the concentration of CO and CO₂ almost not deposited on the surface of the earth and goes down to about a height of 180 m.

Figure 34. Comparison of profiles of the CO₂ mass fraction at specified points for two various heights of chimneys (300.0 and 330.0 m).

Figure 35. Comparison of profiles of the NO₂ mass fraction at specified points for two various heights of chimneys (300.0 and 330.0 m).

Figure 36. Comparison of profiles of the NO₂ mass fraction at the specified points for two various heights of chimneys (300.0 and 330.0 m).

Figure 37. Comparison of profiles of the HNO₃ mass fraction at the specified points for two various heights of chimneys (300.0 and 330.0 m).

However, for the concentration of NO and NO₂ it can be noted that these concentrations descend a little below 100 m despite the fact that NO and NO₂ are lighter than air and than concentrations of CO and CO₂, this is due to the fact that these effects can be observed only in remote sections, since there is observed both the transfer process and the diffusion process. This process is due to turbulent vortex dissipation, which leads to the fact that the concentration of NO and NO₂ distributed in all directions evenly.

However, for concentration HNO₃, it can be seen the opposite situation, that most of the concentration settles on the ground surface and not too far from the Ekibastuz SDPP-1 chimney itself, and only a small part of the concentration HNO₃ remains on the atmospheric boundary layer. This rapid deposition of the most part HNO₃ could be explained that the density of the concentration HNO₃ is many times greater than the density of air. So from the obtained results it can be concluded that nitric acid (HNO₃) pollutes not only near the surface of the earth, but also the surface of the atmospheric layer. In the first case, the nearby territories of Ekibastuz SDPP-1 are polluted, and in the second case, a small part of the concentration of nitric acid (HNO₃) not only pollutes the air, but also lagging behind the small concentration in the atmosphere together with sulfur oxides (SOx) lead to cause and risk in the future the formation of acid rain.

In all the obtained pollution profiles, both the diffusion process and transfer process are seen. Because of the diffusion reason by turbulent vortex dissipation, the gas concentration sets to be higher at ground level by 20.0 m (where diffusion has smoother values). From the computational values, it could be noticed that the concentration rate that is releasing from two chimneys performs the significant role in the impurity distribution. Accordingly, having determined ground-based measurement points of concentrations CO₂, NO₂, HNO₃ in the regions, will be of significant interest in next study.

From the obtained computational values, it could also be noted that behind the chimneys as such, flow vortex regions of different sizes are not formed. Since with this movement of air, inertial forces predominantly prevail, which leads to the fact that the ejected concentration is not concentrated in the areas around the chimneys. And further from the chimneys, a clear dispersion of the allowed concentration is seen. So, in general, flow vortices of various sizes are visible in the concentration plume. Due to these processes, it can be seen from the results that the concentration of emissions is dissipated and expand downstream.

9. Conclusion

This paper presents the results obtained by numerical simulation using the software ANSYS FLUENT for the spread of pollutants formed during the combustion of fuel in an thermal power plant, and their chemical reaction in the atmosphere. The numerical method were tested using three test problems. Using several test problems for determination of the optimal turbulent model (k−ϵ and k–ω). Moreover, the analysis of the received computational values with the data of modeling other authors and measured data was done. As a result, the k – ω turbulent model was defined as the appreciate for the transfer to actual size problem.

Ekibastuz SDPP-1 was used to simulate an actual problem. A remarkable feature of this TPP is that the difference between the chimneys allows us to study the affect of the source height on the dispersion of pollution. The aim of the study was to study the dynamics of the emissions distribution and their chemical reaction with oxygen in the air. Using the example of a real thermal power plant, dispersions of CO, NO, H₂O, NO₂ and products CO₂, NO₂, HNO₃, were simulated during a chemical reaction with oxygen. From the values of the computation it can be concluded that the altitude of the chimneys essentially affects the gases spreading. In the analysis, it was found that it is necessary to simulate the process of locating a source that spreads harmful impurities from a TPP to the atmospheric boundary layer in such a way as to minimize air pollution and surrounding areas. The construction of higher chimneys for a TPP is most suitable for environmental safety.

It has to be noticed that in this work there are some limitations. The major limitation is the size of the computational grid. Computer resources are limited in the size of the computational grid, while a very huge computational grid is necessary to obtain an accurate simulation result. In addition, a further increase in computer resources will be followed by the requirements of researchers to increase the resolution of the grid and the extra physical parameters inclusion that will lead to actual size problem. The next limitation of the study is the complicacy of the analysis of computational studies on the distribution of chemical compounds in the atmosphere from the TPPs activities (size, shape, environmental factors). As well as the third limitation of this work is an accurate description of chemical reaction processes, since under real conditions the number of chemical reactions exceeds more than hundred reactions, and in modeling, a set of simplified chemical mechanisms is usually used to save computer resources. This work is useful for whose interested in the spreading of chemical compounds in the atmosphere from the TPPs activities.

Acknowledgements

This work is supported by grant from the Ministry of education and science of the Republic of Kazakhstan.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by grant from the Ministry of Education and Science of the Republic of Kazakhstan.