ANSYS FLUENT tools were used as part of a standard turbulence k-ε model to simulate the air flow around a number of typical obstacles (a solid cube, a solid hemisphere, and a 2D hill) which form a potential terrain in the NPP emission dispersion area and roughly correspond to the geometry of the buildings and structures within this area. For reproducibility, a non-uniform spatial grid is plotted in the computational region which condenses near the obstacle surface and the outer boundaries. The dimensions and the positions of the obstacles were chosen such that to ensure their best possible coincidence with the conditions of the published experiments. The result of simulating the velocity and direction of the air flow as the whole shows a good agreement with the data from the wind tunnel experiments in the areas in front of and over the obstacle, as well as in its air shadow. Typical accelerated flow, vortex, and reverse flow areas are reproduced reliably. There are variances observed only in the local heavy turbulence regions in the obstacle’s air shadow near the ground surface. All this indicates that it is possible to model in full scale the dispersion of the NPP emissions taking into account the peculiarities of the plant site terrain and the major onsite structures to determine more accurately the personnel and public exposure dose.

NPP gas and aerosol emissions, turbulent diffusion modeling, ANSYS FLUENT

The dispersion of the NPP gas and aerosol emissions depends primarily on the wind direction and speed. The diffusion of emitted products, transversely with respect to the wind, is connected primarily with natural fluctuations of the wind direction, turbulent mixing of air masses caused by the state of the atmosphere (upward flows of warm air), and their own viscosity and friction on the underlying surface. The dispersion of the emission in conditions of a complex terrain (hills, ravines, etc.), as well as in the presence of buildings and structures leads to an additional source of turbulence as the result of the air flowing around such obstacles. In this case, one can hardly expect that it is possible to estimate the concentration of radioactive substances in the air based on simple Gaussian models recommended by the IAEA (Safety Series No. 50-SG-S3 1980).

Among a variety of methods to investigate the air flow in the atmosphere, such as wind tunnel experiments and full-scale ground measurements, extensive use is made of computational fluid dynamics (CFD) (Leelőssy et al. 2018; Yoshihide et al. 2008; Gorle et al. 2009; Ai and Mak 2013), which has proved to be a reliable and cost-effective tool of simulation for turbulent flow problems capable to identify the peculiarities of turbulent flows around obstacles including air shadows, rupture zones, and vortex paths (Stupin and Overko 2006).

The ANSYS FLUENT package tools (ANSYS Fluent Theory Guide 2013) were used in the study to simulate the air flow around a number of typical obstacles (a solid cube (Yu et al. 2017), a solid hemisphere (Zhenqing et al. 2019; Takeshi et al. 1999), a 2D hill (Ferreira et al. 1991; Kim et al. 1997)), which form the potential terrain of the NPP emission dispersion area, and roughly corresponds to the geometry of the buildings and structures within this area. The shapes and dimensions of the obstacles were chosen based on the available publications on wind tunnel experiments (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010), which allowed comparing the numerical simulation results with direct measurement data and determining so the extent to which the computational model used was reliable. The actual purpose of the study is to make certain that it is possible to simulate, in an adequate way, the dispersion of gas and aerosol emissions in conditions of a complex terrain to allow calculating the actual deformation of the radionuclide concentration fields in the subsurface air layer for the NPPs deployed in the respective terrain.

Reynolds-averaged Navier-Stokes (RANS) equations are used as the computational model in the framework of the standard k-e-model of turbulence (Juretic and Hrvoje 2013; Richards and Norris 2011). The latter has been broadly and successfully employed for problems of impurity spreading in the atmosphere (Xing et al. 2013; Kiša and Jelemenský 2009). The key parameters of the model are the kinetic turbulence energy k and the turbulence energy dissipation rate e, which are specified in this paper through the friction velocity u^* as

$\begin{array}{l}k=\frac{{\left(u^{*}\right)}^2}{\sqrt{C_{\mu }}}\\ \epsilon =\frac{{\left(u^{*}\right)}^3}{ky}\end{array}$

where C_m = 0.09 is the dimensionless empirical constant, and y is the vertical coordinate.

The computational region is defined by the type of the obstacle under consideration.

The obstacle in the form of a 2D hill (Fig. 1) is given parametrically as

Figure 1. 

Diagram of the computational region for an obstacle in the form of a 2D hill.

$x=\frac{1}{2}\xi \left[1+\frac{a^2}{{\xi }^2+m^2\left(a^2-{\xi }^2\right)}\right]$ , (1)

$y=\frac{1}{2}m\sqrt{a^2-{\xi }^2}\left[1-\frac{a^2}{{\xi }^2+m^2\left(a^2-{\xi }^2\right)}\right]$ , (2)

where x is the coordinate along the wind direction, and y is the vertical coordinate.

The hill foot dimension is equal to 2a. The parameter x used in formulas (1) and (2) is variable in the limits of |x| £ a. The value m is determined through the average hill slope n = H/a as m = n + (n² + 1)^1/2. For the calculations, the hill height H is assumed to be equal to 0.117 m, and the average slope angle is assumed to be 26°.

The problems of solid hemisphere- and cube-type obstacles flown around were handled using rectangular regions the dimensions of which with the positions of the obstacles relative to the inlet plane are shown in Table 1.

For the reproducibility of the results, a nonuniform spatial mesh was given in the computational region which becomes denser near the obstacle surface and the outside boundaries. Fig. 2 presents an example of such hexahedral mesh for a cubic obstacle.

Figure 2. 

A 3D hexahedral computational mesh for a cubic obstacle.

The law of the wind inlet speed variation with the height u (y) was defined differently for different obstacles to ensure the comparability with the experimental data.

For the hemisphere, the wind speed profile was determined by the power function similar to that used in (Tavakol et al. 2010):

$\frac{u}{U_0}={\left(\frac{y}{H}\right)}^n$ , (3)

where n = 0.135, and U₀ corresponds to the wind speed at the obstacle height H.

For the 2D hill, following the recommendations in (Richards and Hoxey 1993), a logarithmic speed profile was used which is close to that detected in wind tunnel experiments:

$u\left(y\right)=\frac{u^{*}}{K}\left[\ln \left(\frac{y+y_0}{y_0}\right)\right]$ , (4)

where u^* is the friction velocity; y₀ is the roughness height in m; and K is the Karman constant which was assumed to be equal to 0.41. And the free flow velocity value of U₀ = 4 m/s was achieved on the upper plane of the region under consideration.

For the calculations of the flow around the cube, the inlet wind speed was assumed to be constant with the height, that is, u (y) = U₀ where U₀ = 0.6 m/s.

The actual investigation results were calculations of the wind speed values and directions near obstacles of different forms as compared with the wind tunnel experiment data (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010). For illustration, Figs 3 through 5 present the vertical profiles of the flow velocity longitudinal component (along the x axis) in front of the obstacle and over and behind it on the obstacle symmetry axis (y = 0) for the x-coordinates given in Table 2. The estimated velocity at a particular height u (y) was reduced to the value of the inlet flow velocity U₀. The origin of the x axis coincides with the obstacle center.

Figure 3. 

Vertical profiles of the wind speed longitudinal component in front of the obstacles: solid line – calculation; diamonds – experimental data (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010).

The flow distortion is minor in front of the obstacle with the distances from same being of the order of the obstacle height as such (Fig. 3). The reduced speed deviation from unity is explained by the influence of the computational region’s upper and lower boundaries. For the cube, the boundary conditions for both planes define the zero velocity (the wall condition). And the region under consideration is only twice as high as the obstacle. For the latter, however, the flow deceleration in the area y < H is noticeable as the distribution of the considered and measured velocities is asymmetrical relative to the midpoint y = H. Symmetry conditions were set for the calculations of the hill and the hemishpere flown around on the upper boundary due to which there is no deceleration on it. There is a low back flow of u(y) < 0 near the very ground, caused by the swirling in the lower part of the obstacle, found in the calculations and experiments for the hemisphere. On the whole, the coincidence of the estimated and experimental data can be regarded as fairly sastifactory.

The profiles of the wind speed longitudinal component obtained immediately over the obstacle center are shown in Fig. 4. These also exhibit a good agreement between the calculations and the experiments. Besides, both the calculations and the experiments show an extensive vortex and back flow area on the upper face of the obstacle least flown around (cube). And there is also an accelerated flow observed over this area.

Figure 4. 

Vertical profiles of the wind speed longitudinal component over the obstacles: solid line – calculation; diamonds – experimental data (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010).

The distribution of velocities behind the obstacle and, especially, in its air shadow (Fig. 5) is characterized by the flow deceleration and swirling in the reverse direction – it is the larger, the less flown around the obstacle is. This effect is most profound for the cube. On the whole, simulation gives an overestimation of the velocity in the back flow area as compared with the experiment. The maximum differences in the absolute velocity values are observed in the area with the highest turbulence (immediately behind the obstacle near the ground surface). Such overestimation has been most likely caused by an assumption of isotropic viscosity in the numerical model which is not exact, specifically near the ground surface where the intensity of turbulence varies greatly in space. Unfortunately, the sources we used (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010) do not contain data on the experimental data uncertainties.

Figure 5. 

Vertical profiles of the wind speed longitudinal component behind the obstacles: solid line – calculation; diamonds – experimental data (Trombetti et al. 1991; Martinuzzi and Tropea 1993; Tavakol et al. 2010).

The specific properties of the flow around the cubic obstacle are clearly seen in Fig. 6 which presents the distribution of the flow velocity vectors. The vertical lines in the figure correspond to the x-coordinates in which the velocity profiles for the cube shown in Figs 3 through 5 were calculated.

The cube that simulates the buildings and structures of the NPP as such is the obstacle most difficult to flow around. The downwind length of the turbulence area (see Fig. 6) is the maximum one as compared with the hill and the hemisphere. Unlike these, specific to the cube is also the formation of a turbulent area with the flow breakdown over the roof’s upper plane, as well as along the vertical side planes where a pair of oppositely whirling vortices form in the horizontal plane. The flow starts to break down in front of the structure and develops further from the latter’s upper front face and on the side walls. Due to these peculiarities, the problem is of a special interest in terms of testing current CFD calculation algorithms.

Figure 6. 

Flow velocity field for the cubic obstacle flown around.

Comparing the simulation results for the wind flow around typical irregularities of terrain against direct wind tunnels experiments makes it possible to conclude that models based on RANS equations, despite a moderate computational effort, provide for a substantially satisfactory coincidence with the experiment. The exceptions are local areas of heavy turbulence in the air shadow near the ground surface and immediately over the horizontal roofs of the buildings.

The presented comparison gives ground to believe that simulating the dispersion of the NPP gas and aerosol emissions in conditions of an irregular terrain and in the presence of buildings and structures, based on RANS equations and using the ANSYS FLUENT package, is capable to yield adequate results for calculating exposure doses for personnel and locally residing population.