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

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

APPENDIX A 447

A third formula for (

∂φ

∂

can be obtained by rearranging equation (A.9) as

=+O(∆x

) (A.10)

Equation (A.10) uses values at E and W to evaluate the gradient at the mid-

point P, and is called a central difference formula. The central differencing

formula is second-order accurate. The quadratic dependence of the error on

grid spacing means that after grid reﬁnement the error reduces more quickly

in a second-order accurate differencing scheme than in a ﬁrst-order accurate

scheme.

In the ﬁnite volume discretisation procedure developed in section 4.2 the

gradient at a cell face, e.g. at ‘e’, was evaluated using

== (A.11)

By comparing formulae (A.10) and (A.11) it can be easily recognised that

(A.11) evaluates the gradient at the midpoint between P and E through a

central difference formula at point ‘e’: therefore its accuracy is second-order

for uniform grids.

It is relatively straightforward to demonstrate the third-order accuracy of

the QUICK differencing scheme for the convective ﬂux at a midpoint cell

face in a uniform grid of spacing ∆x. The QUICK scheme calculates the

value

at the east cell face of a general node as

−

(A.12)

Taylor series expansion about the east face gives

+∆x

+ O(∆x

) (A.13)

−∆x

+−∆x

+ O(∆x

) (A.14)

−∆x

+−∆x

+ O(∆x

) (A.15)

If we add together 3/8 × (A.13) + 6/8 × (A.14) − 1/8 × (A.15) we obtain

−

+ O(∆x

) (A.16)

The terms involving orders ∆x and ∆x

cancel out in this uniform grid and

the QUICK scheme is a third-order accurate approximation.

If necessary formulae involving more points (in backward or forward

directions) can be derived which have higher-order accuracy: see texts such

as Abbott and Basco (1989) or Fletcher (1991) for further details.

∂

∂φ

∂

∂φ

∂

∂φ

∂

−

2(∆x/2)

−

∆x

∂φ

∂

−

2∆x

∂φ

∂

ANIN_Z01.qxd 29/12/2006 10:06 AM Page 447

For the sake of simplicity the worked examples have focused on uniform

grids of nodal points. However, the derivation of discretised equations in

Chapters 4 and 5 used general geometrical dimensions such as

etc. and is also valid for non-uniform grids. In a non-uniform grid the faces

e and w of a general node may not be at the midpoints between nodes E

and P, and nodes W and P, respectively. In this case the interface values of

diffusion coefﬁcients Γ are calculated as follows:

= (1 − f

)Γ

+ f

(B.1a)

where the interpolation factor f

is given by

= (B.1b)

and

= (1 − f

)Γ

+ f

(B.1c)

where

= (B.1d)

For a uniform grid these expressions simplify to equations (4.5a–b): since

= 0.5, f

= 0.5, we have Γ

= (Γ

+Γ

)/2 and Γ

= (Γ

+Γ

)/2.

Basically there are two practices used to locate control volume faces in

non-uniform grids (Patankar, 1980):

Practice A: Nodal points are deﬁned ﬁrst and the control volume faces

located mid-way between the grid points. This is illustrated in Figure B.1.

Appendix B Non-uniform grids

Figure B.1

ANIN_Z02.qxd 29/12/2006 10:06 AM Page 448

APPENDIX B 449

Here the faces of a control volume are not at the midpoint between the

nodes. The evaluation of gradients obtained through a linear approximation

is unaffected because the gradient remains the same at any point between the

nodes in question, but the values of diffusion coefﬁcient Γ need to be evalu-

ated using interpolation functions (B.1).

It is very important to note that central difference formulae for the

calculation of gradients at cell faces and the QUICK scheme for convective

ﬂuxes are only second- and third-order accurate respectively when the con-

trol volume face is mid-way between nodes. In practice A a control volume

face, e for example, is mid-way between nodes P and E, so the differencing

formula used to evaluate the gradient (

∂φ

∂

is second-order accurate.

A further advantage of practice A is that property values Γ

, Γ

etc. can be

easily evaluated by taking the average values. The disadvantage of practice A

is that the value of the variable

at P may not necessarily be the most repres-

entative value for the entire control volume as point P is not at the centre of

the control volume. In practice B the value of

at P is a good representative

value for the control volume as P lies at the centre of the control volume, but

the discretisation schemes lose accuracy. A thorough discussion of these two

practices can be found in Patankar (1980), to which the reader is referred for

further details.

Figure B.2

Practice B: Locations of the control volume faces are deﬁned ﬁrst and the

nodal points are placed at the centres of the control volumes. This is illus-

trated in Figure B.2.

ANIN_Z02.qxd 29/12/2006 10:06 AM Page 449

Source terms of the discretised equations are evaluated using prevailing

values of the variables. Since a staggered grid is employed in the calculation

procedures, interpolation is required in the calculation of velocity gradient

terms which often appear in source terms. For example, consider the two-

dimensional u-momentum equation

(

uu) + (

vu) =

−+S

(C.1)

where

A small part of the solution grid is shown in Figure C.1; we use the standard

notation for the backward staggered velocity components introduced in

Chapter 6.

∂

Appendix C Calculation of source terms

Figure C.1

ANIN_Z03.qxd 29/12/2006 10:07 AM Page 450

APPENDIX C 451

The discretised form of equation (C.1) for the u-control volume centred

at (i, J ) is

i, J

=∑a

−∆V

+ D

∆V

(C.2)

where

is the width of the u-control volume and ∆V

is the volume of the

u-control volume.

The source term for the u-velocity equation at the cell shown is evaluated

∆V =

Cell

∆V

−

= ∆V

−

(C.3)

Source terms of other equations are calculated in a similar manner.

I, j

− v

I−1, j

∂

I, j+1

− v

I−1, j+1

∂

i, J

− u

i−1, J

∂

i+1, J

− u

i, J

∂

I, J

− P

I−1, J

)

ANIN_Z03.qxd 29/12/2006 10:07 AM Page 451

1 Van Leer limiter (Van Leer, 1974):

(r) = (D.1)

Figure D.1 shows this limiter on a r−

diagram. The limiter function

stays in the TVD region and goes through the point (1, 1), so it is

therefore a second-order accurate TVD limiter. The limiter tends

= 2 for large values of r and is symmetrical.

r +|r |

1 + r

Figure D.1 Van Leer limiter

2 Van Albada limiter (Van Albada et al, 1982):

(r) = (D.2)

This function is shown in Figure D.2. This limiter also goes through

(1, 1) and stays well below the upper limit for TVD schemes. As can be

seen from Figure D.2, this limiter tends to

= 1 as r →∞. It is

straightforward to check that the limiter is symmetrical.

r + r

1 + r

For details of references see Chapter 5.

Appendix D Limiter functions used in

Chapter 5

ANIN_Z04.qxd 29/12/2006 10:08 AM Page 452

APPENDIX D 453

Figure D.2 Van Albada limiter

Figure D.3 Min–Mod limiter

4 SUPERBEE limiter of Roe (1983):

(r) = max[0, min(2r, 1), min(r, 2)] (D.4)

Figure D.4 shows that this limiter is also a combination of linear

expressions. The limiter is symmetrical and follows the upper limit

of the TVD region on the r−

diagram.

3 Min–Mod limiter:

(r) =

min(r, 1) if r > 0

0 if r ≤ 0 (D.3)

Figure D.3 shows that this limiter is piecewise linear and follows the

lowest boundary for a second-order accurate TVD scheme on a r−

diagram. The limiter is symmetrical.

ANIN_Z04.qxd 29/12/2006 10:08 AM Page 453

454 APPENDIX D

Figure D.4 SUPERBEE limiter

Figure D.5 Sweby limiter

6 QUICK limiter of Leonard (Leonard, 1988): A monotonic TVD version

of the QUICK scheme can be achieved by the following function

(see Lien and Leschziner, 1993):

(r) = max[0, min(2r, (3 + r)/4, 2)] (D.6)

5 Sweby limiter (Sweby, 1984): This limiter represents a generalisation of

the Min–Mod and SUPERBEE limiters by means of a single parameter

. The limiter function is as follows:

(r) = max[0, min(

r, 1), min(r,

)] (D.5)

The limiter is symmetrical. We consider only the range of values

1 ≤

≤ 2, and note that Sweby’s limiter becomes the Min–Mod

limiter when

= 1 and the SUPERBEE limiter of Roe when

= 2.

The limiters span the whole TVD region between its upper and

lower limits. Figure D.5 shows the r−

diagram for Sweby’s limiter

when

= 1.5.

ANIN_Z04.qxd 29/12/2006 10:08 AM Page 454

APPENDIX D 455

Figure D.6 QUICK limiter

Figure D.7 UMIST limiter

Figure D.6 shows the r−

diagram for this limiter, which is not

symmetric.

7 UMIST limiter (Lien and Leschziner, 1993): UMIST, which stands

for Upstream Monotonic Interpolation for Scalar Transport, is a

symmetrical form of the QUICK limiter. The function is given by

(r) = max[0, min(2r, (1 + 3r)/4, (3 + r)/4, 2)] (D.7)

Figure D.7 shows this limiter on the r−

diagram.

ANIN_Z04.qxd 29/12/2006 10:08 AM Page 455

We develop one-dimensional equations for frictionless, incompressible ﬂow

by application of the mass conservation and x-momentum equations from

Chapter 2 to a small control volume ∆V of length ∆x in the ﬂow direction.

Mass conservation:

+ div(

V ) = 0 (2.4)

Since the ﬂow is steady we have: div(

V) = 0 (E.1)

We integrate over the control volume and apply Gauss’ divergence theorem

(2.41) to convert the resulting volume integral into an integral over the

bounding surfaces of the control volume:

div(

V)dV =

V . ndA +

V . ndA = 0 (E.2)

There is no mass ﬂow across the side walls, so the last integral over ∆A

equal to zero.

The velocity distribution is uniform and the velocity only changes as

a function of the streamwise x-co-ordinate. The velocity vector V

is in

the positive x-direction, so V

=+u

i, where u

is the magnitude of the

x-velocity at A

; similarly V

=+u

i, where u

is the magnitude of the x-

velocity at A

. The outward normal vector n

is in the negative x-direction,

so n

=−i. The outward normal vector n

on the other hand is in the posi-

tive x-direction, so n

=+i. Thus,



∆A



∆V

∂ρ

∂

Appendix E Derivation of one-dimensional

governing equations for steady,

incompressible flow through

a planar nozzle

Figure E.1

ANIN_Z05.qxd 29/12/2006 10:08 AM Page 456