Steady turbulent flow over a backward-facing step

2.1 Steady turbulent flow over a backward-facing step

In this example we shall investigate steady turbulent flow over a backward-facing step. The problem description is taken from one used by Pitz and Daily in an experimental investigation [**] against which the computed solution can be compared. This example introduces the following OpenFOAM features for the first time:

  • generation of a mesh using blockMesh using full mesh grading capability;
  • steady turbulent flow.

2.1.1 Problem specification

The problem is defined as follows:

Solution domain
The domain is 2 dimensional, consisting of a short inlet, a backward-facing step and converging nozzle at outlet as shown in Figure 2.1.
\relax \special {t4ht=

Figure 2.1: Geometry of backward-facing step


Governing equations
 
  • Mass continuity for incompressible flow
    ∇ ∙U =  0
                                                                          

                                                                          
           \relax \special {t4ht=
    (2.1)

  • Steady flow momentum equation
      ∙          ∙
∇  (UU ) + ∇  R =  - ∇p
           \relax \special {t4ht=
    (2.2)

    where p  \relax \special {t4ht= is kinematic pressure and (in slightly over-simplistic terms) R =  ν   ∇U
      eff  \relax \special {t4ht= is the viscous stress term with an effective kinematic viscosity νeff  \relax \special {t4ht=, calculated from selected transport and turbulence models.

Initial conditions
U  = 0 m ∕s  \relax \special {t4ht=, p = 0 Pa  \relax \special {t4ht= — required in OpenFOAM input files but not necessary for the solution since the problem is steady-state.
Boundary conditions
 
  • Inlet (left) with fixed velocity U  = (10, 0,0)  \relax \special {t4ht= m/s;
  • Outlet (right) with fixed pressure p = 0 Pa  \relax \special {t4ht=;
  • No-slip walls on other boundaries.
Transport properties
 
  • Kinematic viscosity of air ν = μ∕ρ =  18.1 × 10 -6∕1.293 = 14.0 μm2 ∕s  \relax \special {t4ht=
Turbulence model
 
  • Standard k -  ϵ  \relax \special {t4ht=;
  • Coefficients: Cμ = 0.09;C1  = 1.44;C2 = 1.92;αk =  1;αϵ = 0.76923  \relax \special {t4ht=.
Solver name
simpleFoam: an implementation for steady incompressible flow.
Case name
pitzDaily, located in the $FOAM_TUTORIALS/incompressible/simpleFoam directory.

The problem is solved using simpleFoam, so-called as it is an implementation for steady flow using the SIMPLE algorithm. The solver has full access to all the turbulence models in the incompressibleTurbulenceModels library and the non-Newtonian models incompressibleTransportModels library of the standard OpenFOAM release.

2.1.2 Mesh generation

We expect that the flow in this problem is reasonably complex and an optimum solution will require grading of the mesh. In general, the regions of highest shear are particularly critical, requiring a finer mesh than in the regions of low shear. We can anticipate where high shear will occur by considering what the solution might be in advance of any calculation. At the inlet we have strong uniform flow in the x  \relax \special {t4ht= direction and, as it passes over the step, it generates shear on the fluid below, generating a vortex in the bottom half of the domain. The regions of high shear will therefore be close to the centreline of the domain and close to the walls.

The domain is subdivided into 12 blocks as shown in Figure 2.2.


\relax \special {t4ht=


Figure 2.2: Blocks in backward-facing step


The mesh is 3 dimensional, as always in OpenFOAM, so in Figure 2.2 we are viewing the back plane along z = - 0.5  \relax \special {t4ht=. The full set of vertices and blocks are given in the mesh description file below:

1  /*--------------------------------*- C++ -*----------------------------------*\
2  | =========                 |                                                 |
3  | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
4  |  \\    /   O peration     | Version:  plus                                  |
5  |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
6  |    \\/     M anipulation  |                                                 |
7  \*---------------------------------------------------------------------------*/
8  FoamFile
9  {
10      version     2.0;
11      format      ascii;
12      class       dictionary;
13      object      blockMeshDict;
14  }
15  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16  
17  convertToMeters 0.001;
18  
19  vertices
20  (
21      (-20.6 0 -0.5)
22      (-20.6 3 -0.5)
23      (-20.6 12.7 -0.5)
24      (-20.6 25.4 -0.5)
25      (0 -25.4 -0.5)
26      (0 -5 -0.5)
27      (0 0 -0.5)
28      (0 3 -0.5)
29      (0 12.7 -0.5)
30      (0 25.4 -0.5)
31      (206 -25.4 -0.5)
32      (206 -8.5 -0.5)
33      (206 0 -0.5)
34      (206 6.5 -0.5)
35      (206 17 -0.5)
36      (206 25.4 -0.5)
37      (290 -16.6 -0.5)
38      (290 -6.3 -0.5)
39      (290 0 -0.5)
40      (290 4.5 -0.5)
41      (290 11 -0.5)
42      (290 16.6 -0.5)
43      (-20.6 0 0.5)
44      (-20.6 3 0.5)
45      (-20.6 12.7 0.5)
46      (-20.6 25.4 0.5)
47      (0 -25.4 0.5)
48      (0 -5 0.5)
49      (0 0 0.5)
50      (0 3 0.5)
51      (0 12.7 0.5)
52      (0 25.4 0.5)
53      (206 -25.4 0.5)
54      (206 -8.5 0.5)
55      (206 0 0.5)
56      (206 6.5 0.5)
57      (206 17 0.5)
58      (206 25.4 0.5)
59      (290 -16.6 0.5)
60      (290 -6.3 0.5)
61      (290 0 0.5)
62      (290 4.5 0.5)
63      (290 11 0.5)
64      (290 16.6 0.5)
65  );
66  
67  blocks
68  (
69      hex (0 6 7 1 22 28 29 23) (18 7 1) simpleGrading (0.5 1.8 1)
70      hex (1 7 8 2 23 29 30 24) (18 10 1) simpleGrading (0.5 4 1)
71      hex (2 8 9 3 24 30 31 25) (18 13 1) simpleGrading (0.5 0.25 1)
72      hex (4 10 11 5 26 32 33 27) (180 18 1) simpleGrading (4 1 1)
73      hex (5 11 12 6 27 33 34 28) (180 9 1) edgeGrading (4 4 4 4 0.5 1 1 0.5 1 1 1 1)
74      hex (6 12 13 7 28 34 35 29) (180 7 1) edgeGrading (4 4 4 4 1.8 1 1 1.8 1 1 1 1)
75      hex (7 13 14 8 29 35 36 30) (180 10 1) edgeGrading (4 4 4 4 4 1 1 4 1 1 1 1)
76      hex (8 14 15 9 30 36 37 31) (180 13 1) simpleGrading (4 0.25 1)
77      hex (10 16 17 11 32 38 39 33) (25 18 1) simpleGrading (2.5 1 1)
78      hex (11 17 18 12 33 39 40 34) (25 9 1) simpleGrading (2.5 1 1)
79      hex (12 18 19 13 34 40 41 35) (25 7 1) simpleGrading (2.5 1 1)
80      hex (13 19 20 14 35 41 42 36) (25 10 1) simpleGrading (2.5 1 1)
81      hex (14 20 21 15 36 42 43 37) (25 13 1) simpleGrading (2.5 0.25 1)
82  );
83  
84  edges
85  (
86  );
87  
88  boundary
89  (
90      inlet
91      {
92          type patch;
93          faces
94          (
95              (0 22 23 1)
96              (1 23 24 2)
97              (2 24 25 3)
98          );
99      }
100      outlet
101      {
102          type patch;
103          faces
104          (
105              (16 17 39 38)
106              (17 18 40 39)
107              (18 19 41 40)
108              (19 20 42 41)
109              (20 21 43 42)
110          );
111      }
112      upperWall
113      {
114          type wall;
115          faces
116          (
117              (3 25 31 9)
118              (9 31 37 15)
119              (15 37 43 21)
120          );
121      }
122      lowerWall
123      {
124          type wall;
125          faces
126          (
127              (0 6 28 22)
128              (6 5 27 28)
129              (5 4 26 27)
130              (4 10 32 26)
131              (10 16 38 32)
132          );
133      }
134      frontAndBack
135      {
136          type empty;
137          faces
138          (
139              (22 28 29 23)
140              (23 29 30 24)
141              (24 30 31 25)
142              (26 32 33 27)
143              (27 33 34 28)
144              (28 34 35 29)
145              (29 35 36 30)
146              (30 36 37 31)
147              (32 38 39 33)
148              (33 39 40 34)
149              (34 40 41 35)
150              (35 41 42 36)
151              (36 42 43 37)
152              (0 1 7 6)
153              (1 2 8 7)
154              (2 3 9 8)
155              (4 5 11 10)
156              (5 6 12 11)
157              (6 7 13 12)
158              (7 8 14 13)
159              (8 9 15 14)
160              (10 11 17 16)
161              (11 12 18 17)
162              (12 13 19 18)
163              (13 14 20 19)
164              (14 15 21 20)
165          );
166      }
167  );
168  
169  mergePatchPairs
170  (
171  );
172  
173  // ************************************************************************* //

A major feature of this problem is the use of the full mesh grading capability of blockMesh that is described in section ?? of the User Guide. The user can see that blocks 4,5 and 6 use the full list of 12 expansion ratios. The expansion ratios correspond to each edge of the block, the first 4 to the edges aligned in the local x1   \relax \special {t4ht= direction, the second 4 to the edges in the local x2   \relax \special {t4ht= direction and the last 4 to the edges in the local x3   \relax \special {t4ht= direction. In blocks 4, 5, and 6, the ratios are equal for all edges in the local x1   \relax \special {t4ht= and x3   \relax \special {t4ht= directions but not for the edges in the x2   \relax \special {t4ht= direction that corresponds in all blocks to the global y  \relax \special {t4ht=. If we consider the ratios used in relation to the block definition in section ?? of the User Guide, we realize that different gradings have been prescribed along the left and right edges in blocks 4,5 and 6 in Figure 2.2. The purpose of this differential grading is to generate a fine mesh close to the most critical region of flow, the corner of the step, and allow it to expand into the rest of the domain.

The mesh can be generated using blockMesh from the command line and viewed as described in previous examples.

2.1.3 Boundary conditions and initial fields

Edit the case files to set the initial and boundary fields for velocity U  \relax \special {t4ht=, pressure p  \relax \special {t4ht=, turbulent kinetic energy k  \relax \special {t4ht= and dissipation rate ε  \relax \special {t4ht=. The boundary conditions can be specified as: the upper and lower walls are set to Wall, the left patch to Inlet and the right patch to Outlet. These physical boundary conditions require us to specify a fixedValue at the inlet on U  \relax \special {t4ht=, k  \relax \special {t4ht= and ε  \relax \special {t4ht=. U  \relax \special {t4ht= is given in the problem specification, but the values of k  \relax \special {t4ht= and ϵ  \relax \special {t4ht= must be chosen by the user in a similar manner to that described in section 1.1.8.1 of the User Guide. We assume that the inlet turbulence is isotropic and estimate the fluctuations to be 5%  \relax \special {t4ht= of U  \relax \special {t4ht= at the inlet. We have

 ′     ′    ′   -5--
Ux = U y = Uz = 100 10 = 0.5 m ∕s
\relax \special {t4ht=
(2.3)

and

     3     2          2  2
k =  2(0.5) = 0.375 m  ∕s
\relax \special {t4ht=
(2.4)

If we estimate the turbulent length scale l  \relax \special {t4ht= to be 10%  \relax \special {t4ht= of the width of the inlet then

    C0.75k1.5        0.75     1.5
ε = --μ------=  -0.09---0.375---- = 14.855m2 ∕s3
        l       0.1 × 25.4 × 10 -3
\relax \special {t4ht=
(2.5)

At the outlet we need only specify the pressure p = 0Pa  \relax \special {t4ht=.

2.1.4 Case control

The choices of fvSchemes are as follows: the timeScheme should be steadyState; the gradSchemes and laplacianSchemes should be set as default to Gauss; and, the divSchemes should be set to upwind to ensure boundedness.

Special attention should be paid to the solver settings of the fvSolution dictionary. Although the top level simpleFoam code contains only equations for p  \relax \special {t4ht= and U  \relax \special {t4ht=, the turbulence model solves equations for k  \relax \special {t4ht=, ε  \relax \special {t4ht= and R  \relax \special {t4ht=, and tolerance settings are required for all 5 equations. A tolerance of 10-5   \relax \special {t4ht= and relTol of 0.1 are acceptable for all variables with the exception of p  \relax \special {t4ht= when 10- 6   \relax \special {t4ht= and 0.01 are recommended. Under-relaxation of the solution is required since the problem is steady. A relaxationFactor of 0.7 is acceptable for U  \relax \special {t4ht=, k  \relax \special {t4ht=, and ε  \relax \special {t4ht= but 0.3 is required for p  \relax \special {t4ht= to avoid numerical instability.

Finally, in the controlDict dictionary, the time step deltaT should be set to 1 since in steady state cases such as this is effectively an iteration counter. With benefit of hindsight we know that the solution requires 1000 iterations reach reasonable convergence, hence endTime is set to 1000. Ensure that the writeInterval is sufficiently high, e.g. 50, that you will not fill the hard disk with data during run time.

2.1.5 Running the case and post-processing


PIC\relax \special {t4ht=
(a) Velocity vectors after 50 iterations
PIC\relax \special {t4ht=
(b) Velocity vectors at 1000 iterations
PIC\relax \special {t4ht=
(c) Streamlines at 1000 iterations

Figure 2.3: Development of a vortex in the backward-facing step.


Run the case and post-process the results. After a few iterations, e.g. 50, a vortex develops beneath the corner of the step that is the height of the step but narrow in the x  \relax \special {t4ht=-direction as shown by the vector plot of velocities is shown Figure 2.3(a). Over several iterations the vortex stretches in the x  \relax \special {t4ht=-direction from the step to the outlet until at 1000 iterations the system reaches a steady-state in which the vortex is fully developed as shown in Figure 2.3(b-c).