Study of Surface Roughness Effect on a Bluff Body—The Formation of Asymmetric Separation Bubbles

Turbulent flows around bluff bodies are present in a large number of aeronautical, civil, mechanical, naval and oceanic engineering problems and still need comprehension. This paper provides a detailed investigation of turbulent boundary layer flows past a bluff body. The flows are disturbed by superficial roughness effect, one of the most influencing parameters present in engineering applications. A roughness model, recently developed by the authors, is here employed in order to capture the main features of these complex flows. Starting from subcritical Reynolds number simulations (Re = 1.0 × 105), typical phenomena found on critical and supercritical flow regimes are successfully captured, like non-zero lift force and its direction change, drag crisis followed by a gradual increase on this force, and separation and stagnation points displacement. The main contribution of this paper is to present a wide discussion related with the temporal history of aerodynamic loads of a single rough circular cylinder capturing the occurrence of asymmetric separation bubbles generation. The formation of asymmetric separation bubbles is an intrinsic phenomenon of the physical problem, which is successfully reported by our work. Unfortunately, there is a lack of numerical results available in the literature discussing the problem, which has also motivated the present paper. Previous study of our research group has only discussed the drag crisis, without to investigate its gradual increase and the change on lift force direction. Our results again confirm that the Lagrangian vortex method in association with Large-Eddy Simulation (LES) theory enables the development of two-dimensional roughness models.


Introduction
The technologic and science enhancement in many fields of knowledge is because of deep analysis and consequent comprehension of flow fundamentals around bodies of arbitrary shape. Many of the advances observed in ocean, civil and wind engineering have been possible because of studies conducted over time in the field of aerodynamics. The importance of knowing and mastering this subject can be easily recognized when it is observed that, most of bodies present in situations of practical interest for engineering, are exposed to air or water flows. These flows are typical examples of fluid-structure interaction problems, where the transition to turbulence is undeniably a complex phenomenon of nonlinear hydrodynamics. Such phenomenon arouses great scientific interest, has a Re = 3.4 × 10 5 of the base pressure coefficient; during this phase, the lift coefficient, C L , calculated from the pressure distribution, was about 1.3. The separation bubble on the other side formed at Reynolds number flow of Re = 3.8 × 10 5 , and the lift coefficient was C L =0.0, being measured as soon as this happens and remains in this way, as the Reynolds number increases.
Guven et al. [7] presented in one same diagram the results obtained in their experimental tests with other results available in the literature [3,8,9] for circular cylinders with different roughness heights. It was showed scatter results with up to 60%, even for the same roughness and Reynolds numbers. This happens because, although the Reynolds number governs this kind of flow, it is so difficulty to guarantee that two experimental tests have exactly the same roughness surface; besides that, there are some influencing parameters which can cause a premature initiation of a flow regime or even its modification.
According to Zdravkovich [10], the above referred influencing parameters are: (a) free stream turbulence, described by its intensity and scale; (b) surface roughness effect, described by its mean relative surface roughness, ε * /D * , and its texture; (c) wall confinement (blockage effect); (d) aspect ratio; (e) presence of other boundary near the body; (f) end condition of the body; (g) body oscillation effect. The free stream turbulence and surface roughness influence are the most common disturbances in all engineering applications. The wall confinement and the aspect ratio are important influencing parameters that interfere in the comparison between experimental and numerical results, especially when the surface roughness effect is an influencing parameter too [17,18]. These features help us to understand the scatter in the experimental data presented by Guven et al. [7].
In a recent paper, Bimbato et al. [14] developed a Lagrangian vortex method using two-dimensional roughness model, which was based on the physics involved in turbulent boundary layer flows. The roughness surface effect was simulated without make changes on geometry of the body, what means that the investigated circular cylinder surface remains smooth. The key idea is that a point set strategically located close to the body surface injects momentum in its boundary layer aiming to represent surface roughness effect. Their strategy successfully shows that the separation of boundary layer is delayed. They studied the incompressible flow around both smooth and rough circular cylinders at Re = 1.0 × 10 5 , and it was only reported the drag crisis. Based on aerodynamic loads computations, it was concluded that their methodology, starting from a typical subcritical Reynolds number flow, was able to capture supercritical flow patterns. Overall, the drag coefficient computed just fall and it did not present the gradual increase, which is a typical feature of supercritical flows. Also, it was not presented results about the lift force, which would be interesting since a change on the lift force direction is expected [2].
The present paper, therefore, contributes to the literature using the same roughness model developed by Bimbato et al. [14] in order to capture physical details concerning the complex flow around a single rough circular cylinder. Put in other words, the main contribution here states on physical discussions involving the drag and lift forces connected with the occurrence of asymmetric separation bubbles generation. The generation of asymmetric separation bubbles is an intrinsic phenomenon to the physics of the problem [2]. In Section 4, it will be reported that our numerical results have captured that phenomenon, which has rarely been discussed by other numerical investigations. The results also confirm that two-dimensional Lagrangian vortex method with both LES and roughness models is able to capture not only the drag crisis, but also the non-zero lift force, which is typical of critical Reynolds number flows, and the gradual increase on drag coefficient, which is a typical feature of supercritical flows. Furthermore, our results show a delay in the boundary layer separation and also a stagnation point displacement.
Finally, it can be said that numerical simulations using Lagrangian vortex method have been employed to study problems in several areas of engineering. It is possible to mention a large number of studies, from the classics e.g., [19][20][21][22] to the most current ones e.g., [23,24]. Despite the high Reynolds numbers involved in practical engineering problems, the use of two-dimensional numerical simulations is necessary not only to develop and validate computational codes, but also to help Energies 2020, 13, 6094 4 of 20 understand the complex physical phenomena involved in such flows, and also to help conservative designs for engineering applications. A wide range of two-dimensional studies can be found in literature e.g., [14,17,18,[25][26][27][28][29][30][31] which includes relevant contributions of our research group. Figure 1 illustrates the two-dimensional fluid domain, Ω, to be studied, and defined by surface S=S 1 ∪S 2 (S 1 refers to the circular cylinder boundary and S 2 is the boundary far from the body). In the same figure, U * is the incompressible inlet flow, D * defines the outer cylinder diameter, θ denotes upper angular position, x * and y * are the global coordinate system, and the symbol * means dimensional quantities. The flow to be studied is assumed turbulent and, thus, a LES modeling is employed to disconnect the large eddies from the smaller ones. As practical result, the continuity and the momentum equations can be filtered resulting in Equations (1) and (2), such as, respectively:

Mathematical Formulation
In Equations (1) and (2), the i-th filtered velocity component of the flow field is identified by u * i (being u * i = u * i + u * i with u * i representing the velocity fluctuation), ρ defines the density of fluid, p * represents the pressure filtered field (p * = p * + p * ; p * is related to the pressure fluctuation), υ defines the fluid kinematic viscosity, υ t defines the local eddy viscosity coefficient, and S * ij characterizes the deformation tensor of the filtered field [32].
The present methodology generates vortex blobs around the body surface during each time stepping aiming, thus, to solve the large eddy phenomena through the Equations (1) and (2). On the other hand, the small eddy phenomena are simulated using the concept of eddy viscosity coefficient (a typical LES approach). According to Lesieur and Métais [33], the local eddy viscosity coefficient at a point x * in the flow field, and at an instant t * , can be calculated via the local kinetic energy spectrum: being C k = 1.4 the constant of Kolmogorov and F * 2 (x * , ∆ + * , t * ) the local second-order velocity structure function of the filtered field; the latter is defined as follows [33]: In Equation (4), and for the three dimensional space, u * (x * , t * ) − u * (x * + r * , t * ) represents the average speed differences between the center of a sphere of radius |r * | = ∆ + * , located at x * , and points Energies 2020, 13, 6094 5 of 20 placed on the sphere surface (∆ + * separates the scales that are computed by the numerical method of the ones that should be modeled). In this approach the small scales are assumed to be homogeneous and isotropic, and the point of the flow field where one wants to quantify the turbulent activity is assumed as the center of the sphere [33].
The impermeability boundary condition on the cylinder surface (at S 1 in Figure 1) is established such as: u * n − v * n = 0 on S 1 (5) being u * n the fluid normal velocity and v * n the normal velocity of the body surface. The no-slip boundary condition on the cylinder surface (at S 1 in Figure 1) is imposed through the following equation: being u * τ the fluid tangential velocity and v * τ the tangential velocity of the body surface. The far away boundary condition is given by: It is very practical to make dimensionless all the quantities previously presented, aiming to obtain generality gain from the numerical simulations. Thus, in the above equations, D * and U * are assumed as length and velocity scales, respectively. As consequence, the non-dimensional time is defined by t = t * U * /D * .

The Lagrangian Vortex Method
In order to solve numerically the mathematical problem presented on Section 2, the first step is to discretize the circular cylinder surface shown in Figure 1. This is done through the so-called panel method [34], which consists on discretizing any solid boundary in panels, over which singularities are distributed in order to guarantee one of the boundary conditions presented in Section 2 at the "pivotal point" (center point) of each panel ( Figure 2). In the present approach, the cylinder surface is discretized and represented by NP flat panels with source distribution of constant density. The sources generation is necessary to ensure the impermeability boundary condition Equation (5) in association with the mass conservation of the problem.  The Lagrangian vortex method is nothing more than a boundary layer model, which consists on discretizing a flow property (in this case, the vorticity filtered field) using vortex blobs (or vortex particle elements), such that [35]: In Equation (8), ω (ω = ∇ × u) represents the vorticity filtered field; NV is the total number of Lamb discrete vortices necessary to simulate the viscous wake; and Γ k defines the strength of the k-th vortex blob placed at a pivotal point of a flat panel. The discrete vortices generation is necessary to ensure the no-slip boundary condition Equation (6) in association with the global circulation conservation of the problem. In addition, σ 0 defines the Lamb vortex core, which is measured such as [36,37] (see also Figure 3): In Equation (9), the time stepping, ∆t, is computed from an estimate of the velocity and advective length of the flow.
As previously described, the potential problem is solved using the panel method through constant-density source distributions [34]. On the other hand, the vorticity effect is included into the problem by taking the curl of Equation (2), which, in association with Equation (1), results in the well-known vorticity transport equation [38]. In two dimensions, the mentioned equation is scalar and the pressure term is eliminated, such as: Energies 2020, 13, 6094 7 of 20 being Re c a local Reynolds number including the local eddy viscosity coefficient, and defined as [25]: Equation (10) simulates the vorticity filtered field dynamic, being solved to each vortex blob during each time stepping. Figure 3 sketches the typical generation process of vortex blobs used for a smooth surface. Each vortex blob is positioned tangent to the pivotal point of the panel to ensure the no slip boundary condition of the problem Equation (6), as already commented.
The present roughness model modifies the generation process of vortex particle elements without change the body surface. It is known that the roughness effect on flows is to stimulate the turbulence development. Thus, the mentioned generation process is altered, in order to add an instantaneous inertial effect to each nascent vortex blob. This additional inertial effect works as an injection of momentum into the boundary layer and, as consequence, the flow can support the unfavorable pressure gradient present in the boundary layer for longer, in agreement with the physics involved in this kind of complex flows. Since roughness and turbulence are intimately related, the momentum injection imposed by the present roughness model is obtained from turbulence modeling developed by Lesieur and Métais [33], which was originally implemented in a two-dimensional Lagrangian manner by Alcântara Pereira et al. [25], and with improvements made by Bimbato et al. [37].
Thus, the generation process of vortex blobs is modified considering the turbulent activity around "panel's shedding point" (Figure 2), what is not taken into account for smooth surfaces ( Figure 3). The turbulent activity in the vicinity of the shedding point of the i-th panel is computed adopting a semicircle, with radius b = 2ε − σ 0 , and centered on i-th shedding point (ε = ε * /D * is the mean relative surface roughness). According to Bimbato et al. [14], a set of NR points is placed on the semicircle and the average speed differences, required to compute the second-order velocity structure function of the filtered field Equation (4), are calculated between the center of the semicircle (the shedding point) and the points on it ( Figure 4). The second-order velocity function is obtained through the following adapted expression [14]: In Equation (12), u t i (x i , t) is the total velocity on shedding point of the i−th panel (which is located at x i ) at time t, u t w (x i + b, t) is the total velocity at each one of the NR points, defined around the i−th shedding point (that is, on the semicircle of radius b = 2ε − σ 0 ) at time t, and (1 + ε) is idealized to simulate the kinetic energy gain because of the roughness interference. Figure 4 gives us an idea of the turbulent activity computation around the shedding points, which is the velocity fluctuation in the boundary layer.
It is noted from Equation (12) that, at each time stepping, there is a different velocity fluctuation around the shedding point of the i−th panel. Therefore, there is also different local eddy viscosity coefficients associated with these velocity fluctuations, Equations (3) and (4). Thus, the eddy viscosity coefficient at i−th shedding point, at time t, is computed such as: In Equation (13), σ 0 k is the vortex core of the k−th vortex blob placed at i−th shedding point Equation (9).
Since the eddy viscosity is associated with each time and with each panel shedding point, a local Reynolds number must be computed at each time stepping, as defined in Equation (11). Therefore, Re c i (t) represents the Reynolds number modified locally (at i−th shedding point) at time t, and this modification only is computed if the surface roughness effect cannot be neglected, that is, for υ t i (t) 0.
Equation (9) shows how the Reynolds number of the flow and the Lamb vortex core are related. Thus, if the surface roughness effect is present in the boundary layer flow (if υ t i (t) 0), the Reynolds number is modified locally Equation (11). The numerical effect is to change the core radius, σ 0 , of each nascent vortex blob by satisfying Equation (14) instead of Equation (9). Therefore, the generation process of discrete vortex elements, illustrated in Figure 3, is now modified to include an additional inertial effect to each nascent vortex blob (it can also be interpreted by an increase of strength ∆Γ in Figure 5). The new Lamb vortex core assumes the following form [14]: Energies 2020, 13, 6094 9 of 20 In Equation (14), σ 0c k (t) is the vortex core of k−th vortex blob, placed at i−th shedding point at time t, and modified by surface roughness effect (compare Figure 3 with Figure 5). The roughness model is based only on physics of turbulent flows.
The procedure above described to represent the body surface (sources generation) and the vorticity field (creation of nascent vortex blobs) is mathematically defined by two different kinds of linear algebraic equation systems; their solution is iterative and they are solved using the method of least squares, and can be easily linked with the roughness model, when activated.
The set of discrete vortex elements, positioned at shedding points, is necessary to simulate the vorticity filtered field dynamics according to vorticity Equation (10). With this purpose, Chorin [19] idealized an algorithm that splits the advective-diffusive operator of Equation (10) such as, respectively: Equation (15) implies that the vortex cloud is advected as a set of material particles of the flow in a typical Lagrangian manner.
The solution of the advection problem Equation (15) is given by integrating each vortex blob path equation, and using an explicit Euler scheme, in the following form: In Equation (17), u t k (x k , t) represents the velocity filtered field computed at position occupied by k−th vortex blob. Thus, the velocity field of the problem is composed by three contributions, i.e., (i) the incident flow (or the uniform flow) velocity, ui(x, t); the cylinder surface (using the source panels distribution), ub(x, t); the viscous wake (using the vortex cloud), uv(x, t).
For the incident flow, shown in Figure 1, its components take the form: The velocity induced by NP flat panels with constant-density source distribution at position occupied by k−th vortex blob is computed in the following form [34]: In Equation (20), σ i is the source strength per panel length and c n ki [x k (t) − x i ] represents the n−th component of the velocity induced at the k−th vortex blob by the i-th source flat panel.
The third contribution for the velocity filtered field computation is because of the vortex cloud (i.e., the vortex-vortex interaction, obtained through the Biot-Savart law [35]), which takes the following form: In Equation (20), Γ j is the strength of the j−th vortex blob and c n kj x k (t) − x j (t) represents the n−th component of the velocity induced at the k−th vortex blob by other vortex blob placed at x j . Therefore, the velocity filtered field,u t (x, t), is obtained by the summation of Equations (18)-(20) and its computation is used: (a) to obtain the position of vortex blobs at the next time stepping, x k (t + ∆t), by using the material derivative concept (Dω/Dt); (b) by roughness model to compute the average speed differences between shedding points and the points on the semicircles (Equation (12) and Figure 4); (c) by turbulence LES modeling (as will be described below).
Once advective equation Equation (15) has been solved, to properly simulate the vorticity field dynamics, it is necessary to solve, within the same time stepping, the diffusion equation (Equation (16)) to include molecular viscous effect. In the present numerical method, Equation (16) is solved through the random walk method [21], which was originally developed by Chorin [19]. The integral solution of the random walk method is [19]: being: Therefore, the diffusive displacement of the k−th vortex blob is computed, such as: In Equation (24), P and Q are random numbers, which may assume values in the range: 0 < P < 1 and 0 < Q < 1.
The presence of the local Reynolds number (Re c ) in Equation (23) indicates that turbulence manifestations must be computed into the diffusion process [25]. The essence of turbulence modeling used in this work consists on calculating the turbulent activity in points of the flow through the average speed differences existing around such points (Equation (4)). Instead of compute these differences between a center of a sphere and points on its surface, as established by Lesieur and Métais [33] for three-dimensional flows, Alcântara Pereira et al. [25] proposed to compute the average speed differences between a center of an annulus and points located between its internal and external radius. Both the center of the annulus and the points between its inner and outer radius are defined for each vortex blob presents in the computational domain simulating the vorticity field dynamics.
The adaptations proposed by Alcântara Pereira et al. [25] let to compute the second-order velocity structure function of the filtered field in the following form: In Equation (25), u t is the total velocity induced on points of interest ( u t = ui + ub + uv), N is the number of discrete vortices inside the annulus and r j defines the distance between the center of annulus (vortex blob under analysis, i.e., the k−th particle) and the points located between the internal and the external radius of annulus (the j−th particle placed inside it).
A statistical analysis performed by Bimbato et al. [37] concluded that proper internal and external radius of annulus is, respectively: r int = 0.1σ 0 k and r ext = 4.0σ 0 k (for smooth surfaces) or r int = 0.1σ 0c k and r ext = 4.0σ 0c k (for rough surfaces).
The use of Lagrangian vortex method is justified by the ease of work with high Reynolds number simulations (there is no numerical instabilities, which are associated with Eulerian techniques). Furthermore, there is no mesh in the Lagrangian technique and it is kind simple to implement the turbulence LES modeling using velocity differences instead of derivatives Equation (24), which must be closely related to the roughness model Equation (12).
It remains to calculate the aerodynamic coefficients (drag and lift forces). With this purpose, it is used the stagnation pressure definition, such as: In Equation (26), p * defines the static pressure and u * represents the velocity at any point into the fluid domain Ω (Figure 1). Then, the methodology to compute aerodynamic loads starts from a Poisson equation for the pressure, which is solved through the following integral formulation presented by Shintani and Akamatsu [39]: In Equation (27), the static pressure can be computed at i−th calculation point, being H = 1.0 valid for the fluid domain and H = 1/2 valid for solid boundary, Ξ represents a fundamental solution of Laplace equation, and e n defines the unit vector normal to each flat panel used to discretize the cylinder surface.
Using the static pressure obtained from Equation (27), the drag and lift coefficients are evaluated, such as, respectively: In Equations (28) and (29), p ∞ is the reference pressure far from the solid boundary, ∆S i is the length of the i−th panel, and β i is the angle of the i−th panel.
Finally, the algorithm implemented essentially consists of eight steps in the following sequence: (i) simultaneous generation of source panels and nascent vortex blobs (the roughness model is used, when activated); (ii) velocity field computation at each vortex blob; (iii) static pressure computation at each pivotal point followed by integrated aerodynamic loads computation; (iv) advection of the vortex cloud; (v) vorticity diffusion (the LES modeling is used); (vi) reflection of vortex blobs, which eventually migrate into the cylinder; (vii) velocity field computation at each pivotal point to restore the boundary conditions of the problem for new generation of source panels and nascent vortex blobs; (viii) advance by time ∆t. In the last years, the authors of this paper have made an effort to develop the in-house code by using FORTRAN programming language.

Results and Discussion
This Section presents the numerical results obtained for the problem of the two-dimensional, incompressible and unsteady flow past a circular cylinder. The body surface was represented by NP = 300 source flat panels. That chosen number for the panels ensures the convergence for the potential solution of the problem and, as consequence, a refined discretization for the vorticity field, which is reflected on the aerodynamic loads computation. In the case of flow around smooth or rough cylinders for the Reynolds number value studied here (Re = 1.0 × 10 5 ), the drag coefficient is computed only by form (pressure) drag contribution. According to Achenbach [3] measurements, the friction drag contribution does not exceed 2-3% to the total drag. Thus, the friction drag computation has been omitted in this work.
In the present approach, it has been used NR = 21 points on each semicircle defined around each panel shedding point in order to compute the turbulent activity around it ( Figure 4); this number is enough to obtain a reasonable value for the average speed differences Equation (12) in consonance with past study reported by Bimbato et al. [14].
In order to obtain good results, the time increment used is ∆t = 0.05 and an explicit Euler time-marching scheme is used Equation (17). All numerical simulations run until dimensionless time of t = 75, which is enough to reach a statistical equilibrium; the mean coefficient values (pressure distribution, drag and lift) are computed between 37.5 ≤ t ≤ 75.0. As previously presented by Oliveira et al. [18], the adopted criterion to compute aerodynamic loads has sufficiently been refined attaining the saturation state of the numerical simulations; this behavior is typical of a Lagrangian technique. Table 1 presents a summary of some typical experimental results published in the literature and the numerical ones obtained using the present methodology. There is a remarkable scatter among the experimental results, even those referring to the smooth circular cylinder. The scatter among experimental results is because of influencing parameters, especially by the low aspect ratio and high blockage. The cylinder tested by Achenbach [3] presented an aspect ratio of L * /D * = 3.3 (L * is the cylinder length) and the wind tunnel blockage was D * /B * = 0.17 (B * is the distance between tunnel walls). The cylinder tested by Achenbach and Heinecke [4] presented an aspect ratio of L * /D * = 6.75 and the wind tunnel blockage was D * /B * = 0.17, while Zhou et al. [12] used a configuration of cylinder with aspect ratio of L * /D * = 10. Thus, the comparison among experimental results and the present numerical results is not so fair, as previously discussed in Section 1. Unfortunately, there is a lack of numerical data of aerodynamic forces for the flow past a rough cylinder at upper subcritical Reynolds number of Re = 1.0 × 10 5 .  Table 1). The roughness model used here is able to predict supercritical flow features from a subcritical Reynolds number flow simulation. This conclusion comes from the comparison between Figure 6a,b.
The present numerical results show that, for small superficial roughness (0.00 < ε ≤ 0.10 %), the drag coefficient suffers a very small reduction (about 2.9%), which occurs in the subcritical regime (Figure 6a). With the progressive increase in surface roughness (0.20% ≤ ε ≤ 0.45%), the roughness model injects momentum on the boundary layer in a manner that flows with higher Reynolds numbers are simulated; it is remarkable the drag coefficient sudden drops until reach a minimum value (Figure 6b), that is a typical feature of critical flows. The mean drag reduction is about 16.2%. For even greater roughness (ε > 0.45 %), the drag coefficient grows up gradually (about 4.5%), which is the main feature of the supercritical flow regime (Figure 6a). Despite this, it is undeniable that the present numerical results have the same tendency as the experimental results of Achenbach [3] (Table 1). The roughness model used here is able to predict supercritical flow features from a subcritical Reynolds number flow simulation. This conclusion comes from the comparison between Figure 6a The present numerical results show that, for small superficial roughness ( % 0.10 ε 0.00 ≤ < ), the drag coefficient suffers a very small reduction (about 2.9%), which occurs in the subcritical regime (Figure 6a). With the progressive increase in surface roughness ( 0.45% ε 0.20% ≤ ≤ ), the roughness model injects momentum on the boundary layer in a manner that flows with higher Reynolds numbers are simulated; it is remarkable the drag coefficient sudden drops until reach a minimum value (Figure 6b), that is a typical feature of critical flows. The mean drag reduction is about 16.2%. For even greater roughness ( % 0.45 ε > ), the drag coefficient grows up gradually (about 4.5%), which is the main feature of the supercritical flow regime (Figure 6a). Although the percentage drop in the drag coefficient measured by Achenbach [3] is slightly more than double that calculated in the present work, most of the features of transitional flows are qualitatively captured by present methodology. Table 2 shows that boundary layer separation point agree well with the drag coefficient behavior pointed out in Table 1 and Figure 6b, what means the separation point moves downstream with drag reduction (critical flow regime) and upstream with drag gradual increase (supercritical flow regime).  Although the percentage drop in the drag coefficient measured by Achenbach [3] is slightly more than double that calculated in the present work, most of the features of transitional flows are qualitatively captured by present methodology. Table 2 shows that boundary layer separation point agree well with the drag coefficient behavior pointed out in Table 1 and Figure 6b, what means the separation point moves downstream with drag reduction (critical flow regime) and upstream with drag gradual increase (supercritical flow regime). Furthermore, Figure 7 indicates that the lift coefficient,C L , oscillates regularly because of the periodic vortex shedding mechanism [41], and the lift force variation of the smooth cylinder is greater than that of the rough cylinder, which agrees with experimental tests conducted by Zhou et al. [12]; the dash-dot lines in Figure 7 indicates the mean amplitude of lift force. In the present work, no analysis regarding the vortex shedding frequency is presented, however, a detailed discussion can be easily recuperated [14,17,18].
Energies 2020, 13, x FOR PEER REVIEW 15 of 24 Furthermore, Figure 7 indicates that the lift coefficient, L C , oscillates regularly because of the periodic vortex shedding mechanism [41], and the lift force variation of the smooth cylinder is greater than that of the rough cylinder, which agrees with experimental tests conducted by Zhou et al. [12]; the dash-dot lines in Figure 7 indicates the mean amplitude of lift force. In the present work, no analysis regarding the vortex shedding frequency is presented, however, a detailed discussion can be easily recuperated [14,17−18].     Figure 7. It can be noticed that the surface roughness effect causes irregular disturbances on instantaneous pressure distribution curves, which is reflected on time history of drag and lift forces. The mentioned disturbances occur because of inertial effect manifesting on the flow, which is provoked by roughness model. Thus, the boundary layer presents more momentum supporting the unfavorable pressure gradient, as expected, and, as consequence, the flow separation is delayed. This was shown in the joint analysis of Figure 6b and Table 2; now it can also be verified with the analysis of Figure 9. In Figure 9a, the non-zero lift force for the smooth cylinder (C L = 0.021) and for the less roughened cylinders (C L = 0.035 for ε = 0.10% and C L = 0.005 for ε = 0.20%) can be attributed to numerical rounding errors (panels method and vortex cloud contributions). It can be noted that the stagnation point of these three tests is located at θ stag = 0.6 • ; this happens because the discretization of the cylinder using the panel method starts with θ = 0 • (Figure 1) and the panel pivotal point is placed at θ = 0.6 • . It is important to emphasize that all the computations are done on pivotal points ( Figure 2). So, this is the reason why θ stag 0 for the three cases mentioned.
On the other hand, it is notorious not only the non-zero lift force for the rougher cylinders but its direction change (C L = +0.100 and θ stag = −1.8 • for ε = 0.45%; C L = −0.132 and θ stag = +1.8 • for ε = 0.70%), see Figure 9a,b. The boundary layer transition is the only explanation of this occurrence on a symmetrical bluff body, like the circular cylinder. The involved physics indicates that first an asymmetric flow is caused because the boundary layer becomes turbulent at one side of the cylinder; furthermore, a separation bubble and a non-zero lift force can be identified. On the other side of the cylinder is generated another separation bubble, when the boundary layer becomes turbulent on that side. According to Kamiya et al. [6], a fully symmetric flow is not reach as soon as the second separation bubble is formed. They measured C L = 0.2 and θ stag = +2 • when the first separation bubble was formed, followed by C L = 0.0 and θ stag = 0 • at higher Reynolds number (Re ≥ 7.0 × 10 5 ). The critical flow regime is so unstable, and the lift force changing direction was also captured by Kamiya et al. [6] by using a smooth cylinder.
One more feature of the complex transitional flow is captured by the present methodology. According to Bearman [5], the drag coefficient is approximately equal to the measured base pressure coefficient for three representative Reynolds numbers, i.e., Re = 2.0 × 10 5 (no separation bubble), Re = 3.7 × 10 5 (one separation bubble) and Re = 4.0 × 10 5 (two separation bubbles). Figure 10 shows the comparison between drag and base pressure coefficients, calculated by the present Lagrangian approach. It is remarkable that the present methodology agrees with the observations given by Bearman [5] (using the L * /D * = 12 and D * /B * = 0.06 configuration) for the critical flow regime. As expected, there is a disagreement between the computed values of base pressure and drag coefficients for the supercritical flow regime, because of increasing in drag coefficient.
In addition, Figure 11 presents the time-averaged pressure distribution for the flow around a cylinder with different surface roughness influences, where the experimental result was reported by Blevins [42] for smooth cylinder. It can be observed only small differences in the base pressure between the less rough cylinder and the smooth cylinder. As can be identified, the numerical result for the rough cylinder (ε = 0.45%) clearly shows a higher increase in the base pressure, which physically agrees with the higher drag reduction (Figure 6b). Finally, Figure 12 illustrates the viscous wake developed downstream the cylinder in the end of numerical simulations (at t = 75). The blue points in the same figure represent instantaneous vortex blob distributions; they indicate that the vortical structures for ε = 0.45% are narrower than the other ones, which is in accordance to what is expected from critical flow regimes. That conclusion is supported by upper and lower separation points behavior for ε = 0.45%, which are 93 • (Table 2) and 85.8 • , respectively. Finally, Figure 12 illustrates the viscous wake developed downstream the cylinder in the end of numerical simulations (at t = 75). The blue points in the same figure represent instantaneous vortex blob distributions; they indicate that the vortical structures for % 0.45 ε = are narrower than the other ones,

Conclusions
This paper describes a two-dimensional Lagrangian vortex method blended with both LES and roughness models to capture changes in a bluff body aerodynamics. It can be observed that the roughness model employed in the present work is classified as a passive method of vortex shedding control. The passive control method implies to modify the body surface, which is reflected on its aerodynamic loads behaviour. In fact, the roughness model presented here acts as a passive method modifying the boundary layer flow by an injection of momentum into it (which occurs on rough surfaces), but without changing the discretization of the body surface.
In the present approach, the vorticity field is discretized and represented by a cluster of vortex blobs, which are generated during each time stepping from the body surface. The potential theory is adapted to compute the circular cylinder surface, which is discretized using flat panels with constant-density source distribution over them. The aerodynamic loads are calculated through an integral formulation derived from a Poisson equation for the pressure, where the instantaneous vorticity field contributes to the computation. The chosen number of flat panels ensures the convergence of the solution for the potential velocity field and, as consequence, establishes the number of nascent vortex blobs. As the non-dimensional time runs, the number of vortex blobs increases and thus satisfies the vorticity transport equation forming the von-Kármán vortex street. In an upper subcritical Reynolds flow, the drag coefficient amplitudes are significantly lesser then the lift ones (Figure 7a). This behaviour indicates that the numerical simulations have captured the vortex shedding mechanism properly; it is important to observe that the saturation state of the

Conclusions
This paper describes a two-dimensional Lagrangian vortex method blended with both LES and roughness models to capture changes in a bluff body aerodynamics. It can be observed that the roughness model employed in the present work is classified as a passive method of vortex shedding control. The passive control method implies to modify the body surface, which is reflected on its aerodynamic loads behaviour. In fact, the roughness model presented here acts as a passive method modifying the boundary layer flow by an injection of momentum into it (which occurs on rough surfaces), but without changing the discretization of the body surface.
In the present approach, the vorticity field is discretized and represented by a cluster of vortex blobs, which are generated during each time stepping from the body surface. The potential theory is adapted to compute the circular cylinder surface, which is discretized using flat panels with constant-density source distribution over them. The aerodynamic loads are calculated through an integral formulation derived from a Poisson equation for the pressure, where the instantaneous vorticity field contributes to the computation. The chosen number of flat panels ensures the convergence of the solution for the potential velocity field and, as consequence, establishes the number of nascent vortex blobs. As the non-dimensional time runs, the number of vortex blobs increases and thus satisfies the vorticity transport equation forming the von-Kármán vortex street. In an upper subcritical Reynolds flow, the drag coefficient amplitudes are significantly lesser then the lift ones (Figure 7a). This behaviour indicates that the numerical simulations have captured the vortex shedding mechanism properly; it is important to observe that the saturation state of the present numerical simulations was attained and it is supported by previous discussions given by Oliveira et al. [18], since the numerical approach is the same.
In regarding of the rough circular cylinder, the numerical results show that the present methodology is able to capture, even using two-dimensional simulations, many attributes of this complex flow such as: (i) the drag crisis followed by a gradual increase on the drag force; (ii) the non-zero lift force, since according to Zdravkovich [10] the low and high pressure on two sides of the cylinder interchange in different runs, because the side at which the separation is turbulent switches from one side to another occasionally; and (iii) the occurrence of asymmetric separation bubbles generation. The asymmetries are part of the physics involved in critical Reynolds number flows around a circular cylinder; that behavior was previously described by Zdravkovich [10], whom dedicated all his research to investigate different aspects of flows around cylinders.
Although the lift force behavior is an important feature of the critical/supercritical flow regimes, studies that focus on this force are scarce, which valorizes the present numerical results. The results for both predictions of boundary layer separation and stagnation point displacement also help us to report the sensitivity of the roughness model. Furthermore, each experimental study available in the literature uses a different technique to represent the surface roughness effect (wires, spheres, small holes, among others), and each technique is subject to different influencing parameters (free stream turbulence, aspect ratio, blockage and surface texture), which produces a great scatter in the measurements. As consequence, in most of situations it is not appropriate to compare two different experimental results or even an experimental and a numerical one.
Finally, the authors also intend to adapt the present vortex code to study other practical engineering situations, such as problems involving wake interference [43], ground effect mechanisms [17,18,27,28], vortex-induced vibrations [13,26] and thermal effect [44,45]. The present paper also aims to lay the theoretical foundation of the present numerical method for further algorithmic extensions to simulations of such fluid flows in three dimensions.