Chimney Height, a Determining Factor in the Dispersion of Pollutants and Their Concentration

Neyva Maria Lopes Romeiro¹, Eliandro Rodrigues Cirilo¹, Paulo Laerte Natti¹, Letícia Mayumi Doy Okamoto², Tainá Julianotti³

¹Department of Mathematics, State University of Londrina and Member of the Simulation and Numerical Analysis Laboratory—LabSan, Paraná, Brazil.
²Graduated in Computer Science, State University of Londrina, Paraná, Brazil.
³Master’s Student, Applied and Computational Mathematics—PGMAC, State University of Londrina, Paraná, Brazil.

DOI: 10.4236/wjet.2021.91013                

Abstract

This work proposes a way of modelling two-dimensional complex meshes using elliptic equations, which the computational grid coincides with the problem geometry, making computational processing more suitable. A multiblock technique was used in order to achieve a better representation of the problem domain. In this way, numerical simulations of the movement and dispersion of pollutant emissions in the atmosphere are presented in the generated domains, using the Navier-Stokes pollutant transport equations. The curvilinear coordinates and the finite difference method are used for the discretization. The model was verified in two tests. In the first test, three cases were proposed, with geometries containing a chimney followed by an obstacle, using different chimney heights, and the obstacle height was fixed. The test aims to verify the vortices appearance, in the blocks, to obtain agreement with as presented in the literature. In the second test, the geometry is described by a chimney and an obstacle that represents one of the mountains in the valley. The performed tests made possible to verify that the height of the chimney can be considered a determining factor to describe the dispersion of pollutants, as well as their concentration in the proximity of industrial areas.

Keywords

Atmosphere, Valley, Multiblock, Navier-Stokes, Transport Equation

1. Introduction

The emission of pollutants into the atmosphere interferes with air chemical composition and can cause both environmental and human health damage [1] [2] [3]. Several sources are responsible for emitting such pollutants and can be classified as natural, such as natural fires and volcanic emissions, or anthropogenic like burning fossil fuels in factories [4] [5] [6]. Understanding how emissions from these sources occur, especially anthropogenic ones, we try to minimize these damages [7] [8] [9]. Therefore, many researchers have been developing their studies using mathematical modeling [7] [10] [11] [12]. Still, several methods have been used to evaluate pollutant dispersion, considering geometries with point sources like chimney or obstacles [13] [14] [15].

In this context, this work will be conducted in simulations that describe the movement and dispersion of pollutants using mathematical modeling. The numerical solution is designated only for a passive scalar transport, the concentration is solved for a given stationary velocity field, and it can deal with only the constant diffusion coefficient, so that it cannot take into account the turbulent diffusion caused by turbulent flows [15].

The computational mesh that represents the physical domain in the tests was modeled, followed by the numerical model described by the Navier-Stokes equations [16] [17] [18] [19] [20] and the pollutant transport equation [21] [22] [23].

The solutions precision obtained by the simulations of fluid dispersion relies on how well computational mesh represents the studied geometry. For problems with complex geometries, it is convenient to use the curvilinear coordinate system instead of cartesian coordinates, due to the fact that cartesian coordinates lead to a poor boundary fit, since the physical domain doesn’t always coincide with the mesh domain [18] [24] [25] [26]. Another reason that leads to the use of the curvilinear coordinate system, in the discretization of the computational mesh, refers to less difficulty in programming computational codes to solve complex problems and the facility in developing generic methodologies [27].

Despite the advantages of using the curvilinear coordinates for the construction of meshes, it is still possible to obtain low-quality elements. A way to circumvent the quality of the mesh elements, without necessarily improving the refinement, because it would increase the number of nodes in the mesh and, consequently, the computational costs with memory and execution time, refers to the use of the multiblock technique [18] [26] [27] [28] [29]. Applying this technique, the domain is divided into parts, called block or sub-grids, and the equations are solved in each block.

With independent block separation, it is possible to increase mesh refinement at important areas or decrease it in less relevant areas, allowing the generation of structured grids to keep simulations costs balanced. Illustrating a comparison between meshes generated by curvilinear coordinates, consider a single block, Figure 1(a), and multiblocks (5 blocks), Figure 1(b).

For the geometry, Figure 1, the meshes were generated with the same number of partitions. Due to the rectangular form of the geometry and the block division structure, the mesh in curvilinear coordinates resulted in square elements.

(a) (b)

Figure 1. Comparison between meshes generated with a single block and multiblock techniques, both with the same number of partitions. (a) Mesh generated with single block (b) Mesh generated with multiblock technique.

Therefore, through the multi-block technique, the quality of the mesh elements can be improved.

Hence, this work proposes a way of modeling two-dimensional meshes using elliptic equations, generating meshes in the curvilinear coordinate system coinciding with the problem geometry. Still, for a better representation of the problem domain, was used the multiblock technique [26] [27] [28] [29]. The velocity field was obtained by simulating the Navier-Stokes equations and then the transport equation is solved, modeling the pollutant dispersion in the mesh. The discretizations are obtained by applying the finite difference method and using curvilinear coordinates. In the Navier-Stokes equations, it is applied the first order upwind scheme in the convective term [17] [26], central differences in the second order terms and progressive differences for the temporal term [18]. In the pollutant transport equation, it uses central differences in terms of first and second order and progressive differences for the temporal term.

2. Methodology

To generate the meshes, the blocks were represented using the notation involving the cardinal points, where the W side represents the west direction, E, east, N, north and S, south direction. Then, to generate the blocks, parameter with the number of elements in ξ and η directions and their coordinates points to each side were passed.

The boundary conditions are important for a fluid dispersion numerical solution, since they define the simulation limits and help to guarantee the method stability and convergence [18]. Besides, correctly define the boundaries of the block is extremely important. The implemented algorithm recognizes the boundary conditions type of the four sides of each mesh block.

To solve the Navier-Stokes equations were used these boundary conditions:

• CNEI: no-slip and impermeability condition of the fluid, indicates the wall fluid’s velocity is null and there isn’t flow through the wall as well;

• CIPR: prescript input condition, indicates that at the domain border it has non-zero velocity;

• CECO: continuous fluid flow condition, indicates the domain fluid output;

• ADJA: adjacency blocks condition, indicates the data transferring between adjacent blocks.

The boundary condition to solve the transport equation can be:

• CIPRCONC: prescript concentrations condition, used at the domain border, indicates the concentration that enters into the domain across the border;

• CCONCCONC: continuous concentration condition, indicates the domain concentration output;

• ADJACONC: used to indicate adjacency between blocks.

A concern while using multiblocks technique is to maintain simulation consistency, for this purpose, a treatment occurs in the boundaries that indicate adjacency between blocks [18] [26] [27] [28] [29].

The two bands of previous/posterior elements of neighboring blocks serve as transfer area, with no solution calculation in them. To guarantee simulation continuity, it is required that the adjacent boundaries points have the same coordinate in both neighboring blocks. Figure 2 shows an example of adjacency boundary, illustrating how was accomplished the data transfer between two blocks.

3. Governing Equations

Assuming by hypothesis a laminar flow, newtonian, isothermal and incompressible, in two-dimensional, without the existence of the source terms, then the Navier-Stokes equations described in the curvilinear coordinate system are written by:

$\begin{array}{l}\underset{\text{temporal term}}{\underset{︸}{\frac{\partial }{\partial \tau }\left(\frac{\rho u}{J}\right)}}+\underset{\text{convective term}}{\underset{︸}{\frac{\partial }{\partial \xi }\left(\rho Uu\right)+\frac{\partial }{\partial \eta }\left(\rho Vu\right)}}\\ =\underset{\text{diffusive term}}{\underset{︸}{\mu \left[\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial u}{\partial \xi }-\beta \frac{\partial u}{\partial \eta }\right)\right)+\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial u}{\partial \eta }-\beta \frac{\partial u}{\partial \xi }\right)\right)\right]}}+\underset{\text{pressure term}}{\underset{︸}{\left[\frac{\partial p}{\partial \eta }\frac{\partial y}{\partial \xi }-\frac{\partial p}{\partial \xi }\frac{\partial y}{\partial \eta }\right]}}\end{array}$ (1)

$\begin{array}{l}\underset{\text{temporal term}}{\underset{︸}{\frac{\partial }{\partial \tau }\left(\frac{\rho v}{J}\right)}}+\underset{\text{convective term}}{\underset{︸}{\frac{\partial }{\partial \xi }\left(\rho Uv\right)+\frac{\partial }{\partial \eta }\left(\rho Vv\right)}}\\ =\underset{\text{diffusive term}}{\underset{︸}{\mu \left[\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial v}{\partial \xi }-\beta \frac{\partial v}{\partial \eta }\right)\right)+\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial v}{\partial \eta }-\beta \frac{\partial v}{\partial \xi }\right)\right)\right]}}+\underset{\text{pressure term}}{\underset{︸}{\left[\frac{\partial p}{\partial \eta }\frac{\partial x}{\partial \xi }-\frac{\partial p}{\partial \xi }\frac{\partial x}{\partial \eta }\right]}}\end{array}$ (2)

Figure 2. Adjacency treatment between blocks.

where the dynamic viscosity μ and the specific mass ρ are constant. For further details of how to describe these equations in this coordinate system see Maliska [18] and Cirilo et al. [26].

The terms U and V, in Equations (1) and (2), are named contravariant components of the velocity vector, these are normal to ξ and η lines respectively, and

defined as $U=\frac{1}{J}\left(u\frac{\partial \xi }{\partial x}+v\frac{\partial \xi }{\partial y}\right)$ and $V=\frac{1}{J}\left(u\frac{\partial \eta }{\partial x}+v\frac{\partial \eta }{\partial y}\right)$.

As the incompressibility hypothesis was assumed, the continuity equation is such that:

$\underset{\text{pressure term}}{\underset{︸}{\frac{\partial U}{\partial \xi }+\frac{\partial V}{\partial \eta }}}=0$ (3)

therefore, the governing equations of our computational model, in curvilinear coordinates, are given by Equations (1)-(3).

From the velocity field, obtained by the Equations (1)-(3) and the pollutant transport model provided by:

$\begin{array}{l}\underset{\text{temporal term}}{\underset{︸}{\frac{\partial }{\partial \tau }\left(\frac{C}{J}\right)}}+\underset{\text{convective term}}{\underset{︸}{\frac{\partial }{\partial \xi }\left(UC\right)+\frac{\partial }{\partial \eta }\left(VC\right)}}\\ =+\underset{\text{diffusive term}}{\underset{︸}{D\left[\frac{\partial }{\partial \xi }J\left(\alpha \frac{\partial C}{\partial \xi }-\beta \frac{\partial C}{\partial \eta }\right)+\frac{\partial }{\partial \eta }J\left(\gamma \frac{\partial C}{\partial \eta }-\beta \frac{\partial C}{\partial \xi }\right)\right]}}-\underset{\text{reactive term}}{\underset{︸}{\frac{KC}{J}}}\end{array}$ (4)

where the trajectory of the pollutant concentration in the generated meshes is determined, using the curvilinear coordinate system, described in all blocks.

In the Equation (4), $C=C\left(\xi ,\eta ,\tau \right)$ is the pollutant concentration in time t, where $t=\tau$ because we are not admitting mesh movement [26], D is the diffusion coefficient and K is the decay coefficient [21] [30].

4. Numerical Procedure

The numerical code used in the present work refers to the same code used in Cirilo et al. [26] and Saita et al. [30]. Numerical tests on problems involving a single block [26], were presented in the following cases: of parallel plates, of a cavity, and of the arteriosclerosis.

At [30], the mesh used was of the multiblock type, and the velocity field was calculated numerically in accordance with data taken from the literature [31] [32].

Emphasizing that the generalization of the computational code, with its convergence guarantee requirements, from a single block to multiblock is natural. Once the converged velocity is obtained for each block in agreement with adjacent blocks, for reasons of numerical analysis, the multiblock solution is assumed.

Still, the numerical results presented in [26] [30] [33] showed agreement with the literature for several Reynolds numbers [33] [34] [35] [36] [37] by the vortex center location. Finally, the article of [30] validates the code that is used in the present work.

In particular, the parallel plates problem is presented to evaluate the code. Consider a rectangular geometry with eight meters length and a meter height [26] [34] where the flow is predominantly laminar and practically without vortex formation.

The fluid with R_e = 100 is injected into the mesh, with velocity u = 1.0 m∙s⁻¹ and initial conditions for velocity and pressure are considered null. Boundary conditions CIPR, CECO and CNEI are applied to the left, right, top and bottom sides of the geometry, respectively. For the analysis, was considered the mesh refinements [26] as described at Table 1.

The simulations were performed until reaching a steady state, with a tolerance error equivalent to 10⁻³. In solving the linear system, the Gauss-Seidel method was used. For there to be convergence, it was considered $\Delta t={10}^{-2}$ for the meshes P₁ - P₃ and $\text{Δ}t=5\times {10}^{-3}$ for P₄ and P₅. Considering the mesh P₅, was observed that the speed reaches values close to 1.5 m∙s⁻¹ at the beginning of the simulations [26], but with the progress of the process, this value remains only in the central region of the geometry, as expected. Cirilo et al. [26] present a convergence analysis, in the center of the domain, where for y= 0.5 the maximum analytical velocity is 1.5 m∙s⁻¹ [34]. For mesh refinements, Table 1, the results are presented in Table 2, where ${\text{Error}}_{y=0.5}=\left|1.5-u_{y=0.5}\right|$, $u_{y=0.5}$ represents the values obtained numerically in y= 0.5, in each mesh.

5. Results and Discussion

To evaluate the two-dimensional meshes obtained in the curvilinear coordinate system, besides the velocity field when solving the Navier-Stokes equations and the dispersion of the pollutant concentration by solving the transport equation, two tests were performed. Observe that the model considers that the problem is infinitely extended in the third dimension.

Table 1. Meshes for simulations to parallel plates problem.

Table 2. Convergence speed of the numerical code [26].

In the first test, three cases have been proposed, with meshes containing a chimney followed by an obstacle, using different chimneys heights and the obstacle height was fixed. The test aims to verify the vortices appearance, in the blocks, in order to obtain agreement with as proved in the literature [14]. The second test describes a geometry that represents a valley, which portrays a test of environmental interest. Particularly, in this test, the movement and dispersion of pollutant emissions into the atmosphere of an industry located in the Peruvian city of La Oroya were verified.

5.1. Meshes Containing a Chimney Followed by an Obstacle

The region of interest considers a rectangular geometry composed of a chimney and an obstacle, in which the chimney represents the outflow of pollutants from an industry. This test intends to verify how the trajectory of the pollutants dispersion is affected when considering different heights for the chimney, keeping the obstacle fixed. The generated mesh, referring to the region of interest, has a 3 meters height, a width of 17.05 meters, contains an obstacle height of 1.15 meters and a width of 1.7 meters. As for the chimney, it has a width of 0.4 meters and heights of 0.75, 1.15 and 1.55 meters.

In this work, the mesh was divided into 8 blocks where, Figure 3 represents the geometry where the chimney height is 0.75 meters and illustrates the distribution of the blocks. The geometries of the cases where the chimney height is 1.55 and 1.15 meters, hold the blocks positions. The partitions number in ξ and η directions, for each block, is presented in Table 3.

The boundaries conditions, in the three cases, were kept the same to allow a comparative analysis of the results. Table 4 exhibits the boundary conditions to solve the Navier-Stokes equations, resulting in the velocity field. Table 5 exhibits the boundary conditions to solve the transport equation, describing a concentration of pollutants in the evaluated region.

According to Table 4 and Table 5, velocity input was made on the left boundary of the geometry and at the chimney output, trying to simulate possible environmental impacts caused to the atmosphere by the industry operation.

The required parameters to solve the Navier-Stokes equations are shown in Table 6 and were maintained the same for the three cases, where u and v are the velocity components in ξ and η directions, respectively.

Air specific mass ρ value was based on the normal conditions of temperature and pressure, equivalent to 0˚C and 101,325 Pa, respectively [38].

Figure 3. Blocks positions where the chimney height is smaller than the obstacle height.

Table 3. Parameters for generating the mesh of the first test.

Table 4. Boundary conditions to solve Navier-Stokes equations.

Table 5. Boundary conditions to solve transport equations.

Table 6. Parameters to solve Navier-Stokes equations.

The dynamic viscosity μ value was calculated by $\mu =\rho L{v}_i/R_e$, where v_i are the velocity components in ξ and η directions and R_e corresponds to Reynolds number and L, to the mesh height [39]. Was consider L = 3 in the presented cases and, in order to reach a laminar flow, the Reynolds number equal to R_e= 10 was used. The resolution of the linear system and the convergence criterion are similar to that presented in the section that describes the numerical procedure.

Table 7 shows the parameters used in the transport equation. The parameters were also kept the same for all three cases, and the concentration input border occurs only at the chimney output, as shown in Table 5.

5.1.1. First Case: Chimney Height Smaller than the Obstacle Height

The mesh and the division of the blocks used for the first case are shown in Figure 3, where it can notice that the chimney height, 0.75 meters, is smaller than the obstacle height, 1.15 meters. The parameters required for its generation, the number of partitions in ξ and η directions, are reported in Table 3.

The simulations of the Navier-Stokes equations were executed, until reaching the steady state, with a tolerance error equivalent to 10⁻³ and Δt= 10⁻². From the simulations, the results of the maximum velocity v_max, obtained in each block, can be seen at Table 8. It is verified that the maximum velocity approaches a reference value in each block, v_ref, in time t = 200 s, reaching values close to 4.6723 m∙s⁻¹, 6.6272 m∙s⁻¹, 7.3562 m∙s⁻¹, 0.1734 m∙s⁻¹, 9.4283 m∙s⁻¹, 9.2369 m∙s⁻¹, 2.6330 m∙s⁻¹ and 9.0891 m∙s⁻¹, respectively.

The errors in the temporal evolution of each block, which is defined by ${\text{Error}}_{\text{block}}=\left|v_{\text{ref}}-v_{\max }\right|$, are illustrated in Table 9, in which approach zero, thus guaranteeing the convergence process.

The velocity field for t = 200 is shown in Figure 4. It is noticed that higher velocities, represented by the red color, are located in blocks 5, 6 and 8, where the highest velocities were obtained in block 8, as shown in Table 8. Also, due to the geometry configuration, the highest velocities do not approach the soil. Notice that the boundary conditions, in the numerical interfaces, adequately transfer the information between the blocks, demonstrating the applicability of the used method.

Table 7. Parameters to solve transport equations.

Table 8. Maximum velocity obtained in each of the blocks, for several simulation times.

Table 9. Temporal evolution errors of each block, ${\text{Error}}_{\text{block}}=\left|v_{\text{ref}}-v_{\max }\right|$.

Figure 4. Velocity field, in the final time 200 s, for the first case.

Still, in Figure 4, the contour lines, which describe the velocity field, is shown over the entire mesh extension. Notice vortexes formation between the chimney and the obstacle, at block 4, at the left of the chimney, at block 5, and at the right of the obstacle, at block 7. The appearance of vortexes, in these blocks of geometry, is in agreement with the locations obtained in the literature [14].

The final velocity field, at 200 s shown at Figure 4, result from the Navier-Strokes equations, Equations (1)-(3), and it is used to solve the transport equation, Equation (4). The outcome is shown in Figure 5.

It is possible to notice an accumulation of pollutants in the vortex located between the chimney and the obstacle, at block 4, represented in Figure 5, approximately 30% of the pollution emitted from the chimney. Also, due to the flow pattern, the counter-clockwise recirculation causes the concentration of pollutants to be higher and pollutant removal is very prejudiced. These results are in agreement with the results obtained in the literature [14]. As the pollution trajectory moves away from the chimney the concentration contour approaches the soil, just after the obstacle, becoming smaller as it approximates the right side of the geometry.

5.1.2. Second Case: Chimney Height Bigger than the Obstacle Height

The second case mesh differs from the mesh shown in Figure 3 just by the chimney height, which in this case is 1.55 meters, and the obstacle height is 1.15 meters. The partitions number in ξ and η directions, for each block, is presented in Table 3.

The velocity field for this test, obtained as a simulation result of the Navier-Stokes equations, Equations (1)-(3), is shown in Figure 6. Higher velocities, represented by the red color, are concentrated above the obstacle and decrease as the dispersion area increases.

In this mesh configuration, also can be seen the formation of vortices between the chimney and the obstacle, at block 4, and after the obstacle, at block 7. Similarly, to the first case, the final velocity field, given at 200 s from the Navier-Stokes equations, is used to solve the transport equation, Equation (4). The result is shown in Figure 7.

Due to the chimney height being bigger than the obstacle height, the dispersion of pollutants spreads more rapidly, towards the velocity field. Besides that, a large part of the pollutant concentration is above the obstacle, decreasing quickly as it approaches the right side of the geometry, Figure 7. Also, the pollutant concentration accumulated between the chimney and the obstacle, at block 4,

Figure 5. Concentration of pollutants, in the final time 200 s, for the first case.

Figure 6. Velocity field, in the final time 200 s, for the second case.

Figure 7. Concentration of pollutants, in the final time 200 s, for the second case.

and the pollutant concentration that approaches the soil, after the obstacle, at block 7, are small, varying between 5% and 15% of the pollution emitted by the chimney, with lower values when approaching the soil. These results are obtained due to the flow pattern obtained in this case, in which the height of the chimney, bigger than the height of the obstacle, did not allow the concentration of pollutants to enter the recirculation of the vortex in block 4, but in the flow of the velocity field.

5.1.3. Third Case: Chimney and Obstacle with the Same Heights

Dimensions and blocks identification of the third mesh maintains the same pattern of the first and the second case, only difference is that now the chimney and the obstacle are the same height, 1.15 meters.

Parameters to generate this mesh, partitions in ξ and η directions, for each block, is presented in Table 3. In this test, the velocity field, is presented in Figure 8. Higher velocities, represented by the red color, are concentrated above the obstacle, but unlike the previous cases, the higher velocities, in this configuration, are also observed close to the soil.

Figure 8 shows the contours of the velocity vector across the mesh, in which it is possible to notice the formation of several vortexes.

The formation of vortexes is highlighted, two small vortexes between the chimney and the obstacle with recirculations in counterclockwise and clockwise directions, at block 4. One in the 6th and 7th blocks, above and to the right of the obstacle, and other at the top right of the 8th block, in addition to other

Figure 8. Velocity field, in the final time 200 s, for the third case.

smaller vortexes, which cannot be completely seen due to limitations in the geometry size. The appearance of vortexes in these blocks of the geometry is in agreement with the locations obtained in the literature [14]. The vortex that is top right of the 8th block makes that the velocity field stays more close to the soil. The pollutant concentration trajectory, obtained at the end of the transport equation simulation, presents that, because the chimney and the obstacle are the same height, an accumulation of pollutants between them is observed, as shown in Figure 9.

As in case one, due to the configuration of the geometry, the pollution accumulation in the vortexes located between the chimney and the obstacle causes the concentration of pollutants to be higher and the renovation becomes impaired, varying by 15% and 35% of the pollution emitted by the chimney.

Also, due to the vortex at 7th and 8th blocks, a large part of the concentration of the pollutants remained to the right side of the obstacle, next to the soil, until reaching the limit right of the geometry.

5.2. Mesh of a Geometry Representing a Valley—La Oroya Case

The region of interest, in this test, describes part of a La Oroya city located in a valley of the Peruvian Andes. Numerous studies detect the presence of various types of contaminants in the environment and in the body of its inhabitants [40] [41] [42] [43]. One of the main causes of this pollution is the chimney of the La Oroya metallurgical complex, illustrated in Figure 10, exposing the population to critical levels of air quality, through toxic smoke emitted by its vapor.

With this test, the movement and dispersion of pollutants emitted by the La Oroya chimney are simulated, in order to verify the behavior of the pollution trajectory in the valley.

The generated mesh, which represents the La Oroya valley, has 280 meters height, a width of 760 meters wide at the top and 280 meters at the bottom. As for the chimney, it has 160 meters high and 35 meters wide, and can be seen at Figure 11. The mesh was divided into 8 blocks, the 4th and 5th blocks represent the chimney output. Unlike the other blocks, they are more refined, guaranteeing more accurate results by simulating the movement and dispersion of pollutants. The parameters to generate this mesh, partitions in ξ and η directions for each block, are reported in Table 10.

Table 11 and Table 12 show the boundary conditions to solve the Navier-Stokes equations and transport equations, Equations (1)-(4) respectively.

According to the boundary conditions CIPR, CECO and CNEI, Table 11, the

Figure 9. Concentration of pollutants, in the final time 200 s, for the third case.

Figure 10. Chimney of the La Oroya metallurgical complex (Photo credit: International Federation for Human Rights—FIDH).

Figure 11. Mesh simulating the conditions of La Oroya city, located in a valley and site of a lead smelter.

Table 10. Parameters for generating the mesh of La Oroya case.

left sides of the 1st, 2nd and 3rd blocks in the west direction, receives flow input, while at the east side of 6th and 7th blocks, there is no flow, representing the obstacle. Only the right side of the 8th block, allows the flow output, the east direction. These boundary conditions permit the model to represent a valley. With this test, it intends to evaluate whether the model describes pollutant dispersion characteristics, Equation (4) considering the obtained velocity field, Equations (1)-(3), and whether the transfer of information between the blocks was realized. For the simulations, the existence of houses was despised, because the two-dimensional model describes a vertical section of the region.

The parameters used to solve the Navier-Stokes and transport equations, Equations (1)-(4), respectively, are shown in Table 13 and Table 14.

The velocity field obtained for this case is present in Figure 12. Due to the

Table 11. Boundary conditions to solve Navier-Stokes equations of La Oroya case.

Table 12. Boundary conditions to solve transport equations of La Oroya case.

Table 13. Parameters to solve Navier-Stokes equations of La Oroya case.

Table 14. Parameters to solve transport equations of La Oroya case.

Figure 12. Velocity field, in the final time 10,000 s, for Oroya case.

geometry, the recirculation of the velocity field in block 1 has a flow characteristic similar as lid-driven flow in a square cavity [26]. Still, it is possible to observe higher velocities, represented by the red color, at the chimney output towards the 8th block. The contours identify the formations of vortices which it is possible to observe a large vortex formed between the 6th and 7th blocks.

The results and its contour for the transport equation, Equation (4), are shown in Figure 13. It is observed that the chimney height is smaller than the obstacle, borders to the right of blocks 6 and 7. Thus, this test presents characteristics similar to the first case, in which the vortex formed by the velocity field in blocks 6 and 7, causes the pollution emitted by the chimney to accumulate, with maximum values close to 20% and 50% in the respective blocks, with smaller values close to the soil. As expected, there is a high pollutant concentration at right of the chimney, following the velocity field and vortex formation.

It is observed that there is a high pollutant concentration at the right of the

Figure 13. Concentration of pollutants, in the final time 10,000 s, for Oroya case.

chimney, following the velocity field and the vortices formation. These pollutant concentrations are restricted on this side, due to the geometry characteristic, in which the right border, blocks 6th and 7th, represents the obstacle. Despite the chimney height being smaller than the obstacle height, there is a large concentration of pollutants circulating in this region, because of the large vortex formed, Figure 13. It is observed through the first test, in the first case, that the formation of this vortex was expected.

6. Conclusions

In this work, it is proposed the use of a curvilinear coordinate system and the multiblock technique, generating meshes that coincide with the physical domain, with good quality elements and with no increase of computational cost. The consistency of the fluid dispersion simulation in the meshes divided into blocks was guaranteed by treating the limits that indicate adjacency between the blocks.

In this way, numerical simulations of the movement and dispersion of the pollutant emissions into the atmosphere were presented in the generated domains. For the discretization of the Navier-Stokes and transport equations were considered the curvilinear coordinates and the finite difference method.

The suggested method was verified in two tests. In the first test, three cases were proposed, with geometries containing a chimney followed by an obstacle, using different chimney heights and the obstacle height was fixed. Qualitatively, the results obtained showed an expected behavior, resulting in the formation of vortexes in the predicted locations of the velocity field, and higher concentration of pollutants accompanying the velocity vectors. Also, in this test, were maintained the same boundary conditions, which allowed us to conclude that the trajectory of the pollutants is affected by the height of the source in relation to the obstacle.

Through the performed tests, were verified that the height of the chimney can be considered a determining factor to describe the dispersion of pollutants, as well as their concentration in the proximity of industrial areas.

Acknowledgements

The authors are grateful to the Department Mathematics of the State University of Londrina and the Araucária Foundation under grant n˚ 39225, for the encouragement and providing facilities to carry out the present study.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Seinfeld, J.H. and Pandis, S.N. (2012) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 3th Edition, John Wiley & Sons, New York
[2]	Maione, M., Fowler, D., Monks, P.S, Reis, S. Rudich, Y, Williams, M.L. and Fuzzi. S. (2016) Air Quality and Climate Change: Designing New Win-Win Policies for Europe. Environmental Science & Policy, 65, 48-57. https://doi.org/10.1016/j.envsci.2016.03.011
[3]	Popoola, L.T., Adebanjo, S.A. and Adeoye, B.K. (2018) Assessment of Atmospheric Particulate Matter and Heavy Metals: A Critical Review. International Journal of Environmental Science and Technology, 15, 935-948. https://doi.org/10.1007/s13762-017-1454-4
[4]	Cavalcanti, P.M.P.S. (2010) Air Quality Management Model—Preventive and Corrective Approach. Ph.D. Thesis, (in Portuguese). Universidade Federal do Rio de Janeiro, Rio de Janeiro.
[5]	McKain, K., Wofsy, S.C., Nehrkorn, T., Eluszkiewicz, J., Ehleringer, J.R. and Stephens, B.B. (2012) Assessment of Ground-Based Atmospheric Observations for Verification of Greenhouse Gas Emissions from an Urban Region. Proceedings of the National Academy of Sciences of the United States of America, 109, 8423-8429. https://doi.org/10.1073/pnas.1116645109
[6]	Singh, K.P. Gupta, S. and Rai, P. (2013) Identifying Pollution Sources and Predicting Urban Air Quality Using Ensemble Learning Methods. Atmospheric Environment, 80, 426-437. https://doi.org/10.1016/j.atmosenv.2013.08.023
[7]	Pfeffer, M.A., Langmann, B., Heil, A. and Graf, H.-F. (2012) Numerical Simulations Examining the Possible Role of Anthropogenic and Volcanic Emissions during the 1997 Indonesian Fires. Air Quality, Atmosphere & Health, 5, 277-292. https://doi.org/10.1007/s11869-010-0105-4
[8]	Zeebe, R.E., Ridgwell, A. and Zachos, J.C. (2016) Anthropogenic Carbon Release Rate Unprecedented during the Past 66 Million Years. Nature Geoscience, 9, 325-329. https://doi.org/10.1038/ngeo2681
[9]	Yue, X.-L. and Gao, Q-X, (2018) Contributions of Natural Systems and Human Activity to Greenhouse Gas Emissions. Advances in Climate Change Research, 9, 243-252. https://doi.org/10.1016/j.accre.2018.12.003
[10]	Zannetti, P. (1988) Air Pollution Modeling: Theories, Computational Methods and Available Software. Springer, US.
[11]	Aggarwal, A., Haritash, A. K. and Kansal. G. (2014) Air Pollution Modelling—A Review. International Journal of Advanced Technology in Engineering and Science, 2, 355-364.
[12]	Issakhov, A. and Mashenkova, A. (2019) Numerical Study for the Assessment of Pollutant Dispersion from a Thermal Power Plant under the Different Temperature Regimes. International Journal of Environmental Science and Technology, 16, 6089-6112. https://doi.org/10.1007/s13762-019-02211-y
[13]	Fatahillah, A., Masyhudi M.A. and Setiawan, T.B. (2020) Numerical Analysis of Air Pollutant Dispersion in Steam Power Plant Area Using the Finite Volume Method. Journal of Physics: Conference Series, 1490, Article ID: 012002. https://doi.org/10.1088/1742-6596/1490/1/012002
[14]	Madalozzo, D.M.S., Braun, A.L., Awruch, A.M. and Morsch, I.B. (2014) Numerical Simulation of Pollutant Dispersion in Street Canyons: Geometric and Thermal Effects. Applied Mathematical Modelling, 38, 5883-5909. https://doi.org/10.1016/j.apm.2014.04.041
[15]	Nguyen, V.T., Nguyen, T.C. and Nguyen, J. (2019) Numerical Simulation of Turbulent Flow and Pollutant Dispersion in Urban Street Canyons. Atmosphere, 10, 683. https://doi.org/10.3390/atmos10110683
[16]	Fortuna, A.O. (2000) Computational Techniques for Fluid Dynamics: Basic Concepts and Applications (in Portuguese). Editora da Universidade de Sao Paulo, Sao Paulo.
[17]	Thomas, J.L. and Walters, W.R. (1985) Upwind Relaxation Algorithms for the Navier-Stokes Equations. 7th Computational Physics Conference, Cincinnati, 15-17 July 1985, 117-128. https://doi.org/10.2514/6.1985-1501
[18]	Maliska, C.R. (1995) Heat Transfer and Computational Fluid Dynamics. LTC Editora, Rio de Janeiro.
[19]	Wang, P. and Mu, H. (2010) Numerical Simulation of Pollutant Flow and Dispersion in Different Street Layouts. International Journal of Environmental Studies, 67, 155-167. https://doi.org/10.1080/00207231003654227
[20]	Anunciacao, M., Pinto, M.A.V. and Neundorf, R. (2020) Solution of the Navier-Stokes Equations Using Projection Method and Preconditioned Conjugated Gradient with Multigrid and ILU Solver. Revista Internacional de Métodos Numéricos para Cálculo y Diseno en Ingeniería, 36, 17. https://doi.org/10.23967/j.rimni.2020.01.003
[21]	Romeiro, N.M.L., Castro, R.G.S., Cirilo, E.R. and Natti, P.L. (2011) Local Calibration of Coliforms Parameters of Water Quality Problem at Igapó I Lake, Londrina, Paraná, Brazil. Ecological Modelling, 222, 1888-1896. https://doi.org/10.1016/j.ecolmodel.2011.03.018
[22]	Zong, Z., Xu, M., Xu, J. and Lv, X. (2018) Improvement of the Ocean Pollutant Transport Model by Using the Surface Spline Interpolation. Tellus A: Dynamic Meteorology and Oceanography, 70, 1-13. https://doi.org/10.1080/16000870.2018.1481689
[23]	Singh, A., Das, S., Ong, S.H. and Jafari, H. (2019) Numerical Solution of Nonlinear Reaction-Advection-Diffusion Equation. Journal of Computational and Nonlinear Dynamics, 14, Article ID: 041003. https://doi.org/10.1115/1.4042687
[24]	Ikohagi, T. and Shin, B. (1991) Finite-Difference Schemes for Steady Incompressible Navier-Stokes Equations in General Curvilinear Coordinates. Computers & Fluids, 19, 479-488. http://www.sciencedirect.com/science/article/pii/004579309190072P https://doi.org/10.1016/0045-7930(91)90072-P
[25]	Nikitin, N. (2006) Finite-Difference Method for Incompressible Navier-Stoks Equations in Arbitrary Orthogonal Curvilinear Coordinates. Journal of Computational Physics, 217, 759-781. https://doi.org/10.1016/j.jcp.2006.01.036
[26]	Cirilo, E.R., Barba, A., Natti, P.L. and Romeiro, N.M.L. (2018) A Numerical Model Based on the Curvilinear Coordinate System for the MAC Method Simplified. Semina: Ciências Exatas e Tecnológicas, 39, 87-98. http://dx.doi.org/10.5433/1679-0375.2018v39n2p87
[27]	Cirilo, E.R. and Bortoli, A.L. (2006) Generation of the Trachea and Bronchial Tubes mesh by Cubic Splines (in Portuguese). Semina: Ciências Exatas e Tecnológicas, 27, 147-155.
[28]	Romé, C. and Glockner, S. (2005) An Implicit Multiblock Coupling for the Incompressible Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 47, 1261-1267. https://doi.org/10.1002/fld.898
[29]	Almeida, J.O., Lobao, D.C., Stampa C.S. and Alvarez, G.B. (2018) Multi-Block Technique Applied to Navier-Stokes Equations in Two Dimensions. Semina: CiênciasExatas e Tecnológicas, 39, 115-124. http://dx.doi.org/10.5433/1679-0375.2018v39n2p115
[30]	Saita, T.M, Natti, P.L., Cirilo, E.R., Romeiro, N.M.L, Candezano, M.A.C., Acuna, R.B. and Moreno, L.C.G. (2018) Numerical Simulation of Fecal Coliform Dynamics in Lake Luruaco, Colombia. Tema, 18, 2017, 435-447. https://doi.org/10.5540/tema.2017.018.03.435
[31]	Chapra, S. (1997) Surface Water-Quality Modeling. McGraw-Hill, New York.
[32]	Pardo, S.R., Natti, P.L., Romeiro, N.M.L. and Cirilo, E.R. (2012) A Transport Modeling of the Carbon Nitrogen Cycle at Igapó I Lake-Londrina, Paraná State. Acta Scientiarum. Technology, 34, 217-226. https://doi.org/10.4025/actascitechnol.v34i2.11792
[33]	Hou, S., Zou, Q., Chen, S., Doolen, G. and Cogley, A. (1995) Simulation of Cavity Flows by the Lattice Boltzmann Method. Journal of Computational Physics, 118, 329-347. https://doi.org/10.1006/jcph.1995.1103
[34]	Bono, G., Lyra, P.R.M. and Bono G.F.F.O. (2011) Numerical Solution of Incompressible Flows with Large Scale Simulation (in Portuguese). Proceedings of Conference Mecánica Computational, Argentina, 1-4 November 2011, 1423-1440.
[35]	Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and Multigrid Method. Journal of Computational Physics, 48, 387-411. https://doi.org/10.1016/0021-9991(82)90058-4
[36]	Gupta, M.M. and Kalit, J.C. (2005) A New Paradigm for Solving Navier-Stokes Equations: Stream Function-Velocity Formulation. Journal of Computational Physics, 207, 52-68. https://doi.org/10.1016/j.jcp.2005.01.002
[37]	Marchi, C.H., Suero, R. and Araki, L.K (2009) The Lid-Driven Square Cavity ow: Numerical Solution with a 1024 × 1024 Grid. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 31, 186-198. http://dx.doi.org/10.1590/S1678-58782009000300004
[38]	Taylor, J. (2013) Classical Mechanics. Bookman Editora University Science Books, Mill Valley.
[39]	Reynolds, O. (1983) An Experimental Investigation of the Circumstances which Determine Whether the Motion of Water Shall Be Direct or Sinuous, and of the Law of Resistance in Parallel Channels. Philosophical Transactions of the Royal Society of London, 174, 935-982. https://doi.org/10.1098/rstl.1883.0029
[40]	Cooke, C.A. and Abbott, M.B. (2008) A Paleolimnological Perspective on Industrial-Era Metal Pollution in the Central Andes, Peru. Science of the Total Environment, 393, 262-272. https://doi.org/10.1016/j.scitotenv.2007.12.034
[41]	Castro-Bedrinana, J. Chirinos-Peinado, D. and Ríos-Ríos, E. (2013) Niveles de plomoen gestantes y neonatos enlaciudad de laOroya, Peru. Revista Peruana de Medicina Experimental y Salud Publica, 30, 393-398.
[42]	Orihuela, J.C. (2014) The Environmental Rules of Economic Development: Governing Air Pollution from Smelters in Chuquicamata and La Oroya. Journal of Latin American Studies, 46, 151-183. https://doi.org/10.1017/S0022216X13001545
[43]	Neumann, P. (2016) Toxic Talk and Collective (In)action in a Company Town: The Case of La Oroya, Peru. Journal of Latin American Studies, 63, 431-446. https://doi.org/10.1093/socpro/spw010