Applicability of URANS and DES Simulations of Flow Past Rectangular Cylinders and Bridge Sections

This paper discusses the results of computational fluid dynamics simulations carried out for rectangular cylinders with various side ratios of interest for many civil engineering structures. A bridge deck of common cross-section geometry was also considered. Unsteady Reynolds-averaged Navier–Stokes (URANS) equations were solved in conjunction with either an eddy viscosity or a linearized explicit algebraic Reynolds stress model. The analysis showed that for the case studies considered, the 2D URANS approach was able to give reasonable results if coupled with an advanced turbulence model and a suitable computational mesh. The simulations even reproduced, at least qualitatively, complex phenomena observed in the wind tunnel, such as Reynolds number effects for a sharp-edged geometry. The study focused both on stationary and harmonically oscillating bodies. For the latter, self-excited forces and flutter derivatives were calculated and compared to experimental data. In the particular case of a benchmark rectangular 5:1 cylinder, 3D detached eddy simulations were also carried out, highlighting the improvement in the accuracy of the results with respect to both 2D and 3D URANS calculations. All of the computations were performed with the Tau code, a non-commercial unstructured solver developed by the German Aerospace Center.


Introduction
The use of computational fluid dynamics (CFD) is becoming more and more frequent in civil engineering applications.In particular, numerical methods are often used to determine the steady and unsteady aerodynamic loads on structures and structural elements due to wind, although at present, the design process of engineering structures rarely does not include wind tunnel tests.Moreover, structures generally represent bluff obstacles for the airflow, producing massive separation and detachment of eddies.In this case, advanced methods to treat turbulence, such as large eddy simulation (LES), are more appropriate, but they are also very expensive from the computational point of view, requiring fine spatial and temporal discretizations, three-dimensional meshes also for two-dimensional problems (slender prismatic elements) and long simulated time records to obtain statistically-converged results.In addition, civil structures are exposed to turbulent flows, and the proper simulation of oncoming turbulence is another difficult task.Finally, modern engineering structures are very light and slender, so that aeroelastic instabilities may represent a serious problem, and fluid-structure coupling simulations are often necessary.
Hence, it is important to understand whether or not it is possible to simulate with reasonable accuracy the flow around certain types of bluff bodies with relatively economic approaches, such as unsteady Reynolds-averaged Navier-Stokes (URANS) equations in combination with advanced statistical closure models.Furthermore, hybrids of statistical and scale-resolving methods, such as detached eddy simulation (DES), may represent an option, although with a cost much higher than URANS.
The particular attention devoted in this paper to rectangular cylinders is dictated by the fact that civil engineering structures or structural elements often present a rectangular cross-section or can be approximately schematized as elongated rectangular prisms.Furthermore, only one parameter, namely the side ratio, defines the rectangular geometry, which, due to its simplicity, clearly highlights interesting fluid dynamics phenomena and allows a detailed study of the aerodynamic features of the generated flow field.Similar flow characteristics can be encountered also in more complex geometries, but in this case, they are usually more difficult to investigate.Finally, elongated rectangular prisms are prone to different types of flow-induced excitation, depending on the side ratio, which are of general interest in civil engineering structures.However, next to rectangular cylinders, in the second part of the paper, a realistic bridge cross-section is also considered.

Solver
The numerical simulations were performed using a noncommercial finite-volume unstructured flow solver, the DLR-Tau code (e.g., [1]), developed at the German Aerospace Center (DLR).The code solves the compressible Navier-Stokes equations using a vertex-centered approach.Viscous terms were treated using a conventional second-order central differencing scheme, while inviscid fluxes were approximated using a second-order central difference scheme stabilized with scalar artificial dissipation [2,3] (matrix dissipation was also used in a few specific simulations, as explained in Sections 2.2 and 3.4).
The advancing in time was performed through the dual time stepping approach.The time derivatives were first discretized with a second-order backward difference formula, and the resulting sequence of nonlinear steady-state problems was solved through an explicit three-stage Runge-Kutta scheme.In this work, all of the simulations were performed at a low Mach number (M = 0.1), without any algorithm for preconditioning.
In case of simulations with forced motion, the heaving and pitching harmonic vibrations were accounted for, formulating the governing equations in a non-inertial body-fitted moving coordinate system allowing for displacement of the mesh [4].A relative velocity then appears in the convection term of the equations, as well as an additional term in the balance of fluxes.

Turbulence Modeling
In the Reynolds-averaged approach, the unknown variables of the Navier-Stokes equations are statistically averaged.The averaging process results in additional unknown terms, the Reynolds stress tensor components, which have to be modeled.In case of unsteady physical problems, the rationale of this ensemble average of the governing equations is related to the assumption that there is a well-defined spectral gap between the primary low-frequency spectral content of variable fluctuations, which is numerically resolved, and the turbulent content, which is modeled.In the case of vortex shedding past bluff bodies or harmonic forced-motion simulations, this assumption is often acceptable, in view of the dominant spectral contribution at the Strouhal or forced-motion frequencies.
Two different RANS closures were adopted herein.The first one is the one-equation model of Spalart and Allmaras [5], an eddy viscosity model based on the linear Boussinesq relation.According to this hypothesis, the Reynolds stress tensor, analogous to the viscous stress tensor, can be directly related to the flow mean rate of strain.The proportionality constant is the eddy viscosity.The Spalart-Allmaras model proposes a transport equation directly for a quantity representing a modified eddy viscosity.The actual version of the turbulence model used herein is the Spalart-Allmaras-Edwards (SAE) model by Edwards and Chandra [6], characterized by a slight modification of the near-wall treatment to enhance numerical stability and convergence.The Spalart-Allmaras model is known to give accurate results in the case of free-shear flows and attached boundary layers.Similar to the Shear Stress Transport (SST) approach [7], it adheres to Bradshaw's hypothesis (e.g., [8]) and, thus, delivers good predictions for boundary-layer shallow flow separation.It also behaves well in the case of boundary layers in adverse pressure gradients and flow past an airfoil with a blunt trailing edge [5].In the so-called SALSA version (strain-adaptive linear Spalart-Allmaras model), it has also easily been adapted to non-equilibrium conditions [9].Nevertheless, in spite of their numerical advantages, turbulence models based on the standard linear Boussinesq viscosity hypothesis are unable to capture complex interaction mechanisms between Reynolds stresses and the mean velocity field.As a consequence, these models are often inaccurate in the case of strong streamline curvature and massive separation, which characterize the flow past bluff bodies [5,8,10,11].
The second approach followed is an explicit algebraic Reynolds stress model (EARSM), which provides an alternative to the linear Boussinesq relation.It introduces the anisotropy tensor, which offers a more accurate representation of all components of the Reynolds stress tensor and captures rotation and curvature effects, without the drawbacks of high numerical effort and reduced robustness of a full differential stress model.In particular, Rung et al. [12] proposed a quadratic model based on the EARSM of Gatski and Speziale [13] with an improved representation of the production-to-dissipation ratio within the stress-strain relation, which does not require any regularization and yields superior predictive capabilities in case of non-equilibrium large-strain-rates conditions.As compared to the Spalart and Allmaras model, where the turbulence production term is based on the second invariant of the vorticity tensor, this EARSM model relies on a formulation based on the second invariant of both strain rate and vorticity rate tensors, which also does not seem to produce excessive turbulence in the stagnation region in front of the body's upstream face.
In this work, the first-order truncation of this nonlinear EARSM model, called the linearized explicit algebraic (LEA) model, is used in a formulation cast in terms of the standard k − ω background model of Wilcox [14].Despite missing some of the physics covered in the full EARSM equations, this approach has a more physical basis than standard eddy viscosity models.The algebraic equations for the anisotropy tensor components are also characterized by an additional term specific for three-dimensional flows [12,15].Rung et al. [12] tested this model for a 2D transonic bump flow and for the transonic flow over an airfoil.Comparisons to experiments showed good accuracy of the results and significant improvement over the Wilcox k − ω model.Lübcke et al. [11] tested the nonlinear version of this EARSM model for flows past a square and a circular cylinder, showing the significant improvement as compared to Boussinesq viscosity models.In addition to its numerical advantages over the full EARSM models, this model is easily implemented in the existing k − ω standard model.
Detached eddy simulation in the case of a rectangular 5:1 cylinder was also employed as a more advanced strategy of turbulence modeling [16].This non-zonal hybrid method combines the Reynolds-averaged Navier-Stokes approach near solid walls and large eddy simulation away from walls, where significant amounts of turbulent kinetic energy can be economically resolved.The Spalart-Allmaras (SA) model [5] is adopted as an additional equation, and the switch between RANS and LES is obtained through a specific definition of the length scale in the turbulent equation.
The additional parameter with respect to pure URANS simulations is the DES constant (C DES ), which defines the length scale in the SA equation when it works as the LES subgrid scale model.Shur et al. [17] calibrated this parameter against decaying homogeneous isotropic turbulence (DHIT) to ensure that the LES mode returns a proper turbulent kinetic energy cascade.C DES = 0.65 was the recommended value, obtained by adopting a centered fourth-order accurate differencing scheme.A smaller value of the DES constant can be used for more dissipative codes (e.g., [18][19][20]).However, the optimal choice of C DES is closely tied to the numerical properties of the adopted numerical scheme.In the case of scalar dissipation in the term stabilizing the second-order central-difference discretization of the convective term of the mean-flow equations, C DES = 0.45 is assumed on the basis of the calibration using the DHIT reported in [21].Nevertheless, in order to resolve as fine as possible turbulent scales, matrix dissipation was also adopted in the simulations, as scalar dissipation is known to damp out high wavenumber turbulent components.In this case, DHIT calibration suggested using C DES = 0.65.
In 2006, Spalart et al. [22] proposed a new version of DES, called delayed detached eddy simulation (DDES), by modifying the length scale definition in the SA model and, therefore, the RANS/LES switch function, in order to avoid the activation of the LES mode in the case of thick boundary layers or shallow separation regions and, in general, to make the DES results less vulnerable to ambiguous grid densities.Slightly afterwords, improved delayed detached eddy simulation (IDDES) [23] was proposed to resolve the mismatch of the logarithmic layer when DES is used for wall-modeled large eddy simulation.In this work, in combination with matrix dissipation, DDES was employed.Examples of the validation of DES in the framework of the Tau solver can be found in [21,24].
Again, with the purpose not to diffuse the turbulent scales allowed by spatial and temporal resolutions, the skew-symmetric inviscid flux operator was also tested for the discretization of the equations, as it is supposed to minimize the aliasing errors [25].In this case, however, it is important to check that the equation residuals are small and that the continuity equation is sufficiently well satisfied in order to ensure that the skew-symmetric operator performs correctly as an energy-conserving operator.This can require an increase in the number of inner iterations and/or a reduction of the Courant-Friedrichs-Lewy (CFL) number in the dual time stepping algorithm.However, it is worth reminding that second-order schemes are supposed to be less sensitive to aliasing errors than higher-order schemes, as the truncation errors are much larger.
Finally, the convective fluxes in the additional (turbulent) equations were discretized by second-order central differencing for the one-equation Spalart-Allmaras model and by the first-order Roe upwind scheme for the two-equation LEA model.

Spatial and Temporal Discretization
In all of the test cases considered, the governing equations were solved on a circular computational domain with a radius equal to 100 • B, B being the streamwise dimension of the body cross-section.The domain was discretized through hybrid meshes, that is with body-aligned quadrilateral cells near the viscous walls and unstructured triangular cells in the remaining part of the domain.The use of unstructured grids presents several drawbacks: discretization of flow fluxes in the unstructured domain is numerically more diffusive than in the structured domain, and usually, finer spatial discretization is necessary in order to obtain grid-independent results.Nevertheless, the unstructured grids allow one to deal easily with complex geometries, refining the mesh only where it is necessary.For these reasons, the use of unstructured grids is a very popular approach in both research and industrial applications.Naturally, the use of hybrid meshes was not a necessity for the simple geometries investigated in the current paper.However, the current investigation should also be considered as a precursor study for more complex civil and industrial engineering applications.As will be clear in the following, all grids were significantly refined in the wake region.This feature is crucial for the correct resolution of the von Kármán vortex street, i.e., the alternating shedding of eddies typical for separated flows around bluff bodies.
The physical time resolution was chosen in order to discretize the expected period of vortex shedding (estimated either from experiments or preliminary numerical computations) with about 500 time steps, with the exception of the bridge cross-section case, where, however, no less than 220 time steps per period of vortex shedding were employed (see Section 4.1).Although not reported here for the sake of brevity, the satisfactory convergence of the solutions with respect to spatial and temporal discretization has always been verified, with the exception of the 3D DES simulations.For the forced motion test cases (Sections 3.3 and 4.2), the time step size was chosen on the basis not only of the vortex shedding, but also of the period of the driven oscillations.In particular, the latter has always been discretized with no less than 500 time steps.

Boundary and Initial Conditions
Far field boundary condition was imposed at the external contour of the domain, where convective fluxes crossing the boundary faces were calculated employing the theory of Whitfield and Janus [26].
A low-Reynolds boundary condition was imposed at the body surface, applying the turbulent equations down to the wall, where the no-slip condition for momentum was prescribed, without any wall function.This required low-Re grids, i.e., characterized by the distance from the wall of the first layer of grid nodes, such that the wall-unit value in the normal direction y + is of the order of unity and by suitable grid stretching ratios in the direction perpendicular to the wall.In all of the test cases discussed in the following, indeed, it was checked that at any time step, the maximum y + along the body contour (or at the body surfaces, for 3D computations) was smaller than unity.
In case of 3D simulations (Section 3.4), periodic boundary conditions were applied at the spanwise ends of the computational domain.
A small value of the turbulent viscosity (10% of the laminar viscosity) was assigned at the far field boundary, so that the turbulence model was active all over the domain.Nevertheless, no resolved turbulence was generated upstream of the rectangular cylinder, so that simulations were performed with an ideal smooth free stream.As the initial condition, free-stream values of the flow variables were imposed at all grid nodes.

Effect of the Side Ratio
The flow past rectangular sections with side ratios B/D (B and D being respectively the streamwise and cross-flow section dimensions) ranging from 0.2 to 5 and perfectly sharp corners was simulated using a 2D URANS approach in combination with the two-equation LEA k − ω model.The Reynolds number, based on D, was in all cases Re = 26,400, In this first set of simulations, the angle of attack was always null.It is known that rectangular cylinders with side ratios in this range show very different flow fields and mechanisms of vortex shedding.In particular, the flow is fully separated for B/D = 0.2 and 1, while it intermittently reattaches on the side faces for B/D = 5.Moreover, the shear layers separating at the upstream sharp corners do not interfere with the downstream corners for B/D = 0.2, while they do for B/D = 1.In particular, for B/D values above 0.6, the wake starts to grow in the downstream direction, i.e., detached eddies are pushed away from the rear face, resulting in increasing base pressure levels and reducing drag values.
The hybrid grids used in the computations are shown in Figure 1.They are composed of about 23,000-32,000 cells (more details can be found in [27]).In all cases, the time step size was chosen in order to discretize each period of vortex shedding with approximately 500 time steps.The convergence study showed that this temporal discretization is small enough, even conservative, to resolve the dominant unsteady flow phenomena [27].1.It can be noticed that the 2D results are satisfactory for the square and the rectangular 5:1 cylinder, although the Strouhal number in the latter case is slightly underestimated.By contrast, the results are far from the experiments in the case of the rectangular cylinder with B/D = 0.2, which is nearly a flat plate perpendicular to the flow.In particular, the strong overestimation of the drag coefficient is apparent.A similar result yielded for a flat plate by a 2D calculation can be found in [28] and, for an airfoil at a 90 • angle of attack, in [29].For this reason, a 3D URANS-LEA simulation (with a spanwise extension of the computational domain equal to 2•D) was carried out, and this allowed a much better estimation of the mean drag, as well as a dramatic reduction of the standard deviation of the fluctuating lift coefficient.The significant difference between 2D and 3D URANS-LEA results relies on the specific term for three-dimensional flows appearing in the additional algebraic equations for the anisotropy tensor components [12,15].The improvement of the results with a three-dimensional URANS calculation is not surprising, as it has already been underscored and discussed in [29].In particular, in this specific case, to neglect the 3D cross-flow effects in the 2D simulations allows the "sweeping effect" on the rear side of the rectangular cylinder due to the forming vortices to occur at much lower side ratios than around B/D = 0.6.This explains the overprediction of mean drag and force fluctuation amplitudes for B/D = 0.2 [27].

Reynolds Number Effects
In view of the interesting vortex-shedding mechanism [38] and of the fact that the rectangular 5:1 cylinder (B/D = 5) is often assumed as a rough schematization of a bluff bridge deck section, an international benchmark study on this geometry was launched in 2008 [39].Consequently, the numerical investigation was significantly extended for this cross-section.In particular, one of the scopes of this study was to understand the origin of the significant Reynolds number effects observed in a high-pressure wind tunnel for this sharp-edged body [37].In fact, as shown in Figure 2, if a small angle of attack is imposed on the cylinder, a non-monotonous increase of the mean lift coefficient of about 70% is observed when the Reynolds number is varied by about two orders of magnitude.A more refined 2D mesh as compared to that used in the previously discussed analyses was adopted in this case.The computational grid, composed of about 94,000 cells, is reported in Figure 3 (see [40] for more details).The 2D URANS equations were solved in conjunction with the LEA turbulence model for the case of an angle of attack of +4 • , and the results are reported in Figure 2 in terms of the mean lift and moment coefficients.It is worth noting that the mean lift coefficient is well predicted at a low Reynolds number (Re = 5200, based on the cross-flow dimension D) and that the numerical approach is able to reproduce the increase of the mean lift with the Reynolds number, although much less pronounced than in the experiments.It is also very important to remark that the Reynolds number dependance of the results was definitely not captured by the SAE, Wilcox k − ω and Menter SST models [40].Figure 4 clarifies that an increase of the Reynolds number (at least in the range of 5200 to 370,000) produces a significant reduction in the length of the recirculation bubble on the lower side of the rectangle.According to the numerical simulations, this length passes from 3.115 • D for Re = 5200 to 2.735 • D for Re = 370,000, but the variation is presumably larger in reality.As shown in Figure 5, this is due to the increased turbulent kinetic energy, in particular in the lower side recirculation bubble, which fosters the earlier shear layer reattachment.The simulations also showed that the vortex-shedding process is much more regular at a high Reynolds number [40], and this feature is clearly confirmed by the experiments.

Forced Vibration Analysis
For the rectangular 5:1 cylinder, the behavior of the body harmonically oscillating in the heaving and pitching degrees of freedom was also studied.In this case, the 2D URANS equations were solved in conjunction with the one-equation SAE eddy viscosity model [6] and the mesh shown in Figure 1 (about 32,000 cells) [41].
According to the model proposed by Scanlan [42], the lift and the moment coefficients for the oscillating body, referred to the chord B, can be written as follows: where h and α are the heaving displacement and the pitching rotation, the coefficients H * i and A * i are the flutter derivatives, K h = Bω h /U = 2π/U Rh and K α = Bω α /U = 2π/U Rα are the reduced frequencies of oscillation, ω h and ω α , respectively, the circular frequencies of heaving and pitching motion, U the undisturbed flow speed, U Rh and U Rα the reduced wind speeds.In the previous equations, the dot denotes the derivative with respect to time t.
The body was forced to oscillate around a zero-degree angle-of-attack mean position with an amplitude of oscillation of 2 • (to be interpreted for the heaving motion as an apparent angle of attack ḣ/U ) at a Reynolds number of 40,000.The center of rotation was always the geometric center of the section, coinciding also with the reduction point to calculate the torque acting on the body.Three reduced flow speeds were considered, namely 2.5, 7.5 and 15.0.
The period of oscillation of the rectangular section was discretized with a number of time steps ranging from 500 to 3000, depending on the actual reduced frequency of oscillation.In this way, the expected period of vortex shedding (referred to the stationary body) was always resolved with about 370 time steps.
In Figure 6, some numerical results are compared to the experimental data by Matsumoto [43] and Righi [44].Experimental flutter derivatives were measured with a one-DoF forced vibration system and Reynolds numbers up to 40,000 in the first case, and with a two DoF free vibration technique and Reynolds numbers up to 23,000 in the second case.The two experimental results agree only for the so-called "direct" flutter derivatives H * 1 , A * 2 and A * 3 .The functions H * 4 and A * 4 are difficult to be measured, but their influence on self-excitation is small, and they can be often neglected.The numerical results are in good agreement mainly with those reported in [43], for low and medium reduced velocities.Conversely, larger discrepancies are evident for the highest reduced flow speed, and this could be due to the strong interference with vortex shedding, which was negligible at lower reduced velocities, as clearly shown in the lift coefficient amplitude plots of Figure 7.However, the lift force component in quadrature with the heaving motion (H * 1 ) is slightly overestimated also for U Rh = 7.5.It is particularly worth remarking that, in agreement with the experiments, the sign of the flutter derivative A * 2 changes from negative to positive for a reduced velocity around U Rα = 5, thus meaning strong proclivity to one-DoF torsional galloping.

Three-Dimensional DES Simulations
For the rectangular 5:1 cylinder, several three-dimensional calculations were also performed.3D URANS-LEA equations were solved mainly as a reference and to partly account for the effect of the non-perfect correlation of pressure fluctuations on the surface of the cylinder.Moreover, 3D DES-SA equations allowed exploring the potential in terms of accuracy of the results of such an advanced strategy of turbulence modeling and to better investigate detailed features of the flow.The expected improvement of the solutions is paid with a significantly higher cost of the computations, which is not due to the equations themselves, but to the required spatial and temporal resolution (although, herein, the previously discussed URANS grids and time steps were overall fine enough also for DES) and, above all, to the need for three-dimensional meshes with large enough spanwise extension and very long running times to obtain statistically-converged results.A test case with a zero-degree angle of attack and a Reynolds number of 26,400 was considered.Extensive analyses of the results can be found in [45,46].
Figure 8 depicts a sketch of the geometry of the problem, while Figure 9 shows the 3D hybrid grid (1,703,585 nodes and 2,957,440 cells) obtained by extruding a 2D grid for a one-chord length in the spanwise z-direction (L1), using 65 equally-spaced nodes to discretize the resulting edges (∆z/D = 0.078), so that the mesh is perfectly structured in the spanwise direction.In addition, 3,380,961 nodes and 5,914,880 cells characterize the grid with double spanwise extension (L2), due to the same spacing of nodes in the z-direction.The mesh was designed according to the principles discussed in [18].The height of the first structured layer was chosen in order to have maximum wall-unit values in the normal direction y + ∼ 1 (y 1 /D = 2.5 × 10 −4 ).The structured part of the mesh consists of 28 prismatic layers with a stretching factor in the wall normal direction ∆y j+1 /∆y j = 1.23, for a total height of 0.353 • D. In the streamwise and spanwise directions, the spacing is such that the maximum x + ∼ z + ∼ 300, with refinement in the neighborhood of the cylinder sharp edges.The remaining unstructured part of the mesh consists of triangular-based prisms, and it is significantly refined in the near-wake region, where the grid is approximately isotropic (as it is in the external layers of the structured part of the mesh), since DES works in the LES mode there [18].The nondimensional time step size is ∆t * = 0.017, where t * = tU ∞ /D is the number of traveled lengths D by a fluid particle in a time unit (U ∞ = 34.0m/s is the free-stream velocity), in order to discretize the expected period of vortex shedding with more than 500 time steps.
As already mentioned in Section 2.2, DES computations with both scalar (SD) and matrix (MD) artificial dissipation were carried out.The main parameter to control the level of artificial viscosity is the fourth-order artificial dissipation coefficient (k 4 ), and the sensitivity of the solution to the chosen value was investigated.Originally, Jameson et al. [2] and Mavriplis et al. [3] recommended k 4 = 1/256.By contrast, in later papers, a value between 1/64 and 1/32 was suggested.It is worth noting that in their study, Camarri et al. [47] adopted a second-order upwind scheme and observed that the flow solution was sensitive to the amount of artificial viscosity in case a fourth-order dissipation was employed, while the results were more stable with a sixth-order dissipation.A comparison between the solutions obtained by discretizing the skew-symmetric (skew) and the usual divergence form (div) of the inviscid fluxes was also outlined.Finally, the number of inner iterations per physical time step (NITER) was increased from 100 to 150 in order to reduce the residual of the continuity equation (the drop of the residual within each time step was respectively by a factor of 4 to 10 and 15 to 50 in the two cases).The high value of sub-iterations employed is due to the high stiffness of the compressible equations at a low Mach number.Table 2 reports some integral quantities (Strouhal number, mean and standard deviation values of lift and drag coefficients) obtained from the resulting flow field computed by changing several numerical parameters.The data of the reference experimental test case [37] are also reported for comparison.A collection of additional experimental results available in the literature can be found in [45].A large number of nondimensional time units was necessary to obtain the convergence of the considered first-and second-order statistical moments.The 3D URANS-LEA results are very close to those obtained with a two-dimensional simulation, apart from a significant reduction of the standard deviation of the drag coefficient [45].The DES computation with scalar artificial dissipation (second row of Table 2) gave results closer to the experiments and, in particular, a significantly larger value of the Strouhal number and a smaller value of C L .The doubling of the spanwise extension of the computational domain does not affect the Strouhal number and the mean drag coefficient, but it implies a significant reduction of lift and drag standard deviations.By contrast, as a consequence of the use of matrix dissipation, the Strouhal number increases while the drag coefficient decreases.Furthermore, the standard deviation of both lift and drag are significantly reduced.The effect of increasing the number of inner iterations is small, but not null, as a slight reduction of the standard deviation values of drag and lift can be appreciated.The use of the skew-symmetric inviscid flux operator seems to imply a smaller value of the Strouhal number and a slight increase of the standard deviation of lift and drag.It is particularly important to note the strong change in the computed quantities (with the exception of the mean drag coefficient) when the fourth-order artificial dissipation coefficient k 4 is reduced from 1/84 to 1/256.In this case, the insufficient damping of the spurious high-frequency numerical oscillations seems to cause a completely different flow solution.Moreover, if the value of the coefficient k 4 is further reduced to 1/1000, the computation becomes unstable.From the analysis of the corresponding power spectral densities (Figure 10), it is evident how the high-frequency fluctuations become more energetic if matrix dissipation is chosen and, in particular, if the k 4 constant is reduced.This effect is even more evident if local flow quantities are considered.In fact, in terms of turbulence resolution, the key role played by the amount of artificial viscosity is evident from the analysis of the vorticity magnitude isocontours reported in Figure 11.With scalar dissipation, the high wavenumber structures are suppressed.By using the matrix dissipation scheme, finer vortical structures can be resolved, and the full potential offered by the chosen computational grid and time step size seems to be better exploited.The effect of reducing the fourth-order dissipation coefficient is evident in Figure 12, which shows an isosurface of the second invariant Q of the velocity gradient tensor.
The time-averaged flow fields depicted in Figure 13 confirm that the flow solutions obtained with k 4 = 1/84 and k 4 = 1/256 are completely different.In the case with matrix dissipation and k 4 = 1/84, the shear layer mean reattachment point moves non-negligibly upstream (x/D = 1.69) with respect to the computation with scalar dissipation (x/D = 2.25).When the fourth-order dissipation is strongly reduced (k 4 = 1/256), the length of the mean recirculation bubble becomes so small (x/D = 0.54) to be definitely incompatible with the pressure measurements available for Reynolds numbers similar to the one considered (Figure 14).It is worth noting that Bruno et al. [48] obtained x/D = 2.165 with an LES simulation.This fact underlines the sensitivity of the numerical results to the amount of artificial viscosity introduced into the solution through the discretization of the inviscid fluxes.Moreover, it is apparent in Table 2 that significant differences in the mean reattachment location produce variations in the Strouhal number and, above all, in the mean drag coefficient of a few percent only, whilst much more evident is their effect on the standard deviation value of the lift coefficient.To sum up, this flow feature seems to be very sensitive to spatial discretization, turbulence modeling and numerical dissipation.2).However, the wind tunnel results discussed in [38] show that also in the experiments, the location of the shear layer mean reattachment is a critical variable, and significantly different results can easily be obtained in different laboratories.In particular, a crucial role in the shear layer reattachment is played by the location of laminar-to-turbulent transition [37] and by the oncoming turbulence content, which is able to modify the curvature and the thickness growth rate of the shear layers [52].In the present simulations, the reduction of artificial numerical dissipation through the coefficient k 4 enables an over-development of turbulence in the free shear layers at scales allowed by the mesh and the time step size, playing a role similar to oncoming free-stream turbulence.
The important effects of the spanwise extension of the computational domain were also discussed in [45] and [53].Figures 15 and 16 show that, by increasing the spanwise period L, the correlation coefficient of pressures ρ p tends to decrease, as the influence of periodic boundary conditions is reduced and three-dimensional vortical structures are more free to develop.This reduction of pressure correlation entails the decrease of the standard deviation of lift and drag coefficients in the case of domain L2 as compared to domain L1, as previously highlighted in Table 2. Nevertheless, it can also be noted in Figures 15 and 16 that, if matrix dissipation is employed, finer turbulent structures are resolved and, consequently, the spanwise correlation of pressures is significantly reduced.This effect is responsible for the smaller standard deviation of lift and drag coefficients with respect to the computation with scalar dissipation (Table 2).

Bridge Section
Encouraged by the satisfactory results obtained for rectangular cylinders, the 2D URANS-LEA approach was also applied to the case of a realistic bridge section.The chosen case study is the geometry of a single-box girder deck section model with lateral cantilevers tested in the CRIACIV (Inter-University Research Center on Building Aerodynamics and Wind Engineering) Boundary Layer Wind Tunnel, Prato, Italy.Details concerning the model and the experimental setup can be found in [55].Figure 17 shows a schematic view of the cross-section considered.
Figure 17.Schematic of the bridge section geometry.The reference system for generalized displacements and forces is also highlighted.This common geometry for bridge decks was inspired by the Sunshine Skyway Bridge (SSB), Tampa, FL, USA, a cable-stayed bridge with a main span of 364 m, that was tested at a scale of 1:80 in the wind tunnel of the University of Western Ontario, London, ON, Canada [56].This section model was very similar to the one tested at CRIACIV, but had a slightly different chord-to-thickness ratio (B/D = 6.43 for the section tested at CRIACIV and B/D = 6.80 for the Sunshine Skyway Bridge) and presented a certain chamfer at the connection between the lateral cantilevers and the inclined walls of the box.In particular, the degree of sharpness of box lower corners was different, since in the Sunshine Skyway Bridge's case, these were made sharp through acrylic blocks.By contrast, in the CRIACIV section model, these corners were slightly rounded due to the manufacturing of the aluminum foil.Given the significant difference between the two sets of experimental results [55,56], the CFD investigation was mainly devoted to explain this discrepancy and to underscore the key role played by the degree of sharpness of the lower corners for both the stationary and harmonically-oscillating section.Numerical simulations were carried out for a Reynolds number of 600,000, based on the section width B.
The computational hybrid grids used are shown in Figure 18

Stationary Case
For the simulations of flow past the stationary bridge section, the time step size was chosen in order to resolve the expected period of vortex shedding with more than 360 time steps, in the case of sharp edges, and 220 time steps, in the case of smooth edges.A grid-convergence and time step sensitivity study was discussed in [57], showing the stability of the solution with respect to the grid refinement and temporal discretization.
Figures 19 and 20 show that the computed Strouhal numbers and the aerodynamic force and pressure coefficients for the angles of attack considered (α = −5 • , 0 • , +5 • ) are much closer to the experimental data for the SSB [56] than those for the CRIACIV section [55], despite the latter being actually considered in the numerical simulations.As previously mentioned, the most relevant difference between the two geometries is the degree of sharpness of the lower corners, which are perfectly sharp in the SSB case, while they are not in the CRIACIV section.In order to investigate the effect of this feature, a small radius of curvature R was introduced in the computed geometry of the bridge profile.Since it is fairly difficult to estimate the actual value of R in the wind tunnel section model, three different cases were considered: R/B = 0.025, 0.05 and 0.075.The experimental results are from [55] for the CRIACIV section (Re = 576, 000) and from [56] for the Sunshine Skyway Bridge (SSB) section (Re = 295, 000).For the Strouhal number, the 68.3% confidence interval (µ St ± σ St ) is indicated (75, 000 ≤ Re ≤ 810, 000).The moment coefficient is referred to the center of pitching rotation P (Figure 18).The analysis of the integral quantities shown in Figure 19 highlights the importance of the degree of sharpness of the bridge deck lower corners for the resulting flow field.The introduction of a radius of curvature R/B = 0.05 makes the numerical results approach the experimental data for the CRIACIV section for α = 0 • and −5 • .In fact, the drag decreases, and the Strouhal number takes higher values, while the mean lift coefficient at α = 0 • becomes negative, although with a smaller absolute value than in the experiments.The lift at α = −5 • and the moment coefficient at both α = 0 • and −5 • match the wind tunnel data well.By contrast, at α = +5 • , the corner roundness does not produce any significant difference as compared to the sharp-edged configuration.For the larger corner radius, R/B = 0.075, the results are substantially the same as for R/B = 0.05.Conversely, for R/B = 0.025, they are half way between those for R/B = 0 and R/B = 0.05.
Figure 21 clarifies the reason behind the large influence of corner roundness: flow separation occurs at the upstream lower corner in the sharp-edge configuration, while the flow is attached until the downstream corner in the case of the model with rounded edges.Therefore, the width of the wake gets smaller and the drag reduces, while the frequency of vortex shedding increases (higher Strouhal number).This effect is also apparent in Figure 22, which reports the mean-flow streamlines.This figure demonstrates that the difference in the lower corner sharpness also affects the leading-edge flow field.In fact, the separation bubble over the upper bridge deck surface is much smaller in the case of rounded corners.The larger separation zones on the upper side and below the rear portion of the bridge deck section also entail a higher positive (nose-up) moment acting on the sharp-cornered section.

Flutter Derivatives
Forced harmonic vibration numerical tests were carried out for the bridge section in pure heaving and pitching with the aim to determine the flutter derivatives and, therefore, the self-excited forces (Equations ( 1) to ( 2)).Similar analyses for slightly different bridge deck geometries can be found, for instance, in [58][59][60][61].The simulations were performed for the configuration with both sharp (R/B = 0) and rounded edges (R/B = 0.05), and the only experimental results available were those for the CRIACIV section [55,62], which, as previously mentioned, presents lower edges with a certain degree of roundness.
Figure 18 shows the origin of the moving coordinate system, located in the point P of coordinates (0.5, −0.067) with respect to the inertial system xz.The point P was also the center of pitching rotation of the section model employed in the reference experimental tests.The numerical simulations were carried out for three reduced velocities for both heaving and pitching harmonic motion: U Rh = U Rα = 3, 6.5 and 10.The Reynolds number was always Re = 600,000, and the mean angle of attack was kept null.A small amplitude of motion was chosen coherently with the reference experiments, namely α 0 = 2 • in pitch and h 0 /B = 0.02 in heave.The effects of changing the motion amplitude, the angle of attack and the Reynolds number were studied in [63,64].
The period of oscillation of the bridge section was discretized with a number of time steps ranging from 750 to 3700, depending on the geometrical configuration and the considered reduced frequency of oscillation [63,64].In this way, the expected period of vortex shedding (referred to the stationary body) was always resolved with more than 220 time steps.
The complete set of computed flutter derivatives is shown in Figure 23.As expected, the flutter derivatives for the smooth-edged case are closer to the experiments, confirming the results obtained for the stationary profile.In particular, the coefficient A * 2 changes sign for U Rα around four to five for the sharp-edge configuration, while it does not for the smooth-edged configuration, remaining always negative.This difference is extremely important, as it suggests that the section with sharp edges is prone to torsional galloping, whereas only two degrees of freedom flutter at higher reduced wind speed can arise if the box girder lower edges are sufficiently rounded.In this case, the A * 2 coefficient takes positive values despite the fact that the slope of the torque coefficient around a null angle of attack is positive [57].A similar behavior was also observed for the rectangular 8:1 and 10:1 cylinders [49].The tendency to low reduced velocity torsional instability in the case of sharp lower edges is confirmed by the dynamic tests performed on a section model of the Sunshine Skyway Bridge [65], although in that case, the detrimental role on the stability threshold played by lateral and central New Jersey barriers was highlighted.[55,62] in the form of 95% confidence intervals.The usual sign convention of bridge aerodynamics was adopted here, i.e., that of Figure 17 but with the lift force positive downward.
In addition, apart from the coefficients H * 4 and A * 4 , which are difficult to be measured and not very relevant for practical applications, the significant discrepancy between experimental and numerical results for the flutter derivative H * 1 is noted.In this case, the numerical result for smooth edges is closer to that predicted by the thin airfoil theory than to experiments.In the case of sharp edges, the H * 1 derivative also shows negative values, but significantly smaller in magnitude than the wind tunnel data, which means that less aerodynamic damping is introduced in a pure heaving motion.

Conclusions
The main conclusions of the present work can be summarized as follows: • The 2D URANS approach combined with the EARSM-LEA turbulence model gives results of reasonable accuracy for the rectangular cylinders and the bridge section considered.In particular, complex phenomena, such as the Reynolds number effects observed in the wind tunnel for the sharp-edged rectangular 5:1 cylinder (B/D = 5), are qualitatively captured, although they result in being less pronounced than in the experiments.• In the particular case of a rectangular 1:5 cylinder (B/D = 0.2), with the long side perpendicular to the flow, 3D simulations are necessary to obtain acceptable results, even with the URANS equations.
• Forced vibration simulations were carried out for the rectangular 5:1 cylinder with 2D URANS equations in combination with the eddy viscosity model of Spalart and Allmaras (SAE).The results are in overall good agreement with experiments, but significant improvements can be expected by adopting an EARSM turbulence model.In particular, the numerical calculations correctly predict the expected threshold of torsional galloping instability.• The 3D DES-SA approach for the rectangular 5:1 cylinder delivers more accurate results than both the 2D and 3D URANS equations.Nevertheless, the study highlights that the results are sensitive to the amount of artificial dissipation introduced to stabilize the central difference discretization of the convective term in the governing equation.• A large spanwise extension of the computational domain is required in a 3D DES simulation to allow the natural loss of the correlation of the fluctuating pressure field.• 2D URANS-LEA numerical simulations of flow past a realistic bridge section give results in satisfactory agreement with experiments and helped to understand the discrepancy between two sets of experimental data available in the literature for two very similar geometries.• The numerical results underscore the key role played by the degree of sharpness of the bridge section lower corners with respect to both static and aeroelastic behavior.In particular, the section is prone to low speed torsional galloping in the case of sharp edges, while higher speed coupled flutter is expected if the lower corners present a small radius of curvature.

Figure 1 .
Figure 1.Near-body views of the computational grids used for the simulations of flow past rectangular cylinders of various side ratios.

Figure 2 .
Figure 2. Mean lift (a) and moment (b) coefficients for the rectangular 5:1 cylinder at a 4 • angle of attack vs. the Reynolds number: comparison of current numerical results and experimental data from [37].

Figure 3 .
Figure 3. Near body (a) and wake refinement (b) close-ups of the grid used for the 2D rectangular 5:1 section.

Figure 7 .
Figure 7. Numerically-simulated lift coefficient vs. apparent angle of attack in the case of heaving harmonic motion for a reduced flow speed equal to 7.5 (a) and 15.0 (b).

2 LFigure 10 .
Figure 10.Comparison of lift power spectral densities obtained with three different DES computations for the rectangular 5:1 cylinder (see Table2).

Figure 15 .Figure 16 .
Figure 15.Longitudinal lines considered for the calculation of the correlation of pressures for the rectangular 5:1 cylinder: a spanwise position close to the leading edge (a) and one close to the trailing edge (b) were chosen.

Figure 18 .
Figure 18.Computational grids for the bridge section: computational domain (a); wake refinement (b); near body mesh (c); close-up of the upstream lower corner for the sharp-edge configuration (d) and for a smooth-edge configuration (e).

Figure 20 .
Figure 20.Mean (a) and standard deviation (b) values of the pressure coefficients for the configuration with perfectly sharp edges.The experimental results for the SSB section model refer to a null angle of attack and to the bare-deck configuration [56].

Figure 21 .Figure 22 .
Figure 21.Streamwise velocity flow field for the bridge section with sharp (a) (R/B = 0) and with rounded lower edges (b) (R/B = 0.05).The angle of attack is null, and the undisturbed flow speed is U ∞ = 34 m/s.The indicated flow velocities range between −20 and 50 m/s.

3 Figure 23 .
Figure 23.Computed and measured flutter derivatives for the bridge section.R/B = 0 and 0.05 denote, respectively, sharp and smooth lower edges.Experimental data from[55,62] in the form of 95% confidence intervals.The usual sign convention of bridge aerodynamics was adopted here, i.e., that of Figure17but with the lift force positive downward.

Table 2 .
Strouhal number and mean (C D , C L ) and standard deviation (C D , C L ) values of drag and lift coefficients for several 3D URANS and DES computations of flow past the rectangular 5:1 cylinder.t * IN T denotes the nondimensional time units considered.NITER, number of inner iterations per physical time step; SD, scalar artificial dissipation; div, divergence; SA, Spalart-Allmaras; DDES, delayed DES; MD, matrix artificial dissipation.