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

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

6.10 WORKED EXAMPLES OF THE SIMPLE ALGORITHM 207

0.17769p′

= 0.11909p′

+ 0.058593p′

+ 0.18279

0.081093p′

= 0.058593p′

+ 0.02249p′

+ 0.25882

Nodes are A and E are boundary nodes so they need special treatment.

• Pressure nodes A and E

The pressure corrections are set to zero for both nodes:

p′

= 0.0

p′

= 0.0

At node E this follows the practice of Example 6.1, because the static

pressure is given at the nozzle exit. If the static pressure p

at inlet was

given this would also apply at node A without reservation. However, in

this problem we are working with a given stagnation pressure, so we

need to be careful. We note that p

is ﬁxed by Equation (6.65) if the

stagnation pressure p

and velocity u

are known. At the stage in the

SIMPLE algorithm where we start to solve the pressure correction

equations, we have available the guessed velocity u*

as a result of solving

the discretised momentum equation. Whilst it is true that this velocity is

constantly updated as the iterations proceed, we may regard that at each

iteration level the static pressure p

is temporarily ﬁxed by p

and the

current value of u*

, thus justifying the use of p′

= 0.0.

Substitution of p′

= 0.0 and p ′

= 0.0 into the pressure correction equations

for internal nodes B, C and D yields the following system of equations:

0.26161p′

= 0.11909p′

− 0.06830

0.17769p′

= 0.11909p′

+ 0.058593p′

+ 0.18279

0.081093p′

= 0.058593p′

+ 0.25882

These three equations can be solved to give the pressure correction at nodes

B, C and D. The resulting solution is

p′

0.0 1.63935 4.17461 6.20805 0.0

Corrected nodal pressures are now calculated using these pressure

corrections:

= p*

+ p′

= 7.5 + 1.63935 = 9.13935

= p*

+ p′

= 5.0 + 4.17461 = 9.17461

= p*

+ p′

= 2.5 + 6.20805 = 8.70805

Corrected velocities at the end of the ﬁrst iteration are

= u*

+ d

(p′

− p′

) = 2.19935 + 0.31670 × [0.0 − 1.63935] = 1.68015 m/s

= u*

+ d

(p′

− p′

) = 3.02289 + 0.34027 × [1.63935 − 4.17461] = 2.16020 m/s

= u*

+ d

(p′

− p′

) = 3.50087 + 0.23437 × [4.17461 − 6.20805] = 3.02428 m/s

= u*

+ d

(p′

− p′

) = 4.10926 + 0.15 × [6.20805 − 0.0] = 5.04047 m/s

We can also calculate the corrected nodal pressure at A using equation (6.65):

= p

−

–

)

= 10 −

–

× 1.0 × (1.68015 × 0.45/0.5)

= 8.85671

ANIN_C06.qxd 29/12/2006 09:59 AM Page 207

208 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING

First, we check whether the velocity ﬁeld satisﬁes continuity. The mass ﬂow

rates

uA calculated at u-nodes are

Continuity check

Node 1234

uA 0.75607 0.75607 0.75607 0.75607

The exact mass ﬂow rate from Bernoulli’s equation is 0.44721 kg/s, so the

error in the computed mass ﬂow rate is +69%. This is not a problem, because

we would not expect the solution after one iteration to be accurate. Never-

theless, the fact that the mass ﬂow rate at all four velocity nodes is exactly the

same highlights an important feature of SIMPLE, which also applies in

more complex multi-dimensional problems: the algorithm aims to supply

a continuity-satisfying velocity ﬁeld at the end of each iteration cycle. The

close link with this key conservation principle has proved to be a major

strength of the SIMPLE algorithm and its variants.

The computed velocity solution at the end of an iteration cycle is not yet

in balance with the computed pressure ﬁeld, i.e. momentum is not yet con-

served. Of course this is due to the fact that the entries in the discretised

momentum equations were computed on the basis of an assumed initial

velocity ﬁeld. Therefore, we need to perform iterations until both continuity

and momentum equations are satisﬁed.

Under-relaxation

In the iteration process the SIMPLE algorithm requires under-relaxation.

For the next iteration we use under-relaxation factors for both velocity and

pressure (say 0.8 for both) and start the next solution cycle with the follow-

ing velocity and pressure ﬁelds:

new

= (1 − 0.8) × u

old

+ 0.8 × u

calculated

new

= (1 − 0.8) × p

old

+ 0.8 × p

calculated

The velocity ﬁeld for the next iteration is

m/s u

m/s

1.78856 2.29959 3.21942 5.36571

As explained in section 6.4, equations (6.36) and (6.37), a

, S

and d values

of the discretised momentum equations are also under-relaxed. Note that,

for illustration purposes, these under-relaxation measures were not included

in the a

, S

and d values of the discretised momentum equations shown

earlier. In practice under-relaxation is used from the start of the calculation:

hence the resulting solution with the above under-relaxation measures would

be slightly different from the one shown.

Iterative convergence and residuals

If we substitute the under-relaxed velocity and pressure ﬁelds into the

discretised momentum equations they will not satisfy the equations unless

by chance we have computed the ﬁnal solution in one iteration (e.g. due to

ANIN_C06.qxd 29/12/2006 09:59 AM Page 208

6.10 WORKED EXAMPLES OF THE SIMPLE ALGORITHM 209

a fortunate choice of initial velocity and pressure ﬁelds). For example, the

discretised momentum equation at node 1 in the next iteration is

1.20425u

= 1.98592

The difference between the left and right hand sides of the discretised

momentum equation at every velocity node is called the momentum

residual. Substituting the current velocity value of u

= 1.78856 yields:

u-momentum residual at node 1 = 1.20425 × 1.78856 − 1.98592 = 0.16795

If the iteration sequence is convergent this residual should decrease to show

an improving balance between the computed velocity and pressure ﬁelds.

Ideally, we would like to stop the iteration process when mass and momen-

tum are exactly balanced in the discretised pressure correction and momen-

tum equations. In practice, the ﬁnite precision of number representation

in computing machinery would make this impossible and, even if it were

possible to compute with very high precision, this would take far too much

computing time. Our aim is to truncate the iterative sequence when we are

sufﬁciently close to exact balance that further improvement is of limited

practical value.

We calculate momentum residuals at all velocity nodes and monitor the

sum of absolute values of the residuals as an indication of satisfactory progress

of the calculation sequence. We note that residuals can be positive as well

as negative. Using the sum of absolute values prevents spurious indication

of convergence due to cancellation between positive and negative contribu-

tions. We accept the solution when the sum of absolute residuals is less than

a predetermined small value (say 10

−5

). It should be noted that this is a weak

test to determine the point where the iterative sequence can be truncated.

A decreasing sum of residuals could be due to residuals that decrease by

roughly the same amount at every node or due to a small number of decreas-

ing residuals in conjunction with others that do not decrease much at all. In

a grid with a large number of nodes a few increasing residuals might even

be hidden amongst a much larger number of strongly decreasing residuals.

Nevertheless, summed residuals calculations are routinely used as conver-

gence indicator in ﬂuid ﬂow calculations. For a further discussion of the use

of residuals and iterative convergence the reader is referred to Chapter 10.

Application of under-relaxation factors of 0.8 for both u and p and allowing

a maximum sum of absolute momentum residuals of 10

−5

yields convergence

in 19 iterations. The solution is given in the table below

Converged pressure and velocity ﬁeld after 19 iterations

Pressure (Pa) Velocity (m/s)

Node Computed Exact Error (%) Node Computed Exact Error (%)

A 9.22569 9.60000 −3.9 1 1.38265 0.99381 39.1

B 9.00415 9.37500 −4.0 2 1.77775 1.27775 39.1

C 8.25054 8.88889 −7.2 3 2.48885 1.78885 39.1

D 6.19423 7.50000 −17.4 4 4.14808 2.98142 39.1

E0 0 –

Solution

ANIN_C06.qxd 29/12/2006 09:59 AM Page 209

210 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING

Figure 6.12 Predicted mass ﬂow

rate with different grids

The converged mass ﬂow rate for this ﬁve-node grid is 0.62221 kg/s, which

is 39% higher than the exact value. By reﬁning the grid we can get progres-

sively closer solution to the exact solution. Using grids with 10, 20 and

50 nodal points yields mass ﬂow rates of 0.5205, 0.4805 and 0.4597 kg/s,

respectively. This illustrates how the error in the solution can be reduced by

systematic grid reﬁnement. If the grid is further reﬁned to say 200, 500 and

1000 grid nodes the computed mass ﬂow rate converges to the exact solution

of 0.44721 kg/s. This is graphically illustrated in Figure 6.12.

Some comments on exit boundary conditions

Finally, we brieﬂy discuss the issue of the downstream boundary condition.

In Example 6.1 the exit pressure was set to p

= 0. Solution of the pressure

correction equations yields the pressures at nodes other than D as gauge

pressures (relative to p

). Say the absolute pressure at D had a non-zero

reference value p

Abs,D

= p

ref

at this location. The absolute pressure ﬁeld

at node N can be found by adding the absolute pressure at D to the gauge

pressure at N: p

Abs,N

= p

Ref

+ p

Gauge,N

. If the absolute pressure is known at

some other reference location R in the computational domain (p

Abs,R

= p

Ref

)

the absolute pressure at node N can be obtained by means of p

Abs,N

= p

Ref

Gauge,N

– p

Gauge,R

). In constant-property ﬂows the actual value of p

ref

is imma-

terial, since pressure differences appear in the source terms in the discre-

tised momentum equations. When we solve problems involving ﬂuids with

properties that depend on the absolute pressure (e.g. compressible ﬂows) we

modify the SIMPLE algorithm by including one more iterative structures to

update the ﬂuid properties as new absolute pressures become available.

As we have seen in section 2.10 it is also possible to use an outﬂow

boundary condition at the downstream boundary in conjunction with a given

inlet velocity. The outﬂow boundary condition requires the gradient of the

velocity to be zero at the downstream boundary. This can be implemented

by equating the velocities at the two nodes that straddle the boundary, i.e. by

setting

= u

(6.70)

ANIN_C06.qxd 29/12/2006 09:59 AM Page 210

6.11 SUMMARY 211

In Example 6.1 we have seen that a ﬁxed pressure boundary condition is

implemented by setting the pressure correction to zero, which reduces the

original system of pressure correction equations by one equation. Equation

(6.70) cannot provide a pressure, so if this zero velocity gradient boundary

condition was used in Example 6.1 we would have only three pressure

correction equations for nodes A, B and C, but four (unknown) pressure

corrections p′

, p′

. Thus, it appears as if we are one equation short.

To overcome this problem we note again that pressure differences are

important in the discretised momentum equation and use the above device

of setting an arbitrary reference pressure at the exit plane: p

= p

ref

. For the

sake of expediency it is easiest to set p

= p

ref

= 0. Having ﬁxed the pressure

we may set the pressure correction equal to zero, solve the pressure correc-

tion equation as in the above examples and obtain the pressure solution as a

gauge pressure.

The most popular solution algorithms for pressure and velocity calculations

with the ﬁnite volume method have been discussed. They all possess the

following common characteristics:

• The problems associated with the non-linearity of the momentum

equations and the coupling between transport equations are tackled by

adopting an iterative solution strategy.

• Velocity components are deﬁned on staggered grids to avoid problems

associated with pressure ﬁeld oscillations of high spatial frequency.

• In the staggered grid arrangement velocities are stored at the cell faces

of scalar control volumes. The discretised momentum equations are

solved on staggered control volumes whose cell faces contain the

pressure nodes.

• The SIMPLE algorithm is an iterative procedure for the calculation of

pressure and velocity ﬁelds. Starting from an initial pressure ﬁeld p* its

principal steps are:

– solve discretised momentum equation to yield intermediate velocity

ﬁeld (u*, v*)

– solve the continuity equation in the form of an equation for pressure

correction p′

– correct pressure and velocity by means of

I, J

= p*

I, J

+ p′

I, J

i, J

= u*

i, J

+ d

i, J

(p′

I−1, J

− p′

I, J

)

I, j

= v*

I, j

+ d

I, j

(p′

I, J −1

− p′

I, J

)

– solve all other discretised transport equations for scalars

– repeat until p, u, v and

ﬁelds have all converged.

• Reﬁnements to SIMPLE have produced more economical and stable

iteration methods.

• The steady state PISO algorithm adds an extra correction step to

SIMPLE to enhance its performance per iteration.

• It is unclear which of the procedures is the best for general-purpose

computation.

• Under-relaxation is required in all methods to ensure stability of the

iteration process.

Summary6.11

ANIN_C06.qxd 29/12/2006 09:59 AM Page 211

In the previous chapters we have discussed methods of discretising the

governing equations of ﬂuid ﬂow and heat transfer. This process results in

a system of linear algebraic equations which needs to be solved. The com-

plexity and size of the set of equations depends on the dimensionality of the

problem, the number of grid nodes and the discretisation practice. Although

any valid procedure can be used to solve the algebraic equations, the avail-

able computer resources set a powerful constraint. There are two families

of solution techniques for linear algebraic equations: direct methods and

indirect or iterative methods. Simple examples of direct methods are

Cramer’s rule matrix inversion and Gaussian elimination. The number of

operations to the solution of a system of N equations with N unknowns by

means of a direct method can be determined beforehand and is of the order

of N

. The simultaneous storage of all N

coefﬁcients of the set of equations

in core memory is required.

Iterative methods are based on the repeated application of a relatively

simple algorithm leading to eventual convergence after a – sometimes large

– number of repetitions. Well-known examples are the Jacobi and Gauss–

Seidel point-iterative methods. The total number of operations, typically of

the order of N per iteration cycle, cannot be predicted in advance. Stronger

still, it is not possible to guarantee convergence unless the system of equations

satisﬁes fairly exacting criteria. The main advantage of iterative solution

methods is that only non-zero coefﬁcients of the equations need to be stored

in core memory.

The one-dimensional conduction example in Chapter 4, section 4.3, led

to a tri-diagonal system – a system with only three non-zero coefﬁcients per

equation. When QUICK differencing is applied to a convection–diffusion

problem this gives rise to a penta-diagonal system that has ﬁve non-zero

coefﬁcients, which is somewhat more complex to deal with. Nevertheless,

the ﬁnite volume method usually yields systems of equations, each of which

has a vast majority of zero entries. Since the systems arising from realistic

CFD problems can be very large – up to 100 000 to 1 million equations – we

ﬁnd that iterative methods are generally much more economical than direct

methods.

Thomas (1949) developed a technique for rapidly solving tri-diagonal

systems that is now called the Thomas algorithm or the tri-diagonal matrix

algorithm (TDMA). The TDMA is actually a direct method for one-

dimensional situations, but it can be applied iteratively, in a line-by-line

fashion, to solve multi-dimensional problems and is widely used in CFD

programs. It is computationally inexpensive and has the advantage that it

Chapter seven

Solution of discretised equations

Introduction7.1

ANIN_C07.qxd 29/12/2006 04:48PM Page 212

7.2 THE TDMA 213

requires a minimum amount of storage. In this chapter the TDMA is

explained in detail in sections 7.2 to 7.5.

The Jacobi and Gauss–Seidel methods are general-purpose point-iterative

algorithms that are easily implementable, but their convergence rate can be

slow when the system of equations is large. They were not initially consid-

ered suitable for general-purpose CFD procedures. However, more recently

multigrid acceleration techniques have been developed that have improved

the convergence rates of iterative solvers to such an extent that they are now

the method of choice in commercial CFD codes. Point-iterative techniques

and multigrid accelerators will be described in sections 7.6 and 7.7. This

chapter closes with a brief discussion of alternative methods.

Consider a system of equations that has a tri-diagonal form

= C

(7.1a)

−

+ D

−

= C

(7.1b)

−

+ D

−

= C

(7.1c)

−

+ D

−

= C

......= ..

−

n−1

+ D

−

n+1

= C

(7.1n)

n+1

= C

n+1

(7.1n+1)

In the above set of equations

and

n+1

are known boundary values. The

general form of any single equation is

−

j−1

+ D

−

j+1

= C

(7.2)

Equations (7.1b–n) of the above set can be rewritten as

+ (7.3a)

+ (7.3b)

+ (7.3c)

n+1

n−1

+ (7.3n−1)

These equations can be solved by forward elimination and back-substitution.

The forward elimination process starts by removing

from equation

(7.3b) by substitution from equation (7.3a) to give

++C

+ (7.4a)

−

The TDMA7.2

ANIN_C07.qxd 29/12/2006 04:48PM Page 213

If we adopt the notation

= and C ′

+ (7.4b)

equation (7.4a) can be written as

+ (7.4c)

If we let

= and C ′

equation (7.4c) can be recast as

= A

+ C ′

(7.5)

Formula (7.5) can now be used to eliminate

from (7.3c) and the procedure

can be repeated up to the last equation of the set. This constitutes the

forward elimination process.

For the back-substitution we use the general form of recurrence

relationship (7.5):

= A

j+1

+ C ′

(7.6a)

where

= (7.6b)

C ′

= (7.6c)

The formulae can be made to apply at the boundary points j = 1 and j = n + 1

by setting the following values for A and C ′:

= 0 and C ′

n+1

= 0 and C ′

n+1

In order to solve a system of equations it is ﬁrst arranged in the form of

equation (7.2) and

, D

and C

are identiﬁed. The values of A

and C ′

are

subsequently calculated starting at j = 2 and going up to j = n using (7.6b–c).

Since the value of

is known at boundary location (n + 1) the values for

can be obtained in reverse order (

n−1

n−2

,...,

) by means of the

recurrence formula (7.6a). The method is simple and easy to incorporate into

CFD programs. A Fortran subroutine for the TDMA is given in Anderson

et al. (1984).

In the above derivation of the TDMA we assumed that boundary values

and

n+1

were given. To implement a ﬁxed gradient (or ﬂux) boundary

condition, e.g. at j = 1, the coefﬁcient

in equation (7.1b) is set to zero and

the ﬂux across the boundary is incorporated in source term C

. The actual

value of the variable at the boundary is now not directly used in the formu-

lation. The absence of the ﬁrst or the last value does not pose a problem in

applying the TDMA, as will be illustrated in examples below.

C ′

j−1

+ C

−

j−1

−

j−1

C ′

+ C

−

C′

+ C

−

214 CHAPTER 7 SOLUTION OF DISCRETISED EQUATIONS

ANIN_C07.qxd 29/12/2006 04:48PM Page 214

7.4 APPLICATION OF THE TDMA TO 3D PROBLEMS 215

The TDMA can be applied iteratively to solve a system of equations for

two-dimensional problems. Consider the grid in Figure 7.1 and a general

two-dimensional discretised transport equation of the form

= a

+ a

+ b (7.7)

Figure 7.1 Line-by-line

application of the TDMA

To solve the system the TDMA is applied along chosen, e.g. north–south

(n–s), lines. The discretised equation is rearranged in the form

−a

+ a

− a

= a

+ a

+ b (7.8)

The right hand side of (7.8) is assumed to be temporarily known. Equation

(7.8) is in the form of equation (7.2) where

≡ a

, D

≡ a

and C

≡

+ a

+ b. Now we can solve along the n–s direction of the chosen line

for values j = 2, 3, 4,..., n as shown in Figure 7.1.

Subsequently the calculation is moved to the next north–south line. The

sequence in which lines are moved is known as the sweep direction. If we

sweep from west to east the values of

to the west of a point P are known

from the calculations on the previous line. Values of

to its east, however,

are unknown so the solution process must be iterative. At each iteration cycle

is taken to have its value at the end of the previous iteration or a given

initial value (e.g. zero) at the ﬁrst iteration. The line-by-line calculation pro-

cedure is repeated several times until a converged solution is obtained.

For three-dimensional problems the TDMA is applied line by line on a

selected plane and then the calculation is moved to the next plane, scanning

the domain plane by plane. For example, to solve along a n–s line in the x–y

plane of Figure 7.2, a discretised transport equation is written as

−a

+ a

− a

= a

+ a

+ b (7.9)

Application of

the TDMA to two-

dimensional problems

7.3

Application of

the TDMA to

three-dimensional

problems

7.4

ANIN_C07.qxd 29/12/2006 04:48PM Page 215

216 CHAPTER 7 SOLUTION OF DISCRETISED EQUATIONS

Figure 7.3 The grid for

Example 7.1

The values at W and E as well as those at B and T on the right hand side of

Equation (7.9) are considered to be temporarily known. Using the TDMA,

values of

along a selected north–south line are computed. The calculation

is moved to the next line and subsequently swept through the whole plane

until all unknown values on each line have been calculated. After completion

of one plane the process is moved on to the next plane.

Figure 7.2 Application of the

TDMA in a three-dimensional

geometry

In two- and three-dimensional computations the convergence can often

be accelerated by alternating the sweep direction so that all boundary

information is fed into the calculation more effectively. To solve along an

east–west line in the present three-dimensional case the discretised equation

is rearranged as follows:

−a

+ a

− a

= a

+ a

+ b (7.10)

We consider the one-dimensional steady state conductive/convective heat

transfer from a bar of material discussed ﬁrst in Example 4.3 of section 4.3. The

geometry is shown in Figure 7.3. The temperature on the left hand boundary

is taken to be 100°C and the right hand side boundary is insulated so the heat

ﬂux across it is zero. Heat is lost to the surroundings by convective heat

transfer. Solve the matrix equation (4.49) for this problem using the TDMA.

The matrix equation found in section 4.3 was

G20 −5000JG

JG1100J

H−515−500KH

KH100K

H 0 −515−50KH

K = H 100K (4.49)

H 00−515−5KH

KH100K

I 000−510LI

LI100L

The general form of the equation used in the TDMA is

−

j−1

+ D

−

j+1

= C

(7.2)

Examples7.5

Solution

Example 7.1

An illustration

of the TDMA in

one dimension

ANIN_C07.qxd 29/12/2006 04:48PM Page 216