An Arbitrary Hybrid Turbulence Modeling Approach for Efficient and Accurate Automotive Aerodynamic Analysis and Design Optimization

Abstract

The demand in solving complex turbulent fluid flows has been growing rapidly in the automotive industry for the last decade as engineers strive to design better vehicles to improve drag coefficients, noise levels and drivability. This paper presents the implementation of an arbitrary hybrid turbulence modeling (AHTM) approach in OpenFOAM for the efficient simulation of common automotive aerodynamics with unsteady turbulent separated flows such as the Kelvin–Helmholtz effect, which can also be used as an efficient part of aerodynamic design optimization (ADO) tools. This AHTM approach is based on the concept of Very Large Eddy Simulation (VLES), which can arbitrarily combine RANS, URANS, LES and DNS turbulence models in a single flow field depending on the local mesh refinement. As a result, the design engineer can take advantage of this unique and highly flexible approach to tailor his grid according to his design and resolution requirements in different areas of the flow field over the car body without sacrificing accuracy and efficiency at the same time. This paper presents the details of the implementation and careful validation of the AHTM method using the standard benchmark case of the Ahmed body, in comparison with some other existing models, such as RANS, URANS, DES and LES, which shows VLES to be the most accurate among the five examined. Furthermore, the results of this study demonstrate that the AHTM approach has the flexibility, efficiency and accuracy to be integrated with ADO tools for engineering design in the automotive industry. The approach can also be used for the detailed study of highly complex turbulent phenomena such as the Kelvin–Helmholtz instability commonly found in automotive aerodynamics. Currently, the AHTM implementation is being integrated with the DAFoam for gradient-based multi-point ADO using an efficient adjoint solver based on a Sparse Nonlinear optimizer (SNOPT).

1. Introduction

The automobile industry is still developing, and the number of cars is increasing day after day [1]. With the growing demand for automobiles, there appears a question of designing the vehicles with better aerodynamic characteristics, resulting in better fuel economy, lower noise levels and better drivability. Principally, the turbulent flow at the rear part of the car has a strong influence on the above due to the wake that occurs because of the flow separation in that region. Sivaraj, G., et al. [2] define flow separation as a condition due to the lack of energetic flow and inability of flow to move over sharp edges as shown in Figure 1. In turn, it causes shear layer instability or the so-called Kelvin–Helmholtz instability, and thus noise generation.

There have been many studies that have analyzed that concept using Reynolds-Averaged Navier–Stokes (RANS), Unsteady Reynolds-Averaged Navier–Stokes (URANS), Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) models, which are normally validated by comparison of their predictions with the corresponding experimental data.

A typical benchmark case widely used for flow model validation in the automotive industry is the experimental study by Ahmed et al. [4], where they performed an experiment in a wind tunnel using the Ahmed body as a simplified vehicle model. Igali et al. [5] used the Ahmed body as a test case to compare the performance and efficiency of different turbulence models for automotive aerodynamic application and concluded that the RANS approach with the k-ω SST model is still sufficiently accurate in calculating the time-averaged parameters in contrast to the earlier general claim of the model’s inability to capture the Kelvin–Helmholtz effect adequately by Menter [6], the inventor of the model.

Ahmed et al. [4] described their simplified car body as a model that consists of three parts, such as a forebody, midsection, and rear end. The edges are rounded in order to avoid flow separation. The midsection is a sharp-edged box that has a rectangular cross-section. The slant angle in the rear can be configured from 0 degrees to 40 degrees with a 5-degree step. However, in this study, the Ahmed body with a 25-degree slant angle was investigated.

According to Ahmed et al. [4], the experimental value for the drag coefficient is equal to 0.298, and the wake structure was obtained and analyzed behind the body.

At this slant angle, the flow has a transient and complex topology. As can be seen in Figure 2, two vortices are produced on the side edges (red lines). In addition, the flow separates from the slant surface at the sharp corner and then reattaches near the downstream end of the slant surface of the body in Figure 2, where there is an unstable shear layer on top of the created separation zone above the slant surface. This results in the Kelvin–Helmholtz instability of the shear layer there. Similar instability can also happen on top of the separation zone behind the rear end.

This study aimed to implement the Arbitrary Hybrid Turbulence Modeling (AHTM) approach based on the Very Large Eddy Simulation (VLES) in OpenFOAM/DAFoam to explore its flexibility for efficient analysis and design with automotive applications in mind. This could help the design engineers to taylor their his mesh for different requirements, such as design optimization for minimum drag, low level of noise and good drivability, in different parts of the flow field to achieve the required accuracy and resolution with sufficient efficiency at the same time for rapid turn-around time and improved productivity.

2. Literature Review

There are many published studies in the literature that present the results of simulations aimed to provide insight into the flow around the Ahmed body. For example, two recent studies showed the wake structures of the vehicle in detail by Igali et al. [5] and Sadykov et al. [8], respectively, with the former using the RANS approach and the latter the URANS one. Their usefulness is in the comparison of results obtained by the RANS and URANS approaches with different turbulence models, where Igali et al. [5] achieved 2.474% of error in comparison with experimentally measured drag coefficient using the RANS approach, while Sadykov et al. [8] obtained 1.68% error using the URANS method. Moreover, each of the papers has a detailed flow visualization of the wake region at the back of the vehicle model.

However, there is also another approach called AHTM, based on the Very Large Eddy Simulation (VLES), which can show more details and greater resolution of the turbulent eddies and structures in the forms of streamlines and vorticity contours in the rear of the Ahmed body. This approach can arbitrarily combine, the RANS, URANS, LES and DNS models, depending on the local mesh density; therefore, it has the accuracy of LES and computational efficiency of URANS because the user can optimally control the mesh density at critical areas of the flow field to improve simulation resolution and accuracy without sacrificing the efficiency of computation (see Figure 3).

As a result, it can capture the turbulent eddies near and at the rear of the Ahmed Body, which is crucial for properly capturing the Kelvin–Helmholtz instability of shear layers there. However, as the VLES model is a relatively new approach, there are no crucial studies that have been performed to evaluate its performance and efficiency for its potential application in Aerodynamic Design Optimization (ADO) for automotive vehicles.

The turbulence modeling approaches described above, such as RANS and URANS, are successful at predicting the flow behavior for the slant angle of 35°; however, these models are not capable of producing decent results at 25 degrees of slant angle [10]. Furthermore, they cannot resolve the turbulent eddies over the body and in the wake for noise reduction design and optimization. This can be explained by the complexity of the flow structure related to the Kelvin–Helmholtz instability. At this angle, the partial detachment of the flow on the slant develops larger scale unsteadiness in the turbulent flow, thus increasing the complexity of predicting the flow behavior in its numerical simulation. The RANS simulation, for the most part, does not properly describe the separation at the rear part of Ahmed Body [10] due to the time-averaging used for all turbulent eddies; therefore, new turbulence modeling should be considered and implemented to obtain more accurate results. Considering the aforementioned conclusion, LES is thus considered to be a more suitable approach for a body with a 25-degree slant angle, according to [11], which managed to produce results close to the experimental ones. However, LES simulation required extremely fine meshes to resolve the separated turbulent boundary layers around the Ahmed Body. It has serious problems resolving the turbulent boundary layers in the near-wall region. As such, the Detached Eddy Simulation (DES) [12] was introduced to overcome these LES difficulties. However, DES, as a hybrid zonal approach, has a very rigid two-zone demarcation in the flow field introduced before the flow simulation, which makes it very inflexible and difficult to use for engineering applications. In addition, DES is still a largely LES approach in most parts of the flow field. As a result, Both LES and DES have very high computation costs, which hinders the practical application of DES and LES for high Reynolds number flows around Ahmed Body [13] for ADO. It is also noteworthy that Igali et al. [5] and Sadykov et al. [8] showed that the RANS and URANS models could generate time-averaged results that are comparable to those of LES. This is a challenge to the established assumption that the LES/DES is always superior to the more efficient RANS/URANS approaches if the detailed resolution of the turbulent eddies is not required in engineering design.

In order to achieve a compromise between accuracy and computational efficiency, the VLES was developed by Speziale [14]. In this method, the turbulent subgrid stresses were obtained by damping the RANS-based Reynolds stresses in regions where the grid spacing is small and of the order of the Kolmogorov length scale. This method possesses the advantage of continuously transiting between DNS, LES, URANS and RANS computations. However, it excessively damps the Reynolds stresses, and its damping function does not include the turbulence integral scale, which can easily revert to RANS/URANS mode at very large Reynolds numbers even with very fine meshes. Afterward, Johansen et al. [15] developed a model that introduced a filtering factor based on the ratio of modeled length scales of turbulent eddies to local mesh sizes. This method was then used to dampen not the Reynolds stresses as in Speziale’s model [14] but directly dampen the turbulent eddy viscosity. On the other hand, Johansen’s model failed to include the Kolmogorov scale and the ratio of modeled and resolved turbulence energy in the damping of the turbulent viscosity, which is physically more meaningful [16]. VLES models were found to better capture turbulent eddies and the separation around the body more accurately than RANS/URANS [17]. Therefore, it can be thought of as an arbitrary hybrid RANS/URANS/LES/DNS turbulence model in one simulation.

The latest VLES model introduces the resolution control function F_r, which modifies the turbulent viscosity from the RANS modeling [18]. This resolution control factor is the ratio of sub-grid turbulent stress to RANS/URANS turbulent stress, which can also roughly represent the ratio of modeled turbulent energy to total turbulent energy. It is responsible for smooth transitioning between RANS/URANS/LES/DNS modes depending on local mesh density in comparison with turbulence integral and Kolmogorov length scales. The following equation describes the general form of the resolution control function F_r, which is established from the ratio of unresolved turbulence energy to the total turbulence energy [19,20]. The resolution control function represents the ratio between the unresolved turbulent energy to the total turbulent energy [20]. Afterward, it was modified to the following form, which gives the minimum value between 0.0 and 1.0 by adopting the original Speziale model [14] model to calculate the ratio:

$$F_r=min\left[1.0,{\left(\frac{1.0-e^{-\beta \ast L_c/L_k}}{1.0-e^{-\beta \ast L_i/L_k}}\right)}^2\right]$$

(1)

The resolution control function can take values between 0.0 and 1.0, depending on the mesh dimensions, turbulence integral and Kolmogorov length scales, which helps it to avoid the disadvantages of Speziale’s VLES model as discussed above, and can be proven by using the Taylor series to analyze F_r in the extreme condition of L_k tending to zero. Therefore, the VLES model can function in different modes from RANS, URANS, LES and DNS in turbulence modeling, depending on the relative values of the three length scales. For example, control function F_r is equal to 1 when the mesh scale is not capable of resolving the turbulence integral scale, which is the average turbulence length scale; thus, the RANS/URANS simulation is adopted. On the other hand, when the mesh is fine enough to resolve all turbulence scales, the control function reaches 0, thus resulting in DNS modeling [20]. In the middle of these two extremes of the RANS and DNS modeling, the control function takes values of between 0.0 and 1.0, the VLES simulation mode is maintained.

As a result, compared with the standard k-ε RANS/URANS model, the VLES model adjusts the calculation on turbulent eddy viscosity, which is in the following form [17]:

$${\mu }_t=F_r\ast \rho \ast C_{\mu }\ast \frac{k^2}{\epsilon }$$

(2)

Therefore, the VLES model can be developed based on the RANS models such as the k-ω and k-ω SST models. In this study, the VLES based on the k-ω SST turbulence modeling was used for the simulation of aerodynamics over the Ahmed body, whose results were also compared with those of the conventional RANS/URANS k-ω SST turbulence model.

3. Mathematical Formulations

The viscous incompressible flow with constant properties is governed by the Navier–Stokes’ equations [21]:

$$\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+v\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$$

(3)

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

(4)

Since, in most applications, the effect on the characteristics of the mean flow is important, it is convenient to study turbulent flow by dividing it into two parts, which is named as Reynolds Decomposition based on the Reynolds time averaging, where velocity, u, and pressure, p, are expressed as sums of the mean and fluctuating parts [21]:

$$u_i=\overline{u_i}+{u^{\prime }}_i$$

(5)

$$p=\overline{p}+p\prime$$

(6)

Then after inserting expressions above into the Navier–Stokes’s equations, the Reynolds Averaged Navier–Stokes’s equations can be given as follows [22] after dropping all the overbars for convenience:

$$\frac{\partial u_j}{\partial x_j}=0$$

(7)

$$\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+v\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}-\frac{\partial }{\partial x_j}\left\{\overline{u_i{u}_j}\right\}$$

(8)

The Reynolds stresses in the last term on the right-hand side of Equation (8) are modeled by different turbulence models.

The k- SST model

The shear stress transport (SST) k-$\omega$ model is a low Reynolds-number two-equation turbulence model [23]. The current model appears to be the combination of k-$\omega$ model in the viscous boundary sublayer and k-$\epsilon$ model in the regions that are away from walls. The model has two additional transport equations for turbulence kinetic energy, k and specific dissipation rate, $\omega$ [23]:

$$\rho \frac{\partial k}{\partial t}+\rho \frac{\partial \left(u_{i}k\right)}{\partial x_i}=P_k-\beta \rho k\omega +\frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_i}\right]$$

(9)

$$\rho \frac{\partial \omega }{\partial t}+\rho \frac{\partial (u_i\omega )}{\partial x_i}=\frac{\gamma }{v_t}{P}_k-\beta \rho {\omega }^2+\frac{\partial }{\partial x_i}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial \omega }{\partial x_i}\right]+2\left(1-F_1\right)\frac{\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_i}\frac{\partial \omega }{\partial x_i}$$

(10)

$$\begin{gathered} -\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\tau }_{ij}={\mu }_t\left(2S_{ij}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)-\frac{2}{3}\rho k{\delta }_{ij} \\ {\mu }_t=\frac{a_1\rho k}{\max \left(a_1\omega ,\mathsf{\Omega }{F}_2\right)} \\ {P}_k={\tau }_{ij}\frac{\partial u_i}{\partial x_j} \\ {S}_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\right) \\ {F}_1=tanh\left(ar{g}^4\right) \\ arg=min\left[max\left(\frac{\sqrt{k}}{\beta \omega d_w},\frac{500\nu }{d_w^2\omega }\right),\frac{4{\sigma }_{\omega 2}k}{C{D}_{k\omega }{d}_w^2}\right] \\ C{D}_{k\omega }=max\left(\frac{2\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_i}\frac{\partial \omega }{\partial x_i},{10}^{-10}\right) \end{gathered}$$

(11)

Equation (11) represents the Boussinesq analogy used for turbulence model closure. Constants of the equation above have the following values: $\beta =0.09$, ${\sigma }_k=0.85,$${\sigma }_{\omega }=0.5, {\sigma }_{\omega 2}=0.856, a_1=0.31$, $d_w$ is the distance to the wall.

The k-epsilon (k-ε) model

The k-ε model is a semi-empirical model that is based on the transport equations for k and $\epsilon$, which are turbulence kinetic energy and its dissipation rate, respectively [24]:

$$\rho \frac{\partial k}{\partial t}+\rho \frac{\partial \left(u_{i}k\right)}{\partial x_i}=P_k-\rho \epsilon +\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]$$

(12)

$$\rho \frac{\partial \epsilon }{\partial t}+\rho \frac{\partial (u_i\epsilon )}{\partial x_i}=C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]$$

(13)

$${\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon }$$

(14)

Values of model constants are as follows [23]: $C_{\mu }=0.09$, $C_{\epsilon 1}=1.44$, $C_{\epsilon 2}=1.92$, ${\sigma }_k=1.0, {\sigma }_{\epsilon }=1.3$.

The URANS model

The URANS model used is a revised RANS model that describes unsteady turbulent flows with the largest scale above the integral one. Léonard et al. [25] and Sadykov et al. [8] claimed that the URANS method was a good model that could provide a qualitative understanding of the key physics in transient turbulent flows.

In the URANS model, Equation (5) becomes

$$u=\overline{u}+u^{\prime }=\langle \overline{U}\rangle +u^{\prime }$$

(15)

where, $\langle \overline{U}\rangle$ is the transient but time-averaged velocity. The same is also true for pressure.

The VLES model

The following equation describes the general form of control function F_r, which is established from turbulence energy [19,20]:

$$F_r=\frac{{\int }_{L_k}^{L_c}E\left(L\right)dL}{{\int }_{L_k}^{L_i}E\left(L\right)dL}$$

(16)

Here, $L_c=C_{VLES}{({\Delta }_x{\Delta }_y{\Delta }_z)}^{1/3}$ is a cut-off length scale, $L_i=\frac{k^{3/2}}{\epsilon }$ the integral length scale and $L_k=\frac{v^{3/4}}{{\epsilon }^{1/4}}$ is the Kolmogorov length scale. Moreover, the ${\Delta }_x{\Delta }_y{\Delta }_z$ are the mesh dimensions in different directions and the laminar kinematic viscosity is υ. According to Equation (17), the resolution control function represents the ratio between the unresolved turbulent energy to the total turbulent energy [17]. It was subsequently modified to the following form, which adopts the minimum value between 1.0 and modified model [14]:

$$F_r=min\left[1.0, {\left(\frac{1.0-e^{-\beta \ast L_c/L_k}}{1.0-e^{-\beta \ast L_i/L_k}}\right)}^2\right]$$

(17)

The wall functions used in the OpenFOAM simulation are nutKWallFunction, omegaWallFunction and kqRWallFunction.

nutWallFunction

The nutWallFunction provides a wall constraint on various fields, such as turbulent viscosity, i.e., nut, or turbulent kinetic energy dissipation rate, i.e., epsilon, for low- and high-Reynolds number turbulence models.

$$\begin{gathered} v_t=v_{t_{vis}} \mathrm{if} y^{+} <=y_{lam}^{+} \\ {v}_t=v_{t_{log}} \mathrm{if} y^{+} > y_{lam}^{+} \end{gathered}$$

With

$v_t$ = Turbulent viscosity [m²/s]

$v_{t_{vis}}=v_t$ computed by the viscous sublayer assumptions [m²/s]

$v_{t_{log}}$=$v_t$ computed by the inertial sublayer assumptions [m²/s]

$y^{+}$ = Estimated wall-normal distance of the cell centre in wall units

$y_{lam}^{+}$ = Estimated intersection of the viscous and inertial sublayers in wall units

omegaWallFunction

The omegaWallFunction boundary condition provides a wall constraint on the specific dissipation rate, i.e., omega, and the turbulent kinetic energy production contribution, i.e., G, for low- and high-Reynolds number turbulence models.

$$\begin{gathered} {\omega }_{vis}=\frac{6v_{\omega }}{{\beta }_1{y}^2} \\ {\omega }_{log}=\frac{\sqrt{k}}{c_{\mu }Ky} \\ \mathrm{G}=w(v_{t_w}+v_w)\left|n.{\left(\nabla u\right)}_f\right|c_{\mu }^{1/4}\frac{\sqrt{k}}{Ky} \mathrm{if} y^{+} > y_{lam}^{+} \end{gathered}$$

where,

ω = Specific dissipation rate [1/s]

ω_vis = ω computed by the viscous sublayer assumptions [1/s]

ω_log = ω computed by the inertial sublayer assumptions [1/s]

w = Cell-corner weights [-]

k = Turbulent kinetic energy [m²/s²]

y = Wall-normal distance [m]

c_µ = Empirical model constant [-]

v_w = Kinematic viscosity of fluid near wall [m²/s]

$v_{t_w}$ = Turbulent viscosity near wall [m²/s]

n = Face unit normal vector [-]

u = Velocity [m/s]

K = von Kármán constant [-]

ω = ${\omega }_{vis}$ if $y^{+}$ <= $y_{lam}^{+}$

ω = ${\omega }_{log}$ if $y^{+}$ > $y_{lam}^{+}$

where,

ω = ω at $y^{+}$

$y^{+}$ = Estimated wall-normal distance of the cell centre in wall units

kqRWallFunction

The kqRWallFunction boundary condition provides a simple wrapper around the zero-gradient condition, which can be used for the turbulent kinetic energy, i.e., k, square-root of turbulent kinetic energy, i.e., q (e.g., in qZeta turbulence model), and Reynolds stress tensor fields, i.e., R (e.g., in LRR turbulence model), for the cases of high Reynolds number flow using wall function.

4. Model Setup

The mesh around the Ahmed Body was generated using the SnappyHexMesh method. Following Figure 4 demonstrates the expansion layers around Ahmed Body with the first layer thickness of 0.0000084 m. Overall, 10 expansion layers with expansion ratio were generated, and at a distance of 0.001 m from Ahmed Body, there is a refinement region with 4th level. Figure 6 demonstrates refinement boxes near the Ahmed Boy to develop the smooth transition from the boundary layer mesh to coarse fluid domain mesh. The nearest to Ahmed Body refinement box is refined up to the 3rd level, whereas middle and farthest ones are refined up to 2nd and 1st level, respectively. Moreover, the refinement boxes are longer in the direction of flow to capture the wake regions behind the Ahmed Body (Figure 5).

The dimensions of the Ahmed body are shown below (Figure 6):

(1)

W (x-axis) = 389 mm

(2)

L (y-axis) = 1040 mm

(3)

H (z-axis) = 288 mm

The boundary conditions of the model are as listed in

Table 1. Boundary conditions.

Boundary Conditions	k	nut	omega	p	u
Inlet	FixedValue; Uniform I;	Calculated; Uniform 0;	FixedValue; Uniform 10300;	zeroGradient;	FixedValue Uniform (40 0 0);
Outlet	zeroGradient;	Calculated; Uniform O·	zeroGradient;	FixedValue; Uniform O·	zeroGradient”
Left	zeroGradient;	Calculated; Uniform O;	zeroGradient;	zeroGradient;	Slip;
Right	zeroGradient;	Calculated; Uniform O;	zeroGradient;	zeroGradient;	Slip;
Top	zeroGradient;	Calculated; Uniform O;	zeroGradient;	zeroGradient;	Slip;
Ground	kqRWall Fun ction; Uniform I;	nutkWallFu nction; Uniform O;	omegaWallF unction; Uniform 10300·	zeroGradient;	FixedYalu e; Uniform(0 0 0)·
Ahmed	kqRWallFun ction; Uniform I;	nutkWal!Function; Uniform 0;	omegaWallF unction; Uniform 10300·	zeroGradient;	FixedValu e; Uniform (0 0 0)·

:

(1)

Inlet fluid velocity: 40 m/s;

(2)

Outlet: zero gradient velocity.

The top, ground, left, and the right side are assigned as walls. Figure 7 shows the boundary conditions discussed above and the Ahmed body inside it. The inlet velocity values were given at the initial time step; therefore, they were saved in the 0 folder in the OpenFOAM folder.

The box is 20 times larger along the Y-axis to fully develop and stabilize the fluid flow before meeting the Ahmed body. In addition to the above-mentioned walls and ground dimensions, other dimensional parameters are written in the blockMeshDict file for simulation.

The Fluid Properties are given based on the experimental conditions:

(1)

Density (ρ) of fluid (air) at 25 °C is 1.225 kg/m³;

(2)

Kinematic viscosity (υ) of fluid is 1.5 × 10⁻⁵ m²/s;

(3)

Reynolds number based on the height of Ahmed Body: $\frac{UH}{\vartheta }=\frac{40\ast 0.288}{1.5\ast {10}^{-5}}=768000.$

The above fluid properties and others, such as density and kinematic viscosity, were prescribed and stored in the OpenFOAM controlDict file.

5. Implementation of Vles in OpenFOAM

The VLES model was implemented in the OpenFOAM open-source CFD code by developing a new library for its resolution control factor F_r. Firstly, the k-ω SST turbulence model and base files were used to modify it for VLES implementation. This is because the VLES uses a resolution control factor F_r on top of the RANS/URANS model. As a result, the k-ω SST turbulence model was copied to a new folder, which was named newkOmegaSST. The k-ω SST turbulence model contains two files, which are kOmegaSST.H and kOmegaSST.C, and two base files with F1, F2 and other functions declared. Therefore, the same files were copied for the VLES implementation with a new name defined and declared as newkOmegaSST.C, newkOmegaSST.H, newkOmegaSSTBase.C and newkOmegaSSTBase.H. (Figure 8).

Secondly, the files and options directories from the Make folder and turbulentTransportModels.C file were copied to the newkOmegaSST folder. The purpose of Make files is to connect the library with the new turbulence model. Afterward, the wmake libso command was written to compile the library, which can be subsequently used as a new library for simulations.

The resolution control factor F_r was included in the newkOmegaSSTBase.C, as shown in Figure 9 below. Equation (17) (Resolution Control Function) was used to rescale the turbulence viscosity of the k-ω SST turbulence model, as well as the known coefficients, which were declared in the newkOmegaSSTBase.H file. As a result, the new VLES model was implemented in a modified source file of the k-ω SST turbulence model.

In order to use the compiled library with the Ahmed Body case, the controlDict file should be modified to connect the simulation with the VLES turbulence model. Therefore, the libs sub-dictionary with the new library assigned as libmyincompressibleTurbulenceModels.so was written in the controlDict file. Afterward, the turbulenceProperties file in the constant file was changed from kOmegaSST to newkOmegaSST. Thus, the VLES model was fully set up to solve the airflow around the Ahmed Body.

6. Results and Discussion

The mesh convergence study was performed with a drag coefficient as the monitoring parameter using the k-ω SST solver. The results are summarized in

Table 2. Mesh Convergence Study.

Mesh #	Mesh Cell Number	Dimension	Drag Coefficient	Difference %
1	1,786,575	0.14	0.4038	-
2	2,608,663	0.12	0.3693	8.54%
3	4,195,405	0.1	0.3508	4.53%
4	8,070,025	0.08	0.3451	1.65%

.

In order to generate the mesh refinement of the domain, the blockMesh file was modified at each step. For example, firstly, the mesh size of the fluid domain was 0.14 m in all directions; afterward, it was reduced to 0.12 m. As a result, the mesh cell number was doubled each time approximately. This was performed until the mesh convergence was completed; in other words, the simulation results are independent of the mesh size. Table 1 demonstrates the mesh convergence study with drag coefficient comparison. As can be seen from Table 2, the relative difference between the 1st and 2nd mesh is 8.54%; therefore, further meshing was required. Eventually, the difference between 4th and 3rd was calculated to be 1.65%, which shows that 4th mesh is the optimal one, and further mesh refinement is not required. This mesh was used to perform the URANS and VLES simulations, which were found to generate sufficient resolution for the wake structures. It should be noted that, generally speaking, one can only use time-averaged parameters to check mesh convergence and validate the results with measurements as the experimental measurements are all time-averaged, which means that RANS simulation is sufficient in the convergence study. Theoretically, it is also impossible to compare two 4D transient flows apart from their corresponding time-averaged parameters. Further mesh refinement will improve the resolution of the turbulent eddies, and the time-averaged quantities of flow variables remain the same.

The Y+ values near the walls were calculated to check the quality of the mesh, particularly the ability of the turbulence models to capture the viscous sublayers of the flow on the walls of the Ahmed Body, which are shown in Figure 10. The average Y+ value is 4.96, which is lower than the threshold value of 5 to capture the viscous sublayer. Moreover, as can be seen from Figure 10, the maximum values of Y+ are at the corners of the front and backside is less than 100, which shows that the boundary layers at the corners are not accurate. However, this will not affect the drag coefficient dramatically, as most of the surface mesh is still around 5 with a minimum value of 0.022.

The VLES steady-state simulation was tested before running the unsteady simulation. In other words, the steady-state VLES is similar to the RANS simulation, which ignores the transient behavior of the flow. Interestingly, the VLES steady simulation results show richer flow structures and higher resolution, as expected, compared with the RANS simulation, as demonstrated in Figure 11 and Figure 12. Therefore, the resolution control function of VLES demonstrated positive improvements over the RANS simulation. For comparison, Figure 11 and Figure 12 also show the URANS and VLES results in the forms of velocity and pressure contours, which demonstrate the progressively richer transient 4D turbulent flow structures from RANS, URANS to VLES as expected. The stagnation of flow occurs at the front of the Ahmed body, which results in the maximum pressure there. Moreover, as air flows over the front part of the Ahmed Body, the pressure decreases, which can be seen in the blue regions at the corners. Furthermore, over the slant, the turbulent boundary layer (shear layer) separation occurs; therefore, the associated pressure is illustrated graphically as a blue region, and the Kelvin–Helmholtz instability can be vividly captured partially by URANS and fully by VLES.

The vorticity field around the Ahmed body and near the slant was analyzed by implementing Q and vorticity functions. Figure 13 shows the vorticity contours by URANS and VLES, where the two counter-rotating vortices on the two slant edges are captured. Obviously, the VLES has a far better resolution of the turbulent eddies and their vorticity fields than the URANS counterpart. Moreover, over the slant, the flow separates and then reattaches to the surface in an unsteady manner, similar to the dynamic stall phenomenon for airfoils. The front side of Ahmed Body also experiences flow separation at the low corner, and the recirculating flow then emanates from the corner and spreads out to the two sides of the car body, as can be seen from Figure 13c,d, but the URANS results fail to show these phenomena in Figure 13a,b. The implication is that one cannot use URANS for noise analysis for the design of car bodies, although it can generate sufficiently accurate results for the largest flow structures. It can also be noted that two rotating vortices blocks flow originating from the separation zone in front of Ahmed Body, which increases the pressure difference at the front and back side, in turn, increases the energy consumption and results in higher pressure drag.

Figure 14 illustrates the flow field streamlines at the back of the Ahmed Body by URANS and VLES, which illustrates the recirculation zones at the rear end after slant angle surface, as well as the counter-rotating vortices over the two sides of the slant. The rear recirculation zone happens due to the separation of flow at the rear, which develops the low fluid energy zone at the back of the Ahmed Body. Moreover, the VLES simulation generates far richer flow structures here for the complex 3D turbulent separated flow around the Ahmed Body and is confirmed by the experimental observation, which is shown in Figure 2.

The validation of VLES simulations by comparison with the experimental drag coefficient was also performed. Here the drag coefficients of different models, including RANS (k-ω SST), URANS, VLES, DES and LES, were compared with the experimental measurement (0.298) in

Table 3. Drag coefficients by RANS, URANS, VLES, DES and LES and their comparison with experimental data.

Turbulence Models	RANS (k-ω SST)	URANS	VLES	DES-SST [10]	LES-SVV [10]
Drag coefficient	0.345	0.329	0.3127	0.343	0.431
Experimental data	0.298	0.298	0.298	0.298	0.298
Error (%)	14.8	10.4	4.93	15.1	44.6

. As with the observations by Igali et al. [5] and Sadykov et al. [8], the performance of RANS is comparable to those of DES and LES, with DES being the worst inaccuracy, and URANS better than the three. Finally, the VLES is the best in accuracy among the five models.

Figure 15 and Figure 16 represent the predicted drag coefficients as functions of time by the k-ω SST URANS and VLES models, respectively. The fluctuations of drag coefficients develop due to the transient nature of the flow field because of shear layer instability, unsteady flow separation and flow recirculation. Therefore, as can be seen from Figure 15, the time-averaged drag coefficient by URANS was approximately 0.329, which differs from experimental data by 10.4%, as shown in Table 3. The VLES simulation was performed, and the predicted drag coefficient as a function of time was plotted in Figure 16. The average drag coefficient was approximately 0.3127, which was 4.93% higher than the experimental data. It is thus further confirmed that the VLES turbulence model is found to be the best among RANS, URANS, VLES, DES and LES models.

Figure 17 below plots the X-velocity profiles at different locations over the slant surface by different turbulence models. In addition to the simulated results, the experimental velocity profiles for 25-degree slant angle were also plotted in the graph for validation of the models [26]. All data were obtained along the centerline (of the Z-axis) of the Ahmed Body. The X-velocity component was evaluated in the five different distances of the X-axis starting from the tail of the model every 0.06 m in negative X-direction.

The graph shows the agreement of the simulated velocity profiles with the experimental ones, where the best match is observed to be the results of the VLES model. Moreover, it is noted that the velocity profiles between RANS k-ω SST and k-ε models were very similar. Nevertheless, the profiles predicted by the RANS k-ω SST model were closer to the experimental data than the k-ε model in the separation zone, as expected. The notable point is seen in the first two results, where the experimental profiles did not reach the surface of the Ahmed Body completely. Hence, the simulation results extend beyond the experimental results [27].

For the sake of comparing the computational efficiencies of the URANS, VLES, DES and LES models, their CPU times used for the Ahmed body simulations are shown in

Table 4. Comparison of CPU Times of URANS, VLES, DES and LES.

	Time Step Size (s)	Time Period (s)	CPU Time (h)
URANS	3 × 10−5	0.39	175 DELL 4-Core PC
VLES	1.1 × 10−4	0.39	675 DELL 4-Core PC
DES-SST [10]	2.57 × 10−4	0.92	20000 IBM P4 (64 procs)
LES-SVV [10]	6 × 10−5	0.29	500 NEC SX8 (8 procs)

. While the VLES time step size is greater than the URANS’, it was found that VLES needs more iterations to solve velocity, pressure, as well as k and omega per time step as greater fluctuations in the flow field requires more interactions in the PISO scheme to ensure solution convergence. Thus, the VLES solver is four times slower than the URANS one, but it generates a much greater resolution of the turbulent eddies and vorticity field. The VLES was found to be much more efficient than and DES and LES, given that the LES was run on a 4-core PC while the DES and LES were run on parallel supercomputers.

7. Conclusions

This paper presents the development of an AHTM approach based on VLES and its implementation DAFoam/OpenFOAM for the advanced study of aerodynamic flow over an idealized car model (the Ahmed Body) and the evaluation of its potential as an accurate and efficient solver in an ADO tool such as DAFoam. The Ahmed Body with a slant angle of 25 degrees was studied with RANS, URANS and VLES with an extra effort on the study of the Kelvin–Helmholtz instability effect, and their performances were compared among themselves and with those of DES and LES models as well. Firstly, a steady RANS simulation was performed to predict the turbulent flow around the Ahmed Body and calculate its drag coefficient. The URANS and VLES models were then employed to simulate the unsteady turbulent flow to study the Kelvin–Helmholtz instability effect in the shear layers over the body. The VLES was observed to resolve the turbulent structures well with the visualization of the flow field, in particular the vorticity field. The VLES was found to be the best model among the five models examined in this study and by careful validation with experimental measurements. As a result, it can be concluded that the main objective of this study was successfully achieved, which was to implement the VLES model in DAFoam/OpenFOAM for better capturing the counter-rotating vortices, Kelvin–Helmholtz instability and the separation of flow at the back of the Ahmed body. This shows that the AHMT approach based on VLES offers the automotive designer a highly flexible yet accurate and efficient tool for ADO for drag and noise reduction and drivability enhancement for cars based on gradient-based efficient adjoint solver and Sparse Nonlinear optimizer (SNOPT) such as the DAFoam.