The location chosen for this study is a neighborhood in Niigata city in Japan, a building complex with simple building shapes. This location is used in the Architectural Institute of Japan (AIJ) benchmark case E, and its geometric model and data from wind tunnel experiments are provided at Architectural Institute of Japan, (2021). Other published CFD studies on this location available in the literature include (Tominaga et al., 2005; Mochida and Lun, 2008; Caniot et al., 2011; Tominaga, 2012).

The BPGs have been used as a minimum requirement when sizing the computational domain. The guidelines are based on box-shaped domains. They are only suitable if wind at the inlet has one specific direction, which is inward parallel to the domain’s sidewalls.

A cylindrical domain is used as the intention of this study was to have a flexible mesh suitable for simulating wind from any horizontal direction. Therefore, the recommendations of the BPGs had to be applied and adapted to the cylindrical domain.

The routine of using a cylindrical computational domain when performing urban wind simulations was introduced by Kastner and Dogan, (2020). They reported that a cylindrical domain yielded comparable results to a traditional box-shaped domain in terms of accuracy and convergence behavior. Figure 1 presents the computational domain, and the dimensions are listed on the right side of the figure.

The BPGs suggest that the top of the domain should be 5 H _max from the top of the tallest building. H _max, as defined by Franke, (2006), is the height of the tallest building and is essential in sizing the computational domain. Here, H _max is 60 m. As there is no elevation in the terrain patch, the height of the domain is set to 6 H _max. Franke et al. (2007) suggest a distance of 15 H _max downstream of the buildings to the outlet of the domain. As the cylindrical domain is suitable for wind coming in from any direction, the radial distance from the buildings to the domain’s sides is set to 15 H _max.

The computational grid was developed according to the BPGs, and Figure 2 presents a partial view of the grid structure from the side. The final grid is produced with snappyHexMesh, the hex-dominant unstructured mesh generator of OpenFOAM. snappyHexMesh requires a base mesh and 3D geometry in the form of a triangulated surface, here represented by the building model and the terrain model. Hågbo et al. (2019) concluded that an unstructured mesh generated with snappyHexMesh could be sufficient in urban wind simulations where only the general flow features are of interest.

The base mesh was generated with blockMesh, a meshing utility in OpenFOAM suitable for simple structured meshes with grading and curved edges. The grid is stretched vertically in the base mesh with an overall expansion ratio of three. It means that the height of the top cells is three times the height of the bottom cells.

The grid is refined successively in a 2:1 ratio near the terrain and buildings in the final mesh, with eight grid cells between each refinement level. Two refinement levels are used on the outermost part of the terrain patch, while three levels are used within a radius of 500 m from the center of the patch. Four refinement levels are applied near the buildings.

The grid consists of about 12.1 million cells, and key features of the mesh are listed in Table 1.

The height of the cells adjacent to the buildings and on the terrain near the buildings is 0.34 m. On the rest of the terrain within the built area, the cell height is 0.68  m, ensuring the elevation height of 2 m to be within a minimum of the third cell from the ground. These cell heights give the average non-dimensional wall distances y ⁺ value of about 7.0 × 10³ for the building patch and 9.0 × 10³ for the terrain patch in the built region.

These are relatively high values that significantly exceed the recommended values of (30–500) (Franke et al., 2007), meaning that the first cell center is well outside the log-law region. However, standard logarithmic wall functions are typically also used in CFD simulations of atmospheric boundary layer wind flow when y ⁺ is well above the recommended upper limit without reducing the accuracy of the simulated wind velocity field. Examples using y ⁺ values of the same order of magnitude as in this study include the studies cited herein (Blocken et al., 2007; van Hooff and Blocken, 2010a; van Hooff and Blocken, 2010b; van Hooff et al., 2011).

The high Reynolds number flow studied here makes it challenging to be within the log-law region. The most important reason for using these high y ⁺ values is that the recommended range is based on smooth walls and would yield unnecessarily small near-wall cells for our rough terrain wall. It should also be noted that the log-law profile is not necessarily valid given the complex geometry of the terrain.

Apart from the y ⁺-values, the resolution is sufficiently high according to the BPGs, as it satisfies the recommendation of a minimum of 10 cells per building side and in the passage between buildings, as well as at least 3 cells from the ground to the elevation height (Franke et al., 2007; Yoshie et al., 2007; Tominaga et al., 2008). In this study, the elevation height was set to 2 m above ground level, a reasonable height for pedestrian wind comfort evaluation.

An outline of the simulation setup applied is provided at Table 2. The wind simulations were performed by solving the incompressible, three-dimensional steady Reynolds-averaged Navier-Stokes (RANS) equations with the finite volume method. Turbulent closure was provided by the realizable k-ϵ turbulence model developed by Shih et al. (1995) and the SIMPLEC algorithm handled the pressure-velocity coupling ensuring that the continuity equation is satisfied. The realizable k-ϵ turbulence model was chosen, as it is recommended by Blocken et al. (2003) and Blocken et al. (2004) in the application of pedestrian wind environment simulations, as it proved to perform substantially better than the standard k-ϵ model. Second-order discretization schemes were used for the spatial discretization, including the convective terms.

The iterative convergence criteria were that the scaled residuals should drop four orders of magnitude for k (turbulent kinetic energy), ϵ (turbulent kinetic energy dissipation rate), and each Cartesian component of U(velocity), as recommended for all flow variables by the BPGs. The convergence criterion for p (pressure) was relaxed to 5.0 × 10^–4. This was done to speed up the simulations, as reaching a residual drop of four orders of magnitude required significantly longer simulation time. The relaxation of the convergence criterion did not impact the wind velocity magnitude at the pedestrian level in the built area investigated in this paper. These criteria were met for all simulations presented.

The simulations were conducted using the OpenFOAM (version 6) toolbox, which is an open-source CFD software package (Weller et al., 1998).

An overview of the boundary conditions (BCs) used is presented in Table 3. Combined inlet/outlet boundary conditions are used on the computational domain’s sides to allow for wind in any direction. The inlet profiles were set to form a homogeneous atmospheric boundary layer (ABL) with neutral stratified conditions and a logarithmic velocity profile, according to the BPGs. The wind speed at a reference height of 10 m, U _ref10m , was set to 5 m s⁻¹. Outlet conditions were used with constant pressure and zero gradient for the remaining variables. The direction of the flux automatically determines the inlet and outlet regions.

No-slip conditions with wall functions were used for both the terrain and the buildings layer. The logarithmic law for smooth walls was used for the buildings. A rough wall condition was applied for the terrain (Hargreaves and Wright, 2007). The roughness was set to z ₀ = 3.00 m to replicate the profiles used in the experiment. The slip condition was applied to the top boundaries.

The instantaneous regional wind conditions at any urban site continuously change with time. They include quantities such as wind velocity and direction, turbulence, humidity, and temperature, and how these parameters vary with height. The varying conditions result from changing external weather patterns, which depend on other factors such as time of day, season, atmospheric stability, formation of low- and high-pressure systems, and many more.

The pedestrian wind comfort maps are produced by combining wind simulations, local statistical meteorological data in the form of a wind rose, and a set of defined comfort criteria. There are numerous sets of wind comfort criteria, and in this study, a version of the Lawson wind comfort criteria (Lawson, 1978) is used, see Table 4. The Lawson criteria classify areas with “comfortable” activity types based on their associated accepted average wind speed intervals and exceedance frequency. The calculation of wind comfort criteria was also considered in Blocken and Carmeliet, (2004), Blocken et al. (2012) and Vita et al. (2020a). Note that only comfort is considered in this work, not pedestrian distress (Vita et al., 2020b).

The number of simulated wind directions used in the production of the pedestrian wind comfort maps in this study is:

1) 4, increments of 90°

Not all wind speeds are simulated. Instead, a representative wind field for each Lawson category is found by a linear scaling of the simulated wind velocity. The maximum value in each wind speed range is used as the scaling factor to be conservative. Linear scaling is justified as the wind is well within the turbulent flow regime for all the wind speed intervals under consideration, and hence similar flow distributions are expected.

The procedure of producing pedestrian wind comfort maps in this study is as follows. First, the wind is simulated from 64 evenly distributed horizontal wind directions. A set of 4, 8, 16, 32, or all 64 simulations of evenly distributed wind directions are used further.

Then, the pedestrian wind comfort is evaluated in all the computational cells using the algorithm outlined in Figure 3.

Since the Lawson criteria are based on exceedance probability, it needs to be determined how often the local wind velocity in each computational cell is above the given criteria. Since we assume equal probabilities for each wind direction and each wind class, the frequency for each wind class is P class = 1 n sector ⋅ n class (1)

In general, this frequency would be obtained from the wind rose for the location. Next, we need to find the local velocity of the cell for each wind direction and wind class. We only have a single simulation for each wind direction. For each wind class, we calculate a representative local velocity in each cell for this wind class by doing a linear scaling the simulated wind velocity, U local = U class U ref U simulated (2)

Here, U _ref is the reference velocity at 10 m height for the simulated case. For U _class, we use the upper bound of the class to be conservative. If this local velocity for the wind class exceeds any of the comfort criteria, we accumulate the frequency in this cell by P _class. After all the wind directions and wind classes have been considered, we determine if the total accumulated frequency has been exceeded for each category in the Lawson criteria, and assign the appropriate comfort value to the cell.

The algorithm assigns a comfort class to all the cells within the computational domain. Finally, a pedestrian wind comfort map is produced by reading the comfort of the cells located 2 m above the terrain level and projecting its representing color onto the terrain patch below.

In this case study, the local wind statistics were arbitrarily chosen such that the wind was to always be 5 m s⁻¹ at 10 m height and that the frequency of each wind sector was identical. In that way, the importance of each wind simulation is equally weighed in the pedestrian comfort maps.

The wind tunnel experiments were conducted at a 1/250 scale (Yoshie et al., 2007). A boundary layer with properties similar to an atmospheric boundary layer was generated in the wind tunnel upstream of the test section of dimensions 1.8 m length, 1.8 m width, and 1.8 m height. The vertical profile of the horizontal wind velocity was developed to follow a function with a power-law exponent of 0.25 when reaching the test section. Scalar wind velocities at 0.008 m above the wind tunnel floor (2 m above the ground surface in real scale) were measured by multi-point thermistor anemometers.

Figure 4 shows the normalized, simulated wind velocity magnitude at 2 m above ground, with the wind from N (0°). The colored dots are located in the same position as the 80 probes in the wind tunnel experiments, and their color represents the wind velocity magnitude difference, ΔU _val, between the wind tunnel experiment and the simulation, where ΔU _val is: Δ U val = U exp − U sim (3)

Blue-colored dots indicate that the wind velocity magnitude of the simulation is greater than in the wind tunnel experiments, and red-colored dots mean that the wind velocity magnitude of the simulation is lower than in the wind tunnel experiments. The experimental values compared against are averaged over time.

Figure 4 clearly shows the expected speed-up effect at the ground level near the side corners of the high-rise building relative to the overall wind direction. The similarly expected wake regions of low wind velocity magnitudes are seen leeward of the buildings.

When comparing the wind velocity magnitudes at the probes in the simulations to the results from the wind tunnel experiments from 16 wind directions, we see a slight trend that the simulations are overpredicting the effects of the buildings. More specifically, the simulated values in relatively high wind speed areas are generally too high compared to wind tunnel results. Likewise, the simulated values are too low compared to wind tunnel results in low wind speed regions. This trend is slightly noticeable in Figure 4, where red dots are somewhat located in blue low-velocity areas, and blue dots in red high-velocity areas. Note that the calculated differences, represented by the colored dots, are highly sensitive to the location of the probes as they are positioned within regions of high flow gradients.

Similar observations were made on the same benchmark case by Tominaga et al. (2005), Mochida and Lun, (2008), and Caniot et al. (2011). All three studies reported that the CFD simulations generally agree with the wind tunnel experiments, but that the simulations tended to underestimate the wind speed in the wake regions of buildings. Tominaga et al. (2005) note that the disagreement is also seen for other comparable benchmark tests and that it is partly related to differences in the definition of the average wind velocity magnitude in the two methods. While the wind tunnel tests register time averages of the instantaneous scalar velocity, the simulations use the velocity vector’s time averages (Tominaga et al., 2005). Mochida and Lun (2008) argue that it is related to the used turbulence models’ inability to reproduce the vortex shedding from tall buildings. Caniot et al. (2011) emphasize that the points of high relative difference are positioned near the edges of the wakes where high flow gradients occur, making the values highly location-sensitive.

Table 5 lists the average relative difference of wind velocity magnitude, and the Pearson correlation coefficient (PCC) between wind simulations and wind tunnel experiments for 16 wind directions at the 80 probes. The experimental values compared against the simulations are averaged over time. The table quantifies that the wind velocity magnitude of the simulations is overall lower at the ground level than in the wind tunnel experiments. For the 16 wind directions, the average relative difference of wind velocity magnitude is in the range from −1% to −35%, and the Pearson correlation coefficient is in the range 0.71–0.84.

The overprediction of the extent and significance of the wake regions of low wind velocity magnitudes, leeward of the buildings, agrees with other RANS simulations in the literature. Generally, it is because the transient nature of the flow is difficult to fully represent in time-averaged simulations, leading to less turbulent mixing in the simulations. It is also evident in Table 5 where the negative values of the average relative difference of wind velocity magnitude between wind simulations and wind tunnel experiments are because many of the probes are located within the wake regions in the simulations.

Figure 5 presents the wind velocity magnitude at 2 m above ground for 16 simulations representing wind from a 90° wide sector around the reference wind direction N (0°). Notice that there are generally always high wind speeds in the vicinity of any two of the corners in the high-rise building, and which of the corners is affected is highly dependent on the wind direction. For example, the high wind speed corners in the 349° simulation are the West and the North corners, but are the West and the East corners for the 0° simulation, where the wind direction is only shifted by 11°. Similar jumps of high wind velocity areas occur throughout all wind sectors.

The variance within the set is best viewed in Figure 6. It illustrates the wind velocity magnitude difference within the North sector, ΔU _Nsec, to reference wind direction simulation at 2 m above ground, where ΔU _Nsec is: Δ U Nsec = U ref0 ° − U 321 ° − 45 ° (4)

Negative values, represented by blue, correspond to areas where the reference case has a lower wind velocity magnitude than the compared case and vice versa. It is clear that the differences to the reference wind direction simulation grow with increasing wind direction shift from it, in terms of both magnitude and spatial extent. The differences are especially pronounced near the corners of the high-rise building. Notice that the difference is represented by alternating signs/colors going around the high-rise building. This alternation is related to the high sensitivity of the location of the “high wind speed corners,” seen in Figure 5. Quantification of the variances seen in Figure 6 is provided in Table 6.

Each wind sector is equally spaced. The width of the sectors depends on the number of wind directions included in the production of the pedestrian wind comfort maps. Therefore, each wind simulation represents a 90° wide sector when 4 wind directions are used. Wind coming from the North sector ranging from 315° (NW) to 45° (NE) is represented by only one wind simulation of wind coming from 0° (N).

Table 6 shows that out of the 16 wind directions simulated in this sector, about a minimum of 88.6% of the area has an absolute value difference in wind velocity magnitude of less than 0.5 m s⁻¹ in the wind simulation with the highest discrepancy to the reference wind simulation. The remaining 11.4% is the maximum area in the reference wind simulation that is more than 0.5 m s⁻¹ off compared to one of the wind direction simulations it represents in the 90° wide sector.

Doubling the number of wind simulations used in the pedestrian wind comfort maps from 4 to 8 reduces the sector width each wind simulation represents to 45°. It reduces the area to 8.5%. Likewise, 16 wind directions further reduce the wind sector increments to 22.5° and the area to 4.4%; 32 wind directions give increments of 11.25° and an area of 1.0%.

Specially detailed information of the wind velocity magnitude distribution at pedestrian-level obtained from wind simulations is essential when producing the pedestrian wind comfort maps. The number of simulated wind directions decides the equally spaced wind sector increments each wind simulation represents. Fewer wind directions causing larger increments are naturally associated with higher variance in wind velocity distribution within this wind sector. This relation is evident in Figures 5, 6 and Table 6.

Figure 7 shows the pedestrian wind comfort maps at 2 m above ground produced using either 4, 8, 16, 32, or 64 evenly distributed horizontal wind directions. Only a section of the computational domain, 100 m in each horizontal direction of the high-rise building is presented.

Ground areas in the vicinity of the high-rise building, especially near the corners, are classified with the worst (“uncomfortable”) or the second-worst (“business walking”) wind comfort classification in all the maps. In these areas, the frequent strong wind is caused by a wind pull-down effect often seen near high buildings, particularly ones with sharp corners. The wind is considerably weaker at the ground level close to the other lower buildings. The areas with frequent high wind are more evenly distributed around the high-rise building, as more wind directions are included in the production of the maps, especially compared to the map produced with wind simulations from only 4 wind directions.

The differences in the pedestrian wind comfort maps are best viewed in Figure 8. The figure gives the classification difference to the reference case, ΔC, where: Δ C = C i − C j (5)

Ci is the pedestrian wind comfort class of a given location using 64 wind directions. Cj is the classification at the same location using a different number of wind directions. The map produced from the most wind directions, here 64, is likely the most accurate and is used as the reference case.

Negative values, represented by blue, correspond to areas where the reference case has a lower classification number, meaning an improved wind comfort class, than the compared case and vice versa. Notice that there is a relatively large total area in red in Figure 8A. The reference case has a higher classification number in these red areas, meaning a worse wind comfort class, than the case using only 4 wind directions.

Quantification of the differences is provided in Table 7. It shows the percentage of area classified differently compared to maps using all the 64 simulated wind directions. The total area classified differently from the reference map clearly reduces in size as more wind directions are included; 79% of the area within 100 m in each horizontal direction of the high-rise building is classified with the same pedestrian wind comfort class as the map produced using 64 wind directions. This number increases to 92% using 8 wind directions, 96% using 16 wind directions, and finally 99% using 32 wind directions. Similar values were observed when including the whole building area of radius 200 m.