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

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

7.7 MULTIGRID TECHNIQUES 237

Step 3: Prolongation

The next step is to go backwards transferring error vectors from each coarse

grid to the next ﬁne grid level. This is called the prolongation process. Linear

interpolation or any other interpolation scheme could be used to construct

ﬁne grid values from the coarse grid values. Using linear interpolation some

sample values are

e′

= (0.75e

)

e′

= (0.75e

+ 0.25e

) (7.41)

e′

= (0.25e

+ 0.75e

)

and so on.

It should be noted that we use a prime to indicate the prolonged error vec-

tor e′

on Grid 2 to distinguish it from the error vector e

. Furthermore, for

the nodes closest to the boundary we have used the fact that the value of the

problem variable is known so the error on the boundary is zero. The calcu-

lations of (7.41) yield the following values for the prolonged error vector e′

Ge′

JG17.556J

He′

KH31.514K

He′

KH47.726K

H . K

H . KH. K

He′

KH24.062K

Ie′

LI12.261L

The prolonged error vector e′

is now used to correct the original error vec-

tor e

(7.40) on Grid 2:

corrected

= e

+ e′

(7.42)

This yields

JG19.156JG17.556JG36.713J

KH58.310KH31.514KH89.825K

KH96.049KH47.726KH143.775K

H . K

H . KH . KH. KH . K

KH158.591KH24.062KH182.654K

LI55.351LI12.261LI67.612L

At this stage it is usual to do some smoothing iterations before transferring

this error vector to the grid above this level. First we perform two Gauss–

Seidel smoothing iterations and obtain the following corrected and smoothed

error vector on Grid 2:

JG32.639J

KH95.749K

KH152.494K

H . K

H . KH. K

KH188.283K

LI65.248L

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

This result is next used to compute the prolonged error vector e

on Grid 1

by linear interpolation using the process outlined in (7.41), but replacing

superscripts 2h by h and 4h by 2h:

JG24.479J

KH48.416K

KH79.971K

H . K

H . KH. K

KH96.007K

LI48.936L

Step 4: Correction and ﬁnal iterations

Finally, the prolonged error vector e

is used to compute the corrected inter-

mediate solution y on Grid 1:

corrected

= y + e

(7.43a)

Thus,

JG116.755JG24.479JG141.235J

KH141.994KH48.416KH190.411K

KH160.427KH79.971KH240.399K

H . K

(7.43b)

H . KH . KH. KH.K

KH394.392KH96.007KH490.399K

LI468.130LI48.936LI517.067L

Comparison of the corrected solution (7.43b) with intermediate solution

(7.38) and the TDMA solution in Table 7.14 shows that the multigrid pro-

cedure has considerably reduced the error. Substitution of the corrected

solution into r = b − A.ygives an r.m.s. residual of 8.786, which is lower

than the previous r.m.s. residual on Grid 1, which was 14.951. Since inter-

polation errors are involved in the restriction and prolongation processes we

cannot expect to achieve the true solution in one multigrid cycle. In order to

improve the solution further we do more iterations on the ﬁne mesh (say

two) and repeat the ‘ﬁne grid–coarse grid’ procedure until convergence is

achieved. We proceed by using the three-grid procedure and go backwards

and forwards as many times as is needed to reduce the r.m.s. residual to

−6

. This multigrid cycle is called a three-grid V-cycle. The process steps

are illustrated in Figure 7.8 along with annotations of the number of itera-

tions at each level inside the circles. The diagram also reveals the origin of

the term ‘V-cycle’.

238 CHAPTER 7 SOLUTION OF DISCRETISED EQUATIONS

Figure 7.8 Schematic of the V-

cycle multigrid procedure used

in the example

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

7.7 MULTIGRID TECHNIQUES 239

The pattern of convergence achieved by repeating the V-cycle of 2 ﬁne

grid iterations, 10 and 10 coarse grid iterations is shown in Figure 7.9. The

multigrid procedure is fast and effective, since it converges in 60 ﬁne grid

iterations, which compares favourably with the ordinary Gauss–Seidel

method, which takes 664 iterations to achieve the same residual value. Even

after allowing for extra computational effort due to the coarse grid iterations,

the order-of-magnitude improvement of the convergence rate by the multi-

grid procedure is clearly beneﬁcial. When multigrid acceleration techniques

are applied to 2D and 3D problems the convergence gains obtained are very

attractive, which explains their popularity among CFD users.

Figure 7.9 Residual

reduction pattern with

ordinary Gauss–Seidel iterations

and multigrid Gauss–Seidel

iterations

7.7.3 Multigrid cycles

Multigrid techniques can be used in conjunction with any iterative technique.

In our simple example we have illustrated the main concepts of the multigrid

methods. In practical CFD calculations the multigrid transfer process is

more sophisticated and different cycles of coarsening and reﬁnement are

used with special schedules of restriction and prolongation at different

reﬁnement levels. Common choices of multigrid cycles are the so-called V-,

W- and F-cycles, which are illustrated in Figure 7.10.

The simple V-cycle shown in Figure 7.10a consists of two legs. The cal-

culation starts at the ﬁnest level. Iterations at any level are called relaxation.

After a few relaxation sweeps on the ﬁnest level the residuals are restricted

to the next coarse level and after relaxation on that level the residuals are

passed on to the next coarse level, and so on until the coarsest level is

reached. After ﬁnal relaxation on the coarsest level the prolongation steps are

performed on the upward leg of the V-cycle until the ﬁnest level is reached.

In the W-cycle additional restriction and prolongation sweeps are used at

coarser levels to obtain better reduction of long-wavelength errors. A typical

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

pattern is illustrated in Figure 7.10b. The ﬂexible cycle or ‘F-cycle’ is very

similar to the W-cycle, but has a different pattern of coarse-level sweeps as

illustrated in Figure 7.10c.

In the technique known as the full multigrid (FMG) method the calcula-

tions do not start at the ﬁnest grid, but instead at the coarsest level. Solutions

are transferred to successively ﬁner grid levels and on the ﬁnest level the

prolonged solution is used as the initial guess for start of the iterative pro-

cess. This solution process could be accelerated further using any of the

cycling procedures. For example, V-cycles could be used at each successive

reﬁnement level, as illustrated in Figure 7.11.

240 CHAPTER 7 SOLUTION OF DISCRETISED EQUATIONS

Figure 7.10 Illustration of

different multigrid cycle

strategies: (a) V-cycle,

(b) W-cycle and (c) F-cycle

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

7.7 MULTIGRID TECHNIQUES 241

7.7.4 Grid generation for the multigrid method

As illustrated in the above example, grid generation is required to create the

coarse grids. The most straightforward method is to combine control vol-

umes or regenerate the mesh using half the number of nodes of the mesh

above that level. For 2D structured grids, such as the Cartesian grid shown

in Figure 7.12, coarse grids can be readily generated by deleting alternate

grid lines. Thus one coarse grid control is constructed from every four ﬁne

grid control volumes. This can easily be extended to 3D meshes using eight

ﬁne grid control volumes per coarse grid control volume.

Figure 7.11 Cycling strategy

used in full multigrid method

Figure 7.12 A 2D Cartesian

mesh – the coarse grid is

constructed by deleting alternate

grid lines or combining groups of

four control volumes

In the above example we computed the coarse grid system matrix and

other required quantities using actual geometrical properties of the coarse

mesh (Table 7.12). This type of multigrid procedure is called a geometric

multigrid procedure. In the other variation to this method, the coefﬁcients

are not recomputed from the grid geometry to save calculation effort, but

approximated as linear combinations of coefﬁcients of the ﬁne grid equations.

Such multigrid methods are called algebraic multigrid and are widely

used in commercial CFD solvers. The technique known as the additive

correction multigrid (ACM) strategy of Hutchinson and Raithby (1986) is

also a popular multigrid method used in many CFD procedures.

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

Multigrid acceleration of the Gauss–Seidel point-iterative method is cur-

rently the solution algorithm of choice for commercial CFD codes. The rate

of convergence of this procedure can be optimised by speciﬁc choices for:

(i) interpolation of the residual vector and coefﬁcient matrix from ﬁne to

coarse meshes during restriction, (ii) interpolation of the error vector from

coarse to ﬁne meshes during prolongation, (iii) cycles of coarsening and

reﬁnement with special schedules of restriction and prolongation at different

reﬁnement levels. For further details of more advanced multigrid procedures

the reader should consult the appropriate literature (see e.g. Wesseling,

1992; Briggs, 1987). There are also several excellent learning resources avail-

able on the Internet for multigrid methods: see for example the multigrid

network MGNET at http://www.mgnet.org/.

Several other solution algorithms are available for CFD problems with

discretised equations that contain a large number of contributions from sur-

rounding nodes. The Strongly Implicit Procedure (SIP) due to Stone (1968),

in particular with the improvements suggested by Schneider and Zedan

(1981), is more suitable in this case. Details are not presented here in the

interest of brevity and the interested reader is referred to Anderson et al.

(1984). Another solution procedure which is being used in CFD calculations

is the conjugate gradient method (CGM) of Hestenes and Stiefel (1952).

This method is based on matrix factorisation techniques. Improvements

by Reid (1971), Concus et al. (1976) and Kershaw (1978) ensure accelerated

convergence in the CFD calculations. The CGM requires greater storage

than other iterative methods described earlier. Further details of the method

can also be found in Press et al. (1992).

242 CHAPTER 7 SOLUTION OF DISCRETISED EQUATIONS

Summary7.8

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

Having ﬁnished the task of developing the ﬁnite volume method for steady

ﬂows we are now in a position to consider the more complex category of

time-dependent problems. The conservation law for the transport of a scalar

in an unsteady ﬂow has the general form

(

ρφ

) + div(

) = div(Γ grad

) + S

(8.1)

The ﬁrst term of the equation represents the rate of change term and is zero

for steady ﬂows. To predict transient problems we must retain this term in

the discretisation process. The ﬁnite volume integration of equation (8.1)

over a control volume (CV) must be augmented with a further integration

over a ﬁnite time step ∆t. By replacing the volume integrals of the convective

and diffusive terms with surface integrals as before (see section 2.5) and

changing the order of integration in the rate of change term we obtain

(

ρφ

) dt dV + n.(

)dA dt

= n.(Γ grad

)dA dt + S

dV dt (8.2)

So far we have made no approximations but to make progress we need tech-

niques for evaluating the integrals. The control volume integration is essen-

tially the same as in steady ﬂows and the measures explained in Chapters 4

and 5 are again adopted to ensure successful treatment of convection, diffusion

and source terms. Here we focus our attention on methods necessary for the

time integration. The process is illustrated below using the one-dimensional

unsteady diffusion (heat transfer) equation and is later extended to multi-

dimensional unsteady diffusion and convection–diffusion problems.

Unsteady one-dimensional heat conduction is governed by the equation

c = k + S (8.3)

∂



t +∆t



t +∆t



t +∆t



∂

t +∆t



∂

Chapter eight The finite volume method for

unsteady flows

Introduction8.1

One-dimensional

unsteady heat

conduction

8.2

ANIN_C08.qxd 29/12/2006 04:38PM Page 243

In addition to the usual variables we have c, the speciﬁc heat of the material

(J/kg.K).

Consider the one-dimensional control volume in Figure 8.1. Integration

of equation (8.3) over the control volume and over a time interval from t to

t+∆t gives

c dV dt = k dV dt + S dV dt (8.4)

This may be written as

c dt dV = kA

− kA

+ D∆V dt (8.5)

In equation (8.5), A is the face area of the control volume, ∆V is its volume,

which is equal to A∆x, where ∆x =

is the width of the control volume,

and D is the average source strength. If the temperature at a node is assumed

to prevail over the whole control volume, the left hand side can be written as

c dt dV =

c(T

− T

)∆V (8.6)

In equation (8.6) superscript ‘o’ to refers to temperatures at time t; tempera-

tures at time level t+∆t are not superscripted. The same result as (8.6)

would be obtained by substituting (T

− T

)/∆t for

∂

t, so this term has

been discretised using a ﬁrst-order (backward) differencing scheme. Higher-

order schemes, which may be used to discretise this term, will be discussed

brieﬂy later in this chapter. If we apply central differencing to the diffusion

terms on the right hand side equation (8.5) may be written as

c(T

− T

)∆V = k

A − k

A dt

+ D∆V dt (8.7)

t +∆t



− T

t +∆t



∂

t +∆t



t +∆t



∂

t +∆t



∂

t +∆t



t +∆t



∂



t +∆t



∂



t +∆t



244 CHAPTER 8 THE FINITE VOLUME METHOD FOR UNSTEADY FLOWS

Figure 8.1

ANIN_C08.qxd 29/12/2006 04:38PM Page 244

8.2 ONE-DIMENSIONAL UNSTEADY HEAT CONDUCTION 245

To evaluate the right hand side of this equation we need to make an assump-

tion about the variation of T

, T

and T

with time. We could use tempera-

tures at time t or at time t+∆t to calculate the time integral or, alternatively,

a combination of temperatures at time t and t+∆t. We may generalise the

approach by means of a weighting parameter

between 0 and 1 and write

the integral I

of temperature T

with respect to time as

= T

dt = [

+ (1 −

]∆t (8.8)

Hence

0 1/2 1

∆t

–

+ T

)∆tT

∆t

We have highlighted the following values of integral I

: if

= 0 the tem

perature at (old) time level t is used; if

= 1 the temperature at new time

level t+∆t is used; and ﬁnally if

= 1/2, the temperatures at t and t+∆t

are equally weighted.

Using formula (8.8) for T

and T

in equation (8.7), and dividing by A∆t

throughout, we have

c ∆x =

−

+ (1 −

) −+D∆x (8.9)

which may be rearranged to give

c +

+ T

= [

+ (1 −

] + [

+ (1 −

]

c − (1 −

) − (1 −

) T

+ D∆x (8.10)

Now we identify the coefﬁcients of T

and T

as a

and a

and write

equation (8.10) in the familiar standard form:

= a

[

+ (1 −

] + a

[

+ (1 −

]

+ [a

− (1 −

+ b (8.11)

where

+ a

) + a

∆x

∆t

∆x

∆t

− T

)

− T

)

− T

)

− T

)

− T

∆t

t +∆t



ANIN_C08.qxd 29/12/2006 04:38PM Page 245

246 CHAPTER 8 THE FINITE VOLUME METHOD FOR UNSTEADY FLOWS

and

with

D∆x

The exact form of the ﬁnal discretised equation depends on the value of

When

is zero, we only use temperatures T

, T

and T

at the old time

level t on the right hand side of equation (8.11) to evaluate T

at the new time

and the resulting scheme is called explicit. When 0 <

≤ 1 temperatures at

the new time level are used on both sides of the equation and the resulting

schemes are called implicit. The extreme case of

= 1 is termed fully

implicit and the case corresponding to

= 1/2 is called the Crank–

Nicolson scheme (Crank and Nicolson, 1947).

8.2.1 Explicit scheme

In the explicit scheme the source term is linearised as b = S

. Now

the substitution of

= 0 into (8.11) gives the explicit discretisation of the

unsteady conductive heat transfer equation

= a

+ a

+ [a

− (a

+ a

− S

)]T

+ S

(8.12)

where

= a

and

The right hand side of equation (8.12) only contains values at the old time

step so the left hand side can be calculated by forward marching in time. The

scheme is based on backward differencing and its Taylor series truncation

error accuracy is ﬁrst-order with respect to time. As explained in Chapter 5,

all coefﬁcients need to be positive in the discretised equation. The coefﬁcient

of T

may be viewed as the neighbour coefﬁcient connecting the values

at the old time level to those at the new time level. For this coefﬁcient to be

positive we must have a

− a

> 0. For constant k and uniform grid

spacing,

=∆x, this condition may be written as

∆x

∆t

∆x

∆t

ANIN_C08.qxd 29/12/2006 04:38PM Page 246