Simulation of a ground-mounted prism in ABL flow using LES: on overview of error metrics and distribution

The accuracy of wind loading predictions using Large Eddy Simulation (LES) is usually influenced by numerous model parameters, which can influence the obtained results. The validation of numerical simulations with traditional Wind Tunnel Test (WTT) is still an important task, necessary to increase our a priori knowledge of possible inaccuracies and set up mitigation strategies. In this study, LES is used to simulate the wind fields around an isolated model high-rise building, under seven wind attack angles and validated with WTT results. The influence of various settings and parameters on the model performance is studied. For the angle of attack showing higher inaccuracy, different mesh refinement strategies and turbulence models are tested. Results indicate that LES can accurately predict the mean and local maxima of the pressure coefficients for both perpendicular and skew wind attack angles, as well as reproduce global forces and their envelopes with very good accuracy. Conversely, pronounced errors are found in the prediction of the pressure coefficient standard deviation and the local minima. The highest deviations between LES and WTT are found close to the leading edge in correspondence of flow separations which are observed in WTT and not in LES for skew flows. The addition of boundary layer cells and the use of different subgrid models have very limited effectiveness in modifying the obtained results for the analysed case.


Introduction
Large Eddy Simulation, LES, is known to be a powerful tool for the investigation of wind effects on structures such as high-rise buildings [1][2][3][4]. However, comparison with traditional Wind Tunnel Test (WTT) sometimes shows discrepancies, so that their validation and the individuation of critical aspects which might affect results is still an important research activity. Numerous studies can be found in the literature validating LES results in terms of wake flow characteristics [5,6], global aerodynamic forces [7,8] and pressure distribution [9][10][11]. The usefulness of such validations mainly consists in the individuation of critical aspects which can be sometimes moderated by appropriate modelling choices. When moderation strategies are not found or lead to excessive computational Xing Advances in Aerodynamics (2023) 5:9 burden, their individuation anyway allows to build up knowledge, which can help in correctly interpreting results obtained by numerical models when WTT results are not available.
From the studies above and analogous works [12,13], it can be found that the model performance of LES is usually good but, in some cases, sensitive to various model parameters. Among them, the accurate reproduction of flow characteristics in the Atmospheric Boundary Layer (ABL) is known to be crucial for the evaluation of wind effects on structures [14,15]. Nevertheless, numerous other factors can influence the simulation [16][17][18], which include, as pointed out in [19], ground roughness, subgrid scale model, domain size, near-wall treatment and sensitivity to mesh and timestep.
Numerous contributions have been devoted to investigating the influence of some of the aforementioned model parameters on the model accuracy. Unfortunately, due to the variability of the adopted computational setup and analysed flows, contradictory results are sometimes found.
In particular, Wang et al. [20] assessed the effects of the mesh type and boundary mesh on the time-averaged and fluctuating wind characteristics around an isolated high-rise building standing in ABL. They found that the case with hexahedral cells has the best agreement with experiments. They also found that the boundary layer mesh does not improve the numerical accuracy in any circumstance, despite the fact that the nondimensional wall distance ( y + ) is often considered to be a very important quantity to be controlled. Wijesooriya et al. [21] showed that the choice of the Sub-Grid Scale (SGS) model does impact the flow field. The wall adapting local eddy viscosity (WALE) model was found to be the most suitable among the tested models including Smagorinsky-Lilly model, Dynamic model, kinetic energy sub grid scale model and wall modeled LES. Liu et al. [22] examined the effect of mesh resolution and different SGS models on the prediction of the pedestrian level wind fields around building arrays. They observed that the dynamic Smagorinsky-Lilly method gave the best model performance and recommended appropriate values for time step and sampling period.
The aim of the present paper is to assess the numerical accuracy of LES in predicting wind loads on a high-rise building and to study the influence of some important model parameters on the accuracy of the obtained results. The high-rise building model is an isolated flat-topped box with an aspect ratio of 1:3:5, standing in an ABL. It is one of the high-rise building models from the public wind tunnel database constructed by Tokyo Polytechnic University (TPU), Japan, which have been considered also in some of the already available contributions [19,[23][24][25].
Firstly, the numerical results are compared with WTT data in the form of several model performance metrics for 7 angles of attack ranging from 0 • to 90 • . Then, for the worst case, corresponding to a skew flow impinging at 45 • , the spatial distribution of the pressure coefficient statistics (i.e. time-average, standard deviation and peak values) and the corresponding Prediction Error (PE), defined as the difference between numerical and experimental results, are analysed. For such worst case, different numerical setups are evaluated varying mesh configuration and testing three SGS models, namely, WALE, k − ω SSTSAS and Smagorinsky.
The paper is organized as follows. The descriptions of WTT and LES are provided in the next section, along with the procedure used to generate inflow conditions. Then, the numerical results are validated in Section 3. Subsequently, the sensitivity studies of the PE to the mesh configuration and the SGS model are provided in Section 4. Finally, conclusions are drawn in Section 5.

WTT description
The WTT were performed in the Boundary Layer Wind Tunnel (BLWT) of TPU, Japan.
It is an open circuit tunnel with 2.2 m width and 1.8 m height. The ABL in the experiments was generated through spires and square blocks acting as roughness elements. It corresponds to a terrain category IV in the Architectural Institute of Japan (AIJ) standard [26]. The wind speed at building height is around 11.11 m/s, and the stream-wise turbulence intensity ( I u_WTT ) is approximately 15%.
The tested model is a flat-topped prismatic box with 0.5 m height (H) and 0.3 m (L) * 0.1 m (W) rectangular cross-section. The length scale reduction factor of the experiments is 1:400, leading to a building height equal to 200 m in real scale. 480 wind pressure taps are positioned on the four side faces, sampling synchronously with a frequency of 1000 Hz for a duration of 32.8 s. The origin of the coordinate system is the center of the building base. The geometry of the building is presented in Fig. 1 (a).
Incidence angles ranging from 0 • to 100 • with an increment of 5 • are considered in WTT. The 0 • is defined as the condition in which the flow is moving along the +x direction and it is perpendicularly hitting the long building side, as indicated in Fig. 1 (b). The complete database and detailed description of the WTT are publicly accessible [27].

LES description
LES are designed to reproduce the WTT using the open source software OpenFoam v6.0. A 3D view of the computational domain is presented in Fig. 2 (a). The domain has the dimensions of 3 m, 2.2 m and 1.8 m in stream-wise, lateral and vertical directions, respectively. The distance between the inlet patch and the building location is 1 m, which equals the double of the building height H. The blockage ratio of the simulations is approximately 3.0% for 0 • and is 1.0% for 90 • , which are both lower than recommendations [28]. The turbulent flow at the inflow patch is generated using a synthetic turbulence generation technique called PRFG 3 . It is an extended version of the Prescribed-wavevector Random Flow Generator (PRFG) [29], aiming at giving control over three-dimensional spectral densities, in other words, explicitly targeting all integral length and time scales of turbulence. Interested readers are invited to the references [29,30] for details regarding the method. The top and two sides of the domain are treated as symmetry conditions, mimicking the effect of wind tunnel walls. The domain bottom and the building surfaces are treated as walls. In particular, the ground of the domain is modeled as a rough wall with a roughness height of 0.015 m (model scale) and a roughness constant of 0.5. The building surfaces are modeled as smooth walls. An inlet-outlet condition is used for the outlet patch.
An overview of the mesh is shown in Fig. 2 (b). The snappyHexMesh is used for meshing. A rotor-stator method is adopted. Specifically, the rotor is a cylinder of radius 0.25 m (reduced scale) with an axis oriented along the z-direction passing through the building center. It is rigidly rotated when studying the wind effects at different attack angles without re-meshing the geometry. The stator is the rest of the domain, as indicated in Fig. 2 (b) and (c). As can be seen, a hexahedral mesh is used for the whole domain and cell expansion gradients are used, leading to the fact that the mesh size at the domain top and outlet is around two times bigger than at the domain center. The mesh size is around 0.016 m at the inflow patch and it is around 0.004 m up to a distance of 0.1 m from the building surfaces. This yields the number of cells for the building H, L and W to be 125, 75 and 25, respectively, which fulfills the requirements in [18]. The total number of cells is around 1.8 million.
As for the numerical schemes utilized for the simulations, the pressure velocity coupling is imposed using the PISO algorithm. The time discretization is performed using the Crank Nicolson scheme with a blending coefficient of 0.85. The face fluxes for all the quantities are calculated using the bounded Gauss linear scheme except for the velocity which is calculated using the second-order LUST scheme. The k-Equation is used as the LES subgrid scale turbulence model. The non-dimensional time step ( �t * (U r /H) ) for the calculation is 0.01, yielding the Courant number of the calculation to be around 3.0 on average and only attain 6.0 in some small size cells.

Inflow
As mentioned before, synthetic turbulence is applied at the inflow. In particular, PRFG 3 is used for generating the ABL flow. It targets a Category IV profile in the AIJ standard code, which is in agreement with the one adopted in WTT. The mean velocity and the along-wind turbulence intensity, I u , distribution along height are set following the power law in the code prescription, being 0.27 the exponential and 550 m the reference height of ABL. The turbulence intensities in the y and z directions are set to be 0.75 and 0.5 times of the turbulence intensity in the x direction, respectively. We mention that in order to obtain the target wind field characteristics at the building location, the turbulence intensity of each velocity component has been increased of about 20% at the inflow to compensate for the energy dissipation between the inflow and the building location. The reference velocity, U r , is chosen as the time-averaged wind speed at the reference height, H r , which is the height of the building.
In order to check the performance of the turbulence generator, an empty domain with the synthesized turbulent inflow condition is tested first. The wind profiles at the location where the building will be placed are shown in Fig. 3. As can be seen, the mean wind speed profiles from the AIJ code, WTT and LES agree well with each other. In the same figure the turbulence intensity profiles are reported, being I u_LES , I v_LES and I w_LES at H r equal to 12%, 10% and 8%, respectively. Furthermore, the simulated turbulence integral length scales at H r , evaluated from velocity time-histories and assuming frozen turbulence, L ux , L vx and L wx , nondimensionalized to the WTT length scale, are 0.375, 0.2 and 0.125, respectively.
The power spectral density of the velocity components u, v and w at the building height are presented in Fig. 4 (a), (b) and (c), respectively, showing a good agreement with the von Karman (indicated as "Tar. Karman" in the figures) wind spectrum up to the cut-off frequency, equal to approximately 65 Hz. Overall, the simulated turbulent inflow shows a satisfactory agreement with WTT and with the targeted profiles. 3 Results and discussion

Mesh
As is well known, it is ideally necessary to ensure that results obtained using numerical simulations are independent from the adopted mesh size [28]. When using LES, it is actually more appropriate to state that results of interest shall be independent or vary within acceptable ranges with the mesh size, as a complete independency is actually not expected. We start by considering the case at 0 • . For such case, we use three meshes, namely, coarse mesh (CM), medium mesh (MM) and fine mesh (FM), to check the grid dependency of the numerical solution. The number of cells of the three meshes is 0.1, 1.8 and 6.6 million, respectively. In other words, the ratio between the cell size in each direction between two consecutive meshes is approximately 1.5, following the suggestion in [31]. A summary of the characteristics of the three meshes is provided in Table 1. The computational resources needed to perform each simulation are also reported. In order to characterize the obtained results, the time-average, rms, maximum and minimum values of C p , respectively denoted as C p , C ′ p , Ĉ p and Č p , are considered. Peak values are calculated by fitting to the two-minute extreme using a Gumbel cumulative probability distribution function and then, adopting the well-known shifting property of the Gumbel distribution [32], to extract the 10-minute extremes associated with a nonexceedance probability equal to 80%, in agreement with the Cook and Mayne approach [33]. Figure 5 presents the scatter plots of the C p statistics from simulations using the aforementioned three meshes. As shown in Fig. 5 (a) -(c), no remarkable differences in terms of C p , C ′ p and Ĉ p can be observed, except from a very slight underestimation of C p in CM (the difference is approximately 0.1). Regarding Č p , results obtained from CM only arrive at -2.5. However, they attain approximately -4.0 in the simulations with MM and FM, showing no apparent variation between the two meshes. Some localized points have more severe extreme suctions in MM and FM than in CM, while the results from CM are closer to WTT, which seem nonintuitive. However, as indicated in an analogous work [34], a finer grid model might underestimate the turbulence dissipation rate and give a higher turbulent kinetic energy prediction, yielding overall better performance with coarser meshes. In all, results highlight that CM is probably excessively coarse, as it substantially differs from the other two. Results appear to be quite similar with respect to all analyses quantities for MM and FM, so that MM is chosen for the next analyses.

Pressure coefficients statistics
After choosing the mesh size to be used, we proceed at simulating seven attack angles, equally spaced between 0 • and 90 • . The LES results are compared with WTT results for all the simulated wind attack angles in the form of three error metrics, namely, the mean error (ME), the standard deviation of error (SE) and the coefficient of determination ( R 2 ) [12], as reported below.
where the subscript i represents the pressure tap number and the total number of monitors, N, is equal to 480. The WTT (or LES ) represents the time-average and σ WTT (or σ LES ) represents the standard variation of the WTT (or LES) dataset. The error metrics  Table 2.
Notice that ME for C p is not reported, as it might be simply related to the choice of the reference pressure, so being irrelevant. Looking at Table 2, an overall correspondence in terms of C p is obtained between WTT and LES, with the average R 2 around 0.99. Similarly, good agreement could also be found regarding the Ĉ p , whose average R 2 is 0.93.
Considering all the wind attack angles, we observe that no clear trend emerges for the error metrics when the angle is varied. For instance, R 2 of C ′ p for 0 • is 0.72, lower than some skew angles such as 0.86 of 30 • and 0.76 of 75 • . However, at 90 • we obtain the highest value of R 2 among all the cases, equal to 0.89. Substantially, a good matching is obtained when the flow is orthogonal to the short side, leading to a reattached flow. In other conditions, results do not follow simple trends. The prediction of fully detached flow expected at 0 • is accurate, but it is not more accurate than other skew angles. Figures 6 and 7 show the local peak values of C p considering at the same time wind incident angles from 0 • to 90 • , i.e. the extremes from all attack angles are enveloped together, as for design purposes. In particular, we show scatter plots reporting in abscissa WTT results and in ordinates LES results. In order to allow to individuate taps belonging to different building faces, pressure taps are grouped as face I , face II , face III and face IV , following the indication in Fig. 1 (a). The results in terms of the peak Ĉ p appear to be characterized by much higher accuracy with respect to suctions. The extreme negative pressures near the building edges on face II , face III and face IV reach values of -4.0 in WTT but only attain -3.5 in LES. Conversely, LES overestimates the negative extremes near the center bottom of face IV , where the Č p arrives at -5.0 in LES but reaches only -3.5 in WTT.

Global forces and moments
The global forces for all simulated wind attack angles are reported in dimensionless form C Fx , C Fy , C Mx and C My , representing the force and moment components in alongwind (x) and cross-wind (y) directions, as shown in Figs. 8 and 9. Specifically, the time series of C p is integrated over the building surfaces to obtain time histories of the nondimensional force coefficients. Then, for each incidence angle, the peak values (including the maximum and minimum) of forces and moments are calculated following the same procedure used for the pressure peaks explained in Section 3.1. Again, for design purposes, the peaks from all the angles are enveloped and indicated as LES env. and WTT env. in the figures.  in LES are about 25% lower than in WTT. Again, discrepancies can be found for 45 • , for which the peak values of forces (and moments) have smaller values in LES than in WTT. Similar results are also found in another work [35], which declared that the wind attack angle of 15 • approaches the critical angle or glancing angle and has the minimum drag force coefficient and maximum mean lift magnitude, and the Strouhal number (St) is maximized.

Prediction errors for 45 • case
As the 45 • case appears to be the one presenting higher discrepancies between WTT and LES, we further investigate it in this section in more detail. Hereafter the Prediction Error (PE) is used as an error metric and is defined as PE S = LES S − WTT S , where the subscript S represents the C p statistics, i.e., C p , C ′ p , Ĉ p and Č p . Figure 10 presents the scatter plots of C ′ p and Č p for the 45 • case. The scatter plots of C p and Ĉ p are not reported here, since they are satisfactorily accurate. It can be seen that C ′ p on face II shows systemic underestimations from the numerical model. In fact, monitors with C ′ p approaching 0.5 in WTT are lower than 0.3 in LES. As expected, underestimations also regarding Č p can be found on these probes, with the negative extremes in LES being 60% less pronounced than in WTT. Besides, some localized differences regarding Č p can be found on face III . The values of Č p of these points are around -3.0 in WTT but are lower than -2.0 in LES. Figures 11 and 12 present the surface distributions of C ′ p and Č p as well as their prediction errors, i.e., PE C ′ p and PEČ p . More specifically, it can be clearly seen that the underestimations mentioned above are located near the leading edge of face II . Looking at Fig. 11 (a) and (b), the higher values of C ′ p are also located on this region. As for the independent spots near the top edge of face III where the highest PEČ p appears, one possible explanation can be found in a relative study [23]: the extreme value is recorded downstream of the leading edge, from which strong vortices are expected to be shed.  Figure 13 presents the frequency distributions of PE C ′ p and PEČ p . It can be seen that the PE C ′ p concentrates around -0.05. It can also be found that deviations between LES and WTT in terms of extreme suctions ( PEČ p ) are mainly inside the range -1.0 and 0.0, which indicates a tendency of the simulation to underestimate such quantity and provides a measure of the expected underestimations.

Sensitivity study
In this section, the influence of parameters including mesh setting and LES subgrid scale (SGS) model on numerical results for the 45 • case is investigated, attempting to ameliorate the previously obtained prediction errors.

Mesh
We firstly consider the possibility to ameliorate results by refining the mesh. In Section 3, we observed that the selected meshing parameters were able to guarantee a low sensitivity of the obtained results to the cell size. Nevertheless, the study was performed at 0 • , so that it cannot be automatically extended to 45 • . We here consider two mesh refinement strategies: in the first one we refine all the fluid volume around the building, while in the second one we add boundary layer cells. Then, the two strategies are combined.
Consequently, three adjusted meshes based on MM are obtained. The details of the tested meshes are listed below: (I) boundary layer cells, which have 8 layers and the minimum thickness equaling 0.0004 m and scratching ratio of 1.15, are added to the building surfaces, leading to the mesh named BL; (II) the mesh size is refined to 0.004 m up to a distance of 0.2 m from the building and to 0.002 m to a distance of 0.025 m from the building, resulting in the mesh named BR; (III) the two aforementioned strategies are applied at the same time, yielding the mesh named BL_BR. The mesh zones far away from the building are kept unchanged. An overview of the newly considered meshes in the proximity of the building is provided in Fig. 14 and their characteristics are provided in Table 3. It is worth noting that the addition of boundary   layer cells dramatically increases (double) the computational costs, though a relatively lower (half ) value of y + can be obtained at the building surfaces.
The error metrics of the three cases are presented in Table 4. Unexpectedly, the mesh with boundary layer cells, i.e., BL, has worse performance in terms of C ′ p prediction compared to the original mesh, MM. Specifically, the R 2 of C ′ p for mesh MM is 0.60 while for the mesh BL it decreases to 0.54. Differently, the BR mesh shows better performance in terms of C ′ p , reaching R 2 equal to 0.65. As it can be seen, despite such variations, the overall error does not change substantially, especially if we consider the high increase in computational resources needed to make the refinements. We also show that the combination of the two refinement strategies BL_BR does not show a more accurate prediction of C ′ p . Now we see the simulated results in terms of Č p . Again, slightly better performance can be found from BR than from BL, while the results from BR and BL_BR are close to each other. Figures 15 and 16 show the surface distributions of the prediction errors PE C ′ p and PEČ p for the three meshes. The range of PE C ′ p for BL is close to the original simulation, MM, still showing high error values on face II , which seems to indicate a detachment of the flow at the edge for WTT not predicted by LES. This can be deduced also by Fig. 11 (a) and (b) in which C ′ p is shown. The distribution of PEČ p for all the analysed meshes is shown in Fig. 16, which does not show any major difference between the three.

Subgrid scale model
It has been seen that refining the mesh did not yield strong improvements of the results, despite its high increase of computational costs. We thus here investigate the effect of changing the turbulence model. Table 5 reports error metrics of C p statistics for simulations on mesh MM with different turbulence models (namely, WALE, k − ω SSTSAS and Smagorinsky). Overall, none of them shows better model

The prediction errors of C ′
p and Č p for 45 • are not sensitive to the addition of boundary layer cells for the near-wall treatment. Some marginal improvements are obtained by refining the mesh in the surroundings of the building, but this will dramatically increase the computational time; 4. No major improvement is obtained by changing the adopted turbulence model.
Overall, it is difficult to individuate with certainty the cause of the observed discrepancies. Actually, many possible causes can be individuated. Such causes range from small deviations of the geometry with respect to WTT to differences in the incoming turbulence, despite the fact that a good matching was obtained in terms of turbulence intensity, and length scales were prescribed according to usual practice. Surely, we observe that the considered case is not particularly sensitive to the adopted turbulence model and mesh size, once an appropriate mesh has been initially selected. It must be remarked, that the addition of boundary layers proved, beside not ameliorating results, had only a very limited effect on the obtained results.