Abstract
The incompressible high Reynolds number flow around the rectangular cylinder with aspect ratio 5 : 1 has been extensively studied in the recent literature and became a standard benchmark in the field of bluff bodies aerodynamics. The majority of the proposed contributions focus on the simulation of the flow when a smooth inlet condition is adopted. Nevertheless, even when nominally smooth conditions are reproduced in wind tunnel tests, a low turbulence intensity is present together with environmental disturbances and model imperfections. Additionally, many turbulence models are known to be excessively dissipative in laminar-to-turbulent transition zones, generally leading to overestimation of the reattachment length. In this paper, Large Eddy Simulations are performed on a 5 : 1 rectangular cylinder at non-null-attack angle aiming at studying the sensitivity of such flow to a low level of incoming disturbances and compare the performance of standard Smagorinsky-Lilly and Kinetic Energy Transport turbulence models.
1. Introduction
Thanks to the increase in the computer power, the research in the field of Wind Engineering is more and more making use of Computational Fluid Dynamics techniques in order to assess the effects of wind loading on structures.
In fact, a number of flows usually encountered around structures relevant to Civil Engineering applications, like bridge decks and buildings, are characterized by the presence of shear layers and separation bubbles which render their prediction by means of CFD an extremely challenging task.
Due to their special geometric simplicity, prismatic shapes have been deeply investigated in the literature and have been often adopted as a prototype of fully detached and reattached flows. A number of experimental studies identified the main mechanisms involved in the definition of the flow field around such simple, nominally two-dimensional, shapes and highlighted the fundamental role played by the stability conditions of the shear layers detached from the leading edges, especially for shapes whose aspect ratio is higher than three, which approximately corresponds to the threshold separating fully detached and reattached flows [1–4]. In particular, Nakamura and Ohya focused on the effects of incoming turbulence, showing that the reattachment point migrates upstream with increasing turbulence levels [1]. Subsequently, it was clarified that Kelvin-Helmholtz instabilities appear in the detached shear layers due to incoming disturbances, destabilizing them and, thus, causing the reattachment point upstream migration. This result was also highlighted by Hillier and Cherry [5] and Kiya and Sasaki [6] that investigated the effects of free-stream turbulence on the topology of separation bubbles. Their experiments showed that the incoming flow turbulence deeply influences pressure distributions in terms of both mean and standard deviation. Due to these considerations, it appears that the incoming flow turbulence should be accurately taken into account when dealing with wind loading problems, for both experimental and numerical investigations.
With regard to numerical studies, several research works have been proposed aiming at reproducing the flow field around bluff bodies. Among them, Yu and Kareem performed Large Eddies Simulations (LES) around rectangular cylinders with varying aspect ratio at null-attack angle in a two- and three-dimensional framework [7]. Sohankar et al. [8] proposed simulations of rectangular cylinders at non-null-attack angle focusing on very low Reynolds numbers, for which the flow is expected to be mainly laminar. Such simulations showed good results in terms of bulk parameters and highlighted that when using the standard deviation of the lift coefficient as an indicator, the simulation appears to be very sensitive to the adopted numerical setups. Shimada proposed a detailed validation of a two-layer model with a modification on the turbulent kinetic energy production term by analyzing the flow field around a wide selection of rectangular cylinders at null-attack angle [9] while Noda and Nakayama proposed LES including the effect of incoming turbulence (a turbulence intensity equal to 5% was considered) showing that the standard Smagosinsky model is suitable for reproducing the effects induced by incoming turbulence for engineering purposes [10]. Moreover, a wide number of experimental and computational studies have been also proposed in order to study the aeroelastic behavior of such bodies in smooth [11] and turbulent conditions [12]. In particular, Daniels et al. [13] performed LES on an elastically mounted rectangular cylinder in smooth and turbulent inflow conditions, showing that increasing the turbulence intensity can diminish significantly the amplitude of oscillations.
In all, it appears that while all computational models are able to approximately reproduce some characteristics of the flow field, their accuracy strongly depends on the adopted numerical setups and on the quantities taken into consideration. Indeed, if on one hand the mean flow characteristics are generally reproduced by numerical simulations in a satisfactory way, on the other hand the standard deviation of velocity and pressure fields is often predicted with much less accuracy [14]. It also emerged that, due to the large number of parameters which might affect the flow topology, a detailed mapping of the prisms aerodynamic behavior has not been yet achieved.
Recently, the BARC project focused on the simulation of the turbulent flow around a 5 : 1 rectangular cylinder and became a standard benchmark for comparison between experimental data and numerical simulations [15]. The overview of the first four years of activity of the BARC project [16] highlighted that even experimental results extracted by various research groups show a remarkable variability whose causes are currently not completely clear. When numerical simulations are considered, Bruno et al. [16] observed that results are even more dispersed showing great sensitivity to the adopted numerical schemes, turbulence model, and mesh size while their relative importance is difficult to be evaluated at the current stage. In this context, uncertainty studies can play a fundamental role. In particular, Witteveen et al. [17] showed that the flow field around the 5 : 1 rectangular cylinder is very sensible to setup parameters as small variations of the angle of attack, the incoming turbulence intensity, and the incoming turbulence length scales.
It is also noticed that while disturbances are always present in experimental conditions, when numerical simulations are performed, a fictitious environment is produced where imperfections are introduced only by the differential problem discretization and by the numerical solution itself. Beside the absence of incoming disturbances, many turbulence models are known to be excessively dissipative in laminar-to-turbulent transition zones [18]. Thus, numerical simulations can be considered as a limit case, expected to produce excessively stable shear layers and, thus, generally overestimating the reattachment length.
In this paper, the performances of the standard Smagorinsky-Lilly [19, 20] model and the Kinetik Energy Transport [21] model are tested when the incoming flow is characterized by a very low turbulence intensity. To this purpose, the flow around a fixed 5 : 1 rectangular cylinder is considered at 4° attack angle. In fact, the mismatch between experimental and numerical results has been found to increase with the attack angle when perfectly smooth inlet conditions are adopted together with a standard SM turbulence model [22].
The paper is organized as follows. In Section 2 the experimental setup adopted in the wind tunnel tests is briefly described, while in Section 3 the numerical features of the computational model are discussed. Numerical results obtained from simulations are presented in Section 4 and compared with the experimental data. Finally, in Section 5 some conclusions are drawn.
2. Experimental Setup
In this section, the setup used to extract the experimental results used for comparison with numerical simulations is described. The experimental study was performed in the open-circuit boundary layer wind tunnel of CRIACIV laboratory, located in Prato, Italy [23].
An aluminium model characterized by width, , equal to 300 mm, depth, , equal to 60 mm, and length, , equal to 2380 mm was employed in the experimental tests and equipped with pressure taps monitored by two 32-channel PSI miniaturized piezoelectric scanners at a sampling rate of 500 Hz [24]. The wind tunnel test section is 2.42 m large and 1.60 m high, so that the resulting blockage ratio was 3.75%. The Reynolds number based on was varied from to , while incoming turbulence intensity was varied from 0.7% to 13.6%.
In the present study, the experimental configuration with nose-up angle equal to 4° at Reynolds number based on equal to has been used as a reference for the numerical simulations.
3. Computational Model
In this section, the characteristics of the computational model are reported and compared with other studies, together with the description of the adopted boundary conditions and numerical setups.
According to the aim of this study, the lowest turbulence intensity recorded in the available wind tunnel tests, representing the experimental approximation of the smooth flow condition, is adopted, so that the incoming turbulence intensity is set to be equal to 0.7%. The inflow turbulent field is synthetically generated by using the Modified Discretizing and Synthesizing Random Flow Generator method proposed by Castro et al. [25, 26], which guarantees the divergence-free condition and is able to satisfy prescribed turbulence spectrum and spatial correlations. It is worth noting that if the fluctuation field is introduced at the inlet boundary, large pressure fluctuations are generated due to the violation of the Navier-Stokes equations. In order to avoid this phenomenon, fluctuations are introduced in the computational domain by modifying the velocity-pressure coupling algorithm following the procedure reported by Kim et al. [27].
The computational domain size and the grid resolution have been defined according to the guidelines provided by Bruno et al. [16] and the BARC main setup [15]. As showed in Figure 1, the domain dimensions are such that and , while . The resultant blockage ratio is 0.67%, and therefore blockage ratio effects can be neglected. In order to avoid boundary effects on the solution, the distance between the prism leading edge and the inlet boundary is set to be equal to .
As showed in Figure 2(a), a structured mesh is adopted in the boundary layer close to the wall, with along-wind dimension and crosswind dimension . The maximum recorded during all simulations is found to be equal to , while the mean one equals 2.02. Outside the boundary layer, the mesh is unstructured and quad dominated and its size is slowly coarsened up to approximately in the wake (see Figure 2(b)), where the cells aspect ratio is approximately one. In order to limit the numerical dissipation and to propagate the fluctuation field minimizing the energy loss, the mesh sizing in front of the rectangular cylinder has dimensions and is maintained almost constant and structured until the inlet boundary is reached.
(a)
(b)
In direction, the cell dimension is . The final resultant mesh is composed of about finite volumes. A comparison between the domain size and the mesh resolution adopted in the present study and similar ones is reported in Tables 1 and 2.
With respect to turbulence modeling, two subgrid scales models have been considered, namely, the standard Smagorisnky-Lilly (SM) model and the Kinetic Energy Transport (KET) model. The fluid flow governing equations are reported inwhere the overbar denotes the spatially filtered quantities and represents the velocity in the th direction, while is the pressure, is the resolved viscous stress tensor, and is the subgrid stress tensor. The resolved viscous stress tensor can be written aswhere represents the resolved strain rate tensor, whose expression is reported below:By assuming Boussinesq’s hypothesis, the subgrid stress tensor can be written aswhere is the subgrid kinetic energy:and is the turbulent eddy viscosity. When the SM model is adopted, is expressed aswhere is the Smagorinsky constant, set to be equal to 0.12, while is the local grid spacing. Conversely, if the KET model is adopted, the transport equation for the subgrid kinetic energy is considered in addition to equations reported in (1):and the turbulent eddy viscosity is calculated aswhere and are model constants set to be equal, respectively, to 1.05 and 0.094.
With respect to the numerical scheme, a centered second-order accurate scheme is adopted to discretize the spatial derivatives, exception made for the nonlinear convective term. For this term, the Linear Upwind Stabilized Transport scheme is used, which proved to be particularly successful for LES in complex geometries [28].
Time integration is performed by using the the two-step second-order Backward Differentiation Formulae, in accordance with Bruno et al. [29]. The solution at each time step is obtained by means of the well known Pressure Implicit with Splitting of Operator algorithm. The adopted nondimensional time step (based on ) is , leading to a maximum Courant number in all the simulations equal to .
Symmetry boundary conditions are imposed at the top and bottom surfaces, while periodic conditions are adopted on the faces normal to the spanwise direction (see Figure 1). At the outlet boundary, zero pressure is imposed, while, at the inlet, the null normal gradient of pressure is prescribed.
The generated fluctuation field is introduced in the computational domain in a plane shifted with respect to the inlet boundary of approximately , while the mean velocity is prescribed at the inlet. The power spectral density of the velocity components along , , and directions, indicated, respectively, as , , and , is reported in Figure 3 for two different points located near the inlet boundary at and just upstream the rectangular cylinder at . As it can be seen, the high frequency content of the spectra appears to be slightly damped proceeding from the inlet towards the rectangular cylinder. With respect to the along-wind turbulent length scale , it is found to be equal to for both of the aforementioned locations.
The prism has been equipped with 2500 pressure monitors and data are acquired at each time step, leading to 200 samples for one nondimensional time unit. All simulations have been run by using the open source finite volume software on 160 CPUs at CINECA on the Galileo cluster (516 nodes, 2 eight-core Intel Haswell 2.40 GHz processors with 128 GB RAM per node).
4. Numerical Results
In this section, numerical results obtained for each turbulence model with smooth and turbulent inflow condition are reported and systematically compared with experimental data. In particular, the flow bulk parameters obtained from each simulation together with a discussion of the resulting flow topology are reported in Section 4.1. The pressure coefficient statistics on the central section are showed in Section 4.2, while Section 4.3 focuses on the spanwise-averaged pressure distributions. Then, in Section 4.4, correlations in the spanwise direction are discussed.
Numerical results have been obtained by considering a simulation time of , being the nondimensional time unit. The postprocessing of data has been carried out considering only the last 700 nondimensional time units in order to avoid flow initialization effects [30].
In order to check the convergence of the recorded statistics, the time-history of the lift coefficient has been subdivided into ten segments and first- and second-order statistics were extracted by incrementally extending the part of the signal considered in the postprocessing. Despite the relatively long simulation time (longer than the minimum requirements [30]), in the worst case, a plateau has been reached for second-order statistics only when 90% of the total time-history has been used.
4.1. Flow Topology and Bulk Parameters
In order to have a qualitative view of the turbulent structures obtained from the numerical analyses, the flow topology in terms of isosurfaces of the invariant [31] coloured with pressure is reported in Figure 4 for the two investigated turbulence models, in smooth and turbulent inflow conditions. The instantaneous isosurfaces are plotted to correspond with maximum lift. As it can be seen, when the inflow turbulence is considered, the topology showed by the SM model appears to be quite in accordance with that obtained by means of the KET model. Furthermore, both of them comply well with the topology showed by the KET model in smooth inflow conditions, while in this case the topology obtained by using the SM model shows a narrower wake.
The different behavior of the two turbulence models is also reflected in Table 3, which shows that the reattachment point at the bottom surface moves from to when turbulence is introduced and the SM model considered, while no significant differences between turbulent and smooth inlet are observed if the KET model is analyzed.
The effect of incoming turbulence can be also appreciated by observing the streamlines of the time-averaged flow fields reported in Figure 5. In particular, the SM model appears to be significantly affected by the presence of incoming turbulence while the KET model appears to be rather stable. It is also noticed that all considered models, apart from the SM in smooth flow, predict the creation of a small elongated vortex at the leading edge, in correspondence with the top side, whose core is located at approximately .
The statistics of the flow bulk parameters are reported in Table 4 (results are made nondimensional with respect to ) and referred to the global reference system. In particular, and are the drag and the lift coefficients, respectively, while their root mean square values are indicated as and , respectively. Again, the results obtained by the two turbulence models are in good agreement when a low level of incoming turbulence is considered while remarkable differences are observed in perfectly smooth flow. Furthermore, it should be noticed that and , together with the Strouhal number (St), obtained by taking into account the incoming turbulence agree quite well with experimental measurements. In order to provide a clearer picture of the obtained results, in the following, pressure field statistics distributions are analyzed.
4.2. Central Section Statistics
In this section, the pressure coefficient statistics on the prism central section are analyzed. Data are presented by adopting the curvilinear abscissa as reported in Figure 1. The upwind face of the prism is identified by , and the along-wind surfaces are identified by while the downwind face is identified by .
Figure 6 shows the distribution of the pressure statistics obtained by means of the present simulations and comparison with experimental results. As it can be noticed, if the top surface is considered, the time-averaged pressure coefficient, , distribution predicted by the KET model is almost in perfect accordance with the experimental data for both smooth and turbulent inflow conditions. Contrarily, the SM model approaches experimental results only when turbulent inlet conditions are adopted. Observing on the same surface, the KET model always complies with experimental measurements, with the peak position and the experimental trends being reproduced with good accuracy. Generally, the SM model appears to be unable to correctly reproduce the shape of the pressure recovery not catching the change of steepness that occurs at about , in correspondence with the separation between the vortex shedding and the vortex coalescing zones as individuated in [30].
When distribution on the bottom surface is analyzed, the SM model appears to be very sensitive to small turbulence intensities in the incoming flow. In fact, the distribution obtained in smooth flow conditions clearly changes when the inflow turbulence is considered: the reattachment point migrates upstream and the distribution shifts, approaching the experimental data. The KET model, instead, shows again a lower sensitivity to the incoming turbulence, results obtained in both configurations being approximately the same. The SM model is very accurate when turbulence is taken into account and shows better performances if compared with the KET model, which slightly overestimates the reattachment length and underestimates suctions at the leeward edge. The improvement of the SM model with turbulent inlet observed for mean values can be noticed also by considering distribution on the bottom surface. In this case the model correctly predicts the peak position and the trend observed experimentally, even if the peak value is clearly overestimated. The KET model correctly predicts the peak position and the trend in particular at the trailing edge, where an increment in rms is experimentally observed.
The improvement of performances showed by the SM model when the incoming turbulence is taken into account can be due to the fact that the SM model is too dissipative in laminar and transitional regions [18]. Indeed, the SM model predicts a nonvanishing eddy viscosity when the flow is laminar and this results in more stable shear layers when the incoming flow is smooth. Conversely, the KET model is able to adjust the eddy viscosity basing on the subgrid kinetic energy, leading to an overall less dissipative behavior. When the incoming flow turbulence is considered, the SM model performs better. In fact, the laminar-to-turbulent transition is driven by disturbances introduced by the explicitly simulated incoming flow turbulence. This different behavior in perfectly smooth and low turbulence inflow conditions is not observed when the KET model is considered. Indeed, the less dissipative behavior of the KET model with respect to the SM model allows correctly predicting the destabilization of the shear layers without triggering it explicitly by introducing disturbances.
4.3. Spanwise-Averaged Quantities
distributions of , denoted as , are reported in Figure 7 for each case. Data concerning the central section statistics (see Figure 6) are repeated in this figure in order to better highlight differences with respect to the corresponding values. Focusing on the top surface, the SM model shows that the gap between and increases proceeding from the leading edge to the trailing one, where the peak value of is radically decreased with respect to the central alignment for both smooth and turbulent inlet condition. The flow two-dimensionality is dominant close to the leading edge, where statistics are identical to the central section ones, while the flow three-dimensional mechanisms prevail downwind. This trend is confirmed by the KET model, even if in this case the gap between and is reduced if compared with the SM model, suggesting that the predicted flow is more correlated in the spanwise direction.
Considering results of the SM model on the bottom surface, the reduction between and central section statistics is higher when the inflow turbulence is considered. In fact, while for the smooth inlet condition peak value is roughly 66% of the central section one, when turbulence is taken into account, the peak of is only 44% of the corresponding value. Therefore, the SM model shows that the inflow turbulence contributes to decreasing the flow correlation in the spanwise direction to a great extent.
In agreement with the previously reported results, the KET model does not show such a high sensitivity, differences between and in smooth and turbulent condition being almost identical. Qualitatively, by observing the rms distributions in terms of both peak amplitude and position, results obtained by using the KET model appear to be intermediate between the ones obtained in smooth and turbulent flow by using the SM model.
4.4. Correlations
This section reports correlations along the spanwise direction for three different sections, located at , and 4.25. Data acquired on probes close to the edges in the spanwise direction can be affected by the prescribed boundary conditions; therefore results are presented disregarding part of the rectangular cylinder near to the boundary, ranging from to . Figure 8 shows the correlation of , indicated as , obtained for the SM model on top and bottom surfaces for smooth and turbulent inlet, while results obtained by using the KET model are reported in Figure 9.
Interestingly, in this case the correlation appears to be higher when a small incoming turbulence level is considered. The authors conjecture that this behavior might be related to the fact that when perfectly smooth conditions are considered, the flow impinges on the trailing edge as it can be deduced from time-averaged streamlines reported in Figure 5. Contrarily, when incoming turbulence is considered, a higher level of entrapment of the separation bubble is observed, probably leading to higher along-span correlations. Regarding the bottom surface and focusing on results obtained by means of the SM model, an opposite behavior is recorded. This fact might be due to the development of spanwise vortical structures that decrease the flow correlation (see, e.g., Sasaki and Kiya [36]).
Also, in this case, results predicted by using the KET model are qualitatively intermediate between the ones obtained by using the SM model in smooth and turbulent conditions and a good stability with respect to the incoming turbulence level is observed.
5. Conclusions
In this paper, Large Eddy Simulations have been performed aiming at studying the unsteady flow field around a rectangular cylinder with aspect ratio 5 : 1 at attack angle, when a small level of incoming turbulence is taken into account. Two different turbulence models have been considered, namely, the classical Smagorinsky-Lilly model and the Kinetic Energy Transport model. Both of them have been studied by adopting perfectly smooth and turbulent inlet conditions, with turbulence intensity set to be equal to . The inflow turbulence has been synthetically generated with the MDSRFG method and introduced in the computational domain by modifying the velocity-pressure coupling algorithm. Results in terms of flow bulk parameters, time-average, and rms of the pressure coefficient distributions in correspondence with the prism central section have been compared to available experimental data. It appears that when the Smagorinsky-Lilly model is adopted, the modeling of a realistic turbulent inflow is important in order to obtain accurate results, since even small values of the incoming flow turbulence intensity can deeply affect results. The Kinetic Energy Transport model proved to be less sensitive to the low inflow turbulence and provided results qualitatively intermediate between the Smagorinsky-Lilly models adopted with perfectly smooth and turbulent inflow conditions. The different behavior of the two considered models can be due to the fact that the Smagorinsky-Lilly model is not able to make the turbulent eddy viscosity null, appearing to be too dissipative and failing in accurately predicting the laminar-to-turbulent transition. On the contrary, the Kinetic Energy Transport model can adjust the turbulent eddy viscosity depending on the subgrid kinetic energy, showing satisfactory results for both smooth and turbulent inlet conditions. Indeed, the less dissipative behavior of the Kinetic Energy Transport model with respect to the Smagorinsky-Lilly model allows the shear layer instabilities to develop even without directly triggering them with incoming disturbances. In all, when incoming turbulence is considered, a good agreement between numerical results and experimental data is observed, in particular in terms of time-averaged pressure distributions. Moreover, the differences between the flow fields obtained by the two considered models significantly reduce.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The authors are grateful to CRIACIV for providing the experimental results and CINECA for providing the HPC facilities needed to accomplish the present research work. In particular, they would like to thank Ivan Spisso for his invaluable constant support.