Urban Interference Effects on the Aerodynamic Loads of a CAARC High-Rise Building: A Comparative RANS and DDES Study

1. Introduction

Urbanization and technological advances in construction have transformed cities worldwide, increasing population density and promoting the “vertical city” concept through the proliferation of skyscrapers. This evolution has become a major concern in engineering design, as highlighted by Caswell et al. [1], who demonstrated the influence of urban geometry on thermal and visual comfort in high-rise districts. Beyond representing economic development, tall buildings also optimize land use and encourage public transportation in dense urban centers [2]. Nevertheless, increasing building heights and urban densification affect the urban environment by modifying wind flow patterns and pollutant dispersion, while also influencing thermal comfort, urban energy efficiency and intensifying urban heat island effects, as reported by Siddique et al. [3].

The Atmospheric Boundary Layer (ABL) constitutes the region adjacent to the Earth’s surface where the flow is directly influenced by surface friction, thermal stratification resulting from solar radiation, and the planet’s rotational effects, whereas within cities, research specifically addresses the Urban Canopy Layer (UCL), which extends from the ground level to a region immediately above the building tops, thereby encompassing the zone of influence characterized by flow separation and the formation of turbulent wakes. The differences between the ABL and the UCL, along the three primary flow regions are illustrated in Figure 1.

According to Chen et al. [4], ventilation within the UCL plays a key role in pollutant dispersion and thermal regulation, being influenced by geometric aspect ratios, packing density, and wind direction. These factors also drive incoming turbulence and aerodynamic interference, modifying the wind loads acting on building facades. Therefore, Wind Engineering in urban environments requires a multidisciplinary approach that accounts for the strong interaction between geometric and aerodynamic parameters.

The assessment for Wind Engineering studies necessitates the application of specific methodologies, which can be broadly classified into three primary approaches:

• In situ measurements, that utilize environmental monitoring devices. Characterized by high fidelity but constrained by high operational costs, time investments required to gather meaningful datasets [5,6].
• Physical testing in reduced-scale wind tunnels allow complex aerodynamic effects to be reproduced under controlled laboratory conditions but it limited by relatively high costs [7].
• Numerical simulations via Computational Fluid Dynamics (CFD), which offer operational flexibility, but their predictive accuracy remains dependent on the appropriate selection of turbulence models, the experimental validation of numerical results, and available High-Performance Computing (HPC) resources [7,8].

CFD modeling has consolidated as a tool for investigating urbanization effects on natural ventilation, atmospheric flow and wind-induced actions on structures, employing diverse approaches that range from simplified models with idealized geometries and urban canyons to simulations of complex and expansive urban areas:

• Hui et al. [9] investigated the influences of trees planted on building facades on the aerodynamics of tall buildings using Large-Eddy Simulations (LES) turbulence model, finding out that vegetation compresses the separation bubble and reduces mean and fluctuating pressures on the side walls, concluding that planting trees on balconies attenuates lift fluctuations and reduces wind loads on the building.
• Silva et al. [10] assessed Pedestrian Wind Comfort (PWC) within an array of nine buildings characterized by variations in height and wind incidence angle using the Reynolds-Averaged Navier-Stokes (RANS) model with the classic $k-\epsilon$ formulation, finding that wind configurations aligned with the canyons generate dangerous velocities, whereas oblique angles mitigate peak values.
• Villalobos-García et al. [11] investigated the effects of surrounding structures on the wind loads acting upon a high-rise building using a RANS model with the $k-\omega$ Shear Stress Transport (SST) formulation, finding that nearby configurations reduce global pressure coefficients by an average of 27% but generate localized suction fields on windward faces.
• Ibarra-Hernández et al. [12] investigated turbulent flow patterns within a real low-rise urban landscape in Monterrey (Mexico) using the RANS model with the Renormalization Group (RNG) $k-\epsilon$ closure equations, finding that mean velocity fields remain independent of the Reynolds number and, consequently, their evaluation demonstrates that main avenues concentrate high velocities and promote turbulence generation.

The selection of the modeling strategy depends on the equilibrium between representative fidelity and computational feasibility, given that large-scale urban simulations can demand meshes with a large number of cells, which renders the utilization of conventional computational resources unfeasible [13].

CFD simulations enable the evaluation of parameters related to wind comfort and safety, provided that appropriate turbulence models and boundary conditions are implemented. While Large-Eddy Simulations (LES) models exhibit high capability in reproducing complex urban flows, hybrid approaches such as Detached-Eddy Simulations (DES) combine the advantages of RANS and LES formulations, providing improved predictions of boundary layer separation, vortex formation and aerodynamic forces [14]. In Brazil, the evaluation of wind effects is based on the Associação Brasileira de Normas Técnicas (ABNT) NBR 6123/2023 [15] standard code, yet CFD has consolidated as an alternative capable of estimating aerodynamic pressures and forces with accuracy comparable to wind tunnel testing, while offering enhanced flexibility and lower costs.

This work is organized as follows: Section 2 presents a review of the literature concerning aerodynamic forces and turbulence modeling. Section 3 describes the adopted methodology, including the main equations, turbulence models, numerical methods, computational domain, boundary conditions and mesh discretization procedure. Section 4 provides numerical validation through the Michelstadt urban neighborhood benchmark, where the predicted velocity fields obtained with different turbulence models and numerical schemes are compared against the experimental data. Section 5 presents and discusses the aerodynamic coefficients for the investigated urban configuration with the two analyzed turbulence models. Finally, Section 6 summarizes the main conclusions, highlighting the limitations of the present study, and outlines recommendations for future research.

2. Literature Review

Property transport within a fluid, including momentum, heat and scalar concentrations, occurs via convection associated with the mean velocity field and through diffusion of either molecular or turbulent nature, whereas turbulence itself is characterized by chaotic, three-dimensional, and unsteady motions that constitute a property of the flow rather than the fluid. The flow regime is determined by the Reynolds number:

$$Re={\displaystyle \frac{\rho \dot{u}L}{\mu }}$$

(1)

where $\rho$, $\dot{u}$ and $\mu$ are, respectively, the fluid specific mass, velocity and dynamic viscosity and $L$ is the characteristic length. Generally, for internal flows in pipes or ducts, particularly circular ones, values below 2300 indicate laminar flow, meaning it is predominantly laminar with perturbations damped by viscosity, whereas for values greater than 4000, the flow becomes turbulent, leading to perturbation amplification and multi-scale vortex formation. In this turbulent regime, the elevated diffusivity enhances the transport of mass, momentum and energy, while Turbulent Kinetic Energy (TKE) is transferred from large to small scales through the energy cascade process until it is dissipated into heat by molecular viscosity.

2.1. Turbulence Modelling

CFD provides a numerical approach to simulate fluid flow, heat transfer, and related phenomena through the resolution of equations of motion. In urban applications, this field has expanded the comprehension of flow patterns across complex geometries while overcoming limitations inherent to experimental techniques [14]. Its engineering applications encompass the assessment of wind potential, pollutant dispersion, urban microclimates, pedestrian comfort, and wind-induced loads on structures. Each of these areas presents challenges related to geometric complexity, multiscale interactions, and turbulence generation induced by buildings, vegetation, and topography, thereby rendering turbulence modeling an important aspect for simulation accuracy.

The primary turbulence strategies are categorized into: (i) RANS formulations, which are grounded in the time-averaging of flow equations; (ii) LES, which explicitly resolve large-scale vortices while modeling sub-grid scales; and (iii) hybrid frameworks, such as DES, which combine RANS formulations in near-wall regions with LES treatments in detached flow zones to achieve an equilibrium between computational cost and physical fidelity. Each modeling approach involves a direct trade-off between predictive accuracy and computational resource allocation, such that RANS formulations remain the most widely implemented in engineering practices due to their computational efficiency and capacity to adequately represent mean flow characteristics, although they depend on closure equations to account for turbulent transport effects [16]. The main closure models are based on the $k-\epsilon$ and $k-\omega$ formulations, which differ with respect to their transport equations, treatment of turbulence anisotropy, and near-wall effect modeling [17,18,19].

The RANS equations describe turbulent flows through the Reynolds decomposition, which separates each instantaneous quantity into a mean and a fluctuating component [20], a formulation that introduces the Reynolds stresses responsible for the closure problem and consequently requires additional equations to represent their transport effects:

$$\rho \dot{{\overline{u}}_j}{\displaystyle \frac{\partial \dot{{\overline{u}}_i}}{\partial x_j}}=\rho {\overline{f}}_i+{\displaystyle \frac{\partial }{\partial x_j}}\left[-\overline{p}{\delta }_{ij}+\mu \left({\displaystyle \frac{\partial \dot{{\overline{u}}_i}}{\partial x_j}}+{\displaystyle \frac{\partial \dot{{\overline{u}}_j}}{\partial x_i}}\right)-\dot{\rho \overline{u_i^{\text{'}}{u}_j^{\text{'}}}}\right]$$

(2)

where $\dot{\overline{u}}$ are the mean velocity components, $x_j$ are the cartesian coordinates, ${\overline{f}}_i$ are the mean body forces, $\overline{p}$ is the mean pressure, ${\delta }_{ij}$ is the Kronecker delta, and $\dot{\rho \overline{u_i^{\text{'}}{u}_j^{\text{'}}}}$ are the Reynolds stresses. For Wind Engineering purposes, $\rho$ is adopted as 1.225 kg/m³, representing air density at 15 °C at sea level [21].

Among these additional equations, the $k-\omega$ SST scheme combines the advantages of the standard $k-\epsilon$ and $k-\omega$ formulations [22,23] by solving two transport equations: one for the TKE ($k$) and another for the specific dissipation rate ($\omega$):

$${\displaystyle \frac{D}{Dt}}\left(\rho \omega \right)=\mathbf{\nabla }·\left(\rho D_{\omega }\mathbf{\nabla }\omega \right)+{\displaystyle \frac{\rho \gamma G}{\nu }}-{\displaystyle \frac{2}{3}}\rho \gamma \omega \left(\mathbf{\nabla }·\dot{\mathit{u}}\right)-\rho \beta {\omega }^2-\rho \left(F_1-1\right)C{D}_{k\omega }+{\mathsf{\Phi }}_{\omega }$$

(3)

$${\displaystyle \frac{D}{Dt}}\left(\rho k\right)=\mathbf{\nabla }·\left(\rho D_k\mathbf{\nabla }k\right)+\rho G-{\displaystyle \frac{2}{3}}\rho k\left(\mathbf{\nabla }·\dot{\mathit{u}}\right)-\rho {\beta }^{*}\omega k+{\mathsf{\Phi }}_k$$

(4)

where $D/Dt$ represents the material derivative operator represents both local temporal variation and convective transport, $\nu$ denotes the kinematic viscosity, ${\nu }_t$ is the turbulent kinematic viscosity, $\mathsf{\Phi }$ is the invariant of the fluid strain rate tensor, $G$ is the TKE production term arising from mean velocity gradient shear, $D_k$ and $D_{\omega }$ represent the effective diffusivity coefficients for TKE and specific dissipation, respectively, $C{D}_{k\omega }$ is a cross-diffusion term enabling the transition between $k-\omega$ and $k-\epsilon$ formulations, and ${\mathsf{\Phi }}_k$ and ${\mathsf{\Phi }}_{\omega }$ are additional source terms. The empirical constants $\beta$, ${\beta }^{*}$ and $\gamma$ serve model closure and calibration purposes, whereas $F_1$ and $F_{23}$ function as blending factors. The kinematic turbulent viscosity is formulated to limit shear stresses according to:

$${\nu }_t={\displaystyle \frac{a_{1}k}{\mathrm{máx}\left(a_1\omega ,b_1{F}_{23}\mathsf{\Phi }\right)}}$$

(5)

where $a_1$ and $b_1$ are empirical constants and $F_{23}$ is a blending function. This formulation enables ${\nu }_t$ to limit shear stresses and to improve predictive capabilities in regions exhibiting adverse pressure gradients.

The DES formulation combines the advantages of RANS and LES models within a unified approach, implementing RANS in near-wall regions and LES in separated flow zones [24], a transition achieved by replacing the wall distance with a modified length scale:

$$d=\min \left(d,C_{DES}\mathsf{\Delta }\right)$$

(6)

Which is dependent on grid cell dimensions. $C_{DES}$ is a calibration constant and $\mathsf{\Delta }$ denotes the maximum dimension of the mesh cell, a modification that causes the model to behave as a RANS near solid boundaries ($d<C_{DES}\mathsf{\Delta }$) and analogously to a Smagorinsky LES model in far-field regions, thereby making scale dependency a local property.

The original DES formulation can exhibit premature transition to LES mode when excessive grid refinement occurs within the boundary layer, artificially reducing turbulent viscosity and inducing spurious flow separation. To overcome this limitation, the Delayed Detached-Eddy Simulations (DDES) variant was developed by Spalart et al. [25], introducing a shielding function to preserve RANS behavior across wall-bounded attached flow regions.

Among these variants, the $k-\omega$ Shear Stress Transport Delayed Detached-Eddy Simulations (SSTDDES) model stands out by utilizing the turbulence length scale of the $k-\omega$ SST model to govern the transition between RANS and LES modes, providing enhanced predictive fidelity for vortex shedding and separated flows at lower computational cost than full LES simulations [26].

The LES formulation in far-field regions relies on Kolmogorov’s similarity theory [27], according to which large-scale vortices depend directly on flow geometry, whereas smaller scales exhibit universal isotropic behavior. Consequently, LES explicitly resolves large turbulent structures while modeling only sub-grid scales via a Sub-Grid Scale (SGS) model. The scale separation is executed through a spatial filtering operation, computationally represented by the grid resolution itself, yielding the filtered mass and momentum conservation equations:

$${\displaystyle \frac{\partial {\overline{\dot{u}}}_i}{\partial x_i}}=0$$

(7)

$${\displaystyle \frac{\partial {\overline{\dot{u}}}_i}{\partial t}}+{\displaystyle \frac{\partial \left({\overline{\dot{u}}}_i{\overline{\dot{u}}}_j\right)}{\partial x_j}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+\nu {\displaystyle \frac{{\partial }^2{\overline{\dot{u}}}_i}{\partial x_j\partial x_i}}-{\displaystyle \frac{\partial {\tau }_{ij}}{\partial x_j}}$$

(8)

where ${\overline{\dot{u}}}_i$ are the filtered velocity components, and ${\tau }_{ij}$ is the sub-grid stress tensor, which arises to introduce the closure problem. In the adopted SSTDDES formulation, the SGS stress tensor is modeled by the Smagorinsky scheme [28].

2.2. Wind Effects in Urban Engineering

Wind actions induce loads on tall and slender buildings, since they are driven by asymmetric pressure distributions resulting from geometry, different wind incidence angles, and surrounding aerodynamic interference effects, rendering the analysis of aerodynamic and aeroelastic effects essential, particularly regarding transverse and torsional forces that frequently determine structural dimensioning and serviceability criteria [29].

Design standards differ in their definitions of basic wind velocity and implemented correction factors, given that while the Brazilian ABNT NBR 6123 [15] and American ASCE 7-22 [30] standards utilize a 3-second gust velocity, Eurocode 1 [31] considers a 10-minute mean velocity and the Chinese GB 50009/2012 [32] code employs specific reference pressure and height parameters. Nevertheless, despite these distinct methodological approaches, prescriptive codes can be insufficient for high-rise buildings or complex geometries, thereby necessitating complementary evaluation via wind tunnel testing or CFD frameworks.

Within urban regions, surrounding buildings modify the flow field and consequently wind-induced actions, where despite the introduction of neighborhood coefficients by the standard codes [15,30,31,32], its application remains limited to prescriptive estimates that do not adequately account for complex flows or pedestrian comfort parameters, thereby requiring specific investigations that evaluate the influence of urban geometry, relative building placement, and wind incidence angles [33]. According to Blessmann [34], the primary interference phenomena encompass buffeting, which is induced by periodic vortex shedding mechanisms; the Venturi effect, responsible for flow acceleration between adjacent structures; and wake turbulence, which alters the wind loads acting upon downstream configurations through recirculation regions and vortex formation.

CFD has consolidated as a tool for building aerodynamics analysis, being implemented in the evaluation of wind loads, natural ventilation, urban microclimates, and the optimization of urban layouts for energy efficiency [35], while additionally enabling the investigation of urban geometry influence on flow patterns and pedestrian safety [36]. To guarantee the reliability of the numerical results, the implementation of best practices guidelines is required. These include adequate mesh resolution, the appropriate selection of the turbulence model, and the systematic validation of numerical data against experimental benchmarks.

3. Methodology

This section details the methodology, which encompasses the evaluation of aerodynamic coefficients and their computational formulation within OpenFOAM, the simulation domain, boundary conditions, near-wall treatments, and the implemented numerical schemes. All simulations were carried out in OpenFOAM v2206, using HPC resources from Centro Nacional de Desempenho de Alto Desempenho em São Paulo (CENAPAD-SP) at Universidade Estadual de Campinas (UNICAMP). The specific HPC resources employed for the simulations consisted of two AMD EPYC 7662 processors based on the Zen 2 architecture, each featuring 64 cores, for a total of 128 cores, with 512 GB of available RAM and 328 TFLOPS distributed across 58 compute nodes. Since each simulation utilizes a single compute node, approximately 5.65 TFLOPS are employed for each case.

3.1. Aerodynamic Coefficients

The quantification of wind forces on buildings is performed through aerodynamic coefficients, which represent dimensionless parameters relating forces and moments to flow properties and structural geometry. The mean pressure coefficient $C_p$ describes the surface pressure relative to the freestream static and dynamic pressures $p$, according to:

$$C_p={\displaystyle \frac{p-p_{\infty }}{0.5\rho {\dot{u}}_{\infty }^2}}$$

(9)

where $p_{\infty }$ = 0 Pa and ${\dot{u}}_{\infty }$ = 1 m/s are, respectively, the reference pressure and velocity of the undisturbed freestream flow.

Integrating these pressure coefficients over the surface of the main building yields the overall aerodynamic forces and moments acting on the building:

$${\mathit{F}}_p=\sum _i{\rho }_i{\mathit{s}}_{f,i}\left(p-p_{\infty }\right)$$

(10)

$${\mathit{F}}_v=\sum _i{\mathit{s}}_{f,i}·\left(\mu {\mathit{R}}_{dev}\right)$$

(11)

where ${\mathit{F}}_p$ and ${\mathit{F}}_v$ are, respectively, the normal pressure and the tangential viscous forces, ${\mathit{s}}_{i,f}$ is the face area vector, and ${\mathit{R}}_{dev}$ corresponds to the deviatoric stress tensor. These quantities are expressed as the drag ($C_d$), side ($C_s$), and lift ($C_l$) forces, as torsional ($C_t$) and overturning ($C_{mx}$ and $C_{my}$) moments:

$$C_d={\displaystyle \frac{F_d}{0.5\rho {\dot{u}}_{\infty }^2{S}_n}}$$

(12)

$$C_s={\displaystyle \frac{F_s}{0.5\rho {\dot{u}}_{\infty }^2{S}_m}}$$

(13)

$$C_l={\displaystyle \frac{F_l}{0.5\rho {\dot{u}}_{\infty }^2{S}_n}}$$

(14)

$$C_t={\displaystyle \frac{T_o}{0.5\rho {\dot{u}}_{\infty }^2{S}_{h}b}}$$

(15)

$$C_{mx}={\displaystyle \frac{M_x}{0.5\rho {\dot{u}}_{\infty }^2{S}_{m}h}}$$

(16)

$$C_{my}={\displaystyle \frac{M_y}{0.5\rho {\dot{u}}_{\infty }^2{S}_{n}h}}$$

(17)

where $S_n=bw$, $S_m=bh$ and $S_h=wh$ represent the reference areas (where the normal area to the force direction is considered), and $L_{ref}=h,w$ or $b$ donates the reference length associated with the rotational axis of the moment, as illustrated by Figure 2 showing the effects of all these coefficients on a representative building, in which $b$, $w$ and $h$ correspond to depth, width and height dimensions, respectively.

3.2. Simulation Domain

The computational domain and boundary conditions were defined based on the works of Aboshosha et al. [37], Elshaer et al. [38,39], and Liu et al. [40], adopting dimensions of 20$L$×10$L$×6$L$, with the building arrangement positioned at 5$L$ from the inlet, 15$L$ from the outlet, 5$L$ from the lateral boundaries, and 6$L$ from the top boundary, adhering to the best practices guidelines for CFD simulations in urban environments proposed by Franke et al. [41] and Tominaga et al. [42]. Figure 3 shows a perspective view of this computational domain.

The mesh is generated through a two-step process, wherein the blockMesh utility initially builds the background mesh with progressive refinement toward the building region to concentrate elements near the model and ground, followed by snappyHexMesh, which executes local refinement and geometric fitting via castellation and snapping stages.

The building arrangement selected for this study corresponds to the Oklahoma City (USA) downtown cluster dataset provided by the University of Hamburg (Germany) through the European Cooperation in Science and Technology (COST) framework [43]. The selection of this specific urban arrangement is primarily justified by its adoption as a benchmark model in urban aerodynamics and microclimate studies, since high-fidelity three-dimensional geometric datasets of real cities are often proprietary or require costly topographical surveys. The model features an open space between structures near the central location, where two distinct scenarios are evaluated regarding the aerodynamic analysis: the isolated configuration, and the insertion in a previously empty space of the Commonwealth Advisory Aeronautical Research Council (CAARC) ($b$ = 45.72 m; $w$ = 30.48 m; $h$ = 182.88 m) standard tall building model [44]. Figure 4 details the mesh condition in the model before and after the insertion of the high-rise structure. This comparison aims to identify how aerodynamic coefficients and, consequently, the structural performance of a high-rise building are influenced by neighborhood interference effects across different wind incidence angles.

3.3. Boundary Conditions

For the RANS $k-\omega$ SST model, the boundary conditions were specified to reproduce a neutral ABL, in accordance with Hargreaves and Wright [45], Yang et al. [46], and Richards and Norris [47]. The friction velocity is calculated by:

$${\dot{u}}^{*}\left(z\right)={\displaystyle \frac{u_{\infty }\kappa }{\ln \left(z/z_0\right)}}$$

(18)

Where $\kappa$ = 0.41 is the von Kárman constant, $z$ is the height above ground level and $z_0$ = 0.03 m is the atmospheric roughness parameter [48]. The logarithmic profile along with the distributions of TKE ($k$) and specific dissipation rate ($\omega$) are defined by:

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

(19)

$$k={\displaystyle \frac{{\left(u^{*}\right)}^2}{\sqrt{C_{\mu }}}}\sqrt{C_1\ln \left({\displaystyle \frac{z+z_0}{z_0}}\right)+C_2}$$

(20)

$$\omega ={\displaystyle \frac{u^{*}}{k\sqrt{C_{\mu }}\left(z+z_0\right)}}$$

(21)

Where $C_{\mu }$ = 0.09, $C_1$ and $C_2$ are empirical constants to adjust the functions to the logarithmic profile. At the inlet, velocity is prescribed as Eq. (19), $k$ is set to 0.00375 m²/s², $\omega$ = 0.061 s^-1 and ${\nu }_t$ is fixed as zero. At the outlet, zero-gradient conditions are imposed on velocity and transport variables, while pressure is fixed at zero. The top and lateral boundaries are modeled as symmetry planes, thereby eliminating blockage effects and the development of artificial boundary layers.

For the SSTDDES model, the steady-state solution obtained from the RANS simulation serves as the initial condition, maintaining identical boundary conditions for scalar transport quantities. The distinction lies in replacing the inlet velocity condition with the Divergence-Free Synthetic Eddy Method (DFSEM) [49], implemented in OpenFOAM as the turbulentDFSEMInlet utility, generating synthetic turbulence by convecting eddies through a virtual inlet volume, utilizing the mean fields from the prior RANS simulation as a foundation for turbulent fluctuations. This approach yields a transient inlet condition that represents atmospheric turbulence with higher physical fidelity, albeit at a high computational expense relative to RANS simulations.

Finally, to ensure the reliability of the aerodynamic evaluations, appropriate convergence criteria and simulation periods were established. For the steady-state RANS formulation, the solver was executed for 900 pseudo-time iterations. A preliminary analysis indicated that 600 iterations were required to achieve a stable flow profile, so the final 300 iterations were utilized for data extraction. To maintain consistency for transient evaluations, the SSTDDES model utilized this converged steady-state flow field as its initial condition, being subsequently executed for 5 minutes (300 seconds) of actual physical time to adequately capture the transient turbulent phenomena.

3.4. Numerical Schemes

For the RANS $k-\omega$ SST model, the simpleFoam solver is employed, which is suitable for steady-state incompressible flows, using relaxation factors of 0.3 for pressure and 0.7 for all other variables, alongside a convergence tolerance of 10^-4. Temporal integration is handled by the steadyState scheme, whereas gradients, divergence and Laplacian terms are discretized using, respectively, Gauss linear with cellLimited and linearUpwind for velocity; limitedLinear for turbulent transport variables; and Gauss linear corrected. Pressure is solved using the Geometric-Algebraic Multi-Grid (GAMG) solver, while velocity, TKE and specific dissipation rate utilize smoothSolver with Gauss-Seidel smoothing.

For the SSTDDES model, the pimpleFoam solver is used for transient flows, incorporating two PIMPLE outer correctors per time step. Temporal integration employs the second-order implicit backward scheme, gradients are computed using least squares, convective terms rely on the Gauss limitedLinear scheme, and Laplacian terms are discretized via Gauss linear orthogonal. Pressure is solved using GAMG with DICGaussSeidel preconditioner, whereas transport equations utilize smoothSolver paired with a symGaussSeidel preconditioner.

A constant pseudo-time step (1 iteration = 1 second) is adopted for simpleFoam to drive steady-state convergence, whereas pimpleFoam utilizes an actual physical time step of $\mathsf{\Delta }s$ = 0.01 s. The selection of these numerical schemes aligns with established practices in the literature for building aerodynamics and urban flow simulations in OpenFOAM [50,51,52].

3.5. Near-Wall Treatment

To ensure proper modeling of the turbulent boundary layer without excessively refined mesh at the wall, wall functions are applied to the TKE, specific dissipation rate, and turbulent viscosity fields. For the velocity, a no-slip condition ($\dot{u}$ = 0 m/s) is imposed on all solid boundaries. For TKE, the kqRWallFunction boundary condition is specified, acting as a zero-gradient wrapper in high-Reynolds-number flows. Conversely, for ${\nu }_t$ and $\omega$, a continuous transition treatment is performed based on the dimensionless wall distance $y^{+}$.

The omegaWallFunction computes the value of $\omega$ near-wall cells by considering its asymptotic contributions, where its value within the viscous layer is given by:

$${\omega }_{vis}={\displaystyle \frac{6\nu }{{\beta }_1{y}^2}}$$

(22)

Where ${\beta }_1$ = 0.075 is a model constant in $k-\omega$ closures and $y$ is the cell-center normal distance to the wall. In the logarithmic sublayer, equilibrium between turbulence production and dissipation is given by the equation below.

$${\omega }_{log}={\displaystyle \frac{k^{1/2}}{C_{\mu }^{1/4}\kappa y}}$$

(23)

Regarding the turbulent viscosity, the nutkWallFunction condition is applied, which evaluates this parameter based on local TKE and $y^{+}$:

$$y^{+}={\displaystyle \frac{u^{*}y}{\nu }}={\displaystyle \frac{C_{\mu }^{1/4}y{k}^{1/2}}{\nu }}$$

(24)

$${\nu }_t=\nu \left({\displaystyle \frac{y^{+}\kappa }{\ln \left(E_v{y}^{+}\right)}}-1\right)$$

(25)

where $E_v$ = 9.8 denotes the roughness parameter for smooth walls. If the cell lies within the viscous layer, the turbulent viscosity is set directly to zero, yielding no additional turbulent contribution to the effective viscosity at the node. The target $y^{+}$ value should fall within the range of 30 to 300, preferably closer to 30 along solid boundaries, to preserve the physical fidelity of the numerical models [53].

4. Validation

To evaluate model fidelity regarding flow field behavior, the numerical predictions from the RANS $k-\omega$ SST and SSTDDES models are compared against experimental data from the Michelstadt dataset, provided by COST Action ES1006 framework [54]. As illustrated in Figure 5, the experimental setup models an urban environment consisting of structures of varying heights with sharp edges at all vertices, designed to emulate prominent boundary layer separation within the UCL.

In the wind tunnel experiment, physical prober and LIDAR technology were used to measure the velocity field within a central 0.7×0.7 m area of the array (Figure 6a) across five elevations ($z$ = [0.89; 4; 8; 12; 13.33] cm). Results from the RANS and SSTDDES simulations were extracted at the corresponding locations for validation, but due to the implementation of wall functions, data points for which $\dot{u}/{\dot{u}}_{\infty }$ < 0.2 were excluded because they were associated with spurious near-wall predictions along solid surfaces, thereby restricting the quantitative comparison to street canyons (Figure 6b) and regions above the roof level. The Michelstadt dataset is widely implemented for the validation of numerical models predicting flow fields and pollutant dispersion in Architectural and Wind Engineering applications [55,56,57].

To compare the experimental and numerical data, five statistical metrics were evaluated: the coefficient of determination (R²), the fraction of predictions within a factor of two observations (FAC2), the fractional bias (FB), the root-mean-square error (RMSE) and the mean absolute error (MAE). To indicate a good agreement with benchmark literature, R² and FAC2 should be close to 1, while FB, RMSE and MAE should approach zero. Results are presented in scatter plots for both models across three mesh refinement levels: M1, M2 and M3 (near-wall mesh refinement details are illustrated in Figure 7), containing approximately 1.82, 2.76 and 3.96 million finite volumes, respectively, at a Reynolds number of approximately $Re$ = 143,000.

Scatter plots comparing numerical and experimental velocity fields results from the RANS $k-\omega$ SST model are presented in Figure 8, with corresponding statistical metrics and $y^{+}$ values summarized in Table 3. For the SSTDDES model, the comparison is presented in Figure 9, with corresponding metrics and values for $y^{+}$ summarized in Table 4. Computational execution times were normalized relative to the RANS-M3 baseline, which reached 1 hour and 39 minutes of processing time.

For the RANS $k-\omega$ SST model, analysis of the metrics and scatter plots demonstrates a consistent improvement in accuracy as the mesh is refined. The model exhibited positive sensitivity to increased grid resolution, as mesh M3 achieved the best performance, with a coefficient of determination R² = 0.991, indicating excellent agreement with experimental data and surpassing obtained on coarser grids.

This performance gain is further supported by the reduction in error metrics, where RMSE decreased from 0.087 to 0.046, and MAE dropped by approximately 47.9%. Visually, the scatter plots follow this trend: while M1 and M2 exhibit scatter relative to the ideal line, the data points for M3 are closely concentrated along the bisector line, especially when $\dot{u}/{\dot{u}}_{\infty }$ > 0.5. Thus, mesh refinement reduced the near-zero discrepancies close to solid boundaries and provided a more faithful representation of the flow field within the street canyons. $y^{+}$ values show that the data remained below the maximum value of 300 and, even though the average values are lower than 30, the statistical coefficients showed that the model was still able to capture wind flow phenomena comparable to the experimental setup, with good agreement in all evaluated cases.

For the SSTDDES model, M1 exhibited the poorest statistical performance among all tested configurations, recording the lowest determination coefficient R² = 0.958 and highest error metrics, with RMSE = 0.105 and MAE = 0.086. This performance degradation occurs because coarser meshes in hybrid formulations lack sufficient spatial resolution in street canyons and shear layers, delaying the development of resolved sub-grid turbulent fluctuations. Furthermore, the M1 scatter plot exhibits higher mean velocity values than the other cases due to the numerical divergence induced by the coarser mesh. Conversely, M2 provided the best overall predictive accuracy among all evaluated cases, achieving R² = 0.996, alongside lowest error indices RMSE = 0.035 and MAE = 0.025, outperforming M3 by a small margin. Such non-monotonic convergence behavior is well documented in computational wind engineering and hybrid modeling literature, where in DES-type formulations, excessive localized grid refinement without full DNS-level resolution can cause the grid to enter an “ambiguous grid density” regime, potentially triggering subtle grid-induced separation or log-layer mismatch near solid boundaries, where RANS shielding and LES sub-grid viscosity transition [25,58,59]. Thus, mesh M2 offers the optimal balance for solving shear layer dynamics without inducing numerical errors. Regarding the $y^{+}$ parameter, values followed the same trend observed in the RANS simulations, with the upper threshold remaining unexceeded, while the statistical metrics confirm strong agreement of the numerical model.

The computational execution time for the SSTDDES simulations was, on average, approximately 35 times higher than that of the RANS cases, while yielding similar results regarding the mean velocity field in both approaches. However, it should be emphasized that although the time-averaged fields are equivalent, the SSTDDES model inherently resolves turbulent fluctuations and localized flow phenomena near structural elements, making it more relevant when instantaneous flow dynamics are required for engineering design applications. The validation confirms that both proposed models capture the velocity field across densely built urban environments, establishing their capability to perform urban-scale simulations for predicting neighborhood interference effects.

5. Results and Discussion

The aerodynamic coefficients calculated at the different incidence angles were grouped into two categories: mean values, obtained by averaging over the steady-state iterations or physical simulation time, and absolute peak values. Each coefficient in each category was presented in a separate graph, encompassing all wind incidence angles. The drag, lateral force, lift, pitching moment of the Y-axis, rolling moment of the X-axis, and torsional moment coefficients are presented, respectively, in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 for both mean and peak values. For the lift and torsional moment coefficients, the reference area and reference length were defined as the building base area (with a 1:1.5 aspect ratio) and its smaller dimension, respectively. For the remaining aerodynamic coefficients, the reference area was defined as the area of the smaller building façade (with a 1:6 aspect ratio), while the building height was adopted as the reference length. The nominal values of the coefficients, along with a time-domain plot for each one of them, are presented in Appendix A of this study.

The coefficients computed for the building reveal a dependency of structural wind load on incidence angles, neighborhood morphology, and the selected turbulence modeling approach. Because the CAARC model features a rectangular base with an aspect ratio of 1:1.5, the wind incidence angle dictates the exposed area and the magnitude of the resulting forces. When the wind strikes at angles of 0° and 180°, it acts directly upon the smaller façade, whereas incidence angles of 90° and 270° impact the larger façade directly and, consequently, these latter orientations generate the highest mean values for both drag and the overturning moment in the Y-axis.

However, the analysis of peak aerodynamic loads demonstrates that maximum forces do not always align with the largest exposed surface areas, since the highest maximum drag coefficient recorded across all simulations occurs at the oblique incidence angles $\theta$ = 135° and 315° under the isolated SSTDDES configuration ($C_d$ = 3.82). Oblique angles also determine the torsional behavior of the building, since in isolated scenarios for both turbulence models, the mean torsion coefficient remains near zero when the wind is perpendicular to the facades. In these cases, the lack of neighborhood interference means the direct flow does not induce tangential forces along the building faces. Conversely, when the flow approaches at oblique angles, the wind acts simultaneously on two adjacent faces, creating pressure imbalances and flow separations that result in more pronounced mean torsional effects. This behavior is consistent with the wind engineering literature, since experimental findings of Guzmán-Solis et al. [60] show that critical torsion occurs at oblique wind directions rather than when wind is normal to the faces. It should be noted, however, that the SSTDDES peaking torsional values over the RANS $k-\omega$ SST model may be a case-specific result from the dataset, and indicates that further data with different building shapes and aspect ratios should be analyzed.

Furthermore, a comparative analysis between the RANS $k-\omega$ SST and the SSTDDES models underscore disparities inherent to turbulence modeling, particularly regarding transient aerodynamic phenomena, While RANS provides computational efficiency and adequately represent time-averaged mean flow characteristics, they inherently filter out transient turbulent fluctuations. The convergence plots in Appendix A show that, while RANS model exhibits numerical variation in the aerodynamic coefficients at the beginning of the calculation, the forces generally stabilize after approximately 300 iterations as the solver reaches a steady-state response. In contrast, the SSTDDES plots show frequent changes in the behavior of the aerodynamic forces, with significant fluctuations, particularly in the drag and lateral forces, causing the RANS $k-\omega$ SST model to slightly overestimate certain mean coefficients when compared to the SSTDDES approach. Specifically, the $k-\omega$ SST calculates a mean drag that is 17.91% higher in the isolated configuration and 8.23% higher on average in the neighborhood configuration when compared to the SSTDDES mean predictions. Mean lift and moment Y also exhibit similar overestimations in the RANS model across configurations.

In contrast, the SSTDDES formulation explicitly resolves large-scale vortices and sub-grid turbulent fluctuations in separated flow zones, allowing it to capture the extreme aerodynamic peaks. For example, Figure A3 and Figure A4 in Appendix A show that drag and moment Y have high peaks compared to the time history for all turbulence models around $t$ = 250 seconds.

In 50. 08% greater on average than the RANS prediction, and this amplification is observed across all parameters: on average, 202.15% higher for side forces, 51.87% higher for moment Y, 226.80% higher for moment X, and 204.73% higher for torsion compared to $k-\omega$ SST peak values. These divergences persist in the neighborhood configuration, where SSTDDES predicts a peak drag 48.68% higher and peak side forces 143.18% higher on average than $k-\omega$ SST. These results are consistent with previous studies showing that hybrid approaches improved predictive fidelity for boundary layer separation and vortex formation and, for structural design purposes, relying solely on RANS models, can underestimate important peak and fluctuating load effects, especially in separated-flow regions, reinforcing the necessity of models that can capture turbulent fluctuations to evaluate peak structural forces [25,61,62].

The introduction of a surrounding neighborhood alters the wind flow patterns, inducing aerodynamic interference effects on the structure. In terms of mean aerodynamic forces, the urban canopy generally acts as a shielding mechanism, creating a protective zone that mitigates the magnitude of most global coefficients: considering the average values of all evaluated wind incidence angles, mean drag is 17.91% higher in the isolated $k-\omega$ SST scenario and 20.43% higher in the isolated SSTDDES scenario when compared to their respective neighborhood configurations. Similarly, the mean lift force is 23.58% greater in the isolated RANS case and 20.26% greater in the SSTDDES case, demonstrating that upstream structures reduce the direct aerodynamic loads acting on the facades of the target building.

However, despite this overall reduction in the mean coefficients, the neighborhood configuration introduces some interference phenomena that amplify transient peak forces. While the mean lift experiences a reduction, the peak lift values are amplified by the surrounding urban geometry, measuring 16.47% higher in the RANS configuration compared to the isolated case. Torsional behavior further exemplifies this wake-induced phenomenon, since in the RANS isolated models, maximum torsion values remain close to zero. Yet, the presence of the neighborhood introduces asymmetric flow separation and turbulent wakes that amplify torsional effects across all angles, resulting in peak torsion values that are 21.52% higher than those in the isolated setup, indicating that the adjacent buildings induce twisting loads that do not exist when the building stands alone. This type of behavior, which is specific to each region, topography, and morphology of the building distribution within a city, cannot be determined in a generalized manner through normative approaches, demonstrating that numerical and experimental tests with specific setups may be necessary for urban planning and the design of new buildings.

In the isolated configuration, the peak torsional coefficient shows different angle-dependent trends in the RANS and SSTDDES results. This divergence is consistent with modeling strategies: steady RANS responds primarily to the mean asymmetric pressure field induced by oblique wind incidence, whereas SSTDDES with turbulent inflow also captures instantaneous vortex shedding and fluctuating pressure asymmetry. For this reason, torsional peaks in SSTDDES do not need to follow the same monotonic variation with wind angle seen in RANS, as the literature supports linking this behavior to both turbulence-model formulation and incoming turbulence, since changes in angle of direct incidence and inflow turbulence can modify separation, shedding mode, and fluctuating torsional loading [60,63,64].

The most notable deviation from the general shielding effect is observed in the lateral dynamics which, contrary to drag and vertical forces, both mean and maximum side forces are amplified by the neighborhood across the majority of evaluated angles, peaking in the SSTDDES neighborhood configuration at an incidence angle of $\theta$ = 180°, where the wind strikes directly the smaller facade, reaching a maximum coefficient of $C_s$ = 2.59. At this specific angle, upstream wake interactions and localized negative pressure fields generated by adjacent structures likely exacerbate lateral shedding and crosswind instability.

Ultimately, these results indicate that while urban densification may provide a general shielding effect against mean longitudinal wind loads, it concurrently generates highly turbulent, localized wakes that amplify transient peak forces in both evaluated turbulence models, with side forces and torsional moments showing the largest differences between isolated and neighborhood configurations. These amplifications underscore the necessity of modeling specific urban surroundings to ensure structural safety and adequately capture wind effects for the potential structural design of the building.

Furthermore, to assess a deeper understanding of the numerical simulations, the aerodynamic coefficients obtained for the isolated configuration were compared with data from literature. However, a particularity adopted in the present study regarding reference quantities should be highlighted, as the reference area used to calculate the drag coefficient was kept constant and equal to the smaller face of the building (proportions 1:6). Since literature frequently uses the larger projected face of the model as the reference area (proportions 1.5:6), the values extracted from previous studies were multiplied by an adjustment factor of 1.5 to allow for an equitable comparison with the results of the present numerical modeling.

When evaluating the scenarios with wind perpendicular to the building faces, a good agreement between the numerical and experimental results was observed, since for $\theta$ = 0°, the numerical models yielded a mean $C_d$ = 1.33 for the $k-\omega$ SST and $C_d$ = 1.22 for the SSTDDES, with wind tunnel experiments reported in Alminhana et al. [65] reporting approximately $C_d$ = 0.80 based on the larger face. After applying the area conversion factor, the adjusted value is 1.20, indicating that the RANS and SSTDDES formulations yielded a relative difference of approximately 9.77% and 1.64%, respectively, compared with the experimental data, demonstrating good predictive agreement in capturing the mean flow topology.

For $\theta$ = 90°, the simulations yielded $C_d$ = 2.14 for RANS and $C_d$ = 2.01 for SSTDDES. The literature indicates that experimental values are sensitive to the intensity of the incoming turbulence, with Elshaer et al. [38] reporting that the transition from uniform flow to a boundary-layer profile increased $C_d$ from 1.49 to 1.80. Yan and Li [66], using the same larger face as the reference area, reported values of 1.50 for smooth flow and 1.80 for turbulent flow, while Alminhana et al. [65] obtained values of 1.45 and 1.24 for smooth and turbulent flow, respectively. By converting the most conservative range in the literature (1.24 < $C_d$ < 1.50) to the reference area adopted in this study, the expected values range from 1.86 to 2.25, with both models agreeing well with the turbulent-flow conditions.

At $\theta$ = 45°, the values were $C_d$ = 1.95 for RANS and $C_d$ = 1.75 for SSTDDES, while experimental data from Alminhana et al. [65] showed a value of $C_d$ = 1.65 after adjusting for the reference area. Thus, the RANS and SSTDDES formulations exhibited, respectively, a deviation of approximately 15.38% and 5.71%, demonstrating that both models can capture the effects observed in experimental tests, although the DDES model is more consistent with what is expected in this case.

Recent studies suggest parameterizing the drag coefficient as a function of building geometry to improve urban canopy models and predictions. Chen et al. [67] demonstrated that cross-sectional geometry alters recirculation zones, indicating that macroscopic estimates based solely on frontal area density ${\lambda }_f$ or plan area density ${\lambda }_p$ may exhibit discrepancies. When evaluating the parametric formulation proposed by Santiago and Martilli [68], which relates drag to plan area density for sparse flows, the $C_d$ values for isolated rectangular geometries typically range from 1.08 to 1.35, with a limiting upper value of 1.85 for dense arrays. The results obtained in the present study are in good agreement with this range of parametric estimates, since for a $\theta$ = 0° incidence angle, corresponding flow impinging on the smaller facade, the models yielded $C_d$ = 1.22 for SSTDDES and $C_d$ = 1.33 for RANS. At $\theta$ = 90°, after correcting the coefficient based on the frontal area effectively projected to the wind, the actual aerodynamic shape coefficients are 1.34 for SSTDDES and 1.42 for RANS, confirming the consistency of the simulations with the parameterization ranges for tall buildings. Furthermore, inspired by the pronounced variations in drag observed for different morphologies, the presented results allow a directional parameterization to be proposed for the CAARC configuration: variations in the wind incidence angle modify the apparent “aerodynamic shape” and recirculation wake, yielding a drag distribution curve that can be directly incorporated into future mesoscale predictive models to represent buildings with a 1:1.5 aspect ratio more accurately than the use of a single static coefficient.

6. Conclusions

This study presented a numerical investigation of the aerodynamic loads on a CAARC high-rise building placed within a dense urban neighborhood in Oklahoma City and, by using CFD simulations, the influence of the surrounding urban geometry and multiple wind incidence angles were examined. The analysis evaluated aerodynamic coefficients, while comparing the RANS $k-\omega$ SST and the SSTDDES turbulence models.

The results demonstrated a sensitivity of aerodynamic loads to urban geometry and wind direction. Analysis revealed that wind incidence angles of 90° and 270° generated the highest mean drag and overturning moments, as the airflow directly impacts the larger facade of the building. Conversely, oblique incidence angles amplified torsional effects and, while the surrounding neighborhood generally provided a shielding effect that reduced mean structural loads, it concurrently induced complex wake interactions that amplified localized transient peak forces, particularly lateral side forces and peak torsion.

Furthermore, the comparison between the turbulence models revealed that steady-state formulations might misrepresent the nature of certain peak loads. While RANS predicted torsional peaks driven by mean geometric asymmetries at oblique angles, the SSTDDES captured instantaneous alternating vortex shedding at normal incidence angles ($\theta$ = 0° and 90°), suggesting that resolving transient fluctuations is important, as mean-flow symmetries might fail to represent the shedding mechanics that drive extreme bluff body twisting. Additionally, the evaluation of the drag coefficients corroborated literature parameterization theories, and for future work, the proposal of directional coefficient distribution curves could offer a more precise and consistent input for future mesoscale predictive models representing buildings with the same base proportions as the CAARC building and other geometries.

Overall, the findings suggest that aerodynamic interactions in high-density urban environments are complex and oversimplified analytical methods or steady-state assumptions may result in the underestimation of extreme peak loads. The spatial arrangement of the buildings, combined with their orientation relative to prevailing winds, must be considered during early-stage urban design and new buildings development. This study underscores the value of CFD-based wind engineering analysis to ensure structural safety.

Limitations of the Research and Recommendations for Future Work

Despite the findings provided by this investigation, certain limitations should be acknowledged and addressed in future research. Regarding turbulence modeling, this study focused on the comparison between RANS and SSTDDES formulations; however, a full Large-Eddy Simulation (LES) model was not explored. While RANS is adequate for predicting mean velocity fields, it is inherently deficient in capturing turbulence intensity and gustiness. Future studies should conduct a comprehensive comparative analysis involving RANS, DDES and LES frameworks, confronting the computational processing time of each approach with its accuracy and computational viability for structural projects. Also, the current research concentrated on a scenario where the central high-rise prominently stands out above the surrounding urban canopy. Future investigations should examine cases where the target building possesses a height equal or lower than its neighboring structures.

The current study investigates the CAARC model, which features a standard rectangular cross-section with sharp corners, while modern architectural designs could employ aerodynamic shaping, such as chamfered, recessed or round corners to mitigate wind loads. Future research could evaluate how these geometric modifications perform when subjected to complex, highly turbulent wakes of a dense neighborhood, it may alter the aerodynamic effectiveness of such mitigation strategies.

While this research assumed a neutral ABL with isothermal conditions, densely built urban areas are subjected to Urban Heat Island effects, which introduces thermal stratification. Future CFD simulations should incorporate thermal buoyancy effects to evaluate how temperature-driven turbulence modulates boundary layer separation, wake recovery and the resulting aerodynamic forces on high-rise buildings.

Finally, the building model was considered perfectly rigid and, given the transient peak forces captured by the SSTDDES model, particularly the crosswind and torsional loads, future investigations should couple the CFD routine with a Computation Structural Dynamics module. Evaluating fluid-structure interaction could allow the assessment of aeroelastic phenomena, such as vortex-induced vibrations and aerodynamic damping, which are important for fatigue and serviceability analysis of slender buildings.