In theory numerical results may be obtained that are indistinguishable from the ‘exact’ solution of the transport equation when the number of computational cells is infinitely large, irrespective of the differencing method used. However, in practical calculations we can only use a finite – sometimes quite small – number of cells, and our numerical results will only be physically realistic when the discretisation scheme has certain fundamental properties.
- Finite size of control volume introduce numerical issues.
- In order to analyze numerical errors, the numerical discretization scheme are tested for the following three properties.
Integration of the convection–diffusion equation over a finite number of control volumes yields a set of discretised conservation equations involving fluxes of the transported property φ through control volume faces. To ensure conservation of φfor the whole solution domain the flux of φ leaving a control volume across a certain face must be equal to the flux of φentering the adjacent control volume through the same face. To achieve this the flux through a common face must be represented in a consistent manner – by one and the same expression – in adjacent control volumes.
- To ensure conservation of φ for the whole solution domain, the flux of φ leaving a control volume across a certain face must be equal to the flux of φ entering the adjacent control volume through the same face.
- For example, consider the one-dimensional steady state diffusion problem without source terms shown below:
The fluxes across the domain boundaries are denoted by qA and qB. Let us consider four control volumes and apply central differencing to calculate the diffusive flux across the cell faces. The expression for the flux leaving the element around node 2 across its west face is Γw2 (φ2 −φ1)/δx and the flux entering across its east face is Γe2(φ3−φ2)/δx. An overall flux balance may be obtained by summing the net flux through each control volume, taking into account the boundary fluxes for the control volumes around nodes 1 and 4:
Since Γe1=Γw2, Γe2=Γw3 and Γe3=Γw4 the fluxes across control volume faces are expressed in a consistent manner and cancel out in pairs when summed over the entire domain. Only the two boundary fluxes qA and qB remain in the overall balance, so above equation expresses overall conservation of property φ. Flux consistency ensures conservation of φ over the entire domain for the central difference formulation of the diffusion flux.
Inconsistent flux interpolation formulae give rise to unsuitable schemes that do not satisfy overall conservation. For example, let us consider the situation where a quadratic interpolation formula, based on values at 1, 2 and
3, is used for control volume 2, and a quadratic profile, based on values at points 2, 3 and 4, is used for control volume 3.
- Flux consistency ensures conservation of φ over the entire domain for the central difference formulation of the diffusion flux.
- Inconsistent flux interpolation schemes are pron to numerical errors such as quadratic interpolation formula.
As shown in Figure below, the resulting quadratic profiles can be quite different.
Consequently, the flux values calculated at the east face of control volume 2 and the west face of control volume 3 may be unequal if the gradients of the two curves are different at the cell face. If this is the case the two fluxes do
not cancel out when summed and overall conservation is not satisfied. The example should not suggest to the reader that quadratic interpolation is entirely bad. Further on we will meet a quadratic discretisation practice, the so-called QUICK scheme, that is consistent.
The discretised equations at each nodal point represent a set of algebraic equations that needs to be solved. Normally iterative numerical techniques are used to solve large equation sets. These methods start the solution
process from a guessed distribution of the variable φ and perform successive updates until a converged solution is obtained. Scarborough (1958) has shown that a sufficient condition for a convergent iterative method can be expressed in terms of the values of the coefficients of the discretised equations:
- For iterative solvers, the matrix must be diagonally dominant.
- All coefficients of the discretized equations should have the same sign (usually all positive).
Here a′P is the net coefficient of the central node P(i.e. aP−SP), and the summation in the numerator is taken over all the neighbouring nodes (nb). If the differencing scheme produces coefficients that satisfy the above criterion
the resulting matrix of coefficients is diagonally dominant. To achieve diagonal dominance we need large values of net coefficient (aP−SP) so the linearisation practice of source terms should ensure that SP is always negative. If this is the case −SP is always positive and adds to aP.
Diagonal dominance is a desirable feature for satisfying the ‘boundedness’ criterion. This states that in the absence of sources the internal nodal values of property φ should be bounded by its boundary values. Hence
in a steady state conduction problem without sources and with boundary temperatures of 500°C and 200°C, all interior values of T should be less than 500°C and greater than 200°C. Another essential requirement for boundedness is that all coefficients of the discretised equations should have the same sign(usually all positive). Physically this implies that an increase in the variable φ at one node should result in an increase in φ at neighbouring nodes. If the discretisation scheme does not satisfy the boundedness requirements it is possible that the solution does not converge at all, or, if it does, that it contains ‘wiggles’.
The transportiveness property of a fluid flow (Roache, 1976) can be illustrated by considering the effect at a point P due to two constant sources of φ at nearby points W and E on either side as shown in Figure below. We define the non-dimensional cell Peclet number as a measure of the relative strengths of convection and diffusion:
- Relative strength of convection and diffusion is measured by Peclet number.
The lines in Figure below indicate the general shape of contours of constant φ (say φ=1) due to both sources for different values of Pe. The value of φ at any point can be thought of as the sum of contributions due to the two
- Consider the effect at a point P due to two constant sources of φ at nearby points W and E on either side
Let us consider two extreme cases to identify the extent of the influence at node P due to the sources at Wand E:
• no convection and pure diffusion (Pe→0)
• no diffusion and pure convection (Pe→∞)
In the case of pure diffusion the fluid is stagnant (Pe→0) and the contours of constant φ will be concentric circles centred around W and E since the diffusion process tends to spread φ equally in all directions. Figure (a) shows that both φ=1 contours pass through P, indicating that conditions at this point are influenced by both sources at W and E. As Pe increases the contours change shape from circular to elliptical and are shifted in the direction of the flow as shown in Figure (b). Influencing becomes increasingly biased towards the upstream direction at large values of Pe, so, in the present case where the flow is in the positive x-direction, conditions at P will be mainly influenced by the upstream source at W. In the case of pure convection (Pe→∞) the elliptical contours are completely stretched out in the flow direction. All of property φ emanating from the sources at W and E is immediately transported downstream. Thus, conditions at P are now unaffected by the downstream source at E and completely dictated by the upstream source at W. Since there is no diffusion φP is equal to φW. If the flow is in the negative x-direction we would find that φP is equal to φE. It is very important that the relationship between the directionality of influencing and the flow direction and magnitude of the Peclet number, known as the transportiveness, is borne out in the discretisation scheme.
- An Introduction to Computational Fluid Dynamics, THE FINITE VOLUME METHOD, Second Edition
H K Versteeg and W Malalasekera.
- Computational Fluid Dynamics-I, Lecture-9 by Dr. Tariq Talha, College of EME, NUST, Pakistan.