At this point after using a first order approximation () for the nonlinear terms our equations become as the following (equations 149 and 150)
One important feature of these equations is that, the equations are solved individually. That means we first solve the streamfunction equation (167) and then after this we solve the vorticity equation (168). And, most importantly, both equations are unconditionally stable.
When we finite difference equations (167) and (168) using three point central differencing, the operators on the LHS will form a coefficient matrix and we will obtain a mathematical system as the following
where each line in the coefficient matrix denotes diagonal elements which are non zero. Except these diagonal elements, the rest of the elements in the coefficient matrix are all equal to zero.
At this point we have two unconditionally stable equations. This is really what we want. On the other hand these equations are in fully implicit form, that is, each equation requires the solution of a large banded coefficient matrix, which is not computationally very efficient. Previously we have only worried about the numerical stability, however the equations we obtain (equations 167 and 168) are not practical to solve especially when we use large number of grid points. Now we have to worry about the computational efficiency also. Therefore instead of solving equations (167) and (168) in implicit form, more efficiently, we spatially factorize the equations such that
Here we split the single LHS operators into two operators where each only contains derivatives in one direction. Note that equations (169) and (170) are unconditionally stable also. The only difference between equations (167 and 168) and (169 and 170) are the terms produced by the factorization. The advantage of this process is that the equations now require the solution of tridiagonal systems which is numerically more efficient.
Now we have to discuss some issues about spatial factorization. The argument about factorization is that the terms due to factorization are equal or less than the order in the truncation error in time stepping and therefore have a negligible affect on the solution. Let us examine this argument. As an example for simplicity we consider the streamfunction equation (167) only. The pseudo time derivative is approximated by an implicit Euler time step such that
In this approximation the leading truncation error is (). At first it appears that the terms due to factorization have the same order as the leading truncation error and therefore can be neglected. However our equations are steady, and we are seeking a steady solution by iterating in this pseudo time domain. At the steady state the pseudo time derivative is
Therefore at the steady state the truncation errors in equation (171) become zero. However the terms due to factorization will not become zero and will remain in the final solution. To illustrate this further, equation (169) can be written in explicit form as the following
As the solution converges to a steady state, converges to and they cancel out each other. Therefore at convergence equation (174) appears as the following
Note that the streamfunction equation is
In equation (175) there is no guarantee that the RHS, which is the consequence of the factorization, will be equal to zero. Only a very small time step () will ensure that the RHS of equation (175) will be close to zero. However when Dt is very small then the convergence will be very slow. In addition, the cross derivative term may be significant and may be (O ) or higher in many flow problems. Therefore there is problem in equations (169) and (172). To overcome this problem we add similar cross derivative terms at the previous time step to RHS of equations and obtain the following equations
If we analyse, equation (177) would give us the following
As the solution converges to a steady state, converges to and also converges to , so that they cancel out and the final solution converges to the correct physical form where
We would like to emphasize that the added terms to RHS in equations (177) and (178) are NOT artificial numerical dissipation terms. Such terms are often used to stabilize the numerical schemes so that convergence can be achieved. Here, the added terms on the RHS are not meant to stabilize the numerical scheme. The only reason for this approach is to kill the terms introduced by factorization so that the equations are the correct physical representation at the steady state.
One of the most popular spatially factorized schemes is undoubtedly the Beam & Warming method. In that method, the equations are formulated in delta () form. In a delta formulation, there will be no second order terms due to the factorization in the solution at the steady state.
At this point one may suggest to use a delta () form in our streamfunction and vorticity equations (167 and 168) and then factorize the resulting equations so that we would not have any second order terms at the steady state. If we have done this, we would obtain the following equations
where and . Note that when the solution converges and converges to and respectively, therefore and converges to zero. For this reason in this formulation the terms produced by factorization have no effect on the solution at the steady state, since these terms together with every term on the LHS in equations (181) and (182) become homogeneous at the steady state. This is good but let us analyse these two equations in terms of stability. In previous pages, we have worked on the stability of several model equations and as a result we deduced two rules for our numerical method for numerically stability. Let us see if we violate any of the two rules in equations (181) and (182). In equation (182) we have dissipation terms on the RHS. According to Rule 1, even though we do not want any dissipation terms on RHS this is still OK. In equation (182) we have dissipation and most importantly convection terms on the RHS. According to Rule 2, we definitely do not want to have ( order) convection terms on the RHS. Therefore by using a delta formulation () we have violated the rules we started with in the first place.
Equations (181) and (182) is nothing but the Beam & Warming method applied to the streamfunction and vorticity equation. (with and for those who are familiar with this method). Even though Beam & Warming method is always said to be unconditionally stable, it requires artificial numerical dissipation for stability. This is contradictory. One thing I would like to note, I have never seen a study that uses Beam & Warming method applied to N-S equations, that does not need artificial numerical dissipation terms for numerical stability. I have written a code that solves (181) and (182) as it is shown above (with out any additional artificial dissipation), and I had very difficulty in getting convergence even at small Reynolds numbers. When I add numerical dissipation to equations (181) and (182) I was again having difficulty in even moderate Reynolds numbers. The Beam & Warming method applied to streamfunction and vorticity equations was not successful.
This was the reason why I get into developing this numerical method, cause I wanted to find a numerical scheme that is REALLY unconditionally stable that DOES NOT require any artificial numerical dissipation terms for stability. I was never comfortable in using artificial numerical dissipation because as the name implies they are ARTIFICIAL. When artificial numerical dissipation is used then the solution obtained is not actually the solution of the original equations.
Following the rules we deduced and then factorizing the equations we have obtained unconditionaly stable equations but the terms produced by the factorization remained in the converged steady state solution. Then to overcome this problem we have added the same amount of terms to RHS also to kill the terms produced by the factorization. The resulting equations converge to the correct physical form of the N-S equations. Now we have to look at the stability and see if we have violated any of the rules we have. Let us rewrite the final form of our resulting equations
or writing the RHS in explicit form these equations look like the following
We see that we have some cross derivative terms on the RHS. These terms will have, of course, an effect on the numerical stability. But one thing is very important, these terms are second order (). This fact is very important and this is the main difference between this numerical method and the Beam & Warming method. If is small than will be even smaller, therefore the effect of these terms may be actually very small. The only term that may be dangerous is the last term with . As Reynolds number increase the effect of this last term may be significant on numerical stability. To apply a Von Neumann stability analysis to equations (183) and (184) is very difficult. Anyhow when I wrote the code that solves equations (183) and (184), I saw that as I experiment with the code, the scheme has a very good stability characteristics. In fact the numerical stability of the equations (183) and (184) had better numerical stability than the delta formulation shown in equation (181) and (182) even with additional fourth order artificial numerical dissipation terms on the RHS and second order artificial numerical dissipation terms on the LHS as it was suggested by Beam & Warming.
Therefore as conclusion, the numerical method we developed for the solution of 2-D Steady Incompressible Navier-Stokes equations in streamfunction and vorticity formulation takes the final form as the following
The solution methodology of these two equations is quite simple. First we solve for equation (187). In order to solve this equation first we define a new variable, f, such that
then equation (187) can be written as
In equation (190) the only unknown is f. We can easily solve this equation for this new variable f. After this, when we obtain the values of f, then we can easily solve equation (189) for . We note that both equations (190 and 189) require the solution of tridiagonal systems which is straightforward with the Thomas algorithm.
After solving the streamfunction equation we, then, solve for the vorticity equation (188). In order to do this, as we did in solving the streamfunction equation, we define a new variable, g, such that
then equation (188) can be written as
In equation (192) the only unknown is g. After solving this equation for this new variable g, then we solve equation (191) for . Again, we note that, both equations (192 and 191) require the solution of tridiagonal systems.
The numerical method we described here is quite efficient and stable even at very high Reynolds numbers. Most importantly, since this method require the solution of tridiagonal systems, which is computationally very efficient and easy to solve, this method allows us to use very large number of grid points compared to a fully implicit method. This numerical method is used successfully in the following studies: