The calculation of some laminar flows using various discretisation schemes

COMPUTER

METHODS

NORTH-HOLLAND

IN APPLIED

PUBLISHING

IMECHANICS

AND

ENGINEERING

35 (1982)

293-313

COMPANY

THE CALCULATION

OF SOME LAMINAR FLOWS DISCRETISATION SCHEMES

USING VARIOUS

Andrew POLLARD Mechanical

Engineering

Department,

Queen’s

University,

Kingston,

Ontario,

Canada

Alan L.-W. SIU Cowlputer Modelling

Group, Calgary. Alberta,

Received Revised

For various schemes

formulation

schemes

Two

new

the upwind

be

stable,

obtained install

it can, the

one

the solutions

using

validation;

into existing

form obtained

the upwind

further

schemes

and/or

14 May

the upwind,

hybrid

schemes

generate

grid distribution and data from

seems using

to

be

these

the hybrid

1982

and the quadratic

are used in their

both

both

formulation,

found

negative

of

which sources.

sensitive

to

normal

and positive

the

generate

These size

of

the

are equal that

are provided

these

two

coefficients.

positive

influence

to the results

two new schemes

It is concluded

new formulations

The quadratic than

influence

by reference

new formulations

finite-difference

forms.

to be less stable

always

and tested

other

schemes.

and to this end, these

computer

1982

received

has in the past been

on occasion,

for a uniform

and hybrid

although

flows,

and hybrid

scheme)

quadratic

are presented

Furthermore, require

of

upwind

QUICK

because

forms

coefficients, to

two-dimensional The

(Leonard’s

other

from

laminar,

are evaluated.

4 January

manuscript

Canada

obtained are found

finite-difference

to or better these

than

grid. those

new formulations

in a form that is easy to

codes.

0. Nomenclature A?

area for 4-equation in i-direction

L M:,

M;

(i = e, w, n, s); a Bf’

i-spacing for plates; finite-difference coefficients of @-variable in i-direction

N Pe P

(i = E, W, N, S, P); Cf d, D Df 1

I,

convection coefficient of 4 at i-interface; diameters; diffusion coefficient of 4 at i-interface; length; re-attachment length;

004S-7825/82/0000-0000/$02.75

@ 1982 North-Holland

Re %I, S” T

u

length; sign determinator

(see

(2.11)); number of iterations; Peclet number; pressure; radial coordinate direction ; Reynolds number; source terms in U, V momentum equations computational time (CPU sets); axial velocity component ;

A. Pollard, A.L.- W. Siu, Various discretisations for laminar jlows

294

V X

y

transverse velocity component; axial coordinate direction; transverse coordinate direction;

Greek

4E,

&‘, dh &G 6P

N, S, P; * 5

p

density; dynamic viscosity; scalar variable; exchange coefficient 4;

e, w, n, s

of

stream function; vorticity;

Subscripts 4

,S 4 r,

values of scalar variable at grid points E, W,

any variable; compass point directions with respect to central point p, also grid interface locations.

1. Introduction In the realm of computational fluid mechanics, false diffusion is a major obstacle to the achievement of accurate and economical solutions to the governing equations; see, for example, [l] _ False-diffusion arose when practitioners of computational fluid mechanics replaced central-differencing (a scheme that achieves stable, non-false diffusion affected results through the use of fine grid distributions, such that the grid Peclet number, the ratio of convection to diffusion across the grid, is less than 2) with one-sided differencing schemes (which effects stable solutions on coarse grids). The introduction of these schemes (commonly referred to as ‘upwind’) afforded economical solutions but reduced the accuracy of the solutions because they introduced a truncation error; this error is commonly referred to as false or numerical diffusion. Recently, Leonard [2-61 introduced a scheme that tends to eliminate false-diffusion while still maintaining the economical assets of one-sided differencing schemes. Leonard’s scheme (known by the acronym QUICK: Quadratic Upstream Interpolation for Convection Kinematics) proposes that, instead of using linear interpolation for the convection terms as used in standard one-sided differencing schemes, a three-point upstream weighted quadratic interpolation should be used. Indeed, it has been shown that this approach provides solutions to the governing equations that are more accurate than those achieved by other one-sided differencing schemes (i.e., upwind-differencing) (7-101. An attribute of Leonard’s scheme is its apparent simplicity and ease of use. However, as noted by previous researchers [ll, 121, this attribute is open to question. Han et al. [ll] attempted a direct application of Leonard’s QUICK scheme, apparently in a manner similar to [8,9] (in [8,9] details of how the scheme was applied are not given). As a result of this direct application, stable solutions to the equations could not be secured unless a pseudo-source term was introduced; a feature not noted by Leonard [2-61 or in [8,9]. The differencing formulae were then re-formulated and stable solutions could then be obtained. However, when extended to high Reynolds number flows (i.e., turbulent), even these (reformulated) formulae required the introduction of the pseudo-source term to secure stable solutions. The need for a pseudo-source term to secure stable solutions is, to these authors, most distressing.

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

295

The purpose of this paper is to present a redevelopment of the QUICK scheme that effects stable solutions without the need for extra terms, thereby providing a common thread from which further advancements in this area can be achieved. The redeveloped QUICK scheme is tested here by comparison against other schemes (upwind and hybrid-differencing) and data from other sources for three laminar flows of varying complexity. The rest of the paper is now outlined. The next section provides the equations in a general form that govern laminar flow. These equations are then discretised and the approximations made to model the convective terms are provided for the upwind, hybrid and QUICK schemes. Two new formulations of this latter scheme are given; these are called QUICKE and QUICKER. The four schemes are then applied to three laminar flows in Section 4. The paper closes by summarizing the results and provides conclusions and recommendations for future work. 2. The mathematical

formulation

2.1. The equations The equations governing the following forms. The continuity equation:

steady,

incompressible

laminar

flow can be conveniently

cast into

(2.1) The momentum x-direction:

or Navier-Stokes

equations:

(2.2) r-direction

:

(2.3)

where U, V are the velocities in the axial (x) and radial (r) coordinate directions, p is the density, p is the pressure, p is the dynamic viscosity, S, and S, are the source terms in the 17 and V momentum equation respectively. Here, S, equals zero, and S, contains a term generated from transforming the equations from the Cartesian to the cylindrical coordinate system. Indeed, the above equations describe motions in Cartesian coordinates when & and r are set equivalent to dy and one, respectively. For future reference, the above equations can be cast into the following form (2.4)

296

A. Poftard, A.L.- W. Siu, Various discretisations for laminar flows

where C$ denotes any scalar variable. When &, and S, equal zero, and # equals one, (2.4) becomes the continuity equation. When # stands for U (or V>, and &, stands for ,u3 the corresponding momentum equation in the x (or r, or y) direction is obtained. S, now contains the pressure gradient term. 2.2. The finite diflerence discretisation Equations (2.4) are solved, with their appropriate boundary conditions, by first integrating them over finite-difference control volumes that overlay the physical domain considered; an example of the grid structure is shown in Fig. 1 where it is noticed that the grid is staggered so that velocities are located mid-way between the grid points. The centrally located grid points serve as locations for all other general scalar variables (e.g., pressure and viscosity). Now, consider a single control volume for 4, as shown in Fig. 2; integration of (2.4) over this control volume yields

where the A’s represent the areas of the cell faces in four-compass-point directions (12,e, s, w) located mid-way between the grid nodes. The next step in formulating a finite difference equation is the assumption of the 4 profiles between any two grid points. The diffusion terms are formulated using the central difference approximation; and, since this is common practice, further attention will not be placed upon them. Attention is directed to the convection terms (e.g., the rpU+ terms), for it is the

L NUMBER CONTROL VOLUME TYPE

GRID LINES

\-GRIO 4

I

POINTS 5 U

BOlJiOARY

Fig. 1. Control

AND COt$JU-

volume

BOUtkIARY

specification.

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

N

p-

‘I Ax-,

\w

,.

r--------

A

w

P

--,e

4.1

L :

rv ru,

r

. T

Ar

..E ..

I

L_____s__ _.A 8X,

.,S

..

sxe

A

x

Fig. 2. Control

297

volume

--c

. A

for a scalar variable.

approximations that are used for these terms that can lead to the generation of false-diffusion. Furthermore, in what follows, it is considered that the schemes that are used to approximate the convective terms are only applied to the convected variable (i.e., 4 in (2.5)); the convecting velocity is discretised using the central differencing scheme. Also, for convenience, only uniformly distributed control volumes are used for discretisation purposes; the extension to non-uniform grids, while in some instances is not straightforward, can be, for the present, dispensed with. (a) The Central Difference Scheme (CDS). In this scheme, the value of 4 existing at the control volume interface is taken as the average value of the grid points that lie on either side of the control volume interface. For velocities that are low enough, this scheme is recommended; however, it has been found (see, for example, [13]) that when the grid Peclet number (LJAxlr) is greater than 2, the coefficient matrix becomes non-diagonally dominant. As a result, a semi-implicit-type numerical scheme, as is normally used, becomes unstable. (b) The Upwind-DifSerencing Scheme (UDS). The stability problem associated with the use of the CDS can be overcome by the introduction of one-way influences: in evaluating the value of 4 at an interface, the 4 value takes into account the direction of the local velocity vector located at the interface. In other words, we have

(2-6) Unfortunately, by instituting this approximation, numerical diffusion is also introduced (the method by which this is introduced into the finite difference equations can be found in, for example, [3]); and, at the same time, it is presumed that regardless of how large the velocities are, diffusion effects are still present. Furthermore, at very low velocities, such that Pe < 2, the upwind formulation is used, even though CDS would be preferable. (c) The Hybrid-Diffrencing Scheme (HDS). In 1972, Spalding [14] introduced the HDS which permitted the +-value at the control volume interface to be evaluated by letting the CDS prevail when Pe < 2 and UDS for Pe 2 2; and for the latter case (i.e., high velocities) the diffusion term in the general equation is assumed ineffectual (i.e., zero). This latter assumption serves the role of allowing, in part, the false or numerical diffusion to act in place of the now neglected real diffusion.

298

A. Pollard, A.L.- W. Siu, Various discretisati~~s for laminar flows

(d) The Quadratic U~~trea~z Interrogation for convective Ki~ze~at~~.~ (QUK’K) Scheme. In an effort -to the aforementioned generating schemes UDS and Leonard [2-61 and widely his QUICK This scheme a three-point weighted quadratic for each As a when inserted (2.5) his provides a difference equation is third accurate (i.e., diffusion, a order truncation is absent). reference to 3, the formulae in dimension are:

(2.7)

as shown

formulae can Fig. 4.

extended

to

dimensions

by

a 9-point

star configuration

ww

Fig. 3. QUICK

interpolation for interfacial values.

Fig. 4. Nine-point star configuration scheme.

EE

for the QUICK

2.3. The finite difierence equations

inserting

the

formulations

for the

(2.8) where the 61’s are the coefficients made up of contributions from diffusion and convection, the latter being obtained with the aid of the differencing scheme chosen (i.e., CDS, UDS, HDS or QUICK). The B$ t erm is the sum of all the other B’s, and the S’s are the source terms containing all terms resulting from transforming the partial differential equations into their finite difference counterparts and not explicitly shown in (2.8). Obviously, the B’s and S’s are uniquely formulated for each differencing scheme. It is to these formulations that our attention is now turned. Before doing so, however, note that the CDS is not presented here for, because of the extremely fine grid distribution needed to ensure that locally Pe < 2, it finds little direct use in computational fluid mechanics.

A. Pollard,A.L.- W. Siu, Various discretisations for laminarflows

2%

(a) The Uphill ~i~ere~ci~g Scheme. The derivation of the B’s and S’s for this scheme (and the HDS that follows) will not be given in detail: they can be found in, for example, 1151; here the results of manipulating the equations ate given. With reference once again to Fig. 2, the forms of the B’s and S’s are

Also S: = r,Ar(Pp - PE), L.PAx r,

s;: = i0

S1: = rtiAx(Pp - PN) I Sg = 0,

for cylindrical-polar

coordinates

only,

otherwise .

The brackets Q1 denote that the B’s are assigned the maximum of the values contained within them. It is seen that the B coefficients and S terms are always positive, thus making the fractional coefficients in (2.8) positive. This type of finite-difference approximation is known to be unconditionally stable; see, for example, [ 16). (b) The ~y~r~~-~i~ere~c~~g Scheme. The source terms for the HDS are the same as those for UDS. The B’s have the following different definitions, even though they share the same convective and diffusive coefficients as before,

(2.10)

where, once again, it is seen that ail the B’s and S’s are positive, implying therefore an unconditionally stable system of equations. (c) The QUrCK Scheme. Using the same definitions for the convective and diffusive coefficients, as noted in (2.9), the QUICK scheme gives different sets of B coefficients and source terms. The first attempt at formulating the B’s and S’s utilizing the QUICK scheme was not concerned with formulating the coefficients into a neat package as indicated for UDS/HDS. The B’s were formulated so as to provide insight into the effects of convection and diffusion, and the ‘curvature terms’ were inserted into the S’s. The ‘curvature terms’ are defined here as

300

A. Pollard, A.L.- W. Siu, Various discretisations for iaminar flows

those terms outside the standard convection and diffusion formulation normally used in UDS/HDS. This practice appears to be common among those who have used QUICK [9, I I J; the word ‘appears’ is used intentionally since in these references, details as to the actual distribution of all the terms {in both the coefficients and the source terms) are not given. As an example of this first attempt, the coefficients and the source terms for the Umomentum equation are given below. They are

where

and

An examination of the B coefhcients, shown in (2.11), indicates that, if the convection effects are strong enough, they can assume negative values. As a result, stable solutions cannot be guaranteed. A similar situation could occur if the denominator (B$ - St) in (2.8) becomes negative. The above formulations were not inserted into a computer code because of the above mentioned limitations on the formutation of the coeflicients and source terms. However, it is interesting to note that Han et al. [l l] apparently instituted (2.1 I) into their code but convergence could not be obtained unless a ‘pseudo-source’ term, based seemingly upon intuition, was introduced; this term, in effect, damped the system of equations that permitted solutions to be had only after a considerable expenditure of computer time. In an effort to eliminate the need for this ‘pseudo-source’ term, these authors re-formulated the interfacial term was more influential; for exampie, if distributions of (i.e., (2.7)) so that the convection U, > 0, sbp = 0.125(6& + 4&) - 0.125(& + &). However, in doing so, the latter expression had to be evaluated using +-value calculated (and stored) from the previous iteration cycle, otherwise “convergence was destroyed or greatly impaired”. Furthermore, when calculating turbulent flows, these authors found it necessary to introduce, once again, their ‘pseudosource’ term to achieve converged solutions. Even so, for future users of their formulation, no guidance is provided as to what should be the magnitude of this source term. To achieve an ordered approach to the formulation of the coefficients, the finite difference equations were re-arranged in such a way that the fractional coefficients are always positive. Rest&s of such effort, which also provides the coefficients in a neat package for direct

301

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

insertion

into a finite difference

code, are now presented.

They are

;c:>,(D$ - $C$- AC”,),(DZ - $CZ),(Df + xl? - $C”,)], B$/ = I[(D”w- :ct>,(D$ + $c”w+ &?), (D”w+ tct), (lx -SC”, + fct)n, Bfi = up: + :c:>,(0” - :cP: - C), (0: - $C$),(Ix + 32 - ic>g, I?$ = [(of - $C$),(Df + :c: +kc:), (D$ + gc:>, (lx - :c: +ic)~,

B-g = [(D$ +

I?$=l?~+l?~+BgJ+B~, + M:C:

S,” = rUAr(Pp - PE) - i{M:C,U U,, -~{M:c,UU~+ S,“= rVAx(P,

~;c,Uu~

U,, - M,C,“U,,

-~;c:u~ + M;C:V,,

-P,)-${MX,VVw

- M,C:

U.N}

-~;cYu~j,

(2.12)

- M;C:VEE

- M,C,VV,,}

-~{M:c~~+M;:C,UV,-M;~,U~~-M;CYV~}, S: = ${M:C;

+ M:C,U - M,C:

- M;C,U}

+;{M;C:+M;C,U-M,C:-M,C,LI}, S~~=d{M:CY+M~C,V-M,C,V-M,C,V}+~{MdCY+M:c,V-M,cY-M,c,V), Cartesian S&=

coordinate,

- Gz-$+,&

for cylindrical

coordinate

.

set of B coefficients gives positive values regardless of the magnitude and direction of the convective effect. Equations (2.12) constitute and what will be referred to as, the QUICK Extended form, or QUICKE. Note that the source terms S$’ in (2.12) can assume positive values, and may render a negative value to the denominator of (2.8). To avoid such happenings, the source terms are revised as follows. This

new

S$* = S$ + St4*,

and

St’ = 0

where +*p is the current (‘in-store’) value of &. This re-adjustment to the QUICKE formulation will be hereafter referred to as the QUICK Extended and Revised form, or QUICKER. Note that as a result of this revision, the semi-implicitness of (2.8) is reduced. Indeed, by virtue of this, we can expect that the QUICKE formulation should provide solutions in time periods shorter than the corresponding times for QUICKER; however, QUICKE may not be as stable as QUICKER. Note that details regarding the derivation of (2.12) can be found in [17].

3. The solution procedure The finite-difference equations were solved Spalding [18]; this algorithm was incorporated

using the SIMPLE into the general

algorithm of Patankar and two-dimensional numerical

302

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

code, called 2/E/FIX, of Pun and Spalding [19]. In this code the primitive variables are computed in a semi-implicit line-by-line manner over a staggered finite-difference grid (see Fig. 1). Because of the semi-implicitness of the code, under-relaxation parameters are employed. Unless stated to the contrary, these parameters were set to 0.5 for all variables and for all flow situations and for all schemes. Since concern is primarily directed towards the outcome of the computations using this code and not the code itself, readers interested in the details of the 2/E/FIX code can refer to [19] or, for example, [20].

4. The results of the calculations 4.1. The scope In this section of the paper, the results of some calculations are presented. Consideration is given to three laminar flows: (i) flow between parallel plates; (ii) axisymmetric flow in tubular sudden expansions; and (iii) flow in a square cavity with a moving lid. It is noticed that advancing in numerical order also increases the flow complexity; and, while it is acknowledged here that case (iii) is the only flow situation that will exhibit strongly the effects of false diffusion, the two other flow situations are used to gain insight into the relative performance of the upwind, hybrid, QUICKE and QUICKER schemes. Furthermore, laminar flows are chosen so as to establish that the present formulations of QUICK (i.e., QUICKE, QUICKER) are indeed stable and accurate before additional complexities, such as models of turbulence, are introduced. The three physical situations considered are depicted in Fig. 5, together with the appropriate physical quantities characterizing each flow situation. In what follows, each flow

(a)

‘TOP

Fig. 5. Physical

situations

considered.

A. Pollard, AL-W.

Siu, Various discretisations for laminar flows

303

situation will be examined separately. Within each flow situation examined, the UDS, HDS, QUICKE and QUICKER will be evaluated in terms of accuracy and computational expense; this latter term implying the number of iterations required by the computer code to reach a pre-determined level of convergence (i.e., summing all the terms in each equation over all computational control volumes and normalizing with respect to the influx of the quantity being considered). 4.2. Flow between parallel plates

42.1.

The physical situation considered

The physical situation is shown in Fig. 5(a). Fluid enters the channel with a uniform velocity profile. Downstream of the point where the boundary layers intersect, the flow becomes fully developed. There is a lack of good experimental data for this flow situation. However, a number of publications dealing with the numerical solutions to the complete Navier-Stokes equations can be found in the literature, see, for example, 1211. The work by McDonald et al. 1221, who used the stream-function vorticity forms of the equations together with non-false diffusion generating central-differencing, was selected here as the basis for comparison. Note that McDonald et al. used a uniform mesh of 21 x 201 to cover a plate length of 20 channel widths. 4.2.2. Some computational details Since the flow is symmet~c about the mid-plane of the channel, computations were performed only in one-half the channel width. Convergence criteria of lop3 for each equation was found sufficient to ensure that a solution to the governing equations had been obtained. Calculations were performed at a Reynolds number (Re = pU(2a)/p) of 150 since at this Re the information provided by McDonald et al. is quite complete. The boundary conditions for this problem were supposed as follows: (i) inlet: specified uniform velocity profile; (ii) outlet: au/dx = 0 and V = 0; (iii) wall boundary: U = 0 and V = 0; (iv) axis of symmetry: ~~/~y = 0, V = 0. 4.2.3. Comparison of results Fig. 6 depicts the calculated centre-line velocities as a function of the distance downstream of the plate entrance for Re = 150. The curve in the upper right portion was calculated using plate lengths of L/2a = 11.5. However, because of the uniformly spaced 20 x 48 grid (the maximum used) was found to be too coarse for accurate results in the near entrance region, the calculations were repeated using a domain one quarter of the previous channel length (i.e., L/Za = 2.875) while still using a 20 x 48 grid. The computed results have been plotted in the lower left portion of Fig. 6. It can be seen that the calculated results of the UDS, HDS and QUICKE/ER schemes agree very well with the results of McDonald et al. For the case of the shorter plate length, the exit boundary condition is quite probably doubtful; however, checks of the development of the centre-line velocity and pressure gradients indicated no discontinuities in their slopes in this region. It is believed, therefore, that the boundary condition imposed has little effect on the calculations.

A. Pollard,

304

A.L.-

W. Siu, Variow

discretisations

IO

_____“DS --HRS

r

for lamirlar j?ows

---__ -W

.75 x/(Za)

= 2.875

li

@f-t 1 2r

a

.50 1 ----(JDS

-HDS - - - QUICKE/ER 0 MCDONALD

-

.2q-

I I.0

--

10-I

I.0

1 +-PI

I’

I

$=

/E!f

A‘

I

Ir

eta1

::

,

5,

1.0

1.2

b

,

i

-.-I

I

IO

IO2

/n Fig. 6. Laminar Bow between parallel plates: ment of centre-line velocity at Re = 150.

0.0 0.0

I

/

I

I

0.2

0.4

0.6

0.8

U( yl/U

Y

develop-

Fig. 7. Laminar flow between profiles at x/a = 1.0.

parallel

plates:

velocity

In Fig. 7 the velocity profiles near the entrance (at x/a = 1.0) are displayed; these results are from using the shorter plate lengths. The velocity profiles are characterized by overshoots: the physical reason for which can be found in, for example, [21]. With the present 20 x 48 grid, it is evident that the QUICKE/ER forms provide slightly better predictions of the velocity profiles than do the UDS and the HDS. It is interesting to note that near the wall, the UDSIHDS provide slightly smaller velocity gradients than do QUICKE/ER, although the curves in this region seem coincident. A smaller gradient here is characteristic of flows at lower Reynolds numbers (see [23]). This discrepancy seems to be caused by the appearance of false diffusion generated by the UDS/HDS, since the velocity vectors in this region are directed at a slight angle relative to the grid lines. The velocity profiles at x/a = 4.0 are shown in Fig. 8; again, these profiles result from using the shorter plate length. Here, the three curves almost lie on top of one another. It is obvious that the flow is ~1~7~~~fully developed at this location. Hence, the velocity vectors are inclined at much smaller angles to the grid lines. It would appear, then, that the false diffusion generated in the UDS/HDS becomes so small that it no longer affects the solution in this region.

.O -

QUICKE/ER MCOONALD

‘! etol

4 I

0.0 0.0

Fig. 8 Laminar

I

I

25

.50

flow between

parallel

I

.75 UlyVU

plates:

A

I

I

I .oo

1.25

velocity

profiles

I 1.50

at x/a = 4.0.

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

305

It is necessary to emphasize that the solutions obtained by QUICKE and QUICKER are exactly the same. However, as seen from Table 1, the revised form, QUICKER, takes about 2.5 times the number of iterations, and therefore 2.5 times the computational time required by the extended form QUICKE. The reason being that the QUICKER form is less implicit than the QUICKE form in the treatment of source-like terms. Furthermore, Table 1 shows that QUICKE takes about 50% more computational time per node when compared to the UDS/HDS. However, the better interpolating accuracy acquired through the use of the QUICK scheme requires fewer iterations to reach a converged solution. Consequently, in terms of overall computational time, the QUICKE form shows savings of 23% and 13% respectively over the UDS and the HDS. Tabte 1 Comparison of convergence 20 x 48 grid

$

= 2.875

& = 11.5

rates for plane flow at Re = 150,

Scheme

Tl Tugs

N

TIN”

UDS QUICKE HDS

1 0.771 0.889

139 121 73

2.49 3.65 2.54

QUICKER

1.969

191

3.56

UDS

1

1.58

2.54

QUICKE HDS QUICKER

0.853 1.006 1.868

158 9.5 209

2.55 3.60 3.585

“T/N = Actual computer time (secs)/number of iterations.

4.3. Axisymmetric

fiow through a pipe sudden-expansion

4.3.1. physical situation Fig. 5(b) depicts the physical and flow geometry of the pipe sudden-expansion flow. The diameter ratio (d/D) used in this situation was fixed at 0.5. As seen from the diagram, the flow is assumed symmetric about the centre-line axis. Just behind the step, and bounded by the larger diameter pipe, there is a recirculating eddy where the velocity vectors can orientate themselves at very large angles to the grid lines. In the core region, initially bounded by the inner surface of the eddy, the fluid gradually adjusts itself until it reaches the fully developed condition further downstream, i.e., at about 3 times the eddy Iength away from the inlet. Some computational details length of the calculation domain in the stream-wise direction was taken as four times the reattachment length which was found to be of sufficient length such that the exit boundary conditions (as noted below) had no effect upon the details of the flow. Convergence criteria were again chosen as IO-” for each equation. The following boundary conditions were incorporated into the computer program: (i) inlet: specified parabolic velocity profile; 4.3.2.

ne

A. Pollard, A.L.- W. Siu, Various discretisations for laminar j?~ws

306

(ii) outlet: au/&~ = 0 and V= 0; (iii) wall boundary: U = 0, V = 0; (iv) axis of symmetry: au/Jr = 0, V= 0. 4.3.3. Presentation and discussion of results In this flow situation, Reynolds numbers of 50, 100, 150 and 200 were employed so as to determine the ability of the UDS, HDS, QUICKE and QUICKER schemes to provide physical realism and recirculation zone lengths in accord with those obtained by Macagno and Hung [24]. A variety of grid densities were used to calculate these flows; however, it was found that a 16 x 25 grid, distributed in the radial and axial coordinate directions respectfully was sufficient to obtain essentially grid independent results. Note that Pollard [25] found that a 28 x 57 grid was needed to achieve fully grid independent results using the HDS, however, the present results differ from [25] only by 2%. The re-attachment lengths computed via the use of the four discretization schemes are listed in Table 2 together with some pertinent computational details. Table 2 Comparison axisymmetric numbers

of computational details and re-attachment lengths for Iaminar flow in sudden expansions, using various Reynold’s

Reynold’s number

Scheme 16 x 25 grid UDS WTLJDS iv

Md

HDS TITUDS N L/d T/Tuos N Ud QUICKER T/Tuns N Ud 1,/dof Macagno

QUICKE

and Hung [24]

50

100

150

113 2.15 0.96 107 2.26 1.22 93 2.20 l.Y8 151 2.20

1.18 91 4.32 2.02 156 4.32

1 108 6.26 0.96 106 6.61 1.16 89 6.45 2.06 159 6.45

2.2

4.3

6.5

1

111 4.20 0.96 107 4.43

200 1 111 8.34 0.84 100 8.82 1.19 91 8.6 2.06 157 8.6

8.8

From Table 2 it is noticed that all the calculated re-attachment lengths (Ud) are in very close agreement with the experimental values (which represent for Re < 150 the mean values of the experimental data of [24]). For Re > 150 no data were given by [24]. From the entries in Table 2 it is noticed that the computational times for QUICKE are approximately equal to those of UDS/HDS; however, the QUICKER schemes required twice these times to secure a converged solution. The fact that the calculated re-attachment lengths are very close to one another seems to

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

307

confirm the findings of Leschziner [S] that errors arising from false diffusion are small for this flow situation. To test further the performance of these schemes, a more severe test case is required: one that provides for large magnitudes of velocities crossing grid lines at large angles. It is for this reason that attention is turned to the flow in a square cavity with a moving lid.

4.4. ~u~~nu~~~~

in a square cavity with a sliding lid

4.4-f. The physical situation The physical situation considered is shown in Fig. 5(c) where it is noticed that the lid of the cavity slides with a constant velocity UTop above the cavity, the sides of which are fixed and of equal length L. A characteristic parameter governing the flow is the Reynolds number (p&&/~). The flow in the cavity, induced by the lid motion, was calculated at many values of Re; however, as a representative case, only that for Re = 400 is presented here. 4.4.2. Some computational details For this flow situation, various grid sizes were used; but only the results obtained using a 30 x 30 uniformly distributed grid will be concentrated upon. The convergence criteria used were more stringent that those previously applied because of the highly non-linear nature of the problem; that is, convergence criteria of lo-’ were used. The boundary conditions for this problem are: (i) top wail: U - UroP, V = 0; (ii) side walls: U = V = 0. 4.4.3. ~ese~tation and discussion of results This flow situation has become a standard ‘test-case’ for evaluating numerical procedures; although no standard solution has yet to be acknowledged (although many use 1261, see, for example, [ll]). In an effort to supply a standard solution, some of the participants at the First International Conference on Numerical Methods in Laminar and Turbulent Flows 1271 provided solutions for this problem using various finite-difference and finite-element techniques. The results of these calculations were collected, evaluated and presented by Olson [28]. The ‘best’ solutions (subjectively assessed with no experimental confirmation) was argued by Olson to be those of Cliffe et al. [29] who used a finite-element procedure on a 57 x 57 grid with extensive corner refinement. This, of course, is not to say that Cliffe’s et ai. solution is the true solution; for even within the results compiled by Olson, there are some solutions that compare very favourably with those obtained by Cliffe et al., and some others that did not compare well. Because the collection of results by Olson appears to be the most up-to-date, they are used for comparison purposes; but, keeping with previous investigators, Burggraf’s solution [26] is also used. Fig. 9 compares, with the numerical results gathered by Olson, the U-velocity profiles calculated at the cavity centre-line using a 30 X 30 grid. Since Olson presented a cluster of curves rather than a single line, it is the band width of these curves at all critical locations that are depicted in the figure. It is observed that only the QUICKEIER schemes provide

308

A. Pollard.

A.L.- W. Siu. Various

discretisutinns

--_-___ -+ i_

25

for Iuminar frmvs

UDS QIJICKUER HDS OLSON (1979)

“‘“TOP

Fig. 0. Laminar flow in a square cavity with a moving lid: U-velocity Re = 400.

profile along the vertical centre-plane.

30 X 30 grid,

calculations that are in close agreement with Olson’s results; and that due to the effect of false diffusion, neither the UDS nor the HDS were able to reproduce Olson’s results. Fig. 10 displays the vertically directed velocities along the horizontal centre-line of the cavity. Based upon the results presented in Fig. 9, it is expected that the QUICKE/ER schemes would out-perform the UDS and HDS; and indeed this is noted to be the case. From these two figures, it is clear that the QUICKE/ER schemes are best when compared to the results of Olson; providing, of course, that these results are indeed the true solutions. It is interesting to note that Olson did not consider as a base from which to work the results of Burgraff [26]. In fact, when compared to the results of [26] both Han et al. (using their version of QUICK) and the present QUICKE/ER results provide solutions that are imperceptibly different. This of course raises the question as to which solutions are indeed correct. All that can be said is that it is probably not those obtained using UDS/HDS. The static pressure distribution on the surface of the sliding lid is also provided by Olson. 0.4-

---

0.0

LiDS

0.2

0.6

0.4

0.8

1.0

X/L

Fig. 10. Laminar flow in a square 30 x 30 grid, Re = 400.

cavity with a moving

lid: V-velocity

profiles

along the horizontal

centre-plane,

A. Pollard, A.L.- W. Siu, Various discretisations

for laminar ftows

309

Results from two different research groups [29,30] were presented. These two pressure curves are shown in Fig. 13. It is obvious that the results obtained using the QUICKE/ER forms follow closely the pressure curve predicted by Cliffe et al. On the other hand, both the UDS and the HDS gave excellent predictions at both ends of the pressure curve, but the predicted pressure was too high in the region between x/L= 0.4 and x/L= 0.7.

50 fo)CLIFFE

Reb

eta,

(1978)

0

-50

-100

Ib)

c II

OLSON,

TUANN

x o-

CWICKE/ER HDS

A -

UDS

(19791

X/L

Fig. 11. Laminar flow in a square cavity with a moving lid: static pressure distribution Cliffe et al. and Olson and Tuann taken from [28].

along sliding lid; curves of

Table 3 tabulates the parameters for the primary vortex at Re = 400. The entries in Table 3 are defined in the following manner. The stream function r,Gis defined as:

u = -aqyay

(4.1)

Y= a@Jax

(4.2)

and

where, at the walls, J/ is arbitrarily set to zero. The vorticity < is defined as (X.Jfay - ~V/~x);

Table 3 Calculated

x

: t

Rep

parameters

for the primary vortex at Re = 400

UDS

HDS

0.628 -0.085 0.665

0.580 -0.096 0.620

2.0

-24.1

2.0

-31.4

QUICKE/ER 0.560 -0.107 0.607 2.2

-39.9

Olson

Burgraff

0.555 -0.114 0.605

0.565 -0.1017 0.615

2.3

2.142

-44.0

-

A. Pollard, A.L.- W. Siu, Various discretisations for larninar flows

310

and hence, is positive for clockwise rotation. The calculated l value can be treated as a measure of the rate of fluid rotation around the vortex centre. Finally, the dimensionless static pressure ji is given by (p - p,)lU & where pr is the reference pressure on and at the centre of the bottom wall. Note that the dimensionless group (Rep) is tabulated here rather than the dimensionless pressure. A local minimum of this pressure p (or Rep) denotes the x and y locations of the vortex centre. From the table, it is seen that QUICKWER give results which are in very close agreement with the results of both Olson and Burggraf. On the other hand, the results obtained by the UDS/HDS indicate less fluid is flowing across the surfaces between the vortex centre and the walls; that is, the fluid rotates at a slower rate. This behaviour is reflected by the lower vorticity values. As the fluid rotates slower around the vortex centre, the centrifugal force acting on it would be consequently weaker. Smaller pressure gradients are then required to maintain this circulatory motion. Hence, higher than usual static pressures appear at the vortex centre in the computed results of the UDS/HDS. Table 4 details the computational times and the number of iterations required to obtain converged results using UDS, HDS, QUICKE and QUICKER. It is seen that the QUICKE/ER schemes require about four times the number of iterations taken by the UDS and HDS. However, the results obtained by these latter two schemes are of little benefit since they are too far away from what are considered here as the acceptable solutions of Olson and Burggraf. In this test case, the UDS/HDS and the QUICKER forms were very stable even with the relaxation factors set at 0.8. However, the QUICKE form was relatively unstable. The solution diverged for a 10 x 10 grid, and convergence was secured for the finer grid only when the relaxation factors were lowered to 0.3.

Table 4 Comparison Re=400

of convergence

Schemes

10 x 10 grid

20 x 20 grid

30

x

30 grid

UDS HDS QUICKEb QUICKER UDS HDS QUICKE’ QUICKER UDS HDS QUICKE’ QUICKER

rates

for

square

TITUDS

cavity

flow

N

TIN”

1 0.83

69 87

0.145 0.139

3.94

198 181 254 577 586 392 472 1071 926

0.200 0.534 0.544 0.919 0.918 1.306 1.314 2.245 2.244

1.43 5.5 5.6 1.21 4.70 4.06

at

“T/N = actual computer time (secs)/number of iterations. bNo converged solutions were possible even with the relaxation factors equal to 0.3. ‘All relaxation factors equal to 0.3.

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

311

The computed results in this test case indicate that the QUICKE/ER formulations perform much better than the UDS/HDS. ~though the QUICKE form suffered from instability, its revised counterpart QUICKER performed well. 4._5. Summary and further discussion This section of the paper has presented the results of applying two newly developed forms derived from the QUICK scheme, namely the QUICKE and the QUICKER, in flows of varying complexity, Their performances have been checked with respect to the UDS and the X-IDS.QUICKE/ER were observed to give accurate and identical solutions in all the flow cases tested. The QUICKE form was found to be very efficient; even though it suffered from instability when the flow rate was high and the grid spacing was large (i.e., high grid Peclet numbers). The QUICKER form, on the other hand, was found to be very stable under all flow situations tested; but, because of the treatment of its source-like terms, the QUICKER form has been found to require more iterations in order to reach a converged solution. In the last two flow situations, where there are regions of recirculating fluid, it was noted that the number of iterations required by the QUICK formulations increased drastically over those needed to secure convergence using UDS and HDS. This behaviour, we believe, is probably caused by the solution procedure: it marches in a single direction, and during each sweep, only new information upstream of a line currently under consideration can be transmitted downstream. It is therefore easy to see that fewer iterations are required for the parallel plate problem. Coupled with this line-by-line approach, it is noted that for each line visited, the QUICKE/ER formulations need two points upstream or downstream as opposed to the UDS or HDS which needs only one point. This extra piece of information, which itself is not yet stabilized, we believe contributes significantly to the QUICKE/ER formulations slow rate of convergence. Before closing, and as a footnote to the calculations presented here, all the test cases have been repeated for a Reynolds number of lo”. The motivation for these calculations was provided in part by the fact that Han et al. [ll] could not obtain solutions for their QUICK formulated coefficients at this Reynolds number. These authors argue that “this is not surprising since it is unlikely that a stable, steady-flow solution. . . exists at such a high Reynold’s number”. Moreover, when calculating the flow in a sudden expansion they did obtain a converged solution for Re = 500 and 1000 when it is known that at Re > 2.50 a stable, steady-flow probably does not exist (see, for example, 1311).The present authors believe that a converged, stable steady-flow solution is theoretically possible, even though physically such a flow situation may not exist. The results of calculations at Re = 10’ showed that solutions (on coarse grid distributions) can be achieved for all flow situations considered here, using the UDS, HDS and QUICKER schemes but not for the QUICKE scheme. During the course of these calculations, it became clear that whether a solution could be reached seemed to depend upon the local grid Peclet number ratio (i.e., PeJPe,). For example: if this ratio was ~2.5, solutions to all schemes could not be obtained for the parallel plate case; if PeJPe, < 10 solutions for the sudden expansion flow case could not be obtained; whereas for the cavity problem, solutions could be obtained for the QUICKER scheme even if this ratio was less than 2. This behaviour is most disturbing and investigations that will hopefully clear up this matter are continuing.

312

5. Conclusions

A. Pollard,

A.L.- W. Siu, Various

discretisations

for laminar flows

and recommendations

The conclusions that can be drawn from this study are as follows. (i) Two new finite difference forms of the QUICK scheme have been derived; they are referred to here as the QUICK extended (QUICKE), and the QUICK extended and revised (QUICKER) forms respectively. (ii) These two forms have been successfully validated by applying the corresponding finite difference equations to solve for 2-D laminar flow problems of varying complexity. (iii) The QUICKE and QUICKER forms have been found to give solutions that are equal to or more accurate than the UDS and the HDS, especially in those situations where the flow is highly re-circulatory. (iv) In general, it has been found that in those situations where the attributes of the QUICKE/ER schemes are a most beneficial, computational effort exceeds that required for the UDS/HDS. (v) The QUICKE form was found to suffer from instability at high flow rates and coarse grids, while the QUICKER form has proved itself very stable for those flow situations considered. The finite difference QUICKER formulation of the QUICK scheme developed in this study has proved to be accurate and stable. However, this formulation did not give economical solutions in flow situations where there are zones of re-circulating fluid. It seems probable that this slow convergence is a result of the solution procedure. Therefore, it would be beneficial to repeat the tests with a whole field solutions approach, (i.e., use a procedure like the will known TEACH code [32]). In doing so, more computer storage locations would be required, but this solution method may improve substantially the iterative efficiency.

Acknowledgment The authors are glad to acknowledge the financial support provided by the Natural Science and Engineering Research Council of Canada under grant number A-1486. Also, thanks to the Department of Mechanical Engineering, University of Calgary for providing facilities while one of the authors (AP) was there on staff and the other (AS) was there as a graduate student. Thanks also to Miss J. Hill for typing the original manuscript.

References [l] G.D. Raithby, Skew upstream differencing schemes for problems involving fluid flow, Comput. Meths. Appl. Mech. Engrg. 9 (1975) 153-164. [2] B.P. Leonard, M. Leschziner and J. McGuirk, Third-order finite difference method for steady two-dimensional convection, Proc. Numerical Methods in Laminar and Turbulent Flows, Swansea, 1978. [3] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Meths. Appl. Mech. Engrg. 12 (1979) 59-98. [4] B.P. Leonard, Adjusted quadratic upstream algorithms for transient incompressible convection, Proc. AIAA Fourth Computational Fluid Dynamics Conf., Williamsburg, VA, 1979. [5] B.P. Leonard, A Survey of Finite Differences of Opinion on Numerical Mudling of the Incomprehensible Defective Confusion Equation (ASME WAM, New York, 1979).

A. Pollard, A.L.- W. Siu, Various discretisations for laminar flows

313

[6] B.P. Leonard, News flash: upstream parabolic interpolation, Proc. 2nd GAMM Conf. on Numerical Methods in Fluid Mechanics, Koln, Germany, 1977. [7] J. Militzer, W.B. Nicoll and S.A. Alpay, Some observation on the numerical calculation of the recirculation region of twin parallel symmetric jet flow, Proc. Symp. Turbulent Shear Flows, Penn State, 1977. [8] M.A. Leschziner, Practical evaluation of three finite difference schemes for the computation of steady-state recirculating flows, Comput. Meths. Appl. Mech. Engrg. 23 (1980) 293-312. [9] M.A. Leschziner and W. Rodi, Calculation of annular and twin parallel jets using various discretization schemes and turbulence model variants, ASME JFE 103 (1980) 352-360. [lo] I.P. Castro, The numerical prediction of recirculating flows, in: C. Taylor, K. Morgan, C.A. Brebbia, eds., Proc. Numerical Methods in Laminar and Turbulent Flows (Pentech Press, Plymouth, 1978). [l l] T. Han, J.A.C. Humphrey and B.E. Launder, A comparison of hybrid and quadratic upstream differencing in high Reynolds number elliptic flows, Rept. FM-81-1, University of California, 1981. Also: Comput. Meths. Appl. Mech. Engrg. 29 (1981) 81-95. [ 121 M.A. Leschziner, Private communication, 1980. [13] P.J. Roache, Computational Fluid Dynamics (Hermosa Publ. Albuquerque, NM, 1976). [14] D.B. Spalding, A novel finite-difference formulation for differential expressions involving both first and second derivatives, Internat. J. Numer. Meth. Engrg. 4 (1972) 551-559. [15] S.V. Patankar, Numerical Heat Transfer and Fluid Flow (McGraw-Hill, New York, 1980). [16] G.E. Forsythe and W.R. Wasaw, Finite-Difference Methods for Partial Differential Equations (Wiley, New York, 1960). [17] A. L.-W. Siu, Numerical calculation of some laminar flows using various discretization schemes, M.A.Sc. Thesis, Dept. of Mech. Engrg., University of Calgary, 1981. [18] S.V. Patankar and D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Internat. J. Heat Mass Transfer 15 (1972) 1787-1806. [19] W.M. Pun and D.B. Spalding, A general computer program for two-dimensional elliptic flows, Rept. HTS/76/2, Imperial College, London, 1976. [20] A. Pollard, Numerical computation of flow in tee-junctions, Physicochemical Hydrodynamics 2 (2/3) (1981) 203-227. [21] R.K. Shah and A.L. London, Laminar flow forced convection in ducts, Academic Press Adv. in Heat Transfer Series, Supp. 1 (1978). [22] J.W. McDonald, V.E. Denny and A.R. Mills, Numerical solutions of the Navier-Stokes equations in inlet regions, ASME J. Appl. Mech. 39 (1972) 873878. [23] H. Morihara and R.T. Cheng, Numerical solution of the viscous flow in the entrance region of parallel plates, J. Comp. Phys. 11 (1973) 550-572. [24] E.O. Macagno and T.K. Hung, Computational and experimental study of a captive annular eddy, J. Fluid Mech. 28 (1967) 43-61. [25] A. Pollard, A contribution on the effects of inlet conditions when modelling stenoses using sudden expansions, J. Biomechanics 14 (1981) 349-355. [26] O.R. Burggraf, Analytical and numerical studies of the structure of steady separated flows, J. Fluid Mech. 24 (1966) 113-151. [27] C. Taylor, K. Morgan and C.A. Brebbia, Eds., Numerical Methods in Laminar and Turbulent Flows (Pentech Press, Plymouth, 1978). [28] M.D. Olson, Comparison problem No. 1 recirculating flow in a square cavity, Structural Research Series Rept. 22, University of British Columbia, 1979. [29] K.A. Cliffe, C.P. Jackson, J. Rae and K.H. Winters, Finite element flow modelling using velocity and pressure variables, AERE Harewell Rept. No. R-9202, 1978. [30] M.D. Olson and S.Y. Tuann, New finite element results for the square cavity, Comput. Fluids 7 (1979) 123135. [31] L.H. Back and E.J. Roschke, The influence of upstream conditions on flow re-attachment lengths downstream of an abrupt circular channel expansion, J. Biomechanics 9 (1976) 481-483. [32] A.D. Gosman and W.M. Pun, The TEACH program, HTS Course Notes, Imperial College, London. 1974.