Versteeg H., Malalasekra W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method

Подождите немного. Документ загружается.

5.3 THE CENTRAL DIFFERENCING SCHEME 137

Figure 5.2

Figure 5.3 The grid used for

discretisation

additional terms to account for convection. To solve a one-dimensional

convection–diffusion problem we write discretised equations of the form

(5.14) for all grid nodes. This yields a set of algebraic equations that is solved

to obtain the distribution of the transported property

. The process is now

illustrated by means of a worked example.

A property

is transported by means of convection and diffusion through

the one-dimensional domain sketched in Figure 5.2. The governing equation

is (5.3); the boundary conditions are

= 1 at x = 0 and

= 0 at x = L. Using

ﬁve equally spaced cells and the central differencing scheme for convection

and diffusion, calculate the distribution of

as a function of x for (i) Case 1:

u = 0.1 m/s, (ii) Case 2: u = 2.5 m/s, and compare the results with the

analytical solution

= (5.15)

(iii) Case 3: recalculate the solution for u = 2.5 m/s with 20 grid nodes and

compare the results with the analytical solution. The following data apply:

length L = 1.0 m,

= 1.0 kg/m

, Γ=0.1 kg/m.s.

exp(

ux/Γ) − 1

exp(

uL/Γ) − 1

−

Example 5.1

Solution

The method of solution is demonstrated using the simple grid shown in

Figure 5.3. The domain has been divided into ﬁve control volumes giving

x = 0.2 m. Note that F =

u, D =Γ/

x, F

= F

= F and D

= D

everywhere. The boundaries are denoted by A and B.

The discretisation equation (5.14) and its coefﬁcients apply at internal

nodal points 2, 3 and 4, but control volumes 1 and 5 need special treatment

since they are adjacent to the domain boundaries. We integrate governing

equation (5.3) and use central differencing for both the diffusion terms and

the convective ﬂux through the east face of cell 1. The value of

is given

at the west face of this cell (

= 1) so we do not need to make any

approximations in the convective ﬂux term at this boundary. This yields the

following equation for node 1:

(

) − F

= D

(

−

) − D

(

−

) (5.16)

For control volume 5, the

-value at the east face is known (

= 0). We

obtain

− (

) = D

(

−

) − D

(

−

) (5.17)

ANIN_C05.qxd 29/12/2006 04:36PM Page 137

138 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS

The matrix form of the equation set using

= 1 and

= 0 is

G 1.55 −0.45 0 0 0 J G

JG1.1J

H−0.55 1.0 −0.45 0 0 K H

KH0 K

H 0 −0.55 1.0 −0.45 0 K H

K = H0 K (5.19)

H 00−0.55 1.0 −0.45K H

KH0 K

I 000−0.55 1.45L I

LI0 L

The solution to the above system is

JG0.9421J

KH0.8006K

K = H0.6276K (5.20)

KH0.4163K

LI0.1579L

Table 5.1

Node a

= a

+ a

− S

1 0 0.45 1.1

−1.1 1.55

2 0.55 0.45 0 0 1.0

3 0.55 0.45 0 0 1.0

4 0.55 0.45 0 0 1.0

5 0.55 0 0.9

−0.9 1.45

Rearrangement of equations (5.16) and (5.17), noting that D

= 2Γ/

= 2D and F

F, gives discretised equations at boundary nodes of the

following form:

= a

+ a

+ S

(5.18)

with central coefﬁcient

= a

+ a

+ (F

− F

) − S

and

Node a

10 D − F /2 −(2D + F)(2D + F )

2, 3, 4 D + F /2 D − F /2 0 0

5 D + F /2 0 −(2D − F )(2D − F )

To introduce the boundary conditions we have suppressed the link to the

boundary side and entered the boundary ﬂux through the source terms.

(i) Case 1

u = 0.1 m/s: F =

u = 0.1, D =Γ/

x = 0.1/0.2 = 0.5 gives the coefﬁcients as

summarised in Table 5.1.

ANIN_C05.qxd 29/12/2006 04:36PM Page 138

5.3 THE CENTRAL DIFFERENCING SCHEME 139

Comparison with the analytical solution

Substitution of the data into equation (5.15) gives the exact solution of the

problem:

(x) =

The numerical and analytical solutions are compared in Table 5.2 and in

Figure 5.4. Given the coarseness of the grid the central differencing (CD)

scheme gives reasonable agreement with the analytical solution.

2.7183 − exp(x)

1.7183

Figure 5.4 Comparison of

numerical and analytical

solutions for Case 1

(ii) Case 2

u = 2.5 m/s: F =

u = 2.5, D =Γ/

x = 0.1/0.2 = 0.5 gives the coefﬁcients as

summarised in Table 5.3.

Table 5.3

Node a

= a

+ a

− S

10−0.75 3.5

−3.5 2.75

2 1.75 −0.75 0 0 1.0

3 1.75 −0.75 0 0 1.0

4 1.75 −0.75 0 0 1.0

5 1.75 0 −1.5

1.5 0.25

Table 5.2

Node Distance Finite volume Analytical Difference Percentage

solution solution error

1 0.1 0.9421 0.9387 −0.003 −0.36

2 0.3 0.8006 0.7963 −0.004 −0.53

3 0.5 0.6276 0.6224 −0.005 −0.83

4 0.7 0.4163 0.4100 −0.006 −1.53

5 0.9 0.1579 0.1505 −0.007 −4.91

ANIN_C05.qxd 29/12/2006 04:36PM Page 139

Comparison of numerical and analytical solution

The matrix equations are formed from the coefﬁcients in Table 5.3 by the

same method used in Case 1 and subsequently solved. The analytical solu-

tion for the data that apply here is

(x) = 1 +

The numerical and analytical solutions are compared in Table 5.4 and shown

in Figure 5.5. The central differencing scheme produces a solution that

appears to oscillate about the exact solution. These oscillations are often

called ‘wiggles’ in the literature; the agreement with the analytical solution is

clearly not very good.

1 − exp(25x)

7.20 × 10

140 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS

Figure 5.5 Comparison of

numerical and analytical

solutions for Case 2

Table 5.4

Node Distance Finite volume Analytical Difference Percentage

solution solution error

1 0.1 1.0356 1.0000 −0.035 −3.56

2 0.3 0.8694 0.9999 0.131 13.05

3 0.5 1.2573 0.9999 −0.257 −25.74

4 0.7 0.3521 0.9994 0.647 64.70

5 0.9 2.4644 0.9179 −1.546 −168.48

(iii) Case 3

u = 2.5 m/s: a grid of 20 nodes gives

x = 0.05, F =

u = 2.5, D =Γ/

= 0.1/0.05 = 2.0. The coefﬁcients are summarised in Table 5.5 and the

resulting solution is compared with the analytical solution in Figure 5.6.

Table 5.5

Node a

= a

+ a

− S

1 0 0.75 6.5

−6.5 7.25

2–19 3.25 0.75 0 0 4.00

20 3.25 0 1.5

−1.5 4.75

ANIN_C05.qxd 29/12/2006 04:36PM Page 140

5.4 PROPERTIES OF DISCRETISATION SCHEMES 141

The agreement between the numerical results and the analytical solution

is now good. Comparison of the data for this case with the one computed on

the ﬁve-point grid of Case 2 shows that grid reﬁnement has reduced the F/D

ratio from 5 to 1.25. The central differencing scheme seems to yield accurate

results when the F/D ratio is low. The inﬂuence of the F/D ratio and the

reasons for the appearance of ‘wiggles’ in central difference solutions when

this ratio is high will be discussed below.

The failure of central differencing in certain cases involving combined con-

vection and diffusion forces us to take a more in-depth look at the properties

of discretisation schemes. 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 inﬁnitely large, irrespective of the

differencing method used. However, in practical calculations we can only use

a ﬁnite – sometimes quite small – number of cells, and our numerical results

will only be physically realistic when the discretisation scheme has certain

fundamental properties. The most important ones are:

• Conservativeness

• Boundedness

• Transportiveness

5.4.1 Conservativeness

Integration of the convection–diffusion equation over a ﬁnite number of

control volumes yields a set of discretised conservation equations involving

ﬂuxes of the transported property

through control volume faces. To ensure

conservation of

for the whole solution domain the ﬂux of

leaving a

control volume across a certain face must be equal to the ﬂux of

entering

the adjacent control volume through the same face. To achieve this the ﬂux

through a common face must be represented in a consistent manner – by

one and the same expression – in adjacent control volumes.

For example, consider the one-dimensional steady state diffusion prob-

lem without source terms shown in Figure 5.7.

Figure 5.6 Comparison of

numerical and analytical

solutions for Case 3

Properties

of discretisation

schemes

5.4

ANIN_C05.qxd 29/12/2006 04:36PM Page 141

The ﬂuxes across the domain boundaries are denoted by q

and q

. Let us

consider four control volumes and apply central differencing to calculate the

diffusive ﬂux across the cell faces. The expression for the ﬂux leaving the

element around node 2 across its west face is Γ

(

−

x and the ﬂux

entering across its east face is Γ

(

−

x. An overall ﬂux balance may be

obtained by summing the net ﬂux through each control volume, taking into

account the boundary ﬂuxes for the control volumes around nodes 1 and 4:

− q

+Γ

−Γ

+Γ

−Γ

+ q

−Γ

= q

− q

(5.21)

Since Γ

=Γ

, Γ

=Γ

and Γ

=Γ

the ﬂuxes 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 ﬂuxes q

and q

remain in the

overall balance, so equation (5.21) expresses overall conservation of property

. Flux consistency ensures conservation of

over the entire domain for the

central difference formulation of the diffusion ﬂux.

Inconsistent ﬂux 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 proﬁle, based on values at

points 2, 3 and 4, is used for control volume 3.

As shown in Figure 5.8, the resulting quadratic proﬁles can be quite

different.

(

−

)

(

−

)

(

−

)

(

−

)

(

−

)

(

−

)

142 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS

Figure 5.7 Example of

consistent speciﬁcation of

diffusive ﬂuxes

Figure 5.8 Example of

inconsistent speciﬁcation of

diffusive ﬂuxes

Consequently, the ﬂux 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

ANIN_C05.qxd 29/12/2006 04:36PM Page 142

5.4 PROPERTIES OF DISCRETISATION SCHEMES 143

two curves are different at the cell face. If this is the case the two ﬂuxes do

not cancel out when summed and overall conservation is not satisﬁed. 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.

5.4.2 Boundedness

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 sufﬁcient condition for a convergent iterative method

can be expressed in terms of the values of the coefﬁcients of the discretised

equations:

! ≤ 1 at all nodes

< 1 at one node at least

(5.22)

Here a ′

is the net coefﬁcient of the central node P (i.e. a

− S

), and the

summation in the numerator is taken over all the neighbouring nodes (nb). If

the differencing scheme produces coefﬁcients that satisfy the above criterion

the resulting matrix of coefﬁcients is diagonally dominant. To achieve

diagonal dominance we need large values of net coefﬁcient (a

−

) so the

linearisation practice of source terms should ensure that S

is always

negative. If this is the case −S

is always positive and adds to a

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 bounded-

ness is that all coefﬁcients 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 neighbour-

ing 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’. This is powerfully illustrated by the results

of Case 2 of Example 5.1. In all other worked examples we have developed

discretised equations with positive coefﬁcients a

and a

, but in Case 2 most

of the east coefﬁcients were negative (see Table 5.3), and the solution con-

tained large under- and overshoots!

5.4.3 Transportiveness

The transportiveness property of a ﬂuid ﬂow (Roache, 1976) can be illus-

trated 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 5.9. We deﬁne

∑|a

|a′

ANIN_C05.qxd 29/12/2006 04:36PM Page 143

144 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS

Figure 5.9 Distribution of

in the vicinity of two sources

at different Peclet numbers:

(a) pure convection, Pe → 0;

(b) diffusion and convection

the non-dimensional cell Peclet number as a measure of the relative strengths

of convection and diffusion:

Pe == (5.23)

where

x = characteristic length (cell width)

The lines in Figure 5.9 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

sources.

Γ/

Let us consider two extreme cases to identify the extent of the inﬂuence

at node P due to the sources at W and E:

• no convection and pure diffusion (Pe → 0)

• no diffusion and pure convection (Pe →∞)

In the case of pure diffusion the ﬂuid 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 5.9a

shows that both

= 1 contours pass through P, indicating that conditions at

this point are inﬂuenced by both sources at W and E. As Pe increases the

contours change shape from circular to elliptical and are shifted in the direc-

tion of the ﬂow as shown in Figure 5.9b. Inﬂuencing becomes increasingly

biased towards the upstream direction at large values of Pe, so, in the present

case where the ﬂow is in the positive x-direction, conditions at P will be

mainly inﬂuenced by the upstream source at W. In the case of pure convec-

tion (Pe →∞) the elliptical contours are completely stretched out in the ﬂow

direction. All of property

emanating from the sources at W and E is imme-

diately 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

is equal to

. If the ﬂow is in the

negative x-direction we would ﬁnd that

is equal to

. It is very important

that the relationship between the directionality of inﬂuencing and the ﬂow

direction and magnitude of the Peclet number, known as the transportive-

ness, is borne out in the discretisation scheme.

ANIN_C05.qxd 29/12/2006 04:36PM Page 144

5.5 ASSESSMENT OF THE CENTRAL DIFFERENCING SCHEME 145

Assessment of

the central

differencing scheme

for convection---

diffusion problems

Conservativeness: The central differencing scheme uses consistent expres-

sions to evaluate convective and diffusive ﬂuxes at the control volume faces.

The discussions in section 5.4.1 show that the scheme is conservative.

Boundedness:

(i) The internal coefﬁcients of discretised scalar transport equation (5.14) are

+ D

− a

+ a

+ (F

− F

)

A steady one-dimensional ﬂow ﬁeld is also governed by the

discretised continuity equation (5.10). This equation states that

− F

) is zero when the ﬂow ﬁeld satisﬁes continuity. Thus the

expression for a

in (5.14) becomes equal to a

= a

+ a

. The

coefﬁcients of the central differencing scheme satisfy the Scarborough

criterion (5.22).

(ii) With a

= D

− F

/2 the convective contribution to the east coefﬁcient

is negative; if the convection dominates it is possible for a

to be

negative. Given that F

> 0 and F

> 0 (i.e. the ﬂow is unidirectional),

for a

to be positive D

and F

must satisfy the following condition:

= Pe

< 2 (5.24)

If Pe

is greater than 2 the east coefﬁcient will be negative. This

violates one of the requirements for boundedness and may lead to

physically impossible solutions.

In the example of section 5.3 we took Pe = 5 in Case 2 so condition (5.24) is

violated. The consequences were evident in the results, which showed large

‘undershoots’ and ‘overshoots’. Taking Pe less than 2 in Cases 1 and 3 gave

bounded answers close to the analytical solution.

Transportiveness: The central differencing scheme introduces inﬂuencing

at node P from the directions of all its neighbours to calculate the convective

and diffusive ﬂux. Thus the scheme does not recognise the direction of the

ﬂow or the strength of convection relative to diffusion. It does not possess

the transportiveness property at high Pe.

Accuracy: The Taylor series truncation error of the central differencing

scheme is second-order (see Appendix A for further details). The require-

ment for positive coefﬁcients in the central differencing scheme as given

by formula (5.24) implies that the scheme will be stable and accurate only if

Pe = F/D < 2. It is important to note that the cell Peclet number, as deﬁned

by (5.23), is a combination of ﬂuid properties (

and Γ), a ﬂow property (u)

and a property of the computational grid (

x). So for given values of

and

Γ it is only possible to satisfy condition (5.24) if the velocity is small, hence

in diffusion-dominated low Reynolds number ﬂows, or if the grid spacing is

small. Owing to this limitation central differencing is not a suitable discretisa-

tion practice for general-purpose ﬂow calculations. This creates the need for

discretisation schemes which possess more favourable properties. Below we

discuss the upwind, hybrid, power-law, QUICK and TVD schemes.

5.5

ANIN_C05.qxd 29/12/2006 04:36PM Page 145

One of the major inadequacies of the central differencing scheme is its inabil-

ity to identify ﬂow direction. The value of property

at a west cell face is

always inﬂuenced by both

and

in central differencing. In a strongly

convective ﬂow from west to east, the above treatment is unsuitable because

the west cell face should receive much stronger inﬂuencing from node W

than from node P. The upwind differencing or ‘donor cell’ differencing

scheme takes into account the ﬂow direction when determining the value at

a cell face: the convected value of

at a cell face is taken to be equal to the

value at the upstream node. In Figure 5.10 we show the nodal values used to

calculate cell face values when the ﬂow is in the positive direction (west

to east) and in Figure 5.11 those for the negative direction.

146 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS

Figure 5.10

Figure 5.11

The upwind

differencing

scheme

5.6

When the ﬂow is in the positive direction, u

> 0, u

> 0 (F

> 0, F

> 0),

the upwind scheme sets

and

(5.25)

and the discretised equation (5.9) becomes

− F

= D

(

−

) − D

(

−

) (5.26)

which can be rearranged as

+ D

+ F

)

= (D

+ F

)

+ D

to give

[(D

+ F

) + D

+ (F

− F

)]

= (D

+ F

)

+ D

(5.27)

When the ﬂow is in the negative direction, u

< 0, u

< 0 (F

< 0, F

< 0), the

scheme takes

and

(5.28)

Now the discretised equation is

− F

= D

(

−

) − D

(

−

) (5.29)

ANIN_C05.qxd 29/12/2006 04:36PM Page 146