Application of the Numerical Method to Cavity Flow

The following figure shows the schematics of the lid driven cavity flow. The top wall is moving with a speed of U.

The flow inside a driven cavity is governed by the 2-D Steady Incompressible Navier-Stokes equations. These equations in streamfunction and vorticity formulation are given as

(193)

(194)

As explained in previous pages, the numerical method that Erturk et. al. have used to solve these equations are given as

(195)

(196)

First, initial values are assigned to the streamfunction and vorticity variables. We preferred to use a homogeneous initial guess such that

(197)                      ,

In order to decrease the number of iterations for convergence, if a Reynolds number solution is already obtained then this solution could be used for a different Reynolds number as an initial guess, instead of starting the iterations from a homogeneous solution. This numerical method converges to a unique solution no matter what the initial conditions are. Therefore any solution can be used as an initial condition.

After assigning the initial guesses to the variables we solve the streamfunction equation (195). This equation (195) forms the following mathematical system

(198)

where the LHS matrices are defined as

(199)

such that

(200)                      ,          ,

and

(201)

such that

(202)                      ,          ,

and the RHS vector is

(203)

where  and  and also I is the identity matrix.

In order to solve equation (195), a new variable f is introduced. First the following equation is solved for f at every grid point

(204)

Note that the LHS operator has only x-derivatives. Then equation (204) forms the following tridiagonal system

(205)

In order to solve this equation we need boundary condition values of f at the left and right wall boundaries of the cavity. The variable f is defined as

(206)

Therefore using equation (206), the value of  on the left wall and  on the right wall is calculated as

(207)

Since on the wall boundaries the value of  is equal to zero, the value of f on the left and right wall becomes zero, such as

(208)                  ,

With using these boundary conditions, the system in equation (205) is solved using Thomas algorithm easily. After we obtain the values of f at every interior grid point, then we solve equation 206. This equation forms the following tridiagonal system

(209)

With using the value of  at the upper and lower wall boundaries being equal to zero, the system in equation (209) is solved with using Thomas algorithm.

At this point we have updated the values of  at every grid point. Therefore we can calculate the values of  on the wall boundaries. We have used Thom's method such that

(210)

where the subscript 0 denotes the points on the wall and 1 denotes the points adjacent to wall, and  is the velocity of the wall such that it is equal to 1 at the upper wall and is equal to zero at the other three walls. The equation (210) is used to calculate the values of vorticity on the walls.

Before going into solving the vorticity equation, let us talk about how we can discretize the convection terms. We have used three point second order () central difference for the second derivatives (dissipation terms), such as

(211)

For the first derivatives we can either use three point central difference or we can use an upwind difference with deferred correction. We can two different ways of upwind difference with deferred correction.

Ø     We can use three point central difference

(212)

Ø     We can use the discretization suggested by Rubin & Khosla (Computers and Fluids, Vol 9, pp 138-158, 1981)

(213)

where  if  and  if . Note that as the solution converges to steady state where approaches to , the above solution is identical to three point central differencing.

Ø     We can use the following discretization

(214)

where  if  and  if . Again note that as the solution converges to steady state, the above solution is also identical to three point central differencing.

We will formulate our equations considering the three discretizations.

Now we solve the vorticity equation (196). Equation (196) forms the following mathematical system

(215)

where the LHS matrices are defined as

(216)

such that

(217)                      ,          ,

and

(218)

such that

(219)                      ,          ,

and the RHS vector is

(220)

where  and  and also I is the identity matrix.

Ø     If we use three point central differencing for the first order derivatives then

(221)

Ø     If we use the discretization suggested by Rubin & Khosla then

(222)

where the values of  and  are equal to

(223)

(224)

Ø     If we use the third type of discretization then

(225)

where again the values of  and  are equal to

(226)

(227)

In order to solve equation (196), a new variable g is introduced. First the following equation is solved for g at every grid point

(228)

This equation (228) forms the following mathematical system

(229)

In order to solve this equation we need boundary condition values of g at the left and right wall boundaries of the cavity. The variable g is defined as

(230)

Therefore using equation (230), the value of  on the left wall and  on the right wall is calculated as

(231)

Note that on the left and right wall . Using these boundary conditions, the system in equation (229) is solved using Thomas algorithm and we obtain the values of g at every interior grid point. After this we solve the system in equation (230) such that

(232)

Using the values of  on the upper and lower wall, we calculated using equation (210), as boundary conditions, we easily solve the system in equation (232) using Thomas algorithm.

At this point we have updated the values of  at every point, and one iteration cycle is finished. Now we check the residuals to see if the numerical solution has converged to steady state. If the residuals are not below a certain level that we set, then we continue with the next iteration cycle and solve equation (204).

My extension tests on the numerical method have showed that using a three point central difference for the first derivatives works the best. The third way of discretization is the worst, the convergence is slower and some times the code diverges. The second way of discretization (suggested by Rubin & Khosla) works also great, moreover with this way of discretization sometimes the code converges faster than the regular central differencing

Residuals

At the end of each iteration, we have to calculate the residuals in order to see if the numerical solution has converged to a steady state or not. I have used three residuals for each variable (streamfunction and vorticity). These residuals help to decide whether or not the convergence is achieved.

The first residual is simply the residual of the streamfunction and vorticity equations (equations 193 and 194), such that

(233)

(234)

We are solving the streamfunction and vorticity as shown in equation (193 and 194). In the limit, ie. at the steady state,  and  should be equal to zero. The definition residuals in  and  is rather uncommon in CFD, however I think this is the most important residual we should look at.

We have discretized the mathematical differential equations (193) and (194) using three point central difference, therefore we are seeking a solution that satisfies the LHS of equations (223) and (224) at every grid point. The residuals  and  actually shows how far our numerical solution is from the mathematical  solution (for a three point second order () discretization). Therefore, the residuals defined in (223) and (224) have a mathematical significance. In my computations I considered that the convergence to steady state is achieved when both  and  are less than , that is to say when the magnitude of the residuals of streamfunction and vorticity equation is less than  every grid point in the domain. Such a low value is actually not necessary, however it certainly shows the ability of the presented numerical method to converge to very low residuals.

The second residual we monitor during the iterations is defined as

(235)

(236)

This type of residual definitions is very common in CFD. It simply shows how much the variables change in magnitude in an iteration step. In our numerical solutions, when  and  are less than , the values of  and  are in the order of  and  respectively. Again such low values demonstrate the efficiency of the numerical method.

The third residual we monitor during the iterations is defined as

(237)

(238)

These definitions of residuals are similar to previous ones. These show the percentage of the change of magnitude in a variable in an iteration step. WATCH OUT, if you start the iterations from a homogeneous initial guess, you have to skip calculating these residuals ( and ) in the first iteration step. However you do not have to check the value of  and  at every grid to see if they are equal to zero, becouse it is almost impossible for these variables to get an exact value of zero (0.000...000 with double precision) during iterations numerically. They may get a very very small value but not zero exactly. Also in those grid points, the numerator of  and  is even smaller than the denominator. If you are not comfortable do not calculate  and  in your codes. In my solutions when convergence is achieved to steady state the values of  and  are in the order of  and  respectively. This means that at the steady state, at every grid point the streamfunction and vorticity values are changing less than 0.0000000001% of their values between two iteration steps. Again such low values demonstrate the efficiency of the numerical method.