Numerical solutions of 2D laminar flow over a backwardfacing step
at high Reynolds numbers are presented. The governing 2D steady
incompressible NavierStokes equations are solved with a very
efficient finite difference numerical method which proved to be highly stable even at
very high Reynolds numbers. Present solutions of the laminar flow
over a backwardfacing step are compared with experimental and
numerical results found in the literature.
Fluid flows in channels with flow separation and reattachment of the
boundary layers are encountered in many flow problems. Typical
examples are the flows in heat exchangers and ducts. Among this type
of flow problems, a backwardfacing step can be regarded as having
the simplest geometry while retaining rich flow physics manifested
by flow separation, flow reattachment and multiple recirculating
bubbles in the channel depending on the Reynolds number and the
geometrical parameters such as the step height and the channel
height.
In the literature, it is possible to find many numerical studies on
the 2D steady incompressible flow over a backwardfacing step. In
these studies one can notice that there used to be a controversy on
whether it is possible to obtain a steady solution for the flow over
a backwardfacing step at Re=800 or not. For example, this fact
was discussed in Gresho et al. [17] in detail and
they concluded that the flow is stable and computable at Re=800.
Guj and Stella [19], Keskar and Lyn [22], Gartling
[16], Papanastasiou et al. [28],
Rogers and Kwak [31], Kim and Moin [23], Comini
et al. [6], Barton [3,4], Sani and
Gresho [32], Sheu and Tsai [33], Biagioli
[5], Grigoriev and Dargush [18], Morrison
and Napolitano [27] and also Zang et al.
[40] are just example studies that presented solutions up to
Re=800 obtained with different numerical methods, we pick as
example from the literature. In these studies one can also notice
that, Re=800 was believed to be the limit for obtaining a steady
solution.
Apart from these studies
[17,19,22,16,28,31,23,6,3,4,32,33,5,18,27,40]
and many others found in the literature, Ramsak [30]
presented numerical solutions for the steady flow over a
backwardfacing step for Re=1000. Similarly, Cruchaga [7]
solved the steady backwardfacing step flow using finite element
method and obtained steady numerical solutions up to Re=5500,
although he mainly presented solutions for Re=800 case. These two
studies [30,7] have clearly shown that it is
possible to obtain numerical solutions well beyond Re=800.
In the literature it is also possible to find studies that
investigate the numerical stability of the flow over a
backwardfacing step. For example, Fortin et al. [15] have
studied the stability of the 2D steady incompressible flow over a
backwardfacing step until a Reynolds number of 1600. They stated
that with the grid mesh and the boundary conditions they have used,
no pair of eigenvalues has crossed the imaginary axes such that the
flow over a backwardfacing step is stable up to the maximum
Reynolds number (Re=1600) they considered in their study. As
another example, Barkley et al. [2] have done a
computational stability analysis to the flow over a backwardfacing
step with an expansion ratio of 2. They continued their stability
analysis up to Reynolds number of 1500 and stated that the flow
remains linearly stable to twodimensional perturbations and
moreover shows no evidence of any nearby twodimensional
bifurcation. These two studies, [2,15], are important
in the sense that their twodimensional computational stability
analysis show that the flow over a backwardfacing step is temporally
stable at high Reynolds numbers, that is to say there exist steady
solutions of the flow over a backwardfacing step at high Reynolds
numbers.
Erturk et al. [9] introduced an efficient, fast and
stable numerical formulation for the steady incompressible
NavierStokes equations. Their method solve the streamfunction and
vorticity equations separately, and the numerical solution of each
equation requires the solution of two tridiagonal systems. Solving
tridiagonal systems are computationally efficient and therefore they
were able to use very fine grid mesh in their solution. Using this
numerical formulation, they solved the very well known benchmark
problem, the steady flow in a square driven cavity, up to Reynolds
number of 21000 using a 601×601 fine grid mesh. The numerical
formulation introduced by Erturk et al. [9] proved
to be stable and effective at very high Reynolds number flows
([9,10,12]) and also at flows with
nonorthogonal grid mesh at extreme skew angles ([13]).
In this study, we will present numerical solutions of 2D steady
incompressible flow over a backwardfacing step at high Reynolds
numbers. The 2D steady incompressible NavierStokes equations in
streamfunction and vorticity formulation will be solved using the
efficient numerical method introduced by Erturk et al.
[9] and detailed solutions will be presented.
Armaly et al. [1] have done experiments on the
backwardfacing step flow and have presented valuable experimental
benchmark results to the literature. In their experiments the
expansion ratio, ie. the ratio of the channel height downstream of
the step to the channel height upstream of the step, was equal to
1.942. Lee and Mateescu [24] have also done experiments on the
flow over a backwardfacing step with expansion ratios of 1.17 and
2.0. We note that, these experimental results will form a basis to
compare our numerical solutions for moderate Reynolds numbers. For
this, in this study we will consider backwardfacing step flow with
both 1.942 and 2.0 expansion ratios.
While the flow over a backwardfacing step serve as an interesting
benchmark flow problem for many numerical studies, some studies
stated that the inlet and exit boundary condition used for the model
problem can affect the numerical solution. In the literature, Barton
[3] studied the entrance effect for flow over a backwardfacing facing step. He used the QUICK scheme for numerical solution. He
stated that when using an inlet channel upstream of the step,
significant differences occur for low Reynolds numbers, however,
they are localized in the sudden expansion region. He also stated
that if the Reynolds number is high then shorter reattachment and
separation lengths are predicted. Cruchaga [7] solved
the backwardfacing step flow both with using an inlet channel and
without using an inlet channel and he obtained slightly different
solutions for Re=800 for the two case. Papanastasiou et al.
[28] studied the effect of the outflow boundary
conditions on the numerical solutions in a backwardfacing step
flow. They presented a free boundary condition as an outflow
boundary condition where this condition is inherent to the weak
finite element formulation of the momentum equations such that it
minimizes the energy functional for the unbounded Stokes flow. In
his detailed study on the flow over a backwardfacing step, Gartling
[16] stated the importance of the outflow boundary
condition for the considered flow, also Wang and Sheu [37]
discussed on the use of free boundary conditions specifically
together with a discussion on various outflow boundary conditions.
In this study we will examine the effect of both the inlet channel
and the outflow boundary condition on the numerical solution. We
believe that it is not only the inlet channel and outflow boundary
condition that can affect the numerical solution inside the
computational domain but the location of the outflow boundary is
also very important for the accuracy of the numerical solution. In
this study we will examine the effect of the location of the outflow
boundary on the numerical solution by considering different
computational domains with different locations of the outflow
boundary.
In a very important study Yee et al. [39] studied the
spurious behavior of the numerical schemes. They showed that for
backwardfacing step flow when a coarse grid mesh is used, one can
obtain a spurious oscillating numerical solution. They have clearly
stated that fine grids are necessary in order to avoid spurious
solutions for the backwardfacing step flow. Similar to the comments
of Yee et al. [39], Erturk et al. [9]
have reported that for square driven cavity flow when a coarse grid
mesh was used, they observed that the numerical solution was not
converging to the steady state solution, but it was oscillating at
high Reynolds numbers. They reported that when a finer grid mesh was
used, the oscillating behavior of the numerical solution disappeared
and it was possible to obtain a steady solution. They stated that
when finer grids are used, the Mesh Reynolds number defined as
Re_{m}=[(u Dh)/(n)] decreases and this improves the
numerical stability characteristics of the numerical scheme used
(see [36] and [38]), and allows high Reynolds
numbered solutions computable.
In this study, following Yee et al. [39] and Erturk
et al. [9], we will use a fine grid mesh in order to be
able to obtain steady state numerical solutions at high Reynolds
numbers.
This paper presents highly accurate numerical solutions of the
backwardfacing step flow obtained using; the efficient numerical
method introduced by Erturk et al. [9], an exit
boundary located very far away downstream of the step,
nonreflecting boundary conditions at the exit boundary, an inlet
boundary located very upstream of the step, a fine grid mesh and a
convergence criteria that is close to machine accuracy. In the
following section, the details of the numerical method will be
given. In the next section, the details on the numerical procedure
such as the grid mesh, the inlet and exit boundary conditions and
also the wall boundary conditions used in our computations will be
presented. Then, finally, comparisons of our results with
experimental and numerical solutions found in the literature will be
done and detailed numerical solutions will be presented together
with a brief discussion.
We consider the 2D steady incompressible NavierStokes equations in
streamfunction (y) and vorticity (w) formulation. In
nondimensional form, they are given as
¶^{2}y
¶x^{2}
+
¶^{2}y
¶y^{2}
= w
(1)
1
Re
æ è
¶^{2} w
¶x^{2}
+
¶^{2} w
¶y^{2}
ö ø
=
¶y
¶y
¶w
¶x

¶y
¶x
¶w
¶y
(2)
where, Re is the Reynolds number, and x and y are
the Cartesian coordinates. The Reynolds number is defined as
Re=[UD/(n)] where U is the inlet mean velocity or in other
words two thirds of the maximum inlet velocity and D is the
hydraulic diameter of the inlet channel where it is equivalent to
twice the inlet channel height, ie. D=2h_{i} as shown in Figure 1a.
For the numerical solution of streamfunction and vorticity
equations, we preferred to use the efficient numerical method
proposed by Erturk et al. [9]. In this study we
will only describe the method shortly and the reader is referred to
[9,12] for details of the numerical method.
The NavierStokes equations (1 and 2) are nonlinear equations
therefore they need to be solved in an iterative manner. In order to
have an iterative numerical algorithm we assign pseudo time
derivatives to equations (1 and 2). Using implicit Euler
approximation for the pseudo time, the obtained finite difference
equations in operator notation are the following
æ è
1  Dt d_{xx}  Dt d_{yy}
ö ø
y^{n+1} = y^{n} + Dt w^{n}
(3)
æ è
1  Dt
1
Re
d_{xx}  Dt
1
Re
d_{yy} + Dt y_{y}^{n} d_{x} Dt y_{x}^{n} d_{y}
ö ø
w^{n+1} = w^{n}
(4)
where subscripts denote derivatives and also d_{x}
and d_{y} denote the first derivative and similarly
d_{xx} and d_{yy} denote the second derivative finite
difference operators in x and ydirections respectively, for
example
d_{x} q
=
q_{i+1,j}q_{i1,j}
2Dx
+ O (Dx^{2})
d_{xx} q
=
q_{i+1,j}2q_{i,j}+q_{i1,j}
Dx^{2}
+ O (Dx^{2})
(5)
where i and j are the grid index and q denote
any differentiable quantity. We note that, since we have used three
point central differencing the presented solutions in this study are
second order accurate.
The finite difference equations (3 and 4) are in fully implicit form
where each equation requires the solution of a large banded matrix
which is not computationally efficient. Instead, the fully implicit
equations (3 and 4) are spatially factorized such that
æ è
1  Dt d_{xx}
ö ø
æ è
1  Dt d_{yy}
ö ø
y^{n+1} = y^{n} + Dt w^{n}
(6)
æ è
1  Dt
1
Re
d_{xx} + Dt y_{y}^{n} d_{x}
ö ø
æ è
1  Dt
1
Re
d_{yy} Dt y_{x}^{n} d_{y}
ö ø
w^{n+1} = w^{n}
(7)
We note that, now each equation (6 and 7)
requires the solution of two tridiagonal systems which is
computationally more efficient. However, spatial factorization
introduces Dt^{2} terms into the left hand side of the
equations and these Dt^{2} terms remain in the equations
even at the steady state. In order to have the correct physical form
at the steady state, Erturk et al. [9] add the same
amount of Dt^{2} terms to the right hand side of the
equations such that the final form of the finite difference equations become the following
æ è
1  Dt d_{xx}
ö ø
æ è
1  Dt d_{yy}
ö ø
y^{n+1} = y^{n} + Dt w^{n} +
æ è
Dt d_{xx}
ö ø
æ è
Dt d_{yy}
ö ø
y^{n}
(8)
æ è
1  Dt
1
Re
d_{xx} + Dt y_{y}^{n} d_{x}
ö ø
æ è
1  Dt
1
Re
d_{yy} Dt y_{x}^{n} d_{y}
ö ø
w^{n+1} = w^{n}
+
æ è
Dt
1
Re
d_{xx}  Dt y_{y}^{n} d_{x}
ö ø
æ è
Dt
1
Re
d_{yy} + Dt y_{x}^{n} d_{y}
ö ø
w^{n}
(9)
We note that, at the steady state the Dt^{2} terms
on the right hand side of the equations cancel out the Dt^{2} terms due to the factorization on the left hand side, hence
at the steady state the correct physical representation is
recovered. The solution methodology for the two equations (8 and 9) involves a
twostage timelevel updating. For the streamfunction equation (8), the
variable f is introduced such that
æ è
1  Dt d_{yy}
ö ø
y^{n+1} = f
(10)
and
æ è
1  Dt d_{xx}
ö ø
f = y^{n} + Dt w^{n} +
æ è
Dt d_{xx}
ö ø
æ è
Dt d_{yy}
ö ø
y^{n}
(11)
In equation (11) f is the only unknown. First, this
equation is solved for f at each grid point. Following this, the
streamfunction (y) is advanced into the new time level using
equation (10).
In a similar fashion for the vorticity equation (9), the variable g is
introduced such that
æ è
1  Dt
1
Re
d_{yy}  Dt y_{x}^{n} d_{y}
ö ø
w^{n+1} = g
(12)
and
æ è
1  Dt
1
Re
d_{xx} + Dt y_{y}^{n} d_{x}
ö ø
g = w^{n}
+
æ è
Dt
1
Re
d_{xx}  Dt y_{y}^{n} d_{x}
ö ø
æ è
Dt
1
Re
d_{yy} + Dt y_{x}^{n} d_{y}
ö ø
w^{n}
(13)
As we did the same with f, g is determined at every
grid point using equation (13), then vorticity (w) is advanced into the next time level using equation (12).
The schematics of the backwardfacing step flow considered in this
study are shown in Figure 1. In this study, the inlet boundary is
located 20 step heights upstream of the step, as shown as L_{1} in
Figure 1a. In this inlet channel, L_{1}, we used 500 uniform grids
as also shown in Figure 1b. The exit boundary is chosen as 300 step
heights away from the step, as shown as L_{2} in Figure 1a. In
order to have high accuracy in the vicinity of the step, in
xdirection from the step to a distance of 100 step heights we
used 2500 uniform fine grids as shown as L_{3} in Figure 1b. From
100 step heights distance to the exit boundary, also as shown as
L_{4} in Figure 1b, we used 1250 stretched grid points in order to
be able to have the location of the exit boundary far from the step,
as this approach was also used in Gartling [16]. In
L_{4}, the grids are stretched smoothly starting from the constant
grid spacing used in L_{3} and this was done using Robert's
stretching transformation of the original uniform grid (see
Tannehill, Anderson and Pletcher [36]). The
transformation used is given by
x=L_{3} + L_{4}
(b+1)(b1)[(b+1)/(b1)]^{1[`x]}
[(b+1)/(b1)]^{1 [`x]}+1
(14)
where [`x]=[0,1] represent the original uniformly
spaced grid points and x are the stretched grid points and b is the stretching parameter.
In ydirection, we have used 101 uniform grids as it is shown in
Figure 1b.
Barton [3] studied the entrance effect for flow over a
backwardfacing step and reported that the inlet channel upstream of
the step affects the numerical solution significantly at low
Reynolds numbers. In this study the inlet boundary is located 20
step heights upstream of the step. Here we would like to note that
among the numerical studies of flow over a backwardfacing step
found in the literature, this study utilizes the largest inlet
channel, L_{1}. Cruchaga [7] solved the backwardfacing step both with and without using an inlet channel and
obtained slightly different solutions. In this study we will try to
examine the effect of the inlet channel on the considered flow.
At the inlet boundary we imposed that the flow is fully developed
Plane Poiseuille flow between parallel plates such that the inlet
velocity profile is parabolic.
At the exit boundary we used a nonreflecting boundary condition
such that any wave generated in the computational domain could pass
through the exit boundary and leave with out any reflection back
into the computational domain. For details on the subject the reader
is referred to the study of Engquist and Majda [8] in
which the nonreflecting boundary condition concept is first
introduced in the name of Äbsorbing Boundary Condition", and also
to Jin and Braza [21] for a review of nonreflecting boundary
conditions. Liu and Lin [26] attached a buffer region to the
physical domain to damp erroneous numerical fluctuations. In this
region they [26] added a buffer function to the streamwise
second order derivatives in the momentum equations such that the
reflected outgoing waves from an artificially truncated outlet are
thus absorbed. This approach for the exit boundary condition has
been used in various similar studies ([20,11,12])
successfully.
As an exit boundary condition our approach is; at near the exit
boundary we kill the elliptic (¶^{2}/¶x^{2}) terms
in the governing equations (1 and 2) gradually in a buffer zone. To
accomplish this, these elliptic terms are multiplied by a weighting
factor s. At the beginning of the buffer zone, we set s=1 and at
the end of the buffer zone, it is zero, s=0. In between, the
weighting factor changes as
s(i)=
tanh(4) + tanh(arg)
2 tanh(4)
(15)
where
arg=4
æ è
1
2(iibuf)
(imaxibuf)
ö ø
(16)
where i is the numerical streamwise index, imax is the
numerical index of the last grid point in streamwise direction and
ibuf is the index i of the first grid point at the beginning of
the buffer zone. The buffer zone contains 100 grid points as shown
in Figure 1b.
We note that, sufficiently away from the step, the solution of
the 2D steady incompressible NavierStokes equations should
approach to a fully developed Plane Poiseuille flow between parallel
plates, such that the velocity profile is parabolic and also the
second xderivatives of streamfunction and vorticity variables
should be equal to zero, ie. [(¶^{2} y)/(¶x^{2})]=0 and [(¶^{2} w)/(¶x^{2})]=0. As a
matter of fact, the latter two condition is equivalent of what we
impose at the exit boundary. Therefore, when the exit boundary is
sufficiently away from the step, the nonreflecting boundary
condition we used at the exit boundary physically gives the exact
solution at the exit boundary. We note that, in terms of the
primitive variables, the conditions we apply at the exit boundary is
equivalent to [(¶^{2} y)/(¶x^{2})]=[(¶v)/(¶x)]=0 and also [(¶^{2}w)/(¶x^{2})]=[(¶^{3} u)/(¶x^{2}¶y)]  [(¶^{3} v)/(¶x^{3})] = 0 where u and v are
the velocity components in x and ydirections respectively.
After the step, physically the flow needs some streamwise distance
to adjust and become fully developed. The nonreflecting boundary
condition we apply is equivalent to assuming the flow to be fully
developed. Our extensive numerical tests have shown that, using the
nonreflecting boundary condition we explained above, if we placed
the exit boundary close to the step, such that it was not
sufficiently away from the step, the obtained velocity profile at
the exit plane was not the analytical parabolic profile. The result
was inconsistent with the boundary condition, such that, we applied
that the flow to be parallel and fully developed however the
computed output was contradicting with the boundary condition we
apply and also contradicting with the analytical solution as well.
This inconsistency was due to the fact that the exit boundary was
not sufficiently away from the step and also the flow becomes
parallel and fully developed only when sufficiently away from step.
Leone [25] compared his numerical solutions at the outflow
boundary with the Poiseuille flow solution, in order to test the
accuracy of the numerical solution and also the boundary condition
he used. Following Leone [25], in this study we decided to
use the fact that the velocity profile at the exit boundary should
be parabolic, as a mathematical check on our numerical solution. To
do this, we systematically move the location of the exit boundary
away from the step, ie. that is to say we used a larger
computational domain in xdirection, and solve the flow problem
until the obtained exit velocity profile agrees with the analytical
solution. In our computations, the exit boundary is located
sufficiently away (300 step heights) from the step, such that, even
for the highest Reynolds number considered in this study
(Re=3000), the obtained numerical velocity profile at the exit
boundary agrees with the analytical parabolic profile with a maximum
difference less than 0.05 %. The nonreflecting boundary condition
we have used, the large computational domain with exit boundary
sufficiently away from the step and the mathematical check of the
numerical solution with the analytical solution ensure that the
presented numerical results in this study are indeed very accurate.
The vorticity value at the wall is calculated using Jensen's formula
(see Tannehill, Anderson and Pletcher [36])
w_{0}=
3.5 y_{0}4 y_{1} + 0.5 y_{2}
Dy^{2}
(17)
where subscript 0 refers to the points on the wall and 1 refers to
the points adjacent to the wall and 2 refers to the second line of
points adjacent to the wall and also Dy is the grid spacing.
We note that for the vorticity on the wall, Jensen's formula
provides second order accurate results. We also note that near the
inlet and exit boundaries the flow is fully developed Plane
Poiseuille flow between parallel plates and the velocity profile is
parabolic, and hence the streamfunction profile is third order and
vorticity profile is linear in ydirection. Therefore, Jensen's
formula not only provide us numerical solutions at the boundary with
accuracy matched to the numerical method used inside the
computational domain, it also provides the exact analytical value
near the inlet and exit boundaries.
During our computations as a measure of convergence to the steady
state, we monitored three residual parameters. The first residual
parameter, RES1, is defined as the maximum absolute residual of the
finite difference equations of steady streamfunction and vorticity
equations (1 and 2). These are respectively given as
RES1_{y}
= max
æ è
ê ê
y_{i1,j}^{n+1}2y_{i,j}^{n+1}+y_{i+1,j}^{n+1}
Dx^{2}
+
y_{i,j1}^{n+1}2y_{i,j}^{n+1}+y_{i,j+1}^{n+1}
Dy^{2}
+w_{i,j}^{n+1}
ê ê
ö ø
RES1_{w}
= max
æ è
ê ê
1
Re
w_{i1,j}^{n+1}2 w_{i,j}^{n+1}+w_{i+1,j}^{n+1}
Dx^{2}
+
1
Re
w_{i,j1}^{n+1}2w_{i,j}^{n+1} +w_{i,j+1}^{n+1}
Dy^{2}

y_{i,j+1}^{n+1}y_{i,j1}^{n+1}
2 Dy
w_{i+1,j}^{n+1}w_{i1,j}^{n+1}
2 Dx
+
y_{i+1,j}^{n+1}y_{i1,j}^{n+1}
2 Dx
w_{i,j+1}^{n+1}w_{i,j1}^{n+1}
2 Dy
ê ê
ö ø
(18)
The magnitude of RES1 is an indication of the degree
to which the solution has converged to steady state. In the limit
RES1 would be zero.
The second residual parameter, RES2, is defined as the maximum
absolute difference between two iteration steps in the
streamfunction and vorticity variables. These are respectively given
as
RES2_{y}
=
max
æ è
ê ê
y_{i,j}^{n+1}  y_{i,j}^{n}
ê ê
ö ø
RES2_{w}
=
max
æ è
ê ê
w_{i,j}^{n+1}  w_{i,j}^{n}
ê ê
ö ø
(19)
RES2 gives an indication of the significant digit on
which the code is iterating.
The third residual parameter, RES3, is similar to RES2, except that
it is normalized by the representative value at the previous time
step. This then provides an indication of the maximum percent change
in y and w in each iteration step. RES3 is defined as
RES3_{y}
=
max
æ è
ê ê
y_{i,j}^{n+1}  y_{i,j}^{n}
y_{i,j}^{n}
ê ê
ö ø
RES3_{w}
=
max
æ è
ê ê
w_{i,j}^{n+1}  w_{i,j}^{n}
w_{i,j}^{n}
ê ê
ö ø
(20)
RES3_{y}
=
max
æ è
ê ê
y_{i,j}^{n+1}  y_{i,j}^{n}
y_{i,j}^{n}
ê ê
ö ø
RES3_{w}
=
max
æ è
ê ê
w_{i,j}^{n+1}  w_{i,j}^{n}
w_{i,j}^{n}
ê ê
ö ø
(20)
In our calculations, for all Reynolds numbers we
considered that convergence was achieved when both
RES1_{y} £ 10^{10} and RES1_{w} £ 10^{10} was achieved. Such a low value was chosen to
ensure the accuracy of the solution. At these convergence levels the
second residual parameters were in the order of RES2_{y} £ 10^{17} and RES2_{w} £ 10^{15},
that means the streamfunction and vorticity variables are accurate
to 16^{th} and 14^{th} digit accuracy respectively at a grid
point and even more accurate at the rest of the grids. Also at these
convergence levels the third residual parameters were in the order
of RES3_{y} £ 10^{14} and RES3_{w} £ 10^{13}, that means the streamfunction and vorticity
variables are changing with 10^{12}% and 10^{11}% of their
values respectively in an iteration step at a grid point and even
with less percentage at the rest of the grids. These very low
residuals ensure that our solutions are indeed very accurate.
At the beginning of this study, we considered a computational domain
with uniform fine grids that has the exit boundary located 60 step
heights away from the step as it was also used the same in
[16,25,22,6,7], and the inlet boundary
was located right at the step, ie. there was no inlet channel. Using
the described numerical method and the boundary conditions, we
obtained steady numerical solutions for up to Reynolds number of
1000 and above this Reynolds number our numerical solution was not
converging, but it was oscillating. However, as explained above,
when we look at the computed velocity profile at the exit boundary,
in the Reynolds number range of 700 £ Re £ 1000,
we see that our computed profiles did not match with the analytical
parabolic profile, suggesting that at these Reynolds numbers the
flow actually was not fully developed at the exit plane. Then we
decided to use a larger computational domain with uniform fine grids
with the exit boundary located at 100 step heights away from the
step. This time, at the considered Reynolds numbers, the computed
velocity profile at the exit boundary was very close to the
parabolic profile. Surprisingly, when we used a larger computational
domain, this time we were able to obtain steady solutions for up to
Reynolds number of 1600. This fact suggested that as the outflow
boundary was moved away from the step location, it was possible to
obtain steady numerical solutions of the flow over a backwardfacing
step for higher Reynolds numbers. Encouraged by this, we decided to
increase computational domain again and solve for larger Reynolds
numbers. Continuing this process, as the outflow boundary was moved
away from the step with using uniform fine grid mesh, the
computations became time consuming and inefficient after some point,
since we needed to use more grid points every time. Then we decided
to use a gradually graded grid mesh as it was also used by Gartling
[16]. As described in Section 3, from step to 100 step
heights distance in xdirection we decided to use uniform fine
grid mesh for accuracy and from 100 step heights distance to 300
step heights distance where the exit boundary is located, we used
1250 gradually graded grid mesh. With this grid mesh, we were able
to obtain steady numerical solutions of the 2D backwardfacing step
flow up to Re=3000. We note that we have also tried several
different grid meshes with different grid points with different grid
stretching and exit boundary locations, to make sure that our
numerical solutions are independent of the location of the exit
boundary and also grid mesh used. We also note that, among the
studies on backwardfacing step flow found in the literature, the
present study has utilized the largest computational domain with the
largest number of grid points.
At the exit boundary we have used a nonreflecting boundary
condition. As stated before, the condition we applied at the exit
boundary is equivalent of assuming the flow to be parallel and hence
fully developed. If the flow is fully developed the velocity profile
should be the Plane Poiseuille velocity profile, ie. the parabolic
profile. In a case, if the flow is assumed to be parallel at the
exit boundary however the obtained velocity profile is not the
parabolic profile, then this would indicate that there is an
inconsistency with the input boundary condition and the output
numerical solution. This would only occur when a small computational
domain with the location of the exit boundary that is not
sufficiently away from the step is considered. We decided to use the
Plane Poiseuille velocity profile as a mathematical check on our
velocity profile at the exit boundary, as this was also done by
Leone [25]. Figure 2 shows the computed vorticity,
uvelocity and streamfunction profiles at the outflow boundary for
the lowest and the highest Reynolds number considered in this study
(Re=100 and Re=3000) together with the analytical exact solution
of Plane Poiseuille flow.
This figure shows that, from the lowest to
the largest for the whole range of Reynolds number considered, our
computed numerical solutions are indeed in excellent agreement with
the analytical solution.
As mentioned before, Barton [3] studied the entrance
effect for flow over a backwardfacing step and found that the inlet
channel upstream of the step affects the numerical solution at low
Reynolds numbers. During our numerical experimentation, since we
have noticed that the location of the exit boundary affected the
numerical solution, we decided to use a very long inlet channel to
minimize its possible effect on the numerical solution, because
physically the longer the inlet channel the less likely it could
affect the solution. For this purpose, upstream of the step
location, we decided to use an inlet channel with 20 step heights
length. Figure 3 shows how the velocity profiles changes in
streamwise direction from the inlet boundary to the step location at
various Reynolds numbers between 100 £ Re £ 3000.
In this figure we have plotted the uvelocity profiles at five
streamwise locations such as at x/h=20 being the inlet boundary,
x/h=15, x/h=10, x/h=5 and x/h=0 being the step location.
As seen in Figure 3, at Re=100 the velocity profiles at x/h=15, 10 and 5 match with the parabolic profile at the inlet (x/h=20) perfectly, however the velocity profile at the step location
(x/h=0) deviates from the parabolic profile. From this figure we
see that this deviation is large at small Reynolds number and as the
Reynolds number increases the deviation of the velocity profile at
the step from the inlet parabolic profile gets smaller in accordance
with Barton [3]. This result is not surprising such that,
given the elliptic nature of the NavierStokes equations, the effect
of the step propagates back into the inlet channel and affects the
flow upstream. This effect is large at small Reynolds numbers and as
the Reynolds number increases the flow becomes convectively dominant
and this effect gets smaller. Zang et al. [40] stated
that, in a study Perng [29] has shown that with an entrance
section, the velocity profile right at the expansion deviates from a
parabola and appears to have a "downwash". This in fact is what
we observe at small Reynolds numbers as shown in Figure 3. From
Figure 3, we conclude that in backwardfacing step flow studies
there needs to be an at least 5 step heights long inlet channel
upstream of the step, since in the figure, upstream of x/h=5 there is virtually no difference in velocity profiles between the
parabolic input velocity profile at all Reynolds numbers.
Before presenting our high Reynolds number (Re £ 3000) numerical solutions of the flow over a backwardfacing step, we
would like to present solutions in which we compare with the
experimental and numerical solutions found in the literature in
order to demonstrate the accuracy of our numerical solutions.
First, in order to be able to compare our numerical solutions with
the experimental solutions of Armaly et al. [1], we
considered a backwardfacing step with 1.942 expansion ratio. Using
the numerical procedure described above, we have solved the steady
2D Navier Stokes equations up to Re £ 1500. Figure 4, 5, 6, 7 and 8 show the uvelocity profiles at several xlocations
at Reynolds number of Re=100, 389, 1000, 1095 and 1290
respectively.
In each of the Figures 4, 5, 6, 7 and 8, the top
velocity profiles are scanned from the experimental work of Armaly
et al. [1] and then digitally cleaned, also the
bottom profiles are our computed velocity profiles at the
corresponding xlocations drawn to the same scale for comparison.
As we can see from Figure 4, 5, 6 and 7, at Reynolds numbers of
Re=100, 389, 1000 and 1095, our computed velocity profiles agree
well with that of experimental results of Armaly et al.
[1]. In Figure 8, at Re=1290, our computed velocity
profiles agree well with that of experimental results of Armaly
et al. [1] at xlocations close to the step, and the
agreement is moderate at far downstream xlocations. We note that,
Armaly et al. [1] have stated that, in their
experiments the flow starts to show signs of transition near
Re=1200. Therefore, the apparent difference at far downstream
locations in Figure 8 at Re=1290 between our computed velocity
profiles and experimental velocity profiles of Armaly et al.
[1] is due to the 3D effects observed in experiments.
Figures 4, 5, 6, 7 and 8 show that our 2D steady numerical
solutions agree well with the laminar experimental results of Armaly
et al. [1].
Figure 1c schematically illustrates the recirculating regions occur
in the backwardfacing step flow and Table 1 tabulates our numerical
solutions for the backwardfacing step flow with 1.942 expansion
ratio up to Reynolds number of 1500. For comparison, in Figure 9 we
plotted our computed main recirculation region length normalized by
the step height versus the Reynolds number together with the
experimental results of Armaly et al. [1]. In this
figure we see that at low Reynolds numbers our computed results
agree well with experimental results of Armaly et al.
[1], and at higher Reynolds numbers the agreement is
moderate.
We note that, since we only considered the backwardfacing step with
a 1.942 expansion ratio, for the purpose of comparison with the
experimental results of Armaly et al. [1], we solved
this case only up to Re=1500. Our main focus will be on a
backwardfacing step with a 2.0 expansion ratio as it is done in
most of the numerical studies. For the backwardfacing step with
expansion ratio of 2.0, we have solved the steady incompressible
NavierStokes equation for variety of Reynolds numbers.
Gartling [16] tabulated detailed results at x/h=14 and
x/h=30 streamwise locations. Table 2, 3 and 4 tabulates our
computed profiles of flow variables across the channel at x/h=6,
x/h=14 and x/h=30 respectively, for Re=800 case. Figure 10
compares our numerical results tabulated in Table 3 and 4, with that
of Gartling [16].
In this figure we plotted the
horizontal velocity (u), vertical velocity (v), vorticity
(w), horizontal velocity streamwise gradient ([(¶u)/(¶x)]), horizontal velocity vertical gradient
([(¶u)/(¶y)]), vertical velocity streamwise
gradient ([(¶v)/(¶x)]) profiles at x/h=14 and
30 streamwise locations. In Figure 10 we can see that our computed
u, w and [(¶u)/(¶y)] profiles agree
good with that of Gartling [16] at both x/h=14 and 30
streamwise locations. In this figure we also see that while our
computed v, [(¶u)/(¶x)] and [(¶v)/(¶x)] profiles agree with that of Gartling [16]
at x/h=30, they do not agree at x/h=14 streamwise location. This
was interesting such that for these profiles (v, [(¶u)/(¶x)], [(¶v)/(¶x)]) the results of
this study and Gartling [16] agree well with each other
away from the step however they differ from each other at
locations closer to the step. The explanation to this difference was
evident in Fig. 3 of Cruchaga [7]. Cruchaga
[7] solved the backwardfacing step flow both with using
an inlet channel upstream of the step and also without using an
inlet channel. As we did the same in Figure 10, Cruchaga
[7] compared the same profiles with that of Gartling
[16] at both streamwise locations. While the results of
Cruchaga [7] obtained without using an inlet
channel agreed well with that of Gartling [16], the
results of Cruchaga [7] obtained with using an
inlet channel differed from that of Gartling [16]. We
note that Gartling [16] did not use an inlet channel in
his computations. We also note that our results and the results of
Cruchaga [7] obtained with using an inlet channel agree
well with each other. The results of Cruchaga [7] and
Figure 10 clearly show that for the backwardfacing step flow an
inlet channel is necessary for accurate numerical solutions. It is
evident that, some of the flow quantities such as u, w and
[(¶u)/(¶y)] are not very sensitive to the
absence of an inlet channel, however the flow quantities v,
[(¶u)/(¶x)] and [(¶v)/(¶x)]
are very sensitive to the absence of an inlet channel especially in
the regions close to the step.
In the literature even though there are numerous studies on the flow
over a backwardfacing step, in these studies there are too little
tabulated results. In this study we try to present as many tabulated
results as we can for future references. The reader can find many
more of our tabulated results in [14].
Figure 11 shows streamfunction contours for 100 £ Re £ 3000. This figure exhibit the formation of the
recirculation regions as the Reynolds number increases.
We note
that, the yscale of Figure 11 is expanded in order to be able to
see the details. In Table 5 we have tabulated our numerical results
for the backwardfacing step flow with 2.0 expansion ratio up to
Reynolds number of 3000. Also in Table 6 and Table 7, we have
tabulated the normalized location of the centers of the
recirculating eddies and the streamfunction and vorticity values at
these centers, for future references. In Figure 12 we have plotted
the streamfunction contours for the highest Reynolds number we have
considered, ie. Re=3000, in enlarged view. In this figure, we can
clearly see a recirculating region around x/h=19. As tabulated in
Table 5, this new recirculating region appears in the flow at
Re=2600 and starts to grow as the Reynolds number increase
further. For a better visualization of the flow at Re=3000, in
Figure 13 we have plotted the uvelocity profiles at various
downstream locations. We note that these downstream locations are
shown as dashed lines in Figure 12. In Figure 13, we can clearly see
the stages of the flow developing into a fully developed parabolic
profile towards the outflow boundary.
In Table 5, we see that at Re=3000 a new recirculating region
appears in the flow around x/h=69. Even though it is drawn in
enlarged scale, in Figure 12 it is hard to see this new
recirculating region around x/h=69. In Figure 14, we have plotted
these smaller recirculating regions around x/h=19 and 69 in a more
enlarged scale for a better visualization. In Figure 14, with this
scale we can clearly see the new recirculating region appear in the
flow at Re=3000 around x/h=69.
In Table 8 we have tabulated the length of the main recirculation
region, X1, the detachment and reattachment locations of the first
recirculation region on the upper wall, X2 and X3, and also the
length of this recirculating region on the upper wall, X3X2, at
Reynolds number of 800, together with the same data found in the
literature. The difference in results in Table 8 could be attributed
to especially the different outflow locations and different grid
mesh considered in these studies. We believe that, with a fine grid
mesh of 4251×101 points and an exit boundary located at 300
step heights away from the step, our numerical solutions are more
accurate. In Table 9 we also tabulated the streamfunction (y)
and vorticity (w) values at the center of the main
recirculating region, Eddy X1, and the first recirculating region on
the upper wall, Eddy X3X2, and the locations of eddy centers for
Re=800, together with the results found in the literature. The
results in Table 9 are in good agreement with each other.
Since we have solved the backwardfacing step flow both with 1.942
and 2.0 expansion ratios, in order to see the effect of the
expansion ratio on the flow problem, in Figure 15 we have plotted
the length of the main recirculating region, X1, with respect to the
Reynolds number for both cases. Even though both expansion ratio
numbers, 1.942 and 2.0, are very close to each other in magnitude,
in Figure 15 we can see that up to Reynolds number of 400 the length
X1 of 2.0 expansion ratio is greater than the length X1 of 1.942
expansion ratio. However on the contrary, beyond this Reynolds
number (Re > 400) the X1 length of 1.942 expansion ratio is greater
than the length X1 of 2.0 expansion ratio. It is clear that the step
height has an effect on the flow over backwardfacing step as this
was studied by Thangam and Knight [34,35]. We will
analyze this effect in detail later in Part II of this study.
In Figure 16 we have plotted the tabulated values of the lengths X1,
X2, X3, X4 and X5 in Table 5 as a function of the Reynolds number.
Top figure shows the values as it is. In this top figure, we see
that except around the bifurcation Reynolds numbers at which a new
recirculation region appears in the flow field, the lengths X1, X2,
X3, X4 and X5 change almost linearly with respect to the Reynolds
number. In Figure 16, in the bottom figure we have plotted the same
figure distinguishing this bifurcation regions. In Table 5 we see
that there appears a new recirculating region in the flow at
Reynolds numbers of Re=400, 1700, 2700 and 3000. In Figure 16
bottom figure when we exclude the bifurcation Reynolds number region
around Re=400 and 1700 we can clearly see that the lengths X1, X2,
X3, X4 and X5 behave almost linearly with the Reynolds number. In
the bottom figure the red lines denote the fitted linear lines to
the corresponding points in the figure.
When there appears a new recirculating region in the flow field,
this recirculating region affects the other recirculating region
upstream of it. For example at Re=400 when a new recirculating
region between X2 and X3 appears in the flow field, the upstream
recirculating region X1 starts to grow with a different slope by the
Reynolds number. Also at Re=1700 an other recirculating region
appears between X4 and X5 in the flow field and this mostly affects
the first upstream recirculating region (X2X3), such that beyond
this Reynolds number X2 and X3 starts to grow with a different slope
as the Reynolds number increases furthermore. Although this is not
as clear as the others in Figure 16, we believe that, beyond
Re=1700 this new recirculating region X4X5 affects the slope of
the very upstream recirculating region X1 also, however given the
fact that X4X5 region is far away from the X1 region and also the
Reynolds number is high, the change in the slope of X1 is very
small. We note that, we believe we do not have enough data to
comment on the effects of the X6X7 region that appears in the flow
at Re=2700 and X8X9 region that appears in the flow at Re=3000.
Since the yscale of Figure 16, in which we have plotted the
variation of the location of the detachment and reattachment points
of the recirculating regions with respect to the Reynolds number, is
large, in Figure 17 we decided to plot the variation of the absolute
length of the recirculating regions with respect to the Reynolds
number, that is to say we plot the distance between the reattachment
and detachment points with respect to the Reynolds number. With a
larger yscale, in Figure 17 it is possible to distinguish the
linear regions at Reynolds numbers other than the bifurcation
Reynolds numbers.
We have presented highly accurate numerical solutions of the 2D
steady incompressible backwardfacing step flow obtained using the
efficient numerical method introduced by Erturk et al.
[9]. In our computations the outflow boundary was located
downstream very far away from the step (300 step heights) and also
the inlet boundary was located in an inlet channel very upstream
from the step (20 step heights). Using nonreflecting boundary
conditions at the exit boundary, a fine grid mesh with
4251×101 points and a convergence criteria that is close to
machine accuracy, we were able to obtain numerical solutions up to
very high Reynolds numbers. Detailed results of the backwardfacing
step flow up to Re=3000 are presented. Our results showed that,
for the backwardfacing step flow an inlet channel that is at least
5 step heights long is required for accuracy. Also, the location of
the exit boundary and the outflow boundary condition has an effect
both on the accuracy and on the largest Reynolds number that could
be computed numerically. Our results also indicate that the size of
the recirculating regions grows almost linearly as the Reynolds
number increases.
B.F. Armaly, F. Durst, J.C.F. Pereira, B. Schönung, Experimental
and Theoretical Investigation of BackwardFacing Step Flow,
Journal of Fluid Mechanics127 (1983) 473496.
F. Biagioli, Calculation of Laminar Flows with SecondOrder Schemes
and Collocated Variable Arrangement, Int. J. Numer. Methods
Fluids26 (1998) 887905.
G. Comini, M. Manzan, C. Nonino, Finite Element Solution of the
StreamfunctionVorticity Equations for Incompressible
TwoDimensional Flows, Int. J. Numer. Methods Fluids19
(1994) 513525.
E. Erturk, T.C. Corke, C. Gokcol, Numerical Solutions of 2D Steady
Incompressible Driven Cavity Flow at High Reynolds Numbers,
Int. J. Numer. Methods Fluids48 (2005) 747774.
E. Erturk, C. Gokcol, Fourth Order Compact Formulation of
NavierStokes Equations and Driven Cavity Flow at High Reynolds
Numbers, Int. J. Numer. Methods Fluids50 (2006)
421436.
A. Fortin, M. Jardak, J.J. Gervais, R. Pierre, Localization of Hopf
Bifurcations in Fluid Flow Problems, Int. J. Numer. Methods
Fluids24 (1997) 11851210.
P.M. Gresho, D.K. Gartling, J.R. Torczynski, K.A. Cliffe, K.H.
Winters, T.J. Garratt, A. Spence, J.W. Goodrich, Is the Steady
Viscous Incompressible TwoDimensional Flow Over a BackwardFacing
Step at Re=800 Stable?, Int. J. Numer. Methods Fluids17
(1993) 501541.
G. Jin, M. Braza, A NonReflecting Outlet Boundary Condition for
Incompressible Unsteady NavierStokes Calculations, J. Comp.
Physics107 (1993) 239253.
J. Keskar, D.A. Lyn, Computations Of a Laminar BackwardFacing Step
Flow at Re=800 with a Spectral Domain Decomposition Method,
Int. J. Numer. Methods Fluids29 (1999) 411427.
M. Ramsak, L. Skerget, A Subdomain Boundary Element Method
for HighReynolds Laminar Flow Using Stream FunctionVorticity
Formulation, Int. J. Numer. Methods Fluids46 (2004)
815847.