VERSION 5.

CFD II, Turbulent and High

Mach Number Flow
MINICOURSE
 Solved with COMSOL Multiphysics 5.0

Turbulent Flow Over a Backward

Facing Step
Introduction
The backward facing step has long been a central benchmark case in computational
fluid dynamics. The geometry is shown in Figure 1.

Figure 1: Backstep geometry. Dimensions in SI units.

Fully developed channel flow enters at the domain from the left. When the flow
reaches the step, it detaches and a recirculation zone is formed behind the step.
Because of the expansion of the channel, the flow slows down and eventually

1 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

reattaches. The flow field is displayed in Figure 2.

Figure 2: Resulting flow field.

Though seemingly simple, the flow field is a challenge for turbulence models that
utilize wall functions. The reason is that wall functions are derived by invoking
equilibrium assumptions. Separation and reattachment do not adhere to these
assumptions and it must therefore be asserted by numerical experiments that the wall
functions can give accurate results even if the underlying theoretical assumptions are
not strictly satisfied. The experiment is motivated by the fact that flow with separation
and subsequent reattachment are of central importance in many engineering
applications.

2 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Model Definition
The model data is taken from Ref. 1. The parameters are given in Table 1. The
Reynolds number based on Vinl and the step height, S, is 4.8·104 and the flow is
therefore clearly turbulent.
TABLE 1: MODEL PARAMETERS

Property Value Description

S 0.0381 m Step height
hc 2·S Inlet channel height
H 3·S Outlet channel height
L1 0.3048 m Inlet channel length
L2 1.3335 m Outlet channel length
Vinl 18.2 m/s Velocity at centre of upstream channel
3
ρ 1.23 kg/m Density
-5
μ 1.79·10 Dynamic viscosity

You build the model in two steps:

1 Simulate flow in a long channel of the same height as the inlet to give inlet boundary
conditions for the actual geometry.
2 Simulate the flow over the backward facing step using the inlet boundary condition
from Step 1.

THE INLET CHANNEL

Ref. 1 suggests to simulate a channel that is 100·hc in length. Because the channel is
symmetric around the midplane, the geometry is taken to be a rectangle with lower left
corner at (x, y) = (0, 0) and upper right corner at (x, y) = (100·hc, 0.5·hc). The upper
boundary at y = 0.5·hc is a symmetry plane and the lower boundary at y = 0 is the wall.

Inlet Boundary Conditions

At the inlet x = 0, a plug flow boundary condition with 3% turbulent intensity and a
turbulent length scale according to Table 3-5 in the Theory for the Turbulent Flow
Interfaces in the CFD Module User’s Guide is prescribed. The inflow velocity cannot
be set directly to 18.2 m/s since the resulting centerline velocity at the outlet then
becomes too high. While it is possible to set up an ODE that automatically computes
the appropriate inlet velocity, it is far easier for small models like this one to find it by
trial and error. A few iterations reveal that an inlet velocity of 16.58 m/s gives a
centerline value at the outlet very close to 18.2 m/s.

3 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Outlet Boundary Conditions

Outlet boundary conditions can give local artifacts at the outlet. One possible strategy
is to elongate the channel and extract data some distance before the outlet. That is
however not necessary since fully developed flow in a channel only has a velocity
component tangential to the wall, that is normal to the outflow. By prescribing that
the outflow must have no tangential component, the outlet artifacts can be removed.

THE BACKWARD FACING STEP

There are two aspects of the backward facing step that need special consideration.

Mesh Generation
It is important to apply a fine enough mesh at the separation point to accurately
capture the creation of the shear layer. It must also be remembered that both the flow
field and turbulence variables can feature strong gradients close to the walls and that
the mesh must be fine enough there to represent these gradients.

Solver Settings
The balance between the turbulence transport equations and the Navier-Stokes
equations is rather delicate. If an iteration brings the flow into a state with unphysically
large gradients, there is a considerable risk that the simulation will diverge. It is
therefore advisable to use the parametric solver to gradually increase the Reynolds
number of the flow. The most robust way is to decrease the viscosity which will be done
in this model.

Results and Discussion

As shown in Figure 3, the recirculation length normalized by the step height becomes
6.73. Ref. 2 gives an experimental result of 7.1. The result provided by COMSOL is
well within the range shown by other investigations (see Ref. 1 and Ref. 3). The
separation lengths in Ref. 1 ranges between 6.12 and 7.24. In Ref. 3, recirculation
lengths between 5.4 and 7.1 are obtained. Furthermore, Ref. 3 shows that the

4 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

recirculation length can differ significantly by just changing some implementation

details in the wall functions.

Figure 3: Contour plot of streamwise velocity equal to zero, colored by x/S where S is the
step height.

Finally, note that the recirculation length can shift quite significantly with the mesh
resolution. The current result does not shift much if the mesh is refined, but coarser
meshes can yield very different recirculation lengths. This emphasizes the need to
ensure that the mesh is fine enough.

References
1. 1st NAFEMS Workbook of CFD Examples. Laminar and Turbulent
Two-Dimensional Internal Flows, NAFEMS, 2000.

2. J. Kim, S.J. Kline, and J.P. Johnston, “Investigation of a Reattaching Turbulent

Shear Layer: Flow Over a Backward Facing Step,” Transactions of the ASME, vol.
102, p. 302, 1980.

3. D. Kuzmin, O. Mierka, and S. Turek, “On the Implementation of the k-ε

Turbulence Model in Incompressible Flow Solvers Based on a Finite Element

5 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Discretization,” Int’l J Computing Science and Mathematics, vol. 1, no. 2–4, pp.
193–206, 2007.

Model Library path: CFD_Module/Single-Phase_Benchmarks/

turbulent_backstep

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select physics tree, select Fluid Flow>Single-Phase Flow>Turbulent
Flow>Turbulent Flow, k-ε (spf).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

DEFINITIONS

Parameters
1 On the Model toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

S 0.0381[m] 0.038100 m Step height
hc 0.0762[m] 0.076200 m Inlet channel height
H 0.1143[m] 0.11430 m Outlet channel height
L1 0.3048[m] 0.30480 m Inlet channel length
L2 1.3335[m] 1.3335 m Outlet channel length

6 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Name Expression Value Description

Vinl 16.58[m/s] 16.580 m/s Centerline inlet
velocity
rhof 1.23[kg/m^3] 1.2300 kg/m³ Density
muf 1.79e-5[Pa*s]* 1.7900E-5 Pa·s Dynamic viscosity
fact
fact 1.0 1.0000 Viscosity scaling
factor

GEOMETRY 1

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size section.
3 In the Width text field, type 100*L1.
4 In the Height text field, type hc/2.
5 On the Geometry toolbar, click Build All.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Density rho rhof kg/m³ Basic
Dynamic viscosity mu muf Pa·s Basic

Turbulent Flow, k-ε (spf)

Inlet 1
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 1 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 In the U0 text field, type Vinl.
5 Locate the Turbulence Conditions section. In the LT text field, type 0.07*hc.

7 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

6 In the IT text field, type 0.03.

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundary 3 only.

Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
2 Select Boundary 4 only.
3 In the Settings window for Outlet, locate the Pressure Conditions section.
4 Select the Normal flow check box.

MESH 1

Mapped 1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Mapped.

Distribution 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Mapped
1 and choose Distribution.
2 Select Boundaries 1 and 4 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 From the Distribution properties list, choose Predefined distribution type.
5 From the Distribution method list, choose Geometric sequence.
6 In the Element ratio text field, type 10.
7 In the Number of elements text field, type 40.

Distribution 2
1 Right-click Mapped 1 and choose Distribution.
2 Select Boundary 2 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 From the Distribution properties list, choose Predefined distribution type.
5 In the Number of elements text field, type 250.
6 In the Element ratio text field, type 2.
7 Select the Reverse direction check box.

8 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Distribution 3
1 Right-click Mapped 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Distribution section.
3 From the Distribution properties list, choose Predefined distribution type.
4 Select Boundary 3 only.
5 In the Number of elements text field, type 250.
6 In the Element ratio text field, type 2.
7 In the Model Builder window, right-click Mesh 1 and choose Build All.

STUDY 1
On the Model toolbar, click Compute.

RESULTS

Velocity (spf)
Check that the flow is fully developed by plotting μT along the centerline.

1D Plot Group 4
1 On the Model toolbar, click Add Plot Group and choose 1D Plot Group.
2 On the 1D plot group toolbar, click Line Graph.
3 Select Boundary 3 only. This is the top surface.
4 In the Settings window for Line Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Model>Component
1>Turbulent Flow, k-ε>spf.muT - Turbulent dynamic viscosity.
5 In the Model Builder window, right-click 1D Plot Group 4 and choose Rename.
6 In the Rename 1D Plot Group dialog box, type Turbulent viscosity in the New
label text field.
7 Click OK.

9 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

8 On the 1D plot group toolbar, click Plot.

As can be seen in the resulting plot, the turbulent viscosity has obtained a constant
value well before the outlet.

With the initial simulation step completed, create the backstep model.

ROOT
On the Model toolbar, click Add Component and choose 2D.

ADD PHYSICS
1 On the Model toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Recently Used>Turbulent Flow, k-ε (spf).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1 ×

5 Click Add to Component in the window toolbar.

10 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

6 On the Model toolbar, click Add Physics to close the Add Physics window.

ADD STUDY
1 On the Model toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies>Stationary.
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Turbulent Flow, k-ε (spf) ×

5 Click Add Study in the window toolbar.

6 On the Model toolbar, click Add Study to close the Add Study window.

GEOMETRY 2

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size section.
3 In the Width text field, type L1+L2.
4 In the Height text field, type hc.
5 Locate the Position section. In the x text field, type -L1.
6 In the y text field, type S.

Point 1 (pt1)
1 On the Geometry toolbar, click Primitives and choose Point.
2 In the Settings window for Point, locate the Point section.
3 In the x text field, type -L1.
4 In the y text field, type S+hc/2.

Union 1 (uni1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Union.
2 Select the objects r1 and pt1 only.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size section.

11 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

3 In the Width text field, type L2.

4 In the Height text field, type S.
5 Click the Build Selected button.

Union 2 (uni2)
1 On the Geometry toolbar, click Booleans and Partitions and choose Union.
2 Select the objects uni1 and r2 only.

Mesh Control Edges 1 (mce1)

1 On the Geometry toolbar, click Virtual Operations and choose Mesh Control Edges.
2 On the object fin, select Boundary 7 only.
3 On the Geometry toolbar, click Build All.
4 Click the Zoom Extents button on the Graphics toolbar.

DEFINITIONS

Linear Extrusion 1 (linext1)

1 On the Definitions toolbar, click Component Couplings and choose Linear Extrusion.
2 In the Settings window for Linear Extrusion, locate the Source Selection section.
3 From the Geometric entity level list, choose Boundary.
4 Select Boundary 4 only.
5 Locate the Source Vertices section. Select the Active toggle button.
6 Select Point 3 only.
7 Select the Active toggle button.
8 Select Point 4 only.
9 Click to expand the Destination section. From the Destination geometry list, choose
Geometry 2.
10 Locate the Destination Vertices section. Select the Active toggle button.
11 Select Point 1 only.
12 Select the Active toggle button.
13 Select Point 2 only.

Linear Extrusion 2 (linext2)

1 On the Definitions toolbar, click Component Couplings and choose Linear Extrusion.
2 In the Settings window for Linear Extrusion, locate the Source Selection section.
3 From the Geometric entity level list, choose Boundary.

12 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

4 Select Boundary 4 only.

5 Locate the Source Vertices section. Select the Active toggle button.
6 Select Point 3 only.
7 Select the Active toggle button.
8 Select Point 4 only.
9 Click to expand the Destination section. From the Destination geometry list, choose
Geometry 2.
10 Locate the Destination Vertices section. Select the Active toggle button.
11 Select Point 3 only.
12 Select the Active toggle button.
13 Select Point 2 only.

MATERIALS

Material 2 (mat2)
1 In the Model Builder window, under Component 2 (comp2) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Density rho rhof kg/m³ Basic
Dynamic viscosity mu muf Pa·s Basic

Turbulent Flow, k-ε 2 (spf2)

Inlet 1
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 1 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 Click the Velocity field button.
5 Specify the u0 vector as

comp1.linext1(comp1.u) x
0 y

13 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

6 Locate the Turbulence Conditions section. Click the Specify turbulence variables
button.
7 In the k0 text field, type comp1.linext1(comp1.k).
8 In the ε0 text field, type comp1.linext1(comp1.ep).

Inlet 2
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 3 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 Click the Velocity field button.
5 Specify the u0 vector as

comp1.linext2(comp1.u) x
0 y

6 Locate the Turbulence Conditions section. Click the Specify turbulence variables
button.
7 In the k0 text field, type comp1.linext2(comp1.k).
8 In the ε0 text field, type comp1.linext2(comp1.ep).

Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
2 Select Boundary 7 only.
3 In the Settings window for Outlet, locate the Pressure Conditions section.
4 Select the Normal flow check box.

MESH 2
1 In the Model Builder window, under Component 2 (comp2) click Mesh 2.
2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Element size list, choose Coarse.

Size 1
1 Right-click Component 2 (comp2)>Mesh 2 and choose Edit Physics-Induced Sequence.
2 In the Model Builder window, under Component 2 (comp2)>Mesh 2 right-click Size 1
and choose Build Selected.

14 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

Size 2
1 Right-click Mesh 2 and choose Size.
2 In the Settings window for Size, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 Select Boundary 9 only.
5 Locate the Element Size section. Click the Custom button.
6 Locate the Element Size Parameters section. Select the Maximum element growth rate
check box.
7 In the associated text field, type 1.03.

Size 3
1 Right-click Mesh 2 and choose Size.
2 In the Settings window for Size, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Point.
4 Locate the Element Size section. Click the Custom button.
5 Locate the Element Size Parameters section. Select the Maximum element size check
box.
6 Select Point 5 only.
7 In the associated text field, type 5e-4.

Boundary Layer Properties 1

1 In the Model Builder window, expand the Component 2 (comp2)>Mesh 2>Boundary
Layers 1 node, then click Boundary Layer Properties 1.
2 In the Settings window for Boundary Layer Properties, locate the Boundary Layer
Properties section.
3 In the Thickness adjustment factor text field, type 2.
4 In the Number of boundary layers text field, type 6.
5 Click the Build All button.

STUDY 2

Step 1: Stationary
1 In the Model Builder window, under Study 2 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.

15 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

fact 5 1

6 Click to expand the Values of dependent variables section. Locate the Values of
Dependent Variables section. Select the Values of variables not solved for check box.
7 From the Method list, choose Solution.
8 From the Study list, choose Study 1, Stationary.
9 On the Model toolbar, click Compute.

RESULTS

Velocity (spf2)
Check that the wall lift-off is 11.06 almost everywhere by selecting the Wall Resolution
(spf2) plot group.

Next, reproduce the flow-field plot with the following steps:

Velocity (spf2)
1 In the Model Builder window, under Results right-click Velocity (spf2) and choose
Streamline.
2 In the Settings window for Streamline, locate the Streamline Positioning section.
3 From the Positioning list, choose Uniform density.
4 In the Separating distance text field, type 0.007.
5 Locate the Coloring and Style section. From the Color list, choose White.
6 On the 2D plot group toolbar, click Plot.

Finally, visualize the recirculation length.

2D Plot Group 8
1 On the Model toolbar, click Add Plot Group and choose 2D Plot Group.
2 In the Settings window for 2D Plot Group, locate the Data section.
3 From the Data set list, choose Study 2/Solution 2 (3).
4 Right-click Results>2D Plot Group 8 and choose Contour.

16 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

5 In the Settings window for Contour, click Replace Expression in the upper-right
corner of the Expression section. From the menu, choose Model>Component
2>Turbulent Flow, k-ε 2>Velocity field>u2 - Velocity field, x component.
6 Locate the Levels section. From the Entry method list, choose Levels.
7 Right-click Results>2D Plot Group 8>Contour 1 and choose Color Expression.
8 In the Settings window for Color Expression, locate the Expression section.
9 In the Expression text field, type x/S.
10 In the Model Builder window, right-click 2D Plot Group 8 and choose Rename.
11 In the Rename 2D Plot Group dialog box, type Recirculation length in the New
label text field.
12 Click OK.
13 On the 2D plot group toolbar, click Plot.

17 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P
 Solved with COMSOL Multiphysics 5.0

18 | TU R B U L E N T F L O W O V E R A B A C K W A R D F A C I N G S T E P