Our goal is to solve the streamfunction ( ) and vorticity ( ) formulation of the 2-D Steady Incompressible Navier Stokes equations numerically. In non-dimensional form, the equations are given as
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. For this we assign pseudo time derivatives to these equations (98) and (99) such as
Since a steady solution is desired, the final solution will not be affected by the order of the pseudo time derivative. We use a first order () accurate finite difference approximation for these pseudo time derivatives. Following Rule 1 and Rule 2, we will approximate the dissipation and convection terms on the implicit side such that
At this point of developing our numerical method, we have to ask ourselves two questions. The first one is, how should we linearize the nonlinear convective terms in equation (103), and related with this, should the way we linearize these terms affect the numerical stability and also if the answer is yes then in which way (positive or negative). Second one is, which side (explicit or implicit) should we approximate the vorticity term in equation (102).
As we analyse these two questions, we will see that the answer to second question is not really very important however the first question is a very important question and the way we linearize the non linear terms will have a HUGE effect on numerical stability.
Let us analyse the choices we have for linearization the nonlinear terms and their effects on the numerical stability.
Now the question is how to approximate the following nonlinear terms.
We have to linearize these terms somehow. We can use different linearizations but we have to ask ourselves how the numerical stability of the equations are affected by the approximation of these non-linear terms.
Second order () approximation:
First let us use a second order () Newton linearization for these terms and do a stability analysis. For simplicity let use have an example nonlinear term such as
If we expand each term in Taylor series such as
then the nonlinear term will be
multiplying we get
If we substitute (111) and (112) in (110) then we get
or simplifying we have
Note that this approximation is second order accurate () in time
Therefore equations (104) and (105) will be
If we apply this second order linearization to equation (103) and move the vorticity term in equation (102) to LHS, the equations become as the following
We write these equations in operator notation, such as
Let us do a stability analysis. If we choose
and also substitute these into equation (117) and (118), we can easily obtain that the operator in equation (117) is equivalent to the following
where and and
and the operators in equation (118) are also equivalent to the following
where , , , and
and also and
Let A, B, C and D denote
If we substitute these for the LHS operators in equations (117) and (118) we simply get
Note that in equation (126) we didnot include the following terms on the RHS in equation (118) for a simpler numerical analysis.
These terms are convection terms, and they are on the RHS. After we obtain the numerical stability, we should know that these convection terms on the RHS will decrease the numerical stability of the scheme.
If we solve equations (125) and (126) for and we get
Now, examining these two equations (128) and (129) we see that is a function of and and also is a function of and in numerical stability sense. Therefore for absolute numerical stability, in equation (128) the coefficient of have to be greater than the sum of coefficients of and and also in equation (129) the coefficient of have to be greater than the sum of coefficients of and . Hence the following conditions must be satisfied for numerical stability
If we analyse the terms in these inequalities, from equation (124) we find that
Note that during an iteration process the value of C can be positive or negative with any magnitude. If it is positive, it is good for numerical stability since B is negative, however if it is negative then it is bad. If the sign of C changes, the stability criterias in equations (130) and (131) can not always be satisfied, therefore the scheme may sometimes be numerically unstable depending on C.
Let us go back to equation (115) and see that what would happen if we have approximated the vorticity term in this equation on the explicit side (RHS), then the streamfunction equation (102) would become
Then one can show that equation (133) would become
and also equation (134) would become
Rewriting the vorticity equation (126)
Using equation (135) and equation (136) and solving for and we get
We see from these two equations (137) and (138) that is a function of and and also is a function of and in numerical stability sense. Therefore for absolute numerical stability the following conditions must be satisfied
Again if the sign of C is positive it is good, since B is negative, however if the sign is negative it is bad for numerical stability. Therefore, depending on the sign and magnitude of C, this scheme may sometimes be numerically unstable.
As a result we see that if we use a second order () Newton linearization for the nonlinear terms, the resulting numerical schemes we have do not quarantee numerical stability all the time and might be unstable. The worst thing is that we can not switch from numerical unstability to stability by changing any of the free parameters ( and ). The sign and magnitude of C dictates the numerical stability. As seen in equation (132) the value of C is dependent on and and we have no control on these parameters during the iterations. The only thing we can do is to change the initial contiditons we assing for and , therefore have these variables follow a different path to the steady solution and therefore have different values of and during iterations. And if we are lucky we can find such an initial condition that makes the scheme remain stable until convergence, but this almost never happens. I do not like the idea of the initial guess having an effect on numerical stability. Our aim should be to have a stable numerical scheme regardless of the initial condition.
With second order linearization ( ) of the nonlinear terms in vorticity equation, we did not obtain absolute numerical stability. In the next pages we will look for other options that offers better numerical stability.