Passive Control of the Flow Around a Rectangular Cylinder with a Custom Rough Surface

Abstract

Motivated by existing techniques for implementing roughness on cylinders to control flow disturbances, we performed delayed detached eddy simulations (DDES) at Re = $6\times {10}^6$ that generated unsteady turbulent flow around a rectangular cylinder with a controlled wrinkled surface and a 1:4 aspect ratio. A systematic study of the roughness effect was carried out by implementing different configurations of equally spaced grooves and bumps on the top-surface of the cylinder. Our results suggest that groove geometries reduce energy dissipation at higher rates than the smooth reference case, whereas bumped cylinders produce relative pressures characterized by a sawtooth pattern along the middle-upper part of the cylinder. Moreover, cylinders with triangular bumps increase mean drag and lift forces by up to 8% and 0.08 units, respectively, while circular bumps increase vorticity and pressure disturbances on the wrinkled surface. All of these effects impact energy dissipation, vorticity, pressure coefficients, and flow velocity along the wrinkled surface. Both the surface-manufactured cylinders and the proposed visualization techniques could be replicated in a variety of engineering developments involving flow characterization in the presence of roughness.

1. Introduction

Passive control measures are often employed to improve the aerodynamics around bluff bodies. Some case studies have required complex modeling approaches for the analysis and interpretation of results, but the majority require simple prototypes such as square and rectangular cylinders. This is the case in energy conversion in marine flow environments, the drag reduction of land vehicles, and the characterization of wind flow around buildings.

For ocean energy conversion, the aim is usually to enhance the amplitude of wave oscillation while maximizing the hydro-kinetic energy by varying the location of one or more cylinder structure [1,2]. This technique is called the VIVACE (Vortex Induced Vibration for Aquatic Clean Energy) converter, and was first reported in references [3,4,5,6]. Recent investigations still explore alternatives to optimize the VIVACE device, such as the implementation of secondary structures and channels to redirect the flow [2,7].

In contrast, the design of vehicles with silhouettes that mitigate drag forces commonly targets the delay of the boundary layer separation [8] and aerodynamic optimization through geometrical changes at the front [9,10] and rear of the vehicle [11,12,13]. An example of successful modifications is presented in Ref. [9], who tested a semi-trailer truck model in a wind tunnel set to a Reynolds number (Re) of $\approx 8\times {10}^5$ to quantify the efficiency of different wind deflectors and side fairings. This experiment enabled drag reductions of 19–29%, which could translate into 6–8% fuel savings [14]. Other modifications involved the addition of an ancillary to magnify drag, e.g., through external rear view mirrors which enabled increasing force by 3.23% in a modified truck-trailer, as discussed by Chilbule et al. [13]. These examples illustrate areas of opportunity around passive control mechanisms.

Non-applied studies dealing with geometrical changes explore meaningful geometrical changes in test models. For example, Jin et al. [15] proposed cutting the front corners off a squared cylinder to achieve a drag reduction of up to 48–61%; Mazellier [16] added adaptive porous flaps on the upper and lower sides of a square column to reduce energy absorption imparted by the fluid by up to 17–25%; and Malekzadeh et al. [17] placed barriers before a squared cylinder to reduce the mean drag coefficient and the standard deviation of drag and lift by up to 61%, 87%, and 94%, respectively. The first two studies cited here were performed at Re $\sim {10}^4$, whereas the latter was performed at Re $=5\times {10}^2$. Some of them proved to have potential applications in ocean engineering.

The manufacture of rough surface patterns for passive control has been explored in a variety of bluff bodies [18,19], more specifically on flat surfaces, as in references [20,21], and airfoils [22,23]. These methods have generated encouraging results: Mariotti et al. [18] achieved a significant delay in boundary layer separation and a reduction in pressure drag when implementing a single groove on the lateral surface of a boat-tailed axisymmetric body; the size of the grooves ranged between 1.5% and 3.3% times the characteristic length of the object. In turn, Qi et al. [19] reported disturbances in the wake of a circular cylinder leading to a reduction in drag and lift of up to 34% and 16%, respectively, when implementing triangular grooves with depths of 2–5% the size of the main body. Moreover, Aguirre-López et al. [24] scrutinized the use of triangular and squared patterns that mimic balconies as passive control measures of flow around tall buildings. These CFD simulations revealed that the facade with balconies displays zig-zag patterns in the pressure and vorticity fields along it, while improving the air circulation at the pedestrian level when the induced roughness covers the lateral walls of the building.

The present work builds on previous research focused on flow disturbances around squared and rectangular cylinders with surface modifications [24,25]. The object of study is a rectangular cylinder with a height-to-length ratio of 1:4 and a height-to-width ratio of 1:1. Re = $6\times {10}^6$ was set with the dual purpose of testing new size ratios for ocean engineering purposes while scrutinizing the aerodynamics of commercial buses [26] and other objects with an equivalent geometry and Re. To this end, we completed a systematic study considering bump and groove types to passively modify the aerodynamics of the object. A descriptive analysis of these configurations was carried out, involving pressure, velocity, and energy dissipation variables, while reviewing point, 1D, and 3D flow disturbances. Our numerical simulations were performed in OpenFOAM version 2012 [27], a well-known free open source CFD software that allows for simulating fluid–structure interactions.

The contents of this manuscript are presented as follows. Details of the implementation of geometries, meshing, and the numerical scheme are described in Section 2. The general results of this study are shown in Section 3, while they are discussed in Section 4. The concluding remarks are provided in Section 5.

2. Methodology

2.1. Geometries

The geometry of test cylinders was designed in the open source parametric 3D modeler FreeCAD [28] and subsequently imported into OpenFOAM [27] by means of a formatted stereolithography (STL) file compacted with its surfaceFeatureExtractDict tool. The systematic study consisted of 1 smooth rectangular cylinder and 8 modified prototypes.

All cylinders had the same dimensions, namely length $L=12$ m, height $H=3$ m, and width $D=3$ m, with an aspect ratio of 1:4 from a cross-sectional horizontal perspective and 1:1 for the front elevation (aligned with the main flow direction). The characteristic length of the test prototypes was therefore H.

Figure 1 shows the surface roughness implemented in the modified cylinders. The selected wrinkled geometries consisted of ten equispaced prisms of equal shape distributed across the topside of the cylinder. These protruded from the surface to create bumps while penetrating into it to create grooves. Each prism was contained within an imaginary circle with a radius equal to 5% of H, which replicated in the spanwise direction. The shape of each groove symmetrically or anti-symmetrically mirrored its corresponding numbered bump. Similar geometries for both bumps and grooves were designed with a size of 1%, not shown in the figure. The size of all these prisms was selected following Ref. [18] (1.5–3.3%) and Ref. [19] (2–5%), with the purpose of inducing significant flow disturbances.

2.2. Computational Array and Problem Description

The computational domain is shown in Figure 2. The origin of the global Cartesian system is located at the center of the cylinder and the flow runs parallel to the x-axis. Inflow and outflow boundaries were defined at $x=-6H$ and $x=11H$, respectively. In turn, symmetry boundaries were implemented in the front and back and up and down limits, located at $y=\pm 4.5H$, and $z=\pm 4.5H$, respectively. The walls of the cylinders abided by the non-slip boundary condition.

The simulations started with a free stream velocity equal to $U_{\infty }=10$ m/s in x-direction and a kinematic pressure $p=0$ ${\mathrm{m}}^2$/${\mathrm{s}}^2$ in the internal field, while a kinematic viscosity $\nu =5\times {10}^{-6}$ ${\mathrm{m}}^2$/s was considered. These settings corresponded to Re = $U_{\infty }H/\nu =6\times {10}^6$. Maintaining H and changing to standard atmospheric conditions (at 20 °C), the above Re led to an equivalent characteristic velocity of $U_{\infty }^{*}\approx 30$ m/s, leading to a Mach number (M) of $\mathrm{M}=U_{\infty }^{*}/c_s=$ 0.09, with $c_s$ as the speed of sound, which justifies the incompressibility of the simulated flow [29].

The boundary conditions presented in Figure 2 were computed by fixing $U_{\infty }$ and a zero-gradient condition for p at the inlet. In turn, a fixed value of $p=0$ ${\mathrm{m}}^2$/${\mathrm{s}}^2$ and a composed condition for U were set at the outlet; the latter consisted of switching between a zero gradient when the fluid escaped the domain and a fixed value of $\left|\mathbf{U}\right|=0$ when the fluid flowed into the domain. This was achieved by using the option inletOutlet from OpenFOAM. Finally, no-slip and zero-gradient conditions were set for U and p at the cylinder’s boundary.

The mesh contained a mix of structured and non-structured grids, which can be seen in Figure 3. The coarsest grid denoted a uniform-structured section tagged by Zone A, conformed with quasi-cube cells with an approximate volume of 1 ${\mathrm{m}}^3$, i.e., $H/3$ per side. Zone A was constructed with the blockMesh tool in OpenFOAM, while the rest of the grids were constructed by means of the SnappyHexMesh tool. Zone B also consisted of a uniform-structured grid but with higher refinement for capturing flow effects in the wake more accurately. The grids enclosing Zones C, D, and E were designed to further refine measurements when approaching the cylinder, with all of these planned to be non-structured. Zone E encloses the finest grid consisting of the cells intersected by the cylinder edges, allowing us to capture the geometry of the grooves and the bumps with an approximate resolution of $H/384$; in turn, Zone D was formed of cells that defined the imported STL-geometry, while Zone C was a scaled volume that diminished resolution while spanning from the ending of Zone D to a radial distance $2H$ from the cylinder’s walls; it also matched the refinement of Zone B at $0.5H$ from the walls. The referred software built the mesh in a way that, when two different refinement levels were superposed, the higher one was implemented. A space of one cell (of the lower level) was left at the transition between refinement levels.

2.3. Numerical Scheme

To carry out delayed detached eddy simulations (DDES), we used the Pressure Implicit with Splitting of Operators (PISO) solver in OpenFOAM, named pisoFOAM. This solver allowed us to reproduce the unsteady-incompressible-isothermal flow, taking advantage of the computational speed of Reynolds-average simulations (RAS), and switched to large scales when approaching the cylinder with large eddy simulations (LES), considering a smooth transition between both scales with the delayed-detached feature. Moreover, we merged the DDES scheme with the Spalart–Allmaras (SA) model, which was originally designed for simulating flows around airfoils, based on a equation of turbulent viscosity (${\nu }_t$) [30,31]. The general coupling of the numerical scheme is given by the SA equation:

$${\displaystyle \frac{D}{Dt}}\left(\rho \tilde{\nu }\right)=\nabla ·\left(\rho D_{\tilde{\nu }}\tilde{\nu }\right)+{\displaystyle \frac{g_{b2}}{{\sigma }_{{\nu }_t}}}\rho {|\nabla \tilde{\nu }|}^2+g_{b1}\rho \tilde{\zeta }\tilde{\nu }-g_{w1}{f}_w\rho {\left({\displaystyle \frac{\tilde{\nu }}{\tilde{d}}}\right)}^2+{\zeta }_{\tilde{\nu }},$$

(1)

where ${\nu }_t$ is pulled from a nonlinear function of the modified turbulent viscosity ($\tilde{\nu }$), so that ${\nu }_t=\tilde{\nu }{f}_{\tilde{\nu }1}$, with $f_{\tilde{\nu }1}\equiv f_{\tilde{\nu }1}\left(\tilde{\nu }\right)$ as a damping function depending on the distance to the nearest wall and other factors; $\rho$ is the density of the fluid; $\zeta$ refers to the shear stress; and $\tilde{d}$ defines the length scale. To switch between RAS and LES schemes, we used the following rule:

$$\tilde{d}=max\left[L_{RAS}-f_d,max(L_{RAS}-L_{LES},0)\right],$$

(2)

with the RAS scale ($L_{RAS}$) being a function of the distance from the analyzed cell to the closest solid wall and the LES scale ($L_{LES}$) being a function of the local grid spacing and a modified low Reynolds number. The smoothness of the transition was controlled by the hyperbolic-tangential function $f_d$.

In this context, our solution methodology allowed the discretization of the Navier–Stokes equations and the SA equation via the finite volume method (FVM), taking into account a conditionally stable second-order implicit-backward interpolation for time schemes, Gauss linear interpolation for gradient schemes, Gauss linear limited interpolation for divergence and Laplacian schemes, and linear interpolation to transform the cell-center quantities to face centers.This produces a system of algebraic equations of the form $A\mathbf{x}=\mathbf{b}$, where $A\equiv$ coefficient matrix, $\mathbf{x}\equiv$ vector of variables, and $\mathbf{b}\equiv$ source vector, which is solved by pisoFOAM through an implicit predictor and multiple explicit corrector steps, seeking to obtain close approximations of the exact solution of the difference equations at each time-step [32]. In this part, the algorithm works in a fully implicit pressure scheme with the coupling of the velocity and pressure equations being handled through the iterations [33,34]. We set a maximum of three corrections without extra-correction for mesh non-orthogonality, while the resulting algebraic equations were solved via different methods: a generalized geometric-algebraic multi-grid (GAMG) solver method for pressure, with a Gauss–Seidel method as a smoother; a preconditioned pipelined conjugate residuals (PPCR) solver to override solution tolerance for the final pressure; and a symmetric Gauss–Seidel solver for the rest of the variables.

Regarding the values of parameters, the turbulent viscosity was initialized with ${\nu }_t=0$ in the entire domain and calculated (boundary condition) using the function nutUSpaldingWallFunction at the cylinder’s walls. In turn, the modified turbulent viscosity was set to $\tilde{\nu }=5\nu$ in the entire domain, except for cells near to the cylinder’s wall, in which $\tilde{\nu }=0$ ${\mathrm{m}}^2$/s, as usually followed by the Spalart–Allmaras model [35,36]. For consistency with the high Re we were dealing with, a uniform turbulence intensity of $I=15$% was also defined at the inlet [37], which led to a turbulent energy of $k=(3/2){\left(U_{\infty }I\right)}^2=3.375$ ${\mathrm{m}}^2$/${\mathrm{s}}^2$.

A full definition of all variables in Equation (1) and other implicit parameters in the model is given in references [35,38], including the parameters cited above, the solver mechanisms, the meshing techniques, and the map default values of the motorbike example in OpenFOAM-v2012 (Directory: OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pisoFoam/LES/). The reader is referred to the OpenFOAM tutorials [39] for more details on similar flow simulations, and to our previous works [24,25], in which our current computational–mathematical array is applied and validated by means of a resolution study dealing with a squared cylinder at Re = $2.14\times {10}^4$ and a comparison with experimental results focused on buildings subject to a flow regime between Re = $1.59\times {10}^7$ and Re = $4.01\times {10}^7$, which covers the Re used in our current case study.

2.4. Validation and Nomenclature

With the aim of standardizing data representation, the units in the graphics and results are dimensionless, e.g., $(X,Y,Z)=(x,y,z)/H$, $U\equiv (U_x,U_y,U_z)=(u_x,u_y,u_z)/U_{\infty }$, $\mathsf{\Omega }=\left(\omega H\right)/U_{\infty }$ (vorticity), $N_t={\nu }_t/\left(U_{\infty }H\right)$, and $T=\left(t{U}_{\infty }\right)/H$. Drag $Cd$ and lift $Cl$ dimensionless coefficients were directly obtained from OpenFOAM by referring to the library libforces.so of the function forceCoeffs in the directory controlDict. The relative-pressure coefficient was computed using the known formula

$$Cp={\displaystyle \frac{p-p_{\infty }}{0.5\rho U_{\infty }^2}},$$

(3)

where the freestream pressure $p_{\infty }$ is consistent with the initial pressure condition.

Table 1. Number of cells of each test case.

Geometry	# of Cells
Smooth (benchmark)	1,695,959
B1	2,064,356
B2	2,050,865
B3	2,176,280
B4	2,457,100
G1	2,035,171
G2	2,021,147
G3	2,152,050
G4	2,382,913

shows the number of cells contained within the computational arrays that result from implementing the geometries described in Section 2.1 using the mesh parameters discussed in Section 2.2. According to this, wrinkled cylinders occupy a larger quantity of cells given the refinement needed to shape each bump/groove. In this sense, meshes with circular modifications (B4 and G4) increase cell numbers by up to 41–45% compared to the smooth case. Below these, we see the squared modifications (B3 and G3), with a 27–28% increase, and the two types of triangular modifications (B1, B2, G1 and G2), with a 19–22% increase above the smooth case.

In addition to the validation cases mentioned in Section 2.3, we initialized the present simulations by pre-testing the case B4 with three lower refinements near the object: a coarse mesh (with 55% fewer cells than our smooth case, and a maximum resolution of $H/24$), a medium mesh (with 50% fewer cells and $H/48$ in resolution), and a fine mesh (with 10% fewer cells and $H/96$ in resolution), and a same-level refinement but adding more temporal refinement (from $\mathsf{\Delta }T=1\times {10}^{-3}$ units to $\mathsf{\Delta }T=6.67\times {10}^{-4}$). Coarse, medium, and high meshes apply their corresponding maximum resolution to Grids B, D, and E (see Figure 3); then, the high mesh differs from the extra-high one only in the refinement level for Grid E. Figure 4 shows the convergence of our computational array when calculating the aerodynamic coefficients; Ext. high2 array refers to the extra-high one with $\mathsf{\Delta }T=6.67\times {10}^{-4}$. Although convergence is achieved with the high mesh, we selected the extra-high (Ext. high) mesh as an extra measure of security to adequately capture the wrinkled geometries of the non-smooth cases. From this point on, we will use the term “benchmark” to refer to the wrinkle-free cylinder modeled by using the extra-high mesh.

The validation of the model was performed in terms of the energy contents of the drag and lift pressure fluctuations across a frequency range f for the case B4, using the Ext. high mesh. Figure 5 shows the autoscaled power density function ($PSD$), in which we can observe a broadband energy that characterizes turbulent flow. These results show how the low-frequency fluctuations associated with larger eddies or gusts are followed by the inertial range (highlighted in black), covering the range $f={10}^{-4}$–${10}^{-3}$ Hz. The spectral power distribution satisfies the transference of energy according to Kolmogorov’s scaling law. We can also note that the energy dissipation rate takes place at the highest frequencies, obeying a faster decay due to the turbulence of smaller eddies.

Then, by using the OpenFOAM utility mapFields, we mapped the high-mesh solution at time $T=200$ to start the fine tests (those referred to in Table 1). This enabled us to avoid an increase in the transition interval for the stabilization of the flow, which permitted a larger sampling interval. The total time simulation of the fine tests also matched 200 units. Figure 6 allows us to exemplify this by showing the variation in the drag $Cd$ (a) and lift $Cl$ (b) coefficients obtained for the fine-smooth case in superposition to those calculated with the modified geometry. The sampling interval started from $T=20$ to ensure that the flow was stable for all runs as well as for both the $Cd$ and $Cl$ series, which were obtained by referring to the function forceCoeffs of the library libforces.so, evaluating Equations (4) and (5):

$$Cd={\displaystyle \frac{Fd}{0.5\rho DH{U}_{\infty }^2}},$$

(4)

$$Cl={\displaystyle \frac{Fl}{0.5\rho DH{U}_{\infty }^2}},$$

(5)

where $Fd$ and $Fl$ are the force acting on the building in the streamwise and transverse directions, respectively.

The Courant number was kept below 1 with the corresponding grid and time discretization of the Ext. high mesh. In each time step and each velocity component, the convergence criteria were set at $1\times {10}^{-6}$ for pressure and $1\times {10}^{-7}$ for turbulent kinematic viscosity in their respective units. Each simulation ran in parallel across the 24 subdomains, by using the scotch method from the decomposeParDict tool in OpenFOAM. We used a cluster with Intel^® Xeon^® CPU E5-2680 v3 at 2.50 GHz with 48 Cores and 131.072 GB RAM. Under these conditions, the computational time was around 2 days for the benchmark and up to 3 days for the wrinkled cases discussed in the following sections.

3. Numerical Results of the Surrounding Flow

Figure 7 shows the average velocity maps $U_{mean}$, combined with contour maps of pressure coefficients $C{p}_{mean}$. For all geometries, the sharpness of the main body causes the flow to display a separation at the top-left corner of the cylinder ($X=-2$), followed by its corresponding return halfway up the cylinder ($X=0$). This involves a large zone of vorticity that is reflected by the closed iso-pressure contours that extend up to the position corresponding to the fourth groove/bump ($X\approx -0.5$) and more than one characteristic length D in the upward direction. Outside this zone, the main time-averaged disturbances (large and small) caused by the geometries occur.

On a large scale, all of the modified geometries disturb the flow symmetries seen in the benchmark around the cylinder. In general, the effects produced by bumps are greater than those produced by grooves. Despite the fact that the bumps caused a break in the iso-pressure contour of −0.8 units (sky blue contours) in the vortex zone mentioned above, the flow at the front and bottom of the cylinder showed no significant disturbances in the immediate and middle zones, so we will focus on the upper part of the cylinder.

On a smaller scale, around the modified surface topology, it is also observed that the presence of bumps caused larger, perhaps more interesting, disturbances to the flow than the induced grooves. The latter fold the near-zero contour lines inside the six rear grooves. In contrast to this, it is worth highlighting an up–down effect on pressure that became particularly visible around the last five bumps. For triangular cases (B1 and B2), high-pressure zones (red contours) are located before every bump, while low-pressure zones (blue contours) clearly developed around their rear, which was assumed to be related to a softening flow separation. In turn, for the squared and circular cases (B3 and B4), high pressures share the front of each bump with zero-level pressures (dark contours), while very low pressure contours predominate on the exact position and rear side of the bumps. Contrasting with triangular cases, this observed transition of pressure zones suggests a quick flow separation caused by the bluff shape of the bumps.

4. Discussion and Deeper Analysis

4.1. Point Effects: Flow–Structure Interaction

In Figure 8, we plot the time-averaged $C{p}_{mean}$ versus standard deviation $C{p}_{std}$ coefficient pressure calculated by (3), to highlight the observed small-scale effects. The graphics span half of the cylinder surface, from the front to the rear, as drawn in the picture. For all geometries, the highest pressure $C{p}_{mean}=1$ occurs on the front side, and then reduces to a value $C{p}_{mean}=-0.5$ at the top-left corner where the flow separates. Then, pressure reduction continues until $X=-1.3$, where it reaches different local minimums depending on the geometry. From that point, we see a pressure increase until $C{p}_{mean}=0$ at $X=-0.3$ due to the reattachment of the flow. For the benchmark (a), the pressure stabilizes without vorticity until the rear side ($X=2$), while the associated standard deviation only increases in the transition zone $-1.3<X<0.3$.

As pointed out above, cases G1–G4 only disturb the mean flow at the location of the grooves (small-local changes in $C{p}_{mean}$), but they produce larger variations in $C{p}_{std}$ than the benchmark around the transition zone and beyond. As observed in Figure 7, this performance originates past the first groove where small eddies appear, without altering the overall pressure pattern established by the main flow, which is equivalent to the effect produced in a continuous flow by a small cavity, acting as a controller of flow separation, avoiding a cascade effect. This process is replicated across the grooved region since the main flow dominates, which is the reason why the pressure envelope resembles that of the benchmark.

In turn, the structured contour lines observed with the configurations B1–B4 show a sawtooth pattern of $C{p}_{mean}$ from the middle ($X=0$) to the rear side ($X=2$) of the cylinder. The smooth transition to local highs of the sawtooth pattern takes place over the zones with a flat surface, whereas the transitions to local lows are steeper over the bumps. In agreement with the observations illustrated in Figure 7, the pattern is smoother and lower in amplitude for the triangular cases (B1 and B2), where the maximum disturbances reached are about 0.2 dimensionless units, which fluctuate symmetrically with respect to the zero line, while their standard deviations oscillate between approximately 0 and 0.15. For cases B3 and B4, the largest local decreases in $C{p}_{mean}$ occur with the highest values of $C{p}_{std}$, overall increasing their magnitude with their proximity to the last bump to reach $C{p}_{mean}=-0.6$ and $C{p}_{std}=0.25$ for B3, and $C{p}_{mean}=-1$ and more than $C{p}_{std}=0.3$ for B4. Furthermore, the sharpness of the rough geometries leads to the highest $C{p}_{std}$ values occurring at the $C{p}_{mean}$ peaks (before bumps) for B1–B2, whereas they occur at the $C{p}_{mean}$ valleys (over bumps) for B3–B4.

This performance suggests an accumulation of uncontrolled flow redirection caused by the bluff bumps. However, the separation between bumps is large enough to maintain flow separations in a local or point scale, the reason being that the general sawtooth structure is always in phase according to the equally spaced roughness at the middle-rear side, hence the quasi-harmonic pressure configuration.

4.2. Three-Dimensional Effects: Aerodynamic Coefficients and Results Compilation

Figure 9 shows a different perspective regarding the small and large disturbances observed in the mean aerodynamic coefficients $C{d}_{mean}$ and $C{l}_{mean}$. No geometrical change reduced drag with respect to the benchmark. The grooves caused slight increments of no more than 1%, whereas higher increments ranging from 4% to 8% were reached in the presence of bumps, namely 4% (B3), 5% (B4 and B2), and 8% (B1). Similar effects are seen for the lift coefficient. $C{l}_{mean}$ increases by less than 0.02 units for groove cases, while it increases by 0.06 units for B3 and B4, 0.07 for B2, and 0.08 for B1. In this figure, error bars are also plotted to compare the standard deviation calculated for the four coefficients. These results do not show a dependence on geometrical modifications: all $C{d}_{std}$ values vary between 1.5% and 2% of the $C{d}_{mean}$-benchmark value, while the $C{l}_{std}$ values fluctuate by about 0.1 units. The side force is not shown but there are no significant disturbances there due to the main flow being aligned with the x-direction.

The patterns observed for $C{l}_{mean}$ and $C{d}_{mean}$ are, to an extent, expected given the results shown in the previous sections, where we observed that grooved geometries did not cause large disturbances on a local or large scale. This, added to the fact that the benchmark is fully symmetric, leads to the near-zero effects in $C{l}_{mean}$. On the other hand, the larger pressure effects on the top surface exerted by the bump geometries, Figure 8, differ enough from the benchmark to cause the largest variations.

Regarding $C{d}_{mean}$ disturbances, we stated that when examining Figure 7, Figure 8, Figure 9 and Figure 10, there are no clear pressure disturbances on the front and rear sides of the cylinder, meaning that changes in drag may originate from local effects on the wrinkled surface. This assumption is supported by the logical observation that extra added obstacles (bumps) increase the characteristic length D and consequently drag, such as the ancillary modifications commented on in Ref. [13]. This also agrees with a close inspection of the plots in Figure 8, in which $C{p}_{mean}$ takes relative positive values at the immediate front cells of each bump, whereas it takes near-zero values at their adjacent rear cells. These pressure differences are similar for all the bump geometries, which would explain the low variation in $C{d}_{mean}$ between these modifications, from 4% to 8%, despite the higher pressure differences caused by B3 and B4 in the flat zones that affect the lift. In other words, the inclination of triangular shapes leads to a greater surface area being exposed in the mainstream direction, allowing for a higher drag force than the square and circular bumps.

These deductions seem contradictory to previous reports addressing grooved surfaces with a similar relative size, which reported notable drag and lift reductions at lower Re (see references [18,19] in Section 1). The difference is related to the fact that these modifications disturb the pressure at the rear side. This could happen because (i) a more aerodynamic main body allows a visible effect due to the equispaced wrinkles, even for small modifications (see, for example, Figure 5 in Ref. [18] (p. 525), and Figure 12 in Ref. [19] (p. 14)) or (ii) the wrinkle is located at a position at which it can provoke a large effect (see, for example, the contrast between the two locations of the same wrinkle in Figure 9 in Ref. [18] (p. 529)). In the following subsection, we further discuss the contrast observed when implementing grooves and bumps into the model, while exploring and reporting alternative configurations for special purposes.

4.3. One-Dimensional Effects: Passive Control Along the Wrinkled Side

As described above, the wrinkles on the cylinder configure particular pressure and velocity fields. Therefore, it is of interest to understand how these two variables relate to vorticity. In an attempt to answer this question, this subsection examines the relationship between the pressure and velocity fields through $U_{mean}$, $C{p}_{mean}$, $U_{std}$, $C{p}_{std}$, and ${\mathsf{\Omega }}_{mean}$, in the immediate zones of the top cylinder surface shown in Figure 8. In addition, the relationship between these flow characteristics and the energy dissipation rate $ϵ$ is explored.

Employing Bernoulli’s equation and Formula (3) for incompressible flows, we derived the following relationship between $C{p}_{mean}$ and velocity,

$$Cp=1-U^2.$$

(6)

Then, when the local velocity U increases, $Cp$ decreases.

Figure 10 depicts the joint behavior of $C{p}_{mean}$ and $U_{mean}$ as a function (in color) of the auxiliary variable of distance $\xi$ for each geometry considered in the simulation. This variable represents the minimum distance from the recording point of the mid-top surface to the nearest wrinkle. In this way, the points plotted in cyan are inner grooves or bumps, respectively; points in blue fall within the flat zones of the top surface; and those in black correspond to the front and rear sides. For either configuration, we can see that $U_{mean}=0$ and $C{p}_{mean}=0.84$ on the left end (front) of the cylinder at $(X,Y)=(-2,0)$; this is followed by a rapid increase in $U_{mean}$ and a decrease in $C{p}_{mean}$ as functions of X for values of X around $-2$ until they settle down to what we call the “stabilization” zone. This decay in the value of $C{p}_{mean}$ is also observed in Figure 8. After moving away from the stabilization in the right end of the cylinder and for coordinates in the neighbourhood $X=2$, $C{p}_{mean}$ and $U_{mean}$ exhibit values around $-0.1$ and $0.1$, respectively, according to Figure 7, which was expected since the most dynamic flow–cylinder interactions occur at the left end, and hence we only see residual motion of the fluid at the right end.

Figure 10a, corresponding to the benchmark, exhibits the most regular pattern between $C{p}_{mean}$ and $U_{mean}$; it resembles a sigmoid curve. In this case, $C{p}_{mean}$ decreases as $U_{mean}$ increases and vice versa, except in between the inflection points around $X=-1.2$ (around the second wrinkle) and $X=0$ (between the fourth and fifth wrinkle). The coordinates X of these inflection points match the inflection points of $C{p}_{mean}$ observed in Figure 8a. These points define three distinguishable clusters in Figure 10a: the mean velocity $U_{mean}$ increases before the first inflection point and after the second one, while it decreases at the intermediate X-coordinates. Therefore, $U_{mean}$ inversely mirrors the behavior of $C{p}_{mean}$ for the benchmark. In Figure 10b,d,f,h, corresponding to grooved cylinders, the clusters limited by the inflection points can also be identified. In these cases, the mean velocity decays around the grooves, taking values around $U_{mean}\in (0.1,0.2)$ for G1 and G2 (except in the 10th groove in G2) and around $U_{mean}\in (0.1,0.3)$ for G3 and G4. The patterns observed replicate across the tested groove geometries, and the variability inside each cluster is roughly the same; thus, we conclude that the forms in each of these plots correspond to the sigmoid curve of the benchmark, disarranged by the grooves. This observation is supported by the results observed in Figure 8 for the $C{p}_{mean}$ domain. In the case of the latter, corresponding to Figure 10c,e,g,i, the sigmoid curve is altered to a higher degree. In this case, the clusters induced by the inflection points can also be identified, but they overlap. Despite this, the location of points in blue indicates that the velocity–pressure relationship is similar for all of these geometries in the flat and bumped zones, except for B4. It is worth noting that the mean velocity in the bumps is highly unsteady, with values of $U_{mean}\in (0.1,0.5)$ for $B1$–$B3$ and $U_{mean}\in (0.2,0.9)$ for $B4$. Notably, the case shown in Figure 10i, for the geometry $B4$, indicates that circular bumps induce the most chaotic behavior in the flow. Therefore, in spite of $B4$, the presence of wrinkles produces lower values of $U_{mean}$ compared to the benchmark.

Figure 11 illustrates the relationship between $C{p}_{std}$ and $U_{std}$ as a function of the distance $\xi$ for each geometry considered in the simulation. Figure 11a shows the benchmark. In this plot, we can identify four data clusters around the coordinates X: (i) for X very near $X=2$ and $X=-2$; (ii) X between −2 and −0.75; (iii) X between −0.75 and 0; (iv) X between 0 and 2. The first cluster forms at both ends of the cylinder and relates to the closest point to the origin (with the lowest variation: for the left end, $C{p}_{std}$ and $U_{std}$ mostly take values near 0, while for the right end, $C{p}_{std}\approx 0.025$ and $U_{std}\in (0.025,0.075)$. Note that the data outside of the first cluster describe a curve. The second cluster is located at the low end of the curve and falls further away than the third wrinkle in relation to the point where the curve starts to rise at the $U_{std}$ coordinate and backward at the $C{p}_{std}$ coordinate. The third cluster spans from the end of the second cluster to the point where the curve moving along the $U_{std}$ coordinate starts descending after the position of the fifth wrinkle. The fourth cluster starts at the tail of the curve, and it moves downward and reaches the peak of the curve, where we see an abrupt decay in the velocity’s standard deviation inside this cluster between the position of the fifth and sixth wrinkles; Figure 7 suggests that this is a transition zone between the highly dynamic first half of the cylinder and the less active second half, as pointed out in Section 4.1. Figure 11a also shows that the standard deviation of the pressure coefficient $C{p}_{std}$ is, on average, higher for the first part of the cylinder, which contrasts with the behavior of the second part of the object.

Figure 11b,d,f,h seem to display a distorted version of Figure 11a, although one could argue that the clusters recognized in Figure 11a are also present in these plots, but showing higher local variance in the clusters. Another consequence of the referred distortion is that measurements in the grooved cylinders follow the patterns observed in Figure 11a: the standard deviations $C{p}_{std}$ increase from 1 to 3 before starting to decrease. On the other hand, the standard deviation measured in the neighborhood of the grooves depict a line; for G1–G2, these measurements are below the mean line of the points, while for G3–G4, they are near the mean line. The difference in terms of $U_{std}$ between the first and second half of the cylinder is higher than in the benchmark. Finally, in Figure 11c,e,g,i, it is not possible to identify any clear cluster associated with bump geometries. Most of the data points in these figures fluctuate around a backbone line, except for B4, in which case a line seems to fit well for the cases $C{p}_{std}<0.20$ but not for $C{p}_{std}>0.20$. Measurements in the bumps there are not described by a line. Similarly to the records taken for mean values, Figure 11i suggests that the geometry B4 induces the most turbulent behavior.

Moving on, there are different ways to compute the energy dissipation rate, depending on the restrictions on the flow. Since we used the Spalart–Allmaras model (1), this study employs the turbulent viscosity ${\nu }_t$ as a proxy for such a quantity, which is a non-physical variable that qualitatively characterizes the turbulent flow, by relating $ϵ$ to $\mathsf{\Delta }{\mathbf{U}}_{mean}$ according to

$$ϵ\sim {\nu }_t·{\left(\mathsf{\Delta }{\mathbf{u}}_{mean}\right)}^2,$$

(7)

see Landau and Lifshitz [40]. Another formula to consider is given by Ref. [41],

$$\begin{array}{c}\hfill \tilde{\nu }=\sqrt{{\displaystyle \frac{3}{2}}}{u}_{mean}I{l}_t,{\nu }_t=\tilde{\nu }{f}_{\tilde{\nu }1},\end{array}$$

(8)

where $l_t$ is the turbulent length scale (or the characteristic turbulent length), while I, $\tilde{\nu }$ and $f_{\tilde{\nu }1}$ are defined in Section 2.

In this context, Figure 12 presents the joint behavior of the mean energy dissipation, expressed by $N_{t_{mean}}$, the mean velocity $U_{mean}$, and the mean vorticity ${\mathsf{\Omega }}_{mean}$. The most striking result here is that, in all the cases, the data points $(U_{mean},N_{t_{mean}})$ lie over rays converging to the origin, with a gradient bounded by the ray associated with the benchmark, which suggests that the wrinkles control the dissipation of energy on the top surface. In relation to this, Equation (8) implies that variations in the turbulence intensity I, length scale l, and function $f_{\tilde{\nu }1}$ determine the various gradients shown in Figure A1. Moreover, for the benchmark, most of the measurements follow a line, with the exception of data points that, once tracked, are related to the front $X=-2$ and rear sides $X=2$, plus some other points around the inflection observed at $X=0$, where the reattachment of the flow is discussed in Section 3 and Section 4.1. This linear relationship establishes a direct link between the mean velocity and the energy dissipation rate. Furthermore, we also observe that vorticity (color scale) increases proportionally to the value of the velocity.

Regarding prototypes with grooves, Figure 12b,d,f,h show a more established relationship between $N_{t_{mean}}$ and $U_{mean}$. Additionally, the measurements at the position of the grooves reveal an energy dissipation $N{t}_{mean}<2\times {10}^{-6}$, except for G3 and one value in G2. Likewise, for geometries with bumps, Figure 12c,e,g,i show that the measurements of $(U_{mean},N{t}_{mean})$ also define less dense and better defined rays. This added to the fact that all modified geometries reduce dissipation, and hence vorticity is also affected by the presence of wrinkles, showing higher values at the right ends of the lower rays. This performance is replicated in practically all case studies, at similar rates, only changing with the graphic coordinates, which suggests that vorticity is strongly related to flow velocity and energy dissipation. Such behavior yields a well-organized and stable flow, as seen in Figure 7. Notably, the case of B4 exhibits higher mean velocities (and vorticity) than in the other cases; besides, the data points define two rays, with all the measurements in the bumps being located on the lower ray.

To synthesize the observed performance, we rewrite Equation (8), and complement it with dimensional analysis, so that

$${\nu }_t=c_1{l}_t{u}_{mean},\mathrm{with}{c}_1=I{f}_{\tilde{\nu }1}\sqrt{{\displaystyle \frac{3}{2}}},$$

(9)

which expresses the linear relationship between turbulent viscosity and mean velocity shown in Figure 12, where the gradient of the rays is defined by $c_1{l}_t$. It is worth noting that, since we are plotting the nearest cells to the objects, $f_{\tilde{\nu }1}$ does not change the value of the gradient, whereas Equation (9) suggests that the surface roughness induces small turbulent length scales $l_t$. This also holds for more defined rays determined across grooved geometries and for more disorganized clusters for the bumped cases.

Moreover, Figure 12 shows that the greater the gradient of the ray, the greater the vorticity, which can be inferred from Equation (9), since

$$\omega \sim {\displaystyle \frac{u}{l_t}}$$

(10)

and consequently,

$${\nu }_t\sim c_1{l_t}^2{\omega }_{mean}.$$

(11)

In spite of the identified proportionality between ${\nu }_t$ and $\omega$ present in Equation (11), the cases in Figure 12 do not show such performance, as we observe an apparent inverse proportionality instead. The reason for this seems to be related to the high intensity of vortex generation at low turbulence levels, i.e., $\omega$ is more dependent on ${l_t}^2$ than on ${\nu }_t$, which indicates the presence of eddies instead of large vortexes.

In summary, Equations (10) and (11), inferred from our results, suggest that vorticity could be approximated with a function that depends on turbulent viscosity and flow velocity, while being independent from the wrinkled geometry. Figure 13a, which results from super-positioning the time-averaged values of all the geometries presented in Figure 12, supports the established relationship. Figure 13b shows that such a behavior is replicated for the time-averaged-variation values. The reader is referred to Figure A1 in Appendix A for additional evidence involving all cases.

4.4. Potential Applications and Future Work

The results obtained exemplify the effects that modified geometries could have on the aerodynamic performance of bluff bodies and their further engineering applications, as discussed in Section 1. The pressure patterns taking place on the top side of the cylinder could promote the use of composite structures in ocean engineering that could help to identify the optimum location of soft and tough materials. Further technological improvements could be informed by the patterns observed in the standard deviation of pressure fluctuation and including the vorticity and dissipation of energy charts shown in Figure 12 and Figure 13. Furthermore, bump geometries could be manufactured on strong structures to define the degree of energy that these could absorb, while groove geometries could be useful for minimizing material cost since although they modify the aerodynamic performance, they do keep the force fields close to that of the benchmark. See Figure 8 for an example of the above observations. It is noticeable that grooves induce lower pressure disturbances than bumps over the region of reattachment; see, for example, Figure 8h,i (G4 vs. B4), where the former registers a peak $C{p}_{std}$ of around 10%, while the latter reaches 30% while crossing the zero line sporadically. These pressure fluctuations could induce excessive vibrations on flexible surfaces and membrane-like and cladding panels, none of where were modelled in this study.

Additionally, B1–B4 bump configurations could be considered for the prototyping of bus ancillaries, such as air conditioners, ceiling fans, or escape windows. The extra spaces on the bus would require less than 8% in drag. It should also be noted that our models are coarse approximations of the actual shape of a bus and that the effects produced by the floor wall are not considered. In turn, G1–G4 would be the choice in architectural design when the intention is not to significantly vary lift or drag.

Although the study of the wake is outside the scope of the present work, further research could follow on the positioning of objects in the far wake, given the asymmetric inflow that they would perceive. Future research could also involve big data analysis to characterize the time fluctuations of the flow, for which the addition of non-equispaced wrinkles could help to identify additional controlling parameters and their effect on the data recorded.

5. Conclusions

This paper presents a systematic study of flow around a rectangular cylinder at Re = $6\times {10}^6$, considering eight types of surface roughness on its top-surface. The wrinkled surfaces include prismatic grooves/bumps equispaced in the main flow direction.

Our concluding remarks about this study are summarized as follows:

• All of the proposed rough surface configurations cause local flow disturbances such as separation and reattachment, which at the micro-scale produce a sawtooth pressure pattern of pressure along the wrinkled surface. Therefore, variations in the aerodynamic coefficients with respect to the benchmark derive from mechanisms that determine the accumulation of pressure around the grooves/bumps. The largest increments in the drag are reached through the implementation of triangular bumps in model B1, which amounts to 8% with respect to the benchmark. In turn, the mean lift increases by 0.08 units, whereas no significant changes to the standard deviation were observed for either the drag or lift forces.
• This investigation puts forward an alternative visualization technique to scrutinize time-averaged and time-averaged-variation spaces of pressure–velocity, and turbulent viscosity–velocity corresponding to the wrinkled surface. For the pressure–velocity spaces, we observed two different clusters when plotting data points recorded on geometries with grooves and bumps. For the viscosity–velocity spaces, the curves fit a straight line for the benchmark, but split into multiple rays when introducing any type of roughness, with the grooved geometries G1–G4 being the least affected.
• The differences in flow patterns and aerodynamic coefficients determined in this study could help in tackling problems induced by the flow–structure interaction. We speculate that similar methodologies could benefit other engineering applications such as mechanical, maritime, and beyond.