Airflow over a car

In this example we shall study an airflow over a car. Such a kind of simulations are used in automotive industry to modelize the shape of cars. In this example we shall assume the flow to be steady turbulent.

1 Problem specification

We want to study the motion of air circulating over a car. The experimental set up for studying this problem is a wind tunnel. Instead of considering the full spatial domain the simulation can be carried out over a smaller domain by observing that the problem admits a symmetry, as shown in Fig.1 where is outlined the geometry. The air enters perpendicularly to the inlet at a velocity U = 27.8 m/s. The operating pressure is p = 1 atm = 101325 Pa.


Figure 1: The geometry.

1.1 Governing equations

As the case is steady-state and the fluid is supposed viscous and incompressible, the governing equations to be solved are:
  • Mass conservation for incompressible flow: ∇ ⋅ U = 0
  • Steady-flow momentum conservation: (U ⋅ ∇) U + (∇P/ρ) = ν ∇2U + g
where P, ρ, ν and g are, respectively, presuure, density, kinematic viscosity and gravitational field. We have already talked about the initial velocity above. The initial pressure has been set equal to 1 atm everywhere except at the boundaries where a zero value for the pressure gradient is required.


1.2 Turbulence model

We choose the standard k - ε turbulence model with coefficients Cμ = 0.09, C1 = 1.44, C2 = 1.92, α k = 1, α ε = 7.69. In this model the velocity of the turbulent flow is estimated as a fluctuation around an average value: U = ‹ U › + U'; k = (1/2) |U'|2 represents the turbulent kinetic energy per unit of mass and ε the energy dissipation rate. We set the a value of k = 0.02898 m2/s2 at the inlet, the same value has been set at the outlet for modeling the turbulence in case of backflow. The value of ε is set to ε = 1.8239 m2/s3.

2 Mesh generation

The mesh (about 3600000 cells) was generated by converting a .msh GAMBIT file according to the procedure described in the section Mesh conversion.

Figure 2: Top: surfaces of the car geometry obtained by converting the original GAMBIT mesh file. On the left is shown the external surface, on the right the internal surface. Bottom: on the left is shown the body car mesh, on the right the unstructured mesh generated on the central symmetry plane. The mesh gets finer near the car to capture boundary layer effects.

3 Boundary conditions and initial fields

The case requires initial and boundary conditions settings for all the involved fields: velocity U, pressure P, turbulent kinetic energy k and energy dissipation rate ε. The initial values for the fields are specified in the casename/0/fieldname files. The reader may consult the boundary conditions documentation section (in Italian) for more details. We report the rilevant sections of the fieldname files.



The initial velocity U has been set equal to 27.8 m/s in the y direction at the inlet:

dimensions [0 1 -1 0 0 0 0]; // dimension in m s -1
internalField uniform (0 0 0);
boundaryField
{
inlet // inlet patch
{
type fixedValue;
value uniform (0 27.8 0); // U x = U z = 0, U y = 27.8 m/s
}


The initial turbulent kinetic energy has been set equal to k = 0.02898 m2/s2 at the inlet and outlet (for modeling back flow turbulence):

dimensions [0 2 -2 0 0 0 0]; // dimension in m 2 s -2
internalField uniform 0.02898;
boundaryField
{ inlet // inlet patch
{
type fixedValue;
value uniform 0.02898;
}
outlet // inlet patch
{
type fixedValue;
value uniform 0.02898;
}


The initial dissipation energy rate has been set equal to ε = 1.8239 m2/s3 at the inlet and outlet:

dimensions [0 2 -3 0 0 0 0]; // dimension in m 2 s -3
internalField uniform 1.8239;
boundaryField
{
inlet
{
type fixedValue;
value uniform 1.8239;
}
outlet
{
type fixedValue;
value uniform 1.8239;
}


The initial pressure field P has been set equal to 1 atm everywhere except at the boundaries where a smooth shape condition is required. :

dimensions [0 2 -2 0 0 0 0]; dimension in m 2 s -2
internalField uniform 101325;
boundaryField
{
inlet // inlet type patch (extremity of the conduct)
{
type zeroGradient; // smooth condition
}
..............
}

4 Solver and case control

We used the simpleFoam solver, implemented for steady incompressible flow. The numerical schemes for terms such as derivatives, interpolation procedures etc. are set in the casename/system/fvSchemes file. In this file we specify that the case is steady-state by giving the SteadyState value to the keyword timeScheme. The terms ∇ and ∇ 2 controlled respectively by the keywords gradSchemes and laplacianSchemes have been set to the value Gauss because we adopted the standard finite volume discretization of Gaussian integration which requires the interpolation of values from cell centres to face centres. The divergence term ∇ ⋅, controlled by the keyword divSchemes has been set to UD to ensure boundedness.

5 Results

In Fig. 3 are shown the residuals of the involved fields after 10000 iterations. The case has been run in parallel using 32 cores.
Fig. 3: Residuals of the velocity components, the pressure, the turbulent kinetic energy and the dissipation energy rate.

The same simulation has been run with Fluent. In order to make a comparison we shall analyse the shape of the fields on the central symmetry plane.



Figure 4: Velocity magnitude contour plot for Fluent (top) and OpenFOAM (bottom) on the central symmetry plane.

In Fig. 4 we observe that the magnitude of the velocity fields is of the same order. The velocity field around the car has a similar shape in both the two simulations. In particular, we observe a similar behaviour in the region at the rear of the car where turbulence effects are present. In Fig.5 is shown in more detail the vector velocity flow near the surface of the car.


Figure 5: Velocity vector flow near the surface of the car for Fluent (top) and OpenFOAM (bottom) on the central symmetry plane.

In Fig. 6 is shown the pressure contour plot for Fluent and OpenFOAM on the central symmetry plane. The Fluent picture reports the pressure values (in Pa) referred to 1 atm, thus negative values indicate a pressure lower than 1 atm. On the other hand the OpenFOAM picture reports the absolute pressure values as displayed in the scale. Thus in this case an overall agreement is observed, too.


Figure 6: Pressure contour plot for Fluent (top) and OpenFOAM (bottom) on the central symmetry plane.

File LaTeX documentation..