Application and optimization of drag reduction characteristics on the flow around a partial grooved cylinder by using the response surface method

ABSTRACT It is well known that drag reduction occurs when the flow is passing by a grooved circular cylinder at certain Reynolds numbers, which has been used as a powerful energy saving method in a broad range of circumstances. However, a challenge here is how to evaluate the combined effects of depth, width and location of a given triangular groove set covering half of the cylindrical surface area. A useful approach to quantitatively analyze the influence of these different factors on drag reduction using the response surface methodology is described here. The flow characteristics, including drag coefficient, flow velocity, turbulent kinetic energy and vorticity, were calculated by numerical simulation. The results showed a great drag reduction effect under the legitimate set of groove structure parameters at a super-critical Reynolds number providing a base for optimization process in various engineering applications. The drag reduction mechanism found from this research could extend to other cases and should provide insights into engineering applications like car grilles.


Introduction
With the increase in energy consumption, there is a growing pressure to improve aerodynamic performance in order to achieve the energy saving for different applications, especially through reducing the drag force, a major component in energy consumption. For example, the overall drag costs a significant portion of energy consumption in transportation systems, theoretical as well as experimental researches have been conducted to reduce the vehicle aerodynamic drag. Since the radiator grille caused drag force affects the aerodynamic outcome of a vehicle including the cooling system, it has been considered as a major target for drag reduction as well as efficient fuel utilization. In some cases, the drag reduction can be achieved up to 7% by changing the model of radiator grille (Jama, Watkins, & Dixon, 2006). Therefore, it is attractive to modify the structure of grille leading to the optimal gas flow around the grid that will give rise to the reduction in aerodynamic force.
For some vehicles, the grille is designed in a cylinderlike shape, and therefore, the principles of the flow around the cylinder can be applied to grille design for the goal of reducing the aerodynamic resistance coefficient, Moreover, it can also be used in many other applications such as the design of the natural gas pipeline surface (Luo & Zhang, 2012). Previous studies revealed CONTACT Xiaowen Song songxw@zju.edu.cn that cylindrical surface covered with a non-smooth structure made significant changes in the flow characteristics at the high Reynolds number (Achenbach, 1971). When the boundary layers separate from the sides of a cylinder and form the free shear layers, the Karman vortex street generated behind the cylinder is unstable as the Reynold number of the flow over cylinders correspond to the sub-critical Reynolds number (Doolan, 2010). However, the non-smooth cover on the surface of the cylinder will reduce the coefficient of drag in great efficiency as it causes the critical regime at a lower Reynolds number (Rodriguez et al., 2017). This phenomenon has been verified by different experiments, which tested different non-smooth structures on the cylindrical surface of certain applications such as dimples (Bearman & Harvey, 1993;Tan, Koh, & Ng, 2016), grooves (Yamagishi, Kimura, & Oki, 2013;Yokoi, Igarashi, & Hirao, 2011). Numerous researches have drawn the conclusion, that non-smooth structure can delay the full separation of the boundary layer and optimize the flow condition behind the cylinder as the drag force is directly related to the Karman vortex street and the full separation point of the airflow (Yamagishi & Oki, 2007). In addition, the same non-smooth surface with different designed parameters have also been examined as an important aspect in drag reduction, for example, the influence of single rectangular groove on the flow characteristics around circular cylinder have been reported (Canpolat, 2015;Canpolat, 2017). U-grooved (Alonzo-García et al., 2013), Vgrooved (Alonzo-García et al., 2014) or cactus-inspired grooved cylinder (El-Makdah & Oweis, 2013;Karaki, Abboud, Daher, Osman, & Oweis, 2008) have been designed and tested in research, and the effects of groove number (Yamagishi & Oki, 2005), and sharpness (Yamagishi & Oki, 2007) were also investigated with respect to their influences on the aerodynamic performance. Therefore, the drag reduction characteristics of flow around the cylinder covered with non-smooth surface mimic the condition of a moving automobile at highspeed. And the theoretical base of the engineering application in this situation can be established. Besides, previous studies focused on exploring the drag reduction impacts by individual parameters of non-smooth structure (Aoki, Lee, & Oki, 1998;Yamagishi & Oki, 2005, 2007. In order to achieve a better design method for applications, the effects of various structural parameters on drag reduction were studied in considering these parameters collectively rather than individually. The interactive influences were also explored. By using the response surface methodology, studying the interactive relationship among various factors in a given situation was achieved in this research. Meanwhile, through the establishment of the hierarchy of the effects among these factors, a standard procedure concerning the nonsmooth design for engineering application has been set up. To validate the accuracy of the flow field obtained by numerical simulation, a particle image velocimetry (PIV) test was conducted.

Problem description
Increasing flow characteristics based researches demonstrated that the drag reduction is caused by the delay of full separation point due to the turbulence transition in the boundary layers, which suggests that the groove placements in front of the separation point play a major role in the occurrence of drag reduction. However, in most of the studies, the non-smooth structures were arranged as a whole unit covering the entire surface of the cylinder without considering the regional influence (Yamagishi & Oki, 2004, 2005, 2007. Considering the previous analyses were based on full coverage of the nonsmooth structures on the cylinder surface, alternative experiments should be designed to explore whether the grooves on the back half of cylinder will help to reduce the drag force, and to analyze the differences in flow characteristic between the semi-grooved cylinder and the fully grooved one.

Model description
To reduce the challenges in the manufacturing process while still meeting the requirements for the engineering application of interest, we choose the drag reduction model based on triangular groove structure covering the surface of the cylinder. The analysis cylinder was 0.025 m in diameter and 0.15 m in length resembling the size of most car grilles. The size ratio between groove and cylindrical was setup based on the analysis model in Yamagishi's experiment which gave better drag reduction, and two models with different non-smooth surface coverages were established: a fully grooved model, Model A, in which a 0.025 m diameter cylinder was attached with the grooves that were 1.875 millimeter in width and 0.26 millimeter in depth following the ratio between groove and cylinder in the model described by Yamagishi' s experiment (Yamagishi & Oki, 2005). A total number of 32 grooves were arranged with an angular interval of 11.25°. A semi-grooved model, Model B, was also established, in which half of its exterior surface area was covered by the grooves with the same size as model A (Figure 1).

Calculation domain
The analysis was made in the unsteady three-dimensional turbulent flow with following calculation domain set: the distance from the center of the circle to the upstream boundary surface was 6d (Note: d represented the diameter of the cylinder), whereas the distance to the downstream boundary was 16d. The wall boundary surface on both sides was 6d each away from the center of the circle. The size of the calculation domain is slightly larger  than what is reported by Yamagishi and Oki (2005) to ensure the solution is independent of domain extents. The three-dimensional calculation domain is shown in Figure 2.

The boundary conditions for numerical simulation
When a fluid flows pass a cylinder, the increase in Reynolds number leads to the sudden drop of drag coefficient after sub-critical region, followed by a slight rebound in the super-critical region, resulting a constant value in the trans-critical region (Zdravkovich, 1994). The value of C d will be small in the super-critical region, where the Reynolds number for the flow passing the cylinder with 32 grooves is ranging from about 4 × 10 4 to 6 × 10 4 (Yamagishi & Oki, 2005). The Reynolds number of the model was designed within this range to reduce the drag force. Meanwhile, for studying the effect of grooved surface on drag reduction in high-speed wind, the inlet velocity was set at 35 m/s (78 mph or 126 km/h) modeling the car speed on highway giving the Reynolds number Re d = 5.83 × 10 4 (Re d = (ρUd/μ), ρ: air density, ρ = 1.2 kg/m 3 , μ : dynamic viscosity of air, μ = 1.8 × 10 −5 Pa s), which is the exact conditions of the supercritical region. The outlet boundary was defined under the assumption that the flow is fully developed. One side boundary was set up as the symmetry wall, and others were established as stationary walls (slip wall).
Inlet boundary conditions: a uniform incoming fluid velocity is set at the entrance of the computation region. The speed components in three vector directions were set as follows: Outlet boundary conditions: The outlet boundary was set based on the complete development of the flow field at the back of the calculation domain. The pressure and velocity were set as follows: The establishment of physical model The flow is unsteady and incompressible. The Realizable k-model (Shih, Liou, Shabbir, Yang, & Zhu, 1995) was used in combination with the standard wall function near the cylinder surface as this turbulent model gives a superior performance in predicting the turbulent boundary layer separation at high Reynolds number (Mulvany, Tu, Chen, & Anderson, 2004). And the numerical results calculated by the realizable k-model are close to the experimental results under the same Reynolds number (Chen, 2013). The conservation equations of mass and momentum were also used to control the main flow parameters.
The mass equation is shown below: Momentum equation is as follows: ρ: air density; τ ij : shear stress, in which the i is the normal direction of the surface, whereas the j is the direction of the force; F x : the body force in the direction of x; F y : the body force in the direction of y; and F z : the body force in the direction of z. (The body forces in three directions were set as zeros since the gravity is insignificant.)

The establishment of grid models
The grid model for the rectangular cuboid was generated using unstructured grids by STAR CCM + 12.02, which simulates wake of cylinder accurately (Vermeire, Witherden, & Vincent, 2017). Also, it is reliable in predicting the boundary layer transition (Lin, Raghunathan, Raghunathan, & Mcilwain, 2012). To verify the independence of the computational grid, the drag coefficients of the fully grooved model (Model A) were calculated in STAR CCM + 12.02 under five grids with multiple densities. The calculation results are shown in Table 1. (The drag coefficient C d is defined as (2F d /ρu 2 A), in which the F d represents the drag force that is by definition the force component in the direction of the flow velocity, the ρ stand for the mass density of the fluid, u means the flow speed of the object relative to the fluid, and A is the reference area.) The calculation results proved that when the number of elements in mesh reached 989867, the calculated values of the two parameters remain approximately the same: the difference between C d values calculated by models with 989867 and 1211330 elements is 0.362%. Therefore, the flow solution showed almost no change as the computational mesh was refined in this case. Thus the mesh with 989867 elements was chosen for this study. The total thickness of the boundary layer was 0.8 mm, and the average first-cell y+ value was 33.724, close to the 30 used in a turbulent wall-model. The whole meshes of the calculation domain are shown in Figure 3, and the shapes of grooves and meshes adjacent to the circular cylinder surface are shown in Figure 4.

The verification of the hypothesis
The grid models were incorporated into the versatile fluid analysis software package STAR CCM + 12.02, where the flow characteristics were calculated by finite volume method. Before that, to make sure that the solution does not change within refined time interval with a tolerance less than 5%, the independence of time interval was verified by calculating the drag coefficient of Model A under five different time interval, the computational result is shown in Table 2, the time interval was presented in a non-dimensional form T. (T = (u × t/d), u: the velocity of fluid, t: the time step, d: the diameter of cylinder).
As the value of the drag coefficient did not change beyond the tolerance when the non-dimensional time  interval refined to 0.112, this size was chosen to define the time interval. The measurement of the smooth cylinder was also conducted in parallel for comparison, the simulation result of drag coefficient is 1.21 when the value was averaged over 0.1s, which agrees with the measurement of Wieselsberger (1921) at Reynolds number Re d = 5.83 × 10 4 and suggests the reliability of this drag measurement.
The time series of the instantaneous drag (C d ) and lift (C l ) coefficients for models A and B are presented in Figure 5, the numerical simulation confirms that cylinder surface with semi-grooved attachments oriented toward the incident flow showed more significant effect on drag reduction compared with the cylinder fully covered with grooves, as the drag coefficient for model B is about 0.785, smaller than 0.83 obtained from Model A at Re d = 5.83 × 10 4 . And the amplitude value of the lift coefficient for model B is 0.33, also smaller than the 0.37 for Model A when the amplitude value of both cases tend to be stable, showing the positional effects of the groove on drag reduction.
which is by definition the force component in the direction of the flow velocity; ρ is the mass density of the fluid; u is the flow speed of the object relative to the fluid; A is the reference area. C l = (2L/ρu 2 A), where: L is the lift force; ρ is the mass density of the fluid; u is the flow speed of the object relative to the fluid; A is the reference area.)

Optimization procedure
In our efforts to optimize the model, it is obviously time consuming to perform a great deal of CFD simulations, one way to overcome that is to simplify the grid model and use some advanced turbulence model like modified Spalart-Allmaras model (Pasha, Juhany, & Khalid, 2018) to predict the boundary-layer interactions. Also, it can be solved by an alternative better performed method than the calculation of the CFD tools Fotovatikhah et al., 2018). Some of the advanced algorithms like PSO (IPSO) algorithm and neural network are used for optimization (Zhang, Zhang, Tezdogan, Xu, & Lai, 2018).

Response surface method
The response surface is a method that uses a series of definitive experiments to fit the response surface for simulating the real limit state surface. The basic idea is to set an analytic expression between the limit state function and the basic variable. The advantage of the response surface method is that every level of the experimental parameters can be continuously analyzed during the optimization for the test conditions. The effects of different factors can be studied simultaneously and the research results can be demonstrated by graphics intuitively. The response surface method can obtain an optimal result within a certain range more efficiently by establishing the multivariate quadratic regression model.

Selection of optimal parameters
The comparison experiment between the full-and semigrooved cylinder shows a negative effect on drag reduction by the grooves oriented in the direction against the incident flow. In order to conduct a comprehensive study to evaluate how the groove set near the separation point affects the drag reduction, the layout of the grooves on the cylinder surface was changed by adjusting the angle α between the rearmost groove and the forward stagnation point from 90°to 100°when the cylinder surface was partially covered by grooves ( Figure 6).
And the size of groove, including the width and the depth, is also thought to give a great influence on the characteristics of the flow near the cylinder surface. Thus the three levels' response surface method (RSM) test, by which the influences of different factors can be considered simultaneously for engineering application, was designed to find the influence of multiple factors on the drag reduction. At the same time, the combined effects of different factors on drag reduction can also be studied in this process for future application. According to the aforementioned predication based on the simulation results, the depth, the width, and the angle of the rearmost groove on both side of the cylinder surface were chosen as the independent variables. For an affordable manufacture process, the width and the depth of grooves were changed on the basis of model A and met the requirements of the cylinder size as used in engineering application, since it was proved to give a better effect on drag reduction (Yamagishi & Oki, 2007). The depth and width are fixed for all grooves in a model. The angle α of the rearmost groove was changed from 90°to 100°i n order to discuss the influence of groove in this position, which was located near the separation point with its orientation against the incident flow.
The level-factors and the test plan are presented in Table 3:

Results of the optimization
In order to intuitively establish the corresponding relationship between the studied factors and the responses, as well as to improve test efficiency, the Box-Behnken (Box & Behnken, 1960) method was used to establish the response surface model, in which the final optimization results are evaluated by the reduction of the drag coefficient. Also, the CFD simulations were performed and used as inputs for the optimization procedure for the relative accuracy of the simulation results.
Due to the design requirements of the Box-Behnken method, the experiments with a center value of the three parameters were repeated five times. A series of tests were conducted according to the test plan set up by STAR CCM + 12.02. The test plan and the simulation results of the drag coefficient for every test are shown in Table 4: Compared with the fully grooved model A, the drag coefficient value in test 6 under the partially grooved model is reduced to 0.719. The drag reduction ratio was calculated using the following equation: where the C d1 is drag coefficient of model A and the C d2 is drag coefficient of an optimized model. For the model used in test 6, the value of its R D reaches to 13.37%, which positively correlates with the observed drag reduction.
The test results demonstrated that, by adjusting the sizes and placements of the rearmost grooves, the capability of grooves on cylinder surface drag reduction will be further improved.

Analysis of the simulation results
In order to study the effects of individual factors on the drag reduction performance and explore the crosseffects among them, the analysis software Design-Expert (Montgomery, 2004) was used to establish the response surfaces to draw a limit-state surface and to deal with the conversion between each parameter (input) and the C d value (output). The significance of linear and quadric terms consist of different parameters, including the A, B, C, AB, AC, BC, A 2 , B 2 and C 2 (A = depth, B = width, C = angle α of the rearmost groove) were analyzed by the variance module. To ensure the significant characteristic of each independent variable in the regression equation, the model was optimized based on the initial test result of each term to guarantee the statistical significance of the regression equation.
The results of the variance analysis calculated by ANOVA for the response surface quadratic model (Montgomery, 2004) are shown in Table 5 without insignificant terms.
Based on the results of variance analysis, the model is sure of statistical significance as the p-value is less than 0.05. It means that the model has a high reference value on research. Meanwhile, the significant analysis of independent variables shows the significance of linear terms A, B, C, and quadratic terms AB, A 2 , B 2 , C 2 (p-value < 0.05).
The significance of the quadratic term AB verifies the reciprocal effects between these two factors. It explains why the effect of one independent variable results inconsistent changes of another one at a different level, which should be considered in the process of design optimization for engineering application and cost saving.
The p-value of fit-lack is defined as an x 2 -test to measure how well the model fits the data and whether the model error is equal to the assumed value. The pvalue of fit-lack reaches 0.1479 and proves the reliability of the model containing all the necessary items. The null hypothesis states that the model error mean square is equal to the hypothesized value (pure error), versus the alternative that it is greater than the hypothesized value. When the test p-value is greater than 0.05, the null hypothesis can be accepted in the form without fit. Therefore, the regression equation can be used to replace the real value of reality for analyzing. The equation of quadratic response surface regression was built as follows: (A = depth, B = width and C = angle α of the rearmost groove) The diagnosis of the distribution relationship between the residuals, the predicted value and the computational value are shown in Figure 7: The normal plot of residuals shows a linear correlation, which means there is a normal distribution difference between the actual observations and the regression estimated by regression analysis. The corresponding relationship shows a non-regular distribution between the predicted value and the residual value, but a linear distribution between the computational value and the predicted value. This regularity in distribution represents the high confidence of the predicted value calculated by the regression model.  The contour maps of the influence of the factors on the drag coefficient were used to demonstrate the degree of influence, as well as to express the interactive effect among the factors. The results are presented in Figure 8.
The density of the contour line is a standard used here to evaluate the influential levels of different factors on drag coefficient, namely, when the contours in one certain direction show a tendency to be more denser than the other during the change of independent variables, it means that the alternation of the parameter represented by this direction can lead to a larger change of the drag coefficient translating to a higher impact ability of this factor. Through the analysis of the three contour plots, it can be concluded that depth has been shown to have the greatest impact on the solution, followed by the effect of angle, then width afterwards, which provide a practical reference in the application.
From the three-dimensional response surfaces shown in Figure 9, it can be concluded that the influence of the width and the angle α of the rearmost groove on the drag coefficient are quadratic, which means there is a minimum drag coefficient within the interested range of these two parameters. And the relationship between factor A (depth) and the drag coefficient shows a linear correlation, indicating that the drag force is reduced when the value of depth decreases within the range allowed by grilles size. Therefore, a better groove structure can be designed to reduce the drag coefficient under this condition.

Selection of the optimal model
Through the analysis of the regression equation, an optimal model can be built up to have all these three factors with significant effects on drag reduction. The optimal solutions obtained by the numerical optimization module in the software Design Expert is shown below ( Table 6):  Five different optimized solutions were obtained by the regression equation, and an optimal solution that shows the highest perspective was selected for numerical simulation and experimental validation. To explore the mechanism of the drag reduction by this optimized groove structure, the simulation model was built up with a groove structure of 0.13 mm in depth and 1.886 mm in width. The angle α between the rearmost groove and the forwarding stagnation point was 92.904°, the diameter of the cylinder, the size of the calculation domain, and the settings of the boundary conditions were the same as the simulation of Model A to ensure the reflection of numerical simulation. Figure 10 presents the numerical results of the drag coefficients C d of the three models.
Based on the figure, the value of the drag coefficient finally averages around 0.695, with an error of 0.967%. Compared with the predicted result in the optimal solution designed by regression equation, the numerical simulation is considered to be reliable as the error is less than 5%, which is also verified for the reliability of the response surface. The optimal model shows 12% larger drag reduction compared to model B, and this value is larger than any test in the response surface method. Therefore, it can be concluded that the response surface method can be reliably used for designing the parameters of the non-smooth structure within the established parameter range to set up the optimal model for engineering application. Meanwhile, the drag coefficient value of the optimal model based on the aforementioned solution is smaller than other models with different curved sectional grooves reported in the previously study (Yamagishi & Oki, 2007).

Results and discussion
For achieving a better result in drag reduction when flow passes a cylinder, and expanding the useful field of non-smooth surface in engineering applications, the mechanism of drag reduction was studied through the measurable results of numerical simulation.
While moving around a cylinder, the flow will complete the transition from the laminar flow to the turbulent flow on the cylinder surface, and turns irregular. In some real cases of applications, positive pressure on object surface generally varied greatly, among which a narrow windward tends to suffer higher wind pressures, while a wide one is associated with severer negative wind effects (Mou, He, Zhao, & Chau, 2017). In this case, a large region of inverse pressure will be formed within the Karman vortices behind the cylinder. With the pressure difference between the front and rear area of the cylinder, the total drag force increases significantly, and the pressure drag becomes the main component of it, contributing more than 98% proportion wisely, whereas the skin friction (or viscous) component is responsible for the remaining 1-2% (Achenbach, 1971). Based on this mechanism, the roughness induced drag reduction on the cylinder surface could be interpreted by the convergence of the difference in pressure. This forms a specific aspect in the current research.
It is thought that the roughness on cylinder surface will cause partial separation or reattachment of the flow nearby the edge of the grooves. The flow situation will be changed by these repeated cycle of partial separationattachment, and a closed circulating flow will be formed within the cavity of the groove, leading to the turbulence of the flow in the boundary layer, especially near the cylinder surface. As the energy containing parts of the flow plays an important role in controlling vortices formation in the wake region (Apacoglu, Paksoy, & Aradag, 2011), the disturbance of flow accompanied by the momentum transfer causes higher kinetic energy of the flow at the bottom of the boundary layer. Thus the area of the rear adverse pressure can be changed, and the flow is able to overcome more wall shear force on the cylinder surface and to shift the full separation point towards the downstream side. In this way, the drag force reduction is achieved when the pressure drag becomes smaller.
However, the above function of drag reduction is only applicable by the grooves arranged at the position upstream of the full separation point. When the nonsmooth structure located near the full separation point, the turbulence of flow at this position may speed up the vortex shedding as the flow separating off at the edge of the groove would not be able to reattach to the cylinder surface anymore. Also, the grooves set at the rear surface of the cylinder may not help to reduce the drag force, but increase the friction force. Furthermore, the larger turbulence as well as the higher Reynolds number of flow caused by grooves at this area may lead to the lower pressure, all of which have negative effects on drag reduction.
According to the aforementioned flow characteristics influenced by grooves at different regions of the cylinder surface which is responsible for the main causes of the drag reduction, the parameters of flow near the cylinder surface, the turbulence behind the cylinder, and the vorticity distribution were all discussed. The fully grooved cylinder model A, the semi-grooved cylinder model B and the optimal model C designed by the response surface method were used for analysis when the Reynolds number Re d = 5.83 × 10 4 . By contrasting the different flow parameters of these three models obtained from the numerical simulation, the mechanism of drag reduction when the cylinder was partially grooved with optimal structure was studied deeply. The flow is orientated from left to right in all figures.

Turbulent distribution
The distributions of the turbulent kinetic energy near Model A, Model B and optimal model surfaces with Re d = 5.810 4 are shown in Figure 11. By comparing the distributions of turbulent kinetic energy in model A and model B, it is clear that the groove structure located near the full separation point causes larger flow turbulence at the corner, which can be seen in the figure, and the partial separation at this place may trigger the final separation of the flow early without reattachment anymore. The region of high turbulent kinetic energy is more concentrated and closer to the cylinder at the downstream side in model B. That demonstrates that the grooves have a negative effect on the delay of boundary layer separation if they are too close to the full separation point.
The different distributions of turbulence between model B and optimal model show the advantage of the optimized grooves in drag reduction: the optimal model has a larger turbulent kinetic energy near the cylinder surface compared with Model B. The region with large turbulent kinetic energy adheres to the surface more broadly in the optimal model, leading to the exchanges of the momentum within boundary layer closer to the cylinder surface. Thus the flow at this place has higher kinetic energy to overcome more viscous force as well as adverse pressure, resulting the full separation point shifting towards the downstream side. Figure 12 presents the distributions of turbulent kinetic energy near the surface of cylinders at angle θ = 100°from the stagnation point, revealing the characteristic of this flow parameter near the cylinder surfaces at downstream side. The distribution of turbulent kinetic energy near the cylinder surface at rearward position shows a regular pattern that, the position of maximum turbulent energy in the boundary layer of the three models is located at different heights above the cylinder surface. As the aerodynamic drag decreases, it corresponds to the phenomenon that the position of maximum turbulent kinetic energy is closer to the cylinder surface and the full separation point is much delayed. This simulation results can well explain the aforementioned theory: The large turbulent flow of model A is farther away from the cylinder surface at rearward location compared with model B. Meanwhile, less exchange of flow kinetic energy happens near the cylinder surface in model A. These are adverse effects on the delay of the full separation point. Meanwhile, in model B, the value of turbulent kinetic energy is larger at this position, which corresponds to the visual result of turbulent distribution in Figure 11.
Similarly, because of a smaller geometric structure in the optimal model, the reattachment of flow in optimized groove is more inside compared with that in model B. The position of maximum turbulent kinetic energy for optimal model is about 0.5z/d when it is about 0.5025z/d in model B (z is the height from center of cylinder in radial direction, d is the diameter of cylinder). Thus the optimal model seems to keep the transfer of the momentum closer to the cylinder surface.

Velocity distribution
The computational results of the velocity distributions near the surfaces of model A, model B and optimal model with Re d = 5.83 × 10 4 are shown in Figure 13.
Because of the gas viscosity, velocity gradient of flow exists in the boundary layer. The disturbance of the flow causes the exchange of energy between the high-speed flow and low-speed flow, which increases the velocity of flow at the bottom of the boundary layer. The simulation results clearly verify that: as the flow at the bottom of boundary layer in model B contains more turbulent kinetic energy, the region with high velocity becomes larger and gets closer to the cylinder surface compared with that in model A, which confers the flow near the cylinder surface more kinetic energy to adhere to the surface and to delay the full separation of the boundary layer in the area of reverse pressure. Meanwhile, this characteristic also leads to a consequence that the large velocity of flow is more close to the cylinder compared with model A, and decreases the width of the wake. The same velocity distribution characteristics also appear in the comparison of the visual results between model B and optimal model. To verify the above regularity, the velocity distributions in the boundary layer at angle θ = 100°from the stagnation point near the cylinder surface of models A, B and the optimal model were shown in Figure 14.
From the numerical results of model A and model B, it can be found that at θ = 100°, the value of velocity at the rear positions of the semi-grooved cylinder becomes larger compared with the surface of the full-grooved cylinder. The position where the velocity reaches the maximum value shifts to a low z/d, which is much closer to the cylinder surface.
Meanwhile, at the same distance from the cylinder surface, the optimal model shows the large flow velocity, which is consistent with the turbulence distribution at this location. Thus, the maximum flow velocity increase in turn from model A to model B and optimal model.

The pressure coefficient C p
To confirm the influence of the studied parameters on the backward delay of the full separation point, the full separation points of these different models at Re d = 5.83 × 10 4 need to be precisely determined through numerical simulation. Therefore, the pressure coefficient C p was introduced to determine the location of the full separation point, which was calculated using the following equation: where P is the static pressure on the cylinder surface, P 0 is the static pressure on the inflow boundary surface, ρ is the air density and U is the uniform flow velocity. Yamagishi drew a conclusion that when the partial separation and reattachment of turbulent flow repeated in the area nearby the grooves, the value of the relative pressure coefficient C p changes periodically on the cylinder surface in the upstream side, the pressure becomes large in the valley of the groove and small in the peak of the groove. In additional, the backpressure coefficient in the rear area becomes almost invariable after a certain position point as the flow would not reattach to the surface any more, and this position point can be considered as the point of full separation (Yamagishi & Oki, 2007). Thus the pressure coefficient C p can be used to predict the position where flow separating. The computational results of the pressure distributions on models A, B and optimal model surfaces at Re d = 5.83 × 10 4 are shown in Figure 15. The ordinate shows the pressure coefficients C p and the abscissa shows the angle from the forward stagnation point. Figure 15 demonstrates that the point where the pressure coefficient C p of model B begins to stabilize is at about 105°from the forward stagnation point, larger than the approximal 103°for model A, which is corresponding to the position of the boundary layer separation. And this result is similar to the ones reported in Yamagishi's experiment (Yamagishi & Oki, 2007). Also, the backpressure coefficient of model B becomes larger compared with that of model A near the back surface of the cylinder that is against to incident flow, verifying the smaller pressure drag in model B. Meanwhile, the full separation point of the optimal model is about 110°accompanied by the largest backpressure coefficient, which indicates a better effect of optimized grooves on drag reduction.

The velocity gradient
The velocity gradient near the cylinder surface is affected by the combined effect of the boundary layer as mentioned above, and the value of velocity gradient (du/dz) = 0 corresponds to the partial separation or reattachment of the flow at this location where the reverse flow and the small vortex happened. More precisely, the velocity gradient below 0 corresponds to the  Figure 16. The full separation points determined by this method are also marked in the figure, which is consistent with the distributions of the pressure coefficients.

The vortex shedding
The delay of vortex shedding and the smaller width of the Karman vortex street caused by the retardation of the separation point are the key reasons for drag reduction. For verifying the size and the shape of the shedding vortices, visual simulation analysis based study of flow characteristic behind the cylinder is required. The computational results of the vorticity distributions behind models A, B and optimal model with Re d = 5.83 × 10 4 are shown in Figure 17(a) (b) and (c), which clearly demonstrates that the value of vorticity and the region of wake in model B become smaller compared with that in the models A. The same performance also occurs in the comparison analysis between model B and optimal model. And this characteristic of vortex shedding behind the cylinder is mainly responsible for the reduction of pressure drag.
The computational results of the vorticity at x/d = 1.1 behind the cylinder of the three models at the same time point of T = 0.2 s are shown in Figure 18. The trends of the numerical distribution are in close agreement with the simulation results: the vorticity behind the cylinder of model B becomes smaller at approximately 2.6% compared with that of model A by quantitative analysis, and the position where the value of vorticity reaches the maximum shifts to the low y/d, which means a small region of shedding vortices. The percentage difference in the position of maximum vorticity between model A and model B is about 5.8% in the z-direction. In the case where there exist the optimal grooves on the cylinder surface, the distribution of vorticity tends to be more in the middle of flow field compared with model A and model B, which shows successful optimization on drag reduction as the percentage difference between model B and optimal model is about 4.1%.
In addition, it is also feasible to characterize the width of the wake qualitatively using the distribution of turbulent kinetic energy at a certain position behind the cylinder. The computational results of the distribution of turbulent kinetic energy k/U 2 at x/d = 1.1 behind the cylinder are shown in Figure 19, which are also in agreement with the aforementioned results. The value of turbulent kinetic energy behind the cylinder of model B becomes smaller compared with that behind the cylinder of model A in quantitative analysis, the difference is about 6%, which proved the grooves set at the rear surface of the cylinder may cause larger turbulence. The percentage difference in wake size between model A and model B is about 6.4%, and the difference is about 6.8% between model B and the optimal model.

The backpressure
The computational results of the gauge pressure behind the cylinder of the three models at the same time point of t = 0.2 s are shown in Figure 20. Comparing the backpressure of model A and model B, it is evident that the pressure at the rear of the cylinder in model A is significantly smaller than that in model B, showing a larger area of adverse pressure, and the region of adverse pressure tends to be closer to the surface of cylinder. The different value of the backpressures causes the difference of the pressure drag between models can be explained as follows: the grooves arranged on the back surface of the cylinder also trigger the turbulence of flow and generate the vortices in the near region of the cylinder surface. Therefore, in model A, the velocity of the flow behind the cylinder is larger, and also the backpressure is lower compared with that in model B. That is the main reason why the grooves oriented in the direction against the incident flow have a negative effect on drag reduction.
The larger backpressure of the cylinder in the optimal model is mainly because of the retardation of full separation point caused by the optimized grooves, resulting in the delay of vortex shedding and the reduction of high vorticity area behind the cylinder. Thus the backpressure of the cylinder in an optimal model is larger and the drag force is smaller than that of the cylinder in model B.

The power spectra
As the frequency in generating the vortex is calculated from the fluctuation of the lift coefficient, Figure 21 shows the power spectra of C l for models A, B and optimal model. It is clear that the frequency of the generating vortex of model A becomes smaller compared with that of model B.
Based on the simulation results it can be found that the Strouhal number St (St = (f × d/u), f : frequency, d: diameter of the cylinder, and u: uniform velocity) of model B becomes larger compared with that of Model A, which also exists in the comparison between model B and optimal Model.

PIV tests in the wind tunnel
To validate the reliability of simulation experiment, the PIV technology was used in a wind tunnel to measure the flow characteristics, since the PIV technology can record the information of velocity distribution from a large number of space points at the same transient state.
Considering that the drag reduction of flow around the cylinder is mainly related to the position of the full separation point, the verification was implemented though the comparison of this parameter which is obtained by simulations as discussed earlier. A low-speed, return-flow wind tunnel with the test section of 0.3 m × 0.3 m was adopted in the PIV test. In order to accurately capture the distribution of velocity vector near the cylinder surface, the trace particles made of atomized oil produced by a particle generator were ejected into the wind tunnel and displayed under the laser irradiation. The PIV setup and post-processing related insight package is consist of two resonator laser for light emission, one high resolution CCD camera for image acquisition, one synchronizer for control coordination, and one flow display system. This insight package then combines with the external interface    to give rise to the whole test device. The structure of the wind tunnel and the testing devices are shown in Figure 22. When the velocity of airflow near the cylinder surface decreases to a certain level under the action of the adverse pressure, the flow reverses its direction (compared to the main flow) and thus recirculation occurs. As the two kinds of flow from different directions meet, the flow in the boundary layer separates off the surface, and generates a large number of vortices. Different from the partial separation caused by grooves, this flow separation is more ferocious accompanied with larger vortices at the downstream side, and is not reattaching to the surface anymore. To present this phenomenon clearly, the streamlines visualization is shown in Figure 23.
For verifying the accuracy of this simulation result, an original image of the velocity vector field distribution in the x-z plane of an optimal model by the PIV test is shown in Figure 24.
The position of the final separation point in this experiment is about 110 degrees from the forward stagnation point, which agrees with the simulation result, proving that the high credibility of numerical simulation.

Achievements
The main conclusions and achievements are the following: (1) The mechanism of drag reduction in the flow passing the surface of a grooved cylinder was discussed. As the full separation point is decayed by the nonsmooth surface, the region of adverse pressure will be reduced. This is different from the mechanism of drag reduction when flow passes a plain surface with grooves, as the low-speed roller-bearing-like swirl is formed within the non-smooth unit, which translates sliding friction between the object surface and fluid into rolling friction, leading to skin friction reductions (Song, Lin, Liu, & Zhou, 2017;Song, Zhang, & Lin, 2017).
(2) The semi-grooved surface has more significant effects on drag reduction than the fully grooved pattern, with a maximum ratio of drag reduction up to 5.42% when the fluid speed is 35 m/s. The grooves oriented in the direction against the incident flow show a negative effect on the drag reduction as they generate the turbulence of flow and form the vortices near the rear area of the cylinder, which leads to lower backpressure. (3) An equation of quadratic response surface regression was established to demonstrate the relationship between the drag coefficient (C d ) and the parameters of the grooves, which can be used for describing the different effects of parameters on drag reduction more accurately providing a reference for engineering design. (4) The Depth (factor A) has the greatest impact on the solution, followed by the effect of the Angle α between the rearmost groove and forward stagnation point (factor C), and then the Width (factor B).
The experiments of response surface verify the interactive relationship between the Depth (A) and the Width (B), which explains that the effect of one factor on drag reduction shows inconsistencies at each level of the other one and offers a high reference value for engineering application on drag reduction when a fluid flow around a grooved cylinder. (5) An optimized groove structure was proposed with a better geometric structure to make the circle vortex smaller in the groove, which makes the disturbance of flow and the momentum transfer adhere to the surface of the cylinder more broadly and leads to the higher kinetic energy of flow at this place. The position of the rearmost groove was also discussed with respect to its delay the full separation point at most. (6) The position of full separation point in the optimal model was captured though the velocity vector field distribution in a low-speed and return-flow wind tunnel, and the simulation results were validated by the PIV experiment with respect to its accuracy in predicting the location of separated flow.

Limitations and improvements
The influence of different factors is explored and discussed for model optimization only focus on the effect of groove structure on drag reduction. In a practical application, there are some other factors that may influence the effect of drag reduction such as the interaction between multiple cylinders and the spacing of neighbouringgrilles, the starting angle of grooves, which was not considered. The effect of these factors on drag reduction will be investigated in future research. Besides, the study is limited to the engineering application of highway speed air flow, and it remains to be investigated whether it can be extended to the low or medium speed aeronautical applications. At the same time, the flow characteristics obtained by the study can be further explored in the water flow to extend the scope of potential applications.