Finite Element Analysis Toolbox
Turbulent Flow Over a Backwards Facing Step (OpenFOAM)

Flow over backwards facing steps are classic computational fluid dynamics test problems, and are often used for validation of simulation codes. This test case consists of studying how the flow field reacts to a sudden expansion in a channel. The expansion will cause a break in the flow and a recirculation or separation zone will form. To validate the results the computed length of the recirculation zone is compared with the experimental results of Pitz and Daily [1].

In this example the stationary incompressible Navier-Stokes equations are used to model the fluid with simulation parameters corresponding to a Reynolds number, Re = 18000. The flow is therefore fully turbulent, whereby a turbulence model closure must also be applied. Here the standard k-epsilon turbulence model is used which are available with both the OpenFOAM and SU2 CFD solver integrations built into FEATool Multiphysics.

Tutorial

This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Fluid Dynamics > Turbulent Flow Over a Backwards Facing Step from the File menu, viewed as a video tutorial, or alternatively, follow the step-by-step instructions below. Note that the CFDTool interface differ slightly from the FEATool Multiphysics instructions described in the following. (Note that the external OpenFOAM solver needs to be installed for this model to run correctly.)

  1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.

In this example variations in the z-direction are assumed to be negligible which enables using a planar 2D approximation, and saving computational cost compared to a full 3D simulation.

  1. Press the 2D radio button and select the Navier-Stokes Equations physics mode from the Select Physics drop-down menu. (Note that for CFDTool the physics selection is done in the Equation settings dialog box.)

  2. Press OK to finish the physics mode selection.

The backwards facing step geometry features a slightly tapered outflow region, which is easiest to create by directly specifying a polygon by coordinates.

  1. Select Create Object > Polygon from the Geometry menu, and enter the following data into the Point coordinates table.
1 2 3 4 5 6 7 8
x -0.0206 0 0 0.206 0.29 0.29 0.206 -0.0206
y 0 0 -0.0254 -0.0254 -0.0166 0.0166 0.0254 0.0254


  1. Press OK to finish and close the dialog box. The geometry should now look like the following.

The dimensions are here given in meters. However, note that the toolbox works with any unit system, and it is up to the user to choose consistent units for geometry dimensions, material, equation, and boundary coefficients.

  1. Switch to Grid mode by clicking on the corresponding Mode Toolbar button.

For turbulent flows it is particularly important to ensure a fine resolution near walls, due to viscous boundary layers which feature steep gradients.

  1. The default grid will be too coarse to ensure an accurate solution. To create a finer grid, enter 0.002 into the Grid Size edit field, and press the Generate button.

  2. Switch to Equation mode by clicking on the corresponding Mode Toolbar button.

  3. Since the fluid is air, set the density ρ to 1.293 kg/m3 and viscosity µ to 18.1e-6 kg/m/s, in the corresponding edit fields of the Equation Settings dialog box Press OK to finish and close the dialog box. (Note that the Equation Settings dialog box may look different for CFDTool.)

  4. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.
  5. In the Boundary Settings dialog box, first select all boundaries except for the right outflow and left inflow (wall boundaries 1-4, and 6-7) in the left hand side Boundaries selection list box, and choose the Wall/no-slip boundary condition from the boundary condition drop-down menu.

  6. Then select the leftmost boundary (number 8) in the Boundaries list box, and choose the Inlet/velocity boundary condition from the drop-down menu. Enter 10 m/s in the edit field for the x-velocity coefficient u0.

  7. Finally select the right outflow boundary (number 5), and select the Outflow/pressure boundary condition from the drop-down menu. Finish the boundary condition specification by clicking the OK button.

  8. Now that the problem has been defined, press the Solve Mode Toolbar button to switch to solve mode.

Instead of using the built-in multiphysics solver, the external and dedicated CFD solver OpenFOAM will be used to solve this flow problem. (See the corresponding OpenFOAM solver section in the User's Guide on how to install OpenFOAM.)

  1. Press the OpenFOAM Toolbar button to open the OpenFOAM solver settings and control panel dialog box. and set the Stopping criteria/tolerance for initial residuals to 1e-4.

  2. Select the k-epsilon Turbulence model from the corresponding drop-down menu, and press the Edit button to open a dialog box where the turbulence inlet quantities can be prescribed.

The k and epsilon inlet values can either be estimated from a given turbulence intensity and length scale, or as is done here prescribed directly if the quantities are known.

  1. Enter 0.375 into the edit field for the Turbulent kinetic energy, 14.855 into the edit field for the Turbulent dissipation rate, and press OK.

Advanced users can also use the Edit functionality, to view, fine tune, and modify OpenFOAM dictionaries, as well as export case files for running OpenFOAM simulations separately.

  1. Press the Solve button to start the OpenFOAM solver. The solver control panel will automatically change to show the convergence process for the flow variables.

During the solution process one can switch between the convergence tab and the solver output log. The solver will stop when the residuals for all of the flow variables are below the stopping criteria (here 1e-4), or the maximum number of iterations has been reached.

After the solver has finished, the toolbox will automatically switch to postprocessing mode, and show the resulting velocity field. One can see that a recirculation zone has formed after the step.

  1. To visualize the recirculation zone more clearly, open the Postprocessing settings dialog box, and enter the expression for the normalized recirculation zone length (u<0)*x/25.4e-3 into the Surface Plot expression edit field. The Arrow Plot option can also be turned on to help visualize the flow field.

Note that logical or switch type expressions such as a<b evaluate to either 0 or 1, and are here used to limit the plot to the lower half region for which the u velocity is negative.

The resulting plot shows a recirculation zone length of about 6.25 normalized length units, which is quite close to the reference length of 6.4.

The turbulent flow over a backwards facing step fluid dynamics model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB m-script text file, or GUI script (.fes) file.

CLI Postprocessing

To visualize the full 3D solution from the planar model, the data can be exported and processed on the MATLAB command line interface (CLI) console with the Export Model Data Struct to MATLAB option from the File menu. The postextrude and postplot functions can then be applied to extrude and visualize the data, for example 

fea_extruded = postextrude( fea, 5, 50e-3 );

postplot( fea_extruded, 'surfexpr', 'sqrt(u^2+v^2)', ...
          'parent', figure, 'axis', 'off', 'colorbar', 'off' )
view(25, 25)

Video

References

[1] Pitz RW, Daily JW. Combustion in a Turbulent Mixing Layer Formed at a Rearward Facing Step. AIAA Journal 21, 1983.