## W

FIGURE 6.3 Distribution of <p around its source. (from Versteeg and Malalasekara, 1995. Printed with permission from the Publishers.)

at node P will be influenced not only by upstream conditions at W but also by all the conditions further downstream (node E). As the value of Peclet number increases (more convection), directionality of influence becomes increasingly biased towards the upstream direction. This means that conditions at node E are strongly influenced by those at P but conditions at P will experience only a weak influence from those at node E. At the extreme case of infinite Peclet number, the constant <p contours are completely stretched in the direction of flow and conditions at node E will not influence those at node P. Discretization schemes must respect the transportiveness property (directionality of influence) of flow processes.

Source term linearization and interpolation practices to estimate cell face values are discussed with reference to these desirable properties of the discretization method.

(a) Source term linearization

The linearization of source term, S^P should be a good representation of the S and <p relationship. Depending on the functional form of S^P, there are several different ways of formulating coefficients S^ and Sp of Eq. (6.7). Patankar (1980) discussed various linearization practices and recommended the following method:

where * indicates the guess value or the previous iteration value. This linearization practice is recommended provided that the source term decreases with increasing p. Thus, the coefficients of linearized source terms can be written:

This formulation ensures that the slope of the linearized source term representation (Eq. (6.7)) is the same as the slope of the non-linear source term at node P, and is always negative. It must be noted that linearization with steeper slopes normally leads to slower convergence. For detailed discussion of the effect of source term linearization on convergence, and on the handling of source terms with non-negative slope with p, see Patankar (1980).

(b) Interpolation practices

Some commonly used interpolation practices are discussed here. A simple and straightforward approximation to the value at the CV face center is linear interpolation between the two nearest nodes. At location e on a Cartesian grid (Fig. 6.2), general interpolation coefficients (ft) for such a scheme can be written:

where the linear interpolation factor, Xe is defined as

Interpolation with Eq. (6.25) is second-order accurate. This approximation is called the 'central differencing scheme' (CDS).

It can be seen that the central differencing scheme discussed above is conservative. The coefficients of CDS satisfy the Scarborough criterion. However, for uniform grid (Xe = 0.5), when the Peclet number is higher than 2, the coefficients aE will become negative [Fe > (De/ft)]. This violates the boundedness requirements and may lead to physically unrealistic solutions. At all values of Peclet number, CDS retains the same directionality of influence and, therefore, does not possess the transportiveness property. To gain the transportiveness property, several differencing schemes have been proposed. The simplest is the first-order upwind differencing scheme (UDS), which approximates the value of pe by retaining only the upstream influence (with reference to Eq. (6.10)):

It can be seen that this UDS is conservative and transportiveness is built into the formulation. It is the only discretization scheme which unconditionally satisfies the boundedness criteria and therefore, never leads to 'wiggles'. Unfortunately, the scheme is only first-order accurate and is numerically diffusive. This numerical diffusion is magnified in multidimensional problems if the flow is oblique to the grid. The rapid variations in the variables will be smeared out and since the rate of error reduction is only first order, very fine grids are required to obtain an accurate solution. In order to exploit the higher accuracy of CDS, several combinations of CDS and first-order UDS have been proposed (for example, hybrid differencing schemes and power law differencing schemes). In all these combined schemes, CDS is used as long as Peclet numbers are less than 2. For larger Peclet numbers, either the standard UDS or the power law variant is used (Patankar, 1980). These combined schemes are also unconditionally bounded and therefore widely used in various computational codes. However, the accuracy in terms of Taylor series truncation error is still only of first order.

Several attempts have been made to employ higher order interpolation schemes. One of the most popular schemes is QUICK (quadratic upstream interpolation for convective kinetics), proposed by Leonard (1979). In this scheme, the face value of p is obtained by a quadratic function passing through two bracketing nodes (on each side of the face) and a node on the upstream side (Fig. 6.4). The formulae for

FIGURE 6.4 Quadratic profiles used in QUICK scheme.

FIGURE 6.4 Quadratic profiles used in QUICK scheme.

estimating (e using Eq. (6.10) can be written:

For the uniform grid, the coefficients of the three nodal values involved in the interpolation become 3/8 for the downstream point, 6/8 for the first upstream node and -1/8 for the second upstream node. This scheme is more complex than CDS and it extends the computational molecule by one more node in each direction (the conventional tri-diagonal methods are, therefore, not directly applicable. See the discussion in the following subsection). The scheme has a third-order truncation error and was made popular by Leonard (1979). The transportiveness property is built into the scheme by considering two upstream and one downstream node. However, the main coefficients of the discretized equations are not guaranteed to be positive. This may lead to instability and may lead to unbounded (wiggles) solutions under certain conditions.

Several attempts have been made to reformulate the QUICK and other higher order schemes to ensure the boundedness property. One way to achieve this is to place the troublesome negative coefficients in the source term to alleviate the stability problems (see for example, Hayase et al., 1992 and references cited therein). Other attempts include modifications of higher order schemes either based on flux limiting or based on slope limiting (see Leveque, 1996 for a recent review). It will be instructive to examine the various discretization schemes using the so-called 'normalized variable diagram' to understand and to compare these modifications. Depending on the direction of normal velocity at the face, the locally normalized variable is defined as

where subscript D, U and C denote downstream, upstream and central node, respectively. It can be seen that normalized values of <p at downstream and upstream nodes are 1 and 0, respectively (Fig. 6.5). Because of this, all the discretization schemes can now be written in terms of the normalized variable at the center node. Most of the widely used discretization schemes are shown in Fig. 6.6 in terms of the normalized variable diagram (NVD). Leonard (1988) has shown that any linear NVD characteristic which passes through the second quadrant may produce unphysical oscillations (CDS, QUICK). Characteristics, which pass through the fourth quadrant (below point O), are artificially diffusive. Characteristics, which pass above point

FIGURE 6.5 Node variables in the vicinity of a CV Face (at broken line )

FIGURE 6.5 Node variables in the vicinity of a CV Face (at broken line )

P (1, 1), are oscillatory in two dimensions (second-order upwinding; for details of second-order upwinding schemes, see Shyy et al, 1992). Characteristics passing through point Q (0.5, 0.75) have second-order accuracy (third order, if the slope at Q is 3/4). Thus, NVD can be used to evaluate different discretization schemes as well as devise new ones.

Van Leer (1977) introduced a monotized centered (MC) scheme. The SUPER-BEE scheme is specifically developed to handle discontinuities (Roe, 1985). Gaskell and Lau (1988) proposed a SMART (sharp and monotonic algorithm for realistic transport) scheme, and Leonard (1988) proposed the SHARP (simple high accuracy resolution program) scheme based on modifications of QUICK to preserve bound-edness. These schemes are also shown in Fig. 6.6. All of these schemes try to avoid non-physical oscillations by introducing modifications around the basic QUICK scheme. The SUPERBEE limiter tends to be overcompressive, meaning that it tends to steepen up smooth profiles into discontinuities. For this reason, it is useful for preserving discontinuities but inappropriate for problems with smooth solutions. The desire to use higher order discretization schemes must be balanced against the limitations imposed by the complexity of the problem being solved, the availability of computing resources, the stability and convergence of the solution algorithm and the capacity to tolerate non-physical over- and undershoots. Among the various schemes, QUICK (and its modifications SHARP or SMART) and second-order upwind differencing schemes look more attractive from the point of view of accuracy and ease of implementation. For complex multiphase flows, the flow field may exhibit extreme sensitivity to the gradients of the dispersed phase volume fraction. In such cases, a hybrid or power law scheme advocated by Patankar (1980) may be a better choice, at least while initiating the simulations from an arbitrarily guessed flow field.

Using one of the suitable discretization schemes discussed above, it is possible to relate values of variables and their gradient at CV faces to the node values. It is also necessary to use suitable interpolation schemes to estimate other relevant quantities like effective diffusion coefficients, (r) at required locations. Either algebraic mean or harmonic mean can be used to estimate the value of effective diffusion coefficients at cell faces. For example, the effective diffusion coefficient at face e can be written (for a uniform grid):

## Post a comment