# Numerical Solutions of 2-D Steady IncompressibleFlow Over a Backward-Facing Step,Part I: High Reynolds Number Solutions

### Abstract

Numerical solutions of 2-D laminar flow over a backward-facing step at high Reynolds numbers are presented. The governing 2-D steady incompressible Navier-Stokes 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 backward-facing step are compared with experimental and numerical results found in the literature.

### 1  INTRODUCTION

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 backward-facing 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 2-D steady incompressible flow over a backward-facing 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 backward-facing 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 backward-facing step for Re=1000. Similarly, Cruchaga [7] solved the steady backward-facing 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 backward-facing step. For example, Fortin et al. [15] have studied the stability of the 2-D steady incompressible flow over a backward-facing 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 backward-facing 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 backward-facing 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 two-dimensional perturbations and moreover shows no evidence of any nearby two-dimensional bifurcation. These two studies, [2,15], are important in the sense that their two-dimensional computational stability analysis show that the flow over a backward-facing step is temporally stable at high Reynolds numbers, that is to say there exist steady solutions of the flow over a backward-facing step at high Reynolds numbers.

Erturk et al. [9] introduced an efficient, fast and stable numerical formulation for the steady incompressible Navier-Stokes 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 non-orthogonal grid mesh at extreme skew angles ([13]).

In this study, we will present numerical solutions of 2-D steady incompressible flow over a backward-facing step at high Reynolds numbers. The 2-D steady incompressible Navier-Stokes 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 backward-facing 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 backward-facing 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 backward-facing step flow with both 1.942 and 2.0 expansion ratios.

While the flow over a backward-facing 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 backward-facing 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 backward-facing 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 backward-facing 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 backward-facing 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 backward-facing 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 backward-facing 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 Rem=[(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 backward-facing 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, non-reflecting 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.

### 2  NUMERICAL METHOD

We consider the 2-D steady incompressible Navier-Stokes equations in streamfunction (y) and vorticity (w) formulation. In non-dimensional form, they are given as

 ¶2y ¶x2 + ¶2y ¶y2 = -w
(1)

 1 Re æè ¶2 w ¶x2 + ¶2 w ¶y2 öø = ¶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=2hi 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 Navier-Stokes 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  dxx - Dt  dyy öø yn+1 = yn + Dt  wn
(3)

 æè 1 - Dt 1 Re dxx - Dt 1 Re dyy + Dt  yyn  dx- Dt  yxn  dy öø wn+1 = wn
(4)

where subscripts denote derivatives and also dx and dy denote the first derivative and similarly dxx and dyy denote the second derivative finite difference operators in x- and y-directions respectively, for example

 dx q = qi+1,j-qi-1,j 2Dx + O (Dx2)
 dxx q = qi+1,j-2qi,j+qi-1,j Dx2 + O (Dx2)
(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  dxx öø æè 1 - Dt  dyy öø yn+1 = yn + Dt  wn
(6)

 æè 1 - Dt 1 Re dxx + Dt  yyn dx öø æè 1 - Dt 1 Re dyy -Dt  yxn  dy öø wn+1 = wn
(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 Dt2 terms into the left hand side of the equations and these Dt2 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 Dt2 terms to the right hand side of the equations such that the final form of the finite difference equations become the following

 æè 1 - Dt  dxx öø æè 1 - Dt  dyy öø yn+1 = yn + Dt  wn + æè Dt  dxx öø æè Dt  dyy öø yn
(8)

 æè 1 - Dt 1 Re dxx + Dt  yyn dx öø æè 1 - Dt 1 Re dyy -Dt  yxn  dy öø wn+1 = wn

 + æè Dt 1 Re dxx - Dt  yyn  dx öø æè Dt 1 Re dyy + Dt  yxn  dy öø wn
(9)

We note that, at the steady state the Dt2 terms on the right hand side of the equations cancel out the Dt2 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 two-stage time-level updating. For the streamfunction equation (8), the variable f is introduced such that

 æè 1 - Dt  dyy öø yn+1 = f
(10)

and

 æè 1 - Dt  dxx öø f = yn + Dt  wn + æè Dt  dxx öø æè Dt  dyy öø yn
(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 dyy - Dt  yxn  dy öø wn+1 = g
(12)

and

 æè 1 - Dt 1 Re dxx + Dt  yyn dx öø g = wn + æè Dt 1 Re dxx - Dt  yyn  dx öø æè Dt 1 Re dyy + Dt  yxn  dy öø wn
(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).

### 3  THE BACKWARD-FACING STEP FLOW

The schematics of the backward-facing 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 L1 in Figure 1a. In this inlet channel, L1, 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 L2 in Figure 1a. In order to have high accuracy in the vicinity of the step, in x-direction from the step to a distance of 100 step heights we used 2500 uniform fine grids as shown as L3 in Figure 1b. From 100 step heights distance to the exit boundary, also as shown as L4 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 L4, the grids are stretched smoothly starting from the constant grid spacing used in L3 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=L3 + L4 (b+1)-(b-1)[(b+1)/(b-1)]1-[`x] [(b+1)/(b-1)]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 y-direction, we have used 101 uniform grids as it is shown in Figure 1b.

### 3.1  Inlet boundary condition

Barton [3] studied the entrance effect for flow over a backward-facing 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 backward-facing step found in the literature, this study utilizes the largest inlet channel, L1. Cruchaga [7] solved the backward-facing 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.

### 3.1  Exit boundary condition

At the exit boundary we used a non-reflecting 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 non-reflecting boundary condition concept is first introduced in the name of Äbsorbing Boundary Condition", and also to Jin and Braza [21] for a review of non-reflecting 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/x2) 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(i-ibuf) (imax-ibuf) öø
(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 2-D steady incompressible Navier-Stokes equations should approach to a fully developed Plane Poiseuille flow between parallel plates, such that the velocity profile is parabolic and also the second x-derivatives of streamfunction and vorticity variables should be equal to zero, ie. [(2 y)/(x2)]=0 and [(2 w)/(x2)]=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 non-reflecting 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)/(x2)]=[(v)/(x)]=0 and also [(2w)/(x2)]=[(3 u)/(x2y)] - [(3 v)/(x3)] = 0 where u and v are the velocity components in x- and y-directions respectively.

After the step, physically the flow needs some streamwise distance to adjust and become fully developed. The non-reflecting boundary condition we apply is equivalent to assuming the flow to be fully developed. Our extensive numerical tests have shown that, using the non-reflecting 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 x-direction, 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 non-reflecting 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.

### 3.1  Wall boundary conditions

The vorticity value at the wall is calculated using Jensen's formula (see Tannehill, Anderson and Pletcher [36])

 w0= 3.5 y0-4 y1 + 0.5 y2 Dy2
(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 y-direction. 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.

### 4  RESULTS AND DISCUSSIONS

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

RES1y  = max æè
 êê yi-1,jn+1-2yi,jn+1+yi+1,jn+1 Dx2
 + yi,j-1n+1-2yi,jn+1+yi,j+1n+1 Dy2 +wi,jn+1 êê öø
RES1w  = max æè
 êê 1 Re wi-1,jn+1-2 wi,jn+1+wi+1,jn+1 Dx2
 + 1 Re wi,j-1n+1-2wi,jn+1 +wi,j+1n+1 Dy2
 - yi,j+1n+1-yi,j-1n+1 2 Dy wi+1,jn+1-wi-1,jn+1 2 Dx
 + yi+1,jn+1-yi-1,jn+1 2 Dx wi,j+1n+1-wi,j-1n+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

RES2y =  max æè êê yi,jn+1 - yi,jn êê öø
RES2w =  max æè êê wi,jn+1 - wi,jn êê öø
(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

RES3y =  max æè êê yi,jn+1 - yi,jn yi,jn êê öø
RES3w =  max æè êê wi,jn+1 - wi,jn wi,jn êê öø
(20)

RES3y =  max æè êê yi,jn+1 - yi,jn yi,jn êê öø
RES3w =  max æè êê wi,jn+1 - wi,jn wi,jn êê öø
(20)

In our calculations, for all Reynolds numbers we considered that convergence was achieved when both RES1y £ 10-10 and RES1w £ 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 RES2y £ 10-17 and RES2w £ 10-15, that means the streamfunction and vorticity variables are accurate to 16th and 14th 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 RES3y £ 10-14 and RES3w £ 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 backward-facing 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 x-direction 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 2-D backward-facing 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 backward-facing 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 non-reflecting 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, u-velocity 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 backward-facing 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 u-velocity 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 Navier-Stokes 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 "down-wash". This in fact is what we observe at small Reynolds numbers as shown in Figure 3. From Figure 3, we conclude that in backward-facing 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 backward-facing 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 backward-facing step with 1.942 expansion ratio. Using the numerical procedure described above, we have solved the steady 2-D Navier Stokes equations up to Re £ 1500. Figure 4, 5, 6, 7 and 8 show the u-velocity profiles at several x-locations 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 x-locations 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 x-locations close to the step, and the agreement is moderate at far downstream x-locations. 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 3-D effects observed in experiments. Figures 4, 5, 6, 7 and 8 show that our 2-D 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 backward-facing step flow and Table 1 tabulates our numerical solutions for the backward-facing 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 backward-facing 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 backward-facing step with a 2.0 expansion ratio as it is done in most of the numerical studies. For the backward-facing step with expansion ratio of 2.0, we have solved the steady incompressible Navier-Stokes 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 backward-facing 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 backward-facing 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 backward-facing 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 y-scale 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 backward-facing 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 u-velocity 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, X3-X2, 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 X3-X2, 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 backward-facing 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 backward-facing 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 (X2-X3), 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 X4-X5 affects the slope of the very upstream recirculating region X1 also, however given the fact that X4-X5 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 X6-X7 region that appears in the flow at Re=2700 and X8-X9 region that appears in the flow at Re=3000. Since the y-scale 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 y-scale, in Figure 17 it is possible to distinguish the linear regions at Reynolds numbers other than the bifurcation Reynolds numbers.

### 5  CONLUSIONS

We have presented highly accurate numerical solutions of the 2-D steady incompressible backward-facing 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 non-reflecting 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 backward-facing step flow up to Re=3000 are presented. Our results showed that, for the backward-facing 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.

### References

[1]
B.F. Armaly, F. Durst, J.C.F. Pereira, B. Schönung, Experimental and Theoretical Investigation of Backward-Facing Step Flow, Journal of Fluid Mechanics 127 (1983) 473-496.

[2]
D. Barkley, M.G.M. Gomes, R.D. Henderson, Three-Dimensional Instability in Flow Over a Backward-Facing Step, J. Fluid Mech. 473 (2002) 167-190.

[3]
I. E. Barton, The Entrance Effect of Laminar Flow Over a Backward-Facing Step Geometry, Int. J. Numer. Methods Fluids 25 (1997) 633-644.

[4]
I. E. Barton, Improved Laminar Predictions Using a Stabilised Time-Dependent Simple Scheme, Int. J. Numer. Methods Fluids 28 (1998) 841-857.

[5]
F. Biagioli, Calculation of Laminar Flows with Second-Order Schemes and Collocated Variable Arrangement, Int. J. Numer. Methods Fluids 26 (1998) 887-905.
[6]
G. Comini, M. Manzan, C. Nonino, Finite Element Solution of the Streamfunction-Vorticity Equations for Incompressible Two-Dimensional Flows, Int. J. Numer. Methods Fluids 19 (1994) 513-525.

[7]
M. A. Cruchaga, A Study of the Backward-Facing Step Problem Using a Generalized Streamline Formulation, Commun. Numer. Meth. Engng. 14 (1998) 697-708.

[8]
B. Engquist, A. Majda, Absorbing Boundary Conditions for the Numerical Simulation of Waves, Math. Comput. 31 (1977) 629-651.

[9]
E. Erturk, T.C. Corke, C. Gokcol, Numerical Solutions of 2-D Steady Incompressible Driven Cavity Flow at High Reynolds Numbers, Int. J. Numer. Methods Fluids 48 (2005) 747-774.

[10]
E. Erturk, C. Gokcol, Fourth Order Compact Formulation of Navier-Stokes Equations and Driven Cavity Flow at High Reynolds Numbers, Int. J. Numer. Methods Fluids 50 (2006) 421-436.

[11]
E. Erturk, T.C. Corke, Boundary Layer Leading Edge Receptivity to Sound at Incidence Angles, Journal of Fluid Mechanics 444 (2001) 383-407.

[12]
E. Erturk, O.M. Haddad, T.C. Corke, Laminar Incompressible Flow Past Parabolic Bodies at Angles of Attack, AIAA J. 42 (2004) 2254-2265.

[13]
E. Erturk, B. Dursun, Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Skewed Cavity, ZAMM · Z. Angew. Math. Mech. 87 (2007) 377-392.

[14]
E. Erturk, Gebze Institute of Technology, 2007, < http://www.cavityflow.com >

[15]
A. Fortin, M. Jardak, J.J. Gervais, R. Pierre, Localization of Hopf Bifurcations in Fluid Flow Problems, Int. J. Numer. Methods Fluids 24 (1997) 1185-1210.

[16]
D.K. Gartling, A Test Problem For Outflow Boundary Conditions - Flow Over A Backward-Facing Step, Int. J. Numer. Methods Fluids 11 (1990) 953-967.

[17]
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 Two-Dimensional Flow Over a Backward-Facing Step at Re=800 Stable?, Int. J. Numer. Methods Fluids 17 (1993) 501-541.

[18]
M.M. Grigoriev, G.F. Dargush, A Poly-Region Boundary Element Method for Incompressible Viscous Fluid Flows, Int. J. Numer. Meth. Engng. 46 (1999) 1127-1158.

[19]
G. Guj, F. Stella, Numerical Solutions of High-Re Recirculating Flows in Vorticity-Velocity Form, Int. J. Numer. Methods Fluids 8 (1988) 405-416.

[20]
O.M. Haddad, T.C. Corke, Boundary layer receptivity to free-stream sound on parabolic bodies, Journal of Fluid Mechanics 368 (1998) 1-26.

[21]
G. Jin, M. Braza, A Non-Reflecting Outlet Boundary Condition for Incompressible Unsteady Navier-Stokes Calculations, J. Comp. Physics 107 (1993) 239-253.

[22]
J. Keskar, D.A. Lyn, Computations Of a Laminar Backward-Facing Step Flow at Re=800 with a Spectral Domain Decomposition Method, Int. J. Numer. Methods Fluids 29 (1999) 411-427.

[23]
J. Kim, P. Moin, Application of a Fractional-Step Method to Incompressible Navier-Stokes Equations, J. Comp. Physics 59 (1985) 308-323.

[24]
T. Lee, D. Mateescu, Experimental and Numerical Investigation of 2-D Backward-Facing Step Flow, Journal of Fluids and Structures 12 (1998) 703-716.

[25]
J.M. Leone Jr., Open Boundary Condition Symposium Benchmark Solution: Stratified Flow Over a Backward-Facing Step, Int. J. Numer. Methods Fluids 11 (1990) 969-984.

[26]
C. Liu, Z. Lin, High Order Finite Difference and Multigrid Methods for Spatially-Evolving Instability, J. Comp. Physics 106 (1993) 92-100.

[27]
J.H. Morrison, M. Napolitano, Efficient Solutions of Two-Dimensional Incompressible Steady Viscous Flows. Computers & Fluids 16 (1988) 119-132.

[28]
T. C. Papanastasiou, N. Malamataris, K. Ellwood, A New Outflow Boundary Condition, Int. J. Numer. Methods Fluids 14 (1992) 587-608.

[29]
C.-Y. Perng, Ph.D. dissertation, Dept. Mech. Eng., Stanford University, 1990 (unpublished).

[30]
M. Ramsak, L. Skerget, A Subdomain Boundary Element Method for High-Reynolds Laminar Flow Using Stream Function-Vorticity Formulation, Int. J. Numer. Methods Fluids 46 (2004) 815-847.

[31]
S.E. Rogers, D. Kwak, An Upwind Differencing Scheme for the Incompressible Navier-Stokes Equations,Appl. Numer. Math. 8 (1991) 43-64.

[32]
R.L. Sani and P.M. Gresho, Resume and remarks on the open boundary condition minisymposium, Int. J. Numer. Methods Fluids 18 (1994) 983-1008.

[33]
T.W.H. Sheu, S.F. Tsai, Consistent Petrov–Galerkin Finite Element Simulation of Channel Flows, Int. J. Numer. Methods Fluids 31 (1999) 1297-1310.

[34]
S. Thangam, D.D. Knight, Effect of Stepheight on the Separated Flow Past a Backward Facing Step, Phys. Fluids A 1 (1989) 604-606.

[35]
S. Thangam, D.D. Knight, A Computational Scheme in Generalized Coordinates for Viscous Incompressible Flows, Computers & Fluids 18 (1990) 317-327.

[36]
J.C. Tennehill, D.A. Anderson, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer (2nd edn), Taylor & Francis, 1997.

[37]
M.M.T. Wang, T.W.H. Sheu, Implementation of a Free Boundary Condition to Navier-Stokes Equations, Internat. J. Numer. Methods Heat Fluid Flow 7 (1997) 95-111

[38]
E. Weinan, L. Jian-Guo, Vorticity Boundary Condition and Related Issues for Finite Difference Schemes, J. Comp. Physics 124 (1996) 368-382.

[39]
H.C. Yee, J.R. Torczynski, S.A. Morton, M.R. Visbal, P.K. Sweby, On Spurious Behavior of CFD Simulations, Int. J. Numer. Methods Fluids 30 (1999) 675-711.

[40]
Y. Zang, R.L. Street, J.R. Koseff, A Non-staggered Grid, Fractional Step Method for Time-Dependent Incompressible Navier-Stokes Equations in Curvilinear Coordinates, J. Comp. Physics 114 (1994) 18-33.

File translated from TEX by TTH, version 3.63.