Numerical Method

While I was working on a project a while ago, I had to solve 2-D Steady Incompressible Navier Stokes equations. I was having difficulty in obtaining a converged solution for the problem especially at high Reynolds numbers. I tried every method explained in CFD books and in journal papers, however I still couldn't get a solution. While I was about to give up on the subject, I decided to find out, at least, why the numerical methods at hand do not work at high Reynolds numbers for steady flows. I was going to do a stability analysis for the most commonly used numerical methods such as ADI, Beam and Warming, and some other methods that use spatial factorization, applied to the governing equations. While I was working on numerical stability I realized that each numerical method has its own positive and negative sides in terms of numerical stability. Then I tried to combine all the positive sides in a numerical scheme excluding the negative sides and I obtained this numerical method. As you will see, this numerical method is obtained step by step. At each step, every different numerical posibilities are first discussed then the posibility that offers the best numerical stability is choosen, and at the next stage the same is continued.

For incompressible two-dimensional and axi-symmetric flows, it is convenient to use the streamfunction ( ) and vorticity ( ) formulation of the Navier-Stokes equations. In non-dimensional form, 2-D steady incompressible N-S equations are given as

(1)

(2)

where, Re is the Reynolds number, and x and y are the Cartesian coordinates. These equations are non-linear, and therefore need to be solved in an iterative manner. We assign pseudo time derivatives to these equations such as

(3)

(4)

and iterate on the variables ( and ) in this pseudo time domain until the solution converges. Since we seek a steady solution through iteration in pseudo time domain, the discretization order in pseudo time will not affect the final solution when steady state is achieved. Therefore we will use a first order discretization () in the pseudo time domain.

Before going into the numerical method let us discuss a few points about numerical stability using some model equations.

Statement 1: In terms of numerical stability, it is better to approximate the dissipation terms on the implicit side (LHS) rather than to approximate them on the explicit side (RHS).

In order to explain this statement let us use a simple equation as the following

(5)

As explained earlier, I will use a first order approximation () for the time derivatives. We can approximate the dissipation term either on the explicit side or on the implicict side. If we prefer to approximate the dissipation term on the explicit side, here is what we will simply get

(6)

If we apply a stability analysis to this equation, we will find that for numerical stability, we sould use a such that

(7)

----------------------------------------------------------

Von Neuman Stability Analysis of Explicit 1-D Diffusion Equation:

We have the following equation

(8)

Finite difference equation will be

(9)

or simplifying we get

(10)

where

(11)

Let

(12)

then

(13)

Substituting these into (10) we get

(14)

Cancelling out some terms we get

(15)

Note that

(16)

With this, equation (15) can be written as

(17)

Let

(18)                   and

then

(19)

Note that,  is the amplification factor and it must be less than 1 for numerical stability. Consider

(20)

If  is greater than 1, then the value of  will increase at every time step which will yield to divergence.

From trigonometric relations we get

(21)

Therefore equation (19) can be written as

(22)

For numerical stability

(23)

then

(24)

or

(25)

hence for numerical stability

(26)

or

(27)

Reference:

Tannehill, Anderson and Pletcher, Computational Fluid Mechanics and Heat Transfer, Taylor & Francis, page 87 (first paragraph)

----------------------------------------------------------

On the other hand, if we go back to beginning (equation 5) and approximate the dissipation term on the implicit side we will obtain the following

(28)

Note that equation (28) can also be written in operator notation as the following

(29)

Again if we apply stability analysis, we find that the above approximation is unconditionally stable.

----------------------------------------------------------

Von Neuman Stability Analysis of Implicit 1-D Diffusion Equation:

We have the following equation

(30)

With operator notation this equation looks like

(31)

Finite difference equation will be

(32)

or

(33)

where

(34)

Let

(35)

Then

(36)

Substituting these into (33) we get

(37)

cancelling out some terms

(38)

Let

(39)                   and

then

(40)

or

(41)

Hence

(42)

Therefore the amplification factor  for any value of . This finite diference scheme is unconditionally stable. Unconditionally stable means that the scheme is stable with out any necessary condition, or in other words it means that the scheme is always stable.

Reference:

Tannehill, Anderson and Pletcher, Computational Fluid Mechanics and Heat Transfer, Taylor & Francis, page 87 (example 3.4)

----------------------------------------------------------

One of the reason why I prefer to use operator notation is that in a way it simplifies the stability analysis. Let us rewrite equation (37)

(43)

Rerrranging we get

(44)

Since  and , substituting we get

(45)

Now let us rewrite equation (31)

(46)

From equations (45) and (46), we can consider the LHS operator as a constant obtained after a Fourier expansion. In this case the operator is equivalent to

(47)

We see that if the dissipation term is approximated on the explicit side we end up with a restriction on the time step for numerical stability, however there was no restriction on the time step when we approximate the dissipation term on the implicit side.