A CFD-Based Methodology for Impact Assessment of Industrial Emissions and Evaluation of Mitigation Measures for Regulatory Purposes

Abstract

In a context where air quality has become a global concern, modelling techniques are becoming very popular for analysing pollutant dispersion conditions. While models based on empirical formulations are most widely used for industrial scenarios, singular cases involving complex terrain or large obstacles in the vicinity of emission sources require a more robust approach to evaluate the dispersion conditions. In this research, a computational fluid dynamics (CFD) model is developed to analyse the dispersion of pollutants from an industrial facility whose location and characteristics suggest the occurrence of complex flow features that affect the dispersion patterns. In addition, the variation of the gas exit temperatures of waste heat recovery systems is proposed as a measure to mitigate the impact of the plume. The simulation results show unexpected deviations in the plume path affecting vulnerable areas in a nearby mountainside, but increasing the gas exit temperature is useful to prevent this event. Therefore, the proposed methodology can be considered a decision support tool to find a compromise between the environmental impact and the fuel consumption of the plant.

1. Introduction

The rapid industrial and economic development achieved by many regions over the last decade has led to an important expansion of urban areas as well as the construction of many high-rise buildings to cater for the increasing population. According to one of the United Nations’ latest reports [1], 55.3% of the world’s population lived in urban settlements in 2018, and this proportion is expected to increase to 68% by 2050. However, the proliferation of dense built-up areas poses an important problem, as it is known to cause a reduction in the flow circulation and pollutant removal capacity in urban environments. Hence, air quality has become a major concern for governments and citizens, as a significant proportion of urban population is exposed to air that does not meet legislation air quality standards or World Health Organization (WHO) guidelines [2]. Increasing levels of outdoor pollution are associated with a wide range of detrimental effects on the population’s health and comfort [3].

Governments and environmental agencies all around the world are creating policies to reduce the population’s exposure to air pollution, such as the European Union’s air quality directives adopted throughout the last decades. These regulations are intended to promote the development of strategies for the reduction and control of emissions in urban areas, a scenario that has motivated the use of traditional and new techniques to model the dispersion of pollutants in the atmosphere.

Dispersion modelling uses mathematical equations to describe the behaviour of the atmosphere, as well as dispersion mechanisms of pollutants, and they might also include some certain physical and chemical processes. By applying these to different case studies, it is possible to predict pollutant concentrations at various locations of interest. Of these techniques, computational fluid dynamics models (CFD) have become very popular for plume dispersion analysis in microenvironments and also for local air quality assessments [4,5,6]. They can provide the complex analysis of fluid flow based on the conservation of mass and momentum principles by resolving the Navier–Stokes equation. Indeed, the popularity of CFD techniques extends to other fields such as industrial design [7] or biomedicine [8,9,10,11,12]. One of the main advantages of these models is their ability to provide the entire flow and dispersion field data without suffering from similarity requirements [13]. In addition, CFD simulations can save time and money compared to traditional experiments and testing methods by eliminating the need for physical prototypes.

The widespread use of these models in urban physics enables regulators and urban planners to improve urban geometry designs for optimal ventilation and to enhance air quality indicators in cities. Investigations conducted over the last few decades [5,14,15] show that phenomena such as vortex shedding and recirculation regions are commonplace in urban environments and have a remarkable impact on pollutant dispersion, especially in the presence of high-rise structures. These recirculation and stagnation zones emerge in between the buildings where the wind regularly blows and are generally characterised by low flow speeds that cause pollutants to become trapped within these areas [16,17]. This situation is particularly damaging for the fresh-air intake of dwellings, which can be affected by the emissions from the rooftop stacks of nearby buildings because of the influence of such flow recirculation zones. An appropriate CFD dispersion simulation can be useful in the design and positioning of rooftop chimneys and ventilation systems of the buildings to comply with current standards and to minimise the risk of contamination [18,19]. At larger spatial scales, CFD models are also useful for estimating pollutant transport and dispersion over long distances [20].

Not only urban planners, but also industrial agents, benefit from dispersion models in the design and impact assessment of the emissions. In order to evaluate whether their designs comply with regulatory exigences (e.g., Industrial Emissions Directive requirements), industrial designers usually rely on the so-called operational models to simulate pollutant dispersion. Operational models use semi-empirical equations based on parametrization techniques deduced from experimental data [21]. They only require a small amount of input and can treat a variety of situations with relatively low computational cost [15], so they are frequently employed in decision support tools where fast response times are required. Among these tools, Gaussian formulations are some of the most outstanding examples [14,22].

In most cases, industrial facilities are built on flat and unobstructed terrains, normally far from urban centres. Under these assumptions, operational models can provide good performance results at a very low computational cost. They perform well even in the presence of buildings, such as the case of AERMOD [22] or OML [23], for which downwash algorithms were implemented. Indeed, there is an active modelling community committed to the improvement of AERMOD and other U.S. EPA’s (United States Environmental Protection Agency) air quality models through the introduction of state-of-the-art modelling concepts into them. Moreover, as proof of its performance and wide acceptance, many environmental regulatory bodies recommend or even impose the use of a specific operational model to conduct the dispersion analysis required to request or renew the operating permits for the emission sources. This is the case of the U.S. EPA, which accounts for resources like AERMOD or CTDMPLUS [24] in its list of preferred models. Both of the above examples are Gaussian models designed to estimate hourly concentrations of plume material from elevated point sources. However, the detail and complexity of the meteorological and terrain input data requirements are higher in the case of CTDMPLUS, as it is intended for more complex terrain scenarios.

Nevertheless, the accuracy of empirical models is known to decrease under certain specific conditions. For instance, scenarios involving complex or steep terrains are some of the most problematic situations, leading to the deterioration of the model results, as reported by a comparison study [25]. One of the critical factors explaining this loss in the model performance is that empirical formulations do not account for physical process (i.e., the flow features generated by surface–atmosphere interactions), as opposed to CFD techniques, but are based on expressions that assume crude simplifications [26]. As a result, these models are not suited for predicting the flow and concentration fields in the close surrounding area of the stack, being typically used to evaluate concentrations only at distances of at least several hundred meters from the source onwards [24,27,28]. Despite the efforts accomplished to extend the models’ functionalities to account for the unrepresented physical processes, the limitations found in challenging scenarios have not yet been overcome in most cases [29,30].

Therefore, trusting empirical models alone might not be enough to understand the pollution impact of industrial chimneys when these are placed at complex environments subject to the influence of complex flow phenomena. These singular cases should be tackled with a more advanced modelling technique, such as CFD, as otherwise the impacts of recirculation zones and other flow features could not be appropriately considered, potentially undermining the effectiveness of the chimney design.

In this work, a CFD model is developed to study the pollutant dispersion from an industrial facility whose location meets the singularity and complexity criteria previously explained. The target plant is a chemical facility located in a valley with close mountains exceeding the height of the main emission source. Special attention is paid to the occurrence of flow phenomena affecting the behaviour of the plume and the dilution effectiveness of the chimney under these conditions. In addition, the regulation of waste heat recovery systems is proposed as a measure to mitigate the impact of the pollution plume. To evaluate the effectiveness of this measure, new simulations are conducted, varying the temperature of the exit gases to account for the changes in the behaviour of the plume. The following CFD methodology was first developed and tested on an experimental database, allowing for model validation and the justification of its subsequent application to a real case study.

Few examples of studies on the impact of singular locations on pollutant dispersion mechanisms from industrial sources have been previously discussed in the literature. More specifically, to the authors’ knowledge, no previous cases of studies focusing on (i) recirculation flow phenomena caused by natural terrain and (ii) their effects on the behaviour of industrial plumes have hitherto been conducted. Addressing these two specific features simultaneously through the development of a CFD model represents an innovative contribution of this research.

Therefore, the novelty of this research can be found in the study of the impact of large flow recirculation effects on pollutant dispersion when the cause of this phenomenon is the interaction of the wind with natural terrain features. While this phenomenon has been widely studied in the urban environment, there are hardly any studies on natural terrain scenarios, due to the difficulty of correctly modelling the terrain geometry around the pollutant source, which is the main objective of this research.

To overcome this challenge, another contribution of this research is the design of a methodology to develop robust computational domains for the appropriate representation of complex terrain scenarios based on non-uniform rational B-spline (NURBS) surfaces. To achieve this objective, a systematic review of the literature was carried out to identify the best practices for the development of CFD models, both in terms of traditional guidelines and new approaches emerging from more recent studies.

Finally, another relevant contribution of this study is the use of a CFD model to analyse the effect of modifying an operational parameter of an industrial plant on the behaviour of the exhaust gases. Through this analysis, it is possible to adjust the performance of the heat recovery systems of the plant to mitigate the pollution impact caused by unfavourable weather conditions.

2. Materials and Methods

The CFD model methodology to apply will be first evaluated against tracer observations from the Evaluation of Modelling Uncertainty (EMU) project [31]. This evaluation is intended to serve as a validation step for a real-life case study, as the same computational parameters and settings will be applied to the real-life situation. As stated in [32] for the sub-configuration validation strategy: “it can reasonably be assumed that the same or a similar combination will also provide accurate results for the actual configuration”. Indeed, validation is a crucial step to ensuring the accuracy of computational simulation results, as shown in the literature [33,34,35].

The EMU experiments were conducted in a large and stratified wind tunnel (20 m × 3.5 m × 1.5 m) at the University of Surrey, where continuous jet releases of dense, buoyant and neutrally buoyant gases were simulated in neutral or stable atmospheres. Among the available experiments, case B3, consisting of dense gas release in neutral atmospheric conditions was selected for validation, as these conditions are similar to those in the actual case study.

2.1. Computational Domain and Mesh Generation

The first step in the development of a CFD model is the definition of the geometric domain. In case B3, gas was released from a 1 m diameter hole located in a simple L-shaped building. Another slab-shaped building was located downwind of the former. The geometry and dimensions are schematically depicted in Figure 1.

The upstream wind blows along the x-axis. The emission source is therefore located on the leeward side of the L-shaped building, so it is not directly exposed to the incoming wind but may be exposed to recirculation flow generated behind the building.

The computational domain is discretised with a body-fitted, poly-hexcore mesh based on octree decomposition. The poly-hexcore material is composed of a combination of polyhedra, hexahedra and prisms. To select an appropriate mesh size, a grid independency study was conducted based on three different meshes: a coarse grid (2 million cells), a medium grid (3 million cells) and a fine grid (5 million). The aim of this test was to verify whether the simulation results are sensitive to a further refinement of the computational grid. A threshold value of 5% was taken as the maximum relative difference allowed for the results between one mesh resolution and the following. By computing this indicator for the tracer concentrations, the relative differences between the coarse and the medium meshes were about 6%. Also, this difference is slightly below 5% between the medium and the fine grid. According to these results, the fine mesh (5 million cells) was selected for the simulations.

The grid was refined around the building, with a cell size of 0.4 m on the building’s walls and roof. The stack section has the smallest cell size (0.2 m), and the maximum value (2 m) can be found far away from the building, at the model’s boundaries. On the other hand, an inflation strategy was applied to produce three layers of prisms, adjacent to the solid walls of the domain. The purpose of these inflation layers was not to capture the boundary layer developed on solid surfaces, as they are usually intended for when applying low-Reynolds number modelling. Instead, this strategy was intended to limit the height of the first row of wall-adjacent cells in the polyhedral mesh design, as well as to control the stretching ratio of the surrounding cells in the vertical direction to a maximum of 1.15. In addition, the faces of these prismatic cells are either parallel or perpendicular to the walls, which is advisable for the application of wall functions. Figure 2 shows three views of the mesh obtained.

2.2. Boundary Conditions

The different conditions applying to each of the model boundaries are shown in Figure 3.

It should be noted that the various parameters used to set up the model are defined in the corresponding project report [31].

At the inflow boundaries of the domain, profiles of mean wind velocity, turbulence kinetic energy and turbulence dissipation were imposed following the Expressions (1)–(3).

$$U\left(z\right)=\frac{{\mathrm{u}}^{\mathrm{*}}}{\kappa }\ln \left(\frac{z+z_0}{z}\right)$$

(1)

$$k\left(z\right)=\frac{{\mathrm{u}}^{{\mathrm{*}}^2}}{C_{\mu }^{0.5}}$$

(2)

$$\mathsf{\epsilon }\left(z\right)=\frac{{\mathrm{u}}_{}^{{\mathrm{*}}^3}}{\kappa \left(z+z_0\right)}$$

(3)

where u* is the friction velocity, κ the von Karman constant (0.42), z is the height, z₀ is the aerodynamic roughness length and C_μ is a constant equal to 0.09. In addition, zero species background concentration was defined at the inlet planes.

The pressure outlet condition was selected for the outflow boundaries where a height dependent pressure profile was specified according to Expression (4).

$$P\left(z\right)=\rho gz$$

(4)

where ρ is the air density, g is the gravitational acceleration and z is the height.

The emission source is represented through an inlet condition with uniform velocity of 25 m/s. A mixture of 3.51% by volume of ethylene in a krypton balance is released from the 1 m diameter source in the building’s wall, for which ideal gas properties were considered.

At the top of the domain, fixed values for U, k and ɛ were applied, corresponding to the values of the inlet profiles at this height. Lateral boundaries were assigned symmetry conditions. At the walls, the standard wall function with a roughness modification [36,37] was imposed, for which the aerodynamic roughness length (z₀) equals 0.12 m. It is important to note that this function is based on a different roughness parameter, k_s and a roughness constant C_s, which defaults to 0.5. However, the relationship between k_s and z₀ is shown in (5), which can be derived by the first-order matching of the atmospheric boundary layer (ABL) velocity profile with the wall-function velocity profile in the centre point P of the wall-adjacent cell [38].

$$k_{S,ABL}=\frac{9.793z_0}{Cs}$$

(5)

2.3. Numeric Approach and Turbulence Model

Among the different existing turbulence models reported in the literature for CFD, the Reynolds-averaged Navier–Stokes equations (RANS) and large eddy simulation (LES) are the most outstanding ones. RANS models solve the ensemble averaged Navier–Stokes equations by turbulence modelling, which aims to model the entire spectrum of the turbulent eddies. When flow is assumed to be statistically steady, ensemble averaging is equivalent to time averaging (steady-RANS) [39]. These have traditionally been the first choice to simulate a wide range of atmospheric [6,13,40] and pollution applications [5,41,42], providing accurate solutions at moderate to low computational cost. The availability of best practices and well-developed technical guidelines for RANS models [43,44] is another sign that these models have reached a mature state.

Thus, the RANS approach was selected for the present study, together with realizable k-ɛ model for the closure of the system of equations. Realizable k-ε models are frequently applied in urban studies involving complex flow features [6,40,41,45], where they are shown to yield better results in the prediction of the recirculation zones around buildings than other models, such as the standard k-ε.

All the transport equations were discretised with a second-order scheme. The SIMPLE algorithm is used for pressure–velocity coupling. The gradients were computed based on the Green–Gauss theorem, and initial conditions were computed from inflow profiles. A threshold of 10⁻⁴ was set for all residuals to monitor the progress of the solution, except for the energy equation, where a value of 10⁻⁶ was specified. In addition, it was necessary to analyse the results in order to detect possible undesired oscillations requiring specific actions.

2.4. Validation and Results

Figure 4 shows a comparison between predicted and measured tracer concentrations at two cross-stream planes (x/H = 1 and x/H = 10) and for two different heights (z/H = 1 and z/H = 10) where H is the height of the building. The distances and concentrations in the experiment were normalised to facilitate comparisons between the model in the wind tunnel and the full scale one. For greater clarity, Figure 5 shows the location of the sections where the results were evaluated.

Overall, the realizable k-ɛ captures the general trends and peak concentrations, despite certain overpredictions at certain points of the domain. This evaluation serves as a strong basis for simulating the dispersion of the plume in realistic urban environments.

3. Case Study

The case study selected for this research is a chemical plant located in a mining valley in the Principality of Asturias, a province in the north of Spain. The main activity of the plant is the manufacturing of coke as well as the production of coke oven gas as a subproduct of the process. The facility is surrounded by a highly mountainous area, as shown by Figure 6, and it is located less than 500 m from an urban nucleus.

The Integrated Environmental Authorization (IEA) (Integrated Environmental Authorization revision of May 2019 (BOPA, 2019). URL: https://sede.asturias.es/bopa/2019/05/17/2019-04859.pdf (accessed on 1 June 2023)) of the facility provides relevant information about its features, operating conditions and sources of atmospheric releases. Among the different emission sources identified in the IEA form, the one for the evacuation of the fumes of the gas combustion in the boiler was selected for analysis. The source is a prominent chimney (elevated point source), 56 m high and 2.96 m in diameter, located at the geographic coordinates 43.288439, −5.674461. The proximity of the chimney to the nearby hills that widely exceed the height of the stack is a singular feature to highlight, as a chance exists that the exhaust gases are affected by the wind flow’s interaction with the land surface. The proximity of the urban nucleus reinforces the need to examine the dispersion mechanisms to detect potentially risky pollution events. Indeed, this source recently underwent a dispersion modelling study to verify the dilution effectivity of the chimney’s height, as stated in one of its IEA revisions.

Figure 7a shows a 3D representation of the plant outbuildings surrounded by the natural terrain, seen from the north-west direction. A large hill exists south of the main chimney (red body) whose vertical dimensions far exceed the height of the outlet section (represented by a red line projected onto the surrounding terrain surface). Moreover, these high landforms are significantly close to the source, as indicated by the horizontal distances represented in Figure 7b.

3.1. Computational Domain and Mesh Generation

Terrain modelling in CFD studies has traditionally been an important challenge that has been addressed in different ways in the literature [46,47,48,49,50]. When dealing with terrain elevation data, a common approach is to generate a triangulated irregular network (TIN). Although this technique is advantageous in terms of simplicity and fast results, it often suffers from geometric inaccuracies, especially in the case of very large distances with high spatial resolution. These geometry errors have a significant impact on the quality of the subsequent mesh, thus compromising the reliability of the simulations.

To address this problem, this research proposes an alternative technique based on the definition of Bézier curves or splines from terrain elevation data. A sweeping process is then applied using the splines as cross sections to generate a non-uniform rational B-spline (NURBS) surface.

These surfaces are mathematical representations of 3D geometry that can accurately describe any shape. The geometry of the buildings and the wall of the chimney are also included in this step.

Although this method is more time-consuming than TIN surfaces, the good results obtained in terms of geometric quality justify its use. Moreover, this method allows the generation of a continuous surface throughout the whole case study. This is an advantage because it makes it possible to obtain the air section to be modelled (computational domain) as a single volume composed exclusively of the fluid space encompassing the study area. For this purpose, a rectangular parallelepiped is generated and sliced at the bottom by the NURBS surface. This computational domain format facilitates the subsequent grid generation process, with satisfactory results.

The terrain elevation data of the area of interest were retrieved from the Spanish National Centre for Geographic Information (CNIG) (National Center for Geographic Information (CNIG). Data consultation for “Digital Elevation Models: LIDAR 1st Coverage (2008–2015)” [File: PNOA-2012-LOTE-AST-282-4798-ORT-CLA-COL.LAZ] url: http://centrodedescargas.cnig.es/CentroDescargas (accessed on 3 September 2022)) in LiDAR point cloud format, which offers versatility for the creation of three-dimensional bodies. These data were used to create a continuous NURBS surface, representing the topographical features of the terrain. The region selected for the analysis is a circular area of 1 km radius, with the main emission source of the industrial facility located at its geometric centre. Next, the bodies of the chimney and the buildings located within the study area were added to the model in the form of simplified parallelepipeds and cylinders. The geometry of the buildings was defined based on the layout information from CADMAPPER service, together with the classification properties of LiDAR point clouds, which allow for the exploration of the elevation of the points corresponding to the rooftops. The dimensions of the chimney was extracted from its IEA. Figure 8 shows the final appearance of the modelled terrain surface, together with the plant dependencies and buildings.

When designing a computational domain for the simulation of the atmospheric boundary layer, there are some well-known guidelines to consider [43,44]. According to these, the computational domain must be large enough to ensure homogeneous inlet and outlet conditions, as well as to avoid unintended flow gradients due to the flow contraction effect. For that purpose, the model extents are enlarged through the addition of a horizontal plane constituting the bottom or lower limit of the computational domain. This plane bottom is shaped according to the size guidelines given in [44], ensuring that the blockage ratio of the domain is kept below 3% [43]. The recommended minimum dimensions are defined based on the vertical size of the highest obstacle of the model (H_max). In this case study, the H_max parameter corresponds to the elevation of the highest point of the terrain surface with respect to the flat domain bottom.

As well as following these general recommendations, the study has also incorporated other practices from the more recent literature. One of the key features of the case study that requires special attention is the transition between the flat bottom of the domain and the outer boundary of the terrain, the latter having an irregular contour with strong height gradients typical of a valley area. In order to ensure a suitable and progressive transition of the flow from the flat bottom to the terrain area, a smooth transition surface was generated between these two elements with a maximum slope of 30°, in accordance with the findings of [51].

According to the explained scheme, the domain can be virtually divided into an inner domain (including the terrain area and the transition surface) and an outer domain (corresponding to the expanded flat background). Figure 9 shows the different parts and dimensions of the resulting domain schematically, which was fully developed using the commercial computer-aided design software Autocad Civil 3D.

For the domain discretization, a poly-hexcore mesh of 20 million cells was implemented, following a strategy similar to the one explained for the validation case. In this case, three meshes of 10 million cells (coarse), 20 million cells (medium) and 30 million cells (fine) were selected for the grid independency study. Following the same criteria as in the validation case, the difference in the concentration results between a given mesh resolution and the following falls below 5% for the medium case, which was the one selected for the simulations.

The grid was refined around the emission source, where a cylindrical shape of twice the height of the stack and 6 m diameter was defined to limit the size of the cells inside this region up to a maximum of 1.5 m. This configuration ensured that the jet-dominated part of the plume would develop in a high-resolution mesh region [42]. In addition, the outlet face of the chimney was covered with 36 cells of a maximum size of 0.4 m, as previous studies highlight the importance of refining the mesh on the exhaust section [52].

The cells covering the buildings and terrain surfaces were assigned maximum sizes of 4 and 2.5 m, respectively, with further refinements of up to 2.5 m for terrain to capture proximity to building edges. Away from the buildings, the grid grew gradually to maximum sizes of 10 m within the inner domain and up to 40 m in the top boundary of the outer domain. However, the cells in the bottom of the outer domain were limited to 10 m in the upstream part of the model and to 20 m in the downstream region. The expansion ratio between any two consecutive cells of the domain is at most 1.2, which is key for keeping the truncation error small. Moreover, this ratio decreases locally to 1.1 in critical regions such as the vicinity of the chimney and the buildings’ surfaces. Finally, three inflation layers of prisms were added adjacent to the solid walls of the domain, limiting the stretching ratio to a maximum of 1.5 in the vertical direction. Figure 10 shows the final appearance of the computational grid.

The resulting mesh presents a sufficiently high resolution to capture the structure of boundary layer dynamics and shows high-quality indicators, including a satisfactory level of cell skewness (0.49), which is important to achieve the convergency of results and low discretization errors during the simulations.

3.2. Simulation Scenarios and Boundary Conditions

When conducting modelling studies of the atmosphere, weather information is key to identifying and selecting the most interesting meteorological scenarios for the research objectives. For that purpose, data were retrieved from the Meriñán weather station (Air quality monitoring stations network of the Principality of Asturias. Hourly data consultation for Meriñán station. url: https://asturaire.asturias.es/consultas (accessed on 4 January 2023)), located around 3 km northwest of the facility. Based on this information, the wind frequency roses of each season were generated and are shown in Figure 11.

From the plots, it can be concluded that south-west is the prevailing wind direction most of the year but especially during the cold seasons (autumn and winter), followed by the north-east direction, which prevails during the warm seasons (spring and summer). The wind conditions are mostly calm throughout the year; however, moderate to high wind velocities come mainly from south and south-west directions in the winter and to a lesser extent in spring and autumn. In spring, moderate velocities also occur for the prevailing wind direction of that season (north-east).

This south-west direction is especially interesting for the purpose of this work, as it leads to a situation in which the chimney is located downstream of the nearby hill, potentially exposed to a flow recirculation zone that might develop behind the cliff. This situation might also give rise to the entrapment of the exhaust pollutants, thus negatively affecting the dispersion mechanisms and the evacuation capacity.

In terms of wind velocity, as the focus of this work is on the formation of recirculation zones and swirls, only moderate to high speed conditions, under which these flow features intensify, will be considered for the analysis. Therefore, only wind speed values of at least 6 m/s will be simulated. It should be noted that, for such wind velocities, the atmospheric boundary layer exhibits neutral stratification, and thus thermal effects are either absent or negligible.

According to these criteria, the meteorological scenarios shown in

Table 1. Meteorological scenarios selected for the simulations.

Scenario	Uref at 10 m	Incident Wind Direction	Seasonal Mean Temperature
S1.0	6 m s−1	210° (SSW)	7.5 °C (winter)
S2.0	10 m s−1	210° (SSW)	7.5 °C (winter)

were selected for the simulations. Note that a speed of 6 m/s represents a moderate velocity, whereas 10 m/s is considered a high or strong wind. It should also be noted that the combination of high wind speeds and low air kinematic viscosity results in a high Reynolds number in the region of 10⁵.

The lateral planes of the domain act as either inlets or outlets for the flow, depending on the direction of the incident wind considered in each scenario with respect to the domain boundaries. For a wind flow coming from SW direction, the south and west boundaries of the domain will be treated as inlets, while the north and east boundaries will be the outlets. The outer dimensions applied during the definition phase of the computational domain should be set in accordance with the relative position of each boundary type (inflow or outflow) resulting from each wind scenario.

At the inflow boundaries of the domain, profiles of mean wind velocity, turbulence kinetic energy and turbulence dissipation are imposed based on Richards and Hoxey’s set of equations [53] for a neutral ABL, as shown by Expressions (6)–(8). These profiles ensure a homogeneous ABL and constitute a standard among CFD practitioners for wind engineering studies.

$$U(z)=\frac{u_{ABL}^{*}}{\kappa }\mathit{ln}\left(\frac{z+z_0}{z}\right)$$

(6)

$$k(z)=\frac{u_{ABL}^{{*}^2}}{C_{\mu }^{0.5}}$$

(7)

$$\mathsf{\epsilon }(z)=\frac{u_{ABL}^{{*}^3}}{\kappa \left(z+z_0\right)}$$

(8)

where u^*_ABL is the ABL friction velocity. As in the validation case, zero species background concentration was defined at the inlet planes.

As recommended in aerodynamic studies, an upstream fetch of 5 km was considered for the definition of the atmospheric boundary layer properties [54]. The terrain included within a radial distance of 5 km around the facility was analysed to extract the aerodynamic roughness length z₀, as shown by Figure 12. This roughness parameter was determined based on the updated Davenport roughness classification [54] at every 30° sector, in accordance with the divisions of the wind roses.

For the rest of the model boundaries, namely the top, wall and outlet boundaries, the boundary conditions selected are similar to the respective ones applied to the validation case. The source of study is the large chimney dedicated to the evacuation of the fumes resulting from the gas combustion in the boiler. As stated in the plant’s IEA, the gas emission occurs on a continuous basis 24 h a day and the main pollutants released into the atmosphere are identified together with their correspondent emission limits in its atmospheric emissions annex. These limits were taken as the reference values to set the proportion of each pollutant in the fumes. The use of waste heat recovery systems is envisaged in the IEA of the facility, whose purpose is to reduce the fuel consumption of the plant, increasing the energy efficiency. However, reducing the temperature at which gases are expelled from the chimney comes at the cost of a consequent decrease in the exit velocity, which eventually means a lower effective emission height. The impact of heat recovery systems on the effective emission height and its consequences on pollution levels were already investigated in [55]. In the study, the authors set different stack exit temperatures to analyse the effect on near-source air quality and defined an indicator to quantify the impact of heat recovery systems.

As for the present study, an initial or reference gas exit temperature of 150 °C was taken, which assumes the full utilization of the heat recovery systems. Note that gas temperatures under this value are normally avoided as they can lead to acid condensates inside the chimney’s structure. Then, additional simulations were conducted, increasing the temperature of the gases to 200 °C and 250 °C, assuming the regulation or partial use of the devices. The new simulations are useful for evaluating the impact of waste heat recovery of gases on the dispersion mechanisms, making it possible to set different heat conversion rates to ensure a minimum level of ambient air quality for each weather condition.

On the other hand, exhaust gas velocity is typically in the recommended range of 5–20 m/s. This variable can be calculated for the different temperature conditions according to Equation (9).

$$v_{eg}=C\sqrt{2gH\frac{T_i-T_e}{T_i}}$$

(9)

where C is the discharge coefficient (0.65), g is the gravity acceleration, H is the chimney’s height, $T_i$ is the gas exit temperature and $T_e$ is the outdoor temperature. The resulting velocities for the different gas exit temperatures are shown in

Table 2. Exit gas velocities in m/s for the different temperature conditions.

Ti = 250 °C	Ti = 200 °C	Ti = 150 °C
14.7	13.7	12.5

, considering an outdoor temperature of 7.5 °C, as stated before.

The exit cross-section of the chimney is defined as a velocity inlet with a uniform velocity profile. Turbulence quantities are computed based on the hydraulic diameter, assuming a 10% turbulence intensity.

3.3. Numeric Approach and Turbulence Model

As supported by the tests performed in the validation case, a steady RANS approach was selected, together with the k-ε realizable turbulence model. The values of the model constants are C_1ɛ = 1.44, C_2ɛ = 1.92, C_μ = 0.09 and σ_k = 1.0, σ_c = 1.11, with the latter imposed in accordance with Richards and Hoxey’s inflow profiles [53]. The remaining model settings, such as the discretization scheme, residual thresholds and other numeric parameters, are the same as in the validation case.

4. Results and Discussion

In order to analyse and draw better conclusions from the simulation results, a set of visualization graphs was obtained. Initially, the occurrence of flow recirculation phenomena was confirmed according to the following analysis. Figure 13 shows the velocity magnitude contours on a plane parallel to the incident wind direction passing through the chimney’s centre point. In addition, Figure 14 and Figure 15 show, respectively, the flow path lines coloured by velocity magnitude in Scenario S1.0 and Scenario S2.0.

In all cases, the formation of flow recirculation region starting at the highest point of the hill and evolving along the downhill slope is clearly visible. When it reaches the position of the chimney, a swirling motion can be identified, which clearly impacts the behaviour of the plume. For scenario S1.0 (210°, 6 m/s wind), the speed of the plume is higher than the velocity of the flow in the vicinity of the exhaust section, allowing the gases to rise a few meters before being diverted by the swirl motion. As the exhaust gases lose their velocity, the plume turns slightly southwards due to the upward motion of the vortex, and then starts to bend northwards with a much lower speed. The observed behaviour looks notably different in Scenario S2.0 (210°, 10 m/s wind) since the wind speed at the exhaust section is now comparable to the exit velocity of the gases. In this case, the plume barely rises, and it is soon dragged southwards with a downward movement. Since the background flow velocities around the chimney’s outlet section are now similar or even exceed those of the exhaust plume from a certain height, the latter is much more influenced by the wind.

To assess the influence of the flow dynamics on the pollutant plume, one of the most useful visualization forms that can be obtained from CFD results is the representation of iso-surfaces. The pollutant SO₂ was chosen for the visualization of concentration iso-surfaces for different fractions of the emission levels at the exit cross-section. Figure 16 shows SO₂ iso-concentration surfaces of 1%, 5% and 25% of the source emission levels for the two meteorological scenarios.

The flow recirculation phenomenon and its impact on the plume dispersion can be seen in these representations. In Scenario S1.0, the upward movement of the plume and its subsequent 180° turn can be seen, followed by a downward trajectory on the west side of the chimney. The downward bending motion of the plume causes the pollutants to approach the mountainside to the southwest of the chimney, where a few small houses and farms not explicitly represented in the model geometry are located. This situation may, therefore, have a negative effect on the quality of the air that people living in that area of the mountain breathe. For the lowest concentration fraction, the plume also invades the industrial installations themselves, even if it does not reach ground level. Although the concentration levels in the latter case might not be significant, possible consequences for the employees of the plant resulting from long-term exposure should not be overlooked.

For Scenario S2.0, the plume exhibits a different behaviour as the stronger winds barely allow it to rise so that the pollutant gases no longer approach the mountainside that was previously affected. In contrast, the industrial facility now receives the greatest impact as the plume falls rapidly downwards, with an almost vertical trajectory, reaching the ground at a higher concentration than in the previous scenario.

As a measure to try to mitigate the impact of the pollutant plume, new simulations were run by varying the exit gases temperature based on the regulation of the heat conversion rates of the facility. Figure 17 shows the SO₂ iso-concentration surfaces for 5% of the emission levels recorded at the exit cross-section, considering the two meteorological scenarios and the different gas outlet temperatures (150 °C, 200 °C and 250 °C).

It can be seen from the figure that increasing the exhaust temperature of gases can have a major impact on the behaviour of the plume. For a 6 m/s wind velocity, increasing the exhaust temperature from 150 °C to 200 °C and significantly limiting the downward dragging of the plume due to the swirling motion. On the other hand, if the exhaust gas temperature rises to 250 °C, the plume cannot drift towards the hillside, resulting in a much higher plume lift. These changes in the plume behaviour contribute to a more effective dilution of the pollutants when approaching ground level.

Also, for a wind velocity of 10 m/s, a gas exit temperature of 200 °C is enough to prevent the plume from collapsing towards the ground of the industrial facility, allowing the plume to expand in a north-easterly direction at a sufficiently safe height from rooftops of nearby buildings. Increasing the temperature up to 250 °C further improves the dilution mechanisms but it is not necessary to reach such a high temperature to keep the contaminants away from critical areas.

Running CFD simulations for different scenarios and operating conditions has enabled the identification of critical situations in terms of pollutant concentration conditions. Moreover, this tool has helped to define a mitigation strategy, offering a good balance between impact on air quality and operating conditions, the latter determining the costs of the fuel consumption in the plant. The analysis conducted here can be extended, adding other meteorological conditions to cover as many scenarios as possible. This would allow for the adjustment of the emission parameters in most situations.

5. Conclusions

This research has addressed the challenge of simulating the dispersion of pollutants in a case study which, due to its complex nature in terms of flow phenomena, must be tackled using advanced modelling techniques such as CFD.

The complexity and singularity of the selected case study can be summarised as an elevated industrial source located on steep terrain, with mountains in close proximity to the source with elevations well above the height of the chimney. Furthermore, given the location of the source in relation to the hill and the prevailing winds in the area, significant flow recirculation phenomena can be expected in the area.

To address this problem, an original methodology was developed to provide a robust computational domain for cases of complex scenarios such as the one addressed in this research. This methodology relies on the use of LiDAR elevation data (although other formats of digital terrain elevation models are also valid) to shape the topographic surface of the area, as well as urban planning data to integrate buildings and other civil infrastructure that may impact the flow patterns. This methodology attempts to achieve a clean and smoothed geometry domain, which is one of the main challenges when modelling large areas of terrain with steep characteristics and the presence of multiple urban obstacles.

Failure to generate a sufficiently robust geometry for the CFD model can lead to serious problems in the mesh generation process, compromising the feasibility and reliability of the simulations. In order to facilitate the grid creation process, this research proposes a strategy to generate the geometry domain as a single volume, which consists exclusively of the fluid space encompassing the study area. In addition, not only the most common and traditional guidelines for geometry and mesh generation were followed but other practices identified in the recent bibliography were also utilised to achieve an appropriate geometry and grid design.

With the CFD model in place, the simulations of pollutant dispersion under different meteorological regimes have revealed the detailed plume path in the vicinity of the source. The formation of flow recirculation zones affecting the exhaust section of the chimney were confirmed in the simulations. Also, SO₂ iso-concentration surfaces of different values were computed, revealing unexpected deviations in the path of the plume before it is diluted in the air, while being carried away from the stack section.

On the other hand, the versatility of the CFD models has made it possible to verify the effectiveness of a pollution mitigation measure based on the stack exit temperature, which is obtained via the regulation of the waste heat recovery systems of the plant. In the analysis performed, the effect of three different gas temperatures on the plume dispersion was considered. The simulation results indicate that, under certain meteorological conditions, the plume exhibits changes in its trajectory that cause the pollutants to approach the mountainside to the southwest of the chimney, where some small houses could be exposed to high concentration levels. In this situation, the simulation results also show that increasing the gas exit temperature is useful to keep the plume away from the inhabited area of the mountain during the high exposure events identified in the previous step.

Therefore, the main achievements of the research were a methodology to generate robust computational domains for complex terrain scenarios and the use of the CFD model to assess the effectiveness of a pollution mitigation measure in different weather conditions. These results can be seen as a tool to support the decision-making process in finding a compromise between the environmental pollution impact and the fuel consumption of the plant.

As a continuation of this research, certain future tasks have been identified:

• Testing the described methodology in different locations.
• Carrying out simulations in stratified, stable and unstable atmospheric conditions.
• Comparing the results of the CFD model with those of an operational or fast response formulation.