Validation of multi-fidelity urban flow numerical model for innovative air mobility against wind tunnel experiments in two facilities

Validation of multi-fidelity urban flow numerical model for innovative air mobility against wind tunnel experiments in two facilities César Sánchez-Aguadoa,Pau Varelaa,∗,Pedro Quinteroa,Roberto Navarroa,Sebastián Franchinib, Raúl Manzanares-Bercialb,Mikel Ogueta-Gutiérrezb,Israel Quintanillacand Aurea Gallego-Salgueroc aIUI CMT-Clean Mobility & Thermofluids, Universitat Politècnica de València, Camino de Vera, 46022, Valencia, Spain bInstituto Universitario "Ignacio Da Riva" (IDR/UPM). Universidad Politécnica de Madrid, Plaza Cardenal Cisneros, 3, Madrid, E-28040, Spain cDepartamento de Ingeniería Cartográfica, Geodesia y Fotogrametría, Universitat Politècnica de València, Camino de Vera, 46022, Valencia, Spain ARTICLE INFO Keywords: Computational wind engineering Urban air mobility Large eddy simulations Wind Tunnel Testing ABSTRACT The secure integration of Unmanned Aerial Systems (UAS) in urban environments highlights the need for precise characterization of wind patterns within realistic complex models. This study evaluates the differences in validation data between two wind tunnels for numerical models that support Urban Air Mobility (UAM). Two independent datasets in a LiDAR/DTM-derived complex urban scenario supply velocity and surface-pressure time-averaged measurements to validate multi-fidelity numerical results. Results show that LES improves the RANS time-averaged velocity results, reducing the normalized root mean square error from 24%to less than 10%, minimizing the maximum local error to 30% and outperforming RANS in over 70%of measurement locations. Although LES resolves unsteady wake structures and downstream peaks in turbulent kinetic energy, the employed numerical setup does not fully reproduce experimental velocity fluctuations. Despite inflow sensitivities and differences in turbulence and blockage ratios, normalized surface pressures are highly reproducible across facilities, also showing strong agreement between experiments and CFD simulations. Flow diagnostics reveal distinct zones with different risk signatures depending on local urban geometry and wake interactions. Areas with high-rise, low-density buildings produce strong shear and increase peak velocity and turbulent kinetic energy, which can create hazardous conditions for UAS operations. Results suggest that future UAM studies should validate numerical turbulence spectra and integral length scales at the inflow to ensure realistic reproduction of facility conditions. The obtained dataset enables follow-up work analyses targeting specific UAS mission requirements. 1. Introduction1 As cities continue to densify and with current concerns2 about air quality [1] and inhabitant comfort [2,3], Un-3 manned Aircraft Systems (UAS) are bound to contribute4 to the development of the decarbonized, sustainable, and5 energy-efficient Smart City. UAS are expected to become6 increasingly important in applications related to delivery,7 security, emergency response, and urban mobility [4]. In8 this regard, Urban Air Mobility (UAM) promises to revolu-9 tionize transport and services in cities in the near future; in10 fact, the European Union Aviation Safety Agency (EASA)11 forecasts UAM to become a reality within 5 years [5].12 The need for UAS secure inclusion in urban environments13 introduces a new application in the field of wind engineering,14 as urban wind patterns, characterized by complex flow in-15 teractions, are critical for drones, altering their performance16 and hindering safe operations [6], especially during take-17 off and landing at vertiports [7]. The current trend poses18 ∗Corresponding author (P. Varela) ORCID(s): 0009-0000-4079-0017 (C. Sánchez-Aguado); 0000-0002-7909-4569 (P. Varela); 0000-0003-4373-2079 (P. Quintero); 0000-0003-2587-4954 (R. Navarro); 0000-0002-1476-9174 (S. Franchini); 0000-0001-8987-1694 (R. Manzanares-Bercial); 0000-0001-9183-450X (M. Ogueta-Gutiérrez); 0000-0003-0059-5975 (I. Quintanilla); 0000-0002-0793-8902 (A. Gallego-Salguero) an emerging challenge: understanding the velocity and tur19 bulence of actual urban fields to aid UAS in path planning 20 and optimization. Governments are starting to implement 21 regulations that affect the operability of aerial vehicles [8], 22 and are delineating no-fly zones based on environmental 23 constraints, such as acoustic emissions, obstacle proximity, 24 and strong wind gusts [9]. In turn, these restrictions are gov25 erned by a building-induced wind speed field and turbulent 26 instabilities, which impose additional power requirements 27 and compromise drone maneuverability [10]. 28 Wind tunnel experiments stand out as the most reliable 29 data acquisition method in the field of wind engineering. 30 However, they are constrained to a limited number of dis31 crete measurement points, which cannot resolve the full 32 spatial variability of urban flows required for precise UAS 33 trajectory optimization [11]. On the other hand, numerical 34 approaches offer a faster and more versatile alternative for 35 UAM applications. Giersch et al. [12] emphasized that pre36 computing flow solutions across representative boundary37 layer conditions could enable a rapid evaluation of the flight 38 trajectory based on actual meteorological conditions, for 39 which a location-specific analysis must be developed. It 40 can be inferred that failing to consider all relevant wind 41 directions can lead to inadequate designs and potential risks. 42 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 1 of 15 To study a representative amount of different wind condi-43 tions, and because each wind direction requires a full inde-44 pendent simulation, computationally affordable but under-45 resolving building-scale turbulence simulations should be46 leveraged against higher-fidelity methods. Many authors47 have investigated the strengths and limitations of differ-48 ent turbulence approaches in urban settings in the field49 of Computational Wind Engineering (CWE), emphasizing50 differences in results using LES and RANS. García-Sánchez51 et al. [13] quantified LES uncertainties and concluded its52 superior capabilities regarding high-turbulent flow predic-53 tions. However, Blocken [14] justified the use of RANS54 simulations for routine calculations given the availability of55 best practice guidelines [15] and their lower computational56 requirements. In terms of time-averaged velocity fields, ac-57 cording to Blocken et al. [16], steady RANS simulations can58 provide accurate results with an error margin of around 10%59 in areas where the amplification factor, that is, the increase in60 local wind speed with respect to the wind speed at the same61 location without buildings, is greater than 1. The RANS62 increase in accuracy with the amplification factor is par-63 ticularly important in high-rise building complexes, where64 wind speeds can amplify up to 150%of the free-stream65 velocity [17] due to strong channeling effects that accelerate66 the wind through narrow urban canyons. In contrast, when67 the amplification factor is less than 1, the accuracy of RANS68 simulations deteriorates significantly due to shortcomings69 of the RANS method in how the averaged scalar velocities70 are defined [18]. Results also indicate an underestimation71 of the wind speed at the windward corners of the building72 and the overestimation of the recirculating flow in the wake73 of the buildings in comparison to the results of wind tunnel74 experiments [19].75 As aircraft are bound to fly through the turbulent wake76 of buildings [9], urban flow studies applied to UAM require77 accurate characterization of the turbulent field within the78 microscale urban fabric [20]. According to Adkins [21],79 smaller eddies comparable in size to the UAS are of the80 most significant concern for these aircraft. Resolving these81 high-frequency wake effects requires advanced modeling82 techniques [22] due to the inherent unsteadiness of turbulent83 phenomena. It follows that UAM studies tend to employ high84 spatial and temporal resolution, and as Roseman and Argrow85 [23] note, LES is particularly valuable for capturing the flow86 structures that must be quantified for safe UAS operations87 and risk assessment. García-Sánchez et al. [13] also report88 that LES is more accurate than RANS in capturing turbulent89 kinetic energy (𝑘) in 80%of measurement locations. Despite90 proving to have more precise results, studies on inflow un-91 certainties [18,24] highlight additional difficulties for LES.92 García-Sánchez states that in regions where inflow uncer-93 tainties have a greater effect on the flow, the out-performance94 of LES over RANS in terms of averaged velocity prediction95 drops to 50%. In specific contexts, its high sensitivity to96 inlet conditions and the increased computational cost re-97 quired to resolve subgrid-scale turbulence justifies the use98 of faster and less expensive RANS techniques [16]. With99 hybrid approaches yielding very few successes [18] due to 100 the difficulty of switching between RANS and LES, these 101 conditions underscore the need for a multi-fidelity analysis 102 to assess the effect of turbulence treatment in urban flows. 103 A major limitation across current UAM research frame104 works is the scarce collection of experimentally validated 105 urban flow numerical simulations [11,12]. Risk assess106 ments for UAM operations frequently depend on velocity 107 and turbulence fields from RANS simulations; however, the 108 limited validation against wind-tunnel experiments leaves 109 residual uncertainty in safety and performance assessments 110 for UAS operations [25]. Furthermore, comparative wind111 tunnel studies have shown that facility-specific factors, such 112 as blockage ratios, incorrect reference static pressure in 113 open-circuit tunnels, or differences in the turbulence inte114 gral scale, can produce systematic shifts in mean and peak 115 results [26]. While these discrepancies have been studied 116 in single high-rise buildings, their impact on actual urban 117 environments is yet to be analyzed. Modern drone flight 118 planning optimization tools that rely on high-accuracy CFD 119 wind fields [27] have also been demonstrated on simplified 120 building geometries. Upcoming trajectory planning efforts 121 must expand their scope to more practical and realistic lay122 outs. Thus, the current lack of experimental and high-fidelity 123 numerical urban flow datasets represents a significant barrier 124 to both the accurate validation of UAM risk studies and the 125 future application of flight-planning systems to actual urban 126 scenarios. 127 The primary objective of this study is to assess the im128 pact of facility-specific differences in wind-tunnel environ129 ments on the aerodynamic response of a complex, real-case 130 urban model and to compare those results with CFD predic131 tions. To this end, two different wind tunnels have studied 132 the same 1:375-scale mock-up of the authorized UAS flight 133 zone Urban Lab in Benidorm, Spain. These comparisons 134 ensure the calibration of multi-fidelity CFD simulations, 135 which are more suitable for UAM applications, by assessing 136 whether CFD results lie within empirical bounds. As dis137 cussed above, a specific, validated urban wind database is 138 needed to apply CFD predictions to wind-aware real-time 139 UAS mission planning [28]. To ensure reproducibility of 140 results, this work compiles all the steps necessary to perform 141 the CFD urban flow study and its experimental validation. 142 To the authors’ knowledge, this is one of the first attempts in 143 the literature to perform a cross-validation of experimental 144 results on an actual LiDAR/DTM-derived urban mock-up 145 in two wind tunnels. Furthermore, the present study also 146 contributes to the need to expand the existing urban flow 147 database [29] with a calibrated digital twin of an operational 148 flight zone suitable for on-site drone actuation testing and 149 mission planning. 150 The remainder of this paper is organized as follows. Sec151 tion 2describes the two experimental campaigns, including 152 the wind tunnel facilities, boundary layer generation tech153 niques, and the field geometry employed, as well as an in154 depth description of the data acquisition process. The inflow 155 velocity, turbulence characteristics, and their corresponding 156 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 2 of 15 results are compared. The numerical methodology is seen in157 section 3. This section justifies the computational domain,158 mesh generation strategy, boundary conditions (BC), and159 turbulence modeling approaches employed in the study.160 Section 4validates the performance of different turbulence161 models against experimental data, with a focus on time-162 averaged flow characteristics. Finally, the assessment of163 wind tunnel effects on experimental results is presented.164 2. Experimental setup165 As previously stated, to investigate the effect of wind166 tunnel characteristics and assess the reproducibility of the167 study, experimental tests were performed in the facilities168 of two universities. The first set of experiments was con-169 ducted in the ACLA-16 wind tunnel by researchers from170 the Instituto Universitario “Ignacio Da Riva” (IDR/UPM).171 These studies were validated against tests carried out by the172 CMT-Clean Mobility &Thermofluids Institute (UPV) in the173 Profesor Francisco Payri wind tunnel (from now on, the174 authors refer to two sets of experiments and their results as175 IDR or CMT, respectively). This comparison aims not only176 to validate the employed urban flow experimental procedure177 but also to highlight key features that might alter the flow178 behavior in the region of interest and assess the repeatability179 of results. By using the same model in both experiments,180 differences in results can be attributed to differences in the181 Atmospheric Boundary Layer (ABL) generators, the wind182 tunnel test section, and atmospheric air conditions.183 2.1. Wind tunnels184 The IDR ACLA-16 [30], located in the Montegancedo185 campus premises of the Universidad Politécnica de Madrid,186 is a closed test section, closed circuit wind tunnel that187 features a square cross-section, measuring 2.2 m per side,188 and an overall length of 20 m, subdivided into a flow-189 conditioning chamber and a measurement section. The wind190 tunnel is equipped with 16 fans, each with a total power of191 7.5 kW, which allows a maximum test velocity of 32 m/s192 with a turbulence intensity of approximately 1.8%.193 On the other hand, the CMT Profesor Francisco Payri194 [31] open wind tunnel features a working section length of195 22 m and a square cross-section with a side length of 2.84196 m, trimmed with chamfers at all corners. With nine 45 kW197 fans, the wind tunnel has a maximum velocity of around 32.5198 m/s and a turbulence intensity of 0.2%. This facility is also199 divided into different sections to allow the measurement of200 two independent studies.201 2.2. Atmospheric Boundary Layer configuration202 Flow conditioning was conducted independently in the203 two facilities to reproduce the scaled inflow properties im-204 pinging on the model with the real ABL profile developed205 over Benidorm. This profile is influenced by the effect of206 obstacles surrounding the measurement region. Depending207 on the average roughness of these obstacles, the European208 wind standard specifies five different types of terrain [32],209 which describe a logarithmic velocity and turbulence inten210 sity profiles. The Eurocode 1 [32] establishes a neutral strat211 ification approach to tackle urban airflow analyses, which 212 defines the dominant behavior of neutral ABLs by assuming 213 a null mean vertical wind speed [12]. In that sense, studies 214 on the effect of thermal effects, which usually create vertical 215 instabilities in the wind profile, revealed that thermally216 generated turbulence is secondary compared to mechanical 217 effects [21]. Furthermore, an approach that neglects thermal 218 buoyancy effects has been successfully applied in similar 219 studies [13,33,34]. Due to the terrain’s orography and the 220 Urban Lab’s characteristics, a Eurocode 1 terrain Category 221 III was targeted by both wind tunnels. The velocity profile 222 has been approximated by a power law as per the ASCE 7 223 [35] Category C, 𝑢=𝐶⋅𝑧𝛼, where 𝐶is an arbitrary constant 224 and 𝛼= 0.2.225 ABL generator setups can typically be divided into 226 spires, which primarily increase the turbulence intensity 227 of the flow, and subsequent surface roughness elements, 228 which reduce the velocity of the flow near the ground to 229 create the ABL. Due to the different placement of the model 230 within the wind tunnel and especially the disparity in cross231 sections and turbulence characteristics of the flow, a different 232 turbulence generator layout is designed to transform the 233 uniform-velocity incident flow into the desired ABL for each 234 wind tunnel. The two configurations used in this paper can 235 be seen in Figure 1. The available test section of the IDR set 236 of experiments allows for a more gradual transition of the 237 velocity profile: the IDR experiment incorporates different 238 sizes of cubic roughness elements, ranging from 150 to 70 239 mm per side. The assembly also includes a 300 mm height 240 barrier and 220×1930 mm triangular spires in base width 241 (𝑤) and height (ℎ), respectively. Conversely, CMT’s set of 242 experiments uses 65×90×90 mm rectangular blocks and 243 180×1420 mm 𝑤×ℎspires. Main features are gathered in 244 Table 1.245 A comparison of normalized vertical wind profiles for 246 the two tests against the Eurocode Type III ABL profile 247 is presented in Figure 2, along with the normalized root 248 mean square error (NRMSE) values calculated in the region 249 of interest, that is, the airspace that typically limits UAV 250 flights, from 20 to 120 m in height [10]. Due to the dif251 ference in absolute velocity between the two experiments, 252 non-dimensional vertical profiles are expressed with respect 253 to an arbitrary 𝑧𝑟𝑒𝑓 = 60 m. The NRMSE values of the 254 streamwise velocity are calculated according to the follow255 ing equation: 256 NRMSE = RMSE 𝜙ref =√1 N∑N i=1(𝜙i−𝜙ref i)2 1 N∑N i=1 𝜙ref i ,(1) where 𝜙𝑟𝑒𝑓 represents the reference Eurocode profile. NRMSE 257 results indicate a good agreement with the Eurocode ref258 erence, with the boundary-layer thickness (∼0.40 m at 259 model scale) matching the reference Eurocode 1 Type III 260 within ±5%error. Experiments show a greater difficulty in 261 matching the intensity profile of the Eurocode. A maximum 262 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 3 of 15 Figure 1: Effective test sections for IDR (top) and CMT (bottom) wind tunnels and ABL turbulence generator setups NRMSE of 20%can be seen in Figure 2. This error, however,263 falls within reasonable margins [36], especially considering264 the building-induced dominant turbulent flow regime that265 occurs after the first line of buildings.266 u=uref [!] 0.5 1 1.5 z=zref [!] 0 1 2 3 4 Eurocode Type III ASCE Cat. C (0.78%) IDR (4.09%) CMT (2.70%) Iu[!] 0 0.2 0.4 z=zref [!] 0 1 2 3 4 Eurocode Type III IDR (20.12%) CMT (21.74%) Figure 2: ABL velocity (left) and turbulence intensity (right) comparison. 𝑧𝑟𝑒𝑓 = 60 m. The NRMSE with respect to the Eurocode Type III profiles is seen in parentheses 2.3. Field geometry description and mock-up 267 The study focuses on a designated UAS flight zone of 268 the city of Benidorm, shown in Figure 3. The geometry 269 used in this study is a simplification of the landscape and 270 urban fabric of a 350-meter-radius circular domain around 271 the UAS flight zone. Note that significant discrepancies in 272 measurement points can usually be attributed to geometry 273 simplifications [19,37]. In this case, obstacles below a 274 height of 6 m were removed, while maintaining only relevant 275 buildings that influence the overall flow behavior. To recon276 struct the geometry of the buildings and the terrain, high277 resolution geospatial data were obtained from open-access 278 cartographic databases. The building perimeters were down279 loaded as vector layers (Shapefile format) from the Elec280 tronic Office of the Spanish Cadaster [38]. The raster files 281 corresponding to the Digital Terrain Model (DTM) [39] and 282 the Normalized Digital Surface Model of Buildings (nDSM) 283 [40], both with a 1-meter mesh spacing, were obtained from 284 the Valencian Spatial Data Infrastructure (IDEV). The DTM 285 contained the absolute terrain elevations, while the nDSM 286 provided the height of the buildings relative to the ground. 287 Using ArcGIS Pro software [41], both digital models were 288 combined through map algebra, thereby generating a model 289 of the rooftop elevations. Maximum building height values 290 were extracted using the cadastral footprint as a mask. This 291 information enabled the interpolation and extrusion of the 292 3D building volumes, resulting in an accurate representation 293 of the urban morphology. Additionally, to further increase 294 the geometric fidelity of the model, a simplified represen295 tation of the terrain relief was generated using 6-meter 296 elevation contour lines extracted from the DTM, which were 297 simplified and transformed into vector polygons to represent 298 the stepped elevation layers, later used in the topographic 299 model. The final geometry of the city is seen in Figure 4.300 To ensure compatibility with the available wind tunnel 301 facilities, the experimental model was built at a scale of 302 1:375. This scale was selected to maximize the model size 303 within the constraints of the test section while preserving 304 the dominant flow features of interest. Therefore, the block305 age ratio inside the wind tunnel test chambers is less than 306 5%. Moreover, for bluff-body dominated urban flows, the 307 Reynolds number (Re) independence of the wake character308 istics allows for reliable extrapolation of results from scaled 309 models [42], so long as the Re satisfies a critical threshold 310 dependent on the roughness Reynolds number, Re𝑧0, and the 311 dimensionless height, 𝑧+[43]. This assumption is commonly 312 accepted in urban wind engineering and supports the use of 313 reduced-scale models for aerodynamic studies [44]. Taking 314 into account the velocity profile and the main geometric 315 properties of the buildings, the Re of the model, defined as 316 𝑅𝑒 =𝑈 𝑙 𝜈,(2) where 𝑈is the wind speed, 𝑙is the building width, and 𝜈is 317 the air kinematic viscosity, lies between 2⋅104and 6⋅104in 318 the area of interest. Hence, this experimental setup ensures a 319 wind speed-independent flow regime, as this range is above 320 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 4 of 15 Figure 3: Aerial view of Benidorm city. The North direction corresponds to the upward direction. The authorized flight area for UAS is highlighted in red the critical 𝑅𝑒𝑐𝑟𝑖𝑡 = 1.1⋅104recommended for building wind321 tunnel testing [44].322 A flush-mounted rotating platform accommodates the323 scale model and allows for accurate orientation at different324 angles of incidence. Based on the wind rose of Benidorm325 city, this work analyzes the two principal wind directions:326 N-S (or 𝛽= 0◦), and S-N (or 𝛽= 180◦). The asymmetric327 orography of the terrain forces an elevation between the328 wind tunnel’s floor and the model’s ground. In turn, this step329 introduces a detachment of the flow in the model, which330 would alter the ABL velocity and turbulence profiles from331 their actual behavior. To mitigate the impact of the step332 and ensure a smooth transition between the ABL generation333 zone and the urban model, two transition ramps are included334 on each side of the city in the experimental and numerical335 analyses. A parametric numerical study of the ramp geom-336 etry revealed that steep slope angles introduce a contraction337 in the cross-section, accelerating the flow and inducing a338 contraction of the ABL thickness. Nevertheless, the ramps of339 both studies—between 1◦and 3.5◦, depending on the facility340 and the angle of incidence of the flow—do not significantly341 alter the dimensionless behavior of the ABL [45].342 2.4. Instrumentation and data acquisition343 In IDR’s set of experiments, surface pressure has been344 recorded via 294 taps (24 buildings instrumented + 15 Irwin345 probes) connected to 5 Scanivalve ZOC 33/64PxX2 sen-346 sors, ±10 in H2O (2488 Pa) full scale range, and ±0.15 % full347 scale accuracy. The sampling has been performed at 150 Hz348 sampling rate for a 120 s time period. A Pitot tube mounted349 outside of the boundary layer and any wake-interference350 position has provided the reference dynamic pressure for351 𝑐𝑝computation. Hot-wire anemometry (Dantec 55P91) has352 been used to measure the three velocity components at 21353 heights per probe location up to 0.42 m for 𝛽= 0◦and354 𝛽= 180◦, with sampling at 5 kHz for 104 s. Furthermore,355 Figure 4: Mock-up of a selected region of Benidorm city scaled 1:375. Buildings in gray are not instrumented. Kronos building is painted black as a reference the Irwin probes [46], each consisting of two pressure taps 356 separated by 1 cm in the vertical direction, have been used to 357 calculate wind velocity close to the ground. These data yield 358 both mean and turbulent fluctuations for direct validation of 359 CFD velocity fields and turbulence statistics. 360 CMT’s set of experiments consisted of a total of 124 361 pressure taps, measured using a Chell microDAQ3 pressure 362 scanner. This instrument offers a maximum absolute range 363 of 0.5 kPa to 400 kPa, with an error of 0.02%. The sampling 364 frequency was 400 Hz for a 120-second time period. Al365 though a Pitot tube was used to characterize the ABL, these 366 tests only included pressure measurements in the mock-up. 367 Differences in the acquisition process can be seen in Table 1.368 To help understand the measurements across the study, 369 the top view of the model is displayed in Figure 5. The 370 buildings instrumented with pressure probes are depicted in 371 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 5 of 15 Table 1 Comparison of experimental setup characteristics in both wind tunnels Category Parameter IDR CMT Wind tunnel facilities Cross-section [m] 2.2 ×2.2 2.8 ×2.8, chamfered Length [m] 20 22 ABL generation setup Spires (𝑤×ℎ) [mm] 220 ×1930 180 ×1420 Roughness elements dimensions [mm] 150 to 70 per side 65 ×90 ×90 Barrier height [mm] 300 – Upstream development region [m] 11.1 8.6 Instrumentation Velocity measurement method Hot-wire anemometry – Velocity probe sensor Dantec 55P91 – Velocity sampling frequency [Hz] 5000 – Pressure scanner Scanivalve ZOC 33/64PxX2 Chell microDAQ3 Number of pressure taps 294 124 Pressure sampling frequency [Hz] 150 400 Experimental notes 𝑝𝑎𝑚𝑏 [Pa] 92695 101997 Measurement time [s] 120 120 blue, surrounding the UAS airspace outlined in red. Colored372 dots locate the position of velocity probes: red and black373 dots represent hot-wire vertical measurement line locations374 for 𝛽= 0◦and 𝛽= 180◦, respectively, while green dots375 highlight the ground measurement position for Irwin probes.376 Figure 5: Top view of the city model. In blue, instrumented buildings. Black and red dots represent hot-wire measurement positions for 𝛽= 0◦and 𝛽= 180◦, respectively. Green dots indicate the location of Irwin probes 3. Numerical setup377 The numerical calculations were performed on UPV’s378 “Sirius” cluster using the commercial software Star-CCM+379 (v19.04.007). Both LES and RANS simulations have as-380 sumed incompressible flow, given the low Mach number of381 the flow. The following sections describe the computational382 domain and mesh, the BC, and the numerical settings for383 both stationary and unsteady simulations. Then, the reso384 lution of the LES setup assesses the reach of high-fidelity 385 simulations to capture turbulent phenomena. 386 3.1. Computational domain and mesh 387 The computational domain has been selected to repli388 cate IDR’s experimental procedure and is seen in Figure 6.389 Numerical simulations replicate IDR’s experiment because 390 they are more suitable for UAM studies, given the availabil391 ity of both pressure and velocity probes, which provide more 392 information regarding the flow mean and fluctuating velocity 393 that may affect UAS operations. In this type of study, it is 394 common to create a domain that is sufficiently large to avoid 395 disturbing the flow around the buildings [17]. This time, 396 however, the selected domain replicates the dimensions of 397 the wind tunnel, a procedure successfully employed by the 398 AIJ [19], as this would provide a more accurate assessment 399 of the influence of the walls in the model flow measurements. 400 As described in subsection 2.1, the IDR’s wind tunnel do401 main presents a square section with a width of 2.2m, with the 402 city centered on the floor. The length of the wind tunnel has 403 been shortened upstream of the city to match the inlet of the 404 simulations with the ABL experimental measurement line, 405 thereby avoiding the need for an ABL conditioning region 406 and saving computational costs. In total, the domain extends 407 4.7 m upstream of the model and 15 m downstream of it. 408 The meshing strategy employs Star-CCM+’s trimmed 409 mesher algorithm [47] to generate a predominantly hexahe410 dral grid. As shown in Figure 6, the mesh is progressively 411 refined within a dedicated refinement region enclosing the 412 urban model, ensuring increased spatial resolution in areas 413 of complex flow development. Other refinement areas, as 414 characterized in Table 2, include the ground, the other wind 415 tunnel walls, and the city, with its building wakes. This 416 method yields better results for complex geometries [48] and 417 ensures a smooth transition between coarse and fine cells, 418 which facilitates the convergence of the simulations. 419 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 6 of 15 Ground City WT walls Building wake Refinement region N-S flow direction Figure 6: Domain close-up and differentiated mesh refinement areas: ground, city, and wind tunnel (WT) walls. The mid plane displays the mesh refinement region and cell distribution around buildings Table 2 Mesh parameter description. 𝐻𝑟𝑒𝑓 = 0.7m for RANS and 𝐻𝑟𝑒𝑓 = 0.45 for LES Parameter Domain Wind tunnel walls Ground Refinement region City Cell size (%of 𝐻𝑟𝑒𝑓 ) 20 6.25 1.56 2.21 0.78 Number of PL – 5 7 7 7 PL thickness (%of 𝐻𝑟𝑒𝑓 ) (RANS/LES) – 5.7/2.22 3.6/1.11 1.4/1.11 1.4/0.44 Wake size [m] – – – – 0.3 Research studies such as Murakami’s [49] have stressed420 the difficulty of flow modeling around bluff-body sharp421 edges, characterized by the anisotropic distribution of the422 strain-rate tensor. Based on guidelines for RANS simulations423 regarding the minimum number of cells per building face424 [50,51,52], the cell size was chosen to target at least ten425 cells on each side of a building to accurately capture flow426 separation around the upwind corners. RANS simulations427 also cluster an average of 5 cells within the first 1.5 m428 above the ground, which satisfies AIJ’s guidelines that rec-429 ommend a minimum of 3-4 cells in this region [15]. For430 LES simulations, the described grid resolution is increased431 to resolve smaller turbulence scales. Both grids are described432 in Table 2, where 𝐻𝑟𝑒𝑓 = 0.7m for RANS and 𝐻𝑟𝑒𝑓 = 0.45433 for LES.434 The employed boundary layer refinements (see Table 2)435 yield an average 𝑦+= 0.95 in the region of interest, with436 96.6%of elements exhibiting a 𝑦+<2. The resulting grid437 used for RANS simulations, chosen after conducting a mesh438 independence study, has approximately 12.8 million cells439 and is shown in Figure 7. On the contrary, the LES grid440 was obtained by reducing the RANS grid’s base cell size,441 resulting in a mesh of 59.6 million cells. The quality metrics442 employed to ensure that the grid resolution is sufficient to443 resolve the relevant flow structures are described in subsec-444 tion 3.4.445 Figure 7: Mesh zoom-in of the building region. Surface elements and fluid cells are painted in black and red, respectively 3.2. Boundary Conditions 446 García-Sánchez et al. [13] state the challenge of defining 447 precise BC in LES simulations. The primary difficulty lies in 448 defining fluctuating inlet conditions that reproduce both the 449 time-averaged and turbulent characteristics of the velocity 450 field. In fact, uncertainty related to the inflow conditions 451 can decrease the accuracy of high-fidelity simulations [13]. 452 Regarding inflow BCs, Tominaga [18] establishes different 453 approaches to replicate wind tunnel experiments. In contrast 454 to precursor and recycling approaches, which require mod455 eling the entire wind tunnel domain to predict the ABL de456 velopment, the Synthetic Eddy Method (SEM) [53] applied 457 in this study introduces turbulent structures directly at the 458 inlet. To define the inflow ABL profile, experimental data 459 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 7 of 15 are required for the mean velocity distribution, Reynolds460 stresses, and integral length scale. Synthetic eddies then461 develop over the terrain, replicating the mean flow and fluc-462 tuation characteristics of the experimental measurements.463 This approach eliminates the need for precursor domains and464 associated turbulence-generation regions, thereby reducing465 computational costs.466 In the present LES setup, turbulent statistics have been467 derived from wind tunnel measurements of streamwise (𝑢)468 and vertical (𝑤) velocity components. The transversal com-469 ponent has been statistically approximated by assuming its470 variance is equal to that of the vertical component. The471 choice is physically consistent with the lower order of mag-472 nitude of vertical fluctuations relative to the streamwise473 component, providing a conservative estimate for the 𝑘.474 Moreover, García-Sánchez et al. [13] emphasize that uncer-475 tainties in the inflow turbulence decay downstream, where476 local turbulence is primarily governed by canopy-scale shear477 and building interactions. The integral time scale, 𝜏𝑢(𝑧), has478 been obtained by integrating the normalized autocorrelation479 of 𝑢′up to its first zero crossing. Then, the length scale (𝐿𝑢)480 was computed using Taylor’s frozen turbulence hypothesis,481 following Equation 3.2.482 𝐿𝑢=𝑢⋅𝜏𝑢(3) For the RANS simulations, turbulence quantities have483 been defined in accordance with standard practice and the484 AIJ guidelines [15]. The 𝑘was determined from the mea-485 sured RMS of velocity fluctuations, 𝜎, using equation Equa-486 tion 3.2. The dissipation rate 𝜖was estimated assuming local487 equilibrium of the production term 𝑃𝑘in the 𝑘equation488 and using the Boussinesq approximation, leading to the489 expression seen in Equation 3.2, where 𝑈is the vertical490 velocity profile, and 𝐶𝜇is the model constant (= 0.09)[15].491 𝑘(𝑧) = 𝜎2 𝑢+𝜎2 𝑣+𝜎2 𝑤 2≃𝜎2 𝑢(4) 𝜖(𝑧) ≃ 𝐶1∕2 𝜇⋅𝑘⋅d𝑈(𝑧) d𝑧(5) All solid boundaries were treated as no-slip walls. This492 replicates the physical test-section confinement: the wind493 tunnel walls and building faces were modeled as smooth494 walls with low-Reynolds treatment near the wall. The outlet495 is configured as a zero-gradient, constant-pressure boundary.496 3.3. Solver settings497 The steady model has been employed when using Reynolds-498 Averaged Navier–Stokes (RANS) equations with the realiz-499 able k-𝜖turbulence model, which is seen to provide more500 accurate results in the AIJ’s CWE benchmark [19]. Near-501 wall treatment resolves the viscous sublayer, with the first502 layer 𝑦+≈ 1 in the city region, thereby capturing the near-503 wall velocity gradient and turbulence production explicitly504 rather than relying on wall functions. Resolving the near-505 wall velocity profile also improves the prediction of wall506 pressure coefficients (𝑐𝑝), which are used to compare the two507 experimental studies. While this work, due to the 1:375 scale 508 of the model, resolves the viscous sublayer with a 𝑦+around 509 1 in the region of interest, it is worth noting that this would 510 not be possible with a real scale factor, and wall functions 511 would be required [13]. RANS simulations have employed a 512 segregated solver with SIMPLE pressure-velocity coupling 513 and second-order convection schemes for both momentum 514 and turbulence equations. 515 Unsteady, scale-resolving simulations have used LES 516 with the Wall-Adapting Local Eddy-viscosity (WALE) subgrid-517 scale model [54]. The WALE model assumes a mixing518 length type equation for the subgrid scale viscosity and 519 is less dependent on the value of the model coefficient, 520 𝐶𝑤[55]. The WALE model has been chosen because it 521 provides correct near-wall scaling for the eddy viscosity 522 without requiring a dynamic procedure, making it robust 523 for complex, wall-bounded urban geometries. A bounded524 central differencing scheme has been used for the spatial 525 discretization, while second-order implicit schemes have 526 been applied for temporal discretization. LES simulations 527 have been run with a global time step of Δ𝑡= 2 ⋅10−4s, 528 matching the experimental velocity sampling frequency. 529 The simulations have been run for a total time of 3s ≃530 55𝜏𝑐𝑜𝑛𝑣, where 𝜏𝑐𝑜𝑛𝑣 represents the convective time of LES 531 simulations, based on the maximum height of the buildings 532 (0.4 m) and the free-flow velocity of approximately 7.4 m/s. 533 This interval provides a meaningful sampling period for 534 turbulent statistics, consistent with the amount of time scales 535 utilized in previous LES urban studies [13]. 536 All simulations were performed on UPV’s “Sirius” high537 performance computer cluster, equipped with Intel Xeon 538 8480 56-core processors. Each RANS case required between 539 6 and 8 hours of wall time to converge, corresponding to 540 approximately 100 CPU hours per run. In contrast, each 541 LES simulation required approximately 67 wall-clock hours, 542 nearly 45 000 CPU hours per case. LES simulations em543 ployed six compute nodes, each with 112 cores. 544 3.4. LES mesh quality assessment 545 The adequacy of the LES spatial resolution has been 546 verified by the 𝑘ratio, which assesses the accurate prediction 547 of the energy content of turbulent structures. The 𝑘ratio 548 𝜂, described in Equation 3.4, is applicable for almost any 549 complex flow and has been successfully used in bluff-body 550 studies [55]. It relates the ratio of the 𝑘of the calculated 551 non-filtered eddies, 𝑘𝑐, and the kinetic energy introduced 552 by the subgrid-scale model, 𝑘𝑆𝐺𝑆 . The higher the ratio, 553 the closer the solution is to a Direct Numerical Simulation 554 (DNS) resolution (𝜂= 1); according to Torregrosa et al. 555 [55], a correct resolution of a LES can be assumed from 556 𝜂 > 0.7. Values of 𝜂below this threshold in areas with 557 high turbulence indicate increased sensitivity to the SGS 558 closure and grid resolution, leading to errors in the wake and 559 recirculation regions. 560 𝜂=𝑘𝑐 𝑘𝑇 𝑂𝑇 𝐴𝐿 =𝑘𝑐 𝑘𝑆𝐺𝑆 +𝑘𝑐 (6) Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 8 of 15 The distribution of 𝜂(expressed as a percentage) across561 the urban domain is mapped in Figure 8. More than 99.5% of562 the cells in the city region exhibit a 𝑘ratio higher than 80%,563 indicating that the bulk of the domain is well resolved with564 the used grid. A scarce group of cells present ratios below565 70%, but they only cluster near building corners and sep-566 aration regions, especially at high-rise and windward first-567 line buildings. In the first-row buildings, the high-velocity568 inflow (high local Re) ABL meets sharp building edges, thus569 combining strong local production 𝑃𝑘with thin shear layers.570 Linear shear-layer theory predicts that the most unstable571 Kelvin–Helmholtz modes produce eddies with wavelength572 on the order of the shear-layer thickness, 𝛿, thereby shifting573 the dominant energy-containing scales and the dissipation574 range toward smaller scales. This explains the lower values575 of 𝜂in these areas for a given LES filter width, which is576 proportional to cell size. High-rise buildings suffer similar577 effects; as Roth [56] notes, the flow above the urban canopy578 behaves like a mixing layer between the slower within-579 canopy flow and the faster above-canopy wind, which ex-580 plains why roof-level high shear is produced. Downstream of581 the first row, the interaction between shear-generated eddies582 mixes the flow, increasing the integral scale. As the grid583 filter is maintained, a larger fraction of the energy-containing584 motions is resolved rather than modeled.585 Given that the low-resolution cells constitute a very586 small fraction of the domain and are confined to geometri-587 cally singular regions, the present grid was deemed accept-588 able for the objectives of this study. Since the small-scale589 energy they generate tends to dissipate downstream, these590 cells make a negligible contribution to field mean metrics.591 Furthermore, a refinement based on 𝜂in those localized592 regions would reduce the SGS share but at a substantially593 higher computational cost.594 Flow direction η Figure 8: Calculated (𝑘𝑐) over total (𝑘𝑐+𝑘𝑆𝐺𝑆 )𝑘ratio (𝜂) at plane 𝑧= 0.1m (37.5 m full-scale height). Cells with 𝑘ratio (𝜂 < 0.7) are highligthed in black 4. Results and discussions 595 In this section, experimental data are compared with 596 RANS and LES results, highlighting errors in the nu597 merical approach and the shortcomings of different tur598 bulence models. The analysis of results focuses on the 599 time-averaged velocity (subsubsection 4.1.1) and 𝑘field 600 (subsubsection 4.1.2), as these variables are most relevant 601 for UAS flight performance and stability. Since velocity 602 measurements are only available from IDR’s experiments, 603 these data are compared with the corresponding numerical 604 predictions. As explained in subsection 2.4, these compar605 isons include velocity measurements in vertical probes for 606 𝛽= 0◦and 𝛽= 180◦, although comparisons between 607 CFD and wind tunnel results are only displayed for 𝛽=608 0◦to avoid repetition. Furthermore, the time-averaged 𝑐𝑝609 shown in subsubsection 4.1.3 provides a complementary 610 measure of numerical accuracy. Finally, the reproducibil611 ity of experimental results is discussed in subsection 4.2 612 by comparing the 𝑐𝑝distribution measured in Benidorm’s 613 Urban Lab across the two wind tunnels. 614 4.1. Experimental validation of RANS and LES 615 time-averaged flow fields 616 4.1.1. Comparison of mean velocity field 617 The normalized mean velocity field with respect to the 618 inflow velocity at a representative 37.5 m full-scaled height 619 is shown in Figure 9. The relation between local velocity 620 and the unperturbed inflow velocity is defined as the am621 plification factor [16], and directly affects numerical accu622 racy. Both numerical approaches reproduce the general flow 623 trends; however, RANS underpredicts velocity within the 624 wake region of all buildings, where smaller amplification 625 factor values are obtained. As Leitl et al. [57] suggested, the 626 inherent steadiness of the RANS provokes fewer turbulence 627 structures close to the walls. In turn, less turbulence mixing 628 is produced, resulting in an overestimation of the detached 629 flow in the wake of the buildings. This not only affects the 630 velocity of the wake region but also induces error in the 631 reattachment point of the flow. While RANS usually matches 632 LES results for amplification factors greater than 1, the ab633 sence of vortex shedding also explains the overestimation of 634 mean velocities in flow acceleration corridors, where actual 635 flow intermittently accelerates due to shear-layer instabilities 636 rather than sustained channeling. 637 The flow downstream of buildings exhibits different be638 havior depending on local geometry and wake interactions. 639 Atzori et al. [58] studied flow past simple tandem-obstacle 640 configurations and identified three distinct flow regimes for 641 different building separations. As the separation distance 642 increases, the flow regime is characterized by skimming, 643 where the shear layer passes directly over both obstacles with 644 limited penetration between them, by wake interference of 645 the upstream obstacle with the front face of the downstream 646 one, and finally by an independent wake with minimal in647 teraction. As shown in Figure 9, in this case, flow patterns 648 are influenced by the interaction among multiple obstacles, 649 Sánchez-Aguado et al.: Preprint submitted to Elsevier Page 9 of 15