OpenFOAM 3D Steady Turbulent Simulation (OF-1.4.1)

In this example we shall investigate how to make a 3D simulation with OpenFOAM. The case we want to study is a steady turbulent flow circulating inside a room. This example introduces almost all the basic steps which should be done for studying a CFD problem and explains how these things can be carried out in a quite simple way in OpenFOAM.

1 Problem specification

We want to study the motion of air circulating inside a room. The outline of the geometry is shown in Fig. 1. The room is crossed by a conduct. The air enters perpendicularly to one of the extremities of the conduct (inlet) at a velocity U = 0.5 m/s and flows out the other (outlet). There are two holes in the conduct and four small windows on one of the room walls through which the air can also flow. Two pillars are placed between the conduct and the wall with the windows.


Figure 1: The geometry of the room.

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 = 101325 Pa at all the boundaries except at the windows where its value is 0.995 atm to allow a pressure gradient.


1.2 Turbulence model

We choosed 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 assume the fluctuations at the inlet to be isotropic and 5% of U, thus Ux = Uy = Uz = 0.5 ⋅ 0.05 = 25 ⋅ 10-3 m/s and k = 9.375 ⋅ 10-4 m2/s2. The energy dissipation is related to the turbulence scale l which we estimated to be 10% of 0.5 m, which is the width of the inlet. The value of ε is thus set to ε = (C μ 0.75 ⋅ k 3/2 )/l = 9.53 ⋅ 10 -5 m2/s3.

2 Mesh generation

The mesh was generated by converting a .msh GAMBIT file according to the procedure described in the section Mesh conversion. The domain consists of about 5000000 cells. In Fig. 2 are shown some of the meshed patches of the domain.

Figure 2: Structured patch meshes of the conduct, room walls, inlet and pillars obtained by using the OpenFOAM conversion routines.

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 0.5 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 0.5 0); // U x = U z = 0, U y = 0.5 m/s
}


The initial turbulent kinetic energy has been set equal to k = 9.375 ⋅ 10-4 m2/s2 at the inlet:

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


The initial dissipation energy rate has been set equal to ε = 9.53 ⋅ 10 -5 m2/s3 at the inlet:

dimensions [0 2 -3 0 0 0 0]; // dimension in m 2 s -3
internalField uniform 9.53e-05;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.0016875;
}


The initial pressure field P has been set equal to 1 atm at the inlet and at all the outlets except at the windows where its value is 0.995 atm to allow internal flow circulation. The gradient of the pressure in the direction perpendicular to the the patches of type wall has been set equal to 0:

dimensions [0 2 -2 0 0 0 0]; dimension in m 2 s -2
internalField uniform 101325;
boundaryField
{
wall_ceiling // wall type patch (the ceiling)
{
type zeroGradient;
}
inlet // inlet type patch (extremity of the conduct)
{
type fixedValue;
value uniform 101325;
}
outlet_conduct // outlet type patch (extremity of the conduct)
{
type fixedValue;
value uniform 1;
}
outlets_windows // outlet type patch (the windows)
{
type fixedValue;
value uniform 108800;
}

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 1000 iterations. The case has been run in parallel using 64 cores.

Fig. 3: Residuals of the velocity components, the pressure, the turbulent kinetic energy and the dissipation energy rate.

Even if the maximum pressure gradient is only 0.5%, it is sufficient to cause a back flow with the air entering the room from the outlet as shown in Fig. 4.



Figure 4: Velocity flow on a horizontal plane parallel to the floor. On the left is shown the vertical view looking at the floor from the ceiling of the room. On the right is shown in detail the region near the outlet. The pressure gradient forces the air to enter the room from the outlet, where the velocity reaches its maximum value. The air moves through the conduct and flows towards the wall windows, where the pressure has the lowest value.

The same simulation has been run with Fluent. In order to make a comparison we shall analyse the shape of the fields on two cutting planes, as shown in Fig. 5.

Figure 5: The figure shows the two cutting planes chosen to compare the results of the simulation for OpenFOAM and Fluent.


Figure 6: Velocity magnitude contour plot for Fluent (left) and OpenFOAM (right) on the horizontal cutting plane.

We observe that the magnitude of the velocity fields is of the same order. In the regions far from the two small holes placed on the conduct the velocity magnitude approaches to zero in both the two simulations. In the regions next to the conduct holes the velocity contour has a more diffusive shape in the OpenFOAM case. In Fig. 7 are shown the contour plots for the velocity components. An overall agreement is observed.

Fig. 7: Contour plot comparisons for the velocity components on the horizontal plane.

In Fig. 8 are shown the contour plots for the magnitude velocity field on the vertical plane. On average, the strength of the velocity fields at the conduct section are comparable. There is a nice agreement even for the size of the region where the velocity field doest not dump to zero.


Figure 8: Velocity magnitude contour plot for Fluent (top) and OpenFOAM (bottom) on the vertical cutting plane.

In Fig. 8 are shown the pressure contour plots on the horizontal plane. The Fluent picture reports (in Pa) the pressure gradient with respect to 1 atm = 101325 Pa as reference value, thus the negative values mean a value of pressure lower tham 1 atm. On the other hand the OpenFOAM picture reports (in Pa) the absolute pressure values as shown in the dispayed range. Hence in this case too there is an overall agreement.


Figure 8: Pressure contour plots for Fluent (left) and OpenFOAM (right) on the horizontal cutting plane.

Download the test case.