Parabolic multi-grid method for incompressible viscous flows using a group explicit relaxation scheme

Computers & Fluids Vol. 19, No. 1, pp. 3341, 1991

0045-7930/91 $3.00+ 0.00 Copyright© 1991PergamonPress plc

Printed in Great Britain.All rights reserved

PARABOLIC

MULTI-GRID

INCOMPRESSIBLE

VISCOUS

EXPLICIT

METHOD

FLOWS

RELAXATION

USING

FOR A GROUP

SCHEME

SHIGERU MURATA, ~ NOBUYUKI SATOFUKAt and TADASHIKUSHIYAMA2 ~Department of Mechanical and System Engineering, Kyoto Institute of Technology, Matsugasaki, Sakyo-ku, Kyoto 606, Japan and 2Department of Mechanical Engineering, Osaka Industrial University, Nakagaito, Daito-shi, Osaka 574, Japan (Received 4 April 1990)

A~tract--This paper provides an efficient numerical method on supercomputers for two-dimensional incompressible unsteady viscous flows. The unsteady vorticity~elocity Navier-Stokes equations are discretized by means of the implicit Euler formula and the resulting elliptic system of equations is solved by the parabolic multi-grid method with the group explicit relaxation scheme. Application of the present method to unsteady flowsin a square cavity shows the high computational efficiencyon vector and parallel supercomputers.

INTRODUCTION Recently, in the field of fluid mechanics, it has become possible to analyze practical engineering problems in three-dimensions, owing to the installation of supercomputers. The computational efficiency of a numerical method on supercomputers is very dependent on the algorithm and therefore should adapt to the architectural features of the vector and parallel computer. F r o m this point of view, Evans and Yousif[1] and the authors [2] have proposed a class of new iterative methods for elliptic partial differential equations, which are called group explicit iterative (GEI) methods. Their explicit nature makes them most suitable for vector and parallel computers. Application of a G E I method to the analysis of steady flows in a square cavity, has indicated that it is more efficient on a vector computer than other conventional iterative methods. On the other hand, as an approach to the efficient analysis of unsteady problems on supercomputers, the parabolic multi-grid (PMG) method has been provided by Hackbusch [3] for solving parabolic partial differential equations. This method is characterized by the simultaneous computation of several time levels in one step of the computational process. The computations in each time level can be independently performed by using several processors of a parallel computer. In this paper, a new P M G method has been developed for solving unsteady incompressible viscous flows in a square cavity, by using the G E I method as a relaxation scheme. The computed results, based on the Navier-Stokes equations in vorticity-velocity form, are compared with those by other methods or those based on vorticity-streamfunction equations in order to test the reliability and the computational efficiency of the present method.

GOVERNING EQUATIONS The two-dimensional vorticity-velocity Navier-Stokes equations can be written as follows: l ~, + (~u)x + (~V)y = ~ee (~.,x + ~yy + f ) ,

(1)

uxx + Uyy= -- (~ -- g,

(2)

vx~ + Vyy= ~ -- h,

(3)

where t is the time, x , y are Cartesian coordinates, ( is the vorticity, and u, v are the velocity components in the directions x, y, respectively. Re represents Reynolds number, and f, g, h are the residual functions which are assigned to zero on the finest grid. These equations in vorticityvelocity form are expensive as compared to those in the vorticity-stream function form. But 33

34

SHIGERUMURATAet al.

this formation has many advantages, e.g. the easy implementation of boundary conditions, the applicability to the three-dimensional problem, and so on. Furthermore, the most striking advantage is that the form of governing equations is independent of whether or not the frame of reference is inertial, as pointed out by Speziale [4]. So, with a view to further investigation on the present method, it is very worth while developing the computer code based on this approach and checking the computational efficiency of it. NUMERICAL

PROCEDURE

PMG method We will explain shortly the concept of the P M G method by taking eqn (1) as an example. The parabolic equation may be discretized by means of the implicit Euler formula as follows:

(k _At(k-~

1 (A(m)~k + A!7) (~u)k + A~?')(~v)k = Ree" xx +

A(xm)~k +f(")),

(4)

where At is the time step, and Ax, A>., A.,x and A,, are the difference operators for the derivative terms O/ax, a/#y, O2/Oxz and OZ/ay2, respectively. The superscript k denotes the time level, and the superscript (m) indicates that the functions and operators are defined on the grid level m for 1 ~
INITIAL

CONDITION ~t~ ~

A

I

:K I

ITERATIONNUMBER > °

K=2

", ~(CONVERGED) A; 1st PROCESSOR ~ ~ ~ . ~

~TIME

....

Z:26t ~ - - - - ~ O - ~ O ~ O - - -

Fig. 1. Computational procedure of the PMG method.

Parabolic multi-grid method for incompressibleviscous flows

k=O

k=l

k=K=2 time

3t = KAt FINr GRI

I

35

oDnin~GATION

i

m

,.

pmm+l

JR:-,.~ /

./ / ./ / //, 7,

/

RESTRICT] ON

i/

1/

I (m-1) COARSE GRID

Fig. 2. PMG iteration. approximation storage (FAS) mode of the multi-grid algorithms, proposed by Brandt [5], is employed. If eqn (4) may be represented by

r(=)

(U-,)(m)

(Lk(k)(') = R e +

A------~-= (Fk)(")'

(5)

where L is the elliptic operator, the restriction and prolongation of FAS mode are performed as follows: m~m

-

1

(from coarse to fine grid) ~(,,)<_((m) + p~_

l ( ~ ( m - l) __

R~-

l ~(m)),

(6)

m ~-m + 1 (from fine to coarse grid)

F (') <--L ('~)~(m)+ . ~. ,".

t)7(" +, ~ - +°

- L(" + 0~(m+ 0).

(7)

In eqns (6) and (7), the transfer operators R~_,, P,~+I denote 9-point restriction and prolongation, respectively. These operations, as shown in eqns (6) and (7), are executed according to the following criteria: m~m

- 1

(from coarse to fine grid) (Rs')(")<0.25(Rs") (m+°

m--*m + l

for

l~
(8)

(for fine to coarse grid) (9)

R s " / R s ' - ~ > 0.6,

where Rs" is the Lsresidual of the vorticity at the iteration step n, defined by

" -

(,,j

)/(N

-

1):.

(10)

~i=lj=l

The iteration is terminated when the following criterion is satisfied on the finest grid: (Rs") (M) < 1 x l0 -6.

(1 l)

Group explicit relaxation scheme

A GEI method is employed in the smoothing step of the P M G method. The GEI method was originally proposed as an iterative method for the elliptic problems, such as the steady CAF 19/I--D

SHIGERUMURATA et al.

36 (i)(i,3)

(2)(i+l,j)

(3)(i,j+l)

(4)(i+l,j+l)

I

I

,,(i,j+I)~ (i+l.j+l) (i,]) (i+l,j)

J

Fig. 3. Grid molecules in the G E l method.

Navier-Stokes equations, and the PMG method makes it applicable to the unsteady problems. Discretization of the spatial derivative terms by central differences, reduces eqn (2) for each grid point ( i , j ) to a difference equation. Figure 3 shows the computational molecules for the grid points ( i , j ) , (i + 1,j), ( i , j + 1) and (i + 1,j + 1). Then the system of four difference equations on these four grid points are expressed in matrix notation as follows: n+l

a

c

b

O

Uij

c

a

O

b

lti + I,j

di + I,j

"~- Ui W 2,j

- - b u i + l,j-- 1

b

0

a

c

Ui,j + 1

dij + l

+ Ui - I,j + L

- buij + 2

0

b

c

a

Ui+ I,j+ I

di+l,j+l

"~-//i + 2,3 + I

, (12) -bUi+l,j+2

where a = 1+flZ+p, dij = A x 2 ' gQ -

b = _f12,

c = -1,

(1 + f12 _ p )uij + 0.5. Ax • fl((,j+ l - ¢i~ =).

(13)

In eqn (13), Ax is the spatial grid size in x-direction, fl is the mesh ratio, p is the iteration parameter and superscript n denotes the iteration step. When an odd number of grid points is used in the computational domain, a certain number of ungrouped grid points appear adjacent to boundaries, as shown in Fig. 4. For these ungrouped grid points, the following equations are used: (14)

, a u "+1 " " i,j - ~ - - d i , j + U i + l , j

u"

"at- ~ - l . j - b u i j + l

"

-bunj

y

, -1

(15) •

(16)

The (4 x 4) and (2 x 2) matrices on LHS of eqns (12), (14) and (15) can be analytically inverted, hence the solutions at the next step (n + 1) are expressed explicitly. The present iterative method has four steps per iteration because there are four ways of grouping grid points, as shown in Fig. 4.

Parabolic multi-grid method for incompressibleviscous flows

37

7-J i

o

o----o

O---.-O

o--.-o

o----o

1st s~eD

O----o

o

o

O----O

O---O

O

o----..o

2nd step 31-d steo Fig. 4. Grouping of grid points for the case with N = 7.

qth step

Boundary and initial conditions The present method can be easily applied to the analysis of practical unsteady problems, such as the flow past a circular cylinder or an airfoil. But in these problems, the computation may be carried out under many different conditions on the size of computational region, the type of grid, the boundary condition and so on. Consequently, the computed results with different computational conditions may be inconsistent with each other. In this paper, in order not to mix the discussion of the computational efficiency and the reliability of the present method with that of the effects of the other computational conditions, the method has been applied to the unsteady incompressible viscous flow in a square cavity. The steady state flow is a popular and standard example for testing and comparing numerical methods, and it is very convenient for checking the computed results. The upper wall at y = 1 impulsively moves with unit velocity at t = 0. The initial conditions for all variables are given by the stationary state (=0,

u=0,

v=0

at t = 0 ,

(17)

and the boundary conditions for the velocity components u, v are given by

u(x,O)=u(O,y)=u(1,y)=O,

u(x,l)=l,v=O

for t > 0 .

(18)

Giannattasio and Napolitano [6] have compared four different discretizations of the boundary condition for the vorticity and reported that the Wood's and modified Wood's condition provide reasonably accurate results although the results do not satisfy mass conservation exactly. In this research, the boundary condition for the vorticity is imposed using Wood's condition. It may be written for the four boundaries, with Ay representing the grid size in the y-direction, as follows: ((x, 0) = --((x, Ay) - Ay "g(x, 0) + 2{u(x, 0) - u(x, Ay)}/Ay, ((x, 1)= - ( ( x , l - A y ) + A y

"g(x, 1 ) - 2 { u ( x , 1 ) - u ( x , 1 -Ay)}/Ay,

((0, y) = - ( ( A x , y) + Ax" h(0, y) - 2{v(0, y) - v(Ax, y)}/Ax, ((1, y) = --((1 -- Ax, y) - A x . h(1, y) + 2{v(1, y) + v(1 - Ax, y)}/Ax.

(19)

RESULTS AND DISCUSSIONS First of all, the results are shown for the case in which the total number of time levels per computational step is one, K = 1, in order to check the effect of the total number of grid levels M, and the iteration parameter p, on the computational efficiency. The Reynolds number is 100 and the number of grid points, N, in each direction x, y is 65. Figure 5 shows the correlation between the total number of grid levels, M, and the work units needed for the computation up to t = 21.5. A work unit is defined in such a way that the computational effort on the finest grid corresponds to 1 work unit. It is found from Fig. 5 that the computational efficiency for the case with M = 4, 5 is about 21 times of that of the single-grid method, M = 1. The convergence history is depicted for the first time step, at t = fit, in Fig. 6 which shows a plot of the log of the L2-residual vs. work units. Though the optimum value of the iteration parameter is shown to be 3, the computational efficiency is not sensitive to this parameter. Figure 7 shows the convergence history at the first three time levels, t = 6t, 26t, 36t, for the case with M = 4, p = 3 and K = 2. The constant

38

SHIGERU MURATAet al. 101

105

1 Re = 1 0 0 ,

N=

65

,K=I

, p = 3

3 10-1

z~

;

5

p

100

Re

104 3

10-2

0.5

t

4

1

K

n10-3

65

~J d 103

10. 4

t = 21.5

10 z

= 0.5

, ~t

10"5

I

I

I

I

I

2

3

4

I 5

I

10

T o t a l number of Levels

Work

Fig. 5. Improvement in computational work by the multi-grid method.

100

50 units

Fig. 6. Convergence history of the vorticity.

convergence rate can be obtained regardless of the mesh size, which characterizes multi-grid iteration. Next, the time evolution of the vorticity at the center of the cavity, x = 0.5 and y = 0.5, is shown in Fig. 8. The results of two different time steps, fit, are consistent with each other, hence it can be said that the unsteady flow in a square cavity is accurately obtained by the present method, with this time step. Table 1 shows the comparison of CPU time, required for the computation up to t = 1, on a vector computer the FACOM VP-200. It can be seen from this table that the computational efficiency of the present method for the case with N = 129 and M = 3 is about twice of that of the P M G method with Gauss-Seidel relaxation, because the latter relaxation scheme does not vectorize well due to vector interdependencies between neighbouring grid points. In addition, it is obvious that the computational efficiency for the case with K = 2 is almost equal to that of the conventional method, K = 1, in which eqn (4) is solved time step by time step. If a parallel computer is available, it is expected that the efficiency of the P M G method for the case with K > 1 will increase significantly, because the computation of K time levels can be simultaneously carried out by using K processors of the parallel computer. Table 2 shows the comparison of CPU time 101

100

;

N

=

17

0 ;

N

=

33

& ;

N

=

65

~t

=

1.0

M

=

4

101 I'-L R e =

IO0'N=

65'

M= 4'

P =3'

K = 2

1 0 "1

lO.Z -lo rr

10 .4 lO-Z

10 -s

10.. e

I

0

ZOO Work

units

Fig. 7. M u l t i - g r i d c o n v e r g e n c e .

400

10"3L

0

I

I

10

20

Time

Fig. 8. T i m e e v o l u t i o n o f ~ at x = y = 0.5.

39

P a r a b o l i c m u l t i - g r i d m e t h o d for i n c o m p r e s s i b l e viscous flows Table I. Comparison of CPU time between different relaxation schemes (N = 129 on the vector computer VP-200) Relaxation scheme

K= 1

Present Gauss-Seidel S.O.R.

K= 2

13.3(s) 14.0 21.4 24.7 Diverged

Table 2. Comparison of CPU time between different numbers of grid points (N ~ 65 on the parallel processors of INMOS T-800) Total number of grid points

Ratio

17 x 17 33 x 33 65 x 65

1.05 1.15

K= I 120(s) 437 1445

K~ 2

Ratio

74 267 868

0.62 0.61 0.60

on the two processors of an INMOS T-800. It can be seen from this table that the CPU time is reduced to about 60% by using this parallel computer, for any number of grid points. Furthermore, Table 3 shows the comparison of work units between different numbers of time levels in one computation, K. As K increases, the work units gradually increase. But the excess is about 15% for the case with K = 4, so it can be estimated that the computational time will be reduced to about 30% by using four parallel processors. These results are for the first order approximation in time, the implicit Euler formula, as shown in eqn (4). But more accurate time discretizations, which are required for the analysis of practical problems, can be easily implemented in the present method. Table 4 shows the comparison of work units, needed for the computation up to convergence, between the implicit Euler scheme and the Crank-Nicholson scheme. The latter is accurate to second order in time. The excess of work units of the Crank-Nicholson scheme is about 44% for the case with K = 4, and it is about three times as large as that of the implicit Euler scheme. In the Crank-Nicholson scheme, the iterative solution (ik7 j are obtained by using the solutions at five grid points (i,j), (i + 1,j), (i - 1,j), (i,j + 1) and (i,j - 1) at time step k. On the other hand, in the implicit Euler scheme, the solutions at only one grid point (i,j) are needed for the solution (k+ ~. Because the solutions at time step k have not yet converged for k > 1, the iteration may be more wastefully performed in the Crank-Nicholson scheme than in the implicit Euler scheme. However, the computational time will be reduced to about 36% by using four parallel processors, as shown in Table 4, and therefore it can be said that good convergence is also achieved for higher order implicit scheme. The results for the case with Re = 1000 and N = 257, are shown to check the reliability of the computed results by the present method. Figure 9 shows the constant vorticity lines and the velocity vectors in a square cavity in order to check the flow pattern. The vectors indicate only the flow direction; the length has no meaning. One secondary vortex is visible in each bottom corner of the cavity. Figure 10 shows the velocity profiles for u along the vertical line and v along the horizontal line passing through the geometric center of the cavity, and Fig. 11 shows the vorticity distribution on the moving wall. In these figures, the results by Ghia et al. [7] are included, which are denoted by the squares. The results obtained by the present method are in good agreement with their results which are based on the Navier-Stokes equations in vorticity streamfunction formulation, completely satisfying the equation of continuity. Finally, the contours of the residuals for the equations of continuity and vorticity definition for Re = 1000, N = 129, are shown in Figs 12 and 13, respectively, where the residuals are defined as follows: D,,=

8u

(20)

Ov i,j ,

Ov 8u ~,J" A~,.j = ~ - ~xx + ~

(21)

In both figures, it is found that large residuals concentrate at the corners on the moving wall, where the discontinuity of velocity occurs and the vorticity becomes infinite. So the residuals should decrease for problems with continuous boundary conditions. Table 3. Comparison of work units between different numbers of time levels K (N = 65 and At = 0.25.)

Table 4. Comparison of work units between different implicit schemes (N = 65 and At = 0.25.) Implicit

K

Work units

Ratio

Ratio/K

scheme

I 2 3 4

1454 1515 1627 1677

I 1.04 I.I 2 I. 15

I 0.52 0.37 0.29

Implicit Euler CrankNicholson

K= I

K= 4

Ratio

Ratio/4

1454

1677

1.15

0.29

3624

5240

1.44

0.36

I = `f II e~x ~UTAOtU 0 q l UO u o ! m q ! n s ! p ×

,~l!3!l.IOA "I ] "8!zl

~ ' 0 = x SUOle soItJO.ld ,~1!3o1~A "01 "i~!d

'g'0 = a" p u e

f3 ,~.! 3OlaA O'L

g'o

I

I

0

!,-

g'O = X 5u°1° ' / X _

^+~'~-^

;

:

~

£

=

d

Zg;~ 000¢

< o ,-f OCj-

L= x

oj,; o,o,o/

: N : e~

o

:

"I

S-

n ,<

I u~ot ~ - m )

/

;I.UOSe~d : 00~ -

•£$OlOUq3oI jo olnl!lsu I olo£N 3o 0 0 8 - i SOIAINI aolndsueal JO sJossooo.Id Plle.~ed o~1 oql pue 'nsl!fna .to opera ',~].ISTOA.IH["l 0"10£'~I JO .IDlU0~) .mmdtuo 3 ~ql Jo ZSf-IA[ .t~lndtuo3 .~ele3S zql '00~-dA IAIODVH azlnduao3 aolo~h zql zsn ~m '.reded s!ql u I •[leA~ i~maoua ~ql ol lug3~.l'p~ saouao3 ~ql le ~le.tlu~3uo3 uo!l!uij~p Xl!3!laOh 0ql pue ,~l!nu!luo3 jo suo!lenbz ~ql .Ioj suo!lnlos le3!.mtunu j o slenp!sz.t ~q.L •t u z o j uo!13unjtue&~ls-z~l!3!l.~oh ~ql uo pzseq ~soql ql!m luotuooaSe p o o ~ ut oae tuaoj .'~I.130[0A--/~I.1D.II.IOA0ql uo poseq sllnsoa polnduaoo oq.L •aolndtuo3 [Otleaed oql uo aoq$!q somoooq ,(ouo!3~o ieuo!lmndtuoo oql 'soseo33ul SlOaOI otu!l jo .toqumu ielol oql se 'o3uo H ".mlndmo3 CeleOS e uo ~l.tom Ieuo!lelndmo3 oql uo 13o~ oIll!l seq oanpoooad i'euo!lelndtuo3 ouo u! S[OA01 Ottt!l JO aoqtunu i'elol o q i •aolndtuo3 .IO100A 13 UO 'Otuoqos uo!lexeioa loP!oS-ssneo ql!m poqlom p!~$-!llntu jo leql ueql aoq~!q s! poqlotu luoso.~d oql jo /~31./O.t31Jjo leUO!lemdtuo3 o q z

(17) (f)

(E) (I)

:saxolioJ se poz!.~mutuns oq Aetu sllnso.t oqj. "/~1.IA133o.ll3rlbs 13 HI A~o[j gpeolsun oq:t aoj podolo^Op uooq seq otuoqos uo!leXelOa l!3!idxo dno.t$ e ql!m poqlotu p!a$-!llntu 3!ioqeaed oql 'qo.mosoa s!ql u I SNOISflqDNOD 'SJOlaOA £1!aOla A (q) " u o ! l n q ! J l s ! p £1!3!1ao A (~) "000I = ~ 1 ~oJ u o ! m l O S al~ls £P~OlS "6 "$!rl × '[

×

~i0

0'0

,! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..........................

~. . . . . . . . . . . . . . . . . . .

It

| , ,,,;,,h|h,, .... ~ " " H # II t t t l tt~t--~t

I

II'f'lht .... I IIIIP

.......... ~ .... ~. . . . . . . . . ~

I I

,

r

,

~'0,

J

I

,

r

',00,0

dJI h,,,llll

II

......................... : ................ lhl II. . ..... .... . . . . . .I.l.l. . . . . :~ . .". .. .. ....... .. .. .. . I i [,,1111111,

0

u ~

0

$'0

'<

~1

It

IL_

.................. I I H ............ ,,,,,I~ , | II" , H , ;~,$ IIIIIII ""'" ++ t,,,, ~ ................... ~, # I |111

,,,,,,,,, , ......... - ............... IIIIII " " . . . . . . . . . . . . . . . . . . . . . . . . . .

...................

llllllt%~..

.

.

.

.

.

.

.

+mttSI+ItlIIIIPIIII

Illh, ........ . . . . . . . . . . . . . . . . H,! ............ '. . . . . . . . . . . . . . . . . Ih,, . . . . -....-... . . . . , , , , , , , , , . . . . . . . . . . . . . . . . . . . . . . I. . . . . . . . . . . ,hit

l

, ~ fillll I flirt,Ill-" . . . . . . II H

(q) "/D la YI¥~ICI~ CI~I~gIH S

(o) 0~'

Parabolic multi-grid method for incompressible viscous flows 1,0

41

'

,

o,os

20,

o,oo

0

0'%,0

I

L

I

l

1

0,5 x

I

l

I

I

i,

Fig. 12. Residual contours for the equation of continuity.

0,0

I

O,

I

I

I

I

0,5 x

[

1

I

I

1,0

Fig. 13. Residual contours for the defining equation of vorticity.

REFERENCES I. D. J. Evans and W. S. Yousif, The alternating group explicit method (BLAGE) for the solution of elliptic difference equations. Int. J. comput. Math. 22, 177 (1987). 2. S. Murata, N. Satofuka and T. Kushiyama, Numerical analysis for square cavity flow using a new Poisson solver. In Computational Fluid Dynamics (Edited by G. de Vahl Davis and C. Fletcher), p. 547. North-Holland, Amsterdam (1988). 3. W. Hackbusch, Parabolic multi-grid method. In Computing Methods in Applied Sciences and Engineering VI (Edited by R. Glowinski and J. L. Lions), p. 189. North-Holland, Amsterdam (1984). 4. C. G. Speziale, On the advantages of the vorticity-velocity formulation of the equations of fluid dynamics. J. Comput. Phys. 73, 476 (1987). 5. A. Brandt, Multi-level adaptive solutions to boundary-value problems. Math. Comput. 31, 333 (1977). 6. P. Giannattasio and M. Napolitano, Numerical solutions to the Navier-Stokes equations in vorticity velocity form. Proc. Third Italian Meeting of Computational Mechanics, Palermo, p. 95 (1988). 7. U. Ghia, K. N. Ghia and C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J. Comput. Phys. 48, 387 (1982).