18.4 Differencing Schemes

Having obtained expressions for the discretised form of each of the terms in the conservation equation in the previous section, the discretised form of the full equation is obtained by simply adding together all these contributions. If only the convection and diffusion terms are considered, and using arithmetic averaging in the evaluation of the face value of $\phi$ in the convection term, the discretised equation becomes

$$\nfour \nfour \nfour \nfour \sum_f \, A_f \, \left[ \rho_f\,... ...left( \frac{\phi_P - \phi_A}{d_{AP}} \right) \right] \; = \; 0$$ (18.32)

The quantities F_f and D_f are now introduced where

$$F_f \; = \; A_f\,\rho_f (\underline{u}\,.\,\underline{n})_f \qquad \qquad D_f \; = \; A_f\, (\Gamma_\phi)_f / d_{AP}$$ (18.33)

F_f is the strength of the convection of $\phi$ and D_f is the diffusion conductance. Now () can be rewritten as

$$a_P \phi_P \; = \; \sum_{nb} \, a_{nb} \, \phi_{nb}$$ (18.34)

where the summation is over all neighbouring elements. The equations for the coefficients in equation () are

$$\begin{split} a_{nb}& = \; D_f - (1-\alpha_f) F_f\\ a_P \;... ...f\,F_f ) \\ & = \; \sum_{nb} a_{nb} + \sum_f F_f \end{split}$$ (18.35)

The discretisation techniques can be applied to the steady state continuity equation

$$\text{div}(\rho\,\underline{u}) \; = \; 0$$ (18.36)

which is a special case of the general conservation equation, with $\phi$ equal to 1.0 and $\Gamma_\phi$ equal to 0.0. Substitution of these values of $\Gamma_\phi$ and $\phi$ in equation () gives

$$\sum_f \; F_f \; = \; 0$$ (18.37)

Putting this equality into equation () gives

$$a_P \; = \; \sum_{nb} \, a_{nb}$$

Difficulties arise in obtaining the solution of this equation as the coefficients are not guaranteed to satisfy the condition that they remain positive. This can lead to the solution of the discretised equation becoming unbounded, $\phi_P > \phi_{nb}$ for all neighbours, which is physically unrealistic in the absence of a source (P. 83 Patankar for this and other objections).

In order to avoid the problems observed in the previous formulation a better assumption is sought for the variation of $\phi$ between adjacent nodes. In the upwind scheme, first suggested by Courant, Isaacson and Rees, there is no alteration to the handling of the diffusion term but the convection term uses the upwind nodal value of $\phi$ for the estimated interface value. Thus

$$\begin{split} \nfour \nfour \phi_f \; = \; \phi_P \qquad \q... ...; = \; \phi_A \qquad \qquad & \text{if } F_f < 0.0 \end{split}$$ (18.38)

In this case the coefficients in () are

$$\begin{split} a_{nb} \; & = \; D_f \, + \, \max(-F_f, \,0.0)... ...)] \\ & = \; \sum_{nb} \, a_{nb} + \sum_f \, F_f \end{split}$$ (18.39)

The upwind approach guarantees that all the coefficients remain positive and consequently leads to a bounded solution. Unfortunately the scheme over predicts diffusive effects leading to numerical smearing.

If the exact solution of the one dimensional convection - diffusion problem

$$\frac{d \,}{dx}(\rho u \phi) \; = \; \frac{d \,}{dx} \left( \Gamma_\phi \frac{d\phi}{dx} \right)$$ (18.40)

is sought over a region ranging from 0 to L with the boundary conditions

$$\begin{split} \text{At} \quad x \, = \, 0 \qquad & \phi \, =... ...At} \quad x \, = \, L \qquad & \phi \, = \, \phi_L \end{split}$$ (18.41)

the solution is

$$\frac{\phi \; - \; \phi_0}{\phi_L \; - \; \phi_0} \; = \; \frac{\text{exp}(Px/L) \, - \, 1}{\text{exp}(P) \, - \, 1}$$ (18.42)

where P is the Peclet number

$$P \equiv \frac{ \rho u L}{\Gamma_\phi}$$ (18.43)

which is the ratio of the strength of the convection to the strength of the diffusion. It can be seen that for approximately zero values of the Peclet number the problem reverts to one of pure diffusion. In this case the solution given by () means that the variation of $\phi$ between element centroids is nearly linear. When the flow is in the positive direction, u and hence P being positive, the values of $\phi$ in the domain are influenced more by $\phi_0$, which is the upstream value of $\phi$. Similarly for negative flows the value of $\phi$ is more dependent on its value at L which again is the upstream value. As P increases in either the positive or a negative manner the value of $\phi$ within the domain becomes more strongly influenced by the relevant upstream value of $\phi$. If |P| increases then d$\phi$/dx tends to 0 at the centre of the region indicating that diffusion is nearly absent. The upwind scheme always uses a linear relationship between $\phi$ and x and hence at large |P| values it overestimates the influence of the diffusion term.

Let $\underline{J}$ represent the total flux which is equal to the sum of the diffusion flux and the convective flux. Then

$$\underline{J} \; = \; \rho \underline{u} \phi \; - \; \Gamma_\phi \text{grad}(\phi)$$ (18.44)

and so convection - diffusion equation can be written as

$$\text{div}(\underline{J}) \; = \; 0$$ (18.45)

which after integration over a control volume transforms into the discretised equation

$$\sum_f \, A_f \, (\underline{J}\,.\, \underline{n})_f\; = \; 0$$ (18.46)

In the exponential scheme the exact solution of the equation as shown in equation () is used as a profile between points P and A. The suffices 0 and L are replaced by P and A and d_AP replaces L. Under this assumption the expression for $(\underline{J}\,.\,\underline{n})_f$ is

$$(\underline{J}\,.\,\underline{n})_f \; = \; F_f \left( \phi_P... ...rac{\phi_P \, - \, \phi_A} {\text{exp}(P_f) \, - \, 1} \right)$$ (18.47)

where

$$P_f \; = \; \frac{F_f}{D_f} \; = \; \frac{\rho_f (\underline{u}\,.\, \underline{n})_f} {(\Gamma_\phi)_f}$$ (18.48)

Using equation () the coefficients in the discretised equation become

$$\begin{split} a_{nb} \; & = \; \frac{F_f}{\text{exp}(P_f)-1} \\ a_P \; & = \; \sum_{nb}\,a_{nb} + \sum_f \, F_f \end{split}$$ (18.49)

This scheme produces the exact solution for the steady state one dimensional problem given any value of the Peclet number. The problems with this differencing scheme are that exponentials are expensive in terms of computational time and for greater than one dimensions the scheme is not exact and so the computation time cannot be justified.

If the equations () are considered it can be seen that as the Peclet number increases to positive infinity then

$$a_{nb} / D_f \longrightarrow 0.0$$

and as the Peclet number goes to negative infinity

$$a_{nb} / D_f \longrightarrow -P_f$$

and when the Peclet number is close to zero

$$a_{nb} / D_f = 1.0 \, - \, 0.5 P_f$$

The hybrid scheme, developed by Spalding, uses three straight line sections to approximate these properties. For P_f < -2 it uses

$$a_{nb} / D_f \; = \; -P_f$$

for $\vert P_f\vert \; <$ 2

$$a_{nb}/D_f \; = \; 1.0 \, - \, 0.5 P_f$$

and for P_f > 2

$$a_{nb} / D_f \; = \; 0.0$$

This means that

$$\begin{split} a_{nb} \; & = \; \max \left(-F_f,D_f - \frac{F... ... a_P \; & = \; \sum_{nb} \, a_{nb} + \sum_f \, F_f \end{split}$$ (18.50)

It should be noted that when the Peclet number is in the range -2 to 2 then the hybrid scheme reduces to the central difference scheme and outside this range it uses a modification of the upwind scheme where the diffusion has been set to zero. In this way the shortcomings of the upwind scheme are not shared by the hybrid scheme.

The hybrid scheme provides a good approximation to the exact solution but there is a rather large difference around |P_f| = 2 due to the fact that the diffusion is set to zero at these values. The power law scheme, described by Patankar, provides a better approximation by means of a slightly more complex algorithm which is not that expensive to compute. The power law uses

It can be seen that for |P_f| > 10 the power law and the hybrid schemes are identical. The power law scheme gives an extremely close approximation to the exponential scheme and because of its relatively low computational expense it offers an attractive alternative.

The discretised convection - diffusion equation can now be written in a form independent of the differencing scheme being used

$$\nfour \nfour \sum_f \, \left[ D_f\, A(\vert P_f\vert) + \ma... ...ght] \left( \phi_P - \phi_A \right) + F_f \, \phi_P \; = \; 0$$ (18.51)

where the formulae used for the function A(|P|) for the various first order differencing schemes are given in Table . A considerable amount of effort has been expended by many researchers to build and improve on these schemes.

 

Table: Definition of Differencing Schemes

|t:==:t| Scheme	Formulae for A(|P|)
||-|| Central Differencing $\quad$	1 - 0.5 |P|
||-|| Upwind	1
||-|| Hybrid	max ( 0, 1-0.5 |P| )
||-|| Power Law	max(0, (1 - 0.1|P| )5 ) $\quad$
||-|| Exponential	|P| / exp( |P| ) - 1
|b:==:b|