Aerodynamic Shape Optimization Method of Non-Smooth Surfaces for Aerodynamic Drag Reduction on A Minivan

Abstract

To reduce aerodynamic drag of a minivan, non-smooth surfaces are applied to the minivan’s roof panel design. A steady computational fluid dynamics (CFD) method is used to investigate the aerodynamic drag characteristics. The accuracy of the numerical method is validated by wind tunnel test. The drag reduction effects of rectangle, rhombus and arithmetic progression arrangement for circular concaves are investigated numerically, and then the aerodynamic drag coefficient of the rectangle arrangement with a better drag reduction effect is chosen as the optimization objective. Three parameters, that is, the diameter D of the circular concave, the width W and the longitudinal distance L among the circular concaves, are selected as design variables. A 20-level design of an experimental study using a Latin Hypercube scheme is conducted. The responses of 20 groups of sample points are obtained by CFD simulation, based on which a Kriging model is chosen to create the surrogate-model. The multi-island genetic algorithm is employed to find the optimum solution. The result shows that maximum drag reduction effects up to 7.71% can be achieved with a rectangle circular concaves arrangement. The reduction mechanism of the roof with the circular concaves was discussed. The circular concaves decrease friction resistance of the roof and change the flow characteristics of the recirculation area in the wake of the minivan. The roof with the circular concaves reduces the differential pressure drag of the front and rear of the minivan.

1. Introduction

Due to the energy crisis, the field of automotive engineering is now more concerned with decreasing the fuel consumption of road vehicles. Reducing the aerodynamic drag is crucial for improving fuel efficiency Currently, the improvement of aerodynamic drag characteristics is mainly achieved through methods such as streamlining and local improvement of the body. Wolf-Heinrich Hucho achieved drag reduction by modifying vehicle local body shape [1]. Rong Jianglei put forward a local retrofit method of aerodynamic characteristics based on an iterative method, and verified its validity with the aerodynamic drag improvement practice of a car [2]. However, these methods have some limitations in terms of further aerodynamic drag reductions. Modifications of the model are based on trial and error experience, and the CFD calculations also need further validation. If the modified model is still unsatisfactory, then trials must be repeated for further modification [3]. Therefore, this method is not only time-consuming and less efficient, but also blind and arbitrary; it is very difficult to obtain the best results.

The drag reduction research on non-smooth surface structures has received a great deal of attention [4,5,6]. Drag force can be divided into two different types, that is, pressure drag and friction drag. Different drag reduction methods contain different drag reduction mechanisms. Whitmore and Naughton reported a net drag reduction by introducing roughness to the forebody skin using micro-machined surface overlays in the blunt-based flight vehicle [7]. This can be attributed to microscale structures which alter the near wall turbulence of turbulent boundary layers. As forebody roughness increases, the boundary layer at the model aft thickens and reduces the shearing effect of external flow on the separated flow behind the base region, resulting in reduced base drag. Bullen and Norman measured qualitative and quantitative aspects of the head and body pelage of 23 species of Western Australian bats [8]. He found that rib lets can reduce the skin friction drag of the head and body by up to 10%. Gad-el-Hak pointed out that the proper use of compliant coating as a potential method to favorably interfere with the wall-bounded flow can affect laminar-to-turbulence transition, control boundary layer separation, modulate flow-induced noise, or alter skin friction drag in both laminar and turbulent flows [9]. This can be attributed to hydrophobic surface structures which reduce skin friction drag in liquids. The effect of a non-smooth surface on the flow characteristics near the model surface at specific positions was investigated [10]. The numerical results showed that the circular concave surface was the most promising. The application of a golf dimples surface is the most typical for reducing drag reduction. This can be attributed to the dimples on the golf ball surface which trigger laminar-turbulent transition of the boundary layer. While the laminar boundary layer of a smooth ball is prone to separation in the presence of positive pressure gradients, the turbulent boundary layer remains attached much further downstream.

The successful application of dimples on golf balls have inspired vehicle aerodynamics engineers for a long time. Up to now, the effects of the arrangement of circular concave non-smooth surfaces on vehicle drag reduction have not been well studied. We must note that the aerodynamic shape modifications are severely constrained by the automobile design, and the aerodynamic performances. Only small shape modifications are generally acceptable, and the challenging task for the manufacturer is to select a pertinent modification. In recent years, many design optimization methods have been widely used in aerodynamic area, which can bring great convenience to aerodynamic shape design. The application of a non-smooth surface with optimal physical properties in vehicle drag reduction is very prospective.

The purpose of this study is to investigate the aerodynamic drag improvement effect by applying the Surrogate Model to the circular concave shape optimization of a minivan with a non-smooth surface. In the first section, the numerical method employed for CFD simulation is described, including governing equations, meshing and solution procedure, and then the accuracy of the numerical method is validated by wind tunnel test. In the second section, the aerodynamic drag effects of different arrangements of circular concaves on the roof of the minivan were studied. Then, based on the arrangement with the best drag reduction effect, the relevant design variables were selected, the Latin square design of the experimental method was used to select the sample points, the Kriging surrogate-model was built, and finally global optimization was conducted by a multi-island genetic algorithm. In the third section, through the numerical calculation of an external flow field of models, the flow field, velocity vector, turbulent intensity and other parameters of both the smooth and non-smooth models were compared, and the mechanism of aerodynamic drag reduction of the non-smooth roof with a circular concave structure was analyzed.

2. Numerical Method and Validation

2.1. Governing Equations

The three-dimensional incompressible governing equations are considered in this study. Many studies have demonstrated that the realizable $k$–$\epsilon$ turbulence model is more suitable for the simulation of an automobile external flow field [11,12]. However, the realizable $k$–$\epsilon$ turbulence model adapts wall functions to deal with near-wall treatment. This cannot capture the flow feature of dimples on the roof of the minivan. The $k$–$\omega$ SST turbulence model can capture the shear stress distribution correctly in the boundary layer [13]; so, the shear-stress transport (SST) model with low-Reynolds Number correction was employed in this research. The SST model combines the robustness of the $k$–$\omega$ turbulence model near walls and the capabilities of the $k$–$\epsilon$ model away from the walls as a unified two-equation turbulence model. The definition of the turbulent viscosity is modified to account for the transport of turbulent shear stress. The incompressible continuity equation can be expressed with Cartesian index notation by:

$$\frac{\partial u_i}{\partial x_i}=0,$$

(1)

and the Reynolds momentum equation is:

$$\rho \frac{\partial u_i}{\partial t}+\rho u_i\frac{\partial u_j}{\partial x_i}=\rho g_i-\frac{\partial p}{\partial x_i}+\mu \frac{{\partial }^2{u}_j}{\partial x_i{x}_i}+\frac{\partial }{\partial x_i}\left[{\mu }_t\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)\right]-\frac{2}{3}\rho \frac{\partial k}{\partial x_i}$$

(2)

Where $\rho$ is the fluid density and it remains constant in the differential equations as the flow is assumed to be incompressible, $u_i$ is the velocity component, $p$ is the static pressure, $\mu$ is the dynamic viscosity, $k$ is the turbulent kinetic energy, and ${\mu }_t$ is the scalar dynamic eddy viscosity coefficient.

The specific dissipation rate equation of the $k$–$\omega$ SST model is:

$$\rho \frac{\partial k}{\partial t}+\rho u_i\frac{\partial k}{\partial x_i}={\mu }_t\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)\frac{\partial u_i}{\partial x_j}-\rho {\beta }^{*}\omega k+\frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_i}\right]$$

(3)

and the specific dissipation rate equation of the $k$–$\omega$ SST model is

$$\begin{gathered} \rho \frac{\partial \omega }{\partial t}+\rho u_i\frac{\partial \omega }{\partial x_i}=\gamma \rho \left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)\frac{\partial u_i}{\partial x_j}-\rho \beta {\omega }^2 \\ +\frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial \omega }{\partial x_i}\right]+2\rho \left(1-F_1\right){\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_i}\frac{\partial \omega }{\partial x_i} , \end{gathered}$$

(4)

where ${\mu }_t$ is the scalar dynamic eddy viscosity coefficient [13], which is defined as

$${\mu }_t=\rho {\nu }_t=\rho \frac{a_{1}k}{\max \left(a_1\omega , \left|\mathsf{\Omega }\right|F_2\right)},$$

(5)

where $\left|\mathsf{\Omega }\right|$ is the absolute value of the mean vorticity vector, $a_1=0.3$. Menter’s second blending function is $F_2=\tanh \left[\max {\left(\frac{2\sqrt{k}}{0.09\omega y};\frac{500\nu }{\omega y^2}\right)}^2\right]$, ${\beta }^{*}=0.09$,${\sigma }_k$, ${\sigma }_{\omega }$, $\gamma$, $\beta$ are functions defined by Menter [14] as

$${\sigma }_k={\sigma }_{k1}{F}_1+{\sigma }_{k2}\left(1-F_1\right)=0.85·F_1+\left(1-F_1\right)$$

(6)

$${\sigma }_{\omega }={\sigma }_{\omega 1}{F}_1+{\sigma }_{\omega 2}\left(1-F_1\right)=0.65F_1+0.856\left(1-F_1\right)$$

(7)

$$\gamma ={\gamma }_1{F}_1+{\gamma }_2\left(1-F_1\right)=0.469116F_1+0.4403546\left(1-F_1\right)$$

(8)

$$\beta ={\beta }_1{F}_1+{\beta }_2\left(1-F_1\right)=0.075F_1+0.828\left(1-F_1\right),$$

(9)

where $F_1$ is the blending function, which activates the $k$–$\omega$ model in the near wall region and the $k$–$\epsilon$ model for the rest of the flow. All constants are computed by a blend from the corresponding constants of the $k$–$\epsilon$ and $k$–$\omega$ model via $SST=F_1\left(k-\omega \mathrm{model}\right)+\left(1-F_1\right)\left(k-\epsilon \mathrm{model}\right)$.

$$F_1=\tanh \left\{\max \left[\min \left(\frac{\sqrt{k}}{0.09\omega y};0.45\frac{\omega }{\left|\mathsf{\Omega }\right|}\right);\frac{500\nu }{\omega y^2}\right]\right\}$$

(10)

2.2. Calculation Methods

To reduce the calculation cost, the numerical simulation model of the minivan was appropriately simplified. For example, the door handles, windscreen wiper and rain gutters, and so forth, were ignored. A 1:1 model was developed by using UG software, the length × width × height are 3875 × 1620 × 1850 (mm³), respectively, shown in Figure 1. The computational domain of the whole minivan is a cuboid. The computational domain is seven times the width, eleven times the length and five times the height of the CAD model. ANSYS ICEM CFD software was used to generate an unstructured mesh, and mesh refinement was done on the non-smooth surface to obtain the flow field information more accurately. Meanwhile, to avoid the mesh difference on simulation results, the mesh size of the same part of the model remained unchanged in each CFD steady-state simulation. The total elements for the whole minivan generated in each simulation were more than 5,290,000; three layers of prism mesh were inserted in the vicinity of the minivan body. The thickness of the first layer was 1 mm, shown in Figure 1.

The flow fields simulation was carried out using the computational fluid dynamics code ANSYS FLUENT16.2^®. The solver settings are listed in

Table 1. Solver settings.

Function	Setting
Solver	Segregated
Pressure discretization	Implicit 2nd order
Time-stepping	2nd order implicit
Momentum discretization	Central differencing
Pressure-velocity coupling	PISO
Fluid	Air
Temperature	Normal temperature (24 °C)

. The boundary conditions settings are listed in

Table 2. Boundary Conditions settings.

Boundary	Boundary Condition	Value
Inlet	Constant velocity	30 m/s
	Turbulent kinetic energy k	0.0003375
	Specific dissipation rate $\omega$	0.0596412
Outlet	Gauge pressure	0 Pa (gage)
Top and sides	Free slip	-
Road	Moving wall	30 m/s
Minivan surfaces	No slip wall	-

. A constant velocity boundary condition was chosen for the inlet in order to easily replicate the test velocity in experiments.

The quality of the grid directly affects the accuracy of the simulation results and the time needed to complete the computation. It is very important to check grid independence. Nine kinds of grid numbers were considered to check grid independence, which are shown in Figure 2. With the increase of grid number, the cure of drag coefficient declines. When the grid number increases to about 5.1 million, the cure level goes off, and grid-independent results can be gained. So, 5 million grids have relatively high simulation accuracy in this research.

2.3. Wind Tunnel Experiment

The accuracy of the CFD simulation method is validated by wind tunnel test. The test model is a model of 1:4 by a CNC machining center according to the CAD model so as to ensure the consistency of physical models for testing and the CAD model for numerical simulation. The experimental test is performed on the vehicle test model, which is about 1/4 of the scale of the simulation model. The wind tunnel test model of the original minivan is shown in Figure 3a. The force-measurement test was carried out in the HD-2 wind tunnel of the Wind Engineering Research Center at Hunan University, and the aerodynamic force of the model was tested by six-component and floating box type force balance. Startup of the ground boundary layer suction device can eliminate the effect caused by the ground boundary layer during the wind tunnel test. The wind tunnel is a boundary layer wind tunnel of closed-type, single-circumfluence and double-test segments, and the contour of this wind-tunnel is 53 m in length and 18 m in width. It is an HD-2 type wind tunnel with a rectangular test section area of 7.5, 3 m wide and 2.5 m height with a length of 17 m. To eliminate the blockage effect of the wind tunnel, Mercker’s method was employed to correct the measurement results of HD-2 wind tunnel [15].

The test results of the aerodynamic drag coefficient of the model are shown in Figure 3b. From Figure 3b, when the test wind speed exceeds 20 m/s, the aerodynamic drag coefficient does not change with Reynolds number. The results of the test on the scale model in terms of dimensionless aerodynamic coefficients are the same as for the original minivan if the Reynolds numbers are the same. When the wind speed is 30 m/s, compared with the experiment data, the relative error of the simulation result is 1.43%, which is within the engineering permissible error range of 5%, shown in

Table 3. Comparison of the CFD simulation result and the experiment data.

	Experiment Data	Simulation Result	Relative Error
Cd	0.3570	0.3621	1.43%

, and so the reliability of the CFD numerical simulation is validated.

3. Optimal Design Method

3.1. Non-Smooth Processing Region Selection and Size Settings

Vedrtnam pointed out that the rear part of the vehicle contributes 80% of the total aerodynamic drag [16]. The non-smooth processing region should be on the roof panel, which is able to control the wake region properly so as to reduce the loss of turbulent energy and differential pressure resistance, while the minivan roof panel is the surface which has the greatest impact on wake region. Therefore, the drag reduction effect after the circular concave processing of the roof panel is studied. The non-smooth region and part mesh of the non-smooth surface with triangular cells is shown in Figure 4.

Relevant studies indicate that both the differential pressure resistance caused by airflow separation and the frictional resistance caused by the viscous role of air always relate to the boundary layer thickness [17]. The selection of the non-smooth structure should be related to boundary layer thickness. The size, height or depth of the non-smooth structure should be smaller than the distance between the roof surface and logarithm law area. Consequently, the size of the circular concave is mainly determined according to the thickness of the boundary layer.

The thickness of the flat plate laminar flow boundary layer is calculated as follows:

$$\delta \left(x\right)=5.48x/\sqrt{\frac{\rho Ux}{\mu },}$$

(11)

where $\delta \left(x\right)$ is the thickness of the boundary layer, $x$ is the distance from the leading edge of the plate, the roof panel studied is approximated to be a flat plate with the length of 2.04 m and $\mathrm{Re}\left(l\right)$ is the Reynolds number,

$$\mathrm{Re}\left(l\right)=\frac{\rho Ul}{\mu },$$

(12)

where $\rho$ is density, $U$ is a velocity scale, $l$ is a length scale and $\mu$ is the dynamic viscosity, which is 0.0722 m²/s. The velocity scale is taken as the inlet velocity, taken as 30 m/s in this study.

It is mainly considered that the depth of the circular concave structure should be less than the thickness of the boundary layer during the size selection of the circular concave structure in this study. When calculating, the depth of circular concave structure should not be more than 27.25 mm when circular concave non-smooth processing is performed on the roof panel.

During the CFD simulation, the drag reduction effect is mostly assessed by resistance-reducing ratio ΔC_d%:

$${\mathsf{\Delta }\mathrm{C}}_{\mathrm{d}}\%=\left({\mathrm{C}}_{\mathrm{d}-modified}-{\mathrm{C}}_{\mathrm{d}-original}\right)/{\mathrm{C}}_{\mathrm{d}-original}.$$

(13)

3.2. Size Design and Arrangement of Circular Concave Structure

The size parameter of the circular concave structure was considered: diameter D, the lateral distance W, the longitudinal distance L, and the depth S. To facilitate the design and arrangement, the depth S of the circular concave structure is half of the circular concave diameter D. The size design of the circular concave structure is shown specifically in Figure 5. The value range of D, W, L and S is given as [10, 50], [60, 160], [60, 160], [5, 25] respectively; its numeric type is integer and unit is mm.

According to many bionics experiments, the soil animal, the dung beetle, can move freely in the soil due to the circular concave shape of its non-smooth surface and the arrangement of its circular concave structure. Therefore, the arrangement impact of the circular concave structure should be considered. As shown in Figure 6, three arrangements are chosen: rectangle arrangement, rhombus arrangement and arithmetic arrangement.

The D, W, L and S are selected to be 15 mm, 120 mm,120 mm, 7.5 mm, respectively. CFD simulation is carried out for three arrangements.

Table 4. CFD simulation results of three arrangements.

Arrangements	Cd	Drag Reduction Rate ΔCd%
Original	0.3621	-
Rectangle	0.3531	−2.49%
Rhombus	0.3552	−1.90%
Arithmetic	0.3564	−1.57%

shows that the drag reduction effect of the rectangle arrangement is the best with a drag reduction rate of 2.49%.

3.3. Optimal Design of Circular Concave Non-Smooth Surface

According to the CFD simulation results of the three arrangements, the drag reduction effect of the rectangle arrangement is the best, so the rectangle arrangement of the circular concave non-smooth surface is chosen for optimization, the circular concave diameter D, the lateral distance W and the longitudinal distance L are selected as design variables. The objective of seeking the optimal combinations is to achieve the maximum drag reduction effect, that is, to obtain the minimum C_d value in the rectangle arrangement. The flow chart of the aerodynamic optimization of the circular concave non-smooth surface is shown in Figure 7.

Optimal Latin Hypercube is a modified Latin Hypercube, in which the combination of factor levers for each factor is optimized, rather than randomly combined [18]. Optimal Latin Hypercube is a well-established DOE (Design of Experiment) method. A 20-level DOE study is carried out using the Optimal Latin Hypercube method for the above design variables, and the CFD simulations are conducted and 20 sets of response values are obtained, as shown in

Table 5. Sample points and simulation results.

Design Variables			Simulation Results	Drag Reduction Rate
D (mm)	W (mm)	L (mm)	Cd	ΔCd%
10	144	144	0.3585	−0.99%
12	149	60	0.3534	−2.41%
13	139	113	0.3551	−1.92%
15	71	128	0.3595	−0.73%
16	160	86	0.3565	−1.55%
18	60	65	0.3424	−5.45%
19	123	139	0.3558	−1.73%
21	81	92	0.3457	−4.53%
23	65	76	0.3550	−1.97%
24	76	155	0.3540	−2.25%
26	113	102	0.3532	−2.45%
27	107	123	0.3375	−6.79%
29	134	134	0.3554	−1.84%
31	155	149	0.3567	−1.48%
32	128	97	0.3566	−1.52%
36	102	118	0.3544	−2.12%
38	92	107	0.3541	−2.21%
42	118	81	0.3363	−7.12%
47	97	160	0.3571	−1.37%
50	86	71	0.3373	−6.85%

. The 20 levels of each design variable are divided between the lower and upper bounds of the design space. With these 20 design points, the minivan model is updated, meshed and solved for aerodynamic drag automatically. Then, DOE post-processing of the Pareto analysis is performed. To rank each parameter’s contribution, a Pareto plot for drag coefficient is given in Figure 8.

Figure 8 shows that the longitudinal distance L has the greatest impact on aerodynamic drag coefficient C_d, the effect of which on C_d is about 22%. The diameter D and the lateral distance W have the reverse trend of effects, 2% and 4%, respectively, and so the longitudinal distance L is the major factor affecting the aerodynamic drag coefficient. The interaction effect between design variables is shown in Figure 9. The unit of horizontal coordinate in Figure 9 is mm.

The interaction of L with W on C_d is the most obvious one, followed by D with W; the interaction of D with L is the smallest. Although the impact of W on the aerodynamic drag is not particularly significant, the impact of interaction effects between W and other parameters on C_d cannot be ignored.

Based on the DOE results, a surrogate-model is built. Surrogate-model concepts were introduced in structural design optimization by Schmit and Miura in the late 1970s [19], and were then widely used in the engineering design area. One of the primary reasons for using the surrogate-model is to reduce the computational cost of performing optimization. Several multi-dimensional and non-linear interpolation techniques can be used to construct the surrogate-model. The surrogate-model is built using the DOE sample points in order to construct an analytical relation between the design variables and the simulation responses. Kriging models are originated in the areas of mining and geostatistics that involve spatially and temporally correlated data. Kriging models can provide an accurate surrogate-model. Moreover, it allows the construction of a global surrogate-model which is valid for the entire design space. Therefore, the Kriging Model is chosen to create the surrogate-model in this study.

To verify the fitting precision of the surrogate-model, another three design points are chosen beside the 20 DOE design points. The three design points are used for the CFD simulation and the surrogate-model simulation, respectively; the results are shown in

Table 6. Comparison of C_d between CFD simulation results and surrogate model results.

Design Variables			Comparison
D (mm)	W (mm)	L (mm)	CFD Simulation	Surrogate Model	Relative Error (%)
28	110	100	0.3512	0.3491	0.60
36	90	80	0.3531	0.3505	0.74
44	95	100	0.3501	0.3466	1.00

. The relative error is the ratio of the absolute error and the CFD simulation result. The absolute error is the absolute value of the difference between the surrogate-model result and the CFD simulation result.

From Table 6, the relative error between the CFD simulation result and the surrogate model value is within 2% in each group, which indicates that the established surrogate model can describe the relationship between design variables and the response values. The surrogate model has a higher reliability, and it is feasible as a replacement for the CFD simulation results.

3.4. Optimization Method Using Multi-Island GA

Genetic algorithm (GA) is a global optimization method and is used to select the optimum solution candidates and control the optimum inquiry process. GA is built on the theories of Darwinian biological evolution and Mendelism [20]. It performs genetic operations such as selection, crossover, and mutation on the fitness individual and becomes the new optimum individual. This method of optimization has been applied successfully to many problems in engineering optimization [21]. A more efficient Genetic Algorithm called “multi-island GA” is employed in this study. The main feature of multi-island Genetic Algorithm that distinguishes it from traditional genetic algorithms is the fact that each population of individuals is divided into several sub-populations called islands. Multi-Island Genetic Algorithm, as a kind of pseudo-Parallel Genetic Algorithm, can search the global optimal solution in the optimization area well [22]. The basic parameters used in the multi-island GA are shown in

Table 7. Basic Parameters of multi-island GA.

Parameters	Values
Size of Subpopulation	50
Number of Islands	10
Iterative algebra	100

. Based on the surrogate-model and the multi-island GA, the optimal design variables and objectives are calculated. Finally, the optimal design variables are as follows: W is 110 mm, L is 72 mm and D is 45 mm, respectively.

4. Results and Discussion

The optimal design variables are put into the CFD simulation of the optimization process (shown in Figure 7) again, the optimization minivan CAD model is constructed, the mesh is generated by ICEM CFD and solved for aerodynamic drag coefficient with ANASYS Fluent software; the simulation result for aerodynamic drag coefficient is 0.3342. Compared to the surrogate-model result, the tolerance is only 0.57%, shown in

Table 8. Optimal solution.

Design Variables			Comparison of Aerodynamic Drag Coefficient
D(mm)	W(mm)	L(mm)	CFD Simulation Result	Surrogate-Model	Relative Error (%)
45	110	72	0.3342	0.3323	0.57

. That validates the precision of the surrogate-model again. Relative to the minivan with the original circular concave surfaces, the objective is decreased 7.71%, as shown in

Table 9. Comparison for aerodynamic drag coefficient.

Original Model	Optimized Model	Reduction Rate (%)
0.3621	0.3342	7.71

.

When the air flows through the surface of the circular concave, a low-speed rotating air flow is formed inside the circular concave, as shown in Figure 10. The rotating air flow in the circular concave leads to the gas–gas contact between the airflow inside and outside the circular concave, resulting in the vortex cushion effect and the driving effect. These vortices form “air bearing”, which converts sliding friction into rolling friction. Under the combined action of the vortex cushion effect and the driving effect, the decrease of velocity gradient means that the thickness of the viscous bottom layer increases, which leads to the decrease of viscous shear stress. As Ding pointed out, these vortices act like a kind of fluid roller bearing and may reduce the skin’s friction drag [10]. As a result, the viscous drag is reduced.

The velocity distribution cloud diagram and the isoline of the boundary layer thickness can be seen in Figure 11. The green line represents the isoline of the boundary layer thickness ($\delta =0.99U_{\infty }$). The thickness of the boundary layer of the smooth roof is about 1.2 mm, while the thickness of the boundary layer of the roof with circular concaves is higher than 2.3888 mm, which shows that the circular concave structure does produce the effect, which increases the thickness of the boundary layer, reduces the average velocity gradient of the normal at-the-wall surface, and the surface shear stress is reduced. As Whitmore et al. point out, when the roof roughness increases, the boundary layer on the roof thickens and reduces the shearing effect of external flow on the separated flow behind the base region, resulting in reduced base drag [7].

Some positions between the circular concaves are selected for turbulent intensity analysis. The positions are marked as the distance from the first row of circular concaves. The profiles of the turbulence intensity in the direction of perpendicular to the wall, for the smooth wall and the wall with circular concaves, are shown in Figure 12. The turbulent intensity near the wall of the circular concave structure’s surface appears to be larger than that in the original model. The turbulence intensity is increased by the presence of the circular concaves. When the distance L between circular concaves decreases, the roughness on the roof panel increases, which means that the turbulence intensity increases, so the drag reduction is more prominent when the design parameter L is reduced, shown in Figure 8.

Some points between the circular concaves are selected for mean skin friction coefficient analysis. The positions are marked as the distance from the first row of circular concaves. The mean skin friction coefficient on some point is shown in Figure 13. Compared with the original minivan, the wall friction coefficients decrease after the arrangement of circular concaves is optimized. Drag force can be divided into two different types, that is, pressure drag and friction drag. Skin friction reduction plays an important role in this research. The increase of turbulence intensity at the roof of the minivan leads to the decrease of flow shear force, and thus decreases friction resistance.

Figure 14 shows two vector pictures of the longitudinal symmetry plane of the original minivan and the optimal minivan, respectively. It can be clearly seen that the airflow separates forward on the rear roof panel of the optimized minivan, and there is no eddy behind the optimized minivan. However, there are two dominant eddies behind the original minivan. As Wu et al. point out, the reattachment length of a two-dimensional turbulent Backward-Facing Step flow decreases as the thickness of the incoming boundary layer increases [23]. The thickness of the boundary layer of the roof with the circular concaves is higher than that of the smooth roof (show in Figure 11), which leads to a decrease of the size of the recirculation area in the wake of the optimized minivan (shown in Figure 14).

The pressure drag due to the viscosity of the air causes the differential pressure between the front and rear of the minivan, and is a major component of the aerodynamic drag. To reduce the pressure drag is to reduce the positive pressure area in front of the minivan and the negative pressure zone behind the minivan. Figure 15 shows the wake region pressure contours of the original minivan and the optimized model, respectively. The wake negative pressure area of the optimized car reduces significantly, the positive pressure area increases and thus reduces the differential pressure drag of the front and rear.

5. Conclusions

The aerodynamic drag coefficients of the minivan with and without the non-smooth surfaces are simulated using CFD. Simulation results show that circular concave non-smooth processing on the roof’s surface can reduce the aerodynamic drag of the minivan effectively.

Three kinds of circular concave arrangement methods are discussed, that is, rectangle arrangement, rhombus arrangement and arithmetic arrangement, among which the drag reduction effect of the rectangle arrangement method is the best.

An aerodynamic shape optimization method based on the surrogate-model is presented and applied to the case of the circular concave non-smooth surface. Finally, the optimal design variables are obtained: W = 110 mm, L = 72 mm and D = 45 mm, respectively. The relative error between the CFD simulation result and the surrogate-model result is only 0.57% using the optimal design variables, the aerodynamics drag coefficient is decreased 7.71% compared with the drag coefficient of the minivan with the original non-smooth surface. This method, combining the design of the experiment, the approximate model and the optimization algorithm, can provide certain engineering guidance for the research and optimization of vehicle body circular concave non-smooth surface drag reduction.

The reduction mechanism of the roof with the circular concaves was discussed. Compared with the original minivan with the smooth roof, the increase of turbulence intensity at the roof with the circular concaves leads to the decrease of flow shear force, and thus decreases friction resistance. The thickness of the boundary layer of the roof with the circular concaves is higher than that of the smooth roof, which leads to a decrease of the size of the recirculation area in the wake of the optimized minivan, and there is no dominant eddy behind the optimized minivan. The wake negative pressure area of the optimized car reduces significantly, the positive pressure area increases and thus reduces the differential pressure drag of the front and rear.