Comparison of three finite-volume interpolation schemes

 

 

T.M. Harms^I; T.W. Von Backström^II; J.P. du Plessis^III

^ISenior Lecturer, MSAIMechE
^IIProfessor, Department of Mechanical Engineering, University of Stellenbosch, Private Bag XI, Matieland, 7602 South Africa
^IIIProfessor, Department of Applied Mathematics, University of Stellenbosch, Private Bag XI, Matieland, 7602 South Africa

 

 

---

ABSTRACT

Some fundamental aspects of the finite-volume numerical flow analysis approach, which are generally assumed to be understood but at times still incorrectly interpreted in the literature, are reiterated. In this context an accuracy comparison between linear, upwind quadratic and locally analytic interface interpolation has been carried out. The formulations adopted were chosen with stability and arbitrary non-orthogonal physical grids in mind. Artificial diffusion damping through flux blending has been provided for. A finite-volume quality norm is proposed as a qualitative measure to evaluate different discretization approaches. Three laminar flow problems are employed as test cases, namely, the concentric shear flow, wall driven flow in a square cavity, and the mixing of two scalar streams. The results indicate that in the light of growing computer power, simple linear second order interpolation in conjunction with the ability to introduce controlled artificial diffusion damping is adequately accurate and likely to dominate finite volume thinking in the current decade.

---

 

 

1 Introduction

The segregated SIMPLE numerical finite-volume flow analysis procedure, in conjunction with the pressure correction approach, continues to be extensively applied to every flow problem imaginable. In this context the definition of a fluid by Daugherty and Franzini¹ as anything which deforms easily under stress' is useful because it encompasses modelling applications to plastic deformation 'flows' for example. The phenomenal success of the method can partly be attributed to the clarity of the writing of its early contributors, e.g. Patankar and Spalding,² Patankar,³ and Gosman et al.⁴

During the past ten years the field of computational fluid dynamics has experienced an explosion-like surge of publications, some extensively mathematical and thus often rendering it difficult to appreciate the new insight into the relationship between the numerics and physics offered. The current paper presents an attempt to reiterate, in the context of a comparison between three interpolation schemes, some fundamental aspects of the SIMPLE approach, which in the experience of the authors and their students are still often misunderstood.

 

2 Physical model

For the purposes of this paper two-dimensional, incompressible, steady, constant viscosity laminar flow is considered, which is well described by the Navier-Stokes momentum equations and the continuity equation. In the general form described by Patankar³

u is the velocity vector as a measure of convection, ΓΦ the diffusion coefficient and SΦ the volumetric source of any conserved scalar property Φ. In the case of the momentum equations, Φ becomes the Cartesian velocity component U (or V) as a measure of momentum per unit mass, and SΦ the respective negative Cartesian pressure gradient sources. Since the Cartesian velocity components are transformation invariant in space,⁵ they can be regarded as scalar quantities, which is significant for the analysis in non-orthogonal geometries. The continuity equation will emerge from equations (1) and (2) by setting Φequal to unity and SΦ (and ΓΦ) equal to zero.

 

3 Finite volume methodology

A conservation law becomes useful through the fact that the property under consideration is conserved, but its distribution or concentration in space (and time) is not. An equation can therefore be formulated from the balance of the integrated convection-diffusion transport, J, of the property through the boundary of a control volume and the balancing integrated change of state, SΦ, inside this volume. This classical (computational) thermodynamic approach is directly implemented in the physical finite-volume method. The typical cell-centred nodal arrangement (as opposed to corner or interface nodal positions) shown in Figure 1 reflects the unique association between the control volume, the conservation equation and the location of the property value sought. To save space, a local ξ - n interface coordinate system is also indicated in Figure 1. This is further referred to below.

 

 

The partial derivatives contained in equations (1) and (2) can either be approximated through finite differences or determined exactly from a shape function, which approximates the local property field based on, for example, the nine nodal values shown in Figure 1. In the current work, finite differences have been employed, except where otherwise indicated in the following sections.

Specific to the cell-centred formulation is the requirement to determine the interface convection-diffusion flux, which Schneider and Raw⁶ referred to as the first closure problem of the method. Spalding,⁷ for example, implicitly showed that a stabilizing feature of fluid flows is the fact that, due to the effects of convection, the downstream events in the flow become less influential than upstream events when the local Péclet number exceeds approximately 2. Linear, quadratic or locally analytic interpolation is therefore practical, provided numerical measures are taken to simulate this stability efficiently. Khosla and Rubin⁸ showed that, analogous to Spalding's⁷ idea, this can be substantially and conveniently achieved by transferring mainly downstream influences in the interface fluxes, J, from the left hand divergence side of the discrete integrated form of equation (1), from which coefficients in the matrix equations to be solved are assembled, to the right hand source term side. This idea, previously used for example by Denton⁹ and later by Peric' ,¹⁰ remained surprisingly neglected during the 1980s.

In the current work the central difference or linear interpolation formulation adopted is given by

where m is the discrete gradient defined by the two nodes adjacent to the interface. Substituting this expression in equation (2), and employing a maximum function to heed the stability considerations previously mentioned, results in

The flux blending factor A described by Peric ¹¹ identifies those terms, which, upon integration, are to be transferred to the source term side of equation (1), where they then might be referred to as secondary sources. This factor is normally of unit value, but can be reduced to enhance stability by effectively weakening the downstream influences. When this is necessary, it is an indication that the numerical and/or physical model represent the physics inadequately. This is especially the case in high Reynolds number or turbulent flows. Pulliam¹² showed how the reduction of A is equivalent to generating false or artificial diffusion damping, with of course an associated reduction in accuracy.

Employing the quadratic upstream interpolation for convection kinematics (quick) referred to by Leonard and Drummond,¹³ the coefficients of an interpolation parabola can be determined from the two nodal values adjacent to the interface and the upstream nodal property gradient in the interpolation direction. Such a formulation is suitable for flow analysis based on arbitrary non-orthogonal physical grids, since three collinear nodes are not required. The two resultant equations (see Appendix 1) can be combined to yield

and

For the exceptional case of Pe - 0, equations (3) and (4) replace equations (5) and (6).

Finally, equation (1) can be solved to yield the locally analytic interpolation function (loads) of Wong and Raithby,¹⁴ if multi-dimensional terms on the left hand side are transferred to the right and the resultant source terms on the right hand side are treated as piecewise constant. This idea was previously tested by Prakash,¹⁵ Huang et al.¹⁶ and the present authors¹⁷ who adopted the following form:

Typically, for the U-component, the source Sn would consist of

The exponential weighting functions A, B, C, and D (see Appendix 2) were here approximated as previously described by the authors,¹⁷ although more efficient expressions are now indicated.¹⁸

Since secondary source terms arise in all three schemes as the result of the deferred correction approach, it is suggested that these schemes are computationally competitive among themselves. This assumes that the various interpolation expressions are optimized with regard to computational expense and comparative associated accuracy. The formulations presented have however been chosen such that the coefficients of the discretized conservation equations approximate those resulting from the locally analytic solution of the homogeneous one-dimensional form of equation (1), i.e. the exponential coefficients described by Patankar,³ with the view to maximizing stability. All three schemes are of second or higher order accuracy, depending on discretization options and flow conditions.

 

4 Single grid

When equations for the velocity components U and V based on equation (1) are differentiated with respect to X and y in turn and the result summed, the following Poisson equation for pressure results

This situation is indicated in Figure 2. For an equal order discretization of the Laplacian, it is necessary to obtain a local solution of the momentum equations in the space between nodes. In other words, deriving equation (10) through the classical route³ of enforcing mass conservation over the control volume interfaces, a function is needed which expresses the interface velocity in terms of the interface pressure gradient. The first of three routes to obtain this is the solution of the momentum equations on staggered grids, as described by Patankar.³ Lee and Chiu¹⁹ recently summarized the disadvantages of this initially mathematically attractive route, which must now be rejected as unergonomic, because it complicates the mental task of the analyst. Difficulties generally arise in three-dimensional non-orthogonal physical geometries, especially where the grid structure is complicated.

 

 

An alternative route therefore to obtain a velocity-pressure coupling at the control volume interface is provided by equation (7). Since provision is made for the local pressure gradient in the source term, equation (7) can simply be differentiated with regard to the local pressure gradient to obtain an expression of velocity correction to compile a pressure correction equation.^3,17

A third method was proposed by Rhie and Chow,²⁰ who simulated a staggered grid by linearly interpolating the interface velocities in terms of the neighbouring discretized momentum equations, except for the pressure gradient contained therein. In this respect the approach is identical to the unequal order method described by Patankar and Baliga,²¹ who solved for pressure at alternate grid points only. In the method of Rhie and Chow²⁰ the pressures at four nodes, which define the pressure gradients at the two adjacent control volumes in the direction of the interpolation, are assumed to lie on a parabola, i.e. a linear pressure gradient field is implied. This reduction of the degree of freedom, i.e. from an at least cubic to a quadratic pressure profile, provides the pressure field smoothing necessary to avoid checker board pressure fields.³ It is however also the reason why, when grids are too coarse, the computed interface velocities can produce a checkerboard pattern with regard to the nodal velocities.

 

5 Quality norm

Equation (1) suggests a way to quantify the performance or quality of a finite-volume analysis. If the divergence over a control volume is large, or, alternatively the source term is large, an indication is provided that local changes in the flow field are severe. It is then useful to plot contours of absolute integrated source terms to obtain an immediate visual indication of difficult flow regions, where for example grid refinement might be called for.

Such difficulties would hot only arise due to say severe physical property gradients, but also due to the numerics. It was discussed above that for numerical stability the upwinding divergence remains on the left of equation (1). As a result, if the net effect of the secondary terms transferred to the right becomes relatively significant, this will also be indicated on the contour plots. The same would apply to non-orthogonal interface flux contributions, which are treated explicitly.

Summing the absolute integrated source terms over the entire flow domain may serve as an attempt to quantify the quality of the solution by a single norm. The con-ververgence of such a number will indicate to what degree a grid independent solution has been approached by comparison of the value obtained on successive smaller grids. The number should also add useful information for the comparison of different discretization options. Furthermore, the information is a priori in the sense that it does not depend on a known solution. If the value is normalized by a characteristic global (g) flux of the problem concerned, the resultant number, Q, may give a comparative indication for different flow problems of the relative significance of 'source term activity' therein:

 

6 Implementation

Instead of the classical³ geographic indentification of fluxes over different control volume faces, nodal neighbours (and the associated interfaces) are numbered sequentially and only one interface flux calculation is programmed as part of a loop. This calculation is based on a local coordinate system, which has been indicated in Figure 1 for interface number 2. This approach provides maximum freedom with regard to the shape of the control volume in more complex geometries.

In the current work boundary conditions are treated explicitly, i.e. Neumann boundary conditions are converted to Dirichlet conditions through extrapolation. As a result boundary interface fluxes are treated exactly as internal fluxes, since equations (4), (6), and (8) will recover the boundary values automatically, when whole boundary control volumes are employed, i.e. a Patankar³ type B grid. Steady flow under-relaxation as described by Patankar³ has been applied and the tri-diagonal, line-by-line algebraic equation solver has been employed.

 

7 Test cases

7.1 Concentric shear flow

The configuration of the concentric shear flow problem is shown in Figure 3. This frequently used test case previously employed by the authors¹⁷ has an analytic solution given by e.g. Schlichting.²² It may represent the first configuration towards modelling a mixing vessel for example. A regular square grid implies values of the interpolation ratio a of 0.5 and 1, i.e. the interface lies halfway between nodes or at the location of a global boundary node. Analytical velocities were implemented as Dirichlet boundary conditions and the analytical pressure difference was used as a Neumann condition. Following Prakash²³ the problem can be non-dimensionalized in terms of a Reynolds number given by p(wr₂) (r2 - r₁) /Γ. Laminar flow tests were run at Reynolds number 1000 using 9 x 9, 21 x 21 and 41 x 41 nodes. The three interface interpolation schemes, described by equations (4), (6), and (8), were employed. Relaxation factors for velocities and pressure were 0.5 and 0.2, respectively, and the net secondary (higher order) source contribution in the momentum equations was relaxed by 0.04 on the finest grid. The latter value is arbitrary and not optimized. Some additional relaxation seems necessary when the velocity/grid size ratio becomes large, i.e. the iteration 'time step' required some reduction. The interface velocity-pressure coupling of Rhie and Chow²⁰ was used in all cases. For central difference interpolation on the 9x9 grid some artificial viscosity proved necessary to obtain convergence. This was applied by setting A equal to 0.9. Starting with a constant radial outward flow field (approximately 50% of maximum analytical tangential velocity), the total number of outer iterations allowed varied between 195 on the smallest to 4 182 on the largest grids. In all cases the normalized summation of absolute residuals suggested by Peric' et al²⁴ for momentum and mass were less than 10^-6 when iterations were halted. The normalization is based on the above Reynolds number e.g. momentum related quantities were normalized by p(ωr₂) (ωr₂) (r₂ - r₁). In terms of computational costs such convergence is not normally required. It was however applied to enhance the reliability of the comparison.

 

 

In Figure 4 the average absolute percentage errors of the computed U-velocity components for the three grid sizes and interpolation methods are shown.

 

 

In terms of accuracy on the coarsest grid the higher order schemes perform well. The physically based loads is slightly better than the numerically based quick scheme. On the other hand, as the grid is refined, the second order central difference scheme rapidly approaches the accuracy of the other two schemes. From equations (4), (6), (8), and (9) it follows that the central difference scheme is the simplest to employ. In the context of very large problems, as occur in climatology or oceanography for example, and also hardware limitations, e.g. array sizes, the use of more accurate schemes might on the other hand still be worth while for some analyses.

In Table 1 the final values of the U-momentum quality norm obtained, i.e. as defined by equation,¹² are presented. The values are essentially the same, which is to be expected. All three schemes use similar up winding coefficients, which follow from equations (4), (6), and (8), and the total integrated source should be independent of grid size. Since this is a relatively regular flow field, these figures may serve as a comparison with data obtained in subsequent test cases.

 

 

For this well-behaved tlow, contours of absolute integrated source terms do not highlight any significant features. However in Figure 5 absolute U-velocity component source term contours for the loads scheme, in this case locally normalized, can be compared to the contours of the locally normalized percentage error plot.

 

 

The maximum error is then found in the vicinity of the fixed inner shaft, where the velocities used for normalization tend to zero. Some correlation between the two figures can be detected.

7.2 Wall driven flow in a square cavity

The steady incompressible lid driven flow in a square cavity complements the previous case in as much it is not a through-flow problem, but relies on the presence of viscosity to establish a recirculating flow pattern, to be captured using a non-flow aligned square grid. The essential features of such a resulting flow pattern can for example be observed in a fume extraction hood over a conveyor belt.

Huang et al.¹⁶ did report failure to obtain convergence when applying the loads scheme of Wong and Raithby¹⁴ to the cavity problem. Prakash²³ showed that convergence could be obtained for his implementation of loads, with the help of near boundary grid refinement. This was confirmed by the current authors¹⁷ for the loads scheme. As a result the cavity problem might also provide a test of the quality norm approach proposed.

Thus the shorter side length of half size boundary control volumes of an otherwise regular grid is calculated from H/(n - 3)/2. The values of the interface interpolation ratio a occurring are therefore 0.33, 0.5, 0.67 and 1.0. The Reynolds number in terms of lid speed and side wall length H tested was 400 and 11 xll, 21 x 21 and 41 x 41 node grids were employed. Relaxation factors were the same as for the previous case. Apart from the moving top side velocity of unity, all other boundary nodal velocities were set to zero. Using a unity side wall length results in a normalizing factor of unity. Iterations were halted when normalized summation of absolute residuals of momentum and mass were less than 3-10^-6. This was achieved after allowing between 500 and 6 000 iterations respectively for the grids used. Again the interface velocity-pressure coupling of Rhie and Chow²⁰ was used.

In Figure 6 the computed U-component on the vertical centreline obtained for the various schemes are compared to the results obtained by Winters and Cliffe.²⁴ Their grid consisted of 57 x 57 finite element nodes with additional local grid refinement in the lid side corners. On the coarsest grid the resolution of the flow field must be considered inadequate. The accuracy improves substantially with grid refinement, the central difference scheme proving once again adequate. The same weak hierarchy of accuracy among the three schemes as indicated in the previous test case can be detected.

In Table 2 U-momentum quality norms obtained according to equation (12) are presented. Although the data is too limited to make conclusive statements, the following tentative observations can be added to the experience gained in the previous test case. With regard to the first two schemes a grid independent solution has not yet been obtained. All numbers based on loads tend towards a lower value than the those based on the other schemes. All the numbers appear to converge to a significantly larger value than those obtained in the previous test case, which might be attributed to the effect of normalization. The accuracy results of both test cases seem to suggest that a lower norm value can be associated with a better quality resolution of the flow field. The particular difficulty of the problem can be identified from a typical absolute integrated U-component source term contours of the central difference case shown in Figure 7 (contours for the other interpolation schemes being similar). The finest grid used appears necessary to resolve the severity of the flow field in the top two corners, where relatively large pressure gradients are required to change the direction of the flow.

 

 

 

 

7.3 Scalar diffusion test

Another frequently used test case to examine interpolation schemes is the transport of a scalar step discontinuity resulting from the joining of two uniforms flows transporting the same scalar at different magnitudes and at various directions relative to the grid. Specifying a 'sharp' inlet profile not only allows an examination of false diffusion 'smearing', but alsc the capacity of the method to deal with steep property gradients in the flow. Furthermore this test tends to give a clear indication of the extent of the inherent dispersion associated with higher order schemes.

The typical configuration is shown in Figure 8. The flow direction is described by the elevation (relative velocity components) of the line, at which the discontinuity occurs. In the current case only a 45° fixed (not computed) flow field is employed. The scalar step profile is imposed by the boundary conditions indicated in Figure 8. In order to examine the three schemes in the context of the general formulations adopted in equations (4), (6), and (8), the diffusion coefficient, e.g. in equation (1) was set at 10^-10 and the primary source term was set to zero. The general formulations are not suited for zero diffusion. The value of 10~¹⁰ was arrived at by noticing no further graphical changes in the scalar profiles examined for a coefficient smaller than 10~⁶. The regular grid sizes examined were 11 xll, 21 x 21 and 41 x 41. The initial scalar field value was set to 0.5. In all cases no relaxation was employed and less than a hundred iterations were required to reduce the normalized (by a factor one) scalar residual to below 10^-6. However, a non-optimized value of 0.9 was assigned to the value of A for the central difference interpolation to introduce some artificial diffusion. Without this remedy convergence was either extremely slow or produced 'sawtooth' results on the finest grid.

 

 

In Figure 9 a comparison between the three interpolation schemes is presented in the form of scalar profiles along the line y - 1 - x. In the context of yielding a vertical step face, the QUICK scheme seems to perform best. On the coarsest grid, all three schemes produce almost equal order unboundedness. With increasing grid refinement, the dispersive overshoots in the damped central difference based results rapidly disappear while they increase for the other two schemes. In terms of accuracy the central difference scheme performs therefore best, although, through the introduction of artificial diffusion damping, unboundedness was introduced.

 

 

In Table 3 the computed quality norms associated with the results presented in Figure 9 are provided. The norm will diminish, as the step is approximated more accurately and the regions of 'flatness', where the first order divergence vanishes, expand and the overshoots disappear. The rapid improvement in the damped central difference based results is clearly distinguishable.

 

 

In Figure 10 contours of the source term distribution are provided for the 41 x41 nodes central difference case, the other cases being similar. The severity of the boundary condition imposed on the scalar field at the inlet flow corner is quite recognizable in all cases.

 

 

8 Conclusion

Due to its simplicity, inherent boundedness and the transparent manner described by Peric',¹¹ in which artificial damping or downwind decoupling can be introduced, it is likely that central difference interpolation will be widely employed.

The alternatives are more accurate, but also more complicated and less transparent. In the light of rapidly rising computer power, the latter disadvantage is gaining in significance. The results of the test cases do not indicate that grid convergence is substantially enhanced by non-linear interpolation. However, physical circumstances, where the use of non-linear interpolation is essential, may occur.

The analysis of source term distribution is apt in finite volume methodology. The quality norm proposed is a heuristic quantitative tool. Its formulation is presented as a first attempt to generate a simple but useful index. It is hoped to report on alternatives and their implications in the future.

 

Nomenclature

a fraction describing the interface location ξ - aL

A, B,C, D exponential functions and their

approximations g global term characteristic for the flow problem

H typical domain dimension

J interface convection-diffusion flux vector

J flux component

L distance between nodes

m local property gradient

M neighbour number

n number of nodes in one direction

N node number

Pe local grid Péclet number puLj p

p pressure at a node

Q finite volume quality norm

r₁, r₂ inner and outer cylinder radii

Re Reynolds number

S local constant volumetric source

u local convection velocity

u normal interface velocity

U, V orthogonal velocity components at a node

Γ diffusion coefficient

Ω control volume size

n local interface coordinate

λ flux blending factor

μ dynamic fluid viscosity

ξ local interface coordinate

p fluid density

σ control volume surface

Φ any transported variable

ω absolute value of angular velocity

 

References

1. Daugherty RL & Franzini JB. Fluid mechanics with engineering applications. 7th edn. McGraw-Hill, New York, 1977.

2. Patankar SV & Spalding DB. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. Journal Heat Mass Transfer, 1972, 15, pp. 1787-1806.         [ Links ]

3. Patankar SV. Numerical heat transfer and fluid flow. 1st edn. Hemisphere, Washington, 1980.

4. Gosman AD, Launder BE & Reece GJ. Computer-aided engineering: heat transfer and fluid flow. Hor-wood, Chichester, West Sussex, 1985.

5. Aris R. Vectors, tensors, and the basic equations of fluid mechanics. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1962.

6. Schneider GE & Raw MJ. A new control-volume-based finite element procedure for the numerical solution of the fluid flow and scalar transport equations. Final report, Department of Mechanical Engineering, University of Waterloo, Waterloo, Ontario, N2L 3G1, 1985.

7. Spalding DB. A novel finite difference formulation for differential expressions involving both first and second derivatives. International Journal for Numerical Methods in Engineering, 1972, 4, pp.551-559.         [ Links ]

8. Khosla PK & Rubin SG. A diagonally dominant second-order accurate implicit scheme. Computers and Fluids, 1974, 2, pp.207-209.         [ Links ]

9. Denton JD. The calculation of fully three dimensional flow through any type of turbomachinery blade row. Lecture series No. 140, Advisory Group for Aerospace Research and Development, 1984.

10. Peric' M. Analysis of pressure-velocity coupling on nonorthogonal grids. Numerical Heat Transfer, 1990, Part B, 17, pp.63-82.         [ Links ]

11. Peric' M. A finite volume method for the prediction of three-dimensional fluid flow in complex ducts. PhD thesis, Mechanical Engineering Department, Imperial College, University of London, 1985.         [ Links ]

12. Pulliam TH. Artificial dissipation models for the Euier equation. AIAA Journal, 1986, 24, pp.1931-1940.         [ Links ]

13. Leonard BP &; Drummond JE. Why you should not use 'hybrid', 'power-law' or related exponential schemes for convective modelling - there are much better alternatives. International Journal for Numerical Methods in Fluids, 1995, 20, pp.421-442.

14. Wong H H & Raithby GD. Improved finite-difference methods based on a critical evaluation of the approximation errors. Numerical Heat Transfer, 1979, 2, pp.139-163.         [ Links ]

15. Prakash C. Application of the locally analytic differencing scheme to some test problems for the convection-diffusion equation. Numerical Heat Transfer, 1984, 7, pp.165-182.         [ Links ]

16. Huang PG, Launder BE k Leschziner MA. Discretization of nonlinear convection processes: A broad-range comparison of four schemes. Computer Methods in Applied Mechanics and Engineering, 1985, 48, pp.1-24.         [ Links ]

17. Harms TM, Von Backström TW k Du Plessis JP. Reformulation of the simplen discretization scheme to accommodate noncentralized interfaces. Numerical Heat Transfer, 1991, Part B, 20, pp.127-144.         [ Links ]

18. Harms TM, Von Backström TW k Du Plessis J P. Derivation of a modified hybrid approximation. AIAA Journal, 1994, 32, pp.1541-1543.         [ Links ]

19. Lee D k Chiu J J. Covariant velocity-based calculation procedure with nonstaggered grids for computation of pulsatile flows. Numerical Heat Transfer, 1992, Part B, 21, pp.269-286.         [ Links ]

20. Rhie CM k Chow WL. Numerical study of turbulent flow past an airfoil with trailing edge separation. AIAA Journal, 1983, 21, pp.1525-1532.         [ Links ]

21. Baliga BR k Patankar SV. A control volume finite-element method for two-dimensional fluid flow and heat transfer. Numerical Heat Transfer, 1983, 6, p.245-261.         [ Links ]

22. Schlichting H. Boundary-layer theory. 7th edn. McGraw-Hill, New York, 1979.

23. Prakash C. An improved control volume finite-element method for heat and mass transfer, and for fluid flow using equal-order velocity-pressure interpolation. Numerical Heat Transfer, 1986, 9, pp.253-276.         [ Links ]

24. Peric' M, Kessler R k Scheuerer G. Comparison of finite-volume numerical methods with staggered and collocated grids. Computers and Fluids, 1988, 16, pp.389-403.         [ Links ]

25. Winters KH k Cliffe KA. A finite element study of driven cavity flow in a square cavity. AERE-R9444, Harwell, 1979.

 

 

Received July 1996
Final version January 1997

 

 

Appendix 1

Quadratic upwind interpolation

With reference to Figure 1 consider the interpolation equation

Employing the boundary values Φ_N(ξ = 0) and ΦM (ξ- L) and either their nodal derivatives m_n or m_Min the ξ-direction, depending whether the interface velocity u is positive or negative, results in

These equations can be combined resulting in equation (5).

 

Appendix 2

Exponential weighting functions

A removable singularity appears in these equations when Pe - 0. This is accomodated through the numerical approximation employed.¹⁷